#Commands to reproduce the analyses presented in
#Poleto, F.Z., Singer, J.M. and Paulino, C.D. (2010). Comparing diagnostic tests with missing data. To apper in J. Appl. Stat.

#The functions of the Catdata package are available at http://www.poleto.com/missing.html, or you can load using the command
source("http://www.poleto.com/Catdata.r") #for CCA, ICA, AICA, MCAR and MAR

#If you don't have the package gee installed, you cannot reproduce the DICA (marginal analysis)
require(gee)



#Entire Data in Appendix/Table 7 (0=negative, 1=positive, NA=missing)
Id<-1:219
D<-c(0,0,0,1,1,1,1,0,1,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,1,1,0,1,1,0,1,0,1,1,0,0,0,1,0,0,1,0,0,0,1,0,1,1,1,0,0,0,1,1,0,0,1,0,0,0,0,1,1,1,0,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,1,1,1,0,0,0,0,1,0,1,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,1,0,1,1,1,1,1,0,1,0,1,0,1,0,0,0,1,0,1,0,0,0,0,1,0,1,1,0,1,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,1,1,0,1,1,0,1,1,0,0,0,0,0,1,1,0,1,1,0,0,1,0,1,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0,1,1,1,0,1,1,1,0,0,1,1,1,0,0,0,0,0,1,0,1,1,0,1,0,1)
MR<-c(0,0,0,NA,NA,NA,0,0,0,NA,NA,0,NA,NA,NA,0,NA,NA,0,0,0,0,0,NA,NA,0,NA,NA,0,1,0,1,1,NA,0,0,0,0,NA,0,NA,NA,NA,0,NA,NA,0,0,NA,NA,NA,NA,0,1,0,NA,NA,NA,NA,0,NA,NA,NA,NA,0,0,1,NA,0,1,NA,NA,0,0,NA,0,0,NA,NA,NA,NA,NA,NA,0,NA,NA,NA,NA,1,NA,0,0,0,0,NA,NA,NA,NA,0,1,1,NA,NA,0,NA,NA,NA,NA,0,1,1,NA,0,NA,NA,1,1,NA,NA,NA,1,NA,NA,NA,0,0,1,NA,0,0,NA,NA,NA,NA,NA,0,1,0,0,0,0,0,0,NA,0,0,NA,0,0,0,0,NA,1,0,NA,0,1,0,NA,NA,0,NA,NA,NA,0,0,NA,NA,0,NA,NA,0,NA,NA,NA,NA,NA,NA,0,NA,NA,NA,0,0,0,0,NA,NA,NA,NA,NA,0,0,0,0,1,NA,NA,NA,NA,NA,NA,0,NA,0,NA,NA,NA,1,NA,1,0,0,NA,NA,0,NA,0,NA)
US<-c(0,0,0,1,1,0,1,0,1,0,1,0,0,0,0,0,1,0,1,0,0,1,0,0,0,1,1,0,1,1,0,1,0,1,1,0,0,0,1,0,0,1,0,0,0,0,0,1,1,1,0,0,0,1,1,0,0,1,0,0,0,0,1,1,1,0,1,1,0,1,0,1,0,0,1,0,1,0,0,1,1,1,0,0,1,1,1,1,0,0,0,0,1,0,1,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,1,0,1,1,1,1,1,0,1,0,1,0,1,0,0,0,1,0,1,0,0,0,0,1,0,1,1,0,1,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,1,0,1,1,0,0,1,0,0,0,0,0,1,1,0,1,1,0,0,1,0,0,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0,1,1,1,0,1,1,1,0,0,1,1,1,0,0,0,0,0,1,0,1,1,0,1,0,1)
EC<-c(NA,NA,0,NA,NA,0,NA,NA,NA,NA,NA,NA,1,NA,NA,NA,NA,NA,NA,NA,NA,1,NA,NA,NA,0,NA,0,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0,NA,NA,0,NA,NA,NA,NA,NA,NA,1,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0,1,NA,1,0,NA,NA,NA,0,NA,NA,NA,NA,NA,NA,NA,NA,1,NA,NA,1,NA,NA,0,NA,NA,NA,NA,NA,NA,NA,NA,1,NA,NA,NA,NA,0,NA,NA,1,NA,1,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0,0,NA,NA,NA,1,1,NA,NA,NA,NA,NA,NA,NA,NA,1,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,1,NA,NA,NA,NA,NA,NA,NA,0,NA)



#Contrasts for comparison between MR and EC
contr<-rbind(c(1,-1,0,0),c(0,0,1,-1))



#CCA of MR and EC (Tables 3 and 4)
table(EC,D,MR) #or ftable(MR,EC,D), (complete) data as presented in Table 1 and Table 2 (W=1)

freq.cca<-c(table(D,EC,MR)) #organizing the cells in a different way
freq.cca[freq.cca==0]<-0.001 #substituting zero frequencies by a small value

