#Comandos para reproduzir as análises do Exemplo 2, parte 2 (dados agrupados), pp.102-108, 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,36,62,25, 176,145, 28,22),
          c(120,41,47,30, 103, 83, 31,22))
padroes<-t(rep(1,2))%x%cbind(diag(4),diag(2)%x%rep(1,2), rep(1,2)%x%diag(2))
classes<-rep(1,2)%x%t(c(4,2,2))

A<-diag(2)%x%t(c(1,-1,-1,1))


#Intervalos para o melhor-pior caso apresentados na Tabela 3.6

#Limite inferior para \omega_{(1)}
log(167*25/((36+176+22)*(62+145+28)))
#Limite inferior para \omega_{(2)}
log(120*30/((41+103+22)*(47+83+31)))
#Limite superior para \omega_{(1)}
log((167+176+28)*(25+145+22)/(36*62))
#Limite superior para \omega_{(2)}
log((120+103+31)*(30+83+22)/(41*47))


catdata.acc<-readCatdata(TF=TF[,1:4])
loglin.acc<-funlinWLS(model=c("lin","log"),obj=catdata.acc,A1=A,XL=rep(1,2))
summary(loglin.acc) #Logaritmos das razões de chances observadas e sob o modelo H:\omega_{(s)}=\omega, erros padrões e teste Wald de H, Tabela 3.10, ACC


catdata<-readCatdata(TF=TF,Zp=padroes[,-c(1:4,9:12)],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.7
mcar.ml$yst #Tabela 3.8, MCAR, t=1,2,3
mar.ml$yst #Tabela 3.8, 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.8, 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.8, MAR, Total

loglin.mar<-funlinWLS(model=c("lin","log"),obj=mar.ml,A1=A,XL=rep(1,2))
summary(loglin.mar) #Logaritmos das razões de chances observadas e sob o modelo H:\omega_{(s)}=\omega, erros padrões e teste Wald de H, Tabela 3.10, 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



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

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





#menos log-verossimilhança do modelo MNAR1 (protectivo), saturado
mnar1.mlv<-function(p,o1.1=TF[1,1],o1.2=TF[1,2],o1.3=TF[1,3],o1.4=TF[1,4],o1.5=TF[1,5],o1.6=TF[1,6],o1.7=TF[1,7],o1.8=TF[1,8],
	 o2.1=TF[2,1],o2.2=TF[2,2],o2.3=TF[2,3],o2.4=TF[2,4],o2.5=TF[2,5],o2.6=TF[2,6],o2.7=TF[2,7],o2.8=TF[2,8]){
	# p=(pi11(1),...,pi22(2),a2(11),a2(21),a3(11),a3(21),a2(12),a2(22),a3(12),a3(22))
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1

	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2

	a2.1.1<-p[7]; a2.2.1<-p[8]
	a3.1.1<-p[9]; a3.2.1<-p[10]

	a2.1.2<-p[11];a2.2.2<-p[12]
	a3.1.2<-p[13];a3.2.2<-p[14]

	value<- -(
	 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(pi21.1*(1-a2.1.1-a3.2.1))+o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1-a2.2.1-a3.2.1))+

	 o1.5*log(pi11.1*a2.1.1 + pi12.1*a2.2.1)+
	 o1.6*log(pi21.1*a2.1.1 + (1-pi11.1-pi12.1-pi21.1)*a2.2.1)+

	 o1.7*log(pi11.1*a3.1.1 + pi21.1*a3.2.1)+
	 o1.8*log(pi12.1*a3.1.1 + (1-pi11.1-pi12.1-pi21.1)*a3.2.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(pi21.2*(1-a2.1.2-a3.2.2))+o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1-a2.2.2-a3.2.2))+

	 o2.5*log(pi11.2*a2.1.2 + pi12.2*a2.2.2)+
	 o2.6*log(pi21.2*a2.1.2 + (1-pi11.2-pi12.2-pi21.2)*a2.2.2)+

	 o2.7*log(pi11.2*a3.1.2 + pi21.2*a3.2.2)+
	 o2.8*log(pi12.2*a3.1.2 + (1-pi11.2-pi12.2-pi21.2)*a3.2.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(pi21.1*(1-a2.1.1-a3.2.1))+o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1-a2.2.1-a3.2.1))+

	 o1.5*log(pi11.1*a2.1.1 + pi12.1*a2.2.1)+
	 o1.6*log(pi21.1*a2.1.1 + (1-pi11.1-pi12.1-pi21.1)*a2.2.1)+

	 o1.7*log(pi11.1*a3.1.1 + pi21.1*a3.2.1)+
	 o1.8*log(pi12.1*a3.1.1 + (1-pi11.1-pi12.1-pi21.1)*a3.2.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(pi21.2*(1-a2.1.2-a3.2.2))+o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1-a2.2.2-a3.2.2))+

	 o2.5*log(pi11.2*a2.1.2 + pi12.2*a2.2.2)+
	 o2.6*log(pi21.2*a2.1.2 + (1-pi11.2-pi12.2-pi21.2)*a2.2.2)+

	 o2.7*log(pi11.2*a3.1.2 + pi21.2*a3.2.2)+
	 o2.8*log(pi12.2*a3.1.2 + (1-pi11.2-pi12.2-pi21.2)*a3.2.2)

	),c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a2.1.1","a2.2.1","a3.1.1","a3.2.1","a2.1.2","a2.2.2","a3.1.2","a3.2.2"),
	c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a2.1.1","a2.2.1","a3.1.1","a3.2.1","a2.1.2","a2.2.2","a3.1.2","a3.2.2",
	 "o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8",
	 "o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8")
)

mnar1.fit<-.nlmP(objfunc=mnar1.mlv,params=c(rep(1/4,6),rep(1/3,8)),lower=rep(0,14),upper=rep(1,14),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.7, MNAR1
p<-mnar1.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.8, MNAR1
o1.1<-o1.2<-o1.3<-o1.4<-o1.5<-o1.6<-o1.7<-o1.8<-rowSums(TF)[1]
o2.1<-o2.2<-o2.3<-o2.4<-o2.5<-o2.6<-o2.7<-o2.8<-rowSums(TF)[2]
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1
	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2
	a2.1.1<-p[7]; a2.2.1<-p[8]
	a3.1.1<-p[9]; a3.2.1<-p[10]
	a2.1.2<-p[11];a2.2.2<-p[12]
	a3.1.2<-p[13];a3.2.2<-p[14]
x<-c(
	 o1.1*c(pi11.1*(1-a2.1.1-a3.1.1)),o1.2*c(pi12.1*(1-a2.2.1-a3.1.1)),
	 o1.3*c(pi21.1*(1-a2.1.1-a3.2.1)),o1.4*c((1-pi11.1-pi12.1-pi21.1)*(1-a2.2.1-a3.2.1)),

	 o1.5*c(pi11.1*a2.1.1 , pi12.1*a2.2.1),
	 o1.6*c(pi21.1*a2.1.1 , (1-pi11.1-pi12.1-pi21.1)*a2.2.1),

	 o1.7*c(pi11.1*a3.1.1 , pi21.1*a3.2.1),
	 o1.8*c(pi12.1*a3.1.1 , (1-pi11.1-pi12.1-pi21.1)*a3.2.1),

	 o2.1*c(pi11.2*(1-a2.1.2-a3.1.2)),o2.2*c(pi12.2*(1-a2.2.2-a3.1.2)),
	 o2.3*c(pi21.2*(1-a2.1.2-a3.2.2)),o2.4*c((1-pi11.2-pi12.2-pi21.2)*(1-a2.2.2-a3.2.2)),

	 o2.5*c(pi11.2*a2.1.2 , pi12.2*a2.2.2),
	 o2.6*c(pi21.2*a2.1.2 , (1-pi11.2-pi12.2-pi21.2)*a2.2.2),

	 o2.7*c(pi11.2*a3.1.2 , pi21.2*a3.2.2),
	 o2.8*c(pi12.2*a3.1.2 , (1-pi11.2-pi12.2-pi21.2)*a3.2.2)
   )
round(rbind(c(x[1] ,x[2] ,x[5] ,x[6] , x[9],x[11]),
            c(x[3] ,x[4] ,x[7] ,x[8] ,x[10],x[12]),
            c(x[13],x[14],x[17],x[18],x[21],x[23]),
            c(x[15],x[16],x[19],x[20],x[22],x[24])),2)

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],
	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[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8])

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





