#### This file contains the code used for all the analyses reported in Shepherd, Redman, Ankerst (2008). #### Data is not provided. However, the unadjusted analyses (Figures 2 and 3) are completely reproduced with this code. ###### GBH-type analysis: assuming missing S are completely at random, ###### assuming those without biopsy do not have severe cancer, and ###### assuming monotonicity. N<-10168 Nf<-4951 Np<-5217 NS1f<-821 NS1p<-1194 NS1Y1f<-299 NS1Y1p<-264 z<-c(rep(1,Nf),rep(0,Np)) s<-c(rep(1,NS1f),rep(0,Nf-NS1f),rep(1,NS1p),rep(0,Np-NS1p)) y<-c(rep(1,NS1Y1f),rep(0,NS1f-NS1Y1f),rep(NA,Nf-NS1f),rep(1,NS1Y1p),rep(0,NS1p-NS1Y1p),rep(NA,Np-NS1p)) m1<-mean(y[s*z==1]) ey0<-mean(y[s*(1-z)==1]) TE<-1-(sum(s*z)/sum(z))/(sum(s*(1-z))/sum(1-z)) ### Naive analysis ace.naive<-mean(y[s*z==1])-mean(y[s*(1-z)==1]) se.naive<-sqrt(mean(y[s*z==1])*(1-mean(y[s*z==1]))/sum(s*z) + mean(y[s*(1-z)==1])*(1-mean(y[s*(1-z)==1]))/sum(s*(1-z))) ace.naive ace.naive-1.96*se.naive ace.naive+1.96*se.naive ### GBH Analysis ee<-function (a){ (sum(exp(a+beta*y[s*(1-z)==1])/(1+exp(a+beta*y[s*(1-z)==1])))/sum(s*(1-z)) - (1-TE))^2 } p0hat<-sum(s*(1-z))/sum(1-z) p1hat<-sum(s*z)/sum(z) alphahat<-m0hat<-ace<-varace<-NULL b<-c(-12:12)/2.4 for (j in 1:length(b)) { beta<-b[j] min<-optimize(ee, c(-50,50)) alphahat[j]<-min$minimum m0hat[j]<-(sum(exp(alphahat[j]+beta*y[s*(1-z)==1])/(1+exp(alphahat[j]+beta*y[s*(1-z)==1]))*y[s*(1-z)==1])/sum(s*(1-z)))/ifelse(TE<0,1,1-TE) m1hat<-mean(y[s*z==1]) ace[j]<-m1hat-m0hat[j] ############# Calculating GBH Asymptotic 95% CI U1<-(1-z)*(p0hat-s) U2<-z*(p1hat-s) U3<-(1-z)*s*((1+exp(-alphahat[j]-beta*y))^(-1)-p1hat/p0hat) U4<-(1-z)*s*(m0hat[j]-y*(1+exp(-alphahat[j]-beta*y))^(-1)*p0hat/p1hat) U5<-z*s*(m1hat-y) om11<-sum(U1*U1,na.rm=T)/length(U1) om12<-sum(U1*U2,na.rm=T)/length(U1) om13<-sum(U1*U3,na.rm=T)/length(U1) om14<-sum(U1*U4,na.rm=T)/length(U1) om15<-sum(U1*U5,na.rm=T)/length(U1) om22<-sum(U2*U2,na.rm=T)/length(U1) om23<-sum(U2*U3,na.rm=T)/length(U1) om24<-sum(U2*U4,na.rm=T)/length(U1) om25<-sum(U2*U5,na.rm=T)/length(U1) om33<-sum(U3*U3,na.rm=T)/length(U1) om34<-sum(U3*U4,na.rm=T)/length(U1) om35<-sum(U3*U5,na.rm=T)/length(U1) om44<-sum(U4*U4,na.rm=T)/length(U1) om45<-sum(U4*U5,na.rm=T)/length(U1) om55<-sum(U5*U5,na.rm=T)/length(U1) Omega<-matrix(c(om11,om12,om13,om14,om15, om12,om22,om23,om24,om25, om13,om23,om33,om34,om35, om14,om24,om34,om44,om45, om15,om25,om35,om45,om55),5,5) Gamma<-matrix(0,ncol=5,nrow=5) Gamma[1,1]<-sum(1-z)/N Gamma[2,2]<-sum(z)/N Gamma[3,1]<-sum((1-z)*s*p1hat/p0hat^2)/N Gamma[3,2]<-sum(-(1-z)*s/p0hat)/N Gamma[3,3]<-sum((1-z)*s*exp(alphahat[j]+beta*y)/(1+exp(alphahat[j]+beta*y))^2,na.rm=T)/N Gamma[4,1]<-sum(-(1-z)*s*y*(1+exp(-alphahat[j]-beta*y))^(-1)/(p1hat),na.rm=T)/N Gamma[4,2]<-sum((1-z)*s*y*(1+exp(-alphahat[j]-beta*y))^(-1)*p0hat/(p1hat^2),na.rm=T)/N Gamma[4,3]<-sum(-(1-z)*s*y*exp(alphahat[j]+beta*y)/(1+exp(alphahat[j]+beta*y))^2*p0hat/(p1hat),na.rm=T)/N Gamma[4,4]<-sum((1-z)*s)/N Gamma[5,5]<-sum(z*s)/N varpars<-solve(Gamma)%*%Omega%*%t(solve(Gamma))/N varace[j]<-varpars[5,5]+varpars[4,4]-2*varpars[4,5] } lowerw<-ace-1.96*sqrt(varace) upperw<-ace+1.96*sqrt(varace) ####### Worst-case analysis of HHS (when beta = -/+ infinity) ### Creating artificial ranking of y. This makes computation easier when there are ties, but in no way affects results adj<-c(1:length(y))/1000000 yalt<-y+adj ace.hhs.up<-m1-mean(y[s*(1-z)==1&yalt=quantile(yalt[s*(1-z)==1],TE)]) yp<-y[z==0] yv<-y[z==1] sp<-s[z==0] sv<-s[z==1] Nv<-length(yv) Np<-length(yp) VEboot<-ace.hhs.up.boot<-ace.hhs.low.boot<-NULL for (ii in 1:1000) { bootsampp<-sample(c(1:Np),Np,replace=T) spboot<-sp[bootsampp] ypboot<-yp[bootsampp] bootsampv<-sample(c(1:Nv),Nv,replace=T) svboot<-sv[bootsampv] yvboot<-yv[bootsampv] VEboot[ii]<- 1-(sum(svboot)/Nv)/(sum(spboot)/Np) adjboot<-c(1:length(yp))/1000000 ypalt<-ypboot+adjboot m1.boot<-mean(yvboot[svboot==1]) ace.hhs.up.boot[ii]<-m1.boot-mean(ypboot[spboot==1&ypalt=quantile(ypalt[spboot==1],VEboot[ii])]) } hhsbootpup<-quantile(ace.hhs.up.boot,c(.025,.975)) hhsbootplow<-quantile(ace.hhs.low.boot,c(.025,.975)) hhsbootwup<-c(ace.hhs.up-1.96*sqrt(var(ace.hhs.up.boot)),ace.hhs.up+1.96*sqrt(var(ace.hhs.up.boot))) hhsbootwlow<-c(ace.hhs.low-1.96*sqrt(var(ace.hhs.low.boot)),ace.hhs.low+1.96*sqrt(var(ace.hhs.low.boot))) ## This plot is Figure 2 of Shepherd, Redman, and Ankerst (2008) #postscript("finasteride2a.eps",horiz=T) par(mfrow=c(1,1),lab=c(10,5,1),mar=c(6,4,4,2)) plot(c(b,b,0,0,0),c(upperw,lowerw,0,hhsbootpup[2],hhsbootplow[1]),xlab="",ylab="ACE",type="n") lines(b,ace) lines(b,upperw,lty=2) lines(b,lowerw,lty=2) points(c(-5.25,5.25),c(ace.hhs.up,ace.hhs.low)) points(c(-5.25,5.25),c(ace.hhs.up,ace.hhs.low)) points(c(-5.25,-5.25),hhsbootpup,pch=3) points(c(5.25,5.25),hhsbootplow,pch=3) abline(0,0,lty=3) axis(1,at=c(-5:5),labels=c(0.01,0.02,0.05,0.14,0.37,1,2.7,7.4,20,55,148),line=1.5,tick=FALSE) mtext(expression(beta[0]),side=1,at=-6,line=1) mtext("OR",side=1,at=-6,line=2.5) legend(1.5,.35,c("Estimated ACE", "95% C.I."),lty=c(1,2),bty="n") legend(1.8,.32,legend=c(expression("ACE at " * beta[0] == +-infinity), expression("95% C.I. at " * beta[0] ==+- infinity)),pch=c(1,3),bty="n") #dev.off() ##################################################################################################### ############################ Analysis Relaxing monotonicity assuming MAR ##################################################################################################### rm(list=ls()) N<-10168 Nf<-4951 Np<-5217 NS1f<-821 NS1p<-1194 NS1Y1f<-299 NS1Y1p<-264 # The numbers given below are from an earlier dataset. Unfortunately, Figure 3 in Shepherd, Redman, and Ankerst (2008) used these numbers # instead of the more up-to-date numbers given above (the numbers that were used in all other figures). Results are very similar using both sets of numbers. # However, to identically reproduce Figure 3 in the published paper, the numbers given below need to be used. Unfortunately, some of the difference between # Figures 3 and 4 in the published paper are actually due to differences in the numbers used -- not necessarily differences due to adjusting for covariates. # N<-9060 # Nf<-4368 # Np<-4692 # NS1f<-803 # NS1p<-1147 # NS1Y1f<-280 # NS1Y1p<-237 z<-c(rep(1,Nf),rep(0,Np)) s<-c(rep(1,NS1f),rep(0,Nf-NS1f),rep(1,NS1p),rep(0,Np-NS1p)) y<-c(rep(1,NS1Y1f),rep(0,NS1f-NS1Y1f),rep(NA,Nf-NS1f),rep(1,NS1Y1p),rep(0,NS1p-NS1Y1p),rep(NA,Np-NS1p)) ey1<-mean(y[s*z==1]) ey0<-mean(y[s*(1-z)==1]) p0hat<-(sum(s*(1-z))/sum(1-z)) p1hat<-sum(s*z)/sum(z) phi<- .1 beta0<-0 beta1<-0 sensanalysis<-function(phi,beta0,beta1){ pi<-(1-phi)*p1hat eealpha0<-function (a) { (sum((1-z)*s*((1+exp(-a-beta0*y))^(-1)-pi/p0hat),na.rm=T))^2 } eealpha1<-function (a) { (sum(z*s*((1+exp(-a-beta1*y))^(-1)-pi/p1hat),na.rm=T))^2 } min<-optimize(eealpha0, c(-50,50)) alpha0hat<-min$minimum min<-optimize(eealpha1, c(-50,50)) alpha1hat<-min$minimum eem0<-function (a) { (sum((1-z)*s*(a-y*(1+exp(-alpha0hat-beta0*y))^(-1)*p0hat/pi),na.rm=T))^2 } min<-optimize(eem0, c(0,1)) m0hat<-min$minimum eem1<-function (a) { (sum(z*s*(a-y*(1+exp(-alpha1hat-beta1*y))^(-1)*p1hat/pi),na.rm=T))^2 } min<-optimize(eem1, c(0,1)) m1hat<-min$minimum U1<-(1-z)*(p0hat-s) U2<-z*(p1hat-s) U3<-(1-z)*s*((1+exp(-alpha0hat-beta0*y))^(-1)-pi/p0hat) U4<-z*s*((1+exp(-alpha1hat-beta1*y))^(-1)-pi/p1hat) U5<-(1-z)*s*(m0hat-y*(1+exp(-alpha0hat-beta0*y))^(-1)*p0hat/pi) U6<-z*s*(m1hat-y*(1+exp(-alpha1hat-beta1*y))^(-1)*p1hat/pi) om11<-sum(U1*U1,na.rm=T)/length(U1) om12<-sum(U1*U2,na.rm=T)/length(U1) om13<-sum(U1*U3,na.rm=T)/length(U1) om14<-sum(U1*U4,na.rm=T)/length(U1) om15<-sum(U1*U5,na.rm=T)/length(U1) om16<-sum(U1*U6,na.rm=T)/length(U1) om22<-sum(U2*U2,na.rm=T)/length(U1) om23<-sum(U2*U3,na.rm=T)/length(U1) om24<-sum(U2*U4,na.rm=T)/length(U1) om25<-sum(U2*U5,na.rm=T)/length(U1) om26<-sum(U2*U6,na.rm=T)/length(U1) om33<-sum(U3*U3,na.rm=T)/length(U1) om34<-sum(U3*U4,na.rm=T)/length(U1) om35<-sum(U3*U5,na.rm=T)/length(U1) om36<-sum(U3*U6,na.rm=T)/length(U1) om44<-sum(U4*U4,na.rm=T)/length(U1) om45<-sum(U4*U5,na.rm=T)/length(U1) om46<-sum(U4*U4,na.rm=T)/length(U1) om55<-sum(U5*U5,na.rm=T)/length(U1) om56<-sum(U5*U6,na.rm=T)/length(U1) om66<-sum(U6*U6,na.rm=T)/length(U1) Omega<-matrix(c(om11,om12,om13,om14,om15,om16,om12,om22,om23,om24,om25,om26, om13,om23,om33,om34,om35,om36,om14,om24,om34,om44,om45,om46, om15,om25,om35,om45,om55,om56,om16,om26,om36,om46,om56,om66),6,6) Gamma<-matrix(0,ncol=6,nrow=6) Gamma[1,1]<-sum(1-z)/N Gamma[2,2]<-sum(z)/N Gamma[3,1]<-sum((1-z)*s*(1-phi)*p1hat/p0hat^2)/N Gamma[3,2]<-sum(-(1-z)*s*(1-phi)/p0hat)/N Gamma[3,3]<-sum((1-z)*s*exp(alpha0hat+beta0*y)/(1+exp(alpha0hat+beta0*y))^2,na.rm=T)/N Gamma[4,4]<-sum(z*s*exp(alpha1hat+beta1*y)/(1+exp(alpha1hat+beta1*y))^2,na.rm=T)/N Gamma[5,1]<-sum(-(1-z)*s*y*(1+exp(-alpha0hat-beta0*y))^(-1)/(p1hat*(1-phi)),na.rm=T)/N Gamma[5,2]<-sum((1-z)*s*y*(1+exp(-alpha0hat-beta0*y))^(-1)*p0hat/(p1hat^2*(1-phi)),na.rm=T)/N Gamma[5,3]<-sum(-(1-z)*s*y*exp(alpha0hat+beta0*y)/(1+exp(alpha0hat+beta0*y))^2*p0hat/(p1hat*(1-phi)),na.rm=T)/N Gamma[5,5]<-sum((1-z)*s)/N Gamma[6,4]<-sum(-z*s*y*exp(alpha1hat+beta1*y)/(1+exp(alpha1hat+beta1*y))^2/(1-phi),na.rm=T)/N Gamma[6,6]<-sum(z*s)/N varpars<-solve(Gamma)%*%Omega%*%t(solve(Gamma))/N ace<-m1hat-m0hat varace<-varpars[6,6]+varpars[5,5]-2*varpars[5,6] return(list(pi = pi, ace = ace, varace = varace)) } b0 = -20:20/8 b1 = -20:20/8 #### This process is very time consuming. For this reason, I save the results so that I only had to do it once for all. date() phi<-.01 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } date() save(ace,vara,file="relaxmons-01a.Rdata") phi<-.05 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-05a.Rdata") phi<-.1 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-10a.Rdata") phi<-.2 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-20a.Rdata") phi<-.5 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-50a.Rdata") phi<-1-p0hat ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-inda.Rdata") #### This plot is Figure 3 of Shepherd, Redman, Ankerst (2008) #postscript("relaxmonsa1.eps",horiz=F) par(mfrow=c(3,2),mar=c(4,4,2,3)) b0<-b1<--20:20/8 load("relaxmons-01a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=gray(.8),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=3, line=0, expression(phi==0.99)) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(1.05),log(1.05),log(1.35),log(1.35),log(1.05)),c(log(1.05),log(1.35),log(1.35),log(1.05),log(1.05))) ## Thompson points(log(1.2),log(1.2),pch=4) load("relaxmons-05a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.95)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(2),log(2),log(4),log(4),log(2)),c(log(.25),log(.5),log(.5),log(.25),log(.25))) ### Dr Kristal lines(c(log(1.05),log(1.05),log(1.35),log(1.35),log(1.05)),c(log(1.05),log(1.35),log(1.35),log(1.05),log(1.05))) ## Thompson load("relaxmons-10a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.90)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(2),log(2),log(4),log(4),log(2)),c(log(.25),log(.5),log(.5),log(.25),log(.25))) ### Dr Kristal points(log(3),log(1/3),pch=4) load("relaxmons-20a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.80)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(2),log(2),log(4),log(4),log(2)),c(log(.25),log(.5),log(.5),log(.25),log(.25))) ### Dr Kristal load("relaxmons-50a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.50)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() load("relaxmons-inda.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, quote("S(0),S(1) Independent, " * phi==0.24)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() #dev.off() ###################################################################################### ############## Now assuming MAR instead of MCAR ###################################################################################### rm(list=ls()) #### Reading in the PCPT data from August 2006 pcpt =read.delim("pcpt.txt",header = T,sep = '\t') pcpt$Z = 2 - pcpt$TRNO ##### Reading data from an updated file received in January 2007 ##### The only information I am going to take from this new dataset ##### is with regards to whether they had a prior negative biopsy library(foreign) datastuff<-read.dta("final_20dec06.dta") ##### Merging the two together data2merge<-data.frame(datastuff$patno,datastuff$caind,datastuff$lastbiopdt,datastuff$priorbiopdt) new.pcpt<-merge(pcpt,data2merge,by.x="PATNO",by.y="datastuff.patno",all=TRUE) #### Note: One id was in the pcpt dataframe but not the datastuff #### dataframe. In the lines below, this person is assigned no prior negative biopsy. attach(new.pcpt) pnegbx<-ifelse((!is.na(datastuff.lastbiopdt)&hasbx==0)|!is.na(datastuff.priorbiopdt),1,0) levels(cause) = c('Both','DRE','PSA','Neither','Other','Other Proc') # create indicator of high grade y1<-ifelse(PRVCAGLT >=95,NA,PRVCAGLT) y = 1*(PRVCAGLT >= 7) # if y is missing when s=1, then y=0 y<-ifelse(pca==1&is.na(y),0,y) # id,trno, hasbx, cancer, hg, weights new.data = as.data.frame(cbind(PATNO, Z, hasbx,pca,y)) names(new.data) = c('id', 'z','r','s','y') attach(new.data) N<-length(z) agesq<-age^2 psa0sq<-psa0^2 modlambda.1 <- glm(has7yr ~ Z + psa0 + age + agesq + white + aa + fampca, family = binomial, data = new.pcpt) #### Just looking at the impact of age on having a 7 year visit agejunk<-c(51:86) or.age<-exp(modlambda.1$coeff[4]*(agejunk-mean(age))+modlambda.1$coeff[5]*(agejunk^2-mean(age)^2)) ## The X matrix contains a column of 1s, z, and the covariates ## included in the logistic regression model; a quadratic for age was included X.1 <-cbind(1,Z,psa0,age,agesq,white,aa,fampca) dimX.1 <-ncol(X.1) ## lambda is predicted probability of missing lambda.1<-modlambda.1$fitted ## It's worth looking at lambda. If some lambda's are too close to 0 ## or 1 we may have some problems. We don't seem to have any problems summary(lambda.1) ## These next two steps are used to compute Omega score.1 = (has7yr - lambda.1) * X.1 Omega.1<-matrix(NA,dimX.1,dimX.1) for (i in 1:dimX.1) { for (j in 1:dimX.1) { Omega.1[i,j]<-sum(score.1[,i]*score.1[,j])/N } } ## Second part of weights N.2 = sum(has7yr) log.psalast<-log(psalast) summary(psalast[has7yr==1]) missing.psalast<-id[has7yr==1&is.na(psalast)] modlambda.2 = glm(hasbx ~ Z + psa0 + age + agesq + white + aa + fampca + cause_psa + cause_dre + pnegbx + Z:cause_psa +Z:cause_dre, family = binomial, data = new.pcpt,subset = has7yr == 1) X.2 <-cbind(1,Z,psa0,age,agesq,white,aa,fampca,cause_psa,cause_dre,pnegbx,Z*cause_psa,Z*cause_dre)*has7yr dimX.2 <-ncol(X.2) summary(X.2) ## eta.x is the linear predictor; lambda is predicted probability of missing eta.x.2<-t(modlambda.2$coeff %*% t(X.2)) lambda.2<-exp(eta.x.2)/(1+exp(eta.x.2)) ## It's worth looking at lambda. If some lambda's are too close to 0 ## or 1 we may have some problems. summary(lambda.2) ## These next two loops are used to compute Omega score.2<-matrix(NA,N,dimX.2) r.2 = r[has7yr == 1] for (i in 1:N) { score.2[i,]<-(r[i]-lambda.2[i])*X.2[i,] } Omega.2<-matrix(NA,dimX.2,dimX.2) for (i in 1:dimX.2) { for (j in 1:dimX.2) { Omega.2[i,j]<-sum(score.2[,i]*score.2[,j], na.rm = T)/sum(!is.na(lambda.2)) } } # combine two matrixes into one Omega score = cbind(score.1,score.2) Omega = matrix(0,dimX.1 + dimX.2,dimX.1 + dimX.2) Omega[1:dimX.1,1:dimX.1] = solve(Omega.1) Omega[(dimX.1 + 1):(dimX.1+dimX.2),(dimX.1 + 1):(dimX.1+dimX.2)] = solve(Omega.2) ## The MCAR estimating equations are weighted by w w<-ifelse(hasbx == 0, 0, 1/(lambda.1*lambda.2)) ## Estimating P(S=1|Z=0) and P(S=1|Z=1) p0hat<-sum(s*(1-z)*w,na.rm=T)/sum((1-z)*w,na.rm=T) p1hat<-sum(s*z*w,na.rm=T)/sum(z*w,na.rm=T) ## Now I will need to estimate at fixed values of phi, beta0, and ## beta1. Since these parameters are unknown, I perform a sensitivity ## analysis over a range of values. ## I define a function 'sensanalysis' which performs a sensitivity analysis at an ## assumed value of phi (varies beta0 and beta1) sensanalysis<-function(phi, beta0, beta1) { ## Re-parameterizing pi<-(1-phi)*p1hat ## Estimating alpha0 and alpha1 eealpha0<-function (a) { (sum(w*(1-z)*s*((1+exp(-a-beta0*y))^(-1)-pi/p0hat),na.rm=T))^2 } eealpha1<-function (a) { (sum(w*z*s*((1+exp(-a-beta1*y))^(-1)-pi/p1hat),na.rm=T))^2 } min<-optimize(eealpha0, c(-50,50)) alpha0hat<-min$minimum min<-optimize(eealpha1, c(-50,50)) alpha1hat<-min$minimum eem0<-function (a) { (sum(w*(1-z)*s*(a-y*(1+exp(-alpha0hat-beta0*y))^(-1)*p0hat/pi),na.rm=T))^2 } min<-optimize(eem0, c(0,1)) m0hat<-min$minimum eem1<-function (a) { (sum(w*z*s*(a-y*(1+exp(-alpha1hat-beta1*y))^(-1)*p1hat/pi),na.rm=T))^2 } min<-optimize(eem1, c(0,1)) m1hat<-min$minimum ## Estimating equations (these should sum to 0) V1<-w*(1-z)*(p0hat-s) V2<-w*z*(p1hat-s) V3<-w*(1-z)*s*((1+exp(-alpha0hat-beta0*y))^(-1)-pi/p0hat) V4<-w*z*s*((1+exp(-alpha1hat-beta1*y))^(-1)-pi/p1hat) V5<-w*(1-z)*s*(m0hat-y*(1+exp(-alpha0hat-beta0*y))^(-1)*p0hat/pi) V6<-w*z*s*(m1hat-y*(1+exp(-alpha1hat-beta1*y))^(-1)*p1hat/pi) ## Estimating variance. Gamma = matrix(0,ncol = 6, nrow = 6) Gamma[1,1] = sum(w*(1-z),na.rm=T)/N Gamma[2,2] = sum(w*z,na.rm=T)/N Gamma[3,1] = sum(w*(1-z)*s*(1-phi)*p1hat/p0hat^2,na.rm=T)/N Gamma[3,2] = sum(-w*(1-z)*s*(1-phi)/p0hat,na.rm=T)/N Gamma[3,3] = sum(w*(1-z)*s*exp(alpha0hat+beta0*y)/(1+exp(alpha0hat+beta0*y))^2,na.rm=T)/N Gamma[4,4] = sum(w*z*s*exp(alpha1hat+beta1*y)/(1+exp(alpha1hat+beta1*y))^2,na.rm=T)/N Gamma[5,1] = sum(-w*(1-z)*s*y*(1+exp(-alpha0hat-beta0*y))^(-1)/(p1hat*(1-phi)),na.rm=T)/N Gamma[5,2] = sum(w*(1-z)*s*y*(1+exp(-alpha0hat-beta0*y))^(-1)*p0hat/(p1hat^2*(1-phi)),na.rm=T)/N Gamma[5,3] = sum(-w*(1-z)*s*y*exp(alpha0hat+beta0*y)/(1+exp(alpha0hat+beta0*y))^2*p0hat/(p1hat*(1-phi)),na.rm=T)/N Gamma[5,5] = sum(w*(1-z)*s,na.rm=T)/N Gamma[6,4] = sum(-w*z*s*y*exp(alpha1hat+beta1*y)/(1+exp(alpha1hat+beta1*y))^2/(1-phi),na.rm=T)/N Gamma[6,6] = sum(w*z*s,na.rm=T)/N # G.gamma (B matrix in RRZ) b1<-b2<-b3<-b4<-b5<-b6<-NULL for (jj in 1:dimX.1) { b1[jj]<-sum(V1*(1-lambda.1)*X.1[,jj],na.rm=T) b2[jj]<-sum(V2*(1-lambda.1)*X.1[,jj],na.rm=T) b3[jj]<-sum(V3*(1-lambda.1)*X.1[,jj],na.rm=T) b4[jj]<-sum(V4*(1-lambda.1)*X.1[,jj],na.rm=T) b5[jj]<-sum(V5*(1-lambda.1)*X.1[,jj],na.rm=T) b6[jj]<-sum(V6*(1-lambda.1)*X.1[,jj],na.rm=T) } for (jj in (dimX.1+1):(dimX.1 + dimX.2)) { b1[jj]<-sum(V1*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) b2[jj]<-sum(V2*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) b3[jj]<-sum(V3*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) b4[jj]<-sum(V4*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) b5[jj]<-sum(V5*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) b6[jj]<-sum(V6*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) } G.gamma = rbind(b1,b2,b3,b4,b5,b6)/N # MM matrix temp.MM = matrix(0, nrow = N, ncol = (dimX.1 + dimX.2)^2) for(m in 1:N) { temp.MM[m,] = as.vector(score[m,] %o% score[m,]) } MM = matrix(colSums(temp.MM,na.rm = T), nrow = (dimX.1 + dimX.2)) /N # PSI psi.matrix = - solve(MM) %*% t(score) # adjusted estimating equation g.tilda = rbind(V1,V2,V3,V4,V5,V6) + G.gamma %*% psi.matrix temp.gg = matrix(0, nrow = N, ncol = (6)^2) for(m in 1:N) { temp.gg[m,] = as.vector(g.tilda[,m] %o% g.tilda[,m]) } # adjusted variance of estimating equations GG = matrix(colSums(temp.gg,na.rm = T), nrow = (6)) / N varpars.new = (solve(Gamma) %*% GG %*% t(solve(Gamma)))/N # estimate of and variance of ACE ace <-m1hat-m0hat varace<-varpars.new[6,6]+varpars.new[5,5]-2*varpars.new[5,6] return(list(pi = pi, ace = ace, varace = varace)) } b0 = -20:20/8 b1 = -20:20/8 phi<-.01 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) date() for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } date() save(ace,vara,file="relaxmons-mar-01a.Rdata") phi<-.05 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-mar-05a.Rdata") phi<-.1 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-mar-10a.Rdata") phi<-.2 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-mar-20a.Rdata") phi<-.5 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-mar-50a.Rdata") phi<-1-p0hat ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-mar-inda.Rdata") ########### This code creates Figure 4 of Shepherd, Redman, and Ankerst (2008) #postscript("relaxmons-mara.eps",horiz=F) par(mfrow=c(3,2),mar=c(4,4,2,3)) b0<-b1<--20:20/8 load("relaxmons-mar-01a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=gray(.8),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=3, line=0, expression(phi==0.99)) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(1.05),log(1.05),log(1.35),log(1.35),log(1.05)),c(log(1.05),log(1.35),log(1.35),log(1.05),log(1.05))) ## Thompson points(log(1.2),log(1.2),pch=4) load("relaxmons-mar-05a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.95)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(2),log(2),log(4),log(4),log(2)),c(log(.25),log(.5),log(.5),log(.25),log(.25))) ### Dr Kristal lines(c(log(1.05),log(1.05),log(1.35),log(1.35),log(1.05)),c(log(1.05),log(1.35),log(1.35),log(1.05),log(1.05))) ## Thompson load("relaxmons-mar-10a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.90)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(2),log(2),log(4),log(4),log(2)),c(log(.25),log(.5),log(.5),log(.25),log(.25))) ### Dr Kristal points(log(3),log(1/3),pch=4) load("relaxmons-mar-20a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.80)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(2),log(2),log(4),log(4),log(2)),c(log(.25),log(.5),log(.5),log(.25),log(.25))) ### Dr Kristal load("relaxmons-mar-50a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.50)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() load("relaxmons-mar-inda.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, quote("S(0),S(1) Independent, " * phi==0.24)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() #dev.off() #################################################################################################### ###### Analyses used to create Tables 3, 4, and 5 and Figure 5 using prostatectomy results #################################################################################################### rm(list=ls()) pcpt =read.delim("pcpt.txt",header = T,sep = '\t') pcpt$Z = 2 - pcpt$TRNO names(pcpt) library(foreign) datastuff<-read.dta("final_20dec06.dta") summary(datastuff) ##### Merging the two together data2merge<-data.frame(datastuff$patno,datastuff$caind,datastuff$lastbiopdt,datastuff$priorbiopdt) new.pcpt<-merge(pcpt,data2merge,by.x="PATNO",by.y="datastuff.patno",all=TRUE) attach(new.pcpt) pnegbx<-ifelse((!is.na(datastuff.lastbiopdt)&hasbx==0)|!is.na(datastuff.priorbiopdt),1,0) levels(cause) = c('Both','DRE','PSA','Neither','Other','Other Proc') # create indicator of high grade y1<-ifelse(PRVCAGLT >=95,NA,PRVCAGLT) y = 1*(PRVCAGLT >= 7) # if y is missing when s=1, then y=0 y<-ifelse(pca==1&is.na(y),0,y) new.data = as.data.frame(cbind(PATNO, Z, hasbx,pca,y)) names(new.data) = c('id', 'z','r','s','y') sum(pcpt$V) attach(new.data) N<-length(z) agesq<-age^2 psa0sq<-psa0^2 modlambda.1 <- glm(has7yr ~ Z + psa0 + age + agesq + white + aa + fampca, family = binomial, data = new.pcpt) #### Creating tables table(Z,hasbx) 3015/(2808+3015) 4951/(5217+4951) chisq.test(table(Z,hasbx)) summary(psa0[hasbx==0]) summary(psa0[hasbx==1]) wilcox.test(psa0[hasbx==0],psa0[hasbx==1]) summary(age[hasbx==0]) summary(age[hasbx==1]) wilcox.test(age[hasbx==0],age[hasbx==1]) race<-ifelse(white==1,"white",ifelse(aa==1,"aa","other")) table(race,hasbx) c(242,271,5310)/sum(c(242,271,5310)) c(346,351,9471)/sum(c(346,351,9471)) chisq.test(table(race,hasbx)) table(fampca,hasbx) c(5038,785)/sum(c(5038,785)) c(8473,1695)/sum(c(8473,1695)) chisq.test(table(fampca,hasbx)) table(cause_psa,hasbx) 69/(5754+69) 802/(9366+802) chisq.test(table(cause_psa,hasbx)) table(cause_dre,hasbx) 83/(5740+83) 831/(9337+831) chisq.test(table(cause_dre,hasbx)) table(pnegbx,hasbx) 477/(477+5346) 1348/(1348+8820) chisq.test(table(pnegbx,hasbx)) table(Z[pca==1],V[pca==1]) 596/1484 225/531 chisq.test(table(Z[pca==1],V[pca==1])) summary(psa0[pca==1&V==0]) summary(psa0[pca==1&V==1]) wilcox.test(psa0[pca==1&V==0],psa0[pca==1&V==1]) summary(age[pca==1&V==0]) summary(age[pca==1&V==1]) wilcox.test(age[pca==1&V==0],age[pca==1&V==1]) race<-ifelse(white==1,"white",ifelse(aa==1,"aa","other")) table(race[pca==1],V[pca==1]) table(race[pca==1],V[pca==1])/531 chisq.test(table(race[pca==1],V[pca==1])) table(fampca[pca==1],V[pca==1]) chisq.test(table(fampca[pca==1],V[pca==1])) table(cause_psa[pca==1],V[pca==1]) chisq.test(table(cause_psa[pca==1],V[pca==1])) table(cause_dre[pca==1],V[pca==1]) chisq.test(table(cause_dre[pca==1],V[pca==1])) table(pnegbx[pca==1],V[pca==1]) chisq.test(table(pnegbx[pca==1],V[pca==1])) table(y[pca==1],V[pca==1]) chisq.test(table(y[pca==1],V[pca==1])) ###### Looking at association between covariates and pca exp(glm(pca ~ Z, family=binomial, subset = hasbx==1)$coeff[2]) exp(glm(pca ~ psa0, family=binomial, subset = hasbx==1)$coeff[2]) exp(glm(pca ~ age, family=binomial, subset = hasbx==1)$coeff[2]) exp(glm(pca ~ white+aa, family=binomial, subset = hasbx==1)$coeff[2]) exp(glm(pca ~ white+aa, family=binomial, subset = hasbx==1)$coeff[3]) exp(glm(pca ~ fampca, family=binomial, subset = hasbx==1)$coeff[2]) exp(glm(pca ~ cause_psa, family=binomial, subset = hasbx==1)$coeff[2]) exp(glm(pca ~ cause_dre, family=binomial, subset = hasbx==1)$coeff[2]) exp(glm(pca ~ pnegbx, family=binomial, subset = hasbx==1)$coeff[2]) glm(pca ~ cause_psa, family=binomial, subset = hasbx==1)$coeff[2] glm(pca ~ psa0, family=binomial, subset = hasbx==1&Z==0)$coeff[2] glm(pca ~ age, family=binomial, subset = hasbx==1&Z==0)$coeff[2] glm(pca ~ white+aa, family=binomial, subset = hasbx==1&Z==0)$coeff[2] glm(pca ~ fampca, family=binomial, subset = hasbx==1&Z==0)$coeff[2] glm(pca ~ cause_psa, family=binomial, subset = hasbx==1&Z==0)$coeff[2] glm(pca ~ cause_dre, family=binomial, subset = hasbx==1&Z==0)$coeff[2] glm(pca ~ pnegbx, family=binomial, subset = hasbx==1&Z==0)$coeff[2] exp(glm(pca ~ cause_psa, family=binomial, subset = hasbx==1&Z==0)$coeff[2]) glm(pca ~ psa0, family=binomial, subset = hasbx==1&Z==1)$coeff[2] glm(pca ~ age, family=binomial, subset = hasbx==1&Z==1)$coeff[2] glm(pca ~ white+aa, family=binomial, subset = hasbx==1&Z==1)$coeff[2] glm(pca ~ fampca, family=binomial, subset = hasbx==1&Z==1)$coeff[2] glm(pca ~ cause_psa, family=binomial, subset = hasbx==1&Z==1)$coeff[2] glm(pca ~ cause_dre, family=binomial, subset = hasbx==1&Z==1)$coeff[2] glm(pca ~ pnegbx, family=binomial, subset = hasbx==1&Z==1)$coeff[2] exp(glm(pca ~ cause_psa, family=binomial, subset = hasbx==1&Z==1)$coeff[2]) X.1 <-cbind(1,Z,psa0,age,agesq,white,aa,fampca) dimX.1 <-ncol(X.1) lambda.1<-modlambda.1$fitted summary(lambda.1) score.1 = (has7yr - lambda.1) * X.1 Omega.1<-matrix(NA,dimX.1,dimX.1) for (i in 1:dimX.1) { for (j in 1:dimX.1) { Omega.1[i,j]<-sum(score.1[,i]*score.1[,j])/N } } ## Second part of weights N.2 = sum(has7yr) log.psalast<-log(psalast) summary(psalast[has7yr==1]) missing.psalast<-id[has7yr==1&is.na(psalast)] ### models with age^3 and log.psalast^2 didn't seem to add anything modlambda.2 = glm(hasbx ~ Z + psa0 + age + agesq + white + aa + fampca + cause_psa + cause_dre + pnegbx + Z:cause_psa +Z:cause_dre, family = binomial, data = new.pcpt,subset = has7yr == 1) X.2 <-cbind(1,Z,psa0,age,agesq,white,aa,fampca,cause_psa,cause_dre,pnegbx,Z*cause_psa,Z*cause_dre)*has7yr dimX.2 <-ncol(X.2) summary(X.2) ## eta.x is the linear predictor; lambda is predicted probability of missing eta.x.2<-t(modlambda.2$coeff %*% t(X.2)) lambda.2<-exp(eta.x.2)/(1+exp(eta.x.2)) ## It's worth looking at lambda. If some lambda's are too close to 0 ## or 1 we may have some problems. summary(lambda.2) ## These next two loops are used to compute Omega # dlambda.2<-score.2<-matrix(NA,N.2,dimX.2) score.2<-matrix(NA,N,dimX.2) r.2 = r[has7yr == 1] for (i in 1:N) { # dlambda.2[i,]<-lambda.2[i]*(1-lambda.2[i])*X.2[i,] score.2[i,]<-(r[i]-lambda.2[i])*X.2[i,] } Omega.2<-matrix(NA,dimX.2,dimX.2) for (i in 1:dimX.2) { for (j in 1:dimX.2) { Omega.2[i,j]<-sum(score.2[,i]*score.2[,j], na.rm = T)/sum(!is.na(lambda.2)) } } ## Third part of weights - probability of verification N.3 = sum(pca) modlambda.3 = glm(V ~ Z + y+ psa0 + age + white + aa + fampca + psalast + cause_psa + cause_dre + pnegbx, family = binomial, data = pcpt, subset = pca == 1) summary(modlambda.3) X.3 <-cbind(1,Z,y,psa0,age,white,aa,fampca,psalast,cause_psa,cause_dre,pnegbx)*pca dimX.3 <-ncol(X.3) summary(X.3) ## eta.x is the linear predictor; lambda is predicted probability of missing eta.x.3<-t(modlambda.3$coeff %*% t(X.3)) lambda.3<-exp(eta.x.3)/(1+exp(eta.x.3)) ## It's worth looking at lambda. If some lambda's are too close to 0 ## or 1 we may have some problems. summary(lambda.3[pcpt$V ==1]) summary(lambda.3[pcpt$V ==0]) boxplot(split(lambda.3,pcpt$V)) ## These next two loops are used to compute Omega # dlambda.3<-score.3<-matrix(NA,N.3,dimX.3) score.3<-matrix(NA,N,dimX.3) r.3 = pcpt$V for (i in 1:N) { # dlambda.3[i,]<-lambda.3[i]*(1-lambda.3[i])*X.3[i,] score.3[i,]<-(r.3[i]-lambda.3[i])*X.3[i,] } Omega.3<-matrix(NA,dimX.3,dimX.3) for (i in 1:dimX.3) { for (j in 1:dimX.3) { Omega.3[i,j]<-sum(score.3[,i]*score.3[,j], na.rm = T)/sum(!is.na(lambda.3)) } } # combine two matrixes into one Omega score = cbind(score.1,score.2,score.3) Omega = matrix(0,dimX.1 + dimX.2+dimX.3,dimX.1 + dimX.2+dimX.3) Omega[1:dimX.1,1:dimX.1] = solve(Omega.1) Omega[(dimX.1 + 1):(dimX.1+dimX.2),(dimX.1 + 1):(dimX.1+dimX.2)] = solve(Omega.2) Omega[(dimX.1 + dimX.2 + 1):(dimX.1+dimX.2+dimX.3),(dimX.1 + dimX.2 + 1):(dimX.1+dimX.2+dimX.3)] = solve(Omega.3) ## The MCAR estimating equations are weighted by w w<-ifelse(hasbx == 0 , 0, 1/(lambda.1*lambda.2)) w.3 = ifelse(pcpt$V == 0,0,1/lambda.3) postscript(file="ipweights.eps") par(mfrow=c(1,2)) hist(w[hasbx!=0],xlab=expression(1/lambda),main="") hist(w[pcpt$V!=0]*w.3[pcpt$V!=0],xlab=expression(1/(lambda* lambda[q])),main="") dev.off() summary(w[hasbx!=0]) max(w[hasbx!=0])/min(w[hasbx!=0]) summary(w[pcpt$V!=0]*w.3[pcpt$V!=0]) max(w[pcpt$V!=0]*w.3[pcpt$V!=0],na.rm=TRUE)/min(w[pcpt$V!=0]*w.3[pcpt$V!=0],na.rm=TRUE) summary(w[pcpt$V!=0&z==1]*w.3[pcpt$V!=0&z==1]) summary(w[pcpt$V!=0&z==0]*w.3[pcpt$V!=0&z==0]) summary((w[pcpt$V!=0&z==1]*w.3[pcpt$V!=0&z==1])[w[pcpt$V!=0&z==1]*w.3[pcpt$V!=0&z==1]<49]) summary((w[pcpt$V!=0&z==0]*w.3[pcpt$V!=0&z==0])[w[pcpt$V!=0&z==0]*w.3[pcpt$V!=0&z==0]<46]) hist(w[pcpt$V!=0&z==1]*w.3[pcpt$V!=0&z==1]) hist(w[pcpt$V!=0&z==0]*w.3[pcpt$V!=0&z==0]) hgpt[w[pcpt$V!=0&z==1]*w.3[pcpt$V!=0&z==1]==max(w[pcpt$V!=0&z==1]*w.3[pcpt$V!=0&z==1],na.rm=TRUE)] hgpt[w[pcpt$V!=0&z==0]*w.3[pcpt$V!=0&z==0]==max(w[pcpt$V!=0&z==0]*w.3[pcpt$V!=0&z==0],na.rm=TRUE)] ## redefine y y.old = y y = hgpt table(y[z==0&V==1],y.old[z==0&V==1]) table(y[z==1&V==1],y.old[z==1&V==1]) ## Estimating P(S=1|Z=0) and P(S=1|Z=1) p.hat = function(S,Z,W,NPV) { sum((S + (1-S)*(1-NPV))*Z*W,na.rm=T)/sum(Z*W,na.rm=T) } p0hat<- p.hat(s,1-z,NPV = 1, w) p1hat<- p.hat(s,z,NPV = 1, w) ## Now I will need to estimate at fixed values of phi, beta0, and ## beta1. Since these parameters are unknown, I perform a sensitivity ## analysis over a range of values. ## I define a function 'sensanalysis' which performs a sensitivity analysis at an ## assumed value of phi (varies beta0 and beta1) sensanalysis<-function(phi, beta0, beta1) { ## Re-parameterizing pi<-(1-phi)*p1hat ## Estimating alpha0 and alpha1 eealpha0<-function (a) { (sum(w.3*w*(1-z)*s*((1+exp(-a-beta0*y))^(-1)-pi/p0hat),na.rm=T))^2 } eealpha1<-function (a) { (sum(w.3*w*z*s*((1+exp(-a-beta1*y))^(-1)-pi/p1hat),na.rm=T))^2 } min<-optimize(eealpha0, c(-50,50)) alpha0hat<-min$minimum min<-optimize(eealpha1, c(-50,50)) alpha1hat<-min$minimum eem0<-function (a) { (sum(w.3*w*(1-z)*s*(a-y*(1+exp(-alpha0hat-beta0*y))^(-1)*p0hat/pi),na.rm=T))^2 } min<-optimize(eem0, c(0,1)) m0hat<-min$minimum eem1<-function (a) { (sum(w.3*w*z*s*(a-y*(1+exp(-alpha1hat-beta1*y))^(-1)*p1hat/pi),na.rm=T))^2 } min<-optimize(eem1, c(0,1)) m1hat<-min$minimum ## Estimating equations (these should sum to 0) V1<-w*(1-z)*(p0hat-s) V2<-w*z*(p1hat-s) V3<-w.3*w*(1-z)*s*((1+exp(-alpha0hat-beta0*y))^(-1)-pi/p0hat) V4<-w.3*w*z*s*((1+exp(-alpha1hat-beta1*y))^(-1)-pi/p1hat) V5<-w.3*w*(1-z)*s*(m0hat-y*(1+exp(-alpha0hat-beta0*y))^(-1)*p0hat/pi) V6<-w.3*w*z*s*(m1hat-y*(1+exp(-alpha1hat-beta1*y))^(-1)*p1hat/pi) ## Estimating variance. Gamma = matrix(0,ncol = 6, nrow = 6) Gamma[1,1] = sum(w*(1-z),na.rm=T)/N Gamma[2,2] = sum(w*z,na.rm=T)/N Gamma[3,1] = sum(w.3*w*(1-z)*s*(1-phi)*p1hat/p0hat^2,na.rm=T)/N Gamma[3,2] = sum(-w.3*w*(1-z)*s*(1-phi)/p0hat,na.rm=T)/N Gamma[3,3] = sum(w.3*w*(1-z)*s*exp(alpha0hat+beta0*y)/(1+exp(alpha0hat+beta0*y))^2,na.rm=T)/N Gamma[4,4] = sum(w.3*w*z*s*exp(alpha1hat+beta1*y)/(1+exp(alpha1hat+beta1*y))^2,na.rm=T)/N Gamma[5,1] = sum(-w.3*w*(1-z)*s*y*(1+exp(-alpha0hat-beta0*y))^(-1)/(p1hat*(1-phi)),na.rm=T)/N Gamma[5,2] = sum(w.3*w*(1-z)*s*y*(1+exp(-alpha0hat-beta0*y))^(-1)*p0hat/(p1hat^2*(1-phi)),na.rm=T)/N Gamma[5,3] = sum(-w.3*w*(1-z)*s*y*exp(alpha0hat+beta0*y)/(1+exp(alpha0hat+beta0*y))^2*p0hat/(p1hat*(1-phi)),na.rm=T)/N Gamma[5,5] = sum(w.3*w*(1-z)*s,na.rm=T)/N Gamma[6,4] = sum(-w.3*w*z*s*y*exp(alpha1hat+beta1*y)/(1+exp(alpha1hat+beta1*y))^2/(1-phi),na.rm=T)/N Gamma[6,6] = sum(w.3*w*z*s,na.rm=T)/N # G.gamma (B matrix in RRZ) b1<-b2<-b3<-b4<-b5<-b6<-NULL for (jj in 1:dimX.1) { b1[jj]<-sum(V1*(1-lambda.1)*X.1[,jj],na.rm=T) b2[jj]<-sum(V2*(1-lambda.1)*X.1[,jj],na.rm=T) b3[jj]<-sum(V3*(1-lambda.1)*X.1[,jj],na.rm=T) b4[jj]<-sum(V4*(1-lambda.1)*X.1[,jj],na.rm=T) b5[jj]<-sum(V5*(1-lambda.1)*X.1[,jj],na.rm=T) b6[jj]<-sum(V6*(1-lambda.1)*X.1[,jj],na.rm=T) } for (jj in (dimX.1+1):(dimX.1 + dimX.2)) { b1[jj]<-sum(V1*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) b2[jj]<-sum(V2*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) b3[jj]<-sum(V3*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) b4[jj]<-sum(V4*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) b5[jj]<-sum(V5*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) b6[jj]<-sum(V6*(1-lambda.2)*X.2[,jj-dimX.1],na.rm=T) } #************************ for (jj in (dimX.1+dimX.2+1):(dimX.1 + dimX.2 + dimX.3)) { b1[jj]<-0 b2[jj]<-0 b3[jj]<-sum(V3*(1-lambda.3)*X.3[,jj-dimX.1-dimX.2],na.rm=T) b4[jj]<-sum(V4*(1-lambda.3)*X.3[,jj-dimX.1-dimX.2],na.rm=T) b5[jj]<-sum(V5*(1-lambda.3)*X.3[,jj-dimX.1-dimX.2],na.rm=T) b6[jj]<-sum(V6*(1-lambda.3)*X.3[,jj-dimX.1-dimX.2],na.rm=T) } G.gamma = rbind(b1,b2,b3,b4,b5,b6)/N # MM matrix temp.MM = matrix(0, nrow = N, ncol = (dimX.1 + dimX.2 +dimX.3)^2) for(m in 1:N) { temp.MM[m,] = as.vector(score[m,] %o% score[m,]) } MM = matrix(colSums(temp.MM,na.rm = T), nrow = (dimX.1 + dimX.2+dimX.3)) /N # PSI psi.matrix = - solve(MM) %*% t(score) # adjusted estimating equation g = rbind(V1,V2,V3,V4,V5,V6) temp.g = matrix(0, nrow = N, ncol = (6)^2) for(m in 1:N) { temp.g[m,] = as.vector(g[,m] %o% g[,m]) } g.tilda = g + G.gamma %*% psi.matrix temp.gg = matrix(0, nrow = N, ncol = (6)^2) for(m in 1:N) { temp.gg[m,] = as.vector(g.tilda[,m] %o% g.tilda[,m]) } # adjusted variance of estimating equations # naive EGG = matrix(colSums(temp.g,na.rm = T), nrow = (6)) / N # adjusted GG = matrix(colSums(temp.gg,na.rm = T), nrow = (6)) / N varpars.new = (solve(Gamma) %*% GG %*% t(solve(Gamma)))/N varpars.naive = (solve(Gamma) %*% EGG %*% t(solve(Gamma)))/N # estimate of and variance of ACE ace <-m1hat-m0hat rrc = m1hat / m0hat varace<-varpars.new[6,6]+varpars.new[5,5]-2*varpars.new[5,6] varace.n <-varpars.naive[6,6] +varpars.naive[5,5]-2*varpars.naive[5,6] varrrc = c(rep(0,4), 1/m0hat, -m1hat/m0hat^2) %*% varpars.new %*% c(rep(0,4), 1/m0hat, -m1hat/m0hat^2) return(list(pi = pi, ace = ace, varace = varace, varace.naive = varace.n,rrc = rrc, varrrc = varrrc)) } b0 = -20:20/8 b1 = -20:20/8 phi<-.01 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) date() for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } date() save(ace,vara,file="relaxmons-ver-01a.Rdata") phi<-.05 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-ver-05a.Rdata") phi<-.1 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-ver-10a.Rdata") phi<-.2 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-ver-20a.Rdata") phi<-.5 ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-ver-50a.Rdata") phi<-1-p0hat ace = matrix(NA, nrow = length(b0), ncol = length(b1)) vara = matrix(NA, nrow = length(b0), ncol = length(b1)) for(i in 1:length(b0)) { for(j in 1:length(b1)) { temp <-sensanalysis(phi,beta0 = b0[i], beta1 = b1[j]) ace[i,j] = temp$ace vara[i,j] = temp$varace } } save(ace,vara,file="relaxmons-ver-inda.Rdata") ############# Figure 5 of Shepherd, Redman, and Ankerst (2008) #postscript("relaxmons-vera.eps",horiz=F) par(mfrow=c(3,2),mar=c(4,4,2,3)) b0<-b1<--20:20/8 load("relaxmons-ver-01a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.99)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(1.05),log(1.05),log(1.35),log(1.35),log(1.05)),c(log(1.05),log(1.35),log(1.35),log(1.05),log(1.05))) ## Thompson points(log(1.2),log(1.2),pch=4) load("relaxmons-ver-05a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(1),gray(.8)),xlab="",ylab="",axis=1.2,axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.95)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(2),log(2),log(4),log(4),log(2)),c(log(.25),log(.5),log(.5),log(.25),log(.25))) ### Dr Kristal lines(c(log(1.05),log(1.05),log(1.35),log(1.35),log(1.05)),c(log(1.05),log(1.35),log(1.35),log(1.05),log(1.05))) ## Thompson load("relaxmons-ver-10a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.90)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(2),log(2),log(4),log(4),log(2)),c(log(.25),log(.5),log(.5),log(.25),log(.25))) ### Dr Kristal points(log(3),log(1/3),pch=4) load("relaxmons-ver-20a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.80)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() lines(c(log(2),log(2),log(4),log(4),log(2)),c(log(.25),log(.5),log(.5),log(.25),log(.25))) ### Dr Kristal load("relaxmons-ver-50a.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, expression(phi==0.50)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() load("relaxmons-ver-inda.Rdata") lower<-ace-1.96*sqrt(vara) upper<-ace+1.96*sqrt(vara) rejecth0<-ifelse(lower>0,1,ifelse(upper<0,-1,0)) image(b0,b1,rejecth0,col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="",axes=FALSE) contour(b0,b1,ace,add=T,lty=3,levels=c((0:40)/20-1),labcex=0.8) mtext(side=1, line=2.4, expression("OR"[0]==exp(beta[0]))) mtext(side=2, line=2.4, expression("OR"[1]==exp(beta[1]))) mtext(side=3, line=0, quote("S(0),S(1) Independent, " * phi==0.24)) axis(side=1,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) axis(side=2,at=c(-2:2),labels=c(0.14,0.37,1,2.7,7.4),cex.axis=1.2) box() #dev.off()