#Comandos para reproduzir as análises do Exemplo 1, pp.44-48, 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.

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

#Obs.: Ao organizar os comandos, decidi substituir os algoritmos de otimização usados na época do mestrado, nlm e .nlmP (biblioteca geoR),
#pelo optim e constrOptim, mais atuais. Também reescrevi as verossimilhanças de maneira mais intuitiva e menos massante, mas que não possibilitaram
#facilmente a obtenção das derivadas analíticas, deriv3, para a matriz hessiana, o que me fez utilizar a matriz hessiana numérica.
#Com essas mudanças, os resultados desses comandos, em alguns casos, não coincidem exatamente com os apresentados na dissertação,
#embora, em geral, as diferenças sejam pequenas.


TF<-c(4512,21009,3394,24132,1049,1135,142,464,1224)
padroes<-cbind(diag(4),diag(2)%x%rep(1,2),rep(1,2)%x%diag(2),rep(1,4))
classes<-c(4,2,2,1)

#Intervalo para o melhor-pior caso (Tabela 1.21)
exp(funlinWLS(model=c("lin","log"),obj=readCatdata(TF=c(4512,22522,5895,24132)),A1=c(1,-1,-1,1),XL=1)$beta)
exp(funlinWLS(model=c("lin","log"),obj=readCatdata(TF=c(6927,21009,3394,25731)),A1=c(1,-1,-1,1),XL=1)$beta)

catdata.acc<-readCatdata(TF=TF[1:4])
loglin.acc<-funlinWLS(model=c("lin","log"),obj=catdata.acc,A1=c(1,-1,-1,1),XL=1)
exp(loglin.acc$beta) #Estimativa da razão de chances, ACC, Tabela 1.22
exp(loglin.acc$beta+c(-1,1)*qnorm(0.975)*sqrt(loglin.acc$Vbeta)) #IC(95%,razão de chances), ACC, Tabela 1.22

catdata<-readCatdata(TF=TF[1:8],Zp=padroes[,5:8],Rp=classes[2:3])
mcar.ml<-satMarML(catdata,missing="MCAR")
mar.ml<-satMarML(catdata)
mar.ml #Testes de ajuste do mecanismo MCAR

mcar.ml$yst #Tabela 1.23, MCAR, W=1,2,3
mar.ml$yst #Tabela 1.23, MAR, W=1,2,3
TF[9]*mar.ml$theta #Tabela 1.23, MCAR/MAR, W=4

loglin.mcar<-funlinWLS(model=c("lin","log"),obj=mcar.ml,A1=c(1,-1,-1,1),XL=1)
loglin.mar<-funlinWLS(model=c("lin","log"),obj=mar.ml,A1=c(1,-1,-1,1),XL=1)
exp(loglin.mcar$beta) #Estimativa da razão de chances, MCAR, Tabela 1.22
exp(loglin.mcar$beta+c(-1,1)*qnorm(0.975)*sqrt(loglin.mcar$Vbeta)) #IC(95%,razão de chances), MCAR, Tabela 1.22
exp(loglin.mar$beta) #Estimativa da razão de chances, MAR, Tabela 1.22
exp(loglin.mar$beta+c(-1,1)*qnorm(0.975)*sqrt(loglin.mar$Vbeta)) #IC(95%,razão de chances), MAR, Tabela 1.22

gammatij<-function(pij,lambdat.ij) {
  totij<-length(pij)+1
  tottij<-length(lambdat.ij)+totij
  if(tottij%%totij!=0) {
    stop("comprimento de pij nao compativel com o de lambdat.ij")
  }
  gtij<-numeric(tottij)
  cont<-0
  for(tt in 1:(tottij/totij)) {
    for(ij in 1:totij) {
      cont<-cont+1
      gtij[cont]<-ifelse(ij==totij,1-sum(pij[-totij]),pij[ij])*ifelse(tt==tottij/totij,1-sum(lambdat.ij[seq(ij,(tottij/totij-2)*totij+ij,totij)]),lambdat.ij[cont])
    }
  }
  list(totij=totij,tottij=tottij,gtij=gtij)
}
mveross<-function(p,TF,padroes,classes,estrutura) {
  gtij<-gammatij(pij=p[1:3],lambdat.ij=estrutura(p[-(1:3)]))
  totp<-ncol(padroes)
  if(totp!=sum(classes)){
    stop("classes e padroes nao sao compativeis")
  }
  veross<-0
  for(ind in 1:(gtij$tottij/gtij$totij)) {
    faixa<-(1+ifelse(ind==1,0,sum(classes[1:(ind-1)]))):sum(classes[1:ind])
    veross<-veross+sum(TF[faixa]*log(diag(t(padroes[,faixa])%*%diag(gtij$gtij[(1+(ind-1)*gtij$totij):(ind*gtij$totij)])%*%padroes[,faixa])))
  }
  -veross  
}
expito<-function(eta) {
  1/(1+exp(-eta))
}
b<-c(0,0,0,1)
B<-rbind(diag(3),rep(-1,3))

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

