#Comandos para reproduzir as análises do Exemplo 2, parte 1 (dados originais), pp.94-101, 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.

#OBS.: Os comandos deste arquivo funcionam no R 2.6.2 com o pacote geoR 1.6-21 e sua dependência sp 0.9-25.
#Nas versões seguintes, a função .nlmP do geoR aparentemente não está mais disponível para ser utilizada e, então, esses comandos
#teriam que ser substituídos por outros como constrOptim, optim ou nlm.
#No Windows, basta baixar http://cran.r-project.org/bin/windows/base/old/2.6.2/R-2.6.2-win32.exe e, depois de instalá-lo e abrir o R,
#peça para instalar o pacote geoR (menu "Pacotes" > opção "Instalar pacote(s)..." > escolha algum espelho > OK > geoR > OK),
#que ele automaticamente instalará o geoR 1.6-21 e sua dependência sp 0.9-25.

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


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))
padroes<-t(rep(1,2))%x%cbind(diag(9),diag(3)%x%rep(1,3), rep(1,3)%x%diag(3))
classes<-rep(1,2)%x%t(c(9,3,3))

E<-rbind(c(1,-1,0),c(0,1,-1))
A<-diag(2)%x%(E%x%E)
UL<-diag(2)%x%cbind(rep(0,3),diag(3))


#Intervalos para o melhor-pior caso apresentados na Tabela 3.1 (veja as Tabelas B.1, B.2 e B.3 no Apêndice B.1)

#Alocação A, Limites superiores para (\omega_{11(1)},\omega_{22(1)},\omega_{11(2)},\omega_{22(2)})
log(c(371*35/(10*17),35*144/(10*3),254*16/(8*22),16*106/(12*1)))
#Alocação B, Limites inferiores para (\omega_{11(1)},\omega_{11(2)})
log(c(167*1/(62*203),120*5/(42*133)))
#Alocação C, Limites inferiores para (\omega_{22(1)},\omega_{22(2)})
log(c(1*11/(141*39),5*12/(100*18)))
#Alocação D, Limites inferiores para (\omega_{12(1)},\omega_{21(1)},\omega_{12(2)},\omega_{21(2)})
log(c(17*3/(35*207),10*10/(201*35),22*1/(16*136),8*12/(150*16)))
#Alocação E, Limites superiores para (\omega_{12(1)},\omega_{12(2)})
log(c(203*39/(1*19),133*18/(5*19)))
#Alocação F, Limites superiores para (\omega_{21(1)},\omega_{21(2)})
log(c(62*141/(52*1),42*100/(39*5)))


catdata.acc<-readCatdata(TF=TF[,1:9])
loglin.acc<-funlinWLS(model=c("lin","log"),obj=catdata.acc,A1=A,XL=rep(1,8))
summary(loglin.acc) #Logaritmos das razões de chances observadas e sob o modelo H:\omega_{ij(s)}=\omega, erros padrões e teste Wald de H, Tabela 3.5, ACC
funlinWLS(model=c("lin","log"),obj=catdata.acc,A1=A,UL=UL) #Teste de Wald de H:\omega_{12(s)}=\omega_{21(s)}=\omega_{22(s)}=0,s=1,2, p.101, ACC


catdata<-readCatdata(TF=TF,Zp=padroes[,-c(1:9,16:24)],Rp=classes[,-1])

mcar.ml<-satMarML(catdata,missing="MCAR")
mar.ml<-satMarML(catdata)
mar.ml #Estatística de razão de verossimilhanças, MCAR, Tabela 3.2
mcar.ml$yst #Tabela 3.3, MCAR, t=1,2,3
mar.ml$yst #Tabela 3.3, MAR, t=1,2,3

rbind(mcar.ml$yst$st1.1+mcar.ml$yst$st1.2+mcar.ml$yst$st1.3,mcar.ml$yst$st2.1+mcar.ml$yst$st2.2+mcar.ml$yst$st2.3) #Tabela 3.3, MCAR, Total
rbind(mar.ml$yst$st1.1+mar.ml$yst$st1.2+mar.ml$yst$st1.3,mar.ml$yst$st2.1+mar.ml$yst$st2.2+mar.ml$yst$st2.3) #Tabela 3.3, MAR, Total

loglin.mar<-funlinWLS(model=c("lin","log"),obj=mar.ml,A1=A,XL=rep(1,8))
summary(loglin.mar) #Logaritmos das razões de chances observadas e sob o modelo H:\omega_{ij(s)}=\omega, erros padrões e teste Wald de H, Tabela 3.5, MCAR/MAR*
#*Segundo explicado em (2.21), as matrizes (estimadas) de informação observada sob MCAR e MAR e de informação de Fisher sob MAR são todas iguais sob o modelo estrutural saturado

funlinWLS(model=c("lin","log"),obj=mar.ml,A1=A,UL=UL) #Teste de Wald de H:\omega_{12(s)}=\omega_{21(s)}=\omega_{22(s)}=0,s=1,2, p.101, MCAR/MAR



library(geoR)
sat.lv<-sum(TF*log(TF/rowSums(TF))) #log-verossimilhança sob modelo saturado

b<-kronecker(rep(1,2),c(rep(0,8),1))
B<-kronecker(diag(2),rbind(diag(8),rep(-1,8)))





#menos log-verossimilhança do modelo MNAR1, saturado
mnar1.mlv<-function(p,fr=TF){
	# 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
}

mnar1.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")
)

mnar1.fit<-.nlmP(objfunc=mnar1.mlv,params=c(rep(1/9,16),rep(1/3,12)),lower=rep(0,28),upper=rep(1,28),hessian=TRUE,print.level=0,iterlim=500)
-2*(-mnar1.fit$min-sat.lv) #Estatística de razão de verossimilhanças do ajuste do modelo, Tabela 3.2, MNAR1
p<-mnar1.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.3, MNAR1
	fr<-t(rep(1,15))%x%rowSums(TF)
	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]
round(matrix(c(	 fr[1,1]*(pi11.1*(1-a2.1.1-a3.1.1)),
	 fr[1,2]*(pi12.1*(1-a2.2.1-a3.1.1)),
	 fr[1,3]*(pi13.1*(1-a2.3.1-a3.1.1)),
	 fr[1,4]*(pi21.1*(1-a2.1.1-a3.2.1)),
	 fr[1,5]*(pi22.1*(1-a2.2.1-a3.2.1)),
	 fr[1,6]*(pi23.1*(1-a2.3.1-a3.2.1)),
	 fr[1,7]*(pi31.1*(1-a2.1.1-a3.3.1)),
	 fr[1,8]*(pi32.1*(1-a2.2.1-a3.3.1)),
	 fr[1,9]*(pi33.1*(1-a2.3.1-a3.3.1)),

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

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

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

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

	 fr[2,13]*c(pi11.2*a3.1.2 , pi21.2*a3.2.2 , pi31.2*a3.3.2),
	 fr[2,14]*c(pi12.2*a3.1.2 , pi22.2*a3.2.2 , pi32.2*a3.3.2),
	 fr[2,15]*c(pi13.2*a3.1.2 , pi23.2*a3.2.2 , pi33.2*a3.3.2)
),ncol=3,byrow=T),2) #Cenários t=3 precisam ser transpostos (linhas 7:9 e 16:18)

mnar1.analit<-mnar1.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])

loglin.mnar1<-funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar1.fit$est[1:16])),Vtheta=B%*%solve(attr(mnar1.analit,"hessian")[1,,])[1:16,1:16]%*%t(B),A1=A,X=rep(1,8))
summary(loglin.mnar1) #Logaritmos das razões de chances observadas e sob o modelo H:\omega_{ij(s)}=\omega, erros padrões e teste Wald de H, Tabela 3.5, MNAR1

funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar1.fit$est[1:16])),Vtheta=B%*%solve(attr(mnar1.analit,"hessian")[1,,])[1:16,1:16]%*%t(B),
  A1=A,U=UL) #Teste de Wald de H:\omega_{12(s)}=\omega_{21(s)}=\omega_{22(s)}=0,s=1,2, p.101, MNAR1





