Objetivos:

Ajustar regresiones a datos “del mundo real”.
Escribir funciones de likelihood (modelos de datos) y asignar previas a parámetros.

Polinización de Pomelos

Para este ejemplo vamos a usar datos de Natacha Chacoff que observó visitas de polinizadores a flores de pomelo y luego contó los granos de polen depositados. Su objetivo era ver cómo el número de visitas a las flores se relacionaba con la deposición de polen y eventualmente con la producción de frutos.

Veamos los datos, los cuales podemos cargar desde un repositorio online (Opción 1) o desde un archivo guardado en nuestro directorio de trabajo (Opción 2):

# Opción 1:
apis <- read.table("https://sites.google.com/site/modelosydatos/polen_Apis.csv", 
    header = T, sep = ",")

# Opción 2: apis <- read.table('polen_Apis.csv', header = T, sep = ',')

Para ver qué tenemos dentro del objeto apis que creamos al leer los datos podemos usar la función str() que nos muestra la “estructura” del objeto:

str(apis)

## 'data.frame':    190 obs. of  2 variables:
##  $ visitas: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ granos : int  88 34 0 0 0 0 47 2 0 5 ...

Vemos que se trata de un “dataframe” con dos variables; visitas y granos. Las dos variables son números enteros (int). Estas variables corresponden a datos de una variable predictora (visitas) y a una variable de respuesta (granos). Podemos ver gráficamente cómo se relacionan:

op <- par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(apis$visitas, apis$granos, xlab = "Número de visitas", ylab = "Granos de polen depositados", 
    pch = 21, bg = "gray", cex = 1.3)

par(op)

Para modelar esta relación entre número de visitas de polinizadores y polen depositado en la flor, tenemos que definir una distribución para la variable de respuesta (granos de polen) y conectar el valor esperado de esa distribución con la covariable (número de visitas).

Asumamos que la relación entre estas variables es lineal:

