#Comandos para reproduzir as análises do Exemplo 3, pp.108-115, 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<-c(7,11,2,3,9,5,1e-5,10,4, 8,7,3,0, 0,7,14,7)
Zp<-cbind(rbind( cbind(rep(1,2)%x%diag(3),rep(0,6)),
                 cbind(matrix(0,3,3),rep(1,3)) ),
          rbind( cbind(rep(1,3),matrix(0,3,3)),
                 cbind(rep(0,6),rep(1,2)%x%diag(3)) ) )
Rp<-c(4,4)

kappamodel<-c("add","exp","lin","log","lin","exp","lin","log","lin")
kA1<-rbind(
 c(rep(c(1,0,0,0),2),1),
 rep(1,9),
 kronecker(diag(3),t(rep(1,3))),
 kronecker(t(rep(1,3)),diag(3)) )
kA2<-rbind(
 cbind(diag(2),matrix(0,2,6)),
 cbind(matrix(0,3,2),kronecker(t(rep(1,2)),diag(3))) )
kA3<-cbind( c(1,0),c(1,1),-c(2,1)%*%t(rep(1,3)) )
kA4<-t(c(1,-1))
W1<-c(1,0.75,0,0.75,1,0.75,0,0.75,1) #pesos quadráticos Fleiss e Cohen (1973)
W2<-c(1,0.5,0,0.5,1,0.5,0,0.5,1) #pesos absolutos Agresti (2002)
kw1A1<-rbind(
 t(W1),
 rep(1,9),
 kronecker(diag(3),t(rep(1,3))),
 kronecker(t(rep(1,3)),diag(3)) )
kw2A1<-rbind(
 t(W2),
 rep(1,9),
 kronecker(diag(3),t(rep(1,3))),
 kronecker(t(rep(1,3)),diag(3)) )
kwA2<-rbind(
 cbind(diag(2),matrix(0,2,6)),
 cbind(matrix(0,9,2),  cbind(kronecker(diag(3),rep(1,3)) , kronecker(rep(1,3),diag(3)))  ) )
kw1A3<-cbind( c(1,0),c(1,1),kronecker(-c(2,1),t(W1)) )
kw2A3<-cbind( c(1,0),c(1,1),kronecker(-c(2,1),t(W2)) )
Adif<-rbind(c(0,1,1,-1,0,0,-1,0,0),
            c(0,-1,0,1,0,1,0,-1,0) )





#Intervalos para o melhor-pior caso apresentados na Tabela 3.11 (veja as Tabelas B.4 e B.5 no Apêndice B.2)

#Alocação A
(38-25)/97 #Limite superior para \pi_{1+}-\pi_{+1}
(17-51)/97 #Limite inferior para \pi_{2+}-\pi_{+2}
#Alocação B
(20-25)/97 #Limite inferior para \pi_{1+}-\pi_{+1}
(63-51)/97 #Limite superior para \pi_{2+}-\pi_{+2}
#Alocação C, Limites superiores para \kappa, \kappa_{w1} e \kappa_{w2}
alocC<-readCatdata(TF=c(15,11,2, 10,30,8, 0,10,11))
funlinWLS(model=kappamodel,obj=alocC,A1=kA1  ,A2=kA2 ,A3=kA3  ,A4=kA4,PI1=-1,X=1)$beta
funlinWLS(model=kappamodel,obj=alocC,A1=kw1A1,A2=kwA2,A3=kw1A3,A4=kA4,PI1=-1,X=1)$beta
funlinWLS(model=kappamodel,obj=alocC,A1=kw2A1,A2=kwA2,A3=kw2A3,A4=kA4,PI1=-1,X=1)$beta
#Alocação D, Limites inferiores para \kappa, \kappa_{w1} e \kappa_{w2}
alocD<-readCatdata(TF=c(7,18,5,11,9,12,7,24,4))
funlinWLS(model=kappamodel,obj=alocD,A1=kA1  ,A2=kA2 ,A3=kA3  ,A4=kA4,PI1=-1,X=1)$beta
funlinWLS(model=kappamodel,obj=alocD,A1=kw1A1,A2=kwA2,A3=kw1A3,A4=kA4,PI1=-1,X=1)$beta
funlinWLS(model=kappamodel,obj=alocD,A1=kw2A1,A2=kwA2,A3=kw2A3,A4=kA4,PI1=-1,X=1)$beta




catdata.acc<-readCatdata(TF=TF[1:9])
dif.acc<-funlinWLS(model="lin",obj=catdata.acc,A1=Adif,X=diag(2))
dif.acc #Tabela 3.15, ACC, estimativas e erros padrões para \pi_{i+}-\pi_{+i}, i=1,2
waldTest(dif.acc,diag(2)) #Tabela 3.15, ACC, teste de Wald de H:\pi_{i+}-\pi_{+i}, i=1,2
#... ou funlinWLS(model="lin",obj=catdata.acc,U=Adif[,-9])

#Tabela 3.15, ACC, estimativas e erros padrões para \kappa, \kappa_{w1} e \kappa_{w2}, e testes de Wald de os parâmetros serem iguais a zero
funlinWLS(model=kappamodel,obj=catdata.acc,A1=kA1  ,A2=kA2 ,A3=kA3  ,A4=kA4,PI1=-1,X=1)
funlinWLS(model=kappamodel,obj=catdata.acc,A1=kw1A1,A2=kwA2,A3=kw1A3,A4=kA4,PI1=-1,X=1)
funlinWLS(model=kappamodel,obj=catdata.acc,A1=kw2A1,A2=kwA2,A3=kw2A3,A4=kA4,PI1=-1,X=1)





catdata<-readCatdata(TF=TF,Zp=Zp,Rp=Rp)
satmarml<-satMarML(catdata)
satmarml #TRV do ajuste do mecanismo MCAR (Nota de rodapé 7)
satmarml$yst #Tabela 3.14, MAR, t=1,2,3
round(satmarml$yst$st1.1+satmarml$yst$st1.2+satmarml$yst$st1.3,1) #Tabela 3.14, MAR, Total

dif.acc<-funlinWLS(model="lin",obj=satmarml,A1=Adif,X=diag(2))
dif.acc #Tabela 3.15, MAR, estimativas e erros padrões para \pi_{i+}-\pi_{+i}, i=1,2
waldTest(dif.acc,diag(2)) #Tabela 3.15, MAR, teste de Wald de H:\pi_{i+}-\pi_{+i}, i=1,2

