#Comandos para reproduzir as análises do Exemplo 4, pp.116-126, 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")


##### 1º turno

#Satisfação com a vida x Candidato
SAT.TF<-c(14,217,132,44,21,256,140,18,3,64,41,19,3,55,35,13,1,66,40,8,7,2,0,1)
SAT.Zp<-cbind(rep(1,4)%x%diag(4),diag(4)%x%rep(1,4))
SAT.Rp<-c(4,4)
SAT.cd<-readCatdata(TF=SAT.TF,Zp=SAT.Zp,Rp=SAT.Rp)
SAT.mlcd<-satMarML(SAT.cd,epsilon1=1e-16,epsilon2=1e-16)
funlinWLS(model="lin",obj=SAT.mlcd,A1=cbind(diag(3),rep(0,3))%x%t(rep(1,4)),X=rbind(c(1,1,1),c(1,1,0),c(1,0,0)))
#3º parâmetro é \pi_S-\pi_M



#Renda familiar x Candidato
RF.TF<-c(24,55,82,117,89,31,18,28,79,127,108,64,4,14,35,34,22,11,7,13,13,38,18,14,4,4,18,34,32,17,16,13,7,4)
RF.Zp<-cbind(rep(1,4)%x%diag(6),diag(4)%x%rep(1,6))
RF.Rp<-c(6,4)
RF.cd<-readCatdata(TF=RF.TF,Zp=RF.Zp,Rp=RF.Rp)
RF.mlcd<-satMarML(RF.cd,epsilon1=1e-16,epsilon2=1e-16)
funlinWLS(model="lin",obj=RF.mlcd,A1=cbind(diag(3),rep(0,3))%x%t(rep(1,6)),X=rbind(c(1,1,1),c(1,1,0),c(1,0,0)))



#Grau de instrução x Candidato
INS.TF<-c(101,90,130,93,117,129,135,56,34,32,44,17,21,23,42,21,42,28,36,13)
INS.Zp<-cbind(rep(1,4)%x%diag(4))
INS.Rp<-4
INS.cd<-readCatdata(TF=INS.TF,Zp=INS.Zp,Rp=INS.Rp)
INS.mlcd<-satMarML(INS.cd,epsilon1=1e-16,epsilon2=1e-16)
funlinWLS(model="lin",obj=INS.mlcd,A1=cbind(diag(3),rep(0,3))%x%t(rep(1,4)),X=rbind(c(1,1,1),c(1,1,0),c(1,0,0)))



library(geoR)

#MAR utilizando Rejeição: distribuindo os não-respondentes aleatoriamente na intenção de voto dos candidatos em que não rejeitam
mar.mlv<-function(p){
  pi01<-p[1] #Anai Caproni-29-PCO
  pi02<-p[2] #Ciro-36-PTC
  pi03<-p[3] #Dirceu Travesso-16-PSTU
  pi04<-p[4] #Dra. Havanir-56-PRONA
  pi05<-p[5] #Francisco Rossi-31-PHS
  pi06<-p[6] #João Manuel-27-PSDC
  pi07<-p[7] #José Serra-45-PSDB
  pi08<-p[8] #Luiza Erundina-40-PSB
  pi09<-p[9] #Marta Suplicy-13-PT
  pi10<-p[10] #Osmar Lins-26-PAN
  pi11<-p[11] #Paulinho-12-PDT
  pi12<-p[12] #Paulo Maluf-11-PP
  pi13<-p[13] #Penna-43-PV
  #pi14=1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13 #Professor Walter Canoas-21-PCB

	value<- -(
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi12+pi13)+
#esses 25 não-respondentes rejeitam apenas o Maluf
    25*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi13)+
    3*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