MR.EC.cca<-readCatdata(TF=freq.cca)

MR.EC.cca.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.cca,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),X=diag(4))
MR.EC.cca.sens #est.and s.e. for Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - CCA
round(cov2cor(MR.EC.cca.sens$VFu),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - CCA (Table 4)
MR.EC.cca.sens2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.cca,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),A3=contr,X=diag(2))
MR.EC.cca.sens2 #est.and s.e. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - CCA (Table 3)
round(cov2cor(MR.EC.cca.sens2$VFu),3) #Correlations between the parameter estimates: Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - CCA
waldTest(obj=MR.EC.cca.sens,C=contr) #P-value of the test (H_S,H_E) - CCA

MR.EC.cca.ppv<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.cca,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,0,0,0,1,1,1,1),
  c(0,0,1,1,0,0,1,1),
  c(1,1,1,1,0,0,0,0),
  c(1,1,0,0,1,1,0,0)),
  A2=cbind(diag(4),-diag(4)),X=diag(4))
MR.EC.cca.ppv #est.and s.e. for PPV_{MR}, PPV_{EC}, NPV_{MR}, NPV_{EC} - CCA
round(cov2cor(MR.EC.cca.ppv$VFu),3) #Correlations among the parameter estimates: PPV_{MR}, PPV_{EC}, NPV_{MR}, NPV_{EC} - CCA (Table 4)
MR.EC.cca.ppv2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.cca,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,0,0,0,1,1,1,1),
  c(0,0,1,1,0,0,1,1),
  c(1,1,1,1,0,0,0,0),
  c(1,1,0,0,1,1,0,0)),
  A2=cbind(diag(4),-diag(4)),A3=contr,X=diag(2))
print(MR.EC.cca.ppv2,digits=5) #est.and s.e. for PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - CCA (Table 3)
round(cov2cor(MR.EC.cca.ppv2$VFu),3) #Correlations between the parameter estimates: PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - CCA
waldTest(obj=MR.EC.cca.ppv,C=contr) #P-value of the test (H_P,H_N) - CCA



#ICA of MR and EC (Table 3)
MR.EC.ica<-readCatdata(TF=rbind( c(table(D[is.na(EC)],MR[is.na(EC)])) , c(table(D[is.na(MR)],EC[is.na(MR)])) ))

MR.EC.ica.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.ica,
  A1=diag(2)%x%rbind(matrix(c(1,rep(0,6),1),2),t(rep(1,2)%x%diag(2))),A2=(diag(2)%x%t(c(1,-1)%x%diag(2)))[c(2,4,1,3),],X=diag(4))
MR.EC.ica.sens #est.and s.e. for Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - ICA
round(cov2cor(MR.EC.ica.sens$VFu),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - ICA
MR.EC.ica.sens2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.ica,A3=contr,
  A1=diag(2)%x%rbind(matrix(c(1,rep(0,6),1),2),t(rep(1,2)%x%diag(2))),A2=(diag(2)%x%t(c(1,-1)%x%diag(2)))[c(2,4,1,3),],X=diag(2))
MR.EC.ica.sens2 #est.and s.e. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - ICA (Table 3)
round(cov2cor(MR.EC.ica.sens2$VFu),3) #Correlations between the parameter estimates: Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - ICA
waldTest(obj=MR.EC.ica.sens,C=contr) #P-value of the test (H_S,H_E) - ICA

MR.EC.ica.ppv<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.ica,
  A1=diag(2)%x%rbind(matrix(c(1,rep(0,6),1),2),t(diag(2)%x%rep(1,2))),A2=(diag(2)%x%t(c(1,-1)%x%diag(2)))[c(2,4,1,3),],X=diag(4))
MR.EC.ica.ppv #est.and s.e. for PPV_{MR}, PPV_{EC}, NPV_{MR}, NPV_{EC} - ICA
round(cov2cor(MR.EC.ica.ppv$VFu),3) #Correlations among the parameter estimates: PPV_{MR}, PPV_{EC}, NPV_{MR}, NPV_{EC} - ICA
MR.EC.ica.ppv2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.ica,A3=contr,
  A1=diag(2)%x%rbind(matrix(c(1,rep(0,6),1),2),t(diag(2)%x%rep(1,2))),A2=(diag(2)%x%t(c(1,-1)%x%diag(2)))[c(2,4,1,3),],X=diag(2))
MR.EC.ica.ppv2 #est.and s.e. for PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - ICA (Table 3)
round(cov2cor(MR.EC.ica.ppv2$VFu),3) #Correlations between the parameter estimates: PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - ICA
waldTest(obj=MR.EC.ica.ppv,C=contr) #P-value of the test (H_P,H_N) - ICA



#AICA of MR and EC (Table 3)
MR.EC.aica<-readCatdata(TF=rbind( c(table(D,MR)) , c(table(D,EC)) )) #The following commands are exactly equal to the ICA, except that now they are using a different data set