cat(file = "ml.bug", "
      model  {
           # Función de likelihood
           for( i in 1:n){
               granos[i] ~ dpois(lambda[i])
               lambda[i] <- b0 + b1 * visitas[i]
           }
           
           # Previas
           b0 ~ dnorm(0,0.01)
           b1 ~ dnorm(0,0.01)
    }")

Ahora vamos a definir qué cosas le vamos a pasar a JAGS, qué cosas queremos que guarde, y creamos una función para generar valores iniciales de las cadenas Markovianas. También 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:

granos <- apis$granos
visitas <- apis$visitas
n <- length(granos)

data <- list("granos", "visitas", "n")

params <- c("b0", "b1")

inits <- function() list(b0 = rnorm(1, 20, 1), b1 = rnorm(1, 10, 0.5))

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

Ahora llamamos a JAGS usando la función jags del paquete jagsUI:

library(jagsUI)
ml.sim <- jags(data = data, 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: 190
##    Unobserved stochastic nodes: 2
##    Total graph size: 407
## 
## 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.

Como siempre, chequeamos si convergieron las cadenas y si el n.eff es suficiente.

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.007 minutes at time 2019-03-02 12:32:00.
## 
##              mean    sd     2.5%      50%    97.5% overlap0 f  Rhat n.eff
## b0         14.864 0.361   14.189   14.860   15.565    FALSE 1 1.001  1500
## b1          6.861 0.252    6.365    6.871    7.343    FALSE 1 1.002  1500
## deviance 5418.401 1.853 5416.546 5417.875 5423.259    FALSE 1 1.004   468
## 
## 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 = 1.7 and DIC = 5420.113 
## DIC is an estimate of expected predictive error (lower is better).

Como todas las posteriores tienen Rhat \(\leq 1.1\) podemos ver con confianza qué pasó con los parámetros que nos interesan. Vemos que la pendiente de la recta (b1) es positiva y claramente diferente de cero. Esto lo podemos comprobar mirando los cuantiles de la muestra de la posterior (el de \(2.5\) es de 6.36 y el de \(97.5\)% de 7.34), la columna de overlap0 y la columna f.

Para ver gráficamente cómo es el modelo que ajustamos podemos hacer un plot de los datos y de la recta estimada:

op <- par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(apis$visitas, apis$granos, xlab = "Número de visitas", ylab = "Granos de polen depositados", 
    pch = 21, bg = "gray")
xvisitas <- seq(min(visitas), max(visitas), 1)
lines(xvisitas, ml.sim$mean$b0 + ml.sim$mean$b1 * xvisitas, lwd = 3)

par(op)

¿Cómo podemos ver la incertidumbre alrededor de esta estimación? Para esto, vamos a hacer uso de la muestra de la posterior conjunta que obtuvimos con JAGS. Podemos hacer algo un poco “a lo bestia” graficando las distintas rectas que surgen de las combinaciones de ordenada al origen y pendiente que aparecen el la posterior:

op <- par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(apis$visitas, apis$granos, xlab = "Número de visitas", ylab = "Granos de polen depositados", 
    pch = 21, bg = "gray")
xvisitas <- seq(min(visitas), max(visitas), 1)
for (i in 1:ml.sim$mcmc.info$n.samples) curve(ml.sim$sims.list$b0[i] + ml.sim$sims.list$b1[i] * 
    x, add = TRUE, col = rgb(0.5, 0.5, 0.5, alpha = 0.1))

par(op)

Vemos que no hay demasiaba variabilidad respecto a los valores esperados de granos de polen depositado a medida que aumenta el número de visitas de polinizadores y que la incertidumbre es mayor para valores grandes de número de visitas.

Pero como conocemos nuestro “modelo de datos” podemos usar las muestras de la posterior conjunta para ver algo más interesante; los intervalos de credibilidad para la variable de respuesta (granos de polen).

nsample <- ml.sim$mcmc.info$n.samples  # tamaño de la muestra de la posterior
niter <- 1000  # cantidad de predicciones que vamos a simular
pred <- matrix(NA, niter, length(xvisitas))  # matriz donde voy a guardar las predicciones

for (i in 1:niter) {
    idx <- sample(1:nsample, 1)  # idx es una variable temporal que toma un valor de 1 a nsample 
    pred[i, ] <- rpois(length(xvisitas), lambda = ml.sim$sims.list$b0[idx] + 
        ml.sim$sims.list$b1[idx] * xvisitas)  # predice datos con los valores estimados de los parámetros
}

library(coda)
hpd <- HPDinterval(as.mcmc(pred), prob = 0.95)  # calcula el intervalo de credibilidad con los datos predichos.

# Graficamos el ajuste del modelo a los datos:
op <- par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(apis$visitas, apis$granos, xlab = "Número de visitas", ylab = "Granos de polen depositados", 
    pch = 21, bg = "gray", cex = 1.3)
lines(xvisitas, exp(ml.sim$mean$b0 + ml.sim$mean$b1 * xvisitas), lwd = 2)
lines(xvisitas, hpd[, 1], lwd = 2, lty = 2)
lines(xvisitas, hpd[, 2], lwd = 2, lty = 2)

par(op)

Preguntas

¿Qué les parece la comparación entre los datos y los intervalos de credibilidad?
¿Qué podemos hacer al respecto?

Ejercicios:

En el “modelo de datos” que usamos recién, asumimos una función lineal para el valor esperado de la Poisson. Es decir:

\[ y_i \sim \text{Pois} (\lambda_i) \\ \lambda_i = b_0 + b_1 \cdot x_i \\ \]

Si bien en este caso no tuvimos inconvenientes en ajustar nuestro modelo a los datos, tenemos que tener en cuenta que \(\lambda\) no puede tomar valores negativos y no hay nada en nuestro modelo de datos que impida valores de \(\lambda < 0\). La forma convencional de resolver eso es usar un “log link” y hacer que el modelo lineal sea sobre \(\log(\lambda)\). Es decir

\[ y_i \sim \text{Pois} (\lambda_i) \\ \log(\lambda_i) = b_0 + b_1 \cdot x_i \\ \]

1.1 Modifique el modelo de BUGS para que use un log link, ajusten el modelo usando JAGS y grafiquen los resultados.¿Qué diferencias encuentran?

1.2 Ahora vamos a considerar la posibilidad de que la cantidad de polen depositado no pueda aumentar indefinidamente a medida que aumentan las visitas de polinizadores. Es decir, suponer que la relación en vez de ser lineal es de tipo monomolecular. Busquen en el bestiario de funciones cómo es la forma de esta función. Es importante familiarizarse con distintas funciones ya que en general las relaciones en la naturaleza no son lineales. Escriban la función de likelihood en papel. Si se acostumbran a escribir las funciones de likelihood en papel después es más fácil escribir el modelo en R y en BUGS. Tengan en cuenta que cierto grado de autopolinización es posible, de manera que en ausencia de visitas puede haber polen depositado. ¿Cómo podemos tener esto en cuenta si la función monomolecular comienza en cero?

1.3 Ajusten el modelo con JAGS y hagan un gráfico de los datos con los dos modelos propuestos. Más adelante vamos a ver métodos formales (y no tanto) para comprar la utilidad de modelos alternativos para un mismo set de datos.

Extra! Modifiquen alguno de estos “modelos de datos” para que la distribución de granos de polen sea una Binomial negativa.

Ejercicio 2: Establecimiento de Alisos en el Noroeste de Argentina (NOA)

Vamos a modelar el establecimiento de Alisos (Alnus acuminata) en el NOA en función de la variabilidad en la lluvia. Los datos son de Ezequiel Araoz. Carguemos los datos y grafiquemos cómo cambia el número de establecimientos de Aliso (variable respuesta) en función de la variabilidad en la precipitación (covariable) en distintas cuencas.

# Opción 1:
datos <- read.table("https://sites.google.com/site/modelosydatos/alisos.txt", 
    header = TRUE)

# Opción 2: datos <- read.table('alisos.txt', header = T)

str(datos)

## 'data.frame':    567 obs. of  4 variables:
##  $ yr    : int  81 80 79 78 77 76 75 74 73 72 ...
##  $ lluvia: num  0.32 1.633 0.856 -0.694 -1.402 ...
##  $ ne    : int  0 2 2 2 2 0 0 2 1 0 ...
##  $ cuenca: int  5 5 5 5 5 5 5 5 5 5 ...

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(81 - datos$yr[1:81], datos$lluvia[1:81], type = "l", lwd = 2, xlab = "Año", 
    ylab = "Lluvia (desvíos)")

par(op)

ncuencas <- length(unique(datos$cuenca))

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(datos$lluvia, datos$ne, xlab = "lluvia (desvíos)", ylab = "Número de Establecimientos")

for (i in 1:ncuencas) {
    points(datos$lluvia[datos$cuenca == i], datos$ne[datos$cuenca == i], pch = 21, 
        bg = i + 1)
    
}

par(op)

Los colores de puntos corresponden a las distintas cuencas. Vamos a asumir que no hay diferencias entre cuencas.

2.1 Primero que nada piensen qué tipo de distribución usarían para estos datos.

2.2 Luego, ajusten dos modelos diferentes, asumiendo:

que la relación entre el número de establecimientos y la variabilidad en la precipitación es lineal.
que la relación entre el número de establecimientos y la variabilidad en la precipitación tiene un óptimo (Busquen otra vez ideas en el bestiario de funciones).

2.3 Finalmente, comparar gráficamente el ajuste de ambos modelos a los datos.

Opcional: Regresiones con `Stan` y `brms`

Veamos cómo haríamos el análisis de los datos de pomelos con Stan.

cat(file = "poiss_reg.stan", "
data { 
  int<lower=1> n;
  vector[n] visitas;
  int<lower=0> granos[n];
}

parameters {
  real b0;
  real b1;
}

model{
  b0 ~ normal(0,10);
  b1 ~ normal(0,10);
  
  granos ~ poisson_log(b0 + b1 * visitas);
}
")

Ahora llamamos a Stan

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

reg_dat <- list(n = n, granos = granos, visitas = visitas)

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

Vemos los resultados

print(fit)

## Inference for Stan model: poiss_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
## b0      2.83    0.00 0.02    2.80    2.82    2.83    2.85    2.87   339
## b1      0.19    0.00 0.01    0.18    0.19    0.19    0.20    0.20   427
## lp__ 9634.47    0.05 1.02 9631.64 9634.08 9634.78 9635.19 9635.45   501
##      Rhat
## b0   1.00
## b1   1.00
## lp__ 1.01
## 
## Samples were drawn using NUTS(diag_e) at Sat Mar  2 12:32:33 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).

Para ver los resultados gráficamente

pos <- as.data.frame(fit)
op <- par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(apis$visitas, apis$granos, xlab = "Número de visitas", ylab = "Granos de polen depositados", 
    pch = 21, bg = "gray")
xvisitas <- seq(min(visitas), max(visitas), 1)
lines(xvisitas, exp(mean(pos$b0) + mean(pos$b1) * xvisitas), lwd = 3)

par(op)

Si queremos ver la incertidumbre en la variable de respuesta

nsample <- dim(pos)[1]  # tamaño de la muestra de la posterior
niter <- 1000  # cantidad de predicciones que vamos a simular
pred <- matrix(NA, niter, length(xvisitas))  # matriz donde voy a guardar las predicciones

for (i in 1:niter) {
    idx <- sample(1:nsample, 1)  # idx es una variable temporal que toma un valor de 1 a nsample 
    pred[i, ] <- rpois(length(xvisitas), lambda = exp(pos$b0[idx] + pos$b1[idx] * 
        xvisitas))  # predice datos con los valores estimados de los parámetros
}

library(coda)
hpd <- HPDinterval(as.mcmc(pred), prob = 0.95)  # calcula el intervalo de credibilidad con los datos predichos.

# Graficamos el ajuste del modelo a los datos:
op <- par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(apis$visitas, apis$granos, xlab = "Número de visitas", ylab = "Granos de polen depositados", 
    pch = 21, bg = "gray", cex = 1.3)
lines(xvisitas, exp(ml.sim$mean$b0 + ml.sim$mean$b1 * xvisitas), lwd = 2)
lines(xvisitas, hpd[, 1], lwd = 2, lty = 2)
lines(xvisitas, hpd[, 2], lwd = 2, lty = 2)

par(op)

Para ajustar la versión monomolecular tenemos que definir lambda y usar la distribución de Poisson. También tenemos que usar una previa razonable para el parámetro de valores asintóticos de la monomolecular. Un aspecto importante en este ejemplo es que tenemos granos de polen observados para casos en los que hubo \(0\) visitas. Entonces tenemos que usar una versión de la monomolecular que permita un “ofset”. Por eso agregamos el parámetro “c” para permitir que lambda sea mayor que cero cuando hay cero visitas

cat(file = "poiss_mreg.stan", "
data { 
  int<lower=1> n;
  vector[n] visitas;
  int<lower=0> granos[n];
}

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

model{
  vector[n] lambda;
  b0 ~ normal(100,10);
  b1 ~ normal(0,10);
  lambda = b0 * (1 - exp(-b1 * visitas)) + c;
  granos ~ poisson(lambda);
}
")

Por defecto, Stan usa valores iniciales para los parámetros tomando números al azar entre \(-2\) y \(2\). Pero para la asíntota, necesitamos valores iniciales que tengan sentido. Para eso vamos a escribir una función que genere valores iniciales.

init_fun <- function() {
    list(b0 = runif(1, 10, 130), b1 = runif(1, 0, 1), c = runif(1, 0, 10))
}

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

Ejercicio:

Graficar los intervalos de credibilidad para la variable de respuesta (granos de polen) de este modelo.

Para ajustar el modelo linel con log-link usando brms hacemos

library(brms)

priors <- c(prior(normal(0, 10), class = "b"))

fit_brms <- brm(granos ~ visitas, data = data.frame(granos = granos, visitas = visitas), 
    family = poisson(), prior = priors)

Para ajustar del modelo monomolecular, tenemos que avisarle a brms que ajustamos un modelo no-lineal y definir previas para los parámetros. Además, como usamos una Poisson, tenemos que definir log(lambda) en lugar de lambda…

priors <- prior(normal(100, 10), nlpar = "b0") + prior(normal(0, 2), nlpar = "b1") + 
    prior(normal(0, 2), nlpar = "c")

fit_brms <- brm(bf(granos ~ log(b0 * (1 - exp(-b1 * visitas)) + c), b0 + b1 + 
    c ~ 1, nl = TRUE), data = data.frame(granos = granos, visitas = visitas), 
    family = poisson(), prior = priors)

Ejercicio:

Ajustar el model con distribución binomial negativa

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

Objetivos:

Polinización de Pomelos

Preguntas

Ejercicios:

Ejercicio 2: Establecimiento de Alisos en el Noroeste de Argentina (NOA)

Opcional: Regresiones con Stan y brms

Ejercicio:

Ejercicio:

Opcional: Regresiones con `Stan` y `brms`