#esses 7 não-respondentes rejeitam apenas a Marta
    7*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    2*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    2*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi09+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    3*log(pi01+pi02+pi03+pi04+pi05+pi06+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi08+pi09+pi10+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi06+pi07+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    4*log(pi01+pi02+pi03+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    2*log(pi01+pi02+pi03+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi13)+
    1*log(pi01+pi02+pi03+pi05+pi06+pi07+pi08+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi05+pi06+pi07+pi09+pi10+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi05+pi06+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi03+pi04+pi05+pi06+pi07+pi08+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi03+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi03+pi05+pi06+pi07+pi09)+
    1*log(pi01+pi03+pi06+pi09+pi10+pi11+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi04+pi05+pi06+pi07+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi04+pi05+pi06+pi08+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi04+pi06+pi08+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi02+pi03+pi04+pi05+pi07+pi08+pi09+pi10+pi11+pi12)+
    1*log(pi02+pi03+pi05+pi06+pi07+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi02+pi05+pi08+pi09+pi11)+
    1*log(pi03+pi08+pi09+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi05+pi07+pi08+pi09+pi11)+
    1*log(pi07+pi13)+ #esse não-respondente, rejeita todos os candidatos com exceção do Serra e do Penna
    1*log(pi12)+ #esse não-respondente, rejeita todos os candidatos com exceção do Maluf
#Respondentes
    0.001*log(pi01)+
    2*log(pi02)+
    1*log(pi03)+
    7*log(pi04)+
    19*log(pi05)+
    0.001*log(pi06)+
    414*log(pi07)+
    53*log(pi08)+
    437*log(pi09)+
    2*log(pi10)+
    14*log(pi11)+
    127*log(pi12)+
    7*log(pi13)+
    2*log(1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)
	)
	value
}

mar.der<-deriv3(~-(
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi12+pi13)+
    25*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi13)+
    3*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    7*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    2*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    2*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi07+pi09+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    3*log(pi01+pi02+pi03+pi04+pi05+pi06+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi08+pi09+pi10+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi05+pi06+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi04+pi06+pi07+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    4*log(pi01+pi02+pi03+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    2*log(pi01+pi02+pi03+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi13)+
    1*log(pi01+pi02+pi03+pi05+pi06+pi07+pi08+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi05+pi06+pi07+pi09+pi10+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi02+pi03+pi05+pi06+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi03+pi04+pi05+pi06+pi07+pi08+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi03+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi03+pi05+pi06+pi07+pi09)+
    1*log(pi01+pi03+pi06+pi09+pi10+pi11+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi04+pi05+pi06+pi07+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi04+pi05+pi06+pi08+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi01+pi04+pi06+pi08+pi09+pi10+pi11+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi02+pi03+pi04+pi05+pi06+pi07+pi08+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi02+pi03+pi04+pi05+pi07+pi08+pi09+pi10+pi11+pi12)+
    1*log(pi02+pi03+pi05+pi06+pi07+pi09+pi10+pi11+pi12+pi13+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi02+pi05+pi08+pi09+pi11)+
    1*log(pi03+pi08+pi09+1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)+
    1*log(pi05+pi07+pi08+pi09+pi11)+
    1*log(pi07+pi13)+
    1*log(pi12)+
    0.001*log(pi01)+
    2*log(pi02)+
    1*log(pi03)+
    7*log(pi04)+
    19*log(pi05)+
    0.001*log(pi06)+
    414*log(pi07)+
    53*log(pi08)+
    437*log(pi09)+
    2*log(pi10)+
    14*log(pi11)+
    127*log(pi12)+
    7*log(pi13)+
    2*log(1-pi01-pi02-pi03-pi04-pi05-pi06-pi07-pi08-pi09-pi10-pi11-pi12-pi13)

	),
	c("pi01","pi02","pi03","pi04","pi05","pi06","pi07","pi08","pi09","pi10","pi11","pi12","pi13"),
	c("pi01","pi02","pi03","pi04","pi05","pi06","pi07","pi08","pi09","pi10","pi11","pi12","pi13")
)

mar.fit<-.nlmP(objfunc=mar.mlv,params=c(rep(1/14,13)),lower=rep(0,13),upper=rep(1,13),hessian=T,print.level=0,iterlim=500)
p<-mar.fit$est
mar.analit<-mar.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])
round(rbind(
  mar.fit$est
,
  sqrt(diag(solve(attr(mar.analit,"hessian")[1,,])))
  ),4)

