##### 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()