#### R Code for 'Probability-scale residuals for continuous, discrete, and censored data' ##################################### ##### Exponential linear model ##################################### rm(list=ls()) library(survival) n<-250 set.seed(400) x<-rnorm(n) y<-rexp(n,exp(-x+x^2)) #### Generating data under exponential model x2<-x^2 mod.no.quadratic<-survreg(Surv(y,rep(1,n))~x,dist="exponential") mod.quadratic<-survreg(Surv(y,rep(1,n))~x+x2,dist="exponential") #### correct model ## Probability-scale residuals for each model ry.no.quadratic<-2*pexp(y,1/predict(mod.no.quadratic))-1 ry.quadratic<-2*pexp(y,1/predict(mod.quadratic))-1 ## Observed-minus-expected residuals for each model r.naive.no.quadratic<-residuals(mod.no.quadratic) r.naive.quadratic<-residuals(mod.quadratic) ## Deviance residuals for each model res.devi1<-residuals(mod.no.quadratic,type="deviance") res.devi2<-residuals(mod.quadratic,type="deviance") ## Log(Observed)-Log(Expected) log.omer1<-log(y)-log(predict(mod.no.quadratic)) log.omer2<-log(y)-log(predict(mod.quadratic)) ##### Figure 1 setEPS(horizontal = TRUE) postscript("quadratic-vs-linear-plus4-a.eps",height=4.3,width=8.2) layout(matrix(c(1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5, 1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5, 1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5, 1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5, 6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10, 6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10, 6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10, 6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10), ncol=17, byrow=TRUE)) par(mar=c(4,3,1,0.5)) plot(c(0,1),c(0,1),type="n",xlab="",ylab="",axes=FALSE) mtext("Quadratic Model Fit",side=2,line=1,cex=.8) plot(x,ry.quadratic,ylab="",xlab="",xlim=c(-2.7,2.7),col=1) lines(lowess(ry.quadratic~x,f=1),col=gray(.5),lwd=3) mtext("z", side=1, line=2.2, cex=.65) mtext("PSR", side=2, line=2, cex=.65) plot(x,res.devi2,ylab="",xlab="",xlim=c(-2.7,2.7)) lines(lowess(res.devi2~x,f=1),col=gray(.5),lwd=3) mtext("z", side=1, line=2.2, cex=.65) mtext("Deviance Residual", side=2, line=2, cex=.65) plot(x,r.naive.quadratic,ylab="",xlab="",xlim=c(-2.7,2.7)) lines(lowess(r.naive.quadratic~x,f=1),col=gray(.5),lwd=3) mtext("z", side=1, line=2.2, cex=.65) mtext("OMER", side=2, line=2, cex=.65) plot(x,log.omer2,ylab="",xlab="",xlim=c(-2.7,2.7)) lines(lowess(log.omer2~x,f=1),col=gray(.5),lwd=3) mtext("z", side=1, line=2.2, cex=.65) mtext("log(Observed) - log(Expected)", side=2, line=2, cex=.65) plot(c(0,1),c(0,1),type="n",xlab="",ylab="",axes=FALSE) #text(.5,.5,"Linear Model Fit",srt=90,cex=.8) mtext("Linear Model Fit",side=2,line=1,cex=.8) plot(x,ry.no.quadratic,ylab="",xlab="",xlim=c(-2.7,2.7),col=1) lines(lowess(ry.no.quadratic~x,f=1),col=gray(.5),lwd=3) mtext("z", side=1, line=2.2, cex=.65) mtext("PSR", side=2, line=2, cex=.65) plot(x,res.devi1,ylab="",xlab="",xlim=c(-2.7,2.7)) lines(lowess(res.devi1~x,f=1),col=gray(.5),lwd=3) mtext("z", side=1, line=2.2, cex=.65) mtext("Deviance Residual", side=2, line=2, cex=.65) plot(x,r.naive.no.quadratic,ylab="",xlab="",xlim=c(-2.7,2.7)) lines(lowess(r.naive.no.quadratic~x,f=1),col=gray(.5),lwd=3) mtext("z", side=1, line=2.2, cex=.65) mtext("OMER", side=2, line=2, cex=.65) plot(x,log.omer1,ylab="",xlab="",xlim=c(-2.7,2.7)) lines(lowess(log.omer1~x,f=1),col=gray(.5),lwd=3) mtext("z", side=1, line=2.2, cex=.65) mtext("log(Observed) - log(Expected)", side=2, line=2, cex=.65) dev.off() ################################ #### Biomarker example #### Semi-parametric transformation model; Figure 2. ################################# rm(list=ls()) data=read.csv('AIACK23Study_DATA_2016-01-12_1340.csv') data$study_arm.factor = factor(data$study_arm,levels=c("1","2","3")) levels(data$study_arm.factor)=c("HIV-infected non-obese","HIV-infected obese","Uninfected obese") data$race.factor = factor(data$race,levels=c("0","1","2","3","4")) levels(data$race.factor)=c("caucasian","african american","hispanic","asian","other") library(rms) presid.orm <- function(object){ k <- object$non.slopes L <- object$linear.predictors trans <- object$trans cumprob <- if(length(trans)) trans$cumprob else plogis if(length(L)==0) stop('you did not use linear.predictors=TRUE for the fit') if(length(Y <- object$y) == 0) stop("you did not specify y=TRUE in the fit") if(!is.factor(Y)) Y <- factor(Y) Y <- unclass(Y) - 1 cof <- object$coefficients if(length(X <- object$x)==0) stop("you did not use x=TRUE in the fit") N <- length(Y) px <- 1 - cumprob(outer(cof[1:k], as.vector(X %*% cof[- (1:k)]), "+")) low.x = rbind(0, px)[cbind(Y + 1L, 1:N)] hi.x = 1 - rbind(px, 1)[cbind(Y + 1L, 1:N)] return(low.x - hi.x) } d<-data d$bmi<-data$bmi_screening d$cd4<-data$cd4_enroll d$art.duration<-data$duration_art study.arm<-data$study_arm.factor d$hiv<-ifelse(study.arm=="Uninfected obese",0,1) d$race<-with(d,ifelse(race>1,1,race)) dd<-with(d, datadist(age,sex,race,bmi,cd4,art.duration,smoker)) options(datadist='dd') d$y<-d$oa_akg summary(d$y) postscript("quantile-plot-biomarker.eps",width=10,height=3.5) par(mfrow=c(1,3)) mod1<-lm(y~age+sex+race+bmi+cd4+art.duration+smoker,data=d,subset=hiv==1) psr1<-presid(mod1) qqplot(qunif(ppoints(length(psr1)),-1,1),psr1,xlab="Quantiles for Uniform(-1,1)", ylab="PSR Quantiles", main="Linear Model of Y",cex.lab=1.4,cex.axis=1.4,ylim=c(-1,1)) abline(0,1) d$log.y<-log(d$y) mod2<-lm(log.y~age+sex+race+bmi+cd4+art.duration+smoker,data=d,subset=hiv==1) psr2<-presid(mod2) qqplot(qunif(ppoints(length(psr2)),-1,1),psr2,xlab="Quantiles for Uniform(-1,1)", ylab="PSR Quantiles", main="Linear Model of log(Y)",cex.lab=1.4,cex.axis=1.4,ylim=c(-1,1)) abline(0,1) mod<-orm(y~age+sex+race+bmi+cd4+art.duration+smoker,x=TRUE,y=TRUE,data=d,subset=hiv==1,family="probit") psr<-presid.orm(mod) qqplot(qunif(ppoints(length(psr)),-1,1),psr,xlab="Quantiles for Uniform(-1,1)", ylab="PSR Quantiles", main="Semiparametric Transformation Model",cex.lab=1.4,cex.axis=1.4,ylim=c(-1,1)) abline(0,1) dev.off() postscript("quantile-plot-biomarker-omer.eps",width=10,height=3.5) par(mfrow=c(1,3)) omer1<-residuals(mod1) qqnorm(omer1, main="Linear Model of Y") qqline(omer1) d$log.y<-log(d$y) omer2<-residuals(mod2) qqnorm(omer2, main="Linear Model of log(Y)") qqline(omer2) y<-d$y[d$hiv==1] orm.omer<- -mod$coefficients[sapply(y, FUN=function(x) which(mod$yunique==x))] - (predict(mod, type="lp", kint=1) - mod$coefficients[1]) #### can not compute OMER for the largest y because \alpha(y) is unbounded. orm.omer[which(y==max(y))] <- NA qqnorm(orm.omer, main="Semiparametric Transformation Model") qqline(orm.omer) dev.off() ############################### #### Quantile Regression -- Mixture of Normals; Figure 3. ############################### rm(list=ls()) require(quantreg) set.seed(20160502) n <- 100 beta0 <- 0 beta1 <- -1 beta2 <- 1 mu.outliers <- 0 sigma.outliers <-100 x <- rnorm(n,0,1) s <- rbinom(n, 1, 0.9) y <- s*rnorm(n, beta0+beta1*x+beta2*x^2, 1) + (1-s)*rnorm(n, mu.outliers, sigma.outliers) newRS <- function(m){ apply(m$residuals,1, FUN=function(x){2*mean(m$tau[which(abs(x)==min(abs(x)))])-1}) } ###median QR x2<-x^2 m1 <- rq(y~x, tau=0.5) m2 <- rq(y~x+x2, tau=0.5) taus <- (1:99)/100 f1 <- rq(y~x, tau=taus) f2 <- rq(y~x+x2, tau=taus) library(mixtools) X<-cbind(x,x2) mod.mixture2<-regmixEM(y,X, lambda=c(.95,.05), beta=matrix(c(beta0, beta1, beta2, 0, 0, 0), 3, 2), sigma=c(1,100), addintercept=TRUE) X1<-cbind(1,X) psr.mod.mixture2<-2*(mod.mixture2$lambda[1]*pnorm(y,mod.mixture2$beta[,1]%*%t(X1),mod.mixture2$sigma[1])+ mod.mixture2$lambda[2]*pnorm(y,mod.mixture2$beta[,2]%*%t(X1),mod.mixture2$sigma[2])) - 1 qqplot(qunif(ppoints(length(y)),-1,1),psr.mod.mixture2,xlab="Quantiles for Unif(-1,1)", ylab="PSR with quadratic term") abline(0,1) mod.mixture1<-regmixEM(y,x, lambda=c(.95,.05), beta=matrix(c(beta0, beta1, 0, 0), 2, 2), sigma=c(1,100), addintercept=TRUE) W1<-cbind(1,x) psr.mod.mixture1<-2*(mod.mixture1$lambda[1]*pnorm(y,mod.mixture1$beta[,1]%*%t(W1),mod.mixture1$sigma[1])+ mod.mixture1$lambda[2]*pnorm(y,mod.mixture1$beta[,2]%*%t(W1),mod.mixture1$sigma[2])) - 1 ### PSR only using the information from median regression med.res1<-ifelse(m1$residual<0,-.5,.5) med.res2<-ifelse(m2$residual<0,-.5,.5) ##### Figure 3 setEPS(horizontal = TRUE) f<-.95 postscript("quant-reg-plots-mixture-curiosity-a.eps",height=4.3,width=8.2) layout(matrix(c(1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5, 1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5, 1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5, 1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5, 6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10, 6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10, 6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10, 6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10), ncol=17, byrow=TRUE)) par(mar=c(4,3.1,1,0.5)) plot(c(0,1),c(0,1),type="n",xlab="",ylab="",axes=FALSE) mtext("Quadratic Model Fit",side=2,line=1,cex=.8) plot(x, m2$residuals, ylab="",xlab="") lines(lowess(x,m2$residuals,f=f),col=gray(.6)) mtext("z", side=1, line=2.2, cex=.65) mtext("Observed - Median", side=2, line=2, cex=.65) plot(c(min(x),max(x)),c(-1,1),type="n",ylab="",xlab="") points(x, med.res2) lines(lowess(x,med.res2,f=f),col=gray(.6)) mtext("z", side=1, line=2.2, cex=.65) mtext("PSR (Median Regression)", side=2, line=2, cex=.65) plot(x, newRS(f2), ylab="",xlab="") lines(lowess(x,newRS(f2),f=f),col=gray(.6)) mtext("z", side=1, line=2.2, cex=.65) mtext("PSR (Multiple QR)", side=2, line=2, cex=.65) plot(x, psr.mod.mixture2,ylab="",xlab="") lines(lowess(x,psr.mod.mixture2,f=f),col=gray(.6)) mtext("z", side=1, line=2.2, cex=.65) mtext("PSR (Mixture of Normals)", side=2, line=2, cex=.65) plot(c(0,1),c(0,1),type="n",xlab="",ylab="",axes=FALSE) mtext("Linear Model Fit",side=2,line=1,cex=.8) plot(x, m1$residuals, ylab="",xlab="") lines(lowess(x,m1$residuals,f=f),col=gray(.6)) mtext("z", side=1, line=2.2, cex=.65) mtext("Observed - Median", side=2, line=2, cex=.65) plot(c(min(x),max(x)),c(-1,1),type="n",ylab="",xlab="") points(x, med.res1) lines(lowess(x,med.res1,f=f),col=gray(.6)) mtext("z", side=1, line=2.2, cex=.65) mtext("PSR (Median Regression)", side=2, line=2, cex=.65) plot(x, newRS(f1),ylab="",xlab="") lines(lowess(x,newRS(f1),f=f),col=gray(.6)) mtext("z", side=1, line=2.2, cex=.65) mtext("PSR (Multiple QR)", side=2, line=2, cex=.65) plot(x, psr.mod.mixture1,ylab="",xlab="") lines(lowess(x,psr.mod.mixture1,f=f),col=gray(.6)) mtext("z", side=1, line=2.2, cex=.65) mtext("PSR (Mixture of Normals)", side=2, line=2, cex=.65) dev.off() ################################# ### Count Data Examples ################################# ############## Figure 4 ########## Comparing residuals when data generated from Poisson or Negative Binomial and properly or improperly fit rm(list=ls()) library(MASS) set.seed(20140331) n <- 2000 beta0 <- 0 beta1 <- 1 x <- rnorm(n) mu <- beta0 + beta1*x ##### Poisson y <- rpois(n, lambda=exp(mu)) m1 <- glm(y~x, family=poisson()) s1 <- summary(m1) p.resid <- ppois(y-1, m1$fitted.values) + ppois(y, m1$fitted.values) -1 deviance.resid <- resid(m1, type="dev") pearson.resid <- resid(m1, type="pear") ##### Negative binomial phi <- 3 y3 <- rnbinom(n, mu=exp(mu), size=phi) #### NB2 model is implemented in glm.nb nb2.nb <- glm.nb(y3~x) s.nb2.nb <- summary(nb2.nb) nb2.nb.presid <- pnbinom(y3-1, mu=nb2.nb$fitted.values, size=nb2.nb$theta ) + pnbinom(y3, mu=nb2.nb$fitted.values, size=nb2.nb$theta) -1 nb2.nb.deviance <- resid(nb2.nb, type="deviance") nb2.nb.pearson <- resid(nb2.nb, type="pearson") #### Now fitting wrong model (poisson when it should be negative binomial) nb2.poi <- glm(y3~x, family="poisson") nb2.poi.presid <- ppois(y3-1, nb2.poi$fitted.values) + ppois(y3, nb2.poi$fitted.values) -1 nb2.poi.deviance <- resid(nb2.poi, type="deviance") nb2.poi.pearson <- resid(nb2.poi, type="pearson") ##### Figure 4 b<-3 postscript("pscale-poisson-nb.eps",height=8, width=8) par(mfrow=c(3,3), par(mar=c(3,4,2,0.1)), mgp = c(1.5, 0.5, 0)) #, xpd = NA, oma=c(0,4,2, 0)) plot(x, p.resid, cex=0.4, xlab="z", ylab="") abline(h=0, lty=2) lines(supsmu(x, p.resid, bass=b), col="gray70", lwd=3) mtext("Probability-Scale Residual", side=3, line=0.5) mtext("A", side=2, line=2, at=1, las=1, cex=1.6) plot(x, deviance.resid, cex=0.4, xlab="z", ylab="") abline(h=0, lty=2) lines(supsmu(x, deviance.resid, bass=b), col="gray70", lwd=3) mtext("Deviance Residual", side=3, line=0.5) plot(x, pearson.resid, cex=0.4, xlab="z", ylab="") abline(h=0, lty=2) lines(supsmu(x, pearson.resid, bass=b), col="gray70", lwd=3) mtext("Pearson Residual", side=3, line=0.5) plot(x, nb2.nb.presid, cex=0.4, xlab="z", ylab="") abline(h=0, lty=2) lines(supsmu(x, nb2.nb.presid, bass=b), col="gray70", lwd=3) mtext("B", side=2, line=2, at=1, las=1, cex=1.6) plot(x, nb2.nb.deviance, cex=0.4, xlab="z", ylab="") abline(h=0, lty=2) lines(supsmu(x, nb2.nb.deviance, bass=b), col="gray70", lwd=3) plot(x, nb2.nb.pearson, cex=0.4, xlab="z", ylab="") abline(h=0, lty=2) lines(supsmu(x, nb2.nb.pearson, bass=b), col="gray70", lwd=3) plot(x, nb2.poi.presid, cex=0.4, xlab="z", ylab="") abline(h=0, lty=2) lines(supsmu(x, nb2.poi.presid, bass=b), col="gray70", lwd=3) mtext("C", side=2, line=2, at=1, las=1, cex=1.6) plot(x, nb2.poi.deviance, cex=0.4, xlab="z", ylab="") abline(h=0, lty=2) lines(supsmu(x, nb2.poi.deviance, bass=b), col="gray70", lwd=3) plot(x, nb2.poi.pearson, cex=0.4, xlab="z", ylab="") abline(h=0, lty=2) lines(supsmu(x, nb2.poi.pearson, bass=b), col="gray70", lwd=3) dev.off() ################################## ### Survival Example ################################## rm(list=ls()) library(survival) library(rms) library(mvtnorm) load("d-second-chunk-psresiduals.Rda") #### We do not have permission to include the CCASAnet data used for this analysis. #### Note to Bryan: File is in gabriela directory. d$init.class.cat<-ifelse(d$init.class=="3 NRTI"|d$init.class=="OTHER"|d$init.class=="PI"|d$init.class=="unboosted PI", "Other",as.character(d$init.class)) d$sqrt.cd4<-sqrt(d$nadir.cd4.prehaart) d1<-d[!is.na(d$sqrt.cd4)&d$init.class.cat!="non-HAART"&d$male==0&d$age>50,] dim(d1) median(d1$time)/365.25 sum(d1$death1) sum(d1$death1)/length(d1$death1) mod<-coxph(Surv(time,death1)~rcs(age,4)+aids.enrollment+rcs(sqrt.cd4,4)+rcs(year.init,4)+init.class.cat+strata(site), data=d1) mod1<-coxph(Surv(time,death1)~rcs(age,4)+aids.enrollment+rcs(year.init,4)+init.class.cat+strata(site), data=d1) mresids<-residuals(mod) mresids1<-residuals(mod1) delta<-d1$death1 presids<-1-exp(mresids-delta)-delta*exp(mresids-delta) presids1<-1-exp(mresids1-delta)-delta*exp(mresids1-delta) ############# Extra stuff added in Feb 2016: csresids<-delta-mresids csresids1<-delta-mresids1 fit<-survfit(Surv(csresids,delta)~1) fit1<-survfit(Surv(csresids1,delta)~1) presids.csa<-2*(1-exp(-csresids))-1 presids1.csa<-2*(1-exp(-csresids1))-1 Sy<- exp(-csresids) Sy1<- exp(-csresids1) presids.cs<-ifelse(delta==0,presids-Sy*(1+delta),presids) presids1.cs<-ifelse(delta==0,presids1-Sy1*(1+delta),presids1) psr.norm<-qnorm((presids+1)/2) psr1.norm<-qnorm((presids1+1)/2) dres<-ifelse(mresids>0,1,-1)*sqrt(-2*(mresids+delta*log(delta-mresids))) dres1<-ifelse(mresids1>0,1,-1)*sqrt(-2*(mresids1+delta*log(delta-mresids1))) ##### Big plot with all of them; Figure 5 setEPS(horizontal = TRUE) postscript("time-to-event-big.eps",width=7,height=10) par(mar=c(4,3,.5,.5)) layout(matrix(c(1,1,1,1,1,1,1,1,1, 2,2,2,3,3,3,4,4,4, 2,2,2,3,3,3,4,4,4, 2,2,2,3,3,3,4,4,4, 5,5,5,6,6,6,7,7,7, 5,5,5,6,6,6,7,7,7, 5,5,5,6,6,6,7,7,7, 8,8,8,8,8,8,8,8,8, 9,9,9,10,10,10,11,11,11, 9,9,9,10,10,10,11,11,11, 9,9,9,10,10,10,11,11,11, 12,12,12,13,13,13,14,14,14, 12,12,12,13,13,13,14,14,14, 12,12,12,13,13,13,14,14,14), 14, 9, byrow=TRUE)) plot(c(0,1),c(0,1),type="n",axes=FALSE,xlab="",ylab="") mtext("Model does not include CD4",side=1,line=2,cex=1.2) plot(d1$sqrt.cd4^2,presids1,type="n",xlim=c(0,500),xlab="",ylab="",ylim=c(-1,1)) points(d1$sqrt.cd4[delta==0]^2,presids1[delta==0],pch=1,cex=.5,col=gray(.6)) points(d1$sqrt.cd4[delta==1]^2,presids1[delta==1],pch=3,cex=.5,col=gray(.6)) abline(h=0,lty=2) lines(supsmu(d1$sqrt.cd4^2,presids1,bass=1),col=1,lwd=4) mtext("CD4",side=1,line=2,cex=.7) mtext("PSR",side=2,line=2.5,cex=.7,padj=1) plot(d1$sqrt.cd4^2,psr1.norm,type="n",xlim=c(0,500),xlab="",ylab="") points(d1$sqrt.cd4[delta==0]^2,psr1.norm[delta==0],pch=1,cex=.5,col=gray(.6)) points(d1$sqrt.cd4[delta==1]^2,psr1.norm[delta==1],pch=3,cex=.5,col=gray(.6)) abline(h=0,lty=2) lines(supsmu(d1$sqrt.cd4^2,psr1.norm,bass=1),col=1,lwd=4) mtext("CD4",side=1,line=2,cex=.7) mtext("Normalized PSR",side=2,line=2.5,cex=.7,padj=1) fitpsr1<-survfit(Surv(presids1.cs,delta)~1) plot(summary(fitpsr1)$time,qunif(1-summary(fitpsr1)$surv,-1,1), xlab="", ylab="") abline(0,1,col=gray(.5)) mtext("Quantiles of Unif(-1,1)",side=1,line=2,cex=.7) mtext("Quantiles of Cox-Snell-like PSR",side=2,line=2.5,cex=.7,padj=1) plot(d1$sqrt.cd4^2,mresids1,type="n",xlim=c(0,500),xlab="",ylab="") points(d1$sqrt.cd4[delta==0]^2,mresids1[delta==0],pch=1,cex=.5,col=gray(.6)) points(d1$sqrt.cd4[delta==1]^2,mresids1[delta==1],pch=3,cex=.5,col=gray(.6)) abline(h=0,lty=2) lines(supsmu(d1$sqrt.cd4^2,mresids1,bass=1),col=1,lwd=4) mtext("CD4",side=1,line=2,cex=.7) mtext("Martingale Residual",side=2,line=2.5,cex=.7,padj=1) plot(d1$sqrt.cd4^2,dres1,type="n",xlim=c(0,500),xlab="",ylab="") points(d1$sqrt.cd4[delta==0]^2,dres1[delta==0],pch=1,cex=.5,col=gray(.6)) points(d1$sqrt.cd4[delta==1]^2,dres1[delta==1],pch=3,cex=.5,col=gray(.6)) abline(h=0,lty=2) lines(supsmu(d1$sqrt.cd4^2,dres1,bass=1),col=1,lwd=4) mtext("CD4",side=1,line=2,cex=.7) mtext("Deviance Residual",side=2,line=2.5,cex=.7,padj=1) plot(summary(fit1)$time,qexp(1-summary(fit1)$surv), xlab="",ylab="") mtext("Quantiles of Exp(1)",side=1,line=2,cex=.7) mtext("Quantiles of Cox-Snell Residual",side=2,line=2.5,cex=.7,padj=1) abline(0,1,col=gray(.5)) plot(c(0,1),c(0,1),type="n",axes=FALSE,xlab="",ylab="") mtext("Model includes CD4",side=1,line=2,cex=1.2) plot(d1$sqrt.cd4^2,presids,type="n",xlim=c(0,500),xlab="",ylab="",ylim=c(-1,1)) points(d1$sqrt.cd4[delta==0]^2,presids[delta==0],pch=1,cex=.5,col=gray(.6)) points(d1$sqrt.cd4[delta==1]^2,presids[delta==1],pch=3,cex=.5,col=gray(.6)) abline(h=0,lty=2) lines(supsmu(d1$sqrt.cd4^2,presids,bass=1),col=1,lwd=4) mtext("CD4",side=1,line=2,cex=.7) mtext("PSR",side=2,line=2.5,cex=.7,padj=1) plot(d1$sqrt.cd4^2,psr.norm,type="n",xlim=c(0,500),xlab="",ylab="") points(d1$sqrt.cd4[delta==0]^2,psr.norm[delta==0],pch=1,cex=.5,col=gray(.6)) points(d1$sqrt.cd4[delta==1]^2,psr.norm[delta==1],pch=3,cex=.5,col=gray(.6)) abline(h=0,lty=2) lines(supsmu(d1$sqrt.cd4^2,psr.norm,bass=1),col=1,lwd=4) mtext("CD4",side=1,line=2,cex=.7) mtext("Normalized PSR",side=2,line=2.5,cex=.7,padj=1) fitpsr<-survfit(Surv(presids.cs,delta)~1) plot(summary(fitpsr)$time,qunif(1-summary(fitpsr)$surv,-1,1), xlab="", ylab="") abline(0,1,col=gray(.5)) mtext("Quantiles of Unif(-1,1)",side=1,line=2,cex=.7) mtext("Quantiles of Cox-Snell-like PSR",side=2,line=2.5,cex=.7,padj=1) plot(d1$sqrt.cd4^2,mresids,type="n",xlim=c(0,500),xlab="",ylab="") points(d1$sqrt.cd4[delta==0]^2,mresids[delta==0],pch=1,cex=.5,col=gray(.6)) points(d1$sqrt.cd4[delta==1]^2,mresids[delta==1],pch=3,cex=.5,col=gray(.6)) abline(h=0,lty=2) lines(supsmu(d1$sqrt.cd4^2,mresids,bass=1),col=1,lwd=4) mtext("CD4",side=1,line=2,cex=.7) mtext("Martingale Residual",side=2,line=2.5,cex=.7,padj=1) plot(d1$sqrt.cd4^2,dres,type="n",xlim=c(0,500),xlab="",ylab="") points(d1$sqrt.cd4[delta==0]^2,dres[delta==0],pch=1,cex=.5,col=gray(.6)) points(d1$sqrt.cd4[delta==1]^2,dres[delta==1],pch=3,cex=.5,col=gray(.6)) abline(h=0,lty=2) lines(supsmu(d1$sqrt.cd4^2,dres,bass=1),col=1,lwd=4) mtext("CD4",side=1,line=2,cex=.7) mtext("Deviance Residual",side=2,line=2.5,cex=.7,padj=1) plot(summary(fit)$time,qexp(1-summary(fit)$surv), xlab="",ylab="") abline(0,1,col=gray(.5)) mtext("Quantiles of Exp(1)",side=1,line=2,cex=.7) mtext("Quantiles of Cox-Snell Residual",side=2,line=2.5,cex=.7,padj=1) dev.off() ###################################