Introducción

Repetimos el análisis del modelo para la probabilidad de admisión a la universidad de la última clase. En este caso, hacemos sólo el modelo con el vínculo logit.

rm(list=ls())
df <- read.csv("http://halweb.uc3m.es/esp/Personal/personas/mwiper/docencia/Spanish/Inferencia Bayesiana/R Codes/binary.csv") 
logit <- glm(admit~gre+gpa+as.factor(rank),data=df,family=binomial)
summary(logit)
## 
## Call:
## glm(formula = admit ~ gre + gpa + as.factor(rank), family = binomial, 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6268  -0.8662  -0.6388   1.1490   2.0790  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -3.989979   1.139951  -3.500 0.000465 ***
## gre               0.002264   0.001094   2.070 0.038465 *  
## gpa               0.804038   0.331819   2.423 0.015388 *  
## as.factor(rank)2 -0.675443   0.316490  -2.134 0.032829 *  
## as.factor(rank)3 -1.340204   0.345306  -3.881 0.000104 ***
## as.factor(rank)4 -1.551464   0.417832  -3.713 0.000205 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.98  on 399  degrees of freedom
## Residual deviance: 458.52  on 394  degrees of freedom
## AIC: 470.52
## 
## Number of Fisher Scoring iterations: 4
x <- data.frame(gre=790,gpa=3.8,rank=as.factor(1))
predict(logit,x)
##       1 
## 0.85426

Definiendo los datos en el formato adecuado para R2OpenBUGS

Para implementar la inferencia en R2OpenBUGS, primero definimos las variables

admit <- df$admit
gre <- df$gre
gpa <- df$gpa
rank <- df$rank
n <- length(rank)

Ahora tenemos que definir dummies para los ranks 2 a 4.

rank2 <- rep(0,n)
rank3 <- rep(0,n)
rank4 <- rep(0,n)
rank2[rank==2] <- 1
rank3[rank==3] <- 1
rank4[rank==4] <- 1

Finalmente, ponemos los datos en un .

datos <- list(n=n,admit=admit,gre=gre,gpa=gpa,rank2=rank2,rank3=rank3,rank4=rank4)

Definiendo el modelo BUGS.

Podemos construir la siguiente representación gráfica del modelo logístico en OpenBUGS. Observamos que tenemos que definir el vínculo para \(\mu_i\) como logit y tenemos que representar la distribución de \(admit_i\) como bern.

Generamos el siguiente código que representamos en una función en R. Observamos que hemos dejado las distribuciones a priori poco informativas.

mlg <- function(){
  for( i in 1 : n ) {
        admit[i] ~ dbern(p[i])
        logit(p[i]) <- beta0+gre[i]*betagre+gpa[i]*betagpa+rank2[i]*betarank2+rank3[i]*betarank3 
        +rank4[i]*betarank4
    }
    beta0 ~ dnorm(0.0, 1.0E-6)
    betagpa ~ dnorm(0.0, 1.0E-6)
    betagre ~ dnorm(0.0, 1.0E-6)
    betarank2 ~ dnorm(0.0, 1.0E-6)
    betarank3 ~ dnorm(0.0, 1.0E-6)
    betarank4 ~ dnorm(0.0, 1.0E-6)
}

Definiendo los valores iniciales.

Una opción razonable para definir valores iniciales es suponer que los regresores no son útiles, cuando la probabilidad de admisión sería constante. Estimamos \(\beta_0\) a través de está probabilidad. (Otra opción sería definir las estimaciones frecuentistas de los parámetros como valores iniciales.)

p <- length(admit[admit==1])/n
p
## [1] 0.3175
log(p/(1-p))
## [1] -0.7652847
inits <- function(){
  list(beta0=-0.765,betagpa=0,betagre=0,betarank2=0,betarank3=0,betarank4=0)
}

Corriendo el código

Usamos el comando bugs para correr el código.

library(R2OpenBUGS)
mlgout <- bugs(data=datos,inits=inits,parameters.to.save=c("beta0","betagre","betagpa","betarank2","betarank3","betarank4"),
               model.file=mlg,n.chains = 1,n.burn=10000,n.iter = 20000)
mlgout
## Inference for Bugs model at "C:/Users/mwiper/AppData/Local/Temp/RtmpspKJXA/model2a7c10dd1f20.txt", 
## Current: 1 chains, each with 20000 iterations (first 10000 discarded)
## Cumulative: n.sims = 10000 iterations saved
##            mean  sd  2.5%   25%   50%   75% 97.5%
## beta0      -4.1 1.2  -6.4  -4.9  -4.1  -3.3  -1.8
## betagre     0.0 0.0   0.0   0.0   0.0   0.0   0.0
## betagpa     0.8 0.3   0.2   0.6   0.8   1.1   1.5
## betarank2  -0.7 0.3  -1.3  -0.9  -0.7  -0.5   0.0
## betarank3  -1.4 0.4  -2.1  -1.6  -1.4  -1.1  -0.7
## betarank4  -1.6 0.4  -2.4  -1.9  -1.6  -1.3  -0.8
## deviance  464.7 3.6 459.8 462.1 464.0 466.6 473.4
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 6.1 and DIC = 470.8
## DIC is an estimate of expected predictive error (lower deviance is better).

Vemos que las estimaciones media de los parámetros son similares a las estimaciones frecuentistas.

Análisis de convergencia

En el siguiente código miramos la convergencia del \(\beta_0\). En este caso, parece que la convergencia es bastante razonable. (En principio debemos mirar las demás variables también).

beta0 <- mlgout$sims.list$beta0
ts.plot(beta0,xlab='iters',ylab=expression(beta[0]))

ts.plot(cumsum(beta0)/c(1:length(beta0)),xlab='iters',ylab="E[beta0|datos]")

acf(beta0,main="ACF de beta0")

Densidades a posteriori

Podemos graficar las densidades a posteriori de los parámetros y mostrar estimaciones más precisas de la media y cuantiles.

plot(density(beta0),xlab=expression(beta0),ylab='f',main='Densidad a posteriori de beta0')

mean(beta0)
## [1] -4.07143
quantile(beta0,c(0.025,0.5,0.975))
##     2.5%      50%    97.5% 
## -6.36615 -4.07750 -1.82495

Predicción

Podemos calcular la densidad predictiva de la probabilidad de admisión del nuevo alumno.

betagre <- mlgout$sims.list$betagre
betagpa <- mlgout$sims.list$betagpa
mupred <- beta0+betagre*790+betagpa*3.8
padmit <- exp(mupred)/(1+exp(mupred))
plot(density(padmit),xlab='p',ylab='f',main='Probabilidad de admisión del nuevo alumno')

mean(padmit)
## [1] 0.6997353
quantile(padmit,c(0.025,0.5,0.975))
##      2.5%       50%     97.5% 
## 0.5575696 0.7036618 0.8214144

Regresión probit

Podemos también intentar correr el modelo con el vínculo probit. (En principio, sólo se necesita cambiar logit por probit en la función mlg y cambiar el valor inicial de \(\beta_0\) para corresponder con el probit. No obstante, en este caso, da algunos problemas numéricos del algoritmo.