#menos log-verossimilhança do modelo MNAR2, saturado
mnar2.mlv<-function(p,o1.1=TF[1,1],o1.2=TF[1,2],o1.3=TF[1,3],o1.4=TF[1,4],o1.5=TF[1,5],o1.6=TF[1,6],o1.7=TF[1,7],o1.8=TF[1,8],
	 o2.1=TF[2,1],o2.2=TF[2,2],o2.3=TF[2,3],o2.4=TF[2,4],o2.5=TF[2,5],o2.6=TF[2,6],o2.7=TF[2,7],o2.8=TF[2,8]){
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1

	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2

	a2.1.1<-p[7];a2.2.1<-p[8]
	a3.1.1<-p[9];a3.2.1<-p[10]

	a2.1.2<-p[11];a2.2.2<-p[12]
	a3.1.2<-p[13];a3.2.2<-p[14]

	value<- -(
	 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(pi21.1*(1-a2.2.1-a3.2.1))+o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1-a2.1.1-a3.1.1))+

	 o1.5*log(pi11.1*a2.1.1 + pi12.1*a2.2.1)+
	 o1.6*log(pi21.1*a2.2.1 + (1-pi11.1-pi12.1-pi21.1)*a2.1.1)+

	 o1.7*log(pi11.1*a3.1.1 + pi21.1*a3.2.1)+
	 o1.8*log(pi12.1*a3.2.1 + (1-pi11.1-pi12.1-pi21.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(pi21.2*(1-a2.2.2-a3.2.2))+o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1-a2.1.2-a3.1.2))+

	 o2.5*log(pi11.2*a2.1.2 + pi12.2*a2.2.2)+
	 o2.6*log(pi21.2*a2.2.2 + (1-pi11.2-pi12.2-pi21.2)*a2.1.2)+

	 o2.7*log(pi11.2*a3.1.2 + pi21.2*a3.2.2)+
	 o2.8*log(pi12.2*a3.2.2 + (1-pi11.2-pi12.2-pi21.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(pi21.1*(1-a2.2.1-a3.2.1))+o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1-a2.1.1-a3.1.1))+

	 o1.5*log(pi11.1*a2.1.1 + pi12.1*a2.2.1)+
	 o1.6*log(pi21.1*a2.2.1 + (1-pi11.1-pi12.1-pi21.1)*a2.1.1)+

	 o1.7*log(pi11.1*a3.1.1 + pi21.1*a3.2.1)+
	 o1.8*log(pi12.1*a3.2.1 + (1-pi11.1-pi12.1-pi21.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(pi21.2*(1-a2.2.2-a3.2.2))+o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1-a2.1.2-a3.1.2))+

	 o2.5*log(pi11.2*a2.1.2 + pi12.2*a2.2.2)+
	 o2.6*log(pi21.2*a2.2.2 + (1-pi11.2-pi12.2-pi21.2)*a2.1.2)+

	 o2.7*log(pi11.2*a3.1.2 + pi21.2*a3.2.2)+
	 o2.8*log(pi12.2*a3.2.2 + (1-pi11.2-pi12.2-pi21.2)*a3.1.2)

	),c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a2.1.1","a2.2.1","a3.1.1","a3.2.1","a2.1.2","a2.2.2","a3.1.2","a3.2.2"),
	c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a2.1.1","a2.2.1","a3.1.1","a3.2.1","a2.1.2","a2.2.2","a3.1.2","a3.2.2",
	 "o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8",
	 "o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8")
)

mnar2.fit<-.nlmP(objfunc=mnar2.mlv,params=c(rep(1/4,6),rep(1/3,8)),lower=rep(0,14),upper=rep(1,14),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.7, MNAR2
p<-mnar2.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.8, MNAR2
o1.1<-o1.2<-o1.3<-o1.4<-o1.5<-o1.6<-o1.7<-o1.8<-rowSums(TF)[1]
o2.1<-o2.2<-o2.3<-o2.4<-o2.5<-o2.6<-o2.7<-o2.8<-rowSums(TF)[2]
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1
	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2
	a2.1.1<-p[7];a2.2.1<-p[8]
	a3.1.1<-p[9];a3.2.1<-p[10]
	a2.1.2<-p[11];a2.2.2<-p[12]
	a3.1.2<-p[13];a3.2.2<-p[14]
x<-c(
	 o1.1*c(pi11.1*(1-a2.1.1-a3.1.1)),o1.2*c(pi12.1*(1-a2.2.1-a3.2.1)),
	 o1.3*c(pi21.1*(1-a2.2.1-a3.2.1)),o1.4*c((1-pi11.1-pi12.1-pi21.1)*(1-a2.1.1-a3.1.1)),

	 o1.5*c(pi11.1*a2.1.1 , pi12.1*a2.2.1),
	 o1.6*c(pi21.1*a2.2.1 , (1-pi11.1-pi12.1-pi21.1)*a2.1.1),

	 o1.7*c(pi11.1*a3.1.1 , pi21.1*a3.2.1),
	 o1.8*c(pi12.1*a3.2.1 , (1-pi11.1-pi12.1-pi21.1)*a3.1.1),

	 o2.1*c(pi11.2*(1-a2.1.2-a3.1.2)),o2.2*c(pi12.2*(1-a2.2.2-a3.2.2)),
	 o2.3*c(pi21.2*(1-a2.2.2-a3.2.2)),o2.4*c((1-pi11.2-pi12.2-pi21.2)*(1-a2.1.2-a3.1.2)),

	 o2.5*c(pi11.2*a2.1.2 , pi12.2*a2.2.2),
	 o2.6*c(pi21.2*a2.2.2 , (1-pi11.2-pi12.2-pi21.2)*a2.1.2),

	 o2.7*c(pi11.2*a3.1.2 , pi21.2*a3.2.2),
	 o2.8*c(pi12.2*a3.2.2 , (1-pi11.2-pi12.2-pi21.2)*a3.1.2)
   )
round(rbind(c(x[1] ,x[2] ,x[5] ,x[6] , x[9],x[11]),
            c(x[3] ,x[4] ,x[7] ,x[8] ,x[10],x[12]),
            c(x[13],x[14],x[17],x[18],x[21],x[23]),
            c(x[15],x[16],x[19],x[20],x[22],x[24])),2)

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],
	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[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8])

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