CC<-c(0,0,0,0,0,0,1,0,-1,0,0,0,0)
c(CC%*%p) #\pi_S-\pi_M
c(sqrt(CC%*%solve(attr(mar.analit,"hessian")[1,,])%*%CC))





#menos log-verossimilhança do modelo MNAR1, saturado
mnar1.mlv<-function(p,fr=TF){
  #freqüências e probabilidades j.i, indicando j|i, ou seja, coluna j dado a linha i ("m" para missing)
	o1.1<-fr[1,1];o2.1<-fr[1,2];o3.1<-fr[1,3];o4.1<-fr[1,4];om.1<-fr[1,5]
	o1.2<-fr[2,1];o2.2<-fr[2,2];o3.2<-fr[2,3];o4.2<-fr[2,4];om.2<-fr[2,5]
	o1.3<-fr[3,1];o2.3<-fr[3,2];o3.3<-fr[3,3];o4.3<-fr[3,4];om.3<-fr[3,5]
	o1.4<-fr[4,1];o2.4<-fr[4,2];o3.4<-fr[4,3];o4.4<-fr[4,4];om.4<-fr[4,5]

	pi1.1<-p[1]; pi2.1<-p[2]; pi3.1<-p[3]#; pi4.1<-1-pi1.1-pi2.1-pi3.1
	pi1.2<-p[4]; pi2.2<-p[5]; pi3.2<-p[6]#; pi4.2<-1-pi1.2-pi2.2-pi3.2
	pi1.3<-p[7]; pi2.3<-p[8]; pi3.3<-p[9]#; pi4.3<-1-pi1.3-pi2.3-pi3.3
	pi1.4<-p[10];pi2.4<-p[11];pi3.4<-p[12]#;pi4.4<-1-pi1.4-pi2.4-pi3.4

  a1<-p[13];a2<-p[14];a3<-p[15];a4<-p[16]

	value<- -(
	 o1.1*log(pi1.1*(1-a1))+o2.1*log(pi2.1*(1-a2))+o3.1*log(pi3.1*(1-a3))+o4.1*log((1-pi1.1-pi2.1-pi3.1)*(1-a4))+
	 o1.2*log(pi1.2*(1-a1))+o2.2*log(pi2.2*(1-a2))+o3.2*log(pi3.2*(1-a3))+o4.2*log((1-pi1.2-pi2.2-pi3.2)*(1-a4))+
	 o1.3*log(pi1.3*(1-a1))+o2.3*log(pi2.3*(1-a2))+o3.3*log(pi3.3*(1-a3))+o4.3*log((1-pi1.3-pi2.3-pi3.3)*(1-a4))+
	 o1.4*log(pi1.4*(1-a1))+o2.4*log(pi2.4*(1-a2))+o3.4*log(pi3.4*(1-a3))+o4.4*log((1-pi1.4-pi2.4-pi3.4)*(1-a4))+

	 om.1*log(pi1.1*a1+pi2.1*a2+pi3.1*a3+(1-pi1.1-pi2.1-pi3.1)*a4)+
	 om.2*log(pi1.2*a1+pi2.2*a2+pi3.2*a3+(1-pi1.2-pi2.2-pi3.2)*a4)+
	 om.3*log(pi1.3*a1+pi2.3*a2+pi3.3*a3+(1-pi1.3-pi2.3-pi3.3)*a4)+
	 om.4*log(pi1.4*a1+pi2.4*a2+pi3.4*a3+(1-pi1.4-pi2.4-pi3.4)*a4)

	)
	value
}

