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