#menos log-verossimilhança do modelo MNAR2, saturado
mnar2.mlv<-function(p,fr=TF){
	# 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.2.1))+fr[1,3]*log(pi13.1*(1-a2.3.1-a3.3.1))+
	 fr[1,4]*log(pi21.1*(1-a2.2.1-a3.2.1))+fr[1,5]*log(pi22.1*(1-a2.1.1-a3.1.1))+fr[1,6]*log(pi23.1*(1-a2.2.1-a3.2.1))+
	 fr[1,7]*log(pi31.1*(1-a2.3.1-a3.3.1))+fr[1,8]*log(pi32.1*(1-a2.2.1-a3.2.1))+fr[1,9]*log(pi33.1*(1-a2.1.1-a3.1.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.2.1 + pi22.1*a2.1.1 + pi23.1*a2.2.1)+
	 fr[1,12]*log(pi31.1*a2.3.1 + pi32.1*a2.2.1 + pi33.1*a2.1.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.2.1 + pi22.1*a3.1.1 + pi32.1*a3.2.1)+
	 fr[1,15]*log(pi13.1*a3.3.1 + pi23.1*a3.2.1 + pi33.1*a3.1.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.2.2))+fr[2,3]*log(pi13.2*(1-a2.3.2-a3.3.2))+
	 fr[2,4]*log(pi21.2*(1-a2.2.2-a3.2.2))+fr[2,5]*log(pi22.2*(1-a2.1.2-a3.1.2))+fr[2,6]*log(pi23.2*(1-a2.2.2-a3.2.2))+
	 fr[2,7]*log(pi31.2*(1-a2.3.2-a3.3.2))+fr[2,8]*log(pi32.2*(1-a2.2.2-a3.2.2))+fr[2,9]*log(pi33.2*(1-a2.1.2-a3.1.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.2.2 + pi22.2*a2.1.2 + pi23.2*a2.2.2)+
	 fr[2,12]*log(pi31.2*a2.3.2 + pi32.2*a2.2.2 + pi33.2*a2.1.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.2.2 + pi22.2*a3.1.2 + pi32.2*a3.2.2)+
	 fr[2,15]*log(pi13.2*a3.3.2 + pi23.2*a3.2.2 + pi33.2*a3.1.2)

	)
	value
}

mnar2.der<-deriv3(~-(
	 o1.1*log(pi11.1*(1-a2.1.1-a3.1.1))+o1.2*log(pi12.1*(1-a2.2.1-a3.2.1))+o1.3*log(pi13.1*(1-a2.3.1-a3.3.1))+
	 o1.4*log(pi21.1*(1-a2.2.1-a3.2.1))+o1.5*log(pi22.1*(1-a2.1.1-a3.1.1))+o1.6*log(pi23.1*(1-a2.2.1-a3.2.1))+
	 o1.7*log(pi31.1*(1-a2.3.1-a3.3.1))+o1.8*log(pi32.1*(1-a2.2.1-a3.2.1))+o1.9*log((1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(1-a2.1.1-a3.1.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.2.1 + pi22.1*a2.1.1 + pi23.1*a2.2.1)+
	 o1.12*log(pi31.1*a2.3.1 + pi32.1*a2.2.1 + (1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*a2.1.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.2.1 + pi22.1*a3.1.1 + pi32.1*a3.2.1)+
	 o1.15*log(pi13.1*a3.3.1 + pi23.1*a3.2.1 + (1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*a3.1.1)+

	 o2.1*log(pi11.2*(1-a2.1.2-a3.1.2))+o2.2*log(pi12.2*(1-a2.2.2-a3.2.2))+o2.3*log(pi13.2*(1-a2.3.2-a3.3.2))+
	 o2.4*log(pi21.2*(1-a2.2.2-a3.2.2))+o2.5*log(pi22.2*(1-a2.1.2-a3.1.2))+o2.6*log(pi23.2*(1-a2.2.2-a3.2.2))+
	 o2.7*log(pi31.2*(1-a2.3.2-a3.3.2))+o2.8*log(pi32.2*(1-a2.2.2-a3.2.2))+o2.9*log((1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(1-a2.1.2-a3.1.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.2.2 + pi22.2*a2.1.2 + pi23.2*a2.2.2)+
	 o2.12*log(pi31.2*a2.3.2 + pi32.2*a2.2.2 + (1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*a2.1.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.2.2 + pi22.2*a3.1.2 + pi32.2*a3.2.2)+
	 o2.15*log(pi13.2*a3.3.2 + pi23.2*a3.2.2 + (1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*a3.1.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")
)

mnar2.fit<-.nlmP(objfunc=mnar2.mlv,params=c(rep(1/9,16),rep(1/3,12)),lower=rep(0,28),upper=rep(1,28),hessian=TRUE,print.level=0,iterlim=500)
-2*(-mnar2.fit$min-sat.lv) #Estatística de razão de verossimilhanças do ajuste do modelo, Tabela 3.2, MNAR2
p<-mnar2.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.4, MNAR2
	fr<-t(rep(1,15))%x%rowSums(TF)
	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]
round(matrix(c( 
	 fr[1,1]*(pi11.1*(1-a2.1.1-a3.1.1)),
   fr[1,2]*(pi12.1*(1-a2.2.1-a3.2.1)),
   fr[1,3]*(pi13.1*(1-a2.3.1-a3.3.1)),
	 fr[1,4]*(pi21.1*(1-a2.2.1-a3.2.1)),
   fr[1,5]*(pi22.1*(1-a2.1.1-a3.1.1)),
   fr[1,6]*(pi23.1*(1-a2.2.1-a3.2.1)),
	 fr[1,7]*(pi31.1*(1-a2.3.1-a3.3.1)),
   fr[1,8]*(pi32.1*(1-a2.2.1-a3.2.1)),
   fr[1,9]*(pi33.1*(1-a2.1.1-a3.1.1)),

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

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

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

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

	 fr[2,13]*c(pi11.2*a3.1.2 , pi21.2*a3.2.2 , pi31.2*a3.3.2),
	 fr[2,14]*c(pi12.2*a3.2.2 , pi22.2*a3.1.2 , pi32.2*a3.2.2),
	 fr[2,15]*c(pi13.2*a3.3.2 , pi23.2*a3.2.2 , pi33.2*a3.1.2)
),ncol=3,byrow=T),2) #Cenários t=3 precisam ser transpostos (linhas 7:9 e 16:18)

mnar2.analit<-mnar2.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])

loglin.mnar2<-funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar2.fit$est[1:16])),Vtheta=B%*%solve(attr(mnar2.analit,"hessian")[1,,])[1:16,1:16]%*%t(B),A1=A,X=rep(1,8))
summary(loglin.mnar2) #Logaritmos das razões de chances observadas e sob o modelo H:\omega_{ij(s)}=\omega, erros padrões e teste Wald de H, Tabela 3.5, MNAR2

funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar2.fit$est[1:16])),Vtheta=B%*%solve(attr(mnar2.analit,"hessian")[1,,])[1:16,1:16]%*%t(B),
  A1=A,U=UL) #Teste de Wald de H:\omega_{12(s)}=\omega_{21(s)}=\omega_{22(s)}=0,s=1,2, p.101, MNAR2





#menos log-verossimilhança do modelo MNAR3, 4 graus de liberdade
mnar3.mlv<-function(p,fr=TF){
	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

	a200.1<-p[17];a200.2<-p[18]
	a202<-p[19];a203<-p[20]

	a300.1<-p[21];a300.2<-p[22]
	a302<-p[23];a303<-p[24]

	value<- -(
	 fr[1,1]*log(pi11.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 fr[1,2]*log(pi12.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302))))+
	 fr[1,3]*log(pi13.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 fr[1,4]*log(pi21.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 fr[1,5]*log(pi22.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302))))+
	 fr[1,6]*log(pi23.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 fr[1,7]*log(pi31.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 fr[1,8]*log(pi32.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302))))+
	 fr[1,9]*log(pi33.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303))))+

	 fr[1,10]*log(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi12.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi13.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 fr[1,11]*log(pi21.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi22.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi23.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 fr[1,12]*log(pi31.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi32.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi33.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303))))+

	 fr[1,13]*log(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi21.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi31.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))))+
	 fr[1,14]*log(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi22.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi32.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))))+
	 fr[1,15]*log(pi13.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))) +
	  pi23.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))) +
	  pi33.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))))+

	 fr[2,1]*log(pi11.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 fr[2,2]*log(pi12.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302))))+
	 fr[2,3]*log(pi13.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 fr[2,4]*log(pi21.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 fr[2,5]*log(pi22.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302))))+
	 fr[2,6]*log(pi23.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 fr[2,7]*log(pi31.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 fr[2,8]*log(pi32.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302))))+
	 fr[2,9]*log(pi33.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303))))+

	 fr[2,10]*log(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi12.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi13.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 fr[2,11]*log(pi21.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi22.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi23.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 fr[2,12]*log(pi31.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi32.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi33.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303))))+

	 fr[2,13]*log(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi21.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi31.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))))+
	 fr[2,14]*log(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi22.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi32.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))))+
	 fr[2,15]*log(pi13.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))) +
	  pi23.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))) +
	  pi33.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))))

	)
	value
}