#menos log-verossimilhança do modelo MNAR3, 2 graus de liberdade
mnar3.mlv<-function(p,o1.1=TF[1,1],o1.2=TF[1,2],o1.3=TF[1,3],o1.4=TF[1,4],o1.5=TF[1,5],o1.6=TF[1,6],o1.7=TF[1,7],o1.8=TF[1,8],
	 o2.1=TF[2,1],o2.2=TF[2,2],o2.3=TF[2,3],o2.4=TF[2,4],o2.5=TF[2,5],o2.6=TF[2,6],o2.7=TF[2,7],o2.8=TF[2,8]){
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1

	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2

	a200.1<-p[7];a200.2<-p[8]
	a202<-p[9]

	a300.1<-p[10];a300.2<-p[11]
	a302<-p[12]

	value<- -(
	 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(pi21.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a202)+exp(a300.1+a302))))+

	 o1.5*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))))+
	 o1.6*log(pi21.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))))+

	 o1.7*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))))+
	 o1.8*log(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))))+

	 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(pi21.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a202)+exp(a300.2+a302))))+

	 o2.5*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))))+
	 o2.6*log(pi21.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))))+

	 o2.7*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))))+
	 o2.8*log(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))))

	)
	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(pi21.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a202)+exp(a300.1+a302))))+

	 o1.5*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))))+
	 o1.6*log(pi21.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302))))+

	 o1.7*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))))+
	 o1.8*log(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))))+

	 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(pi21.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a202)+exp(a300.2+a302))))+

	 o2.5*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))))+
	 o2.6*log(pi21.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302))))+

	 o2.7*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))))+
	 o2.8*log(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))))

	),c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a200.2","a202","a300.1","a300.2","a302"),
	c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a200.2","a202","a300.1","a300.2","a302",
	 "o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8",
	 "o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8")
)

mnar3.fit<-.nlmP(objfunc=mnar3.mlv,params=c(rep(1/4,6),rep(0,6)),lower=c(rep(0,6),rep(-Inf,6)),upper=c(rep(1,6),rep(Inf,6)),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.7, MNAR3
p<-mnar3.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.9, MNAR3
o1.1<-o1.2<-o1.3<-o1.4<-o1.5<-o1.6<-o1.7<-o1.8<-rowSums(TF)[1]
o2.1<-o2.2<-o2.3<-o2.4<-o2.5<-o2.6<-o2.7<-o2.8<-rowSums(TF)[2]
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1
	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2
	a200.1<-p[7];a200.2<-p[8]
	a202<-p[9]
	a300.1<-p[10];a300.2<-p[11]
	a302<-p[12]
x<-c(
	 o1.1*c(pi11.1*(1/(1+exp(a200.1)+exp(a300.1)))),
	 o1.2*c(pi12.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302)))),
	 o1.3*c(pi21.1*(1/(1+exp(a200.1)+exp(a300.1)))),
	 o1.4*c((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a202)+exp(a300.1+a302)))),

	 o1.5*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)))),
	 o1.6*c(pi21.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) ,
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a202)/(1+exp(a200.1+a202)+exp(a300.1+a302)))),

	 o1.7*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)))),
	 o1.8*c(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) ,
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302)))),

	 o2.1*c(pi11.2*(1/(1+exp(a200.2)+exp(a300.2)))),
	 o2.2*c(pi12.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302)))),
	 o2.3*c(pi21.2*(1/(1+exp(a200.2)+exp(a300.2)))),
	 o2.4*c((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a202)+exp(a300.2+a302)))),

	 o2.5*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)))),
	 o2.6*c(pi21.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) ,
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a202)/(1+exp(a200.2+a202)+exp(a300.2+a302)))),

	 o2.7*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)))),
	 o2.8*c(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) ,
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))))
   )
round(rbind(c(x[1] ,x[2] ,x[5] ,x[6] , x[9],x[11]),
            c(x[3] ,x[4] ,x[7] ,x[8] ,x[10],x[12]),
            c(x[13],x[14],x[17],x[18],x[21],x[23]),
            c(x[15],x[16],x[19],x[20],x[22],x[24])),2)

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],
	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[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8])

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





#menos log-verossimilhança do modelo MNAR4, saturado
mnar4.mlv<-function(p,o1.1=TF[1,1],o1.2=TF[1,2],o1.3=TF[1,3],o1.4=TF[1,4],o1.5=TF[1,5],o1.6=TF[1,6],o1.7=TF[1,7],o1.8=TF[1,8],
	 o2.1=TF[2,1],o2.2=TF[2,2],o2.3=TF[2,3],o2.4=TF[2,4],o2.5=TF[2,5],o2.6=TF[2,6],o2.7=TF[2,7],o2.8=TF[2,8]){
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1

	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2

	a200.1<-p[7];a200.2<-p[8]
	a220<-p[9];a202<-p[10]

	a300.1<-p[11];a300.2<-p[12]
	a320<-p[13];a302<-p[14]

	value<- -(
	 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(pi21.1*(1/(1+exp(a200.1+a220)+exp(a300.1+a320))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))))+

	 o1.5*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))))+
	 o1.6*log(pi21.1*(exp(a200.1+a220)/(1+exp(a200.1+a220)+exp(a300.1+a320))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a220+a202)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))))+

	 o1.7*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))))+
	 o1.8*log(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a320+a302)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))))+

	 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(pi21.2*(1/(1+exp(a200.2+a220)+exp(a300.2+a320))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))))+

	 o2.5*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))))+
	 o2.6*log(pi21.2*(exp(a200.2+a220)/(1+exp(a200.2+a220)+exp(a300.2+a320))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a220+a202)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))))+

	 o2.7*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))))+
	 o2.8*log(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a320+a302)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))))

	)
	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(pi21.1*(1/(1+exp(a200.1+a220)+exp(a300.1+a320))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))))+

	 o1.5*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))))+
	 o1.6*log(pi21.1*(exp(a200.1+a220)/(1+exp(a200.1+a220)+exp(a300.1+a320))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a220+a202)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))))+

	 o1.7*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))))+
	 o1.8*log(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a320+a302)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302))))+

	 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(pi21.2*(1/(1+exp(a200.2+a220)+exp(a300.2+a320))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))))+

	 o2.5*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))))+
	 o2.6*log(pi21.2*(exp(a200.2+a220)/(1+exp(a200.2+a220)+exp(a300.2+a320))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a220+a202)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))))+

	 o2.7*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))))+
	 o2.8*log(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a320+a302)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))))

	),c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a200.2","a220","a202","a300.1","a300.2","a320","a302"),
	c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a200.2","a220","a202","a300.1","a300.2","a320","a302",
	 "o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8",
	 "o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8")
)

