#  For the assay data, compute OLS and weighted studentized
#  residuals and construct diagnostic plots.  The weighted versions
#  are based on the GLS-PL fit

#  Also, compute standard errors based on the "folklore" theorem
#  and "robust" standard errors using the "sandwich" estimator

#  Assumes you have run a program to compute the OLS and GLS-PL 
#  estimates already

outfile <- "assayresid.Rout"

#  Enter the data

thedat <- scan("assay.dat")
thedat <- matrix(thedat,ncol=3,byrow=T)
y <- thedat[,3]
x <- thedat[,1]

n <- length(x)

p <- 4

#  Enter the estimates 

#  OLS

bols <- bols
sigols <- sigols

#  GLS

bgls <- bgls
theta <-  theta
siggls <- sigma

#  Mean function

ffunc <- function(x,b1,b2,b3,b4){

#  note that we must be careful at x=0 -- the function approaches
#  b2 as x goes to 0

   xx <-(x>0)*x+(x==0)*0.0001
   denom <- 1+exp(b4*(log(xx)-b3))
   f1 <- b1 + (b2-b1)/denom
   fpl <- (x>0)*f1+(x==0)*b2
   fpl
}

#  The gradient matrix X(beta)

fbfunc <- function(x,b1,b2,b3,b4){
 
   xx <-(x>0)*x+(x==0)*0.0001
   denom <- 1+exp(b4*(log(xx)-b3))

   meangrad <- array(0,c(length(x),4),list(NULL,c("b1","b2","b3","b4")))
   meangrad[,"b1"] <- (x>0)*(1-1/denom)+(x==0)*0
   meangrad[,"b2"] <- (x>0)*(1/denom)+(x==0)*1
   meangrad[,"b3"] <- (x>0)*(b4*(b2-b1)*exp(b4*(log(xx)-b3)))/(denom**2)
        +(x==0)*0
   meangrad[,"b4"] <- -(x>0)*(b2-b1)*
          exp(b4*(log(xx)-b3))*(log(xx)-b3)/(denom**2)+(x==0)*0 
    meangrad
}

#  The "weighted" gradient matrix X*(beta)

fbwfunc <- function(x,b1,b2,b3,b4,mut){
    
   w12<-1/mut
   xx <-(x>0)*x+(x==0)*0.0001
   denom <- 1+exp(b4*(log(xx)-b3))

   meangrad <- array(0,c(length(x),4),list(NULL,c("b1","b2","b3","b4")))
   meangrad[,"b1"] <- ((x>0)*(1-1/denom)+(x==0)*0)*w12
   meangrad[,"b2"] <- ((x>0)*(1/denom)+(x==0)*1)*w12
   meangrad[,"b3"] <- ((x>0)*(b4*(b2-b1)*exp(b4*(log(xx)-b3)))/(denom**2)
        +(x==0)*0)*w12
   meangrad[,"b4"] <- (-(x>0)*(b2-b1)*
          exp(b4*(log(xx)-b3))*(log(xx)-b3)/(denom**2)+(x==0)*0 )*w12
    meangrad
}
	
xmat <- as.matrix(fbfunc(x,bols[1],bols[2],bols[3],bols[4]))

#  Compute the "hat" matrix for OLS

xtx <- t(xmat)%*%xmat
ixtx <- solve(xtx)
hmat <- xmat%*%ixtx%*%t(xmat)
seols <- sigols*sqrt(diag(ixtx))

#  The OLS leverages

hjj <- diag(hmat)

#  (a) Predicted values, residuals, studentized residuals, etc

pred <- ffunc(x,bols[1],bols[2],bols[3],bols[4])	
resid <- y-pred
bresid <- resid/(sigols*sqrt(1-hjj))
b23 <- (bresid**2)**(1/3)
logb <- log(abs(bresid))
logr <- log(abs(resid))

#  Make a four-panel plot

postscript("assayolsresid.ps",width=9)

par(mfrow=c(2,2),cex=0.75,oma=c(0,0,2,0),mar=c(4,4,3,2))

plot(pred,resid/sigols,type="p",ylab="Resid.",xlab="Log Pred.",
        cex=0.75,log="x",
        pch=18,mkh=0.06)
	lines(c(0.001,pred+10,800),rep(0,n+2),type="l",lty=1)