MR.EC.aica.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.aica,
  A1=diag(2)%x%rbind(matrix(c(1,rep(0,6),1),2),t(rep(1,2)%x%diag(2))),A2=(diag(2)%x%t(c(1,-1)%x%diag(2)))[c(2,4,1,3),],X=diag(4))
MR.EC.aica.sens #est.and s.e. for Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - AICA
round(cov2cor(MR.EC.aica.sens$VFu),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - AICA
MR.EC.aica.sens2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.aica,A3=contr,
  A1=diag(2)%x%rbind(matrix(c(1,rep(0,6),1),2),t(rep(1,2)%x%diag(2))),A2=(diag(2)%x%t(c(1,-1)%x%diag(2)))[c(2,4,1,3),],X=diag(2))
MR.EC.aica.sens2 #est.and s.e. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - AICA (Table 3)
round(cov2cor(MR.EC.aica.sens2$VFu),3) #Correlations between the parameter estimates: Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - AICA

MR.EC.aica.ppv<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.aica,
  A1=diag(2)%x%rbind(matrix(c(1,rep(0,6),1),2),t(diag(2)%x%rep(1,2))),A2=(diag(2)%x%t(c(1,-1)%x%diag(2)))[c(2,4,1,3),],X=diag(4))
MR.EC.aica.ppv #est.and s.e. for PPV_{MR}, PPV_{EC}, NPV_{MR}, NPV_{EC} - AICA
round(cov2cor(MR.EC.aica.ppv$VFu),3) #Correlations among the parameter estimates: PPV_{MR}, PPV_{EC}, NPV_{MR}, NPV_{EC} - AICA
MR.EC.aica.ppv2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.aica,A3=contr,
  A1=diag(2)%x%rbind(matrix(c(1,rep(0,6),1),2),t(diag(2)%x%rep(1,2))),A2=(diag(2)%x%t(c(1,-1)%x%diag(2)))[c(2,4,1,3),],X=diag(2))
MR.EC.aica.ppv2 #est.and s.e. for PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - AICA (Table 3)
round(cov2cor(MR.EC.aica.ppv2$VFu),3) #Correlations between the parameter estimates: PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - AICA
waldTest(obj=MR.EC.aica.ppv,C=contr) #P-value of the test (H_P,H_N) - AICA



#DICA of MR and EC (Tables 3 and 4)
contr2<-rbind(c(1,-1,0,0),c(0,0,-1,1)) #The contrasts of DICA have different signs for the second row because of the different parametrization

#Rearranging the dataset for the marginal analysis (GEE) / DICA
GEEdata<-data.frame(Id=rep(Id,3),D=rep(D,3),ScreeningTest=rep(c("EC","MR","US"),rep(219,3)),TestResult=c(EC,MR,US),diag(3)%x%D,diag(3)%x%(1-D),ECT1=c(EC,rep(0,2*219)),MRT1=c(rep(0,219),MR,rep(0,219)),UST1=c(rep(0,2*219),US),ECT0=c(1-EC,rep(0,2*219)),MRT0=c(rep(0,219),1-MR,rep(0,219)),UST0=c(rep(0,2*219),1-US))
colnames(GEEdata)[5:10]<-c("ECD1","MRD1","USD1","ECD0","MRD0","USD0")
GEEdata<-na.omit(GEEdata)
GEEdata<-GEEdata[order(GEEdata$Id),]

