##### 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,
###### 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],1-TE)])
ace.hhs.low<-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],1-VEboot[ii])])
ace.hhs.low.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()