Regresiones Simples con JAGS

En este práctico vamos a ajustar regresiones lineales. Para eso tenemos que aprender a escribir las correspondientes funciones de likelihood y aprender a definir las previas de los parámetros.

Objetivos:

• Familiarizarse con la formulación de funciones de likelihood
• Simular datos a partir de distintas funciones de likelihood
• Usar JAGS para estimar parámetros a partir de los datos simulados

Empecemos con una relación lineal entre una variable predictora y una variable respuesta, y asumiendo una distribución normal para la variabilidad. Este es básicamente el modelo detrás de una regresión lineal simple. Más formalmente, podemos escribir el modelo de la siguiente manera:

\[ \begin{aligned} & y_i \sim N(\mu_i, \sigma^{2}) \\ & \mu_i = \beta_0 + \beta_1 \times x_i\\ & \textrm{para } i \textrm{ de } 1:n \\ \end{aligned} \]

Para simular datos de acuerdo con este modelo, primero tenemos que definir los parámetros (pendiente, ordenada al origen y desvío estándar), la variable predictora (también llamada covariable), y el número de observaciones. Noten que la distribución normal tiene dos parámetros, media y desvío. El “truco” para armar un modelo de regresión es hacer que la media dependa de la variable predictora. Como ahora la media es una función lineal de \(x\), el modelo tiene tres parámetros (\(\beta_0\), \(\beta_1\), y \(\sigma\)).

Definamos los parámetros, el número de observaciones y simulemos valores para la variable predictora (en este caso asumimos valores de una distribución uniforme entre \(0\) y \(20\)):

set.seed(2345)
b0 <- 2  # Ordenada al origen (valor de la media cuando la predictora vale cero)
b1 <- 0.2  # Pendiente (cuánto cambia la media cuando la predictora cambia en una unidad)
s <- 3  # Desvío estándar
n <- 10  # Número de observaciones
x.sim <- runif(n, 0, 20)  # covariable simulada

Ahora podemos calcular el valor esperado (\(\mu_i\)) para cada valor de \(x\) y simular datos con una distribución Normal. Esto es hasta ahora igual a lo que hicimos en el práctico de probabilidades

m <- b0 + b1 * x.sim  # media de la distribución normal
y.sim <- rnorm(n, mean = m, sd = s)  # datos simulados

Este modelo se puede ajustar fácilmente en R usando el comando lm, pero vamos a tomarnos el trabajo de escribir la función de likelihood y ajustarlo con JAGS. Para esto último tenemos que escribir el modelo en lenguaje BUGS, especificando la función de likelihood y las previas para todos los valores “no observados”, que en este caso son los parámetros. Una cuestión importante es que en JAGS la distribución normal está formulada en base a la precisión en vez de desvío estándar. La precisión es \(1/\sigma^2\) y normalmente se representa con la letra griega \(\tau\) (tau). Vemos entonces que en el modelo que vamos a definir, pusimos una previa uniforme entre \(0\) y \(100\) sobre el desvío estándar (s) y luego transformamos a tau. Además, las previas que ponemos sobre b0 y b1 equivalen a una normal con media \(0\) y varianza \(=1/0.1 = 10\).