MR.EC.dica.sens<-gee(TestResult~MRD1+ECD1+MRD0+ECD0-1,id=Id,data=subset(GEEdata,ScreeningTest!="US"),family=binomial(link="identity"),corstr="independence",scale.fix=TRUE,scale.value=1)
summary(MR.EC.dica.sens) #est.and s.e. for Sens_{MR}, Sens_{EC}, 1-Spec_{MR}, 1-Spec_{EC} - DICA (note that the last two parameters are different from the other analyses)
round(cov2cor(MR.EC.dica.sens$robust.variance),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, 1-Spec_{MR}, 1-Spec_{EC} - DICA (Table 4)
c(contr2%*%MR.EC.dica.sens$coef) #est. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - DICA (Table 3)
sqrt(diag(contr2%*%MR.EC.dica.sens$robust.variance%*%t(contr2))) #s.e. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - DICA (Table 3)

MR.EC.dica.ppv<-gee(D~MRT1+ECT1+MRT0+ECT0-1,id=Id,data=subset(GEEdata,ScreeningTest!="US"),family=binomial(link="identity"),corstr="independence",scale.fix=TRUE,scale.value=1)
summary(MR.EC.dica.ppv) #est.and s.e. for PPV_{MR}, PPV_{EC}, 1-NPV_{MR}, 1-NPV_{EC} - DICA (note that the last two parameters are different from the other analyses)
round(cov2cor(MR.EC.dica.ppv$robust.variance),3) #Correlations among the parameter estimates: PPV_{MR}, PPV_{EC}, 1-NPV_{MR}, 1-NPV_{EC} - DICA (Table 4)
c(contr2%*%MR.EC.dica.ppv$coef) #est. for PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - DICA (Table 3)
sqrt(diag(contr2%*%MR.EC.dica.ppv$robust.variance%*%t(contr2))) #s.e. for PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - DICA (Table 3)



#MCAR/MAR of MR and EC (Tables 3 and 4)
MisMR<-is.na(MR)
MisEC<-is.na(EC)
table(MisMR,MisEC) #Sample sizes in each of the four missingness patterns
freq.table2<-c( ftable(EC,D,MR) , table(D[MisEC],MR[MisEC]) , table(D[MisMR],EC[MisMR])  , table(D[MisMR & MisEC]) )
freq.table2[freq.table2==0]<-0.001
MR.EC.mar.TF<-readCatdata(TF=freq.table2,
  Zp=cbind(diag(2)%x%rep(1,2)%x%diag(2), rep(1,2)%x%diag(4), rep(1,4)%x%diag(2)),Rp=c(4,4,2))
#Zp informs to the function the patterns of missingness Z_2,...Z_T; See Paulino (1991) or Poleto et al. (2009a) to understand the matrices Z_t - See also the results of summary(MR.EC.mcar.TF)
#Rp has the number of columns of each matrix Z_t, t=2,...,T


#under MCAR
MR.EC.mcar<-satMarML(MR.EC.mar.TF,missing="MCAR",maxit=200) #Default s.e. are computed assuming MAR if missing="MCAR" not used
MR.EC.mcar #P-values of likelihood ratio and Pearson tests of MCAR conditional MAR (first paragraph of Section 4)

MR.EC.mcar.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.mcar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),X=diag(4))
MR.EC.mcar.sens #est.and s.e. for Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - MCAR
round(cov2cor(MR.EC.mcar.sens$VFu),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - MCAR (Table 4)
MR.EC.mcar.sens2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.mcar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),A3=contr,X=diag(2))
MR.EC.mcar.sens2 #est.and s.e. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - MCAR (Table 3)
round(cov2cor(MR.EC.mcar.sens2$VFu),3) #Correlations between the parameter estimates: Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - MCAR
waldTest(obj=MR.EC.mcar.sens,C=contr) #P-value of the test (H_S,H_E) - MCAR

MR.EC.mcar.ppv<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.mcar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,0,0,0,1,1,1,1),
  c(0,0,1,1,0,0,1,1),
  c(1,1,1,1,0,0,0,0),
  c(1,1,0,0,1,1,0,0)),
  A2=cbind(diag(4),-diag(4)),X=diag(4))
MR.EC.mcar.ppv #est.and s.e. for PPV_{MR}, PPV_{EC}, NPV_{MR}, NPV_{EC} - MCAR
round(cov2cor(MR.EC.mcar.ppv$VFu),3) #Correlations among the parameter estimates: PPV_{MR}, PPV_{EC}, NPV_{MR}, NPV_{EC} - MCAR (Table 4)
MR.EC.mcar.ppv2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.mcar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,0,0,0,1,1,1,1),
  c(0,0,1,1,0,0,1,1),
  c(1,1,1,1,0,0,0,0),
  c(1,1,0,0,1,1,0,0)),
  A2=cbind(diag(4),-diag(4)),A3=contr,X=diag(2))
MR.EC.mcar.ppv2 #est.and s.e. for PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - MCAR (Table 3)
round(cov2cor(MR.EC.mcar.ppv2$VFu),3) #Correlations between the parameter estimates: PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - MCAR
waldTest(obj=MR.EC.mcar.ppv,C=contr) #P-value of the test (H_P,H_N) - MCAR


#under MAR
MR.EC.mar<-satMarML(MR.EC.mar.TF,maxit=200) #The following commands are exactly equal to the analysis under MCAR, except that now they are using a different source object
MR.EC.mar #P-values of likelihood ratio and Pearson tests of MCAR conditional MAR (first paragraph of Section 4)

MR.EC.mar.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.mar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),X=diag(4))
MR.EC.mar.sens #est.and s.e. for Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - MAR
round(cov2cor(MR.EC.mar.sens$VFu),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - MAR (Table 4)
MR.EC.mar.sens2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.mar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),A3=contr,X=diag(2))
MR.EC.mar.sens2 #est.and s.e. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - MAR (Table 3)
round(cov2cor(MR.EC.mar.sens2$VFu),3) #Correlations between the parameter estimates: Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - MAR
waldTest(obj=MR.EC.mar.sens,C=contr) #P-value of the test (H_S,H_E) - MAR

MR.EC.mar.ppv<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.mar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,0,0,0,1,1,1,1),
  c(0,0,1,1,0,0,1,1),
  c(1,1,1,1,0,0,0,0),
  c(1,1,0,0,1,1,0,0)),
  A2=cbind(diag(4),-diag(4)),X=diag(4))
