# For the dissolution 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 <- "dissresid.Rout" # Enter the data thedata <- scan("dissolution.dat") thedata <- matrix(thedata,ncol=3,byrow=T) v <- thedata[,1] x <- thedata[,2] y <- thedata[,3] n <- length(x) p <- 6 # Enter the estimates # OLS bols <- bols sigols <- sigols # GLS bgls <- bgls theta <- theta siggls <- sigma # Mean function ffunc <- function(t,W0,W1,beta0,beta1,tau0,tau1,v){ log2 <- log(2) t0 <- 20 f10 <- W0*(1-exp(-log2*(t/tau0)^beta0)) joinpt <- W0*(1-exp(-log2*(t0/tau0)^beta0)) f11 <- (t>t0)*(joinpt + W1*(1-exp(-log2*((t-t0)*(t>t0)/tau1)^beta1))) f1 <- (1-v)*f10 + v*( (t<=t0)*f10 + (t>t0)* f11 ) f1 } # The gradient matrix X(beta) fbfunc <- function(t,W0,W1,beta0,beta1,tau0,tau1,v){ log2 <- log(2) t0 <- 20 stuff <- (t+0.1*(t==20)-t0)/tau1 piece <- log(abs(stuff)) meangrad <- array(0,c(length(t),6),list(NULL,c("W0","W1","beta0","beta1", "tau0","tau1" ))) meangrad[,"W0"] <- ( (1-v)*(1-exp(-log2*(t/tau0)^beta0)) + v*( (t<=t0)*(1-exp(-log2*(t/tau0)^beta0)) + (t>t0)*(1-exp(-log2*(t0/tau0)^beta0)) )) meangrad[,"W1"] <- v*(t>t0)*(1-exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)) meangrad[,"beta0"] <- ( (1-v)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(log(t/tau0))*log2 + v*( (t<=t0)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(log(t/tau0))*log2 + (t>t0)*W0*exp(-log2*(t0/tau0)^beta0)*((t0/tau0)^beta0)*(log(t0/tau0))*log2 ) ) meangrad[,"beta1"] <- v*(t>t0)*W1*exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)*(((t-t0)*(t>t0)/tau1)^beta1)*piece*log2 meangrad[,"tau0"] <- -( (1-v)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(beta0/tau0)*log2 + v*( (t<=t0)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(beta0/tau0)*log2 + (t>t0)*W0*exp(-log2*(t0/tau0)^beta0)*((t0/tau0)^beta0)*(beta0/tau0)*log2 ) ) meangrad[,"tau1"] <- -v*(t>t0)*W1*exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)*(((t-t0)*(t>t0)/tau1)^beta1)*(beta1/tau1)*log2 meangrad } # The "weighted" gradient matrix X*(beta) fbwfunc <- function(t,W0,W1,beta0,beta1,tau0,tau1,v,mut){ w12<-1/mut log2 <- log(2) t0 <- 20 stuff <- (t+0.1*(t==20)-t0)/tau1 piece <- log(abs(stuff)) meangrad <- array(0,c(length(t),6),list(NULL,c("W0","W1","beta0","beta1", "tau0","tau1" ))) meangrad[,"W0"] <- w12*( (1-v)*(1-exp(-log2*(t/tau0)^beta0)) + v*( (t<=t0)*(1-exp(-log2*(t/tau0)^beta0)) + (t>t0)*(1-exp(-log2*(t0/tau0)^beta0)) )) meangrad[,"W1"] <- w12*v*(t>t0)*(1-exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)) meangrad[,"beta0"] <- w12*( (1-v)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(log(t/tau0))*log2 + v*( (t<=t0)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(log(t/tau0))*log2 + (t>t0)*W0*exp(-log2*(t0/tau0)^beta0)*((t0/tau0)^beta0)*(log(t0/tau0))*log2 ) ) meangrad[,"beta1"] <- w12*v*(t>t0)*W1*exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)*(((t-t0)*(t>t0)/tau1)^beta1)*piece*log2 meangrad[,"tau0"] <- -w12*( (1-v)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(beta0/tau0)*log2 + v*( (t<=t0)*W0*exp(-log2*(t/tau0)^beta0)*((t/tau0)^beta0)*(beta0/tau0)*log2 + (t>t0)*W0*exp(-log2*(t0/tau0)^beta0)*((t0/tau0)^beta0)*(beta0/tau0)*log2 ) ) meangrad[,"tau1"] <- -w12*v*(t>t0)*W1*exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)*(((t-t0)*(t>t0)/tau1)^beta1)*(beta1/tau1)*log2 meangrad } xmat <- as.matrix(fbfunc(x,bols[1],bols[2],bols[3],bols[4],bols[5],bols[6],v)) # 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],bols[5],bols[6],v) 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("dissolsresid.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],bgls[5],bgls[6],v) predt <- pred^theta xwmat <- as.matrix(fbwfunc(x,bgls[1],bgls[2],bgls[3],bgls[4],bgls[5],bgls[6],v,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("dissglsresid.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")