#Commands to reproduce the analyses presented in
#Poleto, F. Z., Singer, J. M. and Paulino, C. D. (2012). A product-multinomial framework for categorical data analysis with missing responses. Submitted for publication.

#The functions are available at http://www.poleto.com/missing.html, or you can load directly using the command
source("http://www.poleto.com/Catdata.r")



#Example 1
smoking.TF<-rbind(c(167,17,19,10,1,3,52,10,11, 176,24,121, 28,10,12),
                  c(120,22,19, 8,5,1,39,12,12, 103, 3, 80, 31, 8,14)) #Table of frequencies (Table 1)
smoking.Zp<-t(rep(1,2))%x%cbind(diag(3)%x%rep(1,3), rep(1,3)%x%diag(3)) #[Z_{st},t=2,...,T] matrices combined column-wise for s=1,...,S
smoking.Rp<-rbind(c(3,3),c(3,3)) #[R_{st},t=2,...,T] vectors combined row-wise for s=1,...,S
smoking.catdata<-readCatdata(TF=smoking.TF,Zp=smoking.Zp,Rp=smoking.Rp) #Reading the categorical data with missingness
smoking.catdata #Proportions of the complete data

smoking.satmarml<-satMarML(smoking.catdata) #Application of the theory of Section 4.1
smoking.satmarml
smoking.satmcarwls<-satMcarWLS(smoking.catdata) #Application of the theory of Section 4.2
smoking.satmcarwls

smoking.E<-rbind(c(1,-1,0),c(0,1,-1))
smoking.A<-diag(2)%x%smoking.E%x%smoking.E #Equal to matrix A detailed in (6.1)
smoking.loglin2.ml<-loglinML(smoking.satmarml,A=smoking.A,XL=rep(1,8)) #Fit log-linear model (6.1) by ML methodology
smoking.loglin2.ml
smoking.loglin2.hybrid<-funlinWLS(model=c("lin","log"),obj=smoking.satmarml,A1=smoking.A,XL=rep(1,8)) #Fit log-linear model (6.1) by the hybrid approach
smoking.loglin2.hybrid

smoking.loglin0.ml<-loglinML(smoking.satmarml,U=smoking.A) #Fit log-linear independence model by ML methodology
smoking.loglin0.ml
z2<-smoking.Zp[,1:3]
z3<-smoking.Zp[,4:6]
#Tests of independence conditionally on model (6.1)
#Likelihood ratio
smoking.loglin0.ml$QvH-smoking.loglin2.ml$QvH
1-pchisq(smoking.loglin0.ml$QvH-smoking.loglin2.ml$QvH,smoking.loglin0.ml$glH-smoking.loglin2.ml$glH)
#Pearson
sum(((smoking.loglin2.ml$ystH$st1.1-smoking.loglin0.ml$ystH$st1.1)^2)/smoking.loglin0.ml$ystH$st1.1)+sum(((smoking.loglin2.ml$ystH$st2.1-smoking.loglin0.ml$ystH$st2.1)^2)/smoking.loglin0.ml$ystH$st2.1)+
  sum(((t(z2)%*%(smoking.loglin2.ml$ystH$st1.2-smoking.loglin0.ml$ystH$st1.2))^2)/(t(z2)%*%smoking.loglin0.ml$ystH$st1.2))+sum(((t(z2)%*%(smoking.loglin2.ml$ystH$st2.2-smoking.loglin0.ml$ystH$st2.2))^2)/(t(z2)%*%smoking.loglin0.ml$ystH$st2.2))+
  sum(((t(z3)%*%(smoking.loglin2.ml$ystH$st1.3-smoking.loglin0.ml$ystH$st1.3))^2)/(t(z3)%*%smoking.loglin0.ml$ystH$st1.3))+sum(((t(z3)%*%(smoking.loglin2.ml$ystH$st2.3-smoking.loglin0.ml$ystH$st2.3))^2)/(t(z3)%*%smoking.loglin0.ml$ystH$st2.3))