MR.EC.mar.ppv #est.and s.e. for PPV_{MR}, PPV_{EC}, NPV_{MR}, NPV_{EC} - MAR
round(cov2cor(MR.EC.mar.ppv$VFu),3) #Correlations among the parameter estimates: PPV_{MR}, PPV_{EC}, NPV_{MR}, NPV_{EC} - MAR (Table 4)
MR.EC.mar.ppv2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.mar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,0,0,0,1,1,1,1),
  c(0,0,1,1,0,0,1,1),
  c(1,1,1,1,0,0,0,0),
  c(1,1,0,0,1,1,0,0)),
  A2=cbind(diag(4),-diag(4)),A3=contr,X=diag(2))
MR.EC.mar.ppv2 #est.and s.e. for PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - MAR (Table 3)
round(cov2cor(MR.EC.mar.ppv2$VFu),3) #Correlations between the parameter estimates: PPV_{MR}-PPV_{EC}, NPV_{MR}-NPV_{EC} - MAR
waldTest(obj=MR.EC.mar.ppv,C=contr) #P-value of the test (H_P,H_N) - MAR



#Repeating previous analyses under CCA, MCAR and MAR to the modified dataset (results in the text)
#i.e., moving a single observation from (EC=+) to (EC=-) in (MR=-;D=+;W=1), and another one from (EC=-) to (EC=+) in (MR=+;D=+;W=1)
na.omit(cbind(Id,D,MR,EC)) #complete cases
EC.mod<-EC
EC.mod[subset(Id,EC==1&MR==0&D==1)[1]]<-0
EC.mod[subset(Id,EC==0&MR==1&D==1)[1]]<-1
na.omit(cbind(Id,D,MR,EC.mod)) #complete cases of the modified dataset

#under CCA
freq.mod.cca<-c(table(D,EC.mod,MR)) #organizing the cells in a different way
freq.mod.cca[freq.mod.cca==0]<-0.001 #substituting zero frequencies by a small value
MR.EC.mod.cca<-readCatdata(TF=freq.mod.cca)
MR.EC.mod.cca.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.mod.cca,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),X=diag(4))
MR.EC.mod.cca.sens #est.and s.e. for Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - CCA
round(cov2cor(MR.EC.mod.cca.sens$VFu),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - CCA (text)
MR.EC.mod.cca.sens2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.mod.cca,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),A3=contr,X=diag(2))
MR.EC.mod.cca.sens2 #est.and s.e. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - CCA (text)
round(cov2cor(MR.EC.mod.cca.sens2$VFu),3) #Correlations between the parameter estimates: Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - CCA

#Obviously this modified dataset generates the same frequencies used in the ICA and AICA of the original dataset
rbind( c(table(D[is.na(EC)],MR[is.na(EC)])) , c(table(D[is.na(MR)],EC[is.na(MR)])) )
#is equal to
rbind( c(table(D[is.na(EC.mod)],MR[is.na(EC.mod)])) , c(table(D[is.na(MR)],EC.mod[is.na(MR)])) )
#and
rbind( c(table(D,MR)) , c(table(D,EC)) )
#is equal to
rbind( c(table(D,MR)) , c(table(D,EC.mod)) )

#under DICA
GEEdata.mod<-data.frame(Id=rep(Id,3),D=rep(D,3),ScreeningTest=rep(c("EC","MR","US"),rep(219,3)),TestResult=c(EC.mod,MR,US),diag(3)%x%D,diag(3)%x%(1-D),ECT1=c(EC.mod,rep(0,2*219)),MRT1=c(rep(0,219),MR,rep(0,219)),UST1=c(rep(0,2*219),US),ECT0=c(1-EC.mod,rep(0,2*219)),MRT0=c(rep(0,219),1-MR,rep(0,219)),UST0=c(rep(0,2*219),1-US))
colnames(GEEdata.mod)[5:10]<-c("ECD1","MRD1","USD1","ECD0","MRD0","USD0")
GEEdata.mod<-na.omit(GEEdata.mod)
GEEdata.mod<-GEEdata.mod[order(GEEdata.mod$Id),]