mnar1<-function(a) {
  i<-c(1,1,2,2)
  j<-c(1,2,1,2)
  psi1 <-expito(a[1]+a[4]*(i-1)+a[5]*(j-1))
  psi21<-expito(a[2]+a[4]*(i-1)+a[5]*(j-1))
  psi20<-expito(a[3]+a[4]*(i-1)+a[5]*(j-1))
  c(psi1*psi21,psi1*(1-psi21),(1-psi1)*psi20)
}
fit.mnar1<-optim(par=c(mar.ml$theta[1:3],rep(0,5)),fn=mveross,method="L-BFGS-B",lower=c(rep(0,3),rep(-Inf,5)),upper=c(rep(1,3),rep(Inf,5)),hessian=TRUE,
  control=list(maxit=1000),TF=TF,padroes=padroes,classes=classes,estrutura=mnar1)
#ou
fit.mnar1.ini<-constrOptim(theta=c(mar.ml$theta[1:3],rep(0,5)),f=mveross,method="Nelder-Mead",ui=cbind(B,matrix(0,4,5)),ci=-b+1e-16,
  TF=TF,padroes=padroes,classes=classes,estrutura=mnar1) #a constrOptim não fornece a hessiana numérica, então..
fit.mnar1<-optim(par=fit.mnar1.ini$par,fn=mveross,method="Nelder-Mead",hessian=TRUE,control=list(maxit=1),
  TF=TF,padroes=padroes,classes=classes,estrutura=mnar1)

loglin.mnar1<-funlinWLS(model=c("lin","log"),theta=c(b+B%*%fit.mnar1$par[1:3]),Vtheta=B%*%solve(fit.mnar1$hess)[1:3,1:3]%*%t(B),A1=t(c(1,-1,-1,1)),X=1)
exp(loglin.mnar1$beta) #Estimativa da razão de chances, MNAR1, Tabela 1.22
exp(loglin.mnar1$beta+c(-1,1)*qnorm(0.975)*sqrt(loglin.mnar1$Vbeta)) #IC(95%,razão de chances), MNAR1, Tabela 1.22
round(sum(TF)*gammatij(pij=fit.mnar1$par[1:3],lambdat.ij=mnar1(fit.mnar1$par[-(1:3)]))$gtij,0)  #Tabela 1.23, MNAR1
2*(sat.lv+fit.mnar1$value) #Estatística de razão de verossimilhanças do ajuste do modelo



mnar2<-function(a) {
  i<-c(1,1,2,2)
  j<-c(1,2,1,2)
  psi1 <-expito(a[1]+a[3]*(i-1)+a[4]*(j-1)+a[5]*(i-1)*(j-1))
  psi21<-expito(a[2]+a[3]*(i-1)+a[4]*(j-1)+a[5]*(i-1)*(j-1))
  psi20<-expito(a[2]+a[3]*(i-1)+a[4]*(j-1)+a[5]*(i-1)*(j-1))
  c(psi1*psi21,psi1*(1-psi21),(1-psi1)*psi20)
}
fit.mnar2<-optim(par=c(mar.ml$theta[1:3],rep(0,5)),fn=mveross,method="L-BFGS-B",lower=c(rep(0,3),rep(-Inf,5)),upper=c(rep(1,3),rep(Inf,5)),hessian=TRUE,
  control=list(maxit=1000),TF=TF,padroes=padroes,classes=classes,estrutura=mnar2)