mnar3.der<-deriv3(~-(
	 o1.1*log(pi11.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 o1.2*log(pi12.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302))))+
	 o1.3*log(pi13.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 o1.4*log(pi21.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 o1.5*log(pi22.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302))))+
	 o1.6*log(pi23.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 o1.7*log(pi31.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 o1.8*log(pi32.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302))))+
	 o1.9*log((1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(1/(1+exp(a200.1+a203)+exp(a300.1+a303))))+

	 o1.10*log(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi12.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi13.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 o1.11*log(pi21.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi22.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi23.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 o1.12*log(pi31.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi32.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  (1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303))))+

	 o1.13*log(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi21.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi31.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))))+
	 o1.14*log(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi22.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi32.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))))+
	 o1.15*log(pi13.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))) +
	  pi23.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))) +
	  (1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))))+

	 o2.1*log(pi11.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.2*log(pi12.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302))))+
	 o2.3*log(pi13.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 o2.4*log(pi21.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.5*log(pi22.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302))))+
	 o2.6*log(pi23.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 o2.7*log(pi31.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.8*log(pi32.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302))))+
	 o2.9*log((1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(1/(1+exp(a200.2+a203)+exp(a300.2+a303))))+

	 o2.10*log(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi12.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi13.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 o2.11*log(pi21.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi22.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi23.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 o2.12*log(pi31.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi32.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  (1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303))))+

	 o2.13*log(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi21.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi31.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))))+
	 o2.14*log(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi22.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi32.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))))+
	 o2.15*log(pi13.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))) +
	  pi23.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))) +
	  (1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))))

	),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",
	 "a200.1","a200.2","a202","a203","a300.1","a300.2","a302","a303"),
	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",
	 "a200.1","a200.2","a202","a203","a300.1","a300.2","a302","a303",
	 "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")
)

mnar3.fit<-.nlmP(objfunc=mnar3.mlv,params=c(rep(1/9,16),rep(0,8)),lower=c(rep(0,16),rep(-Inf,8)),upper=c(rep(1,16),rep(Inf,8)),hessian=TRUE,print.level=0,iterlim=500)
-2*(-mnar3.fit$min-sat.lv) #Estatística de razão de verossimilhanças do ajuste do modelo, Tabela 3.2, MNAR3
1-pchisq(-2*(-mnar3.fit$min-sat.lv),4)
p<-mnar3.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.4, MNAR3
	fr<-t(rep(1,15))%x%rowSums(TF)
	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
	a200.1<-p[17];a200.2<-p[18]
	a202<-p[19];a203<-p[20]
	a300.1<-p[21];a300.2<-p[22]
	a302<-p[23];a303<-p[24]
round(rbind(	c( rowSums(TF)[1]*(pi11.1*(1/(1+exp(a200.1)+exp(a300.1)))),
	 rowSums(TF)[1]*(pi12.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302)))),
	 rowSums(TF)[1]*(pi13.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303)))) ),
	c( rowSums(TF)[1]*(pi21.1*(1/(1+exp(a200.1)+exp(a300.1)))),
	 rowSums(TF)[1]*(pi22.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302)))),
	 rowSums(TF)[1]*(pi23.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303)))) ),
	c( rowSums(TF)[1]*(pi31.1*(1/(1+exp(a200.1)+exp(a300.1)))),
	 rowSums(TF)[1]*(pi32.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302)))),
	 rowSums(TF)[1]*(pi33.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303)))) )
,
	 rowSums(TF)[1]*c(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))),
	  pi12.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))),
	  pi13.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303)))),
	 rowSums(TF)[1]*c(pi21.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))),
	  pi22.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))),
	  pi23.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303)))),
	 rowSums(TF)[1]*c(pi31.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))),
	  pi32.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))),
	  pi33.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303))))
,
	 rowSums(TF)[1]*c(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))),
	  pi21.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))),
	  pi31.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1)))),
	 rowSums(TF)[1]*c(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))),
	  pi22.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))),
	  pi32.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302)))),
	 rowSums(TF)[1]*c(pi13.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))),
	  pi23.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))),
	  pi33.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))))
,
	c( rowSums(TF)[2]*(pi11.2*(1/(1+exp(a200.2)+exp(a300.2)))),
	 rowSums(TF)[2]*(pi12.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302)))),
	 rowSums(TF)[2]*(pi13.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303)))) ),
	c( rowSums(TF)[2]*(pi21.2*(1/(1+exp(a200.2)+exp(a300.2)))),
	 rowSums(TF)[2]*(pi22.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302)))),
	 rowSums(TF)[2]*(pi23.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303)))) ),
	c( rowSums(TF)[2]*(pi31.2*(1/(1+exp(a200.2)+exp(a300.2)))),
	 rowSums(TF)[2]*(pi32.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302)))),
	 rowSums(TF)[2]*(pi33.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303)))) )
,
	 rowSums(TF)[2]*c(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))),
	  pi12.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))),
	  pi13.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303)))),
	 rowSums(TF)[2]*c(pi21.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))),
	  pi22.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))),
	  pi23.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303)))),
	 rowSums(TF)[2]*c(pi31.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))),
	  pi32.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))),
	  pi33.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303))))
,
	 rowSums(TF)[2]*c(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))),
	  pi21.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))),
	  pi31.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2)))),
	 rowSums(TF)[2]*c(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))),
	  pi22.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))),
	  pi32.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302)))),
	 rowSums(TF)[2]*c(pi13.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))),
	  pi23.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))),
	  pi33.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))))
),2) #Cenários t=3 precisam ser transpostos (linhas 7:9 e 16:18)

mnar3.analit<-mnar3.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],
	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])

loglin.mnar3<-funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar3.fit$est[1:16])),Vtheta=B%*%solve(attr(mnar3.analit,"hessian")[1,,],tol=1e-50)[1:16,1:16]%*%t(B),A1=A,X=rep(1,8))
summary(loglin.mnar3) #Logaritmos das razões de chances observadas e sob o modelo H:\omega_{ij(s)}=\omega, erros padrões e teste Wald de H, Tabela 3.5, MNAR3

funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar3.fit$est[1:16])),Vtheta=B%*%solve(attr(mnar3.analit,"hessian")[1,,],tol=1e-50)[1:16,1:16]%*%t(B),
  A1=A,U=UL) #Teste de Wald de H:\omega_{12(s)}=\omega_{21(s)}=\omega_{22(s)}=0,s=1,2, p.101, MNAR3