mnar1.der<-deriv3(~-(
	 o1.1*log(pi1.1*(1-a1))+o2.1*log(pi2.1*(1-a2))+o3.1*log(pi3.1*(1-a3))+o4.1*log((1-pi1.1-pi2.1-pi3.1)*(1-a4))+
	 o1.2*log(pi1.2*(1-a1))+o2.2*log(pi2.2*(1-a2))+o3.2*log(pi3.2*(1-a3))+o4.2*log((1-pi1.2-pi2.2-pi3.2)*(1-a4))+
	 o1.3*log(pi1.3*(1-a1))+o2.3*log(pi2.3*(1-a2))+o3.3*log(pi3.3*(1-a3))+o4.3*log((1-pi1.3-pi2.3-pi3.3)*(1-a4))+
	 o1.4*log(pi1.4*(1-a1))+o2.4*log(pi2.4*(1-a2))+o3.4*log(pi3.4*(1-a3))+o4.4*log((1-pi1.4-pi2.4-pi3.4)*(1-a4))+

	 om.1*log(pi1.1*a1+pi2.1*a2+pi3.1*a3+(1-pi1.1-pi2.1-pi3.1)*a4)+
	 om.2*log(pi1.2*a1+pi2.2*a2+pi3.2*a3+(1-pi1.2-pi2.2-pi3.2)*a4)+
	 om.3*log(pi1.3*a1+pi2.3*a2+pi3.3*a3+(1-pi1.3-pi2.3-pi3.3)*a4)+
	 om.4*log(pi1.4*a1+pi2.4*a2+pi3.4*a3+(1-pi1.4-pi2.4-pi3.4)*a4)

	),
	c("pi1.1","pi2.1","pi3.1","pi1.2","pi2.2","pi3.2","pi1.3","pi2.3","pi3.3","pi1.4","pi2.4","pi3.4",
    "a1","a2","a3","a4"),
	c("pi1.1","pi2.1","pi3.1","pi1.2","pi2.2","pi3.2","pi1.3","pi2.3","pi3.3","pi1.4","pi2.4","pi3.4",
    "a1","a2","a3","a4",
 	  "o1.1","o2.1","o3.1","o4.1","om.1","o1.2","o2.2","o3.2","o4.2","om.2",
    "o1.3","o2.3","o3.3","o4.3","om.3","o1.4","o2.4","o3.4","o4.4","om.4")
)

library(geoR)



#Grau de instrução x Candidato
TF<-cbind(c(101, 90,130,93),
          c(117,129,135,56),
          c( 34, 32, 44,17),
          c( 21, 23, 42,21),
          c( 42, 28, 36,13))
sat.lv<-sum(TF*log(TF/rowSums(TF)))
mnar1.fit<-.nlmP(objfunc=mnar1.mlv,params=c(rep(1/4,12),rep(1/2,4)),lower=rep(0,16),upper=rep(1,16),hessian=T,print.level=0,iterlim=500)
-2*(-mnar1.fit$min-sat.lv)
p<-mnar1.fit$est
mnar1.analit<-mnar1.der(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8],p[9],p[10],p[11],p[12],p[13],p[14],p[15],p[16],
	TF[1,1],TF[1,2],TF[1,3],TF[1,4],TF[1,5],TF[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],
  TF[3,1],TF[3,2],TF[3,3],TF[3,4],TF[3,5],TF[4,1],TF[4,2],TF[4,3],TF[4,4],TF[4,5])
round(matrix(mnar1.fit$est[1:12],4,3,byrow=T),4) #Estim.das probs.condicionais de intenção de votar nos candidatos dado o grau de instr.
rowSums(TF)/sum(TF) #Estim.das probs.marginais de se ter o grau de instr.
x<-c((rowSums(TF)/sum(TF))%*%matrix(mnar1.fit$est[1:12],4,3,byrow=T))
x #Estim.das probs.marginais de intenção de votar nos candidatos
x[1]-x[2] #\pi_S-\pi_M
round(mnar1.fit$est[13:16],4) #aloca quase todos os não-respondentes para o Maluf



#Satisfação com a vida (desconsiderando omissos) x Candidato
TF<-cbind(c(14,217,132,44),
          c(21,256,140,18),
          c( 3, 64, 41,19),
          c( 3, 55, 35,13),
          c( 1, 66, 40, 8))