mnar4.fit<-.nlmP(objfunc=mnar4.mlv,params=c(rep(1/4,6),rep(0,8)),lower=c(rep(0,6),rep(-Inf,8)),upper=c(rep(1,6),rep(Inf,8)),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.7, MNAR4
p<-mnar4.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.9, MNAR4
o1.1<-o1.2<-o1.3<-o1.4<-o1.5<-o1.6<-o1.7<-o1.8<-rowSums(TF)[1]
o2.1<-o2.2<-o2.3<-o2.4<-o2.5<-o2.6<-o2.7<-o2.8<-rowSums(TF)[2]
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1
	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2
	a200.1<-p[7];a200.2<-p[8]
	a220<-p[9];a202<-p[10]
	a300.1<-p[11];a300.2<-p[12]
	a320<-p[13];a302<-p[14]
x<-c(
	 o1.1*c(pi11.1*(1/(1+exp(a200.1)+exp(a300.1)))),
	 o1.2*c(pi12.1*(1/(1+exp(a200.1+a202)+exp(a300.1+a302)))),
	 o1.3*c(pi21.1*(1/(1+exp(a200.1+a220)+exp(a300.1+a320)))),
	 o1.4*c((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302)))),

	 o1.5*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)))),
	 o1.6*c(pi21.1*(exp(a200.1+a220)/(1+exp(a200.1+a220)+exp(a300.1+a320))) ,
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a220+a202)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302)))),

	 o1.7*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)))),
	 o1.8*c(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) ,
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a320+a302)/(1+exp(a200.1+a220+a202)+exp(a300.1+a320+a302)))),

	 o2.1*c(pi11.2*(1/(1+exp(a200.2)+exp(a300.2)))),
	 o2.2*c(pi12.2*(1/(1+exp(a200.2+a202)+exp(a300.2+a302)))),
	 o2.3*c(pi21.2*(1/(1+exp(a200.2+a220)+exp(a300.2+a320)))),
	 o2.4*c((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302)))),

	 o2.5*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)))),
	 o2.6*c(pi21.2*(exp(a200.2+a220)/(1+exp(a200.2+a220)+exp(a300.2+a320))) ,
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a220+a202)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302)))),

	 o2.7*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)))),
	 o2.8*c(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) ,
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a320+a302)/(1+exp(a200.2+a220+a202)+exp(a300.2+a320+a302))))
   )
round(rbind(c(x[1] ,x[2] ,x[5] ,x[6] , x[9],x[11]),
            c(x[3] ,x[4] ,x[7] ,x[8] ,x[10],x[12]),
            c(x[13],x[14],x[17],x[18],x[21],x[23]),
            c(x[15],x[16],x[19],x[20],x[22],x[24])),2)

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],
	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[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8])

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





#menos log-verossimilhança do modelo MNAR5, 1 grau de liberdade
mnar5.mlv<-function(p,o1.1=TF[1,1],o1.2=TF[1,2],o1.3=TF[1,3],o1.4=TF[1,4],o1.5=TF[1,5],o1.6=TF[1,6],o1.7=TF[1,7],o1.8=TF[1,8],
	 o2.1=TF[2,1],o2.2=TF[2,2],o2.3=TF[2,3],o2.4=TF[2,4],o2.5=TF[2,5],o2.6=TF[2,6],o2.7=TF[2,7],o2.8=TF[2,8]){
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1

	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2

	a200.1<-p[7];a300.1<-p[8]
	a020<-p[9];a002<-p[10]
	a022<-p[11]
	a200.2<-p[12];a300.2<-p[13]

	value<- -(
	 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(pi21.1*(1/(1+exp(a200.1+a020)+exp(a300.1+a020))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))))+

	 o1.5*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))))+
	 o1.6*log(pi21.1*(exp(a200.1+a020)/(1+exp(a200.1+a020)+exp(a300.1+a020))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))))+

	 o1.7*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))))+
	 o1.8*log(pi12.1*(exp(a300.1+a002)/(1+exp(a200.1+a002)+exp(a300.1+a002))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))))+

	 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(pi21.2*(1/(1+exp(a200.2+a020)+exp(a300.2+a020))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))))+

	 o2.5*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))))+
	 o2.6*log(pi21.2*(exp(a200.2+a020)/(1+exp(a200.2+a020)+exp(a300.2+a020))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))))+

	 o2.7*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))))+
	 o2.8*log(pi12.2*(exp(a300.2+a002)/(1+exp(a200.2+a002)+exp(a300.2+a002))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))))

	)
	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(pi21.1*(1/(1+exp(a200.1+a020)+exp(a300.1+a020))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))))+

	 o1.5*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))))+
	 o1.6*log(pi21.1*(exp(a200.1+a020)/(1+exp(a200.1+a020)+exp(a300.1+a020))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))))+

	 o1.7*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))))+
	 o1.8*log(pi12.1*(exp(a300.1+a002)/(1+exp(a200.1+a002)+exp(a300.1+a002))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022))))+

	 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(pi21.2*(1/(1+exp(a200.2+a020)+exp(a300.2+a020))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))))+

	 o2.5*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))))+
	 o2.6*log(pi21.2*(exp(a200.2+a020)/(1+exp(a200.2+a020)+exp(a300.2+a020))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))))+

	 o2.7*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))))+
	 o2.8*log(pi12.2*(exp(a300.2+a002)/(1+exp(a200.2+a002)+exp(a300.2+a002))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))))

	),c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a300.1","a020","a002","a022","a200.2","a300.2"),
	c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a300.1","a020","a002","a022","a200.2","a300.2",
	 "o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8",
	 "o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8")
)

mnar5.fit<-.nlmP(objfunc=mnar5.mlv,params=c(rep(1/4,6),rep(0,7)),lower=c(rep(0,6),rep(-Inf,7)),upper=c(rep(1,6),rep(Inf,7)),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.7, MNAR5
p<-mnar5.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.9, MNAR5
o1.1<-o1.2<-o1.3<-o1.4<-o1.5<-o1.6<-o1.7<-o1.8<-rowSums(TF)[1]
o2.1<-o2.2<-o2.3<-o2.4<-o2.5<-o2.6<-o2.7<-o2.8<-rowSums(TF)[2]
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1
	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2
	a200.1<-p[7];a300.1<-p[8]
	a020<-p[9];a002<-p[10]
	a022<-p[11]
	a200.2<-p[12];a300.2<-p[13]
x<-c(
	 o1.1*c(pi11.1*(1/(1+exp(a200.1)+exp(a300.1)))),
	 o1.2*c(pi12.1*(1/(1+exp(a200.1+a002)+exp(a300.1+a002)))),
	 o1.3*c(pi21.1*(1/(1+exp(a200.1+a020)+exp(a300.1+a020)))),
	 o1.4*c((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022)))),

	 o1.5*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)))),
	 o1.6*c(pi21.1*(exp(a200.1+a020)/(1+exp(a200.1+a020)+exp(a300.1+a020))) ,
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022)))),

	 o1.7*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)))),
	 o1.8*c(pi12.1*(exp(a300.1+a002)/(1+exp(a200.1+a002)+exp(a300.1+a002))) ,
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a020+a002+a022)/(1+exp(a200.1+a020+a002+a022)+exp(a300.1+a020+a002+a022)))),

	 o2.1*c(pi11.2*(1/(1+exp(a200.2)+exp(a300.2)))),
	 o2.2*c(pi12.2*(1/(1+exp(a200.2+a002)+exp(a300.2+a002)))),
	 o2.3*c(pi21.2*(1/(1+exp(a200.2+a020)+exp(a300.2+a020)))),
	 o2.4*c((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022)))),

	 o2.5*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)))),
	 o2.6*c(pi21.2*(exp(a200.2+a020)/(1+exp(a200.2+a020)+exp(a300.2+a020))) ,
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022)))),

	 o2.7*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)))),
	 o2.8*c(pi12.2*(exp(a300.2+a002)/(1+exp(a200.2+a002)+exp(a300.2+a002))) ,
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a020+a002+a022)/(1+exp(a200.2+a020+a002+a022)+exp(a300.2+a020+a002+a022))))
   )