loglin.mnar2<-funlinWLS(model=c("lin","log"),theta=c(b+B%*%fit.mnar2$par[1:3]),Vtheta=B%*%solve(fit.mnar2$hess)[1:3,1:3]%*%t(B),A1=t(c(1,-1,-1,1)),X=1)
exp(loglin.mnar2$beta) #Estimativa da razão de chances, MNAR1, Tabela 1.22
exp(loglin.mnar2$beta+c(-1,1)*qnorm(0.975)*sqrt(loglin.mnar2$Vbeta)) #IC(95%,razão de chances), MNAR1, Tabela 1.22
round(sum(TF)*gammatij(pij=fit.mnar2$par[1:3],lambdat.ij=mnar2(fit.mnar2$par[-(1:3)]))$gtij,0)  #Tabela 1.23, MNAR1
2*(sat.lv+fit.mnar2$value) #Estatística de razão de verossimilhanças do ajuste do modelo



mveross2<-function(p,TF,padroes,classes,estrutura,pars) {
  gtij<-gammatij(pij=expit(c(sum(p[1:3]),p[1],p[2])),lambdat.ij=estrutura(p[-(1:3)],pars) )
  totp<-ncol(padroes)
  if(totp!=sum(classes)){
    stop("classes e padroes nao sao compativeis")
  }
  veross<-0
  for(ind in 1:(gtij$tottij/gtij$totij)) {
    faixa<-(1+ifelse(ind==1,0,sum(classes[1:(ind-1)]))):sum(classes[1:ind])
    veross<-veross+sum(TF[faixa]*log(diag(t(padroes[,faixa])%*%diag(gtij$gtij[(1+(ind-1)*gtij$totij):(ind*gtij$totij)])%*%padroes[,faixa])))
  }
  -veross  
}
expit<-function(p) {
  exp(p)/(1+sum(exp(p)))
}
mnar3<-function(a,a3) {
  i<-c(1,1,2,2)
  j<-c(1,2,1,2)
  psi1 <-expito(a[1]+a[4]*(i-1)+a[5]*(j-1)+a3*(i-1)*(j-1))
  psi21<-expito(a[2]+a[4]*(i-1)+a[5]*(j-1)+a3*(i-1)*(j-1))
  psi20<-expito(a[3]+a[4]*(i-1)+a[5]*(j-1)+a3*(i-1)*(j-1))
  c(psi1*psi21,psi1*(1-psi21),(1-psi1)*psi20)
}

pnm<-seq(-5,5,0.01) #Faixa de variação do parâmetro de sensibilidade
lpn<-length(pnm)
res<-matrix(0,lpn,4)

date()
for(i in 1:lpn) { #Quase 20min. num Intel Core 2 Duo 1.7GHz
  pn<-pnm[i]
  fit.mnar3<-nlm(p=rep(0,8),f=mveross2,hessian=TRUE,iterlim=10000,TF=TF,padroes=padroes,classes=classes,estrutura=mnar3,
    pars=pn)
  res[i,1]<-exp(fit.mnar3$est[3]) #Estimativa da razão de chances, MNAR3
  res[i,2:3]<-exp(fit.mnar3$est[3]+c(-1,1)*qnorm(0.975)*sqrt(solve(fit.mnar3$hess)[3,3])) #IC(95%,razão de chances), MNAR3
  res[i,4]<-2*(sat.lv+fit.mnar3$min)
}
date()

plot(cbind(pnm,pnm,pnm),res[,1:3],xlab=expression(paste(alpha[3])), ylab="Razão de chances",pch=16,type="n")
lines(pnm,res[,1],lwd=1)
lines(pnm,res[,2],lwd=3)
lines(pnm,res[,3],lwd=3)
abline(1,0,lty=2)

c(min(res[,1]),max(res[,1])) #Intervalo de ignorância
c(min(res[,2]),max(res[,3])) #Intervalo de 95% de incerteza

max(res[,4]) #Maior valor da estatística de razão de verossimilhanças do ajuste do modelo
plot(pnm,res[,4],xlab=expression(paste(alpha[3])), ylab="Estatística de razão de verossimilhanças do ajuste do modelo",pch=16,type="l")