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 (modelos de datos)
• 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(1234)
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). En el modelo que vamos a definir, ponemos como previa para el desvío una uniforme entre \(0\) y \(100\) y luego transformamos a valores de 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/(sigma*sigma) 
           sigma ~ 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. En general, es importante centrar y estandarizar las variables predictoras de una regresión. Esto hace que las pendientes queden en unidades de desvió estándar de las predictoras (covariables) y que la ordenada al origen nos de el valor esperado de la variable de respuesta (y) cuando la predictora (x) está en su valor promedio. Para transformar (escalar) la variabe predictora, usamos la función scale

library(jagsUI)

# defino las variables
x <- as.numeric(scale(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", "sigma")

# 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), sigma = 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..... 
## Adaptive phase complete 
##  
## 
##  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 revisar 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\)) pero en general, valores del orden de \(1000\) son suficientes.

print(ml.sim)

## JAGS output for model 'ml.bug', generated by jagsUI.
## Estimates based on 3 chains of 1000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 500 iterations and thin rate = 1,
## yielding 1500 total samples from the joint posterior. 
## MCMC ran for 0.001 minutes at time 2019-02-28 14:40:04.
## 
##            mean    sd   2.5%    50%  97.5% overlap0     f  Rhat n.eff
## b0        2.804 0.755  1.274  2.808  4.323    FALSE 1.000 1.001  1240
## b1        0.913 0.821 -0.753  0.922  2.619     TRUE 0.886 1.000  1500
## sigma     2.411 0.708  1.451  2.258  4.118    FALSE 1.000 1.000  1500
## deviance 44.122 3.037 40.575 43.367 51.804    FALSE 1.000 1.000  1500
## 
## 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 = 4.6 and DIC = 48.738 
## DIC is an estimate of expected predictive error (lower is better).

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). Noten que tenemos que transformar los valores de los coeficientes porque habíamos centrado y estandarizado la variable predictora.

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

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

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

Figure 1: Comparación entre la relación verdadera (en rojo), la estimación de JAGS (en negro), y la estimación del modelo lineal (en gris)

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. De esta manera, podemos comprobar que nuestro procedimiento de generar muestras de la posterior a partir de cadenas Markovianas hace una buena estimación de la pendiente.

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
hist(ml.sim$sims.list$b1/sd(x.sim), 60, freq = FALSE, main = "", xlab = expression("b"[1]))
abline(v = b1, lwd = 3, lty = 2)  #Valor verdadero de b1.

Figure 2: Histograma de la muestra de la posterior de la pendiente de la regresión (b_1) generada con JAGS. La línea muestra el valor verdadero de la pendiente

par(op)

Por otro lado vemos que con estos datos, tenemos cierta incertidumbre respecto de si la pendiente es poisitiva o negativa porque parte de la posteior es menor que cero. Para ver con más detalle qué porcentaje de la posterior es menor que cero podemos hacer length(which(ml.sim$sims.list$b1 < 0)) / length(ml.sim$sims.list$b1), que en este caso da 0.114.

¿Cómo compara este resultado con lo que obtendríamos si hacemos un análisis frecuentista? Podemos ajustar fácilmente este modelo usando la función lm de R para obtener:

fit_lm <- lm(y ~ x)
summary(fit_lm)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2115 -1.2566 -0.6228  0.8495  3.8785 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   2.9676     0.6405   4.633  0.00168 **
## x             0.9817     0.6751   1.454  0.18400   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.025 on 8 degrees of freedom
## Multiple R-squared:  0.209,  Adjusted R-squared:  0.1102 
## F-statistic: 2.114 on 1 and 8 DF,  p-value: 0.184

Vemos que los coeficientes que estimamos son parecidos a los que encontramos con el análisis Bayesiano y que la pendiente no es significativamente diferente de \(0\) para un alpha de \(0.05\). Estos resultados no deberían sorprendernos porque las previas que usamos en el análisis Bayesiano no eran muy informativas. De hecho podemos tratar de usar los resultados de el modelo lineal para aproximar una medida de incertidumbre alrededor de la pendiente y compararla con la posterior. Para hacer esto, vamos a necesitar estimar el error estándar de la pendiente y suponer que podemos representar la incertidumbre con una distribución normal.

se <- sqrt(diag(vcov(fit_lm)))

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
hist(ml.sim$sims.list$b1/sd(x.sim), 60, freq = FALSE, main = "", xlab = expression("b"[1]))
curve(dnorm(x, fit_lm$coefficients[2]/sd(x.sim), se[2]/sd(x.sim)), add = TRUE, 
    lwd = 2, col = 2)
abline(v = b1, lwd = 3, lty = 2)  #Valor verdadero de b1.

Figure 3: Posterior estimada con MCMC (en negro) y estiamción de incertidumbre asumiendo distribución normal y usando el error estándar de lo estimado con un modelo lineal (en rojo)

par(op)

Ejercicios

1. ¿De qué otra forma podemos estimar la fracción de la posterior que es menor que cero?

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

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

4. 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?

5. 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.

---

Opcional: Regresiones con Stan y brms

Veamos cómo haríamos la regesión en Stan. Como vimos antes, tenemos que especificar varios detalles. Además, Stan usa desvío estándar en lugar de precisión. En este caso vamos a usar una distribución t para la previa del desvío de la normal. Un aspecto relevante es que como definimos los límites de las variables (por ejemplo con real<lower=0> sigma;), no necesitamos truncar distribuciones. De esta forma, como decimos que el desvío es real<lower=0> la previa va a ser automáticamente truncada a valores no negativos.

cat(file = "reg.stan", "
data { 
  int<lower=1> n;
  vector[n] x;
  vector[n] y;
}

parameters {
  real b0;
  real b1;
  real<lower=0> sigma;
}

model{
  vector[n] mu;
  b0 ~ normal(0,10);
  b1 ~ normal(0,10);
  sigma ~ student_t(3, 0, 10);
  
  mu = b0 + b1 * x;
  y ~ normal(mu, sigma);
}
")

Ahora armamos una lista con los datos para Stan y generamos las cadenas Markovianas.

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

reg_dat <- list(n = n, y = y.sim, x = as.numeric(scale(x.sim)))

fit <- stan(file = "reg.stan", data = reg_dat, iter = 1000, thin = 1, chains = 3)

Para ver las estimaciones de parámetros hacemos

print(fit)

## Inference for Stan model: reg.
## 3 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=1500.
## 
##         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## b0      2.93    0.02 0.76   1.45   2.46   2.94   3.38   4.43  1105 1.00
## b1      1.03    0.03 0.86  -0.72   0.53   0.99   1.53   2.83   881 1.00
## sigma   2.39    0.03 0.74   1.42   1.89   2.24   2.71   4.09   645 1.00
## lp__  -12.17    0.07 1.46 -15.74 -12.90 -11.78 -11.09 -10.53   465 1.01
## 
## Samples were drawn using NUTS(diag_e) at Thu Feb 28 14:40:08 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Las estimaciones son un poco diferentes a las de JAGS, veamos cómo comparan en la posterior de la pendiente

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

Para ajustar este modelo usando brms hacemos

library(brms)

priors <- c(prior(normal(0, 10), class = "b"), prior(student_t(3, 0, 10), class = "sigma"))
fit_brms <- brm(y ~ x, data = data.frame(x = scale(x.sim), y), family = gaussian(), 
    prior = priors)

Para ver las estimaciones de parámetros hacemos

print(fit_brms)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ x 
##    Data: data.frame(x = scale(x.sim), y) (Number of observations: 10) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept     2.94      0.80     1.33     4.55       2394 1.00
## x             0.97      0.85    -0.77     2.64       2691 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     2.40      0.76     1.43     4.30       1670 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Podemos graficar fácilmente los resutlados

plot(fit_brms)

Veamos ahora como comparan las tres estimaciones de la pendiente.

sims <- posterior_samples(fit_brms)

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(ml.sim$sims.list$b1/sd(x.sim)), main = "", xlab = expression("b"[1]), 
    lwd = 2)  # JAGS
lines(density(b1_sam$b1/sd(x.sim)), lwd = 2, col = 2)  # Stan
lines(density(sims$b_x/sd(x.sim)), lwd = 2, col = 3)  # brms
abline(v = b1, lwd = 2, lty = 2)  #Valor verdadero de b1.

Figure 4: Tres aproximaciones a la posterior de la pendiente de una regresión lineal. JAGS (en negro), Stan (rojo) y Stan via brms (verde).

---

Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2019-02-28