#Tabela 3.15, MAR, estimativas e erros padrões para \kappa, \kappa_{w1} e \kappa_{w2}, e testes de Wald de os parâmetros serem iguais a zero
funlinWLS(model=kappamodel,obj=satmarml,A1=kA1  ,A2=kA2 ,A3=kA3  ,A4=kA4,PI1=-1,X=1)
funlinWLS(model=kappamodel,obj=satmarml,A1=kw1A1,A2=kwA2,A3=kw1A3,A4=kA4,PI1=-1,X=1)
funlinWLS(model=kappamodel,obj=satmarml,A1=kw2A1,A2=kwA2,A3=kw2A3,A4=kA4,PI1=-1,X=1)





require(geoR)
TF<-TF[-(13:14)] #Daqui para frente vamos eliminar os zeros estruturais, que tinham sido usados apenas como um artifício (vide pp.58-60)
sat.lv<-sum(TF*log(TF/sum(TF))) #log-verossimilhança sob modelo saturado

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





#menos log-verossimilhança do modelo MAR_{red}, 4 g.l.
marred.mlv<-function(p,fr=TF){
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]

	pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
	pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
	pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32

	a2<-p[9]
	a3<-p[10]

	value<- -(
	 o1*log(pi11*(1-a2))+
	 o2*log(pi12*(1-a2))+
	 o3*log(pi13*(1-a2))+
	 o4*log(pi21*(1-a2-a3))+
	 o5*log(pi22*(1-a2-a3))+
	 o6*log(pi23*(1-a2-a3))+
	 o7*log(pi31*(1-a3))+
	 o8*log(pi32*(1-a3))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a3))+

	 o10*log(pi11*a2 + pi21*a2)+
	 o11*log(pi12*a2 + pi22*a2)+
	 o12*log(pi13*a2 + pi23*a2)+

	 o13*log(pi21*a3 + pi31*a3)+
	 o14*log(pi22*a3 + pi32*a3)+
	 o15*log(pi23*a3 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a3)

	)
	value
}

marred.der<-deriv3(~-(
	 o1*log(pi11*(1-a2))+
	 o2*log(pi12*(1-a2))+
	 o3*log(pi13*(1-a2))+
	 o4*log(pi21*(1-a2-a3))+
	 o5*log(pi22*(1-a2-a3))+
	 o6*log(pi23*(1-a2-a3))+
	 o7*log(pi31*(1-a3))+
	 o8*log(pi32*(1-a3))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a3))+

	 o10*log(pi11*a2 + pi21*a2)+
	 o11*log(pi12*a2 + pi22*a2)+
	 o12*log(pi13*a2 + pi23*a2)+

	 o13*log(pi21*a3 + pi31*a3)+
	 o14*log(pi22*a3 + pi32*a3)+
	 o15*log(pi23*a3 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a3)

	),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a2","a3"),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a2","a3",
	 "o1","o2","o3","o4","o5","o6","o7","o8","o9","o10",
	 "o11","o12","o13","o14","o15")
)

marred.fit<-.nlmP(objfunc=marred.mlv,params=c(rep(1/9,8),rep(1/3,2)),lower=rep(0,10),upper=rep(1,10),hessian=TRUE,print.level=0,iterlim=500)
-2*(-marred.fit$min-sat.lv) #Estatística de razão de verossimilhanças do ajuste do modelo, Tabela 3.13, MAR_{red}
1-pchisq(-2*(-marred.fit$min-sat.lv),4) #Valor-P do teste do ajuste do modelo, fim da p.109, MAR_{red}

p<-marred.fit$est
fr<-rep(sum(TF),15)
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]
	pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
	pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
	pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32
	a2<-p[9]
	a3<-p[10]
x<-rbind(matrix(c(
	 o1*c(pi11*(1-a2)),
	 o2*c(pi12*(1-a2)),
	 o3*c(pi13*(1-a2)),
	 o4*c(pi21*(1-a2-a3)),
	 o5*c(pi22*(1-a2-a3)),
	 o6*c(pi23*(1-a2-a3)),
	 o7*c(pi31*(1-a3)),
	 o8*c(pi32*(1-a3)),
	 o9*c((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a3))
),3,3,byrow=TRUE),matrix(c(
   o10*c(pi11*a2 , pi21*a2),
	 o11*c(pi12*a2 , pi22*a2),
	 o12*c(pi13*a2 , pi23*a2)
),2,3),matrix(c(
   o13*c(pi21*a3 , pi31*a3),
	 o14*c(pi22*a3 , pi32*a3),
	 o15*c(pi23*a3 , (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a3)
),2,3) )
round(x,1) #Tabela 3.14, MAR_{red}, t=1,2,3
round(rowSums(rbind(t(x[4:5,]),t(x[6:7,]))),1)
round(rbind(x[1,]+x[4,],x[2,]+x[5,]+x[6,],x[3,]+x[7,]),1) #Tabela 3.14, MAR_{red}, Total
round(rowSums(rbind(x[1,]+x[4,],x[2,]+x[5,]+x[6,],x[3,]+x[7,])),1)
round(colSums(x),1)
sum(x)


marred.analit<-marred.der(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8],p[9],p[10],
	TF[1],TF[2],TF[3],TF[4],TF[5],TF[6],TF[7],TF[8],TF[9],TF[10],TF[11],TF[12],TF[13],TF[14],TF[15])