1-pchisq(sum(((smoking.loglin2.ml$ystH$st1.1-smoking.loglin0.ml$ystH$st1.1)^2)/smoking.loglin0.ml$ystH$st1.1)+sum(((smoking.loglin2.ml$ystH$st2.1-smoking.loglin0.ml$ystH$st2.1)^2)/smoking.loglin0.ml$ystH$st2.1)+
  sum(((t(z2)%*%(smoking.loglin2.ml$ystH$st1.2-smoking.loglin0.ml$ystH$st1.2))^2)/(t(z2)%*%smoking.loglin0.ml$ystH$st1.2))+sum(((t(z2)%*%(smoking.loglin2.ml$ystH$st2.2-smoking.loglin0.ml$ystH$st2.2))^2)/(t(z2)%*%smoking.loglin0.ml$ystH$st2.2))+
  sum(((t(z3)%*%(smoking.loglin2.ml$ystH$st1.3-smoking.loglin0.ml$ystH$st1.3))^2)/(t(z3)%*%smoking.loglin0.ml$ystH$st1.3))+sum(((t(z3)%*%(smoking.loglin2.ml$ystH$st2.3-smoking.loglin0.ml$ystH$st2.3))^2)/(t(z3)%*%smoking.loglin0.ml$ystH$st2.3))
  ,smoking.loglin0.ml$glH-smoking.loglin2.ml$glH)
#Neyman
sum(((smoking.loglin2.ml$ystH$st1.1-smoking.loglin0.ml$ystH$st1.1)^2)/smoking.loglin2.ml$ystH$st1.1)+sum(((smoking.loglin2.ml$ystH$st2.1-smoking.loglin0.ml$ystH$st2.1)^2)/smoking.loglin2.ml$ystH$st2.1)+
  sum(((t(z2)%*%(smoking.loglin2.ml$ystH$st1.2-smoking.loglin0.ml$ystH$st1.2))^2)/(t(z2)%*%smoking.loglin2.ml$ystH$st1.2))+sum(((t(z2)%*%(smoking.loglin2.ml$ystH$st2.2-smoking.loglin0.ml$ystH$st2.2))^2)/(t(z2)%*%smoking.loglin2.ml$ystH$st2.2))+
  sum(((t(z3)%*%(smoking.loglin2.ml$ystH$st1.3-smoking.loglin0.ml$ystH$st1.3))^2)/(t(z3)%*%smoking.loglin2.ml$ystH$st1.3))+sum(((t(z3)%*%(smoking.loglin2.ml$ystH$st2.3-smoking.loglin0.ml$ystH$st2.3))^2)/(t(z3)%*%smoking.loglin2.ml$ystH$st2.3))
1-pchisq(sum(((smoking.loglin2.ml$ystH$st1.1-smoking.loglin0.ml$ystH$st1.1)^2)/smoking.loglin2.ml$ystH$st1.1)+sum(((smoking.loglin2.ml$ystH$st2.1-smoking.loglin0.ml$ystH$st2.1)^2)/smoking.loglin2.ml$ystH$st2.1)+
  sum(((t(z2)%*%(smoking.loglin2.ml$ystH$st1.2-smoking.loglin0.ml$ystH$st1.2))^2)/(t(z2)%*%smoking.loglin2.ml$ystH$st1.2))+sum(((t(z2)%*%(smoking.loglin2.ml$ystH$st2.2-smoking.loglin0.ml$ystH$st2.2))^2)/(t(z2)%*%smoking.loglin2.ml$ystH$st2.2))+
  sum(((t(z3)%*%(smoking.loglin2.ml$ystH$st1.3-smoking.loglin0.ml$ystH$st1.3))^2)/(t(z3)%*%smoking.loglin2.ml$ystH$st1.3))+sum(((t(z3)%*%(smoking.loglin2.ml$ystH$st2.3-smoking.loglin0.ml$ystH$st2.3))^2)/(t(z3)%*%smoking.loglin2.ml$ystH$st2.3))
  ,smoking.loglin0.ml$glH-smoking.loglin2.ml$glH)
#Wald, ML
waldTest(smoking.loglin2.ml,C=1)
#Wald, hybrid
waldTest(smoking.loglin2.hybrid,C=1)

smoking.loglin1.ml<-loglinML(smoking.satmarml,A=smoking.A,XL=rep(1,2)%x%diag(4)) #Fit log-linear model (6.2) by ML methodology
smoking.loglin1.ml
smoking.loglin1.hybrid<-funlinWLS(model=c("lin","log"),obj=smoking.satmarml,A1=smoking.A,XL=rep(1,2)%x%diag(4)) #Fit log-linear model (6.2) by the hybrid approach
smoking.loglin1.hybrid