#menos log-verossimilhança da estimação do modelo MNAR4, saturado
mnar4.mlv<-function(p,fr=TF){
	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

	a200.1<-p[17];a200.2<-p[18]
	a220<-p[19];a230<-p[20];a202<-p[21];a203<-p[22]

	a300.1<-p[23];a300.2<-p[24]
	a320<-p[25];a330<-p[26];a302<-p[27];a303<-p[28]

	value<- -(
	 fr[1,1]*log(pi11.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 fr[1,2]*log(pi12.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302))))+
	 fr[1,3]*log(pi13.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 fr[1,4]*log(pi21.1*(1/(1+exp(a200.1+a220)+exp(a300.1+a320))))+
	 fr[1,5]*log(pi22.1*(1/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))))+
	 fr[1,6]*log(pi23.1*(1/(1+exp(a200.1+a220+a203)+exp(a300.1+a320+a303))))+
	 fr[1,7]*log(pi31.1*(1/(1+exp(a200.1+a230)+exp(a300.1+a330))))+
	 fr[1,8]*log(pi32.1*(1/(1+exp(a200.1+a230+a202)+exp(a300.1+a330+a302))))+
	 fr[1,9]*log(pi33.1*(1/(1+exp(a200.1+a230+a203)+exp(a300.1+a330+a303))))+

	 fr[1,10]*log(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi12.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi13.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 fr[1,11]*log(pi21.1*(exp(a200.1+a220)/(1+exp(a200.1+a220)+exp(a300.1+a320))) +
	  pi22.1*(exp(a200.1+a220+a202)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))) +
	  pi23.1*(exp(a200.1+a220+a203)/(1+exp(a200.1+a220+a203)+exp(a300.1+a320+a303))))+
	 fr[1,12]*log(pi31.1*(exp(a200.1+a230)/(1+exp(a200.1+a230)+exp(a300.1+a330))) +
	  pi32.1*(exp(a200.1+a230+a202)/(1+exp(a200.1+a230+a202)+exp(a300.1+a330+a302))) +
	  pi33.1*(exp(a200.1+a230+a203)/(1+exp(a200.1+a230+a203)+exp(a300.1+a330+a303))))+

	 fr[1,13]*log(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi21.1*(exp(a300.1+a320)/(1+exp(a200.1+a220)+exp(a300.1+a320))) +
	  pi31.1*(exp(a300.1+a330)/(1+exp(a200.1+a230)+exp(a300.1+a330))))+
	 fr[1,14]*log(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi22.1*(exp(a300.1+a320+a302)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))) +
	  pi32.1*(exp(a300.1+a330+a302)/(1+exp(a200.1+a230+a202)+exp(a300.1+a330+a302))))+
	 fr[1,15]*log(pi13.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))) +
	  pi23.1*(exp(a300.1+a320+a303)/(1+exp(a200.1+a220+a203)+exp(a300.1+a320+a303))) +
	  pi33.1*(exp(a300.1+a330+a303)/(1+exp(a200.1+a230+a203)+exp(a300.1+a330+a303))))+

	 fr[2,1]*log(pi11.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 fr[2,2]*log(pi12.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302))))+
	 fr[2,3]*log(pi13.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 fr[2,4]*log(pi21.2*(1/(1+exp(a200.2+a220)+exp(a300.2+a320))))+
	 fr[2,5]*log(pi22.2*(1/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))))+
	 fr[2,6]*log(pi23.2*(1/(1+exp(a200.2+a220+a203)+exp(a300.2+a320+a303))))+
	 fr[2,7]*log(pi31.2*(1/(1+exp(a200.2+a230)+exp(a300.2+a330))))+
	 fr[2,8]*log(pi32.2*(1/(1+exp(a200.2+a230+a202)+exp(a300.2+a330+a302))))+
	 fr[2,9]*log(pi33.2*(1/(1+exp(a200.2+a230+a203)+exp(a300.2+a330+a303))))+

	 fr[2,10]*log(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi12.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi13.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 fr[2,11]*log(pi21.2*(exp(a200.2+a220)/(1+exp(a200.2+a220)+exp(a300.2+a320))) +
	  pi22.2*(exp(a200.2+a220+a202)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))) +
	  pi23.2*(exp(a200.2+a220+a203)/(1+exp(a200.2+a220+a203)+exp(a300.2+a320+a303))))+
	 fr[2,12]*log(pi31.2*(exp(a200.2+a230)/(1+exp(a200.2+a230)+exp(a300.2+a330))) +
	  pi32.2*(exp(a200.2+a230+a202)/(1+exp(a200.2+a230+a202)+exp(a300.2+a330+a302))) +
	  pi33.2*(exp(a200.2+a230+a203)/(1+exp(a200.2+a230+a203)+exp(a300.2+a330+a303))))+

	 fr[2,13]*log(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi21.2*(exp(a300.2+a320)/(1+exp(a200.2+a220)+exp(a300.2+a320))) +
	  pi31.2*(exp(a300.2+a330)/(1+exp(a200.2+a230)+exp(a300.2+a330))))+
	 fr[2,14]*log(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi22.2*(exp(a300.2+a320+a302)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))) +
	  pi32.2*(exp(a300.2+a330+a302)/(1+exp(a200.2+a230+a202)+exp(a300.2+a330+a302))))+
	 fr[2,15]*log(pi13.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))) +
	  pi23.2*(exp(a300.2+a320+a303)/(1+exp(a200.2+a220+a203)+exp(a300.2+a320+a303))) +
	  pi33.2*(exp(a300.2+a330+a303)/(1+exp(a200.2+a230+a203)+exp(a300.2+a330+a303))))

	)
	value
}

mnar4.der<-deriv3(~-(
	 o1.1*log(pi11.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 o1.2*log(pi12.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302))))+
	 o1.3*log(pi13.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 o1.4*log(pi21.1*(1/(1+exp(a200.1+a220)+exp(a300.1+a320))))+
	 o1.5*log(pi22.1*(1/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))))+
	 o1.6*log(pi23.1*(1/(1+exp(a200.1+a220+a203)+exp(a300.1+a320+a303))))+
	 o1.7*log(pi31.1*(1/(1+exp(a200.1+a230)+exp(a300.1+a330))))+
	 o1.8*log(pi32.1*(1/(1+exp(a200.1+a230+a202)+exp(a300.1+a330+a302))))+
	 o1.9*log((1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(1/(1+exp(a200.1+a230+a203)+exp(a300.1+a330+a303))))+

	 o1.10*log(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi12.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi13.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303))))+
	 o1.11*log(pi21.1*(exp(a200.1+a220)/(1+exp(a200.1+a220)+exp(a300.1+a320))) +
	  pi22.1*(exp(a200.1+a220+a202)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))) +
	  pi23.1*(exp(a200.1+a220+a203)/(1+exp(a200.1+a220+a203)+exp(a300.1+a320+a303))))+
	 o1.12*log(pi31.1*(exp(a200.1+a230)/(1+exp(a200.1+a230)+exp(a300.1+a330))) +
	  pi32.1*(exp(a200.1+a230+a202)/(1+exp(a200.1+a230+a202)+exp(a300.1+a330+a302))) +
	  (1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(exp(a200.1+a230+a203)/(1+exp(a200.1+a230+a203)+exp(a300.1+a330+a303))))+

	 o1.13*log(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi21.1*(exp(a300.1+a320)/(1+exp(a200.1+a220)+exp(a300.1+a320))) +
	  pi31.1*(exp(a300.1+a330)/(1+exp(a200.1+a230)+exp(a300.1+a330))))+
	 o1.14*log(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  pi22.1*(exp(a300.1+a320+a302)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))) +
	  pi32.1*(exp(a300.1+a330+a302)/(1+exp(a200.1+a230+a202)+exp(a300.1+a330+a302))))+
	 o1.15*log(pi13.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))) +
	  pi23.1*(exp(a300.1+a320+a303)/(1+exp(a200.1+a220+a203)+exp(a300.1+a320+a303))) +
	  (1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(exp(a300.1+a330+a303)/(1+exp(a200.1+a230+a203)+exp(a300.1+a330+a303))))+

	 o2.1*log(pi11.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.2*log(pi12.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302))))+
	 o2.3*log(pi13.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 o2.4*log(pi21.2*(1/(1+exp(a200.2+a220)+exp(a300.2+a320))))+
	 o2.5*log(pi22.2*(1/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))))+
	 o2.6*log(pi23.2*(1/(1+exp(a200.2+a220+a203)+exp(a300.2+a320+a303))))+
	 o2.7*log(pi31.2*(1/(1+exp(a200.2+a230)+exp(a300.2+a330))))+
	 o2.8*log(pi32.2*(1/(1+exp(a200.2+a230+a202)+exp(a300.2+a330+a302))))+
	 o2.9*log((1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(1/(1+exp(a200.2+a230+a203)+exp(a300.2+a330+a303))))+

	 o2.10*log(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi12.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi13.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303))))+
	 o2.11*log(pi21.2*(exp(a200.2+a220)/(1+exp(a200.2+a220)+exp(a300.2+a320))) +
	  pi22.2*(exp(a200.2+a220+a202)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))) +
	  pi23.2*(exp(a200.2+a220+a203)/(1+exp(a200.2+a220+a203)+exp(a300.2+a320+a303))))+
	 o2.12*log(pi31.2*(exp(a200.2+a230)/(1+exp(a200.2+a230)+exp(a300.2+a330))) +
	  pi32.2*(exp(a200.2+a230+a202)/(1+exp(a200.2+a230+a202)+exp(a300.2+a330+a302))) +
	  (1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(exp(a200.2+a230+a203)/(1+exp(a200.2+a230+a203)+exp(a300.2+a330+a303))))+

	 o2.13*log(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi21.2*(exp(a300.2+a320)/(1+exp(a200.2+a220)+exp(a300.2+a320))) +
	  pi31.2*(exp(a300.2+a330)/(1+exp(a200.2+a230)+exp(a300.2+a330))))+
	 o2.14*log(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  pi22.2*(exp(a300.2+a320+a302)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))) +
	  pi32.2*(exp(a300.2+a330+a302)/(1+exp(a200.2+a230+a202)+exp(a300.2+a330+a302))))+
	 o2.15*log(pi13.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))) +
	  pi23.2*(exp(a300.2+a320+a303)/(1+exp(a200.2+a220+a203)+exp(a300.2+a320+a303))) +
	  (1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(exp(a300.2+a330+a303)/(1+exp(a200.2+a230+a203)+exp(a300.2+a330+a303))))

	),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",
	 "a200.1","a200.2","a220","a230","a202","a203","a300.1","a300.2","a320","a330","a302","a303"),
	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",
	 "a200.1","a200.2","a220","a230","a202","a203","a300.1","a300.2","a320","a330","a302","a303",
	 "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")
)