MR.EC.mod.dica.sens<-gee(TestResult~MRD1+ECD1+MRD0+ECD0-1,id=Id,data=subset(GEEdata.mod,ScreeningTest!="US"),family=binomial(link="identity"),corstr="independence",scale.fix=TRUE,scale.value=1)
summary(MR.EC.mod.dica.sens) #est.and s.e. for Sens_{MR}, Sens_{EC}, 1-Spec_{MR}, 1-Spec_{EC} - DICA (note that the last two parameters are different from the other analyses)
round(cov2cor(MR.EC.mod.dica.sens$robust.variance),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, 1-Spec_{MR}, 1-Spec_{EC} - DICA (text)
c(contr2%*%MR.EC.mod.dica.sens$coef) #est. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - DICA
sqrt(diag(contr2%*%MR.EC.mod.dica.sens$robust.variance%*%t(contr2))) #s.e. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - DICA (text)

#MCAR/MAR frequencies
freq.table2.mod<-c( ftable(EC.mod,D,MR) , table(D[MisEC],MR[MisEC]) , table(D[MisMR],EC.mod[MisMR])  , table(D[MisMR & MisEC]) )
freq.table2.mod[freq.table2.mod==0]<-0.001
MR.EC.mod.mar.TF<-readCatdata(TF=freq.table2.mod,Zp=cbind(diag(2)%x%rep(1,2)%x%diag(2), rep(1,2)%x%diag(4), rep(1,4)%x%diag(2)),Rp=c(4,4,2))

#under MCAR
MR.EC.mod.mcar<-satMarML(MR.EC.mod.mar.TF,missing="MCAR",maxit=200)
MR.EC.mod.mcar.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.mod.mcar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),X=diag(4))
MR.EC.mod.mcar.sens #est.and s.e. for Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - MCAR
round(cov2cor(MR.EC.mod.mcar.sens$VFu),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - MCAR (text)
MR.EC.mod.mcar.sens2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.mod.mcar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),A3=contr,X=diag(2))
MR.EC.mod.mcar.sens2 #est.and s.e. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - MCAR (text)
round(cov2cor(MR.EC.mod.mcar.sens2$VFu),3) #Correlations between the parameter estimates: Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - MCAR
waldTest(obj=MR.EC.mod.mcar.sens,C=contr) #P-value of the test (H_S,H_E) - MCAR

#under MAR
MR.EC.mod.mar<-satMarML(MR.EC.mod.mar.TF,maxit=200) #The following commands are exactly equal to the analysis under MCAR, except that now they are using a different source object
MR.EC.mod.mar.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.mod.mar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),X=diag(4))
MR.EC.mod.mar.sens #est.and s.e. for Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - MAR
round(cov2cor(MR.EC.mod.mar.sens$VFu),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, Spec_{MR}, Spec_{EC} - MAR (text)
MR.EC.mod.mar.sens2<-funlinWLS(model=c("lin","exp","lin","log","lin"),obj=MR.EC.mod.mar,A1=rbind(
  c(0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1),
  c(1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0),
  c(0,1,0,1,0,1,0,1),
  c(1,0,1,0,1,0,1,0)),
  A2=cbind(diag(4),-diag(2)%x%rep(1,2)),A3=contr,X=diag(2))
MR.EC.mod.mar.sens2 #est.and s.e. for Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - MAR (text)
round(cov2cor(MR.EC.mod.mar.sens2$VFu),3) #Correlations between the parameter estimates: Sens_{MR}-Sens_{EC}, Spec_{MR}-Spec_{EC} - MAR
waldTest(obj=MR.EC.mod.mar.sens,C=contr) #P-value of the test (H_S,H_E) - MAR



#Analysis of MR, EC and US under MCAR (Table 5)
freq.table1<-c( ftable(MR,EC,D,US) , table(D[MisEC],MR[MisEC],US[MisEC]) , table(D[MisMR],EC[MisMR],US[MisMR]) ,
  table(D[MisMR & MisEC],US[MisMR & MisEC]) )
freq.table1[freq.table1==0]<-0.001
MR.EC.US.mcar.TF<-readCatdata(TF=freq.table1,
  Zp=cbind(diag(4)%x%rep(1,2)%x%diag(2), diag(2)%x%rep(1,2)%x%diag(4), diag(2)%x%rep(1,4)%x%diag(2)),Rp=c(8,8,4))

MR.EC.US.mcar.mar<-satMarML(MR.EC.US.mcar.TF,missing="MCAR",maxit=200)
MR.EC.US.mcar.mar #P-values of likelihood ratio and Pearson tests of MCAR conditional MAR (first paragraph of Section 4)


MR.EC.US.mcar.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.US.mcar.mar,A1=rbind(
  c(1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0),
  c(1,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0),
  c(0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1),
  c(0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1),
  c(1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0),
  c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1)),
  A2=cbind(diag(6),-diag(2)%x%rep(1,3))[c(5,6,4,2,3,1),],X=diag(6))
MR.EC.US.mcar.sens #est.and s.e. for Sens_{MR}, Sens_{EC}, Sens_{US}, Spec_{MR}, Spec_{EC}, Spec_{US} - MCAR (Table 5)
round(cov2cor(MR.EC.US.mcar.sens$VFu),3) #Correlations among the parameter estimates: Sens_{MR}, Sens_{EC}, Sens_{US}, Spec_{MR}, Spec_{EC}, Spec_{US} - MCAR (Table 5)

C1<-rbind(c(1,-1,0,0,0,0),c(1,0,-1,0,0,0))
C2<-rbind(c(0,0,0,1,-1,0),c(0,0,0,1,0,-1))
waldTest(obj=MR.EC.US.mcar.sens,C=C1) #P-value of the test Sens_{MR}=Sens_{EC}=Sens_{US}
waldTest(obj=MR.EC.US.mcar.sens,C=C2) #P-value of the test Spec_{MR}=Spec_{EC}=Spec_{US}
waldTest(obj=MR.EC.US.mcar.sens,C=rbind(C1,C2)) #P-value of both tests simultaneously