#Tests of independence conditionally on model (6.2)
#Likelihood ratio
smoking.loglin0.ml$QvH-smoking.loglin1.ml$QvH
1-pchisq(smoking.loglin0.ml$QvH-smoking.loglin1.ml$QvH,smoking.loglin0.ml$glH-smoking.loglin1.ml$glH)
#Pearson
sum(((smoking.loglin1.ml$ystH$st1.1-smoking.loglin0.ml$ystH$st1.1)^2)/smoking.loglin0.ml$ystH$st1.1)+sum(((smoking.loglin1.ml$ystH$st2.1-smoking.loglin0.ml$ystH$st2.1)^2)/smoking.loglin0.ml$ystH$st2.1)+
  sum(((t(z2)%*%(smoking.loglin1.ml$ystH$st1.2-smoking.loglin0.ml$ystH$st1.2))^2)/(t(z2)%*%smoking.loglin0.ml$ystH$st1.2))+sum(((t(z2)%*%(smoking.loglin1.ml$ystH$st2.2-smoking.loglin0.ml$ystH$st2.2))^2)/(t(z2)%*%smoking.loglin0.ml$ystH$st2.2))+
  sum(((t(z3)%*%(smoking.loglin1.ml$ystH$st1.3-smoking.loglin0.ml$ystH$st1.3))^2)/(t(z3)%*%smoking.loglin0.ml$ystH$st1.3))+sum(((t(z3)%*%(smoking.loglin1.ml$ystH$st2.3-smoking.loglin0.ml$ystH$st2.3))^2)/(t(z3)%*%smoking.loglin0.ml$ystH$st2.3))
1-pchisq(sum(((smoking.loglin1.ml$ystH$st1.1-smoking.loglin0.ml$ystH$st1.1)^2)/smoking.loglin0.ml$ystH$st1.1)+sum(((smoking.loglin1.ml$ystH$st2.1-smoking.loglin0.ml$ystH$st2.1)^2)/smoking.loglin0.ml$ystH$st2.1)+
  sum(((t(z2)%*%(smoking.loglin1.ml$ystH$st1.2-smoking.loglin0.ml$ystH$st1.2))^2)/(t(z2)%*%smoking.loglin0.ml$ystH$st1.2))+sum(((t(z2)%*%(smoking.loglin1.ml$ystH$st2.2-smoking.loglin0.ml$ystH$st2.2))^2)/(t(z2)%*%smoking.loglin0.ml$ystH$st2.2))+
  sum(((t(z3)%*%(smoking.loglin1.ml$ystH$st1.3-smoking.loglin0.ml$ystH$st1.3))^2)/(t(z3)%*%smoking.loglin0.ml$ystH$st1.3))+sum(((t(z3)%*%(smoking.loglin1.ml$ystH$st2.3-smoking.loglin0.ml$ystH$st2.3))^2)/(t(z3)%*%smoking.loglin0.ml$ystH$st2.3))
  ,smoking.loglin0.ml$glH-smoking.loglin1.ml$glH)
#Neyman
sum(((smoking.loglin1.ml$ystH$st1.1-smoking.loglin0.ml$ystH$st1.1)^2)/smoking.loglin1.ml$ystH$st1.1)+sum(((smoking.loglin1.ml$ystH$st2.1-smoking.loglin0.ml$ystH$st2.1)^2)/smoking.loglin1.ml$ystH$st2.1)+
  sum(((t(z2)%*%(smoking.loglin1.ml$ystH$st1.2-smoking.loglin0.ml$ystH$st1.2))^2)/(t(z2)%*%smoking.loglin1.ml$ystH$st1.2))+sum(((t(z2)%*%(smoking.loglin1.ml$ystH$st2.2-smoking.loglin0.ml$ystH$st2.2))^2)/(t(z2)%*%smoking.loglin1.ml$ystH$st2.2))+
  sum(((t(z3)%*%(smoking.loglin1.ml$ystH$st1.3-smoking.loglin0.ml$ystH$st1.3))^2)/(t(z3)%*%smoking.loglin1.ml$ystH$st1.3))+sum(((t(z3)%*%(smoking.loglin1.ml$ystH$st2.3-smoking.loglin0.ml$ystH$st2.3))^2)/(t(z3)%*%smoking.loglin1.ml$ystH$st2.3))