mnar4.fit<-.nlmP(objfunc=mnar4.mlv,params=c(rep(1/9,16),rep(0,12)),lower=c(rep(0,16),rep(-Inf,12)),upper=c(rep(1,16),rep(Inf,12)),hessian=TRUE,print.level=0,iterlim=500)
-2*(-mnar4.fit$min-sat.lv) #Estatística de razão de verossimilhanças do ajuste do modelo, Tabela 3.2, MNAR4
p<-mnar4.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.4, MNAR4
	fr<-t(rep(1,15))%x%rowSums(TF)
	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
	a200.1<-p[17];a200.2<-p[18]
	a220<-p[19];a230<-p[20];a202<-p[21];a203<-p[22]
	a300.1<-p[23];a300.2<-p[24]
	a320<-p[25];a330<-p[26];a302<-p[27];a303<-p[28]
round(matrix(c( fr[1,1]*(pi11.1*(1/(1+exp(a200.1)+exp(a300.1)))),
	 fr[1,2]*(pi12.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302)))),
	 fr[1,3]*(pi13.1*(1/(1+exp(a200.1+a203)+exp(a300.1+a303)))),
	 fr[1,4]*(pi21.1*(1/(1+exp(a200.1+a220)+exp(a300.1+a320)))),
	 fr[1,5]*(pi22.1*(1/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302)))),
	 fr[1,6]*(pi23.1*(1/(1+exp(a200.1+a220+a203)+exp(a300.1+a320+a303)))),
	 fr[1,7]*(pi31.1*(1/(1+exp(a200.1+a230)+exp(a300.1+a330)))),
	 fr[1,8]*(pi32.1*(1/(1+exp(a200.1+a230+a202)+exp(a300.1+a330+a302)))),
	 fr[1,9]*(pi33.1*(1/(1+exp(a200.1+a230+a203)+exp(a300.1+a330+a303)))),

	 fr[1,10]*c(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) ,
	  pi12.1*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))) ,
	  pi13.1*(exp(a200.1+a203)/(1+exp(a200.1+a203)+exp(a300.1+a303)))),
	 fr[1,11]*c(pi21.1*(exp(a200.1+a220)/(1+exp(a200.1+a220)+exp(a300.1+a320))) ,
	  pi22.1*(exp(a200.1+a220+a202)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))) ,
	  pi23.1*(exp(a200.1+a220+a203)/(1+exp(a200.1+a220+a203)+exp(a300.1+a320+a303)))),
	 fr[1,12]*c(pi31.1*(exp(a200.1+a230)/(1+exp(a200.1+a230)+exp(a300.1+a330))) ,
	  pi32.1*(exp(a200.1+a230+a202)/(1+exp(a200.1+a230+a202)+exp(a300.1+a330+a302))) ,
	  pi33.1*(exp(a200.1+a230+a203)/(1+exp(a200.1+a230+a203)+exp(a300.1+a330+a303)))),

	 fr[1,13]*c(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) ,
	  pi21.1*(exp(a300.1+a320)/(1+exp(a200.1+a220)+exp(a300.1+a320))) ,
	  pi31.1*(exp(a300.1+a330)/(1+exp(a200.1+a230)+exp(a300.1+a330)))),
	 fr[1,14]*c(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) ,
	  pi22.1*(exp(a300.1+a320+a302)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))) ,
	  pi32.1*(exp(a300.1+a330+a302)/(1+exp(a200.1+a230+a202)+exp(a300.1+a330+a302)))),
	 fr[1,15]*c(pi13.1*(exp(a300.1+a303)/(1+exp(a200.1+a203)+exp(a300.1+a303))) ,
	  pi23.1*(exp(a300.1+a320+a303)/(1+exp(a200.1+a220+a203)+exp(a300.1+a320+a303))) ,
	  pi33.1*(exp(a300.1+a330+a303)/(1+exp(a200.1+a230+a203)+exp(a300.1+a330+a303)))),

	 fr[2,1]*(pi11.2*(1/(1+exp(a200.2)+exp(a300.2)))),
	 fr[2,2]*(pi12.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302)))),
	 fr[2,3]*(pi13.2*(1/(1+exp(a200.2+a203)+exp(a300.2+a303)))),
	 fr[2,4]*(pi21.2*(1/(1+exp(a200.2+a220)+exp(a300.2+a320)))),
	 fr[2,5]*(pi22.2*(1/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302)))),
	 fr[2,6]*(pi23.2*(1/(1+exp(a200.2+a220+a203)+exp(a300.2+a320+a303)))),
	 fr[2,7]*(pi31.2*(1/(1+exp(a200.2+a230)+exp(a300.2+a330)))),
	 fr[2,8]*(pi32.2*(1/(1+exp(a200.2+a230+a202)+exp(a300.2+a330+a302)))),
	 fr[2,9]*(pi33.2*(1/(1+exp(a200.2+a230+a203)+exp(a300.2+a330+a303)))),

	 fr[2,10]*c(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) ,
	  pi12.2*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))) ,
	  pi13.2*(exp(a200.2+a203)/(1+exp(a200.2+a203)+exp(a300.2+a303)))),
	 fr[2,11]*c(pi21.2*(exp(a200.2+a220)/(1+exp(a200.2+a220)+exp(a300.2+a320))) ,
	  pi22.2*(exp(a200.2+a220+a202)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))) ,
	  pi23.2*(exp(a200.2+a220+a203)/(1+exp(a200.2+a220+a203)+exp(a300.2+a320+a303)))),
	 fr[2,12]*c(pi31.2*(exp(a200.2+a230)/(1+exp(a200.2+a230)+exp(a300.2+a330))) ,
	  pi32.2*(exp(a200.2+a230+a202)/(1+exp(a200.2+a230+a202)+exp(a300.2+a330+a302))) ,
	  pi33.2*(exp(a200.2+a230+a203)/(1+exp(a200.2+a230+a203)+exp(a300.2+a330+a303)))),

	 fr[2,13]*c(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) ,
	  pi21.2*(exp(a300.2+a320)/(1+exp(a200.2+a220)+exp(a300.2+a320))) ,
	  pi31.2*(exp(a300.2+a330)/(1+exp(a200.2+a230)+exp(a300.2+a330)))),
	 fr[2,14]*c(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) ,
	  pi22.2*(exp(a300.2+a320+a302)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))) ,
	  pi32.2*(exp(a300.2+a330+a302)/(1+exp(a200.2+a230+a202)+exp(a300.2+a330+a302)))),
	 fr[2,15]*c(pi13.2*(exp(a300.2+a303)/(1+exp(a200.2+a203)+exp(a300.2+a303))) ,
	  pi23.2*(exp(a300.2+a320+a303)/(1+exp(a200.2+a220+a203)+exp(a300.2+a320+a303))) ,
	  pi33.2*(exp(a300.2+a330+a303)/(1+exp(a200.2+a230+a203)+exp(a300.2+a330+a303))))
),ncol=3,byrow=T),2) #Cenários t=3 precisam ser transpostos (linhas 7:9 e 16:18)