cat(file = "ml.bug", "
      model  {
           # Función de likelihood
           for( i in 1:n){
                y[i] ~ dnorm (mu[i], tau) 
                mu[i] <- b0 + b1 * x[i]
           }

           # Previas
           b0 ~ dnorm(0, 0.1) 
           b1 ~ dnorm(0, 0.1)
           tau <- 1/(s*s) 
           s ~ dunif(0, 100)
    }
    ")

Pregunta: ¿Qué tan informativas son estas previas?

Como vimos antes, para llamar a JAGS, cargamos el paquete jagsUI, definimos los datos, los parámetros y luego usamos la función jags.

library(jagsUI)

# defino las variables
x <- x.sim
y <- y.sim

# lista con los datos para pasarle a JAGS
datos <- list("x", "y", "n")

# vector con los parámetros para que JAGS guarde los valores de las cadenas
# Markovianas
params <- c("b0", "b1", "s")

# función para generar valores iniciales de los parámetros (no confundir con
# las previas!)
inits <- function() list(b0 = runif(1, 0, 1), b1 = runif(1, 0, 1), s = runif(1, 
    0, 10))

Finalmente, definimos cuántas iteraciones queremos correr por cadena, cuántos valores vamos a ir descartado en una secuencia (“thin”), qué largo tiene el “burn in” y cuántas cadenas queremos correr.

ni <- 1000
nt <- 1
nb <- 500
nc <- 3

Ahora llamamos a JAGS para que genere las cadenas Markovianas

ml.sim <- jags(data = datos, inits = inits, parameters.to.save = params, model.file = "ml.bug", 
    n.chains = nc, n.iter = ni, n.burnin = nb, n.thin = nt)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 10
##    Unobserved stochastic nodes: 3
##    Total graph size: 53
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 500 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 500 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

Ahora veamos los resultados. Recuerden que primero tenemos que fijarnos si las cadenas convergieron (Rhat \(\leq 1.1\)) y además chequear que el n.eff sea suficientemente grande como para estimar con confianza propiedades de las posteriores (con un n.eff del orden de “cientos”, podemos estimar media e intervalos de credibilidad aproximados. Algunas revistas son bastante exigentes y piden n.eff \(\geq 4000\)).

print(ml.sim)

## JAGS output for model 'ml.bug', generated by jagsUI.
## Estimates based on 3 chains of 1000 iterations,
## burn-in = 500 iterations and thin rate = 1,
## yielding 1500 total samples from the joint posterior. 
## MCMC ran for 0.003 minutes at time 2017-06-18 20:05:54.
## 
##            mean    sd   2.5%    50%  97.5% overlap0     f  Rhat n.eff
## b0        1.514 1.945 -2.513  1.645  5.356     TRUE 0.794 1.002  1500
## b1        0.310 0.212 -0.091  0.305  0.740     TRUE 0.931 1.005  1500
## s         3.869 1.350  2.262  3.628  7.244    FALSE 1.000 1.057   159
## deviance 53.652 3.208 50.120 52.781 61.496    FALSE 1.000 1.049   135
## 
## Successful convergence based on Rhat values (all < 1.1). 
## Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## For each parameter, n.eff is a crude measure of effective sample size. 
## 
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
## 
## DIC info: (pD = var(deviance)/2) 
## pD = 5.1 and DIC = 58.727 
## DIC is an estimate of expected predictive error (lower is better).

¿Qué pueden decir de estos resultados?

Podemos graficar los datos y el modelo ajustado (en negro). Y ya que estamos compararlo con el resultado de lm (en gris). También podemos agregar al gráfico la relación verdadera entre \(x\) y \(\mu\) (en rojo)

xvec <- seq(0, 20, length = 100)  # generamos una secuenca de valores de x para hacer el gráfico

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")

plot(x, y, xlim = c(0, 20), pch = 21, bg = "gray")
lines(xvec, ml.sim$mean$b0 + ml.sim$mean$b1 * xvec, lwd = 3)  #jags
lines(xvec, b0 + b1 * xvec, lwd = 3, col = 2)  # la verdad
abline(lm(y ~ x), lty = 2, col = "gray", lwd = 3)  # lm

par(op)

Otra forma de ver qué aprendimos acerca de la relación entre \(x\) e \(y\) es ver la posterior de la pendiente. En este caso, como conocemos el valor verdadero de la pendiente, lo podemos incluir en el gráfico a ver cuánto difiere de lo que dice el análisis.

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(ml.sim$sims.list$b1), main = "", xlab = expression("b"[1]), lwd = 2)  #Posterior para b1
abline(v = b1, lwd = 2, lty = 2)  #Valor verdadero de b1.

Ejercicios

1. Mostrar gráficamente la posterior para la ordenada al origen y su valor verdadero.

2. ¿Qué pasa si para la previa de `b1 usamos una precisión de \(10\)? ¿Y si usamos una de \(0.001\)?

3. Repetir el análisis pero esta vez simulando \(100\) observaciones. ¿Qué puede decir respecto a las posteriores de la ordenada al origen y la pendiente?

4. Hacer simulaciones y análisis como los de arriba pero reemplazando la distribución Normal por una Poisson. Para eso tienen que seguir los mismo pasos que hicimos con la distribución Normal. Presten atención al número de parámetros que tiene la distribución de Poison y qué valores pueden tomar.

---

sessionInfo()

## R version 3.4.0 (2017-04-21)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.5
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] jagsUI_1.4.4    lattice_0.20-35 knitr_1.16     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.11    digest_0.6.12   rprojroot_1.2   grid_3.4.0     
##  [5] backports_1.1.0 formatR_1.5     magrittr_1.5    coda_0.19-1    
##  [9] evaluate_0.10   stringi_1.1.5   rjags_4-6       rmarkdown_1.5  
## [13] tools_3.4.0     stringr_1.2.0   parallel_3.4.0  yaml_2.1.14    
## [17] compiler_3.4.0  htmltools_0.3.6

---

Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2017-06-18