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
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)
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)
}
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)
}
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.
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")
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
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
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.