mnar4.analit<-mnar4.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])

loglin.mnar4<-funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar4.fit$est[1:16])),Vtheta=B%*%solve(attr(mnar4.analit,"hessian")[1,,],tol=1e-50)[1:16,1:16]%*%t(B),A1=A,X=rep(1,8))
summary(loglin.mnar4) #Logaritmos das razões de chances observadas e sob o modelo H:\omega_{ij(s)}=\omega, erros padrões e teste Wald de H, Tabela 3.5, MNAR4

funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar4.fit$est[1:16])),Vtheta=B%*%solve(attr(mnar4.analit,"hessian")[1,,],tol=1e-50)[1:16,1:16]%*%t(B),
  A1=A,U=UL) #Teste de Wald de H:\omega_{12(s)}=\omega_{21(s)}=\omega_{22(s)}=0,s=1,2, p.101, MNAR4





#menos log-verossimilhança da estimação do modelo MNAR5, saturado
mnar5.mlv<-function(p,fr=TF){
	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

	a200.1<-p[17];a300.1<-p[18]
	a020<-p[19];a030<-p[20];a002<-p[21];a003<-p[22]
	a022<-p[23];a023<-p[24];a032<-p[25];a033<-p[26]
	a200.2<-p[27];a300.2<-p[28]

	value<- -(
	 fr[1,1]*log(pi11.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 fr[1,2]*log(pi12.1*(1/(1+exp(a200.1+a002)+exp(a300.1+a002))))+
	 fr[1,3]*log(pi13.1*(1/(1+exp(a200.1+a003)+exp(a300.1+a003))))+
	 fr[1,4]*log(pi21.1*(1/(1+exp(a200.1+a020)+exp(a300.1+a020))))+
	 fr[1,5]*log(pi22.1*(1/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))))+
	 fr[1,6]*log(pi23.1*(1/(1+exp(a200.1+a020+a003+a023)+exp(a300.1+a020+a003+a023))))+
	 fr[1,7]*log(pi31.1*(1/(1+exp(a200.1+a030)+exp(a300.1+a030))))+
	 fr[1,8]*log(pi32.1*(1/(1+exp(a200.1+a030+a002+a032)+exp(a300.1+a030+a002+a032))))+
	 fr[1,9]*log(pi33.1*(1/(1+exp(a200.1+a030+a003+a033)+exp(a300.1+a030+a003+a033))))+

	 fr[1,10]*log(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi12.1*(exp(a200.1+a002)/(1+exp(a200.1+a002)+exp(a300.1+a002))) +
	  pi13.1*(exp(a200.1+a003)/(1+exp(a200.1+a003)+exp(a300.1+a003))))+
	 fr[1,11]*log(pi21.1*(exp(a200.1+a020)/(1+exp(a200.1+a020)+exp(a300.1+a020))) +
	  pi22.1*(exp(a200.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))) +
	  pi23.1*(exp(a200.1+a020+a003+a023)/(1+exp(a200.1+a020+a003+a023)+exp(a300.1+a020+a003+a023))))+
	 fr[1,12]*log(pi31.1*(exp(a200.1+a030)/(1+exp(a200.1+a030)+exp(a300.1+a030))) +
	  pi32.1*(exp(a200.1+a030+a002+a032)/(1+exp(a200.1+a030+a002+a032)+exp(a300.1+a030+a002+a032))) +
	  pi33.1*(exp(a200.1+a030+a003+a033)/(1+exp(a200.1+a030+a003+a033)+exp(a300.1+a030+a003+a033))))+

	 fr[1,13]*log(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi21.1*(exp(a300.1+a020)/(1+exp(a200.1+a020)+exp(a300.1+a020))) +
	  pi31.1*(exp(a300.1+a030)/(1+exp(a200.1+a030)+exp(a300.1+a030))))+
	 fr[1,14]*log(pi12.1*(exp(a300.1+a002)/(1+exp(a200.1+a002)+exp(a300.1+a002))) +
	  pi22.1*(exp(a300.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))) +
	  pi32.1*(exp(a300.1+a030+a002+a032)/(1+exp(a200.1+a030+a002+a032)+exp(a300.1+a030+a002+a032))))+
	 fr[1,15]*log(pi13.1*(exp(a300.1+a003)/(1+exp(a200.1+a003)+exp(a300.1+a003))) +
	  pi23.1*(exp(a300.1+a020+a003+a023)/(1+exp(a200.1+a020+a003+a023)+exp(a300.1+a020+a003+a023))) +
	  pi33.1*(exp(a300.1+a030+a003+a033)/(1+exp(a200.1+a030+a003+a033)+exp(a300.1+a030+a003+a033))))+

	 fr[2,1]*log(pi11.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 fr[2,2]*log(pi12.2*(1/(1+exp(a200.2+a002)+exp(a300.2+a002))))+
	 fr[2,3]*log(pi13.2*(1/(1+exp(a200.2+a003)+exp(a300.2+a003))))+
	 fr[2,4]*log(pi21.2*(1/(1+exp(a200.2+a020)+exp(a300.2+a020))))+
	 fr[2,5]*log(pi22.2*(1/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))))+
	 fr[2,6]*log(pi23.2*(1/(1+exp(a200.2+a020+a003+a023)+exp(a300.2+a020+a003+a023))))+
	 fr[2,7]*log(pi31.2*(1/(1+exp(a200.2+a030)+exp(a300.2+a030))))+
	 fr[2,8]*log(pi32.2*(1/(1+exp(a200.2+a030+a002+a032)+exp(a300.2+a030+a002+a032))))+
	 fr[2,9]*log(pi33.2*(1/(1+exp(a200.2+a030+a003+a033)+exp(a300.2+a030+a003+a033))))+

	 fr[2,10]*log(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi12.2*(exp(a200.2+a002)/(1+exp(a200.2+a002)+exp(a300.2+a002))) +
	  pi13.2*(exp(a200.2+a003)/(1+exp(a200.2+a003)+exp(a300.2+a003))))+
	 fr[2,11]*log(pi21.2*(exp(a200.2+a020)/(1+exp(a200.2+a020)+exp(a300.2+a020))) +
	  pi22.2*(exp(a200.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))) +
	  pi23.2*(exp(a200.2+a020+a003+a023)/(1+exp(a200.2+a020+a003+a023)+exp(a300.2+a020+a003+a023))))+
	 fr[2,12]*log(pi31.2*(exp(a200.2+a030)/(1+exp(a200.2+a030)+exp(a300.2+a030))) +
	  pi32.2*(exp(a200.2+a030+a002+a032)/(1+exp(a200.2+a030+a002+a032)+exp(a300.2+a030+a002+a032))) +
	  pi33.2*(exp(a200.2+a030+a003+a033)/(1+exp(a200.2+a030+a003+a033)+exp(a300.2+a030+a003+a033))))+

	 fr[2,13]*log(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi21.2*(exp(a300.2+a020)/(1+exp(a200.2+a020)+exp(a300.2+a020))) +
	  pi31.2*(exp(a300.2+a030)/(1+exp(a200.2+a030)+exp(a300.2+a030))))+
	 fr[2,14]*log(pi12.2*(exp(a300.2+a002)/(1+exp(a200.2+a002)+exp(a300.2+a002))) +
	  pi22.2*(exp(a300.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))) +
	  pi32.2*(exp(a300.2+a030+a002+a032)/(1+exp(a200.2+a030+a002+a032)+exp(a300.2+a030+a002+a032))))+
	 fr[2,15]*log(pi13.2*(exp(a300.2+a003)/(1+exp(a200.2+a003)+exp(a300.2+a003))) +
	  pi23.2*(exp(a300.2+a020+a003+a023)/(1+exp(a200.2+a020+a003+a023)+exp(a300.2+a020+a003+a023))) +
	  pi33.2*(exp(a300.2+a030+a003+a033)/(1+exp(a200.2+a030+a003+a033)+exp(a300.2+a030+a003+a033))))

	)
	value
}