dif.marred<-funlinWLS(theta=b+c(B%*%marred.fit$est[1:8]),Vtheta=B%*%solve(attr(marred.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model="lin",A1=Adif,X=diag(2))
dif.marred #Tabela 3.15, MAR_{red}, estimativas e erros padrões para \pi_{i+}-\pi_{+i}, i=1,2
waldTest(dif.marred,diag(2)) #Tabela 3.15, MAR_{red}, teste de Wald de H:\pi_{i+}-\pi_{+i}, i=1,2

#Tabela 3.15, MAR_{red}, estimativas e erros padrões para \kappa, \kappa_{w1} e \kappa_{w2}, e testes de Wald de os parâmetros serem iguais a zero
funlinWLS(theta=b+c(B%*%marred.fit$est[1:8]),Vtheta=B%*%solve(attr(marred.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kA1  ,A2=kA2 ,A3=kA3  ,A4=kA4,PI1=-1,X=1)
funlinWLS(theta=b+c(B%*%marred.fit$est[1:8]),Vtheta=B%*%solve(attr(marred.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kw1A1,A2=kwA2,A3=kw1A3,A4=kA4,PI1=-1,X=1)
funlinWLS(theta=b+c(B%*%marred.fit$est[1:8]),Vtheta=B%*%solve(attr(marred.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kw2A1,A2=kwA2,A3=kw2A3,A4=kA4,PI1=-1,X=1)





#menos log-verossimilhança do modelo MNAR1, 2 g.l.
mnar1.mlv<-function(p,fr=TF){
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]

	pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
	pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
	pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32

	a2.1<-p[9]; a2.2<-p[10]
	a3.1<-p[11];a3.2<-p[12]

	value<- -(
	 o1*log(pi11*(1-a2.1))+
	 o2*log(pi12*(1-a2.2))+
	 o3*log(pi13*(1-a2.2))+
	 o4*log(pi21*(1-a2.2-a3.1))+
	 o5*log(pi22*(1-a2.1-a3.1))+
	 o6*log(pi23*(1-a2.1-a3.2))+
	 o7*log(pi31*(1-a3.2))+
	 o8*log(pi32*(1-a3.2))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a3.1))+

	 o10*log(pi11*a2.1 + pi21*a2.2)+
	 o11*log(pi12*a2.2 + pi22*a2.1)+
	 o12*log(pi13*a2.2 + pi23*a2.1)+

	 o13*log(pi21*a3.1 + pi31*a3.2)+
	 o14*log(pi22*a3.1 + pi32*a3.2)+
	 o15*log(pi23*a3.2 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a3.1)

	)
	value
}

mnar1.der<-deriv3(~-(
	 o1*log(pi11*(1-a2.1))+
	 o2*log(pi12*(1-a2.2))+
	 o3*log(pi13*(1-a2.2))+
	 o4*log(pi21*(1-a2.2-a3.1))+
	 o5*log(pi22*(1-a2.1-a3.1))+
	 o6*log(pi23*(1-a2.1-a3.2))+
	 o7*log(pi31*(1-a3.2))+
	 o8*log(pi32*(1-a3.2))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a3.1))+

	 o10*log(pi11*a2.1 + pi21*a2.2)+
	 o11*log(pi12*a2.2 + pi22*a2.1)+
	 o12*log(pi13*a2.2 + pi23*a2.1)+

	 o13*log(pi21*a3.1 + pi31*a3.2)+
	 o14*log(pi22*a3.1 + pi32*a3.2)+
	 o15*log(pi23*a3.2 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a3.1)

	),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a2.1","a2.2","a3.1","a3.2"),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a2.1","a2.2","a3.1","a3.2",
	 "o1","o2","o3","o4","o5","o6","o7","o8","o9","o10",
	 "o11","o12","o13","o14","o15")
)

mnar1.fit<-.nlmP(objfunc=mnar1.mlv,params=c(rep(1/9,8),rep(1/3,4)),lower=rep(0,12),upper=rep(1,12),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.13, MNAR1
1-pchisq(-2*(-mnar1.fit$min-sat.lv),2) #Valor-P do teste do ajuste do modelo, fim da p.109, MNAR1

p<-mnar1.fit$est
fr<-rep(sum(TF),15)
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]
	pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
	pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
	pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32
	a2.1<-p[9]; a2.2<-p[10]
	a3.1<-p[11];a3.2<-p[12]
x<-rbind(matrix(c(
	 o1*c(pi11*(1-a2.1)),
	 o2*c(pi12*(1-a2.2)),
	 o3*c(pi13*(1-a2.2)),
	 o4*c(pi21*(1-a2.2-a3.1)),
	 o5*c(pi22*(1-a2.1-a3.1)),
	 o6*c(pi23*(1-a2.1-a3.2)),
	 o7*c(pi31*(1-a3.2)),
	 o8*c(pi32*(1-a3.2)),
	 o9*c((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a3.1))
),3,3,byrow=TRUE),matrix(c(
	 o10*c(pi11*a2.1 , pi21*a2.2),
	 o11*c(pi12*a2.2 , pi22*a2.1),
	 o12*c(pi13*a2.2 , pi23*a2.1)
),2,3),matrix(c(
	 o13*c(pi21*a3.1 , pi31*a3.2),
	 o14*c(pi22*a3.1 , pi32*a3.2),
	 o15*c(pi23*a3.2 , (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a3.1)
),2,3) )
round(x,1) #Tabela 3.14, MNAR1, t=1,2,3
round(rowSums(rbind(t(x[4:5,]),t(x[6:7,]))),1)
round(rbind(x[1,]+x[4,],x[2,]+x[5,]+x[6,],x[3,]+x[7,]),1) #Tabela 3.14, MNAR1, Total
round(rowSums(rbind(x[1,]+x[4,],x[2,]+x[5,]+x[6,],x[3,]+x[7,])),1)
round(colSums(x),1)
sum(x)


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],
	TF[1],TF[2],TF[3],TF[4],TF[5],TF[6],TF[7],TF[8],TF[9],TF[10],TF[11],TF[12],TF[13],TF[14],TF[15])

dif.mnar1<-funlinWLS(theta=b+c(B%*%mnar1.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar1.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model="lin",A1=Adif,X=diag(2))
dif.mnar1 #Tabela 3.15, MNAR1, estimativas e erros padrões para \pi_{i+}-\pi_{+i}, i=1,2
waldTest(dif.mnar1,diag(2)) #Tabela 3.15, MNAR1, teste de Wald de H:\pi_{i+}-\pi_{+i}, i=1,2

#Tabela 3.15, MNAR1, estimativas e erros padrões para \kappa, \kappa_{w1} e \kappa_{w2}, e testes de Wald de os parâmetros serem iguais a zero
funlinWLS(theta=b+c(B%*%mnar1.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar1.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kA1  ,A2=kA2 ,A3=kA3  ,A4=kA4,PI1=-1,X=1)
funlinWLS(theta=b+c(B%*%mnar1.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar1.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kw1A1,A2=kwA2,A3=kw1A3,A4=kA4,PI1=-1,X=1)
funlinWLS(theta=b+c(B%*%mnar1.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar1.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kw2A1,A2=kwA2,A3=kw2A3,A4=kA4,PI1=-1,X=1)





#menos log-verossimilhança do modelo MNAR2, 2 g.l.
mnar2.mlv<-function(p,fr=TF){
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]

	pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
	pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
	pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32

	a2.1<-p[9]; a2.2<-p[10]
	a3.1<-p[11];a3.2<-p[12]

	value<- -(
	 o1*log(pi11*(1-a2.1))+
	 o2*log(pi12*(1-a2.1))+
	 o3*log(pi13*(1-a2.1))+
	 o4*log(pi21*(1-a2.2-a3.1))+
	 o5*log(pi22*(1-a2.2-a3.1))+
	 o6*log(pi23*(1-a2.2-a3.1))+
	 o7*log(pi31*(1-a3.2))+
	 o8*log(pi32*(1-a3.2))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a3.2))+

	 o10*log(pi11*a2.1 + pi21*a2.2)+
	 o11*log(pi12*a2.1 + pi22*a2.2)+
	 o12*log(pi13*a2.1 + pi23*a2.2)+

	 o13*log(pi21*a3.1 + pi31*a3.2)+
	 o14*log(pi22*a3.1 + pi32*a3.2)+
	 o15*log(pi23*a3.1 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a3.2)

	)
	value
}

mnar2.der<-deriv3(~-(
	 o1*log(pi11*(1-a2.1))+
	 o2*log(pi12*(1-a2.1))+
	 o3*log(pi13*(1-a2.1))+
	 o4*log(pi21*(1-a2.2-a3.1))+
	 o5*log(pi22*(1-a2.2-a3.1))+
	 o6*log(pi23*(1-a2.2-a3.1))+
	 o7*log(pi31*(1-a3.2))+
	 o8*log(pi32*(1-a3.2))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a3.2))+

	 o10*log(pi11*a2.1 + pi21*a2.2)+
	 o11*log(pi12*a2.1 + pi22*a2.2)+
	 o12*log(pi13*a2.1 + pi23*a2.2)+

	 o13*log(pi21*a3.1 + pi31*a3.2)+
	 o14*log(pi22*a3.1 + pi32*a3.2)+
	 o15*log(pi23*a3.1 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a3.2)

	),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a2.1","a2.2","a3.1","a3.2"),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a2.1","a2.2","a3.1","a3.2",
	 "o1","o2","o3","o4","o5","o6","o7","o8","o9","o10",
	 "o11","o12","o13","o14","o15")
)