1-pchisq(sum(((smoking.loglin1.ml$ystH$st1.1-smoking.loglin0.ml$ystH$st1.1)^2)/smoking.loglin1.ml$ystH$st1.1)+sum(((smoking.loglin1.ml$ystH$st2.1-smoking.loglin0.ml$ystH$st2.1)^2)/smoking.loglin1.ml$ystH$st2.1)+
  sum(((t(z2)%*%(smoking.loglin1.ml$ystH$st1.2-smoking.loglin0.ml$ystH$st1.2))^2)/(t(z2)%*%smoking.loglin1.ml$ystH$st1.2))+sum(((t(z2)%*%(smoking.loglin1.ml$ystH$st2.2-smoking.loglin0.ml$ystH$st2.2))^2)/(t(z2)%*%smoking.loglin1.ml$ystH$st2.2))+
  sum(((t(z3)%*%(smoking.loglin1.ml$ystH$st1.3-smoking.loglin0.ml$ystH$st1.3))^2)/(t(z3)%*%smoking.loglin1.ml$ystH$st1.3))+sum(((t(z3)%*%(smoking.loglin1.ml$ystH$st2.3-smoking.loglin0.ml$ystH$st2.3))^2)/(t(z3)%*%smoking.loglin1.ml$ystH$st2.3))
  ,smoking.loglin0.ml$glH-smoking.loglin1.ml$glH)
#Wald, ML
waldTest(smoking.loglin1.ml,C=diag(4))
#Wald, hybrid
waldTest(smoking.loglin1.hybrid,C=diag(4))

#Minus log-likelihood of the MNAR described in the last paragraph of Section 3
mnar.mll<-function(p,fr){
	# p=(pi11(1),...,pi33(2),a2(11),a2(21),a2(31),a3(11),a3(21),a3(31),a2(12),a2(22),a2(32),a3(12),a3(22),a3(32))
	pi11.1<-p[1]; pi12.1<-p[2]; pi13.1<-p[3]
	pi21.1<-p[4]; pi22.1<-p[5]; pi23.1<-p[6]
	pi31.1<-p[7]; pi32.1<-p[8]; pi33.1=1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1

	pi11.2<-p[9]; pi12.2<-p[10];pi13.2<-p[11]
	pi21.2<-p[12];pi22.2<-p[13];pi23.2<-p[14]
	pi31.2<-p[15];pi32.2<-p[16]; pi33.2=1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2

	a2.1.1<-p[17];a2.2.1<-p[18];a2.3.1<-p[19]
	a3.1.1<-p[20];a3.2.1<-p[21];a3.3.1<-p[22]

	a2.1.2<-p[23];a2.2.2<-p[24];a2.3.2<-p[25]
	a3.1.2<-p[26];a3.2.2<-p[27];a3.3.2<-p[28]

	value<- -(
	 fr[1,1]*log(pi11.1*(1-a2.1.1-a3.1.1))+fr[1,2]*log(pi12.1*(1-a2.2.1-a3.1.1))+fr[1,3]*log(pi13.1*(1-a2.3.1-a3.1.1))+
	 fr[1,4]*log(pi21.1*(1-a2.1.1-a3.2.1))+fr[1,5]*log(pi22.1*(1-a2.2.1-a3.2.1))+fr[1,6]*log(pi23.1*(1-a2.3.1-a3.2.1))+
	 fr[1,7]*log(pi31.1*(1-a2.1.1-a3.3.1))+fr[1,8]*log(pi32.1*(1-a2.2.1-a3.3.1))+fr[1,9]*log(pi33.1*(1-a2.3.1-a3.3.1))+

	 fr[1,10]*log(pi11.1*a2.1.1 + pi12.1*a2.2.1 + pi13.1*a2.3.1)+
	 fr[1,11]*log(pi21.1*a2.1.1 + pi22.1*a2.2.1 + pi23.1*a2.3.1)+
	 fr[1,12]*log(pi31.1*a2.1.1 + pi32.1*a2.2.1 + pi33.1*a2.3.1)+

	 fr[1,13]*log(pi11.1*a3.1.1 + pi21.1*a3.2.1 + pi31.1*a3.3.1)+
	 fr[1,14]*log(pi12.1*a3.1.1 + pi22.1*a3.2.1 + pi32.1*a3.3.1)+
	 fr[1,15]*log(pi13.1*a3.1.1 + pi23.1*a3.2.1 + pi33.1*a3.3.1)+

	 fr[2,1]*log(pi11.2*(1-a2.1.2-a3.1.2))+fr[2,2]*log(pi12.2*(1-a2.2.2-a3.1.2))+fr[2,3]*log(pi13.2*(1-a2.3.2-a3.1.2))+
	 fr[2,4]*log(pi21.2*(1-a2.1.2-a3.2.2))+fr[2,5]*log(pi22.2*(1-a2.2.2-a3.2.2))+fr[2,6]*log(pi23.2*(1-a2.3.2-a3.2.2))+
	 fr[2,7]*log(pi31.2*(1-a2.1.2-a3.3.2))+fr[2,8]*log(pi32.2*(1-a2.2.2-a3.3.2))+fr[2,9]*log(pi33.2*(1-a2.3.2-a3.3.2))+

	 fr[2,10]*log(pi11.2*a2.1.2 + pi12.2*a2.2.2 + pi13.2*a2.3.2)+
	 fr[2,11]*log(pi21.2*a2.1.2 + pi22.2*a2.2.2 + pi23.2*a2.3.2)+
	 fr[2,12]*log(pi31.2*a2.1.2 + pi32.2*a2.2.2 + pi33.2*a2.3.2)+

	 fr[2,13]*log(pi11.2*a3.1.2 + pi21.2*a3.2.2 + pi31.2*a3.3.2)+
	 fr[2,14]*log(pi12.2*a3.1.2 + pi22.2*a3.2.2 + pi32.2*a3.3.2)+
	 fr[2,15]*log(pi13.2*a3.1.2 + pi23.2*a3.2.2 + pi33.2*a3.3.2)

	)
	value
}