round(rbind(c(x[1] ,x[2] ,x[5] ,x[6] , x[9],x[11]),
            c(x[3] ,x[4] ,x[7] ,x[8] ,x[10],x[12]),
            c(x[13],x[14],x[17],x[18],x[21],x[23]),
            c(x[15],x[16],x[19],x[20],x[22],x[24])),2)

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],
	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[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8])

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





#menos log-verossimilhança do modelo MNAR6, saturado
mnar6.mlv<-function(p,o1.1=TF[1,1],o1.2=TF[1,2],o1.3=TF[1,3],o1.4=TF[1,4],o1.5=TF[1,5],o1.6=TF[1,6],o1.7=TF[1,7],o1.8=TF[1,8],
	 o2.1=TF[2,1],o2.2=TF[2,2],o2.3=TF[2,3],o2.4=TF[2,4],o2.5=TF[2,5],o2.6=TF[2,6],o2.7=TF[2,7],o2.8=TF[2,8]){
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1

	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2

	a200.1<-p[7];a200.2<-p[8]
	a020.1<-p[9];a020.2<-p[10]

	a300.1<-p[11];a300.2<-p[12]
	a002.1<-p[13];a002.2<-p[14]

	value<- -(
	 o1.1*log(pi11.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 o1.2*log(pi12.1*(1/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))))+
	 o1.3*log(pi21.1*(1/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a020.1+a002.1)+exp(a300.1+a020.1+a002.1))))+

	 o1.5*log(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi12.1*(exp(a200.1+a002.1)/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))))+
	 o1.6*log(pi21.1*(exp(a200.1+a020.1)/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a020.1+a002.1)/(1+exp(a200.1+a020.1+a002.1)+exp(a300.1+a020.1+a002.1))))+

	 o1.7*log(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi21.1*(exp(a300.1+a020.1)/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))))+
	 o1.8*log(pi12.1*(exp(a300.1+a002.1)/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a020.1+a002.1)/(1+exp(a200.1+a020.1+a002.1)+exp(a300.1+a020.1+a002.1))))+

	 o2.1*log(pi11.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.2*log(pi12.2*(1/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))))+
	 o2.3*log(pi21.2*(1/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a020.2+a002.2)+exp(a300.2+a020.2+a002.2))))+

	 o2.5*log(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi12.2*(exp(a200.2+a002.2)/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))))+
	 o2.6*log(pi21.2*(exp(a200.2+a020.2)/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a020.2+a002.2)/(1+exp(a200.2+a020.2+a002.2)+exp(a300.2+a020.2+a002.2))))+

	 o2.7*log(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi21.2*(exp(a300.2+a020.2)/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))))+
	 o2.8*log(pi12.2*(exp(a300.2+a002.2)/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a020.2+a002.2)/(1+exp(a200.2+a020.2+a002.2)+exp(a300.2+a020.2+a002.2))))

	)
	value
}

mnar6.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.1)+exp(a300.1+a002.1))))+
	 o1.3*log(pi21.1*(1/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a020.1+a002.1)+exp(a300.1+a020.1+a002.1))))+

	 o1.5*log(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi12.1*(exp(a200.1+a002.1)/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))))+
	 o1.6*log(pi21.1*(exp(a200.1+a020.1)/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a020.1+a002.1)/(1+exp(a200.1+a020.1+a002.1)+exp(a300.1+a020.1+a002.1))))+

	 o1.7*log(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi21.1*(exp(a300.1+a020.1)/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))))+
	 o1.8*log(pi12.1*(exp(a300.1+a002.1)/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a020.1+a002.1)/(1+exp(a200.1+a020.1+a002.1)+exp(a300.1+a020.1+a002.1))))+

	 o2.1*log(pi11.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.2*log(pi12.2*(1/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))))+
	 o2.3*log(pi21.2*(1/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a020.2+a002.2)+exp(a300.2+a020.2+a002.2))))+

	 o2.5*log(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi12.2*(exp(a200.2+a002.2)/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))))+
	 o2.6*log(pi21.2*(exp(a200.2+a020.2)/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a020.2+a002.2)/(1+exp(a200.2+a020.2+a002.2)+exp(a300.2+a020.2+a002.2))))+

	 o2.7*log(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi21.2*(exp(a300.2+a020.2)/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))))+
	 o2.8*log(pi12.2*(exp(a300.2+a002.2)/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a020.2+a002.2)/(1+exp(a200.2+a020.2+a002.2)+exp(a300.2+a020.2+a002.2))))

	),c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a200.2","a020.1","a020.2","a300.1","a300.2","a002.1","a002.2"),
	c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a200.2","a020.1","a020.2","a300.1","a300.2","a002.1","a002.2",
	 "o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8",
	 "o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8")
)

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

#Freqüências ampliadas esperadas estimadas, Tabela 3.9, MNAR6
o1.1<-o1.2<-o1.3<-o1.4<-o1.5<-o1.6<-o1.7<-o1.8<-rowSums(TF)[1]
o2.1<-o2.2<-o2.3<-o2.4<-o2.5<-o2.6<-o2.7<-o2.8<-rowSums(TF)[2]
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1
	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2
	a200.1<-p[7];a200.2<-p[8]
	a020.1<-p[9];a020.2<-p[10]
	a300.1<-p[11];a300.2<-p[12]
	a002.1<-p[13];a002.2<-p[14]
x<-c(
	 o1.1*c(pi11.1*(1/(1+exp(a200.1)+exp(a300.1)))),
	 o1.2*c(pi12.1*(1/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1)))),
	 o1.3*c(pi21.1*(1/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1)))),
	 o1.4*c((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a020.1+a002.1)+exp(a300.1+a020.1+a002.1)))),

	 o1.5*c(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) ,
	  pi12.1*(exp(a200.1+a002.1)/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1)))),
	 o1.6*c(pi21.1*(exp(a200.1+a020.1)/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))) ,
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a020.1+a002.1)/(1+exp(a200.1+a020.1+a002.1)+exp(a300.1+a020.1+a002.1)))),

	 o1.7*c(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) ,
	  pi21.1*(exp(a300.1+a020.1)/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1)))),
	 o1.8*c(pi12.1*(exp(a300.1+a002.1)/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))) ,
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a020.1+a002.1)/(1+exp(a200.1+a020.1+a002.1)+exp(a300.1+a020.1+a002.1)))),

	 o2.1*c(pi11.2*(1/(1+exp(a200.2)+exp(a300.2)))),
	 o2.2*c(pi12.2*(1/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2)))),
	 o2.3*c(pi21.2*(1/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2)))),
	 o2.4*c((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a020.2+a002.2)+exp(a300.2+a020.2+a002.2)))),

	 o2.5*c(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) ,
	  pi12.2*(exp(a200.2+a002.2)/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2)))),
	 o2.6*c(pi21.2*(exp(a200.2+a020.2)/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))) ,
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a020.2+a002.2)/(1+exp(a200.2+a020.2+a002.2)+exp(a300.2+a020.2+a002.2)))),

	 o2.7*c(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) ,
	  pi21.2*(exp(a300.2+a020.2)/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2)))),
	 o2.8*c(pi12.2*(exp(a300.2+a002.2)/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))) ,
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a020.2+a002.2)/(1+exp(a200.2+a020.2+a002.2)+exp(a300.2+a020.2+a002.2))))
   )