mnar2.fit<-.nlmP(objfunc=mnar2.mlv,params=c(rep(1/9,8),rep(1/3,4)),lower=rep(0,12),upper=rep(1,12),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.13, MNAR2
1-pchisq(-2*(-mnar2.fit$min-sat.lv),2) #Valor-P do teste do ajuste do modelo, fim da p.109, MNAR2

p<-mnar2.fit$est
fr<-rep(sum(TF),15)
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]
	pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
	pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
	pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32
	a2.1<-p[9]; a2.2<-p[10]
	a3.1<-p[11];a3.2<-p[12]
x<-rbind(matrix(c(
	 o1*(pi11*(1-a2.1)),
	 o2*(pi12*(1-a2.1)),
	 o3*(pi13*(1-a2.1)),
	 o4*(pi21*(1-a2.2-a3.1)),
	 o5*(pi22*(1-a2.2-a3.1)),
	 o6*(pi23*(1-a2.2-a3.1)),
	 o7*(pi31*(1-a3.2)),
	 o8*(pi32*(1-a3.2)),
	 o9*((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a3.2))
),3,3,byrow=TRUE),matrix(c(
	 o10*c(pi11*a2.1 , pi21*a2.2),
	 o11*c(pi12*a2.1 , pi22*a2.2),
	 o12*c(pi13*a2.1 , pi23*a2.2)
),2,3),matrix(c(
	 o13*c(pi21*a3.1 , pi31*a3.2),
	 o14*c(pi22*a3.1 , pi32*a3.2),
	 o15*c(pi23*a3.1 , (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a3.2)
),2,3) )
round(x,1) #Tabela 3.14, MNAR2, t=1,2,3
round(rowSums(rbind(t(x[4:5,]),t(x[6:7,]))),1)
round(rbind(x[1,]+x[4,],x[2,]+x[5,]+x[6,],x[3,]+x[7,]),1) #Tabela 3.14, MNAR2, Total
round(rowSums(rbind(x[1,]+x[4,],x[2,]+x[5,]+x[6,],x[3,]+x[7,])),1)
round(colSums(x),1)
sum(x)


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],
	TF[1],TF[2],TF[3],TF[4],TF[5],TF[6],TF[7],TF[8],TF[9],TF[10],TF[11],TF[12],TF[13],TF[14],TF[15])