mnar.fit<-constrOptim(theta=c(rep(1/9,16),rep(1/3,12)),f=mnar.mll,method="Nelder-Mead",ui=diag(28),ci=rep(0,28),
  control=list(maxit=10000),outer.iterations=1000,fr=smoking.TF)

#hessian matrix
mnar.der<-deriv3(~-(
	 o1.1*log(pi11.1*(1-a2.1.1-a3.1.1))+o1.2*log(pi12.1*(1-a2.2.1-a3.1.1))+o1.3*log(pi13.1*(1-a2.3.1-a3.1.1))+
	 o1.4*log(pi21.1*(1-a2.1.1-a3.2.1))+o1.5*log(pi22.1*(1-a2.2.1-a3.2.1))+o1.6*log(pi23.1*(1-a2.3.1-a3.2.1))+
	 o1.7*log(pi31.1*(1-a2.1.1-a3.3.1))+o1.8*log(pi32.1*(1-a2.2.1-a3.3.1))+o1.9*log((1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(1-a2.3.1-a3.3.1))+

	 o1.10*log(pi11.1*a2.1.1 + pi12.1*a2.2.1 + pi13.1*a2.3.1)+
	 o1.11*log(pi21.1*a2.1.1 + pi22.1*a2.2.1 + pi23.1*a2.3.1)+
	 o1.12*log(pi31.1*a2.1.1 + pi32.1*a2.2.1 + (1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*a2.3.1)+

	 o1.13*log(pi11.1*a3.1.1 + pi21.1*a3.2.1 + pi31.1*a3.3.1)+
	 o1.14*log(pi12.1*a3.1.1 + pi22.1*a3.2.1 + pi32.1*a3.3.1)+
	 o1.15*log(pi13.1*a3.1.1 + pi23.1*a3.2.1 + (1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*a3.3.1)+

	 o2.1*log(pi11.2*(1-a2.1.2-a3.1.2))+o2.2*log(pi12.2*(1-a2.2.2-a3.1.2))+o2.3*log(pi13.2*(1-a2.3.2-a3.1.2))+
	 o2.4*log(pi21.2*(1-a2.1.2-a3.2.2))+o2.5*log(pi22.2*(1-a2.2.2-a3.2.2))+o2.6*log(pi23.2*(1-a2.3.2-a3.2.2))+
	 o2.7*log(pi31.2*(1-a2.1.2-a3.3.2))+o2.8*log(pi32.2*(1-a2.2.2-a3.3.2))+o2.9*log((1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(1-a2.3.2-a3.3.2))+

	 o2.10*log(pi11.2*a2.1.2 + pi12.2*a2.2.2 + pi13.2*a2.3.2)+
	 o2.11*log(pi21.2*a2.1.2 + pi22.2*a2.2.2 + pi23.2*a2.3.2)+
	 o2.12*log(pi31.2*a2.1.2 + pi32.2*a2.2.2 + (1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*a2.3.2)+

	 o2.13*log(pi11.2*a3.1.2 + pi21.2*a3.2.2 + pi31.2*a3.3.2)+
	 o2.14*log(pi12.2*a3.1.2 + pi22.2*a3.2.2 + pi32.2*a3.3.2)+
	 o2.15*log(pi13.2*a3.1.2 + pi23.2*a3.2.2 + (1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*a3.3.2)

	),c("pi11.1","pi12.1","pi13.1","pi21.1","pi22.1","pi23.1","pi31.1","pi32.1",
	 "pi11.2","pi12.2","pi13.2","pi21.2","pi22.2","pi23.2","pi31.2","pi32.2",
	 "a2.1.1","a2.2.1","a2.3.1","a3.1.1","a3.2.1","a3.3.1","a2.1.2","a2.2.2","a2.3.2","a3.1.2","a3.2.2","a3.3.2"),
	c("pi11.1","pi12.1","pi13.1","pi21.1","pi22.1","pi23.1","pi31.1","pi32.1",
	 "pi11.2","pi12.2","pi13.2","pi21.2","pi22.2","pi23.2","pi31.2","pi32.2",
	 "a2.1.1","a2.2.1","a2.3.1","a3.1.1","a3.2.1","a3.3.1","a2.1.2","a2.2.2","a2.3.2","a3.1.2","a3.2.2","a3.3.2",
	 "o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8","o1.9","o1.10",
	 "o1.11","o1.12","o1.13","o1.14","o1.15",
	 "o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8","o2.9","o2.10",
	 "o2.11","o2.12","o2.13","o2.14","o2.15")
)

p<-mnar.fit$par
TF<-smoking.TF
mnar.InfObs<-mnar.der(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8],p[9],p[10],p[11],p[12],p[13],p[14],p[15],p[16],p[17],p[18],p[19],p[20],p[21],p[22],p[23],p[24],p[25],p[26],p[27],p[28],
	TF[1,1],TF[1,2],TF[1,3],TF[1,4],TF[1,5],TF[1,6],TF[1,7],TF[1,8],TF[1,9],TF[1,10],TF[1,11],TF[1,12],TF[1,13],TF[1,14],TF[1,15],
	TF[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8],TF[2,9],TF[2,10],TF[2,11],TF[2,12],TF[2,13],TF[2,14],TF[2,15])

b<-smoking.catdata$b #b in (8), i.e., rep(1,2)%x%c(rep(0,8),1)
B<-smoking.catdata$B #B in (8), i.e., diag(2)%x%rbind(diag(8),rep(-1,8))
smoking.loglin2mnar.hybrid<-funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar.fit$par[1:16])),
  Vtheta=B%*%solve(attr(mnar.InfObs,"hessian")[1,,])[1:16,1:16]%*%t(B),A1=smoking.A,X=rep(1,8)) #Fit log-linear model (6.1) under MNAR by the hybrid approach
smoking.loglin2mnar.hybrid
smoking.loglin1mnar.hybrid<-funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar.fit$par[1:16])),
  Vtheta=B%*%solve(attr(mnar.InfObs,"hessian")[1,,])[1:16,1:16]%*%t(B),A1=smoking.A,X=rep(1,2)%x%diag(4)) #Fit log-linear model (6.2) under MNAR by the hybrid approach