waldTest(obj=MR.EC.US.mcar.sens,C=c(1,-1,0,0,0,0)) #P-value of the test Sens_{MR}=Sens_{EC} (Table 5)
waldTest(obj=MR.EC.US.mcar.sens,C=c(1,0,-1,0,0,0)) #P-value of the test Sens_{MR}=Sens_{US} (Table 5)
waldTest(obj=MR.EC.US.mcar.sens,C=c(0,1,-1,0,0,0)) #P-value of the test Sens_{EC}=Sens_{US} (Table 5)

waldTest(obj=MR.EC.US.mcar.sens,C=c(0,0,0,1,-1,0)) #P-value of the test Spec_{MR}=Spec_{EC} (Table 5)
waldTest(obj=MR.EC.US.mcar.sens,C=c(0,0,0,1,0,-1)) #P-value of the test Spec_{MR}=Spec_{US} (Table 5)
waldTest(obj=MR.EC.US.mcar.sens,C=c(0,0,0,0,1,-1)) #P-value of the test Spec_{EC}=Spec_{US} (Table 5)

MR.EC.US.mcar.sens2<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.US.mcar.mar,A1=rbind(
  c(1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0),
  c(1,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0),
  c(0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1),
  c(0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1),
  c(1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0),
  c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1)),
  A2=cbind(diag(6),-diag(2)%x%rep(1,3))[c(5,6,4,2,3,1),],X=rbind(cbind(c(1,1),matrix(0,2,4)),cbind(rep(0,4),diag(4))))
MR.EC.US.mcar.sens2 #Estimates under Sens_{MR}=Sens_{EC} (text)
MR.EC.US.mcar.sens2$beta-qnorm(0.975)*sqrt(diag(MR.EC.US.mcar.sens2$Vbeta)) #(text)
MR.EC.US.mcar.sens2$beta+qnorm(0.975)*sqrt(diag(MR.EC.US.mcar.sens2$Vbeta)) #(text)

Cadhoc1<-rbind(c(-1,1,0,0,0),c(0,0,-1,0,1),c(0,0,1,-1,0))
c(Cadhoc1%*%MR.EC.US.mcar.sens2$beta) #(text)
sqrt(diag(Cadhoc1%*%MR.EC.US.mcar.sens2$Vbeta%*%t(Cadhoc1)))
c(Cadhoc1%*%MR.EC.US.mcar.sens2$beta)-qnorm(0.975)*sqrt(diag(Cadhoc1%*%MR.EC.US.mcar.sens2$Vbeta%*%t(Cadhoc1))) #(text)
c(Cadhoc1%*%MR.EC.US.mcar.sens2$beta)+qnorm(0.975)*sqrt(diag(Cadhoc1%*%MR.EC.US.mcar.sens2$Vbeta%*%t(Cadhoc1))) #(text)





MR.EC.US.mcar.ppv<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.US.mcar.mar,A1=rbind(
  c(1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0),
  c(1,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0),
  c(0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1),
  c(0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1),
  c(1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0),
  c(1,1,1,1,0,0,0,0,1,1,1,1,0,0,0,0),
  c(1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0),
  c(0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1),
  c(0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1),
  c(0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1)),
  A2=cbind(diag(6),-diag(6))[c(5,6,4,2,3,1),],X=diag(6))
MR.EC.US.mcar.ppv #est.and s.e. for PPV_{MR}, PPV_{EC}, PPV_{US}, NPV_{MR}, NPV_{EC}, NPV_{US} - MCAR (Table 5)
round(cov2cor(MR.EC.US.mcar.ppv$VFu),3) #Correlations among the parameter estimates: PPV_{MR}, PPV_{EC}, PPV_{US}, NPV_{MR}, NPV_{EC}, NPV_{US} - MCAR (Table 5)

C1<-rbind(c(1,-1,0,0,0,0),c(1,0,-1,0,0,0))
C2<-rbind(c(0,0,0,1,-1,0),c(0,0,0,1,0,-1))
waldTest(obj=MR.EC.US.mcar.ppv,C=C1) #P-value of the test PPV_{MR}=PPV_{EC}=PPV_{US}
waldTest(obj=MR.EC.US.mcar.ppv,C=C2) #P-value of the test NPV_{MR}=NPV_{EC}=NPV_{US}
waldTest(obj=MR.EC.US.mcar.ppv,C=rbind(C1,C2)) #P-value of both tests simultaneously