dif.mnar2<-funlinWLS(theta=b+c(B%*%mnar2.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar2.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model="lin",A1=Adif,X=diag(2))
dif.mnar2 #Tabela 3.15, MNAR2, estimativas e erros padrões para \pi_{i+}-\pi_{+i}, i=1,2
waldTest(dif.mnar2,diag(2)) #Tabela 3.15, MNAR2, teste de Wald de H:\pi_{i+}-\pi_{+i}, i=1,2

#Tabela 3.15, MNAR2, estimativas e erros padrões para \kappa, \kappa_{w1} e \kappa_{w2}, e testes de Wald de os parâmetros serem iguais a zero
funlinWLS(theta=b+c(B%*%mnar2.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar2.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kA1  ,A2=kA2 ,A3=kA3  ,A4=kA4,PI1=-1,X=1)
funlinWLS(theta=b+c(B%*%mnar2.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar2.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kw1A1,A2=kwA2,A3=kw1A3,A4=kA4,PI1=-1,X=1)
funlinWLS(theta=b+c(B%*%mnar2.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar2.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kw2A1,A2=kwA2,A3=kw2A3,A4=kA4,PI1=-1,X=1)





#menos log-verossimilhança do modelo MNAR3, saturado
mnar3.mlv<-function(p,fr=TF){
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]

	pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
	pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
	pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32

	a1<-p[9]; a2<-p[10];a3<-p[11]
	a4<-p[12];a5<-p[13];a6<-p[14]

	value<- -(
	 o1*log(pi11*(1-a1))+
	 o2*log(pi12*(1-a3))+
	 o3*log(pi13*(1-a5))+
	 o4*log(pi21*(1-a2-a1))+
	 o5*log(pi22*(1-a4-a3))+
	 o6*log(pi23*(1-a6-a5))+
	 o7*log(pi31*(1-a2))+
	 o8*log(pi32*(1-a4))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a6))+

	 o10*log(pi11*a1 + pi21*a2)+
	 o11*log(pi12*a3 + pi22*a4)+
	 o12*log(pi13*a5 + pi23*a6)+

	 o13*log(pi21*a1 + pi31*a2)+
	 o14*log(pi22*a3 + pi32*a4)+
	 o15*log(pi23*a5 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a6)

	)
	value
}

mnar3.der<-deriv3(~-(
	 o1*log(pi11*(1-a1))+
	 o2*log(pi12*(1-a3))+
	 o3*log(pi13*(1-a5))+
	 o4*log(pi21*(1-a2-a1))+
	 o5*log(pi22*(1-a4-a3))+
	 o6*log(pi23*(1-a6-a5))+
	 o7*log(pi31*(1-a2))+
	 o8*log(pi32*(1-a4))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a6))+

	 o10*log(pi11*a1 + pi21*a2)+
	 o11*log(pi12*a3 + pi22*a4)+
	 o12*log(pi13*a5 + pi23*a6)+

	 o13*log(pi21*a1 + pi31*a2)+
	 o14*log(pi22*a3 + pi32*a4)+
	 o15*log(pi23*a5 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a6)

	),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a1","a2","a3","a4","a5","a6"),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a1","a2","a3","a4","a5","a6",
	 "o1","o2","o3","o4","o5","o6","o7","o8","o9","o10",
	 "o11","o12","o13","o14","o15")
)

mnar3.fit<-.nlmP(objfunc=mnar3.mlv,params=c(rep(1/9,8),rep(1/3,6)),lower=rep(0,14),upper=rep(1,14),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.13, MNAR3

p<-mnar3.fit$est
fr<-rep(sum(TF),15)
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]
  pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
  pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
  pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32
  a1<-p[9]; a2<-p[10];a3<-p[11]
  a4<-p[12];a5<-p[13];a6<-p[14]
x<-rbind(matrix(c(
	 o1*c(pi11*(1-a1)),
	 o2*c(pi12*(1-a3)),
	 o3*c(pi13*(1-a5)),
	 o4*c(pi21*(1-a2-a1)),
	 o5*c(pi22*(1-a4-a3)),
	 o6*c(pi23*(1-a6-a5)),
	 o7*c(pi31*(1-a2)),
	 o8*c(pi32*(1-a4)),
	 o9*c((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a6))
),3,3,byrow=TRUE),matrix(c(
	 o10*c(pi11*a1 , pi21*a2),
	 o11*c(pi12*a3 , pi22*a4),
	 o12*c(pi13*a5 , pi23*a6)
),2,3),matrix(c(
	 o13*c(pi21*a1 , pi31*a2),
	 o14*c(pi22*a3 , pi32*a4),
	 o15*c(pi23*a5 , (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a6)
),2,3) )
round(x,1) #Tabela 3.14, MNAR3, t=1,2,3
round(rowSums(rbind(t(x[4:5,]),t(x[6:7,]))),1)
round(rbind(x[1,]+x[4,],x[2,]+x[5,]+x[6,],x[3,]+x[7,]),1) #Tabela 3.14, MNAR3, Total
round(rowSums(rbind(x[1,]+x[4,],x[2,]+x[5,]+x[6,],x[3,]+x[7,])),1)
round(colSums(x),1)
sum(x)


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],
	TF[1],TF[2],TF[3],TF[4],TF[5],TF[6],TF[7],TF[8],TF[9],TF[10],TF[11],TF[12],TF[13],TF[14],TF[15])

dif.mnar3<-funlinWLS(theta=b+c(B%*%mnar3.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar3.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model="lin",A1=Adif,X=diag(2))
dif.mnar3 #Tabela 3.15, MNAR3, estimativas e erros padrões para \pi_{i+}-\pi_{+i}, i=1,2
waldTest(dif.mnar3,diag(2)) #Tabela 3.15, MNAR3, teste de Wald de H:\pi_{i+}-\pi_{+i}, i=1,2

#Tabela 3.15, MNAR3, estimativas e erros padrões para \kappa, \kappa_{w1} e \kappa_{w2}, e testes de Wald de os parâmetros serem iguais a zero
funlinWLS(theta=b+c(B%*%mnar3.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar3.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kA1  ,A2=kA2 ,A3=kA3  ,A4=kA4,PI1=-1,X=1)
funlinWLS(theta=b+c(B%*%mnar3.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar3.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kw1A1,A2=kwA2,A3=kw1A3,A4=kA4,PI1=-1,X=1)
funlinWLS(theta=b+c(B%*%mnar3.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar3.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kw2A1,A2=kwA2,A3=kw2A3,A4=kA4,PI1=-1,X=1)





#menos log-verossimilhança da estimação do modelo MNAR4, saturado
mnar4.mlv<-function(p,fr=TF){
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]

	pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
	pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
	pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32

	a1<-p[9]; a2<-p[10];a3<-p[11]
	a4<-p[12];a5<-p[13];a6<-p[14]

	value<- -(
	 o1*log(pi11*(1-a1))+
	 o2*log(pi12*(1-a3))+
	 o3*log(pi13*(1-a5))+
	 o4*log(pi21*(1-a2-a6))+
	 o5*log(pi22*(1-a4-a4))+
	 o6*log(pi23*(1-a6-a2))+
	 o7*log(pi31*(1-a5))+
	 o8*log(pi32*(1-a3))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a1))+

	 o10*log(pi11*a1 + pi21*a2)+
	 o11*log(pi12*a3 + pi22*a4)+
	 o12*log(pi13*a5 + pi23*a6)+

	 o13*log(pi21*a6 + pi31*a5)+
	 o14*log(pi22*a4 + pi32*a3)+
	 o15*log(pi23*a2 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a1)

	)
	value
}

mnar4.der<-deriv3(~-(
	 o1*log(pi11*(1-a1))+
	 o2*log(pi12*(1-a3))+
	 o3*log(pi13*(1-a5))+
	 o4*log(pi21*(1-a2-a6))+
	 o5*log(pi22*(1-a4-a4))+
	 o6*log(pi23*(1-a6-a2))+
	 o7*log(pi31*(1-a5))+
	 o8*log(pi32*(1-a3))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a1))+

	 o10*log(pi11*a1 + pi21*a2)+
	 o11*log(pi12*a3 + pi22*a4)+
	 o12*log(pi13*a5 + pi23*a6)+

	 o13*log(pi21*a6 + pi31*a5)+
	 o14*log(pi22*a4 + pi32*a3)+
	 o15*log(pi23*a2 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a1)

	),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a1","a2","a3","a4","a5","a6"),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a1","a2","a3","a4","a5","a6",
	 "o1","o2","o3","o4","o5","o6","o7","o8","o9","o10",
	 "o11","o12","o13","o14","o15")
)