plot(pred,bresid,type="p",ylab="Stud. Resid.",xlab="Log Pred.",
        cex=0.75,log="x",
        pch=18,mkh=0.06)
	lines(c(0.001,pred+10,5),rep(0,n+2),type="l",lty=1)
plot(pred,abs(bresid),type="p",ylab="Log Stud. Resid.",xlab="Log Pred.",
        cex=0.75,log="xy",
        pch=18,mkh=0.06)
plot(pred,b23,type="p",ylab="2/3 Root Stud. Resid.",xlab="Log Pred.",
        cex=0.75,log="x",
        pch=18,mkh=0.06)

#  (b) Now do the weighted versions based on GLS-PL

pred <- ffunc(x,bgls[1],bgls[2],bgls[3],bgls[4])
predt <- pred**theta

xwmat <- as.matrix(fbwfunc(x,bgls[1],bgls[2],bgls[3],bgls[4],predt))

wxtx <- t(xwmat)%*%xwmat
iwxtx <- solve(wxtx)
hwmat <- xwmat%*%iwxtx%*%t(xwmat)

hwjj <- diag(hwmat)

resid <- y-pred
wresid <- resid/(pred**theta)
bresid <- wresid/(siggls*sqrt(1-hwjj))
b23 <- (bresid**2)**(1/3)
logb <- log(abs(bresid))
logr <- log(abs(wresid))

postscript("assayglsresid.ps",width=9)

par(mfrow=c(2,2),cex=0.75,oma=c(0,0,2,0),mar=c(4,4,3,2))
plot(pred,wresid/siggls,type="p",ylab="Wt. Resid.",xlab="Log Pred.",
        cex=0.75,log="x",
        pch=18,mkh=0.06)
	lines(c(0.001,pred+10,800),rep(0,n+2),type="l",lty=1)
plot(pred,bresid,type="p",ylab="Stud. Resid.",xlab="Log Pred.",
        cex=0.75,log="x",
        pch=18,mkh=0.06)
lines(c(0.001,pred+10,150),rep(0,n+2),type="l",lty=1)
plot(pred,abs(bresid),type="p",ylab="Log Stud. Resid.",xlab="Log Pred.",
       cex=0.75,log="xy",
        pch=18,mkh=0.06)
plot(pred,b23,type="p",ylab="2/3 Root Stud. Resid.",xlab="Log Pred.",
       cex=0.75,log="x",
        pch=18,mkh=0.06)

graphics.off()

#  Compute OLS results

cat("OLS estimate of beta and usual standard errors",file=outfile,
  append=F,"\n","\n")
cat(round(bols,3),file=outfile,append=T,"\n")
cat(round(seols,3),file=outfile,append=T,"\n","\n","\n")

#  (c) Compute GLS-PL standard errors using the "folklore" theorem

covmat <- siggls^2*iwxtx

se <- sqrt(diag(covmat))

cat("GLS estimate of beta and usual standard errors",file=outfile,
  append=T,"\n","\n")
cat(round(bgls,3),file=outfile,append=T,"\n")
cat(round(se,3),file=outfile,append=T,"\n","\n","\n")

#  (d) Compute "robust sandwich" standard errors

r2 <- resid^2
rmat <- diag(r2,nrow=length(r2))

#  compute the "center" piece -- we need the matrix WX to form
#  the center piece, so construct it from what we already have 
#  computed

wt12 <- 1/predt
wmat12 <- diag(wt12,nrow=length(wt12))
xmatu <- wmat12%*%xwmat

sterm <- t(xmatu)%*%rmat%*%(xmatu)

sandwich <- iwxtx%*%sterm%*%iwxtx
robustse <- sqrt(diag(sandwich))

cat("Robust sandwich standard errors",file=outfile,append=T,"\n","\n")

cat(round(robustse,3),file=outfile,append=T,"\n","\n")

#  estimates, standard errors -- use usual ses, GLS

#  plot with ols and gls fits superimposed

postscript("assaydata.ps",width=9)

xgrid <- seq(0,90,0.05)
fols <- ffunc(xgrid,bols[1],bols[2],bols[3],bols[4])
fgls <- ffunc(xgrid,bgls[1],bgls[2],bgls[3],bgls[4])

plot(x,y,type="p",xlab="concentration", ylab="response",
   xlim=c(0,90))
lines(xgrid,fols,lty=2)
lines(xgrid,fgls,lty=1)

graphics.off()