Regresión lineal con R2OpenBUGS

Regresión lineal

En clase, hemos visto como ajustar un modelo de regresión lineal a los datos de velocidades y distancias necesarias para para unos coches. Ahora ilustramos como usar para hacer el ajuste bayesiano. Primero vemos el ajuste clásico.

data(cars)
modelolineal <- lm(dist ~ speed, data=cars)
print(modelolineal)

## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Coefficients:
## (Intercept)        speed  
##     -17.579        3.932

plot(cars)
abline(modelolineal,col='red',lwd=2)

xnuevo <- data.frame(20)
names(xnuevo) <- "speed"
predict(modelolineal,xnuevo,interval="confidence")

##        fit      lwr      upr
## 1 61.06908 55.24729 66.89088

predict(modelolineal,xnuevo,interval="prediction")

##        fit      lwr      upr
## 1 61.06908 29.60309 92.53507

Modelización en OpenBUGS

Podemos construir un modelo de regresión utilizando un doodle en OpenBUGS.

Para correrlo en R, primero escribimos una función con el código del modelo.

modelobayes <- function(){
  for( i in 1 : n ) {
        dist[i] ~ dnorm(mu[i], phi)
        mu[i] <- beta0 + beta1 * speed[i]
    }
    beta0 ~ dnorm(0.0, 1.0E-6)
    beta1 ~ dnorm(0.0, 1.0E-6)
    phi ~ dgamma(0.001, 0.001)
    sigma2 <- 1 / phi
}

Definición de los datos

Para correr el programa a través de R, es necesario definir los datos en forma de una lista.

datos <- list(n=length(cars$dist),dist=cars$dist,speed=cars$speed)

Valores iniciales

Metemos los valores iniciales en otra lista dentro de una función.

iniciales <- function(){
  list(beta0=0,beta1=0,phi=1)
}

Ejecución del código

Ya llamamos a y usamos el comando para ejecutar el programa. El código corre 10000 iteraciones en equilibrio y 5000 de burnin.

library(R2OpenBUGS)
resultados <- bugs(data = datos, inits = iniciales, parameters.to.save = c("beta0", 
    "beta1", "phi","sigma2"), model.file = modelobayes, n.chains = 1, n.iter = 10000)

Mirando los resultados

Podemos obtener un breve resumen de los resultados.

resultados 

## Inference for Bugs model at "C:/Users/mwiper/AppData/Local/Temp/RtmpWafHvR/model2cac5e14983.txt", 
## Current: 1 chains, each with 10000 iterations (first 5000 discarded)
## Cumulative: n.sims = 5000 iterations saved
##           mean   sd  2.5%   25%   50%   75% 97.5%
## beta0    -18.1  7.1 -32.3 -22.6 -18.0 -13.4  -4.3
## beta1      4.0  0.4   3.1   3.7   4.0   4.2   4.8
## phi        0.0  0.0   0.0   0.0   0.0   0.0   0.0
## sigma2   247.2 53.9 163.1 210.1 240.1 277.5 368.8
## deviance 416.4  2.6 413.4 414.5 415.7 417.5 423.1
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 3.2 and DIC = 419.5
## DIC is an estimate of expected predictive error (lower deviance is better).

En este caso, las estimaciones medias de los regresores son parecidas a las frecuentistas.

Análisis de los datos simulados

Podemos analizar los resultados a través del paquete , pero es casi más fácil hacerlo directamente, sacando los datos generados por el MCMC en R.

names(resultados$sims.list)

## [1] "beta0"    "beta1"    "phi"      "sigma2"   "deviance"

beta0 <- resultados$sims.list$beta0
beta1 <- resultados$sims.list$beta1
phi <- resultados$sims.list$phi
sigma2 <- resultados$sims.list$sigma2
head(sigma2)

## [1] 299.7 304.1 249.4 212.9 177.3 212.6

head(1/phi)

## [1] 299.7602 304.1363 249.3766 212.9019 177.3364 212.5850

Se puede ver que los valores generados de \(\sigma^2\) son los mismos que uno partido por los valores generados de \(\phi\).

Convergencia

Vemos si el algoritmo ha convergido bien.

ts.plot(beta0)

ts.plot(cumsum(beta0)/c(1:length(beta0)))

acf(beta0)

Parece que sí

Densidades a posteriori

Es fácil ver las caracteristicas de las distribuciones a posteriori.

plot(density(beta0),xlab=expression(beta[0]),ylab='f',main='')

mean(beta0)

## [1] -18.05206

quantile(beta0,c(0.025,0.5,0.975))

##      2.5%       50%     97.5% 
## -32.27075 -18.01500  -4.26745

Predicción

También podemos hacer predicciones.

mediapred <- beta0+beta1*20
quantile(mediapred,c(0.025,0.5,0.975))

##     2.5%      50%    97.5% 
## 55.38047 61.10000 67.12000

distpred <- mediapred+rnorm(length(sigma2),0,sqrt(sigma2))
quantile(distpred,c(0.025,0.5,0.975))

##     2.5%      50%    97.5% 
## 30.16503 60.71175 93.87650

Los resultados son parecidos a los resultados clásicos.