mnar4.fit<-.nlmP(objfunc=mnar4.mlv,params=c(rep(1/9,8),rep(1/3,6)),lower=rep(0,14),upper=rep(1,14),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.13, MNAR4

p<-mnar4.fit$est
fr<-rep(sum(TF),15)
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]
  pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
  pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
  pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32
  a1<-p[9]; a2<-p[10];a3<-p[11]
  a4<-p[12];a5<-p[13];a6<-p[14]
x<-rbind(matrix(c(
	 o1*c(pi11*(1-a1)),
	 o2*c(pi12*(1-a3)),
	 o3*c(pi13*(1-a5)),
	 o4*c(pi21*(1-a2-a6)),
	 o5*c(pi22*(1-a4-a4)),
	 o6*c(pi23*(1-a6-a2)),
	 o7*c(pi31*(1-a5)),
	 o8*c(pi32*(1-a3)),
	 o9*c((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a1))
),3,3,byrow=TRUE),matrix(c(
	 o10*c(pi11*a1 , pi21*a2),
	 o11*c(pi12*a3 , pi22*a4),
	 o12*c(pi13*a5 , pi23*a6)
),2,3),matrix(c(
	 o13*c(pi21*a6 , pi31*a5),
	 o14*c(pi22*a4 , pi32*a3),
	 o15*c(pi23*a2 , (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a1)
),2,3) )
round(x,1) #Tabela 3.14, MNAR4, t=1,2,3
round(rowSums(rbind(t(x[4:5,]),t(x[6:7,]))),1)
round(rbind(x[1,]+x[4,],x[2,]+x[5,]+x[6,],x[3,]+x[7,]),1) #Tabela 3.14, MNAR4, Total
round(rowSums(rbind(x[1,]+x[4,],x[2,]+x[5,]+x[6,],x[3,]+x[7,])),1)
round(colSums(x),1)
sum(x)


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],TF[2],TF[3],TF[4],TF[5],TF[6],TF[7],TF[8],TF[9],TF[10],TF[11],TF[12],TF[13],TF[14],TF[15])

dif.mnar4<-funlinWLS(theta=b+c(B%*%mnar4.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar4.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model="lin",A1=Adif,X=diag(2))
dif.mnar4 #Tabela 3.15, MNAR4, estimativas e erros padrões para \pi_{i+}-\pi_{+i}, i=1,2
waldTest(dif.mnar4,diag(2)) #Tabela 3.15, MNAR4, teste de Wald de H:\pi_{i+}-\pi_{+i}, i=1,2

#Tabela 3.15, MNAR4, estimativas e erros padrões para \kappa, \kappa_{w1} e \kappa_{w2}, e testes de Wald de os parâmetros serem iguais a zero
funlinWLS(theta=b+c(B%*%mnar4.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar4.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kA1  ,A2=kA2 ,A3=kA3  ,A4=kA4,PI1=-1,X=1)
funlinWLS(theta=b+c(B%*%mnar4.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar4.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kw1A1,A2=kwA2,A3=kw1A3,A4=kA4,PI1=-1,X=1)
funlinWLS(theta=b+c(B%*%mnar4.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar4.analit,"hessian")[1,,])[1:8,1:8]%*%t(B),
  model=kappamodel,A1=kw2A1,A2=kwA2,A3=kw2A3,A4=kA4,PI1=-1,X=1)





#menos log-verossimilhança do modelo MNAR5, SOBRESATURADO
mnar5.mlv<-function(p,a7=s7,a8=s8,fr=TF){
	o1<-fr[1];o2<-fr[2];o3<-fr[3];o4<-fr[4];o5<-fr[5];o6<-fr[6];o7<-fr[7];o8<-fr[8];o9<-fr[9]
	o10<-fr[10];o11<-fr[11];o12<-fr[12];o13<-fr[13];o14<-fr[14];o15<-fr[15]

	pi11<-p[1]; pi12<-p[2]; pi13<-p[3]
	pi21<-p[4]; pi22<-p[5]; pi23<-p[6]
	pi31<-p[7]; pi32<-p[8]#; pi33=1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32

	a1<-p[9];a2<-p[10];a3<-p[11];a4<-p[12];a5<-p[13];a6<-p[14]

	value<- -(
	 o1*log(pi11*(1-a7))+
	 o2*log(pi12*(1-a8))+
	 o3*log(pi13*(1-a8))+
	 o4*log(pi21*(1-a1-a4))+
	 o5*log(pi22*(1-a2-a5))+
	 o6*log(pi23*(1-a3-a6))+
	 o7*log(pi31*(1-a8))+
	 o8*log(pi32*(1-a8))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a7))+

	 o10*log(pi11*a7 + pi21*a1)+
	 o11*log(pi12*a8 + pi22*a2)+
	 o12*log(pi13*a8 + pi23*a3)+

	 o13*log(pi21*a4 + pi31*a8)+
	 o14*log(pi22*a5 + pi32*a8)+
	 o15*log(pi23*a6 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a7)

	)
	value
}