sat.lv<-sum(TF*log(TF/rowSums(TF)))
mnar1.fit<-.nlmP(objfunc=mnar1.mlv,params=c(rep(1/4,12),rep(1/2,4)),lower=rep(0,16),upper=rep(1,16),hessian=T,print.level=0,iterlim=500)
-2*(-mnar1.fit$min-sat.lv)
p<-mnar1.fit$est
mnar1.analit<-mnar1.der(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8],p[9],p[10],p[11],p[12],p[13],p[14],p[15],p[16],
	TF[1,1],TF[1,2],TF[1,3],TF[1,4],TF[1,5],TF[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],
  TF[3,1],TF[3,2],TF[3,3],TF[3,4],TF[3,5],TF[4,1],TF[4,2],TF[4,3],TF[4,4],TF[4,5])
round(matrix(mnar1.fit$est[1:12],4,3,byrow=T),4) #Estim.das probs.condicionais de intenção de votar nos candidatos dado o grau de satisf.
rowSums(TF)/sum(TF) #Estim.das probs.marginais de se ter o grau de satisf.
x<-c((rowSums(TF)/sum(TF))%*%matrix(mnar1.fit$est[1:12],4,3,byrow=T))
x #Estim.das probs.marginais de intenção de votar nos candidatos
x[1]-x[2] #\pi_S-\pi_M
round(mnar1.fit$est[13:16],4) #aloca quase todos os não-respondentes para Outros candidatos





TF<-c(414,437,234,119)

sat.lv<-sum(TF*log(TF/sum(TF)))

#menos log-verossimilhança do modelo MNAR, sobresaturado
mnar.mlv<-function(p,fr=TF,a1=s1,a2=s2){
	o1<-fr[1];o2<-fr[2];o3<-fr[3];om<-fr[4]
	pi1<-p[1];pi2<-p[2]#;pi3<-1-pi1-pi2
  a3<-p[3]
	value<- -(
	 o1*log(pi1*(1-a1))+o2*log(pi2*(1-a2))+o3*log((1-pi1-pi2)*(1-a3))+
	 om*log(pi1*a1+pi2*a2+(1-pi1-pi2)*a3)
	)
	value
}

mnar.der<-deriv3(~-(
	 o1*log(pi1*(1-a1))+o2*log(pi2*(1-a2))+o3*log((1-pi1-pi2)*(1-a3))+
	 om*log(pi1*a1+pi2*a2+(1-pi1-pi2)*a3)

	),
	c("pi1","pi2","a3"),
	c("pi1","pi2","a3",
	 "a1","a2",
	 "o1","o2","o3","om")
)



CC<-c(1,-1,0)
ss1<-ss2<-seq(0.05,0.20,0.01)
l1<-length(ss1)
l2<-length(ss2)
res<-matrix(0,l1*l2,11)
i<-0
for(i1 in 1:l1) {
 for(i2 in 1:l2) {
  i<-i+1
  res[i,1]<-s1<-ss1[i1]
  res[i,2]<-s2<-ss2[i2]
  mnar.fit<-.nlmP(objfunc=mnar.mlv,params=c(rep(1/3,2),0.1),lower=rep(0,3),upper=rep(1,3),hessian=T,print.level=0,iterlim=500)
  res[i,3]<- -2*(-mnar.fit$min-sat.lv)
  res[i,4:6]<-mnar.fit$est
  p<-mnar.fit$est
  mnar.final<-mnar.der(p[1],p[2],p[3],s1,s2,TF[1],TF[2],TF[3],TF[4])
  res[i,7:9]<-diag(solve(attr(mnar.final,"hessian")[1,,]))
  res[i,10]<-c(CC%*%mnar.fit$est)
  res[i,11]<-c(CC%*%solve(attr(mnar.final,"hessian")[1,,])%*%CC)
 }
}
round( rbind(apply(res[,3:11],2,min),apply(res[,3:11],2,max)) ,4)
c(min(res[,10]),max(res[,10])) #Intervalo de ignorância, p.122