smoking.loglin1mnar.hybrid
waldTest(smoking.loglin1mnar.hybrid,C=diag(4))





#Example 2
obes.TF<-rbind(
  c( 90, 9, 3, 7, 0,1, 1, 8,16, 5,0, 0, 9,3,0,0,129,18,6,13,32, 5,33,11,70,24) ,
  c(150,15, 8, 8, 8,9, 7,20,38, 3,1,11,16,6,1,3, 42, 2,3,13,45, 7,33, 4,55,14) ,
  c(152,11, 8,10, 7,7, 9,25,48, 6,2,14,13,5,0,3, 36, 5,4, 3,59,17,31, 9,40, 9) ,
  c(119, 7, 8, 3,13,4,11,16,42, 4,4,13,14,2,1,4, 18, 3,3, 1,82,24,23, 6,37,14) ,
  c(101, 4, 2, 7, 8,0, 6,15,82, 9,8,12, 6,1,0,1, 13, 1,2, 2,95,23,34,12,15, 3) ,
  c( 75, 8, 2, 4, 2,2, 1, 8,20, 0,0, 4, 7,2,0,1,109,22,7,24,23, 5,27, 5,65,19) ,
  c(154,14,13,19, 2,6, 6,21,25, 3,1,11,16,3,0,4, 47, 4,1, 8,47, 7,23, 5,39,13) ,
  c(148, 6,10, 8,12,0, 8,27,36, 0,7,17, 8,1,1,4, 39, 6,7,13,53,16,25, 9,23, 8) ,
  c(129, 8, 7, 9, 6,2, 7,14,36, 9,4,13,31,4,2,6, 19, 1,2, 2,58,37,21, 1,23,10) ,
  c( 91, 9, 5, 3, 6,0, 6,15,83,15,6,23, 5,0,0,1, 11, 1,2, 3,89,32,43,15,14, 5) ) #Table of frequencies (Table 2)
