##########################################################################################
#
# DADOS MUNICIPAIS BRASILEIROS (ipeadata.gov.br)
#
##########################################################################################
#
#  1. Estado
#  2. Codigo 
#  3. Municipio
#  4. Renda domiciliar per capita - media do 1o quintil (R$ por mes)
#  5. Renda domiciliar per capita - media do 2o quintil (R$ por mes)	
#  6. Renda domiciliar per capita - media do 3o quintil (R$ por mes)
#  7. Renda domiciliar per capita - media do 4o quintil (R$ por mes)	
#  8. Renda domiciliar per capita - media do 5o quintil (R$ por mes)	
#  9. Razao entre a renda dos 10% mais ricos e 40% mais pobres 	
# 10. Mulheres chefes de familia sem conjuge e com filhos menores de 15 anos (%)	
# 11. Medicos residentes (por mil habitantes)
# 12. Enfermeiros residentes com curso superior (%)	
# 13. Alfabetizados - pessoas 15 anos e mais (%) 	
# 14. Indice de Desenvolvimento Humano 	
# 15. Indice de Desenvolvimento Humano - renda 	
# 16. Indice de Desenvolvimento Humano - longevidade 	
# 17. Indice de Desenvolvimento Humano - educacao 
# 18. Mortalidade ate cinco anos de idade (por mil nascidos vivos)
# 19. Probabilidade de sobrevivencia ate 40 anos (%)
# 20. Probabilidade de sobrevivencia ate 60 anos (%)
# 21. Taxa de fecundidade (%)
# 22. Pessoas 65 anos ou mais - morando sozinhas (%)
# 23. Pessoas 10 e 14 anos - mulheres com filhos (%)	
# 24. Professores do fundamental residentes com curso superior (%)
# 25. Esperanca de vida ao nascer
# 26. Mortalidade infantil (por mil nascidos vivos)
# 27. Domicilios - com agua encanada - pessoas (%)
# 28. Domicilios - com servico de coleta de lixo - pessoas (%)
#
##########################################################################################
rm(list=ls())

nomes = c("Renda per capita - media do 2o quintil (LOG)",	
"Renda per capita - media do 5o quintil (LOG)",
"Alfabetizados - 15 ou mais",
"Mortalidade ate 5 anos (por mil vivos)",
"Prob. sobrev. ate 40",
"Prob. sobrev. ate 60",
"Esperanca de vida ao nascer",
"Mortalidade infantil (por mil vivos)",
"Agua encanada (LOGIT)",
"Coleta de lixo (LOGIT)")

data = read.csv("dadosmunicipais-ipea.csv",header=TRUE)
n  = nrow(data)
attach(data)

reg = lm(idh~rendapc1+rendapc2+rendapc3+rendapc4+rendapc5+razao1040+mulherchefe+medico+enfermeiro+alfabetizado+mortal5+prob40+prob60+taxafecund+LXVoumais+maes10a14+professores+esperanca+mortal+agua+lixo)

plot(reg$fit,reg$res)

lrendapc2       = log(rendapc2)
lrendapc5       = log(rendapc5)
agua[agua==0]   = min(agua[agua>0])
agua[agua==100] = max(agua[agua<100])
lixo[lixo==0]   = min(lixo[lixo>0])
lixo[lixo==100] = max(lixo[lixo<100])
lagua           = log(agua/(100-agua))
llixo           = log(lixo/(100-lixo))

pdf(file="dadosmunicipais.pdf",width=12,height=6)

par(mfrow=c(2,5))
plot(lrendapc2,idh,xlab=nomes[1],ylab="IDH")
abline(lm(idh~lrendapc2),lwd=2,col=2)
plot(lrendapc5,idh,xlab=nomes[2],ylab="IDH")
abline(lm(idh~lrendapc5),lwd=2,col=2)
plot(alfabetizado,idh,xlab=nomes[3],ylab="IDH")
abline(lm(idh~alfabetizado),lwd=2,col=2)
plot(mortal5,idh,xlab=nomes[4],ylab="IDH")
abline(lm(idh~mortal5),lwd=2,col=2)
plot(prob40,idh,xlab=nomes[5],ylab="IDH")
abline(lm(idh~prob40),lwd=2,col=2)
plot(prob60,idh,xlab=nomes[6],ylab="IDH")
abline(lm(idh~prob60),lwd=2,col=2)
plot(esperanca,idh,xlab=nomes[7],ylab="IDH")
abline(lm(idh~esperanca),lwd=2,col=2)
plot(mortal,idh,xlab=nomes[8],ylab="IDH")
abline(lm(idh~mortal),lwd=2,col=2)
plot(lagua,idh,xlab=nomes[9],ylab="IDH")
abline(lm(idh~lagua),lwd=2,col=2)
plot(llixo,idh,xlab=nomes[10],ylab="IDH")
abline(lm(idh~llixo),lwd=2,col=2)