mnar5.der<-deriv3(~-(
	 o1.1*log(pi11.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 o1.2*log(pi12.1*(1/(1+exp(a200.1+a002)+exp(a300.1+a002))))+
	 o1.3*log(pi13.1*(1/(1+exp(a200.1+a003)+exp(a300.1+a003))))+
	 o1.4*log(pi21.1*(1/(1+exp(a200.1+a020)+exp(a300.1+a020))))+
	 o1.5*log(pi22.1*(1/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))))+
	 o1.6*log(pi23.1*(1/(1+exp(a200.1+a020+a003+a023)+exp(a300.1+a020+a003+a023))))+
	 o1.7*log(pi31.1*(1/(1+exp(a200.1+a030)+exp(a300.1+a030))))+
	 o1.8*log(pi32.1*(1/(1+exp(a200.1+a030+a002+a032)+exp(a300.1+a030+a002+a032))))+
	 o1.9*log((1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(1/(1+exp(a200.1+a030+a003+a033)+exp(a300.1+a030+a003+a033))))+

	 o1.10*log(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi12.1*(exp(a200.1+a002)/(1+exp(a200.1+a002)+exp(a300.1+a002))) +
	  pi13.1*(exp(a200.1+a003)/(1+exp(a200.1+a003)+exp(a300.1+a003))))+
	 o1.11*log(pi21.1*(exp(a200.1+a020)/(1+exp(a200.1+a020)+exp(a300.1+a020))) +
	  pi22.1*(exp(a200.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))) +
	  pi23.1*(exp(a200.1+a020+a003+a023)/(1+exp(a200.1+a020+a003+a023)+exp(a300.1+a020+a003+a023))))+
	 o1.12*log(pi31.1*(exp(a200.1+a030)/(1+exp(a200.1+a030)+exp(a300.1+a030))) +
	  pi32.1*(exp(a200.1+a030+a002+a032)/(1+exp(a200.1+a030+a002+a032)+exp(a300.1+a030+a002+a032))) +
	  (1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(exp(a200.1+a030+a003+a033)/(1+exp(a200.1+a030+a003+a033)+exp(a300.1+a030+a003+a033))))+

	 o1.13*log(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi21.1*(exp(a300.1+a020)/(1+exp(a200.1+a020)+exp(a300.1+a020))) +
	  pi31.1*(exp(a300.1+a030)/(1+exp(a200.1+a030)+exp(a300.1+a030))))+
	 o1.14*log(pi12.1*(exp(a300.1+a002)/(1+exp(a200.1+a002)+exp(a300.1+a002))) +
	  pi22.1*(exp(a300.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))) +
	  pi32.1*(exp(a300.1+a030+a002+a032)/(1+exp(a200.1+a030+a002+a032)+exp(a300.1+a030+a002+a032))))+
	 o1.15*log(pi13.1*(exp(a300.1+a003)/(1+exp(a200.1+a003)+exp(a300.1+a003))) +
	  pi23.1*(exp(a300.1+a020+a003+a023)/(1+exp(a200.1+a020+a003+a023)+exp(a300.1+a020+a003+a023))) +
	  (1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*(exp(a300.1+a030+a003+a033)/(1+exp(a200.1+a030+a003+a033)+exp(a300.1+a030+a003+a033))))+

	 o2.1*log(pi11.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.2*log(pi12.2*(1/(1+exp(a200.2+a002)+exp(a300.2+a002))))+
	 o2.3*log(pi13.2*(1/(1+exp(a200.2+a003)+exp(a300.2+a003))))+
	 o2.4*log(pi21.2*(1/(1+exp(a200.2+a020)+exp(a300.2+a020))))+
	 o2.5*log(pi22.2*(1/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))))+
	 o2.6*log(pi23.2*(1/(1+exp(a200.2+a020+a003+a023)+exp(a300.2+a020+a003+a023))))+
	 o2.7*log(pi31.2*(1/(1+exp(a200.2+a030)+exp(a300.2+a030))))+
	 o2.8*log(pi32.2*(1/(1+exp(a200.2+a030+a002+a032)+exp(a300.2+a030+a002+a032))))+
	 o2.9*log((1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(1/(1+exp(a200.2+a030+a003+a033)+exp(a300.2+a030+a003+a033))))+

	 o2.10*log(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi12.2*(exp(a200.2+a002)/(1+exp(a200.2+a002)+exp(a300.2+a002))) +
	  pi13.2*(exp(a200.2+a003)/(1+exp(a200.2+a003)+exp(a300.2+a003))))+
	 o2.11*log(pi21.2*(exp(a200.2+a020)/(1+exp(a200.2+a020)+exp(a300.2+a020))) +
	  pi22.2*(exp(a200.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))) +
	  pi23.2*(exp(a200.2+a020+a003+a023)/(1+exp(a200.2+a020+a003+a023)+exp(a300.2+a020+a003+a023))))+
	 o2.12*log(pi31.2*(exp(a200.2+a030)/(1+exp(a200.2+a030)+exp(a300.2+a030))) +
	  pi32.2*(exp(a200.2+a030+a002+a032)/(1+exp(a200.2+a030+a002+a032)+exp(a300.2+a030+a002+a032))) +
	  (1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(exp(a200.2+a030+a003+a033)/(1+exp(a200.2+a030+a003+a033)+exp(a300.2+a030+a003+a033))))+

	 o2.13*log(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi21.2*(exp(a300.2+a020)/(1+exp(a200.2+a020)+exp(a300.2+a020))) +
	  pi31.2*(exp(a300.2+a030)/(1+exp(a200.2+a030)+exp(a300.2+a030))))+
	 o2.14*log(pi12.2*(exp(a300.2+a002)/(1+exp(a200.2+a002)+exp(a300.2+a002))) +
	  pi22.2*(exp(a300.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))) +
	  pi32.2*(exp(a300.2+a030+a002+a032)/(1+exp(a200.2+a030+a002+a032)+exp(a300.2+a030+a002+a032))))+
	 o2.15*log(pi13.2*(exp(a300.2+a003)/(1+exp(a200.2+a003)+exp(a300.2+a003))) +
	  pi23.2*(exp(a300.2+a020+a003+a023)/(1+exp(a200.2+a020+a003+a023)+exp(a300.2+a020+a003+a023))) +
	  (1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*(exp(a300.2+a030+a003+a033)/(1+exp(a200.2+a030+a003+a033)+exp(a300.2+a030+a003+a033))))

	),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",
	 "a200.1","a300.1","a020","a030","a002","a003","a022","a023","a032","a033","a200.2","a300.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",
	 "a200.1","a300.1","a020","a030","a002","a003","a022","a023","a032","a033","a200.2","a300.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")
)

mnar5.fit<-.nlmP(objfunc=mnar5.mlv,params=c(rep(1/9,16),rep(0,12)),lower=c(rep(0,16),rep(-Inf,12)),upper=c(rep(1,16),rep(Inf,12)),hessian=TRUE,print.level=0,iterlim=500)
-2*(-mnar5.fit$min-sat.lv) #Estatística de razão de verossimilhanças do ajuste do modelo, Tabela 3.2, MNAR4
p<-mnar5.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.4, MNAR4
	fr<-t(rep(1,15))%x%rowSums(TF)
	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
	a200.1<-p[17];a300.1<-p[18]
	a020<-p[19];a030<-p[20];a002<-p[21];a003<-p[22]
	a022<-p[23];a023<-p[24];a032<-p[25];a033<-p[26]
	a200.2<-p[27];a300.2<-p[28]