obesR.TF<-obes.TF
obesR.TF[obesR.TF==0]<-1e-6 #Replacing null frequencies by 10^{-6}
obes.Zp<-t(rep(1,10))%x%cbind( diag(4)%x%rep(1,2) ,  diag(2)%x%rep(1,2)%x%diag(2) , rep(1,2)%x%diag(4) ,
  diag(2)%x%rep(1,4) ,  rep(1,2)%x%diag(2)%x%rep(1,2) , rep(1,4)%x%diag(2) ) #[Z_{st},t=2,...,T] matrices combined column-wise for s=1,...,S
obes.Rp<-rep(1,10)%x%t(c(4,4,4,2,2,2)) #[R_{st},t=2,...,T] vectors combined row-wise for s=1,...,S

sum(obes.TF<5)/sum(obes.TF>=0) #Proportion of frequencies less than 5
sum(obes.TF==0) #Number o null frequencies

obes.catdata<-readCatdata(TF=obes.TF,Zp=obes.Zp,Rp=obes.Rp) #Reading the categorical data with missingness
obesR.catdata<-readCatdata(TF=obesR.TF,Zp=obes.Zp,Rp=obes.Rp) #Again, with null frequencies replaced by 10^{-6}
obes.catdata #Proportions of the complete data
obesR.catdata

obes.mar<-satMarML(obes.catdata) #Application of the theory of Section 4.1
obesR.mar<-satMarML(obesR.catdata) #Again, with null frequencies replaced by 10^{-6}
obes.mar
obesR.mar

obes.mcar<-satMcarWLS(obes.catdata) #Application of the theory of Section 4.2
obesR2.mcar<-satMcarWLS(obes.catdata,zeroN=1/8) #Again, with null frequencies replaced by 10^{-6}
obes.mcar
obesR2.mcar

obes.A.marg<-diag(10)%x%t(cbind( diag(2)%x%rep(1,4) ,  rep(1,2)%x%diag(2)%x%rep(1,2) , rep(1,4)%x%diag(2) ))[c(2,4,6),]
obes.X1<-diag(2)%x%cbind(
  c(1,rep(0,14)),
  c(0,1,0,1,rep(0,11)),
  c(0,0,1,0,1,0,1,rep(0,8)),
  c(rep(0,5),1,0,1,0,1,rep(0,5)),
  c(rep(0,8),1,0,1,0,1,0,0),
  c(rep(0,11),1,0,1,0),
  c(rep(0,14),1)) #Model specification matrix (6.3)
obes.age<-c(6,8,10,8,10,12,10,12,14,12,14,16,14,16,18)
obes.X2<-diag(2)%x%cbind(rep(1,15),obes.age,obes.age^2) #Model specification matrix (6.4)
obes.piece<-obes.age-10
obes.piece[obes.age>10]<-0
obes.X3<-cbind(rep(1,2)%x%cbind(rep(1,15),obes.piece),c(0,1)%x%rep(1,15)) #Model specification matrix (6.5)

obes.lin1.ml<-linML(obes.mar,A=obes.A.marg,X=obes.X1) #Fit linear model (6.3) by ML methodology
obes.lin2.ml<-linML(obes.mar,A=obes.A.marg,X=obes.X2) #Fit linear model (6.4) by ML methodology
obes.lin3.ml<-linML(obes.mar,A=obes.A.marg,X=obes.X3) #Fit linear model (6.5) by ML methodology

obesR.lin1.ml<-linML(obesR.mar,A=obes.A.marg,X=obes.X1) #Fit linear model (6.3) by ML methodology (zeros replaced by 1e-6)
obesR.lin2.ml<-linML(obesR.mar,A=obes.A.marg,X=obes.X2) #Fit linear model (6.4) by ML methodology (zeros replaced by 1e-6)
obesR.lin3.ml<-linML(obesR.mar,A=obes.A.marg,X=obes.X3) #Fit linear model (6.5) by ML methodology (zeros replaced by 1e-6)