x<-cbind(res[,10]-1.96*sqrt(res[,11]) , res[,10]+1.96*sqrt(res[,11]) )
round( rbind(apply(x,2,min),apply(x,2,max)) ,4)
c(min(x[,1]),max(x[,2])) #Intervalo de incerteza, p.122

library(lattice)
trellis.par.set(background=list(col="#FFFFFF"))
wireframe(c(res[,10],x)~rep(res[,1],3)*rep(res[,2],3),groups=c(1,2,3)%x%rep(1,l1*l2),scales=list(arrows=FALSE),screen=list(z=40,x=-90,y=15),xlab=expression(paste(alpha["S"])),ylab=expression(paste(alpha["M"])),zlab=expression(paste(pi["S"]-pi["M"])),col.groups=c(0,0,0),zoom=0.75 )
#Figura 3.7

#Últimas linhas da p.122
res[ res[,1]==0.16 & res[,2]==0.05 , ]
x[ res[,1]==0.16 & res[,2]==0.05 , ]



CC<-c(1,-1,0)
var<-seq(-5,5,0.5)
ss1<-ss2<-exp(var)/(1+exp(var))
l1<-length(ss1)
l2<-length(ss2)
res<-matrix(0,l1*l2,11)
i<-0
for(i1 in 1:l1) {
 for(i2 in 1:l2) {
  i<-i+1
  res[i,1]<-s1<-ss1[i1]
  res[i,2]<-s2<-ss2[i2]
  mnar.fit<-.nlmP(objfunc=mnar.mlv,params=c(rep(1/3,2),0.1),lower=rep(0,3),upper=rep(1,3),hessian=T,print.level=0,iterlim=500)
  res[i,3]<- -2*(-mnar.fit$min-sat.lv)
  res[i,4:6]<-mnar.fit$est
  p<-mnar.fit$est
  mnar.final<-mnar.der(p[1],p[2],p[3],s1,s2,TF[1],TF[2],TF[3],TF[4])
  res[i,7:9]<-diag(solve(attr(mnar.final,"hessian")[1,,]))
  res[i,10]<-c(CC%*%mnar.fit$est)
  res[i,11]<-c(CC%*%solve(attr(mnar.final,"hessian")[1,,])%*%CC)
 }
}
round( rbind(apply(res[,3:11],2,min),apply(res[,3:11],2,max)) ,4)
c(min(res[,10]),max(res[,10])) #Intervalo de ignorância, p.123

x<-cbind(res[,10]-1.96*sqrt(res[,11]) , res[,10]+1.96*sqrt(res[,11]) )
round( rbind(apply(x,2,min),apply(x,2,max)) ,4)
c(min(x[,1]),max(x[,2])) #Intervalo de incerteza, p.123





TF<-c(414,437,234,119)
sat.lv<-sum(TF*log(TF/sum(TF)))

#menos log-verossimilhança do modelo MNAR, saturado, apenas os alphas são estimados (post-mortem)
mnarP.mlv<-function(p,fr=TF,pi1=2686396/6167371,pi2=2209264/6167371){
	o1<-fr[1];o2<-fr[2];o3<-fr[3];om<-fr[4]
	#pi1<-p[1];pi2<-p[2];pi3<-p[3]#;pi4<-1-pi1-pi2-pi3
  a1<-p[1]
  a2<-p[2]
  a3<-p[3]
	value<- -(
	 o1*log(pi1*(1-a1))+o2*log(pi2*(1-a2))+o3*log((1-pi1-pi2)*(1-a3))+
	 om*log(pi1*a1+pi2*a2+(1-pi1-pi2)*a3)
	)
	value
}
mnarP.der<-deriv3(~-(
	 o1*log(pi1*(1-a1))+o2*log(pi2*(1-a2))+o3*log((1-pi1-pi2)*(1-a3))+
	 om*log(pi1*a1+pi2*a2+(1-pi1-pi2)*a3)
	),
	c("a1","a2","a3"),
	c("a1","a2","a3",
	 "pi1","pi2",
	 "o1","o2","o3","om")
)