round(rbind(c(x[1] ,x[2] ,x[5] ,x[6] , x[9],x[11]),
            c(x[3] ,x[4] ,x[7] ,x[8] ,x[10],x[12]),
            c(x[13],x[14],x[17],x[18],x[21],x[23]),
            c(x[15],x[16],x[19],x[20],x[22],x[24])),2)

mnar6.analit<-mnar6.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],
	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[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8])

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





#menos log-verossimilhança do modelo MNAR7, saturado
mnar7.mlv<-function(p,o1.1=TF[1,1],o1.2=TF[1,2],o1.3=TF[1,3],o1.4=TF[1,4],o1.5=TF[1,5],o1.6=TF[1,6],o1.7=TF[1,7],o1.8=TF[1,8],
	 o2.1=TF[2,1],o2.2=TF[2,2],o2.3=TF[2,3],o2.4=TF[2,4],o2.5=TF[2,5],o2.6=TF[2,6],o2.7=TF[2,7],o2.8=TF[2,8]){
	# p=(pi11(1),...,pi22(2),a2(11),a2(12),a2(21),a2(22),a3(11),a3(12),a3(21),a3(22))
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1

	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2

	a2.11<-p[7]; a2.12<-p[8]
	a2.21<-p[9]; a2.22<-p[10]

	a3.11<-p[11];a3.12<-p[12]
	a3.21<-p[13];a3.22<-p[14]

	value<- -(
	 o1.1*log(pi11.1*(1-a2.11-a3.11))+o1.2*log(pi12.1*(1-a2.12-a3.12))+
	 o1.3*log(pi21.1*(1-a2.21-a3.21))+o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1-a2.22-a3.22))+

	 o1.5*log(pi11.1*a2.11 + pi12.1*a2.12)+
	 o1.6*log(pi21.1*a2.21 + (1-pi11.1-pi12.1-pi21.1)*a2.22)+

	 o1.7*log(pi11.1*a3.11 + pi21.1*a3.21)+
	 o1.8*log(pi12.1*a3.12 + (1-pi11.1-pi12.1-pi21.1)*a3.22)+

	 o2.1*log(pi11.2*(1-a2.11-a3.11))+o2.2*log(pi12.2*(1-a2.12-a3.12))+
	 o2.3*log(pi21.2*(1-a2.21-a3.21))+o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1-a2.22-a3.22))+

	 o2.5*log(pi11.2*a2.11 + pi12.2*a2.12)+
	 o2.6*log(pi21.2*a2.21 + (1-pi11.2-pi12.2-pi21.2)*a2.22)+

	 o2.7*log(pi11.2*a3.11 + pi21.2*a3.21)+
	 o2.8*log(pi12.2*a3.12 + (1-pi11.2-pi12.2-pi21.2)*a3.22)

	)
	value
}

mnar7.der<-deriv3(~-(
	 o1.1*log(pi11.1*(1-a2.11-a3.11))+o1.2*log(pi12.1*(1-a2.12-a3.12))+
	 o1.3*log(pi21.1*(1-a2.21-a3.21))+o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1-a2.22-a3.22))+

	 o1.5*log(pi11.1*a2.11 + pi12.1*a2.12)+
	 o1.6*log(pi21.1*a2.21 + (1-pi11.1-pi12.1-pi21.1)*a2.22)+

	 o1.7*log(pi11.1*a3.11 + pi21.1*a3.21)+
	 o1.8*log(pi12.1*a3.12 + (1-pi11.1-pi12.1-pi21.1)*a3.22)+

	 o2.1*log(pi11.2*(1-a2.11-a3.11))+o2.2*log(pi12.2*(1-a2.12-a3.12))+
	 o2.3*log(pi21.2*(1-a2.21-a3.21))+o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1-a2.22-a3.22))+

	 o2.5*log(pi11.2*a2.11 + pi12.2*a2.12)+
	 o2.6*log(pi21.2*a2.21 + (1-pi11.2-pi12.2-pi21.2)*a2.22)+

	 o2.7*log(pi11.2*a3.11 + pi21.2*a3.21)+
	 o2.8*log(pi12.2*a3.12 + (1-pi11.2-pi12.2-pi21.2)*a3.22)

	),c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a2.11","a2.12","a2.21","a2.22","a3.11","a3.12","a3.21","a3.22"),
	c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a2.11","a2.12","a2.21","a2.22","a3.11","a3.12","a3.21","a3.22",
	 "o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8",
	 "o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8")
)

mnar7.fit<-.nlmP(objfunc=mnar7.mlv,params=c(rep(1/4,6),rep(1/3,8)),lower=rep(0,14),upper=rep(1,14),hessian=TRUE,print.level=0,iterlim=500)
-2*(-mnar7.fit$min-sat.lv) #Estatística de razão de verossimilhanças do ajuste do modelo, Tabela 3.7, MNAR7
p<-mnar7.fit$est

#Freqüências ampliadas esperadas estimadas, Tabela 3.8, MNAR7
o1.1<-o1.2<-o1.3<-o1.4<-o1.5<-o1.6<-o1.7<-o1.8<-rowSums(TF)[1]
o2.1<-o2.2<-o2.3<-o2.4<-o2.5<-o2.6<-o2.7<-o2.8<-rowSums(TF)[2]
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1
	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2
	a2.11<-p[7]; a2.12<-p[8]
	a2.21<-p[9]; a2.22<-p[10]
	a3.11<-p[11];a3.12<-p[12]
	a3.21<-p[13];a3.22<-p[14]
x<-c(
	 o1.1*c(pi11.1*(1-a2.11-a3.11)),o1.2*c(pi12.1*(1-a2.12-a3.12)),
	 o1.3*c(pi21.1*(1-a2.21-a3.21)),o1.4*c((1-pi11.1-pi12.1-pi21.1)*(1-a2.22-a3.22)),

	 o1.5*c(pi11.1*a2.11 , pi12.1*a2.12),
	 o1.6*c(pi21.1*a2.21 , (1-pi11.1-pi12.1-pi21.1)*a2.22),

	 o1.7*c(pi11.1*a3.11 , pi21.1*a3.21),
	 o1.8*c(pi12.1*a3.12 , (1-pi11.1-pi12.1-pi21.1)*a3.22),

	 o2.1*c(pi11.2*(1-a2.11-a3.11)),o2.2*c(pi12.2*(1-a2.12-a3.12)),
	 o2.3*c(pi21.2*(1-a2.21-a3.21)),o2.4*c((1-pi11.2-pi12.2-pi21.2)*(1-a2.22-a3.22)),

	 o2.5*c(pi11.2*a2.11 , pi12.2*a2.12),
	 o2.6*c(pi21.2*a2.21 , (1-pi11.2-pi12.2-pi21.2)*a2.22),

	 o2.7*c(pi11.2*a3.11 , pi21.2*a3.21),
	 o2.8*c(pi12.2*a3.12 , (1-pi11.2-pi12.2-pi21.2)*a3.22)
   )
round(rbind(c(x[1] ,x[2] ,x[5] ,x[6] , x[9],x[11]),
            c(x[3] ,x[4] ,x[7] ,x[8] ,x[10],x[12]),
            c(x[13],x[14],x[17],x[18],x[21],x[23]),
            c(x[15],x[16],x[19],x[20],x[22],x[24])),2)

