####################################################################################################
#####  R code to obtain the global and local semiparametric estimates of dose-response curve  ######
####################################################################################################

library(lokern)

# The pool-adjacent-violators algorithm (PAVA) to perform isotonic transformation for x, with weights wt.
pava <- function (x, wt=rep(1,length(x)))   
{ 
	n <- length(x) 
	if (n <= 1) return (x) 
	if (any(is.na(x)) || any(is.na(wt))) { 
		stop ("Missing values in 'x' or 'wt' not allowed") 
	} 
	lvlsets <- (1:n) 
	repeat { 
		viol <- (as.vector(diff(x)) < 0) # Find adjacent violators 
		if (!(any(viol))) break 
		i <- min( (1:(n-1))[viol]) # Pool first pair of violators 
		lvl1 <- lvlsets[i] 
		lvl2 <- lvlsets[i+1] 
		ilvl <- (lvlsets == lvl1 | lvlsets == lvl2) 
		x[ilvl] <- sum(x[ilvl]*wt[ilvl]) / sum(wt[ilvl]) 
		lvlsets[ilvl] <- lvl1 
	} 
	x 
} 

### data for example 2, a phase II dose finding study from Bretz, Pinheiro, and Branson (2005, Biometrics, 61, pp. 738)
phase2 <- read.table("phase2.txt", header=T)
x = unique(phase2$dose)   # unique dose levels
m = as.vector(table(phase2$dose))  # the number of subjects at each dose level

## convert the continuous outcome in the original dataset to a binary response variable y
cutoff=0.8 # cutoff to define binary response y
y = NULL
for(j in x) { y = c(y, sum(phase2$resp[phase2$dose==j]>cutoff)) }

y.m = cbind(y, m-y)
ndose = length(x)
grid = seq(min(x), max(x), length.out = 1000)  # define a fine grid over doses

### bootstrap to obtain MSEs and covariance to determine the weights for the semiparametric estimates
nboot = 1000
pa.bt = matrix(nrow=length(grid), ncol=nboot)
np.bt = matrix(nrow=length(grid), ncol=nboot)
y.b = rep(NA, ndose)
for(j in 1:nboot)
{
	for(i in 1:ndose) #bootstrap y independently at each dose level
	{
		resamples.ind <- sample(1:m[i], replace=T);
		y.b[i] <- sum(resamples.ind<=y[i]);
	}
	y.m.b = cbind(y.b, m-y.b)
	
	## parameteric fitting
	pafit.b <- glm(y.m.b~x, family=binomial(link="probit"))
	pa.bt[,j] <- predict(pafit.b, data.frame(x=grid), type="response")
	
	## nonparameteric fitting
	npfit.b <- glkerns(x, y.b/m, x.out=grid); 
	np.bt[,j] <- npfit.b$est;
}
npfit <- glkerns(x, y/m, x.out=grid)
true.est = npfit$est   # using the nonparametric estimate as the estimate of true dose-response curve
mse.pa = rowMeans((pa.bt-true.est)^2, na.rm=T)
mse.np = rowMeans((np.bt-true.est)^2, na.rm=T)
covar = rep(0, length(grid))
for(i in 1:length(grid))
{
	valid = (!is.na(pa.bt[i,])) &  (!is.na(np.bt[i,]))
	covar[i] = cov(pa.bt[i,valid], np.bt[i,valid])
}
covb = covar + rowMeans(pa.bt-true.est, na.rm=T)*rowMeans(np.bt-true.est, na.rm=T)

### calcualte weights
w.g = sum(mse.np-covb)/sum(mse.pa+mse.np-2*covb)  # global weights
w.l = (mse.np-covb)/(mse.pa+mse.np-2*covb)       # local weights

### regularize w when w>1 or <0
if(w.g>1) {w.g=1;}
if(w.g<0) {w.g=0;}
w.l[w.l>1]=1;
w.l[w.l<0]=0;

### dose-response estimates
pafit <- glm(y.m ~ x, family=binomial(link="probit"))  # parametri fit
est.pa = predict(pafit, data.frame(x=grid), type="response") # parametric estimate
est.np = npfit$est  # nonparametric estimate
est.gsp = w.g*est.pa + (1-w.g)*est.np  # global semiparametric estimate
est.lsp = w.l*est.pa + (1-w.l)*est.np  # local semiparametric estiamte


### plot dose-response curves
par(mfrow=c(1, 3))
#### plot (non-isotonic) estimates of dose-response curve
plot(grid, est.pa, type="l", ylim=c(0, 1), xlab="Dose", ylab="Probability of response", main="(a) Dose-response curve")
points(x, y/m)
lines(grid, est.np, lty=3, col=3)
lines(grid, est.gsp, lty=2, col=2)
lines(grid, est.lsp, lty=4, col=4)
legend(0.35, 1, lty=c(1, 3, 2, 4), col=c(1, 3, 2, 4), box.lwd=0, legend=c("Parametric", "Nonparametric", "Global semiparametric", "Local semiparametric"))

###### plot weights in semiparametric estimates #####
wwd=3
plot(grid, est.pa, type="l", ylim=c(0, 1), xlab="Dose", ylab="", main="(c) Weight")
lines(grid, est.np, lty=3)
abline(h=w.g, lwd=3, lty=3) 
lines(grid, w.l, lwd=3)
legend(0.4, 0.25, lty=c(3, 1, 1, 3), box.lwd=0, lwd=c(wwd, wwd, 1, 1), legend=c("Global weight", "Local weight", "Parametric estimate", "Nonparametric estimate"))


### isotonic estimates using PAVA
est.np.iso = pava(est.np);
est.gsp.iso = pava(est.gsp);
est.lsp.iso = pava(est.lsp);

#### plot isotonic estimates of dose-response curve
plot(grid, est.pa, type="l", ylim=c(0, 1), xlab="Dose", ylab="Probability of response", 
main="(b) Isotonic dose-response curve")
points(x, y/m)
lines(smooth.spline(grid, est.np.iso, spa=0.6), lty=3, col=3)
lines(smooth.spline(grid, est.gsp.iso, spa=0.6), lty=2, col=2)
lines(smooth.spline(grid, est.lsp.iso, spa=0.4), lty=4, col=4)
legend(0.35, 1, lty=c(1, 3, 2, 4),  col=c(1, 3, 2, 4), box.lwd=0, legend=c("Parametric", "Nonparametric", "Global semiparametric", "Local semiparametric"))