mnar5.der<-deriv3(~-(
	 o1*log(pi11*(1-a7))+
	 o2*log(pi12*(1-a8))+
	 o3*log(pi13*(1-a8))+
	 o4*log(pi21*(1-a1-a4))+
	 o5*log(pi22*(1-a2-a5))+
	 o6*log(pi23*(1-a3-a6))+
	 o7*log(pi31*(1-a8))+
	 o8*log(pi32*(1-a8))+
	 o9*log((1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*(1-a7))+

	 o10*log(pi11*a7 + pi21*a1)+
	 o11*log(pi12*a8 + pi22*a2)+
	 o12*log(pi13*a8 + pi23*a3)+

	 o13*log(pi21*a4 + pi31*a8)+
	 o14*log(pi22*a5 + pi32*a8)+
	 o15*log(pi23*a6 + (1-pi11-pi12-pi13-pi21-pi22-pi23-pi31-pi32)*a7)

	),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a1","a2","a3","a4","a5","a6"),
	c("pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32",
	 "a1","a2","a3","a4","a5","a6","a7","a8",
	 "o1","o2","o3","o4","o5","o6","o7","o8","o9","o10",
	 "o11","o12","o13","o14","o15")
)

b1<-seq(-0.6,0.4,0.01)
b2<-seq(-0.6,0.4,0.01)
lb1<-length(b1)
lb2<-length(b2)
bv<-expand.grid(b1,b2)
ff<-function(b)
	c((media-b)%*%inf%*%(media-b))
plot(b1,b2,type="n",xlab=expression(paste(pi[paste(1,"+")]-pi[paste("+",1)])),ylab=expression(paste(pi[paste(2,"+")]-pi[paste("+",2)])),pch=16,xlim=c(-0.25,0.35),ylim=c(-0.6,0.4))
abline(h=0,lty=2)
abline(v=0,lty=2)

var<-seq(-5,5,0.5)
ss7<-ss8<-exp(var)/(1+exp(var))
ls7<-length(ss7)
ls8<-length(ss8)
res<-matrix(0,ls7*ls8,14)

i<-0
date() #Cerca de 5min. num Intel Core 2 Duo 1.7GHz
for(i8 in 1:ls8) {
 cat(i8,"\n")
 for(i7 in 1:ls7) {
  i<-i+1
  cat(" ",i7)

  s7<-ss7[i7]
  s8<-ss8[i8]

  mnar5.fit<-.nlmP(objfunc=mnar5.mlv,params=c(rep(1/9,8),rep(1/3,6)),lower=rep(0,14),upper=rep(1,14),hessian=TRUE,print.level=0,iterlim=500)
  p<-mnar5.fit$est
  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],s7,s8,
	TF[1],TF[2],TF[3],TF[4],TF[5],TF[6],TF[7],TF[8],TF[9],TF[10],TF[11],TF[12],TF[13],TF[14],TF[15])

  dif.mnar5<-funlinWLS(theta=b+c(B%*%mnar5.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar5.analit,"hessian")[1,,],tol=1e-323)[1:8,1:8]%*%t(B),
    model="lin",A1=Adif,X=diag(2))
  mnar5.kappa  <-funlinWLS(theta=b+c(B%*%mnar5.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar5.analit,"hessian")[1,,],tol=1e-323)[1:8,1:8]%*%t(B),
    model=kappamodel,A1=kA1  ,A2=kA2 ,A3=kA3  ,A4=kA4,PI1=-1,X=1)
  mnar5.kappaw1<-funlinWLS(theta=b+c(B%*%mnar5.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar5.analit,"hessian")[1,,],tol=1e-323)[1:8,1:8]%*%t(B),
    model=kappamodel,A1=kw1A1,A2=kwA2,A3=kw1A3,A4=kA4,PI1=-1,X=1)
  mnar5.kappaw2<-funlinWLS(theta=b+c(B%*%mnar5.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar5.analit,"hessian")[1,,],tol=1e-323)[1:8,1:8]%*%t(B),
    model=kappamodel,A1=kw2A1,A2=kwA2,A3=kw2A3,A4=kA4,PI1=-1,X=1)

  res[i,]<-c(s7,s8,-2*(-mnar5.fit$min-sat.lv),
	mnar5.kappa$beta  ,mnar5.kappa$beta  -qnorm(0.975)*c(sqrt(mnar5.kappa$Vbeta  )),mnar5.kappa$beta  +qnorm(0.975)*c(sqrt(mnar5.kappa$Vbeta  )),
	mnar5.kappaw1$beta,mnar5.kappaw1$beta-qnorm(0.975)*c(sqrt(mnar5.kappaw1$Vbeta)),mnar5.kappaw1$beta+qnorm(0.975)*c(sqrt(mnar5.kappaw1$Vbeta)),
	mnar5.kappaw2$beta,mnar5.kappaw2$beta-qnorm(0.975)*c(sqrt(mnar5.kappaw2$Vbeta)),mnar5.kappaw2$beta+qnorm(0.975)*c(sqrt(mnar5.kappaw2$Vbeta)),
        dif.mnar5$beta
	)

  media<-dif.mnar5$beta
  inf<-solve(dif.mnar5$Vbeta)
  qstat<-matrix(apply(bv,1,ff),lb1,lb2)
  contour(b1,b2,qstat,levels=qchisq(0.95,2),add=TRUE,drawlabels=FALSE)
 
 }
 cat("\n")
}
date()
# O gráfico final é a Figura 3.5 (a)

#Intervalos de ignorância e de 95% de incerteza para \kappa, \kappa_{w1} e \kappa_{w2}, e de ignorância para \pi_{i+}-\pi_{+i}, i=1,2
round( rbind(apply(res[,3:14],2,min),apply(res[,3:14],2,max)) ,3)

res[res[,3]>100,1:2] #mostra os alfas que tiveram o pior ajuste

kappa      <-matrix(res[,4] ,ls7,ls8)
kappalinf  <-matrix(res[,5] ,ls7,ls8)
kappalsup  <-matrix(res[,6] ,ls7,ls8)
kappaw1    <-matrix(res[,7] ,ls7,ls8)
kappaw1linf<-matrix(res[,8] ,ls7,ls8)
kappaw1lsup<-matrix(res[,9] ,ls7,ls8)
kappaw2    <-matrix(res[,10],ls7,ls8)
kappaw2linf<-matrix(res[,11],ls7,ls8)
kappaw2lsup<-matrix(res[,12],ls7,ls8)

library(lattice)
trellis.par.set(background=list(col="#FFFFFF")) #trellis.par.get("background")