p[p<0.0001]<-0.005
mnar7.analit<-mnar7.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],
	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[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8])

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





#menos log-verossimilhança da estimação do modelo MNAR8, sobresaturado (1 parâmetro a mais)
mnar8.mlv<-function(p,o1.1=TF[1,1],o1.2=TF[1,2],o1.3=TF[1,3],o1.4=TF[1,4],o1.5=TF[1,5],o1.6=TF[1,6],o1.7=TF[1,7],o1.8=TF[1,8],
	 o2.1=TF[2,1],o2.2=TF[2,2],o2.3=TF[2,3],o2.4=TF[2,4],o2.5=TF[2,5],o2.6=TF[2,6],o2.7=TF[2,7],o2.8=TF[2,8],
   a022=alfa022){
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1

	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2

	a200.1<-p[7];a200.2<-p[8]
	a020.1<-p[9];a020.2<-p[10]

	a300.1<-p[11];a300.2<-p[12]
	a002.1<-p[13];a002.2<-p[14]

	value<- -(
	 o1.1*log(pi11.1*(1/(1+exp(a200.1)+exp(a300.1))))+
	 o1.2*log(pi12.1*(1/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))))+
	 o1.3*log(pi21.1*(1/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a020.1+a002.1+a022)+exp(a300.1+a020.1+a002.1+a022))))+

	 o1.5*log(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi12.1*(exp(a200.1+a002.1)/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))))+
	 o1.6*log(pi21.1*(exp(a200.1+a020.1)/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a020.1+a002.1+a022)/(1+exp(a200.1+a020.1+a002.1+a022)+exp(a300.1+a020.1+a002.1+a022))))+

	 o1.7*log(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi21.1*(exp(a300.1+a020.1)/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))))+
	 o1.8*log(pi12.1*(exp(a300.1+a002.1)/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a020.1+a002.1+a022)/(1+exp(a200.1+a020.1+a002.1+a022)+exp(a300.1+a020.1+a002.1+a022))))+

	 o2.1*log(pi11.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.2*log(pi12.2*(1/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))))+
	 o2.3*log(pi21.2*(1/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a020.2+a002.2+a022)+exp(a300.2+a020.2+a002.2+a022))))+

	 o2.5*log(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi12.2*(exp(a200.2+a002.2)/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))))+
	 o2.6*log(pi21.2*(exp(a200.2+a020.2)/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a020.2+a002.2+a022)/(1+exp(a200.2+a020.2+a002.2+a022)+exp(a300.2+a020.2+a002.2+a022))))+

	 o2.7*log(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi21.2*(exp(a300.2+a020.2)/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))))+
	 o2.8*log(pi12.2*(exp(a300.2+a002.2)/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a020.2+a002.2+a022)/(1+exp(a200.2+a020.2+a002.2+a022)+exp(a300.2+a020.2+a002.2+a022))))

	)
	value
}

mnar8.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.1)+exp(a300.1+a002.1))))+
	 o1.3*log(pi21.1*(1/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a020.1+a002.1+a022)+exp(a300.1+a020.1+a002.1+a022))))+

	 o1.5*log(pi11.1*(exp(a200.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi12.1*(exp(a200.1+a002.1)/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))))+
	 o1.6*log(pi21.1*(exp(a200.1+a020.1)/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a020.1+a002.1+a022)/(1+exp(a200.1+a020.1+a002.1+a022)+exp(a300.1+a020.1+a002.1+a022))))+

	 o1.7*log(pi11.1*(exp(a300.1)/(1+exp(a200.1)+exp(a300.1))) +
	  pi21.1*(exp(a300.1+a020.1)/(1+exp(a200.1+a020.1)+exp(a300.1+a020.1))))+
	 o1.8*log(pi12.1*(exp(a300.1+a002.1)/(1+exp(a200.1+a002.1)+exp(a300.1+a002.1))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a020.1+a002.1+a022)/(1+exp(a200.1+a020.1+a002.1+a022)+exp(a300.1+a020.1+a002.1+a022))))+

	 o2.1*log(pi11.2*(1/(1+exp(a200.2)+exp(a300.2))))+
	 o2.2*log(pi12.2*(1/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))))+
	 o2.3*log(pi21.2*(1/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a020.2+a002.2+a022)+exp(a300.2+a020.2+a002.2+a022))))+

	 o2.5*log(pi11.2*(exp(a200.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi12.2*(exp(a200.2+a002.2)/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))))+
	 o2.6*log(pi21.2*(exp(a200.2+a020.2)/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a020.2+a002.2+a022)/(1+exp(a200.2+a020.2+a002.2+a022)+exp(a300.2+a020.2+a002.2+a022))))+

	 o2.7*log(pi11.2*(exp(a300.2)/(1+exp(a200.2)+exp(a300.2))) +
	  pi21.2*(exp(a300.2+a020.2)/(1+exp(a200.2+a020.2)+exp(a300.2+a020.2))))+
	 o2.8*log(pi12.2*(exp(a300.2+a002.2)/(1+exp(a200.2+a002.2)+exp(a300.2+a002.2))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a020.2+a002.2+a022)/(1+exp(a200.2+a020.2+a002.2+a022)+exp(a300.2+a020.2+a002.2+a022))))

	),c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a200.2","a020.1","a020.2","a300.1","a300.2","a002.1","a002.2"),
	c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a200.2","a020.1","a020.2","a300.1","a300.2","a002.1","a002.2",
	 "o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8",
	 "o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8","a022")
)

alfas<- seq(-5,5,0.01)
nalfas<-length(alfas)
z<-qnorm(0.975)

Qvs<-numeric(nalfas)
ws<-numeric(nalfas)
vars<-numeric(nalfas)
pvs<-numeric(nalfas)

date()
for(i in 1:nalfas){ #Quase 6min. num Intel Core 2 Duo 1.7GHz
  alfa022<-alfas[i]
  mnar8.fit<-.nlmP(objfunc=mnar8.mlv,params=c(rep(1/4,6),rep(0,8)),lower=c(rep(0,6),rep(-Inf,8)),upper=c(rep(1,6),rep(Inf,8)),hessian=TRUE,print.level=0,iterlim=500)
  Qvs[i]<- -2*(-mnar8.fit$min-sat.lv)
  p<-mnar8.fit$est
  mnar8.analit<-mnar8.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],
  	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[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8],alfa022)
  fit<-funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%p[1:6])),Vtheta=B%*%solve(attr(mnar8.analit,"hessian")[1,,])[1:6,1:6]%*%t(B),A1=A,X=rep(1,2))
  ws[i]<-fit$beta
  vars[i]<-c(fit$Vbeta)
  pvs[i]<-1-pchisq(fit$QwH,1)
}
date()

paste("Intervalo de ignorância: (",round(min(ws),2),";",round(max(ws),2),")",sep="") #(-1.62;2.78)

winf<-ws-z*sqrt(vars)
wsup<-ws+z*sqrt(vars)
paste("Intervalo de incerteza: (",round(min(winf),2),";",round(max(wsup),2),")",sep="") #(-2;3.16)

plot(alfas,ws,xlab=expression(paste(alpha[paste("022")])), ylab=expression(paste(omega,"'")),pch=16,type="n",ylim=c(min(winf),max(wsup)))
lines(alfas,ws,lwd=1)
lines(alfas,winf,lwd=3)
lines(alfas,wsup,lwd=3)
abline(0,0,lty=2)

