#Comandos para reproduzir as análises do Exemplo 5, pp.126-135, de
#Poleto, F.Z. (2006). Análise de dados categorizados com omissão. Dissertação de mestrado. Instituto de Matemática e Estatística, Universidade de São Paulo.

#As funções estão disponíveis em http://www.poleto.com/missing.html, ou pode-se carregá-las diretamente usando o comando
source("http://www.poleto.com/Catdata.r")



#ACC de RM e EC (Tabelas 3.23 e 3.26)
RM.EC.acc<-readCatdata(TF=c(6,1,1,2, 0.001,1,0.001,2))

RM.EC.acc.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.EC.acc,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))
print(RM.EC.acc.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{EC}, Espec_{RM}, Espec_{EC} - ACC (Tabela 3.23)
print(waldTest(obj=RM.EC.acc.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - ACC (Tabela 3.23)
print(waldTest(obj=RM.EC.acc.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - ACC (Tabela 3.23)
print(waldTest(obj=RM.EC.acc.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - ACC (Tabela 3.23)
round(cov2cor(RM.EC.acc.sens$VFu),3) #Correlações entre as estimativas dos parâmetros: Sens_{RM}, Sens_{EC}, Espec_{RM}, Espec_{EC} - ACC (Tabela 3.26)

RM.EC.acc.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.EC.acc,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))
print(RM.EC.acc.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{EC}, VPN_{RM}, VPN_{EC} - ACC (Tabela 3.23)
print(waldTest(obj=RM.EC.acc.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - ACC (Tabela 3.23)
print(waldTest(obj=RM.EC.acc.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - ACC (Tabela 3.23)
print(waldTest(obj=RM.EC.acc.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - ACC (Tabela 3.23)
round(cov2cor(RM.EC.acc.vpp$VFu),3) #Correlações entre as estimativas dos parâmetros: VPP_{RM}, VPP_{EC}, VPN_{RM}, VPN_{EC} - ACC (Tabela 3.26)

CC<-rbind(c(1,-1,0,0),c(0,0,1,-1))
round(cov2cor(CC%*%RM.EC.acc.vpp$VFu%*%t(CC)),3) #Correlation between linear contrasts mentioned in the paragraph between Tabelas 3.25 e 3.26 - ACC



#ACI de RM e EC (Tabela 3.23)
RM.EC.aci<-readCatdata(TF=rbind(c(51,22,5,13),c(3,5,3,6)))

RM.EC.aci.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.EC.aci,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))
print(RM.EC.aci.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{EC}, Espec_{RM}, Espec_{EC} - ACI (Tabela 3.23)
print(waldTest(obj=RM.EC.aci.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - ACI (Tabela 3.23)
print(waldTest(obj=RM.EC.aci.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - ACI (Tabela 3.23)
print(waldTest(obj=RM.EC.aci.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - ACI (Tabela 3.23)

RM.EC.aci.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.EC.aci,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))
print(RM.EC.aci.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{EC}, VPN_{RM}, VPN_{EC} - ACI (Tabela 3.23)
print(waldTest(obj=RM.EC.aci.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - ACI (Tabela 3.23)
print(waldTest(obj=RM.EC.aci.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - ACI (Tabela 3.23)
print(waldTest(obj=RM.EC.aci.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - ACI (Tabela 3.23)



#ACAI de RM e EC (Tabela 3.23)
RM.EC.acai<-readCatdata(TF=rbind(c(58,25,5,16),c(9,7,4,10))) #The following commands are exactly equal to the ACI, except that now they are using a different data set

RM.EC.acai.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.EC.acai,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))
print(RM.EC.acai.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{EC}, Espec_{RM}, Espec_{EC} - ACAI (Tabela 3.23)
print(waldTest(obj=RM.EC.acai.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - ACAI (Tabela 3.23)
print(waldTest(obj=RM.EC.acai.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - ACAI (Tabela 3.23)
print(waldTest(obj=RM.EC.acai.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - ACAI (Tabela 3.23)

RM.EC.acai.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.EC.acai,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))
print(RM.EC.acai.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{EC}, VPN_{RM}, VPN_{EC} - ACAI (Tabela 3.23)
print(waldTest(obj=RM.EC.acai.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - ACAI (Tabela 3.23)
print(waldTest(obj=RM.EC.acai.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - ACAI (Tabela 3.23)
print(waldTest(obj=RM.EC.acai.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - ACAI (Tabela 3.23)



#Análises de RM e EC (Tabela 3.24)
RM.EC.mar.TF<-readCatdata(TF=c(6,1,1,2, 0,1,0,2, 51,22,5,13, 3,5,3,6, 53,45),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))

#sob MCAR
RM.EC.mcar<-satMarML(RM.EC.mar.TF,,missing="MCAR",maxit=1000)

RM.EC.mcar.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.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))
print(RM.EC.mcar.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{EC}, Espec_{RM}, Espec_{EC} - MCAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mcar.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - MCAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mcar.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - MCAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mcar.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - MCAR (Tabela 3.24)

RM.EC.mcar.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.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))
print(RM.EC.mcar.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{EC}, VPN_{RM}, VPN_{EC} - MCAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mcar.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - MCAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mcar.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - MCAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mcar.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - MCAR (Tabela 3.24)


#sob MAR
RM.EC.mar<-satMarML(RM.EC.mar.TF,maxit=1000)

RM.EC.mar.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.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))
print(RM.EC.mar.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{EC}, Espec_{RM}, Espec_{EC} - MAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mar.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - MAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mar.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - MAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mar.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - MAR (Tabela 3.24)

RM.EC.mar.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.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))
print(RM.EC.mar.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{EC}, VPN_{RM}, VPN_{EC} - MAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mar.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - MAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mar.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - MAR (Tabela 3.24)
print(waldTest(obj=RM.EC.mar.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - MAR (Tabela 3.24)


                         
#Análises de RM e EC - substituindo freqüências nulas por 0.001 (Tabela 3.25)
RM.EC.mar.subs.TF<-readCatdata(TF=c(6,1,1,2, 0.001,1,0.001,2, 51,22,5,13, 3,5,3,6, 53,45),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))

#sob MCAR
RM.EC.mcar.subs<-satMarML(RM.EC.mar.subs.TF,missing="MCAR",maxit=1000)

RM.EC.mcar.subs.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.EC.mcar.subs,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))
print(RM.EC.mcar.subs.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{EC}, Espec_{RM}, Espec_{EC} - MCAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mcar.subs.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - MCAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mcar.subs.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - MCAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mcar.subs.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - MCAR (Tabela 3.25)
round(cov2cor(RM.EC.mcar.subs.sens$VFu),3) #Correlações entre as estimativas dos parâmetros: Sens_{RM}, Sens_{EC}, Espec_{RM}, Espec_{EC} - MCAR (Tabela 3.26)

RM.EC.mcar.subs.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.EC.mcar.subs,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))
print(RM.EC.mcar.subs.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{EC}, VPN_{RM}, VPN_{EC} - MCAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mcar.subs.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - MCAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mcar.subs.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - MCAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mcar.subs.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - MCAR (Tabela 3.25)
round(cov2cor(RM.EC.mcar.subs.vpp$VFu),3) #Correlações entre as estimativas dos parâmetros: VPP_{RM}, VPP_{EC}, VPN_{RM}, VPN_{EC} - MCAR (Tabela 3.26)

CC<-rbind(c(1,-1,0,0),c(0,0,1,-1))
round(cov2cor(CC%*%RM.EC.mcar.subs.vpp$VFu%*%t(CC)),3) #Correlações entre os contrastes lineares mencionados no parágrafo entre as Tabelas 3.25 e 3.26 - MCAR


#sob MAR
RM.EC.mar.subs<-satMarML(RM.EC.mar.subs.TF,maxit=1000)

RM.EC.mar.subs.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.EC.mar.subs,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))
print(RM.EC.mar.subs.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{EC}, Espec_{RM}, Espec_{EC} - MAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mar.subs.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - MAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mar.subs.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - MAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mar.subs.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - MAR (Tabela 3.25)
round(cov2cor(RM.EC.mar.subs.sens$VFu),3) #Correlações entre as estimativas dos parâmetros: Sens_{RM}, Sens_{EC}, Espec_{RM}, Espec_{EC} - MAR (Tabela 3.26)

RM.EC.mar.subs.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.EC.mar.subs,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))
print(RM.EC.mar.subs.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{EC}, VPN_{RM}, VPN_{EC} - MAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mar.subs.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - MAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mar.subs.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - MAR (Tabela 3.25)
print(waldTest(obj=RM.EC.mar.subs.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - MAR (Tabela 3.25)
round(cov2cor(RM.EC.mar.subs.vpp$VFu),3) #Correlações entre as estimativas dos parâmetros: VPP_{RM}, VPP_{EC}, VPN_{RM}, VPN_{EC} - MAR (Tabela 3.26)

CC<-rbind(c(1,-1,0,0),c(0,0,1,-1))
round(cov2cor(CC%*%RM.EC.mar.subs.vpp$VFu%*%t(CC)),3) #Correlações entre os contrastes lineares mencionados no parágrafo entre as Tabelas 3.25 e 3.26 - MAR



#ACC de RM e US (Tabela 3.27)
RM.US.acc<-readCatdata(TF=c(58,1,0.001,24, 4,1,1,15))

RM.US.acc.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.US.acc,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))
print(RM.US.acc.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{US}, Espec_{RM}, Espec_{US} - ACC (Tabela 3.27)
print(waldTest(obj=RM.US.acc.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - ACC (Tabela 3.27)
print(waldTest(obj=RM.US.acc.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - ACC (Tabela 3.27)
print(waldTest(obj=RM.US.acc.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - ACC (Tabela 3.27)

RM.US.acc.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.US.acc,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))
print(RM.US.acc.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{US}, VPN_{RM}, VPN_{US} - ACC (Tabela 3.27)
print(waldTest(obj=RM.US.acc.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - ACC (Tabela 3.27)
print(waldTest(obj=RM.US.acc.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - ACC (Tabela 3.27)
print(waldTest(obj=RM.US.acc.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - ACC (Tabela 3.27)



#ACAI de RM e EC (Tabela 3.27)
RM.US.acai<-readCatdata(TF=rbind(c(58,25,5,16),c(119,6,3,91)))

RM.US.acai.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.US.acai,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))
print(RM.US.acai.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{US}, Espec_{RM}, Espec_{US} - ACAI (Tabela 3.27)
print(waldTest(obj=RM.US.acai.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - ACAI (Tabela 3.27)
print(waldTest(obj=RM.US.acai.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - ACAI (Tabela 3.27)
print(waldTest(obj=RM.US.acai.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - ACAI (Tabela 3.27)

RM.US.acai.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.US.acai,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))
print(RM.US.acai.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{US}, VPN_{RM}, VPN_{US} - ACAI (Tabela 3.27)
print(waldTest(obj=RM.US.acai.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - ACAI (Tabela 3.27)
print(waldTest(obj=RM.US.acai.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - ACAI (Tabela 3.27)
print(waldTest(obj=RM.US.acai.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - ACAI (Tabela 3.27)



#Análises de RM e US (Tabela 3.28)
RM.US.mar.TF<-readCatdata(TF=c(58,1,0.001,24, 4,1,1,15, 57,4,2,52),Zp=rep(1,2)%x%diag(4),Rp=4)

#sob MCAR
RM.US.mcar<-satMarML(RM.US.mar.TF,missing="MCAR",maxit=1000)

RM.US.mcar.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.US.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))
print(RM.US.mcar.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{US}, Espec_{RM}, Espec_{US} - MCAR (Tabela 3.28)
print(waldTest(obj=RM.US.mcar.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - MCAR (Tabela 3.28)
print(waldTest(obj=RM.US.mcar.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - MCAR (Tabela 3.28)
print(waldTest(obj=RM.US.mcar.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - MCAR (Tabela 3.28)

RM.US.mcar.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.US.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))
print(RM.US.mcar.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{US}, VPN_{RM}, VPN_{US} - MCAR (Tabela 3.28)
print(waldTest(obj=RM.US.mcar.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - MCAR (Tabela 3.28)
print(waldTest(obj=RM.US.mcar.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - MCAR (Tabela 3.28)
print(waldTest(obj=RM.US.mcar.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - MCAR (Tabela 3.28)


#sob MAR
RM.US.mar<-satMarML(RM.US.mar.TF,maxit=1000)

RM.US.mar.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.US.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))
print(RM.US.mar.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{US}, Espec_{RM}, Espec_{US} - MAR (Tabela 3.28)
print(waldTest(obj=RM.US.mar.sens,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_S - MAR (Tabela 3.28)
print(waldTest(obj=RM.US.mar.sens,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_E - MAR (Tabela 3.28)
print(waldTest(obj=RM.US.mar.sens,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{SE} - MAR (Tabela 3.28)

RM.US.mar.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.US.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))
print(RM.US.mar.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{US}, VPN_{RM}, VPN_{US} - MAR (Tabela 3.28)
print(waldTest(obj=RM.US.mar.vpp,C=c(1,-1,0,0)),digits=3) #Valor-P do teste H_P - MAR (Tabela 3.28)
print(waldTest(obj=RM.US.mar.vpp,C=c(0,0,1,-1)),digits=3) #Valor-P do teste H_N - MAR (Tabela 3.28)
print(waldTest(obj=RM.US.mar.vpp,C=rbind(c(1,-1,0,0),c(0,0,1,-1))),digits=3) #Valor-P do teste H_{PN} - MAR (Tabela 3.28)



#Análises de RM, EC e US sob MCAR (Tabela 3.29)
RM.EC.US.mcar.TF<-readCatdata(TF=c(6,0.001,1,0.001, 0.001,0.001,0.001,0.001, 0.001,1,0.001,2, 0.001,1,0.001,2,   51,1,4,1, 0.001,21,1,12,   3,1,3,1,0.001,4,0.001,5,   51,2,2,43),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))

RM.EC.US.mcar.mar<-satMarML(RM.EC.US.mcar.TF,missing="MCAR",maxit=1000)
print(RM.EC.US.mcar.mar,digits=3) #Valores-P dos testes de razão de verossimilhanças e de Pearson para o mecanismo MCAR condicionalmente ao MAR (primeiro parágrafo da p.134)


RM.EC.US.mcar.sens<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.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))
print(RM.EC.US.mcar.sens,digits=3) #est.e e.p. para Sens_{RM}, Sens_{EC}, Sens_{US}, Espec_{RM}, Espec_{EC}, Espec_{US} - MCAR (Tabela 3.29)
round(cov2cor(RM.EC.US.mcar.sens$VFu),3) #Correlações entre as estimativas dos parâmetros: Sens_{RM}, Sens_{EC}, Sens_{US}, Espec_{RM}, Espec_{EC}, Espec_{US} - MCAR (Tabela 3.29)
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=RM.EC.US.mcar.sens,C=C1) #Valor-P do teste Sens_{RM}=Sens_{EC}=Sens_{US} (não apresentado)
waldTest(obj=RM.EC.US.mcar.sens,C=C2) #Valor-P do teste Espec_{RM}=Espec_{EC}=Espec_{US} (não apresentado)
waldTest(obj=RM.EC.US.mcar.sens,C=rbind(C1,C2)) #Valor-P de ambos os testes simultaneamente (não apresentado)

print(waldTest(obj=RM.EC.US.mcar.sens,C=c(1,-1,0,0,0,0)),digits=3) #Valor-P do teste Sens_{RM}=Sens_{EC} (Tabela 3.29)
print(waldTest(obj=RM.EC.US.mcar.sens,C=c(1,0,-1,0,0,0)),digits=3) #Valor-P do teste Sens_{RM}=Sens_{US} (Tabela 3.29)
print(waldTest(obj=RM.EC.US.mcar.sens,C=c(0,1,-1,0,0,0)),digits=3) #Valor-P do teste Sens_{EC}=Sens_{US} (Tabela 3.29)

print(waldTest(obj=RM.EC.US.mcar.sens,C=c(0,0,0,1,-1,0)),digits=3) #Valor-P do teste Espec_{RM}=Espec_{EC} (Tabela 3.29)
print(waldTest(obj=RM.EC.US.mcar.sens,C=c(0,0,0,1,0,-1)),digits=3) #Valor-P do teste Espec_{RM}=Espec_{US} (Tabela 3.29)
print(waldTest(obj=RM.EC.US.mcar.sens,C=c(0,0,0,0,1,-1)),digits=3) #Valor-P do teste Espec_{EC}=Espec_{US} (Tabela 3.29)


RM.EC.US.mcar.vpp<-funlinWLS(model=c("exp","lin","log","lin"),obj=RM.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))
print(RM.EC.US.mcar.vpp,digits=3) #est.e e.p. para VPP_{RM}, VPP_{EC}, VPP_{US}, VPN_{RM}, VPN_{EC}, VPN_{US} - MCAR (Tabela 3.29)
round(cov2cor(RM.EC.US.mcar.vpp$VFu),3) #Correlações entre as estimativas dos parâmetros: VPP_{RM}, VPP_{EC}, VPP_{US}, VPN_{RM}, VPN_{EC}, VPN_{US} - MCAR (Tabela 3.29)

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=RM.EC.US.mcar.vpp,C=C1) #Valor-P do teste VPP_{RM}=VPP_{EC}=VPP_{US} (não apresentado)
waldTest(obj=RM.EC.US.mcar.vpp,C=C2) #Valor-P do teste VPN_{RM}=VPN_{EC}=VPN_{US} (não apresentado)
waldTest(obj=RM.EC.US.mcar.vpp,C=rbind(C1,C2)) #Valor-P de ambos os testes simultaneamente (não apresentado)

print(waldTest(obj=RM.EC.US.mcar.vpp,C=c(1,-1,0,0,0,0)),digits=3) #Valor-P do teste VPP_{RM}=VPP_{EC} (Tabela 3.29)
print(waldTest(obj=RM.EC.US.mcar.vpp,C=c(1,0,-1,0,0,0)),digits=3) #Valor-P do teste VPP_{RM}=VPP_{US} (Tabela 3.29)
print(waldTest(obj=RM.EC.US.mcar.vpp,C=c(0,1,-1,0,0,0)),digits=3) #Valor-P do teste VPP_{EC}=VPP_{US} (Tabela 3.29)

print(waldTest(obj=RM.EC.US.mcar.vpp,C=c(0,0,0,1,-1,0)),digits=3) #Valor-P do teste VPN_{RM}=VPN_{EC} (Tabela 3.29)
print(waldTest(obj=RM.EC.US.mcar.vpp,C=c(0,0,0,1,0,-1)),digits=3) #Valor-P do teste VPN_{RM}=VPN_{US} (Tabela 3.29)
print(waldTest(obj=RM.EC.US.mcar.vpp,C=c(0,0,0,0,1,-1)),digits=3) #Valor-P do teste VPN_{EC}=VPN_{US} (Tabela 3.29)

waldTest(obj=RM.EC.US.mcar.vpp,C=rbind(c(1,-1,0,0,0,0),c(0,0,0,1,-1,0))) #Valor-P do teste VPP_{RM}=VPP_{EC} e VPN_{RM}=VPN_{EC} (não apresentado)