obes.Xassoc<-cbind(c(0,1)%x%rep(1,4),rep(1,2)%x%c(0,1)%x%rep(1,2),rep(1,4)%x%c(0,1))
obes.Xassoc<-cbind(diag(10)%x%obes.Xassoc, rep(1,10)%x%cbind(obes.Xassoc[,1]*obes.Xassoc[,2]+obes.Xassoc[,2]*obes.Xassoc[,3],
  obes.Xassoc[,1]*obes.Xassoc[,3],obes.Xassoc[,1]*obes.Xassoc[,2]*obes.Xassoc[,3]) ) #Model specification matrix (6.6)
obesL.mar<-loglinML(obj=obes.mar,X=obes.Xassoc) #Fit log-linear model (6.6) by ML methodology
obesL.mar
waldTest(obesL.mar,C=cbind(matrix(0,2,31),diag(2))) #Wald test of Markovian dependence structure

obesL.lin2.hybrid<-funlinWLS(model="lin",theta=obesL.mar$thetaH,Vtheta=obesL.mar$VthetaH,A1=obes.A.marg,X=obes.X2)
obesL.lin2.hybrid #Fitting of linear model (6.4) by the hybrid approach subsequently to the log-linear model (6.6)

obesR.lin2.hybrid<-funlinWLS(model="lin",obj=obesR.mar,A1=obes.A.marg,X=obes.X2)
obesR.lin2.hybrid #Fitting of linear model (6.4) by the hybrid approach (zeros replaced by 1e-6)

obes.lin2.hybrid<-funlinWLS(model="lin",obj=obes.mar,A1=obes.A.marg,X=obes.X2)
obes.lin2.hybrid #Fitting of linear model (6.4) by the hybrid approach

obesL.lin1.hybrid<-funlinWLS(model="lin",theta=obesL.mar$thetaH,Vtheta=obesL.mar$VthetaH,A1=obes.A.marg,X=obes.X1)
obesL.lin1.hybrid #Fitting of linear model (6.3) by the hybrid approach subsequently to the log-linear model (6.6)

obesL.lin3.hybrid<-funlinWLS(model="lin",theta=obesL.mar$thetaH,Vtheta=obesL.mar$VthetaH,A1=obes.A.marg,X=obes.X3)
obesL.lin3.hybrid #Fitting of linear model (6.5) by the hybrid approach subsequently to the log-linear model (6.6)

waldTest(obesL.lin1.hybrid,t(c(1,-1))%x%diag(7)) #Wald test for gender effect under the model (6.3)-(6.6)
waldTest(obesL.lin3.hybrid,c(0,0,1)) #Wald test for gender effect under the model (6.5)-(6.6)
waldTest(obesL.lin2.hybrid,t(c(1,-1))%x%diag(3)) #Wald test for gender effect under the model (6.4)-(6.6)

waldTest(obesR.lin2.hybrid,t(c(1,-1))%x%diag(3)) #Wald test for gender effect under the model (6.4) (zeros replaced by 1e-6) 
waldTest(obes.lin2.hybrid,t(c(1,-1))%x%diag(3)) #Wald test for gender effect under the model (6.4)

#Fit log-linear model (6.6) by ML methodology after replacing zeros by 1e-6, and then repeat the previous fitting of linear models
obesRL.mar<-loglinML(obj=obesR.mar,X=obes.Xassoc)
obesRL.lin1.hybrid<-funlinWLS(model="lin",theta=obesRL.mar$thetaH,Vtheta=obesRL.mar$VthetaH,A1=obes.A.marg,X=obes.X1)
obesRL.lin1.hybrid
obesRL.lin2.hybrid<-funlinWLS(model="lin",theta=obesRL.mar$thetaH,Vtheta=obesRL.mar$VthetaH,A1=obes.A.marg,X=obes.X2)
obesRL.lin2.hybrid
obesRL.lin3.hybrid<-funlinWLS(model="lin",theta=obesRL.mar$thetaH,Vtheta=obesRL.mar$VthetaH,A1=obes.A.marg,X=obes.X3)
obesRL.lin3.hybrid
waldTest(obesRL.lin1.hybrid,t(c(1,-1))%x%diag(7))
waldTest(obesRL.lin2.hybrid,t(c(1,-1))%x%diag(3))
waldTest(obesRL.lin3.hybrid,c(0,0,1))