plot(alfas,Qvs,type="l",lwd=3)
paste("Min/Max Qvs: ",round(min(Qvs),2),"/",round(max(Qvs),2),sep="") #0/0

plot(alfas,pvs,type="l",lwd=3)
paste("Min/Max pvs: ",round(min(pvs),3),"/",round(max(pvs),3),sep="") #0.204/1





#menos log-verossimilhança do modelo MNAR9, sobresaturado (1 parâmetro a mais)
mnar9.mlv<-function(p,o1.1=TF[1,1],o1.2=TF[1,2],o1.3=TF[1,3],o1.4=TF[1,4],o1.5=TF[1,5],o1.6=TF[1,6],o1.7=TF[1,7],o1.8=TF[1,8],
	 o2.1=TF[2,1],o2.2=TF[2,2],o2.3=TF[2,3],o2.4=TF[2,4],o2.5=TF[2,5],o2.6=TF[2,6],o2.7=TF[2,7],o2.8=TF[2,8],
   a022=alfa022){
	pi11.1<-p[1]; pi12.1<-p[2]
	pi21.1<-p[3] #pi22.1=1-pi11.1-pi12.1-pi21.1

	pi11.2<-p[4]; pi12.2<-p[5]
	pi21.2<-p[6] #pi22.2=1-pi11.2-pi12.2-pi21.2

	a200.1<-p[7];a200.2<-p[8]
	a220<-p[9];a202<-p[10]

	a300.1<-p[11];a300.2<-p[12]
	a320<-p[13];a302<-p[14]

	value<- -(
	 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(pi21.1*(1/(1+exp(a200.1+a220)+exp(a300.1+a320))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a220+a202+a022)+exp(a300.1+a320+a302+a022))))+

	 o1.5*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))))+
	 o1.6*log(pi21.1*(exp(a200.1+a220)/(1+exp(a200.1+a220)+exp(a300.1+a320))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a220+a202+a022)/(1+exp(a200.1+a220+a202+a022)+exp(a300.1+a320+a302+a022))))+

	 o1.7*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))))+
	 o1.8*log(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a320+a302+a022)/(1+exp(a200.1+a220+a202+a022)+exp(a300.1+a320+a302+a022))))+

	 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(pi21.2*(1/(1+exp(a200.2+a220)+exp(a300.2+a320))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a220+a202+a022)+exp(a300.2+a320+a302+a022))))+

	 o2.5*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))))+
	 o2.6*log(pi21.2*(exp(a200.2+a220)/(1+exp(a200.2+a220)+exp(a300.2+a320))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a220+a202+a022)/(1+exp(a200.2+a220+a202+a022)+exp(a300.2+a320+a302+a022))))+

	 o2.7*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))))+
	 o2.8*log(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a320+a302+a022)/(1+exp(a200.2+a220+a202+a022)+exp(a300.2+a320+a302+a022))))

	)
	value
}

mnar9.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(pi21.1*(1/(1+exp(a200.1+a220)+exp(a300.1+a320))))+
	 o1.4*log((1-pi11.1-pi12.1-pi21.1)*(1/(1+exp(a200.1+a220+a202+a022)+exp(a300.1+a320+a302+a022))))+

	 o1.5*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))))+
	 o1.6*log(pi21.1*(exp(a200.1+a220)/(1+exp(a200.1+a220)+exp(a300.1+a320))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a200.1+a220+a202+a022)/(1+exp(a200.1+a220+a202+a022)+exp(a300.1+a320+a302+a022))))+

	 o1.7*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))))+
	 o1.8*log(pi12.1*(exp(a300.1+a302)/(1+exp(a200.1+a202)+exp(a300.1+a302))) +
	  (1-pi11.1-pi12.1-pi21.1)*(exp(a300.1+a320+a302+a022)/(1+exp(a200.1+a220+a202+a022)+exp(a300.1+a320+a302+a022))))+

	 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(pi21.2*(1/(1+exp(a200.2+a220)+exp(a300.2+a320))))+
	 o2.4*log((1-pi11.2-pi12.2-pi21.2)*(1/(1+exp(a200.2+a220+a202+a022)+exp(a300.2+a320+a302+a022))))+

	 o2.5*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))))+
	 o2.6*log(pi21.2*(exp(a200.2+a220)/(1+exp(a200.2+a220)+exp(a300.2+a320))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a200.2+a220+a202+a022)/(1+exp(a200.2+a220+a202+a022)+exp(a300.2+a320+a302+a022))))+

	 o2.7*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))))+
	 o2.8*log(pi12.2*(exp(a300.2+a302)/(1+exp(a200.2+a202)+exp(a300.2+a302))) +
	  (1-pi11.2-pi12.2-pi21.2)*(exp(a300.2+a320+a302+a022)/(1+exp(a200.2+a220+a202+a022)+exp(a300.2+a320+a302+a022))))

	),c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a200.2","a220","a202","a300.1","a300.2","a320","a302"),
	c("pi11.1","pi12.1","pi21.1","pi11.2","pi12.2","pi21.2",
	 "a200.1","a200.2","a220","a202","a300.1","a300.2","a320","a302",
	 "o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8",
	 "o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8","a022")
)

alfas<- seq(-5,5,0.01)
nalfas<-length(alfas)
z<-qnorm(0.975)

Qvs<-numeric(nalfas)
ws<-numeric(nalfas)
vars<-numeric(nalfas)
pvs<-numeric(nalfas)
meth<-numeric(nalfas)

require(MASS)

date()
for(i in 1:nalfas){ #Quase 15min. num Intel Core 2 Duo 1.7GHz
  alfa022<-alfas[i]
  mnar9.fit<-.nlmP(objfunc=mnar9.mlv,params=c(rep(1/4,6),rep(0,8)),lower=c(rep(0,6),rep(-Inf,8)),upper=c(rep(1,6),rep(Inf,8)),hessian=TRUE,print.level=0,iterlim=500)
  Qvs[i]<- -2*(-mnar9.fit$min-sat.lv)
  p<-mnar9.fit$est
  mnar9.analit<-mnar9.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],
  	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[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8],alfa022)
  fit<-funlinWLS(model=c("lin","log"),theta=as.vector(b+c(B%*%p[1:6])),Vtheta=B%*%ginv(attr(mnar9.analit,"hessian")[1,,])[1:6,1:6]%*%t(B),A1=A,X=rep(1,2))
  ws[i]<-fit$beta
  vars[i]<-c(fit$Vbeta)
  pvs[i]<-1-pchisq(fit$QwH,1)
}
date()

sel<-alfas>=-2 & alfas<=2
paste("Intervalo de ignorância: (",round(min(ws[sel]),2),";",round(max(ws[sel]),2),")",sep="")

winf<-(ws[sel]-z*sqrt(vars[sel]))
wsup<-(ws[sel]+z*sqrt(vars[sel]))
paste("Intervalo de incerteza: (",round(min(winf),2),";",round(max(wsup),2),")",sep="")

plot(alfas[sel],ws[sel],xlab=expression(paste(alpha[paste("022")])), ylab=expression(paste(omega,"'")),pch=16,type="n",ylim=c(min(winf),max(wsup)))
lines(alfas[sel],ws[sel],lwd=1)
lines(alfas[sel],winf,lwd=3)
lines(alfas[sel],wsup,lwd=3)
abline(0,0,lty=2)