estados = sort(unique(data[,1]))
ne      = length(estados)
coefs   = array(0,c(ne,10,2))
par(mfrow=c(2,5))
for (i in 1:ne){
  ii    = sigla==estados[i]
  reg1  = lm(idh[ii]~lrendapc2[ii])
  reg2  = lm(idh[ii]~lrendapc5[ii])
  reg3  = lm(idh[ii]~alfabetizado[ii])
  reg4  = lm(idh[ii]~mortal5[ii])
  reg5  = lm(idh[ii]~prob40[ii])
  reg6  = lm(idh[ii]~prob60[ii])
  reg7  = lm(idh[ii]~esperanca[ii])
  reg8  = lm(idh[ii]~mortal[ii])
  reg9  = lm(idh[ii]~lagua[ii])
  reg10 = lm(idh[ii]~llixo[ii])
  plot(lrendapc2[ii],idh[ii],main=estados[i],xlab=nomes[1],ylab="IDH")
  abline(reg1$coef,col=2)
  plot(lrendapc5[ii],idh[ii],xlab=nomes[2],ylab="IDH")
  abline(reg2$coef,col=2)
  plot(alfabetizado[ii],idh[ii],xlab=nomes[3],ylab="IDH")
  abline(reg3$coef,col=2)
  plot(mortal5[ii],idh[ii],xlab=nomes[4],ylab="IDH")
  abline(reg4$coef,col=2)
  plot(prob40[ii],idh[ii],xlab=nomes[5],ylab="IDH")
  abline(reg5$coef,col=2)
  plot(prob60[ii],idh[ii],xlab=nomes[6],ylab="IDH")
  abline(reg6$coef,col=2)
  plot(esperanca[ii],idh[ii],xlab=nomes[7],ylab="IDH")
  abline(reg7$coef,col=2)
  plot(mortal[ii],idh[ii],xlab=nomes[8],ylab="IDH")
  abline(reg8$coef,col=2)
  plot(lagua[ii],idh[ii],xlab=nomes[9],ylab="IDH")
  abline(reg9$coef,col=2)
  plot(llixo[ii],idh[ii],xlab=nomes[10],ylab="IDH")
  abline(reg10$coef,col=2)
  coefs[i,1,] = reg1$coef
  coefs[i,2,] = reg2$coef
  coefs[i,3,] = reg3$coef
  coefs[i,4,] = reg4$coef
  coefs[i,5,] = reg5$coef
  coefs[i,6,] = reg6$coef
  coefs[i,7,] = reg7$coef
  coefs[i,8,] = reg8$coef
  coefs[i,9,] = reg9$coef
  coefs[i,10,] = reg10$coef
}

# Comparing several models
X = cbind(lrendapc2,lrendapc5,alfabetizado,mortal5,prob40,prob60,esperanca,
          mortal,lagua,llixo)
model = NULL
for (i1 in 0:1)
for (i2 in 0:1)
for (i3 in 0:1)         
for (i4 in 0:1)         
for (i5 in 0:1)         
for (i6 in 0:1)         
for (i7 in 0:1)         
for (i8 in 0:1)         
for (i9 in 0:1)         
for (i10 in 0:1)
  model = rbind(model,c(i1,i2,i3,i4,i5,i6,i7,i8,i9,i10))
model = model[2:1024,]
  
SST    = sum((idh-mean(idh))^2)  
ind    = 1:10
nmodel = 1023
R2a    = rep(0,nmodel)
npar   = rep(0,nmodel)
for (i in 1:nmodel){
  k      = sum(model[i,])
  reg    = lm(idh~X[,ind[model[i,]==1]])
  R2     = 1-sum(reg$res^2)/SST
  R2a[i] = 1-(1-R2)*(n-1)/(n-k-1)
  npar[i] = 1+k
}

topmodels = NULL
R2max     = NULL
for (i in 2:10){
  model1 = model[npar==i,]
  R2a1   = R2a[npar==i]
  R2max  = c(R2max,max(R2a1))
  topmodels = rbind(topmodels,model1[R2a1==max(R2a1),])
}
R2max  = c(R2max,R2a[nmodel])
topmodels = rbind(topmodels,model[nmodel,])

nomes1 = c("lrendapc2","lrendapc5","alfabet","mortal5",
           "prob40","prob60","esperanca","mortal","lagua","llixo")

par(mfrow=c(1,1))
plot(npar,R2a,xlab="Numero de regressores (incluindo o intercepto)",ylab="R2 ajustado",
     pch=16,ylim=c(min(R2a),1.05),axes=FALSE)
axis(2);box();axis(1,at=2:11)
for (i in 2:11)
  text(i,1.02,round(100*R2max[i-1],3),col=4)
text(6.5,1.05,"R2a maximo",col=4)
abline(h=1,lty=2)

par(mfrow=c(1,1))
plot(c(1,11),c(0,13),axes=FALSE,col=0,xlab="",ylab="")
for (i in 1:10){
 text(1:10,11-i,topmodels[i,])
 text(11,11-i,round(R2max[i],4))
 text(i,12,nomes1[i])
 text(i,11,paste(round(100*mean(topmodels[,i]),1),"%",sep=""))
}
text(11,12,"R2a")
abline(h=12.5)
abline(h=11.5)
abline(h=10.5)
abline(h=0.5)

reg.r = lm(idh~lrendapc5+alfabetizado+prob60)
summary(reg.r)

reg.f = lm(idh~lrendapc2+lrendapc5+alfabetizado+mortal5+prob40+prob60+esperanca+
          mortal+lagua+llixo)
summary(reg.f)

par(mfrow=c(2,2))
U = max(abs(reg.r$res))
L = -U
xxx = seq(L,U,length=40)
hist(reg.r$res,prob=TRUE,xlab="",ylab="Residuals",main="Restrited model",breaks=xxx)
xxx = seq(L,U,length=1000)
lines(xxx,dnorm(xxx,mean(reg.r$res),sqrt(var(reg.r$res))),col=2)
box()
plot(reg.r$fit,reg.r$res,xlab="Fitted",ylab="Residuals",ylim=c(L,U))
abline(h=0,col=2)

U = max(abs(reg.f$res))
L = -U
xxx = seq(L,U,length=40)
hist(reg.f$res,prob=TRUE,xlab="",ylab="Residuals",main="Full model",breaks=xxx)
xxx = seq(L,U,length=1000)
lines(xxx,dnorm(xxx,mean(reg.f$res),sqrt(var(reg.f$res))),col=2)
box()
plot(reg.f$fit,reg.f$res,xlab="Fitted",ylab="Residuals",ylim=c(L,U))
abline(h=0,col=2)

dev.off()