#
# Modelos lineales generales para datos binarios
#
rm(list=ls())
#
# Datos sobre admisiones de estudiantes en UCLA.  
# Cambiar la siguiente linea con el directorio donde has guardado los datos.
#
setwd("Z:/docencia/Spanish/Inferencia Bayesiana/R Codes")
df <- read.csv("binary.csv")
#
# Los datos son admit (1 = admitido, 0 no), gre and gpa que son medidas de lo bueno del estudiante y rank que es 
# un ranking del colegio donde hicieron sus estudios.
#
str(df)
summary(df)
xtabs(~ admit + rank ,data=df)
#
# Inferencia clásica via regresión logistica.
#
logit <- glm(admit~gre+gpa+as.factor(rank),data=df,family=binomial)
summary(logit)
#
# Predicción.
#
x <- data.frame(gre=790,gpa=3.8,rank=as.factor(1))
predict(logit,x)
#
# Inferencia via regresión probit.
#
probit <- glm(admit~gre+gpa+as.factor(rank),data=df,family=binomial(link="probit"))
summary(probit)
predict(probit,x)
#
# Inferencia bayesiana con MCMC.
#
if (!require("MCMCpack")) install.packages("MCMCpack")
library(MCMCpack)
burnin <- 10000
mcmc <- 20000
logitmcmc <- MCMClogit(admit~gre+gpa+as.factor(rank),data=df, burnin = burnin, mcmc = mcmc,b0=0,B0=0.00001)
summary(logitmcmc)
densplot(logitmcmc) # Ilustrar las densidades a posteriori
intercepto <- logitmcmc[,1]
betagre <- logitmcmc[,2]
betagpa <- logitmcmc[,3]
betarank <- matrix(rep(0,mcmc*4),nrow=mcmc)
betarank[,2:4] <- logitmcmc[4:6]
predmcmc <- 1/(1+exp(-(intercepto+betagre*790+betagpa*3.8)))
plot(density(predmcmc))
mean(predmcmc)
#
# Miramos la convergencia.
#
ts.plot(intercepto)
ts.plot(cumsum(intercepto)/c(1:mcmc))
acf(intercepto)
#
# No se converge muy bien.  Intentamos el mcmc con thinning=30.
#
thinning <- 30
burnin <- burnin*thinning
mcmc <- mcmc*thinning
logitmcmc <- MCMClogit(admit~gre+gpa+as.factor(rank),data=df, burnin = burnin, mcmc = mcmc,b0=0,B0=0.00001,thin=thinning)
summary(logitmcmc)
densplot(logitmcmc) # Ilustrar las densidades a posteriori
intercepto <- logitmcmc[,1]
betagre <- logitmcmc[,2]
betagpa <- logitmcmc[,3]
betarank <- matrix(rep(0,mcmc*4),nrow=mcmc)
betarank[,2:4] <- logitmcmc[4:6]
predmcmc <- 1/(1+exp(-(intercepto+betagre*790+betagpa*3.8)))
plot(density(predmcmc),main="Densidad predictiva",ylab="f")
mean(predmcmc)
#
# Regresión probit.
#
probitmcmc <- MCMCprobit(admit~gre+gpa+as.factor(rank),data=df, burnin = burnin, mcmc = mcmc,b0=0,B0=0.00001,thin=thinning)
summary(probitmcmc)
densplot(probitmcmc) # Ilustrar las densidades a posteriori
intercepto <- probitmcmc[,1]
betagre <- probitmcmc[,2]
betagpa <- probitmcmc[,3]
betarank <- matrix(rep(0,mcmc*4),nrow=mcmc)
betarank[,2:4] <- probitmcmc[4:6]
predmcmc <- pnorm(intercepto+betagre*790+betagpa*3.8)
plot(density(predmcmc),main="Densidad predictiva",ylab="f")
mean(predmcmc)