mnarP.fit<-.nlmP(objfunc=mnarP.mlv,params=rep(0.1,3),lower=rep(0,3),upper=rep(1,3),hessian=TRUE,print.level=0,iterlim=500)
-2*(-mnarP.fit$min-sat.lv) #Estatística de razão de verossimilhanças para o ajuste do modelo, p.123
p<-mnarP.fit$est
mnarP.analit<-mnarP.der(p[1],p[2],p[3],2686396/6167371,2209264/6167371,TF[1],TF[2],TF[3],TF[4])
round(mnarP.fit$est,4) #EMV de \alpha_S, \alpha_M e \alpha_O, p.123
round(sqrt(diag(solve(attr(mnarP.analit,"hessian")[1,,]))),4) #Erros padrões de \alpha_S, \alpha_M e \alpha_O, p.123
CC<-rbind(c(1,-1,0),c(1,0,-1))
#Valor-P do teste de Wald de H:\alpha_S=\alpha_M=\alpha_O
1-pchisq(c( t(CC%*%mnarP.fit$est)%*%solve(CC%*%solve(attr(mnarP.analit,"hessian")[1,,])%*%t(CC))%*%(CC%*%mnarP.fit$est) ),2)





##### 2º turno

TF<-c(980,821,201)
sat.lv<-sum(TF*log(TF/sum(TF)))

#menos log-verossimilhança do modelo MNAR, sobresaturado
mnar.mlv<-function(p,fr=TF,a1=s1){
	o1<-fr[1];o2<-fr[2];om<-fr[3]
	pi1<-p[1]#;pi2<-1-pi1
  a2<-p[2]
	value<- -(
	 o1*log(pi1*(1-a1))+o2*log((1-pi1)*(1-a2))+
	 om*log(pi1*a1+(1-pi1)*a2)
	)
	value
}

mnar.der<-deriv3(~-(
	 o1*log(pi1*(1-a1))+o2*log((1-pi1)*(1-a2))+
	 om*log(pi1*a1+(1-pi1)*a2)
	),
	c("pi1","a2"),
	c("pi1","a2",
	 "a1",
	 "o1","o2","om")
)

ss1<-seq(0.001,0.999,0.001)
l1<-length(ss1)
res<-matrix(0,l1,6)
i<-0
date()
for(i1 in 1:l1) {
 i<-i+1
 res[i,1]<-s1<-ss1[i1]
 mnar.fit<-.nlmP(objfunc=mnar.mlv,params=c(0.5,0.1),lower=rep(0,2),upper=rep(1,2),hessian=T,print.level=0,iterlim=500)
 res[i,2]<- -2*(-mnar.fit$min-sat.lv)
 res[i,3:4]<-mnar.fit$est
 p<-mnar.fit$est
 mnar.analit<-mnar.der(p[1],p[2],s1,TF[1],TF[2],TF[3])
 res[i,5:6]<-diag(solve(attr(mnar.analit,"hessian")[1,,]))
}
date()

round( rbind(apply(res[,2:6],2,min),apply(res[,2:6],2,max)) ,3)
c(min(res[,3]),max(res[,3])) #Intervalo de ignorância, p.125

x<-cbind(res[,3]-1.96*sqrt(res[,5]) , res[,3]+1.96*sqrt(res[,5]) )
round( rbind(apply(x,2,min),apply(x,2,max)) ,3)
c(min(x[,1]),max(x[,2])) #Intervalo de incerteza, p.125

plot(res[,rep(1,3)],cbind(res[,3],x),xlab=expression(paste(alpha["S"])),ylab=expression(paste(pi["S"])),pch=16,type="n")
lines(res[,1],res[,3],lwd=1)
lines(res[,1],x[,1],lwd=3)
lines(res[,1],x[,2],lwd=3)
abline(h=0.5,lty=2)
#Figura 3.8

#\alpha_S de 5% a 20%
abline(v=0.05,lty=2)
abline(v=0.20,lty=2)