waldTest(obj=MR.EC.US.mcar.ppv,C=c(1,-1,0,0,0,0)) #P-value of the test PPV_{MR}=PPV_{EC} (Table 5)
waldTest(obj=MR.EC.US.mcar.ppv,C=c(1,0,-1,0,0,0)) #P-value of the test PPV_{MR}=PPV_{US} (Table 5)
waldTest(obj=MR.EC.US.mcar.ppv,C=c(0,1,-1,0,0,0)) #P-value of the test PPV_{EC}=PPV_{US} (Table 5)

waldTest(obj=MR.EC.US.mcar.ppv,C=c(0,0,0,1,-1,0)) #P-value of the test NPV_{MR}=NPV_{EC} (Table 5)
waldTest(obj=MR.EC.US.mcar.ppv,C=c(0,0,0,1,0,-1)) #P-value of the test NPV_{MR}=NPV_{US} (Table 5)
waldTest(obj=MR.EC.US.mcar.ppv,C=c(0,0,0,0,1,-1)) #P-value of the test NPV_{EC}=NPV_{US} (Table 5)

waldTest(obj=MR.EC.US.mcar.ppv,C=rbind(c(1,-1,0,0,0,0),c(0,0,0,1,-1,0))) #P-value of the test PPV_{MR}=PPV_{EC} and NPV_{MR}=NPV_{EC}

MR.EC.US.mcar.ppv2<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.US.mcar.mar,A1=rbind(
  c(1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0),
  c(1,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0),
  c(0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1),
  c(0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1),
  c(1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0),
  c(1,1,1,1,0,0,0,0,1,1,1,1,0,0,0,0),
  c(1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0),
  c(0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1),
  c(0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1),
  c(0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1)),
  A2=cbind(diag(6),-diag(6))[c(5,6,4,2,3,1),],X=rbind(cbind(c(1,1),matrix(0,2,4)),cbind(rep(0,4),diag(4))))
MR.EC.US.mcar.ppv2 #Estimates under PPV_{MR}=PPV_{EC} (text)
MR.EC.US.mcar.ppv2$beta-qnorm(0.975)*sqrt(diag(MR.EC.US.mcar.ppv2$Vbeta)) #(text)
MR.EC.US.mcar.ppv2$beta+qnorm(0.975)*sqrt(diag(MR.EC.US.mcar.ppv2$Vbeta)) #(text)

Cadhoc2<-rbind(c(-1,1,0,0,0))
c(Cadhoc2%*%MR.EC.US.mcar.ppv2$beta) #(text)
sqrt(diag(Cadhoc2%*%MR.EC.US.mcar.ppv2$Vbeta%*%t(Cadhoc2)))
c(Cadhoc2%*%MR.EC.US.mcar.ppv2$beta)-qnorm(0.975)*sqrt(diag(Cadhoc2%*%MR.EC.US.mcar.ppv2$Vbeta%*%t(Cadhoc2))) #(text)
c(Cadhoc2%*%MR.EC.US.mcar.ppv2$beta)+qnorm(0.975)*sqrt(diag(Cadhoc2%*%MR.EC.US.mcar.ppv2$Vbeta%*%t(Cadhoc2))) #(text)

MR.EC.US.mcar.ppv3<-funlinWLS(model=c("exp","lin","log","lin"),obj=MR.EC.US.mcar.mar,A1=rbind(
  c(1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0),
  c(1,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0),
  c(1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0),
  c(0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1),
  c(0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,1),
  c(0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1),
  c(1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0),
  c(1,1,1,1,0,0,0,0,1,1,1,1,0,0,0,0),
  c(1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0),
  c(0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1),
  c(0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1),
  c(0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1)),
  A2=cbind(diag(6),-diag(6))[c(5,6,4,2,3,1),],X=rbind(cbind(diag(4),rep(0,4)),c(0,0,0,1,0),c(rep(0,4),1)) )
MR.EC.US.mcar.ppv3 #Estimates under NPV_{MR}=NPV_{EC} (text)
MR.EC.US.mcar.ppv3$beta-qnorm(0.975)*sqrt(diag(MR.EC.US.mcar.ppv3$Vbeta)) #(text)
MR.EC.US.mcar.ppv3$beta+qnorm(0.975)*sqrt(diag(MR.EC.US.mcar.ppv3$Vbeta)) #(text)

Cadhoc3<-rbind(c(0,0,0,-1,1))
c(Cadhoc3%*%MR.EC.US.mcar.ppv3$beta) #(text)
sqrt(diag(Cadhoc3%*%MR.EC.US.mcar.ppv3$Vbeta%*%t(Cadhoc3)))
c(Cadhoc3%*%MR.EC.US.mcar.ppv3$beta)-qnorm(0.975)*sqrt(diag(Cadhoc3%*%MR.EC.US.mcar.ppv3$Vbeta%*%t(Cadhoc3))) #(text)
c(Cadhoc3%*%MR.EC.US.mcar.ppv3$beta)+qnorm(0.975)*sqrt(diag(Cadhoc3%*%MR.EC.US.mcar.ppv3$Vbeta%*%t(Cadhoc3))) #(text)