round(matrix(c( fr[1,1]*(pi11.1*(1/(1+exp(a200.1)+exp(a300.1)))),
	 fr[1,2]*(pi12.1*(1/(1+exp(a200.1+a002)+exp(a300.1+a002)))),
	 fr[1,3]*(pi13.1*(1/(1+exp(a200.1+a003)+exp(a300.1+a003)))),
	 fr[1,4]*(pi21.1*(1/(1+exp(a200.1+a020)+exp(a300.1+a020)))),
	 fr[1,5]*(pi22.1*(1/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022)))),
	 fr[1,6]*(pi23.1*(1/(1+exp(a200.1+a020+a003+a023)+exp(a300.1+a020+a003+a023)))),
	 fr[1,7]*(pi31.1*(1/(1+exp(a200.1+a030)+exp(a300.1+a030)))),
	 fr[1,8]*(pi32.1*(1/(1+exp(a200.1+a030+a002+a032)+exp(a300.1+a030+a002+a032)))),
	 fr[1,9]*(pi33.1*(1/(1+exp(a200.1+a030+a003+a033)+exp(a300.1+a030+a003+a033)))),

	 fr[1,10]*c(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) ,
	  pi12.1*(exp(a200.1+a002)/(1+exp(a200.1+a002)+exp(a300.1+a002))) ,
	  pi13.1*(exp(a200.1+a003)/(1+exp(a200.1+a003)+exp(a300.1+a003)))),
	 fr[1,11]*c(pi21.1*(exp(a200.1+a020)/(1+exp(a200.1+a020)+exp(a300.1+a020))) ,
	  pi22.1*(exp(a200.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))) ,
	  pi23.1*(exp(a200.1+a020+a003+a023)/(1+exp(a200.1+a020+a003+a023)+exp(a300.1+a020+a003+a023)))),
	 fr[1,12]*c(pi31.1*(exp(a200.1+a030)/(1+exp(a200.1+a030)+exp(a300.1+a030))) ,
	  pi32.1*(exp(a200.1+a030+a002+a032)/(1+exp(a200.1+a030+a002+a032)+exp(a300.1+a030+a002+a032))) ,
	  pi33.1*(exp(a200.1+a030+a003+a033)/(1+exp(a200.1+a030+a003+a033)+exp(a300.1+a030+a003+a033)))),

	 fr[1,13]*c(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) ,
	  pi21.1*(exp(a300.1+a020)/(1+exp(a200.1+a020)+exp(a300.1+a020))) ,
	  pi31.1*(exp(a300.1+a030)/(1+exp(a200.1+a030)+exp(a300.1+a030)))),
	 fr[1,14]*c(pi12.1*(exp(a300.1+a002)/(1+exp(a200.1+a002)+exp(a300.1+a002))) ,
	  pi22.1*(exp(a300.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))) ,
	  pi32.1*(exp(a300.1+a030+a002+a032)/(1+exp(a200.1+a030+a002+a032)+exp(a300.1+a030+a002+a032)))),
	 fr[1,15]*c(pi13.1*(exp(a300.1+a003)/(1+exp(a200.1+a003)+exp(a300.1+a003))) ,
	  pi23.1*(exp(a300.1+a020+a003+a023)/(1+exp(a200.1+a020+a003+a023)+exp(a300.1+a020+a003+a023))) ,
	  pi33.1*(exp(a300.1+a030+a003+a033)/(1+exp(a200.1+a030+a003+a033)+exp(a300.1+a030+a003+a033)))),

	 fr[2,1]*(pi11.2*(1/(1+exp(a200.2)+exp(a300.2)))),
	 fr[2,2]*(pi12.2*(1/(1+exp(a200.2+a002)+exp(a300.2+a002)))),
	 fr[2,3]*(pi13.2*(1/(1+exp(a200.2+a003)+exp(a300.2+a003)))),
	 fr[2,4]*(pi21.2*(1/(1+exp(a200.2+a020)+exp(a300.2+a020)))),
	 fr[2,5]*(pi22.2*(1/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022)))),
	 fr[2,6]*(pi23.2*(1/(1+exp(a200.2+a020+a003+a023)+exp(a300.2+a020+a003+a023)))),
	 fr[2,7]*(pi31.2*(1/(1+exp(a200.2+a030)+exp(a300.2+a030)))),
	 fr[2,8]*(pi32.2*(1/(1+exp(a200.2+a030+a002+a032)+exp(a300.2+a030+a002+a032)))),
	 fr[2,9]*(pi33.2*(1/(1+exp(a200.2+a030+a003+a033)+exp(a300.2+a030+a003+a033)))),

	 fr[2,10]*c(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) ,
	  pi12.2*(exp(a200.2+a002)/(1+exp(a200.2+a002)+exp(a300.2+a002))) ,
	  pi13.2*(exp(a200.2+a003)/(1+exp(a200.2+a003)+exp(a300.2+a003)))),
	 fr[2,11]*c(pi21.2*(exp(a200.2+a020)/(1+exp(a200.2+a020)+exp(a300.2+a020))) ,
	  pi22.2*(exp(a200.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))) ,
	  pi23.2*(exp(a200.2+a020+a003+a023)/(1+exp(a200.2+a020+a003+a023)+exp(a300.2+a020+a003+a023)))),
	 fr[2,12]*c(pi31.2*(exp(a200.2+a030)/(1+exp(a200.2+a030)+exp(a300.2+a030))) ,
	  pi32.2*(exp(a200.2+a030+a002+a032)/(1+exp(a200.2+a030+a002+a032)+exp(a300.2+a030+a002+a032))) ,
	  pi33.2*(exp(a200.2+a030+a003+a033)/(1+exp(a200.2+a030+a003+a033)+exp(a300.2+a030+a003+a033)))),

	 fr[2,13]*c(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) ,
	  pi21.2*(exp(a300.2+a020)/(1+exp(a200.2+a020)+exp(a300.2+a020))) ,
	  pi31.2*(exp(a300.2+a030)/(1+exp(a200.2+a030)+exp(a300.2+a030)))),
	 fr[2,14]*c(pi12.2*(exp(a300.2+a002)/(1+exp(a200.2+a002)+exp(a300.2+a002))) ,
	  pi22.2*(exp(a300.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))) ,
	  pi32.2*(exp(a300.2+a030+a002+a032)/(1+exp(a200.2+a030+a002+a032)+exp(a300.2+a030+a002+a032)))),
	 fr[2,15]*c(pi13.2*(exp(a300.2+a003)/(1+exp(a200.2+a003)+exp(a300.2+a003))) ,
	  pi23.2*(exp(a300.2+a020+a003+a023)/(1+exp(a200.2+a020+a003+a023)+exp(a300.2+a020+a003+a023))) ,
	  pi33.2*(exp(a300.2+a030+a003+a033)/(1+exp(a200.2+a030+a003+a033)+exp(a300.2+a030+a003+a033))))
),ncol=3,byrow=T),2) #Cenários t=3 precisam ser transpostos (linhas 7:9 e 16:18)

mnar5.analit<-mnar5.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])

loglin.mnar5<-funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar5.fit$est[1:16])),Vtheta=B%*%solve(attr(mnar5.analit,"hessian")[1,,],tol=1e-50)[1:16,1:16]%*%t(B),A1=A,X=rep(1,8))
summary(loglin.mnar5) #Logaritmos das razões de chances observadas e sob o modelo H:\omega_{ij(s)}=\omega, erros padrões e teste Wald de H, Tabela 3.5, MNAR4

funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%mnar5.fit$est[1:16])),Vtheta=B%*%solve(attr(mnar5.analit,"hessian")[1,,],tol=1e-50)[1:16,1:16]%*%t(B),
  A1=A,U=UL) #Teste de Wald de H:\omega_{12(s)}=\omega_{21(s)}=\omega_{22(s)}=0,s=1,2, p.101, MNAR4