max(res[ss1<=0.17,2])
plot(res[,1],res[,2],xlab=expression(paste(alpha["S"])),ylab="TRV",pch=16,type="l")

min(x[ss1>=0.065,1])





ss1<-seq(0.05,0.20,0.001)
l1<-length(ss1)
res<-matrix(0,l1,6)
i<-0
date()
for(i1 in 1:l1) {
 i<-i+1
 res[i,1]<-s1<-ss1[i1]
 mnar.fit<-.nlmP(objfunc=mnar.mlv,params=c(0.5,0.1),lower=rep(0,2),upper=rep(1,2),hessian=T,print.level=0,iterlim=500)
 res[i,2]<- -2*(-mnar.fit$min-sat.lv)
 res[i,3:4]<-mnar.fit$est
 p<-mnar.fit$est
 mnar.analit<-mnar.der(p[1],p[2],s1,TF[1],TF[2],TF[3])
 res[i,5:6]<-diag(solve(attr(mnar.analit,"hessian")[1,,]))
}
date()

round( rbind(apply(res[,2:6],2,min),apply(res[,2:6],2,max)) ,3)
c(min(res[,3]),max(res[,3])) #Intervalo de ignorância, p.125

x<-cbind(res[,3]-1.96*sqrt(res[,5]) , res[,3]+1.96*sqrt(res[,5]) )
round( rbind(apply(x,2,min),apply(x,2,max)) ,3)
c(min(x[,1]),max(x[,2])) #Intervalo de incerteza, p.125





#menos log-verossimilhança do modelo MNAR, saturado, apenas os alphas são estimados (post-mortem)
mnarP.mlv<-function(p,fr=TF,pi1=3330179/6070331){
	o1<-fr[1];o2<-fr[2];om<-fr[3]
  a1<-p[1]
  a2<-p[2]
	value<- -(
	 o1*log(pi1*(1-a1))+o2*log((1-pi1)*(1-a2))+
	 om*log(pi1*a1+(1-pi1)*a2)
	)
	value
}
mnarP.der<-deriv3(~-(
	 o1*log(pi1*(1-a1))+o2*log((1-pi1)*(1-a2))+
	 om*log(pi1*a1+(1-pi1)*a2)
	),
	c("a1","a2"),
	c("a1","a2",
	 "pi1",
	 "o1","o2","om")
)

mnarP.fit<-.nlmP(objfunc=mnarP.mlv,params=rep(0.1,2),lower=rep(0,2),upper=rep(1,2),hessian=T,print.level=0,iterlim=500)
-2*(-mnarP.fit$min-sat.lv) #Estatística de razão de verossimilhanças para o ajuste do modelo, p.126
p<-mnarP.fit$est
mnarP.analit<-mnarP.der(p[1],p[2],3330179/6070331,TF[1],TF[2],TF[3])
round(mnarP.fit$est,4) #EMV de \alpha_S e \alpha_M
round(sqrt(diag(solve(attr(mnarP.analit,"hessian")[1,,]))),4) #Erros padrões de \alpha_S e \alpha_M
CC<-c(1,-1)
round(c(CC%*%mnarP.fit$est),4)
round(sqrt(c(CC%*%solve(attr(mnarP.analit,"hessian")[1,,])%*%CC)),4)
round(c(CC%*%mnarP.fit$est)-1.96*sqrt(c(CC%*%solve(attr(mnarP.analit,"hessian")[1,,])%*%CC)),4)
round(c(CC%*%mnarP.fit$est)+1.96*sqrt(c(CC%*%solve(attr(mnarP.analit,"hessian")[1,,])%*%CC)),4)
c(CC%*%mnarP.fit$est^2)/c(CC%*%solve(attr(mnarP.analit,"hessian")[1,,])%*%CC)
1-pchisq(c(CC%*%solve(attr(mnarP.analit,"hessian")[1,,])%*%CC),1) #Valor-P do teste de Wald de H:\alpha_S=\alpha_M