#Figura 3.6, \kappa
wireframe(c(res[,4:6])~rep(res[,1],3)*rep(res[,2],3),groups=c(1,2,3)%x%rep(1,ls7*ls8),scales=list(arrows=FALSE),screen=list(z=30,x=-75),xlab=expression(paste(alpha[7])),ylab=expression(paste(alpha[8])),zlab=expression(paste(kappa)),col.groups=c(0,0,0) )

#Figura 3.6, \kappa_{w1}
wireframe(c(res[,7:9])~rep(res[,1],3)*rep(res[,2],3),groups=c(1,2,3)%x%rep(1,ls7*ls8),scales=list(arrows=FALSE),screen=list(z=30,x=-75),xlab=expression(paste(alpha[7])),ylab=expression(paste(alpha[8])),zlab=expression(paste(kappa["w1"])),col.groups=c(0,0,0) )

#Figura 3.6, \kappa_{w2}
wireframe(c(res[,10:12])~rep(res[,1],3)*rep(res[,2],3),groups=c(1,2,3)%x%rep(1,ls7*ls8),scales=list(arrows=FALSE),screen=list(z=30,x=-75),xlab=expression(paste(alpha[7])),ylab=expression(paste(alpha[8])),zlab=expression(paste(kappa["w2"])),col.groups=c(0,0,0) )



b1<-seq(-0.6,0.4,0.01)
b2<-seq(-0.6,0.4,0.01)
lb1<-length(b1)
lb2<-length(b2)
bv<-expand.grid(b1,b2)
ff<-function(b)
	c((media-b)%*%inf%*%(media-b))
plot(b1,b2,type="n",xlab=expression(paste(pi[paste(1,"+")]-pi[paste("+",1)])),ylab=expression(paste(pi[paste(2,"+")]-pi[paste("+",2)])),pch=16,xlim=c(-0.25,0.35),ylim=c(-0.6,0.4))
abline(h=0,lty=2)
abline(v=0,lty=2)

s7<-0.01
s8<-0.01
mnar5.fit<-.nlmP(objfunc=mnar5.mlv,params=c(rep(1/9,8),rep(1/3,6)),lower=rep(0,14),upper=rep(1,14),hessian=TRUE,print.level=0,iterlim=500)
p<-mnar5.fit$est
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],s7,s8,
  TF[1],TF[2],TF[3],TF[4],TF[5],TF[6],TF[7],TF[8],TF[9],TF[10],TF[11],TF[12],TF[13],TF[14],TF[15])
dif.mnar5<-funlinWLS(theta=b+c(B%*%mnar5.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar5.analit,"hessian")[1,,],tol=1e-323)[1:8,1:8]%*%t(B),
  model="lin",A1=Adif,X=diag(2))
media<-dif.mnar5$beta
inf<-solve(dif.mnar5$Vbeta)
qstat<-matrix(apply(bv,1,ff),lb1,lb2)
contour(b1,b2,qstat,levels=qchisq(0.95,2),add=TRUE,drawlabels=FALSE)
text(-0.15,0.33,"(0.01,0.01)")

s7<-0.99
s8<-0.01
mnar5.fit<-.nlmP(objfunc=mnar5.mlv,params=c(rep(1/9,8),rep(1/3,6)),lower=rep(0,14),upper=rep(1,14),hessian=TRUE,print.level=0,iterlim=500)
p<-mnar5.fit$est
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],s7,s8,
  TF[1],TF[2],TF[3],TF[4],TF[5],TF[6],TF[7],TF[8],TF[9],TF[10],TF[11],TF[12],TF[13],TF[14],TF[15])
dif.mnar5<-funlinWLS(theta=b+c(B%*%mnar5.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar5.analit,"hessian")[1,,],tol=1e-323)[1:8,1:8]%*%t(B),
  model="lin",A1=Adif,X=diag(2))
media<-dif.mnar5$beta
inf<-solve(dif.mnar5$Vbeta)
qstat<-matrix(apply(bv,1,ff),lb1,lb2)
contour(b1,b2,qstat,levels=qchisq(0.95,2),add=TRUE,drawlabels=FALSE)
text(0.12,-0.06,"(0.99,0.01)")

s7<-0.01
s8<-0.99
mnar5.fit<-.nlmP(objfunc=mnar5.mlv,params=c(rep(1/9,8),rep(1/3,6)),lower=rep(0,14),upper=rep(1,14),hessian=TRUE,print.level=0,iterlim=500)
p<-mnar5.fit$est
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],s7,s8,
  TF[1],TF[2],TF[3],TF[4],TF[5],TF[6],TF[7],TF[8],TF[9],TF[10],TF[11],TF[12],TF[13],TF[14],TF[15])
dif.mnar5<-funlinWLS(theta=b+c(B%*%mnar5.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar5.analit,"hessian")[1,,],tol=1e-323)[1:8,1:8]%*%t(B),
  model="lin",A1=Adif,X=diag(2))
media<-dif.mnar5$beta
inf<-solve(dif.mnar5$Vbeta)
qstat<-matrix(apply(bv,1,ff),lb1,lb2)
contour(b1,b2,qstat,levels=qchisq(0.95,2),add=TRUE,drawlabels=FALSE)
text(-0.11,-0.11,"(0.01,0.99)")

s7<-0.99
s8<-0.99
mnar5.fit<-.nlmP(objfunc=mnar5.mlv,params=c(rep(1/9,8),rep(1/3,6)),lower=rep(0,14),upper=rep(1,14),hessian=TRUE,print.level=0,iterlim=500)
p<-mnar5.fit$est
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],s7,s8,
  TF[1],TF[2],TF[3],TF[4],TF[5],TF[6],TF[7],TF[8],TF[9],TF[10],TF[11],TF[12],TF[13],TF[14],TF[15])
dif.mnar5<-funlinWLS(theta=b+c(B%*%mnar5.fit$est[1:8]),Vtheta=B%*%solve(attr(mnar5.analit,"hessian")[1,,],tol=1e-323)[1:8,1:8]%*%t(B),
  model="lin",A1=Adif,X=diag(2))
media<-dif.mnar5$beta
inf<-solve(dif.mnar5$Vbeta)
qstat<-matrix(apply(bv,1,ff),lb1,lb2)
contour(b1,b2,qstat,levels=qchisq(0.95,2),add=TRUE,drawlabels=FALSE)
text(0.23,-0.51,"(0.99,0.99)")
# O gráfico final é a Figura 3.5 (b)