Introduccion modelos jerárquicos (I)
1 Introducción
Los modelos jerárquicos (mixtos o multinivel), corresponden a un grupo de modelos estadísticos que permiten incorporar explícitamente la jerarquía de la variación en los datos. Esto es, estiman la distribución posterior de los parámetros del modelo a través agrupamientos parciales (partial pooling) de los grupos, donde los parámetros particulares de cada grupo se estima de una distribución de probabilidad de la “poblacional” de grupos — i.e. una distribución de probabilidad de coeficientes. Dado que todos los grupos provienen de una misma distribución de probabilidad, las observaciones de un grupo ayudan a estimar los coeficientes de otro. Lo que resulta ventajoso cuando el esfuerzo de muestreo no está balanceado entre grupos. Veamos un ejemplo, aclaremos la idea de agrupamientos en los datos y sus implicaciones para la estimación de los coeficientes en un modelo jerárquico.
Modelemos la probabilidad de supervivencia de huevos en 50 nidadas de una ave cualquiera en un parche de bosque. Tenemos tres alternativas:
Modelar los datos asumiendo que todos los nidos tienen igual probabilidad \(p\) de ser depredados (pooling), y por lo tanto estimando un único parámetro común para todos. Esto es, asumimos que la variación entre nidos es cero.
Estimar de manera independiente el coeficiente \(p\) para cada nido (no pooling).
Estimar \(p\) para cada nido, pero asumiendo que cada \(p_i\) hace parte de una distribución de probabilidad \(Binomial(p)\) de la población de nidos (partial pooling).
Simulemos los datos
Ajustemos el modelo pooling.
cat(file = 'multilevel_mods/pooling.stan',
'
data{
int N;
array[N] int surv;
array[N] int huevos;
}
parameters{
real alpha;
}
model{
real p;
alpha ~ normal(5, 2.5);
p = alpha;
p = inv_logit(p);
surv ~ binomial(huevos, p);
}
generated quantities{
vector[N] log_lik;
real p;
p = alpha;
p = inv_logit(p);
for (i in 1:N) log_lik[i] = binomial_lpmf(surv[i] | huevos[i], p);
}
'
)file <- paste(getwd(), '/pooling.stan', sep = '')
fit_pooling <- cmdstan_model(file, compile = T)
pooling <-
fit_pooling$sample(
data = dat,
iter_sampling = 2e3,
iter_warmup = 500,
chains = 3,
parallel_chains = 3,
thin = 3,
refresh = 500,
seed = 123
)
trace_plot <- function(fit, par, n_chains) {
d <- as_mcmc.list(fit)
d <- lapply(d, FUN =
function(x) {
x <- as.data.frame(unclass(as.matrix(x)))
x$iter <- 1:nrow(x)
x
})
for (i in 1:n_chains) d[[i]]$chain <- i
plot(d[[1]][, 'iter', drop = T], d[[1]][, par, drop = T],
type = 'l', col = 1, main = '', ylab = par,
xlab = 'Iteration')
for (i in 2:n_chains) {
lines(d[[i]][, 'iter', drop = T], d[[i]][, par, drop = T],
col = i, main = '', ylab = par,
xlab = 'Iteration')
}
}
mod_diagnostics <- function(model, output) {
l <- model$loo()
pk <- l$diagnostics$pareto_k
pareto <- tibble(x = 1:length(pk),
pareto_k = pk)
diags <- na.omit(tibble(x = 1:nrow(output),
rhat = output$rhat,
ess_bulk = output$ess_bulk,
ess_tail = output$ess_tail))
par(mfrow = c(1, 3), mar = c(4, 4, 1, 1))
diags %$% plot(rhat ~ x, pch = 16, col = 'tomato3', xlab = 'Parameters',
ylab = 'Rhat')
abline(h = 1.1, lty = 3, col = 'red')
diags %$% plot(sort(ess_bulk) ~ x, pch = 16, col = 'cyan4', xlab = 'Parameters',
ylab = 'ESS')
diags %$% points(sort(ess_tail) ~ x, pch = 16, col = 'tan1')
abline(h = 1000, lty = 3, col = 'red')
pareto %$% plot(pareto_k ~ x, pch = 16, col = 'purple', xlab = 'Observations',
ylab = 'Pareto-k', ylim = c(-0.5, 1.2))
abline(h = c(0.5, 0.7, 1), lty = 3, col = 'red')
text(x = rep(nrow(pareto)/2, 4), y = c(0.25, 0.6, 0.8, 1.2),
labels = c('Perfect', 'ok', 'high', 'too hight'))
par(mfrow = c(1, 1))
}
(out_pooling <- pooling$summary())
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1))
for (i in 1:2) trace_plot(pooling, out_pooling$variable[i], 3)
Modelo no pooling
cat(file = 'multilevel_mods/no_pooling.stan',
'
data{
int N;
array[N] int surv;
array[N] int huevos;
array[N] int nido;
}
parameters{
vector[N] alpha;
}
model {
vector[N] p;
alpha ~ normal(5, 25); // previa para cada nido
for (i in 1:N) {
p[i] = alpha[nido[i]]; // estimación estratificada por nido
p[i] = inv_logit(p[i]);
}
surv ~ binomial(huevos, p);
}
generated quantities{
vector[N] log_lik;
vector[N] p;
for (i in 1:N) {
p[i] = alpha[nido[i]];
p[i] = inv_logit(p[i]);
}
for (i in 1:N) log_lik[i] = binomial_lpmf(surv[i] | huevos[i], p[i]);
}
')file <- paste(getwd(), '/no_pooling.stan', sep = '')
fit_no_pooling <- cmdstan_model(file, compile = T)
no_pooling <-
fit_no_pooling$sample(
data = dat,
iter_sampling = 2e3,
iter_warmup = 500,
chains = 3,
parallel_chains = 3,
thin = 3,
refresh = 500,
seed = 123
)
(out_no_pooling <- no_pooling$summary())
par(mfrow = c(3, 3), mar = c(4, 4, 1, 1))
for (i in 1:9) trace_plot(no_pooling, out_no_pooling$variable[i], 3)
Modelo partial pooling
cat(file = 'multilevel_mods/par_pooling.stan',
'
data{
int N;
array[N] int surv;
array[N] int huevos;
array[N] int nido;
}
parameters{
vector[N] alpha;
real mu;
real<lower = 0> sigma;
}
model{
vector[N] p;
alpha ~ normal(mu, sigma); // previa para la población de nidos (previa adaptativa)
mu ~ normal(5, 2.5); // Hiper previa (previa de la previa)
sigma ~ exponential(1); // Hiper previa (previa de la previa)
for (i in 1:N) {
p[i] = alpha[nido[i]]; // modelo estratificado por nido
p[i] = inv_logit(p[i]);
}
surv ~ binomial(huevos, p);
}
generated quantities{
vector[N] log_lik;
vector[N] p;
for (i in 1:N) {
p[i] = alpha[nido[i]];
p[i] = inv_logit(p[i]);
}
for (i in 1:N) log_lik[i] = binomial_lpmf(surv[i] | huevos[i], p[i]);
}
')file <- paste(getwd(), '/par_pooling.stan', sep = '')
fit_par_pooling <- cmdstan_model(file, compile = T)
par_pool <-
fit_par_pooling$sample(
data = dat,
iter_sampling = 2e3,
iter_warmup = 500,
chains = 3,
parallel_chains = 3,
thin = 3,
refresh = 500,
seed = 123
)
(out_par_pool <- par_pool$summary())
par(mfrow = c(3, 3), mar = c(4, 4, 1, 1))
for (i in 1:9) trace_plot(par_pool, out_par_pool$variable[i], 3)
El modelo partial pooling recibe el nombre de modelo jerárquico, porque la estimación de los parámetros se realiza en varios niveles. Veamos en notación matemática el modelo anterior:
\[ \begin{align} &Sup.~huevos_i \sim Binomial(huevos_i,~p_i) & ,~Funcion~de~likelihhod \\ &logit(p_i) \sim \alpha[nido_i] &,~Modelo~lineal \\ &\alpha[nido] \sim Normal(\mu_j, \sigma_j) &,~Previa~para~la~población~de~nidos \\ &\mu_j \sim Normal(5, 2.5) &,~Previa~para~cada~nido~(hiperparámetro) \\ &\sigma_j \sim exponential(1) &,~Previa~para~cada~nido~(hiperparámetro) \\ \end{align} \]
\(\mu_j\) y \(\sigma_j\) se conocen como hiperparámetros porque son parámetros de parámetros (la previa poblacional), por la misma razón \(Normal(5, 2.5)\) y \(exponential(1)\) se les llama hiper-previas.
Usemos el looic para comparar los tres modelos en términos de su capacidad predictiva fuera de la muestra:
elpd_diff se_diff
model3 0.0 0.0
model1 -11.9 5.5
model2 -15.4 5.1 El looic del modelo jerárquico es notablemente más bajo comparado con los modelos pooling y no pooling (i.e. tiene un mayor poder predictivo). Veamos qué tan precisa es la estimación del log-odds de los modelos jerárquico y pooling:
log_odds_P <- no_pooling$draws(format = 'df')
log_odds_P <- log_odds_P[, grepl('alpha', colnames(log_odds_P))]
log_odds_parP <- par_pool$draws(format = 'df')
plot(NULL, xlim = c(-5, 5), ylim = c(0, 1.2), xlab = 'log-odds nidos',
ylab = 'Densidad')
for (i in 1:500) {
curve(dnorm(x, log_odds_parP$mu[i], log_odds_parP$sigma[i]),
lwd = 0.1, add = T, col = "red")
lines(density(as.matrix(log_odds_P[i, ])), lwd = 0.1, col = 'tan1')
}
curve(dnorm(x, mean(log_odds_parP$mu), mean(log_odds_parP$sigma)), add = T,
col = 'red', lwd = 2.5)
lines(density(apply(log_odds_P, 2, mean)), col = 'tan1', lwd = 2.5)
lines(density(logit(rbeta(1e3, 3, 5))), col = 'black', lwd = 2.5)
Ambos modelos predicen bien el promedio del log-odd real (negro). Sin embargo, notemos la amplitud de las distribuciones posteriores, la estimación del modelo jerárquico (rojo) es mucho más precisa comparada con la del modelo pooling (naranja).
Ahora comparemos el error de la probabilidad de supervivencia estimada por ambos modelos, con la probabilidad “real” — la simulada a partir de \(Beta(\phi_1 = 3, \phi_2= 5)\).
post_parPool <- par_pool$draws(format = 'df')
post_parPool <- as_tibble(post_parPool[, grepl('alpha', colnames(post_parPool))])
colnames(post_parPool) <- paste('nido', 1:n, sep = '')
post_parPool <- do.call('rbind', apply(post_parPool, 2, simplify = 'list', FUN =
function(x) {
tibble(li = quantile(x, 0.025),
ls = quantile(x, 0.975),
mu = mean(x),
mod = 'Par pooling')
}))
post_parPool$x <- paste('nido', 1:n, sep = '')
post_noPool <- no_pooling$draws(format = 'df')
post_noPool <- as_tibble(post_noPool[, grepl('alpha', colnames(post_noPool))])
colnames(post_noPool) <- paste('nido', 1:n, sep = '')
post_noPool <- do.call('rbind', apply(post_noPool, 2, simplify = 'list', FUN =
function(x) {
tibble(li = quantile(x, 0.025),
ls = quantile(x, 0.975),
mu = mean(x),
mod = 'No pooling')
}))
post_noPool$x <- paste('nido', 1:n, sep = '')
mods <- rbind(post_noPool, post_parPool)
for (i in 1:3) mods[[i]] <- inv_logit(mods[[i]])ggplot() +
geom_errorbar(data = mods, aes(x = x, ymin = li, ymax = ls, color = mod),
position = position_dodge(width = 0.7), width = 0) +
geom_point(data = mods, aes(x = x, y = mu, color = mod),
position = position_dodge(width = 0.7)) +
geom_point(data = tibble(y = p,
x = paste('nido', 1:n, sep = '')),
aes(x, y)) +
geom_hline(yintercept = mean(p), linetype = 2) +
labs(y = 'Probabilidad', x = NULL) +
coord_flip() +
theme_bw() +
theme(panel.grid = element_blank(),
legend.title = element_blank(),
legend.position = 'top')
Los puntos negros corresponden a la probabilidad de supervivencia real. Pareciera que los valores predichos por el modelo jerárquico están, en general, más cercanos a la probabilidad de supervivencia real de los nidos. Calculemos este error de estimación y veamos cómo se comporta con el tamaño de la muestra (i.e. número de huevos observados por nido).
mods$real_p <- rep(p, 2)
mods$sample_size <- c(huevos, huevos)
mods$error <- mods %$% abs(mu - real_p)
ggplot(mods, aes(sample_size, error, color = mod)) +
geom_point() +
geom_smooth(method = 'lm', se = F) +
labs(x = 'Tamaño de la muestra', y = 'Error P(supervivencia)') +
theme_bw() +
theme(panel.grid = element_blank(),
legend.title = element_blank(),
legend.position = 'top')
En ambos modelos el error decrece con el aumento del tamaño de muestra; es natural, más datos por lo general proveen más información y, por lo tanto, una estimación más precisa de los coeficientes. Notemos, sin embargo, que en promedio el modelo jerárquico (par pooling) estima la probabilidad de supervivencia con un menor error, aunque para ambos modelos este valor tiende a converger con el incremento del tamaño de muestra.
2 Modelos con más de un agrupamiento
En el ejemplo anterior modelamos la probabilidad de supervivencia de huevos en 50 nidadas del ave X ubicadas en un parche de bosque cualquiera. Supongamos que la investigadora decidió expandir la investigación, duplicar el número de nidadas analizadas y distribuir su esfuerzo de muestreo en 9 parches de bosques adicionales. Ahora, dado el nuevo diseño, la supervivencia de los huevos no solo varía en función del nido, sino también dado el parche de bosque en el que se encuentra el nido.
Generemos los datos y describamos el modelo en notación matemática:
n_nidos <- 100
n_parches <- 10
set.seed(555)
p <- rbeta(n_nidos, 3, 5)# probabilidad de visita (5 min) sin efecto
set.seed(555)
huevos <- rpois(n_nidos, 10)
set.seed(555)
var_parche <- rnorm(1e3, 0, 0.05)
set.seed(555)
var_parche <- sample(var_parche[var_parche >= 0], n_parches, replace = F)
var_parche <- sapply(var_parche, simplify = 'array', FUN =
function(x) rnorm(n_parches, 0, x))
var_parche <- as.vector(var_parche)
set.seed(555)
huevos_obs <- rbinom(n_nidos,
prob = p + var_parche,
size = huevos)
dat <-
list(
n = huevos,
superv = huevos_obs,
parche = rep(1:10, each = 10),
nido = 1:100,
N = length(huevos),
N_parche = 10
)\
\[ \begin{align} &Sup.~huevos_i \sim Binomial(huevos_i,~p_i) & ,~Funcion~de~likelihhod \\ &logit(p_i) \sim \alpha[nido_i] + \psi[parche] &,~Modelo~lineal \\ &\psi[parche] \sim Normal(0, \tau) &,~Previa variación~entre~parches\\ &\tau \sim exponential(1) &, ~ Hiper-parámetro/previa~de~la~previa~del~parche\\ &\alpha[nido] \sim Normal(\mu_j, \sigma_j) &,~Previa~para~la~población~de~nidos \\ &\mu_j \sim Normal(5, 2.5) &,~Previa~para~cada~nido~(hiperparámetro) \\ &\sigma_j \sim exponential(1) &,~Previa~para~cada~nido~(hiperparámetro) \\ \end{align} \]
El parche corresponde a un factor que altera la supervivencia de los huevos. Es decir, \(p_i\) reduce o aumenta en función del parche. Dichas oscilaciones las incorporamos en el modelo en forma de una distribución \(Normal(0, \tau)\), donde \(\tau\) controla la ‘magnitud’ de los cambios positivos o negativos.
Ajustemos el modelo:
cat(file = 'multilevel_mods/par_pooling2.stan',
'
data{
int N;
int N_parche;
array[N] int n;
array[N] int superv;
array[N] int parche;
array[N] int nido;
}
parameters{
vector[N] alpha;
vector[N_parche] tau;
real muA;
real<lower = 0> sigmaA;
real<lower = 0> sigmaT;
}
model{
vector[N] p;
alpha ~ normal(muA, sigmaA);
tau ~ normal(0, sigmaT);
muA ~ normal(5, 2.5);
sigmaA ~ exponential(1);
sigmaT ~ exponential(1);
for (i in 1:N) {
p[i] = alpha[nido[i]] + tau[parche[i]];
p[i] = inv_logit(p[i]);
}
superv ~ binomial(n, p);
}
generated quantities {
vector[N] log_lik;
vector[N] p;
for (i in 1:N) {
p[i] = alpha[nido[i]] + tau[parche[i]];
p[i] = inv_logit(p[i]);
}
for (i in 1:N) log_lik[i] = binomial_lpmf(superv[i] | n[i], p[i]);
}
')# A tibble: 314 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ -589. -589. 13.1 13.1 -611. -567. 1.00 456.
2 alpha[1] -0.964 -0.969 0.595 0.584 -1.93 -0.00790 1.00 1645.
3 alpha[2] -0.411 -0.427 0.514 0.515 -1.25 0.443 1.00 1352.
4 alpha[3] -0.162 -0.185 0.502 0.489 -0.985 0.684 1.00 1825.
5 alpha[4] 0.878 0.860 0.521 0.526 0.0856 1.77 0.999 1896.
6 alpha[5] -0.454 -0.434 0.633 0.633 -1.49 0.532 1.00 1744.
7 alpha[6] -1.59 -1.57 0.640 0.626 -2.70 -0.583 1.00 1997.
8 alpha[7] -0.655 -0.645 0.480 0.470 -1.44 0.137 1.00 1691.
9 alpha[8] -1.13 -1.11 0.588 0.573 -2.12 -0.234 1.00 1680.
10 alpha[9] -1.20 -1.19 0.553 0.541 -2.15 -0.325 1.00 1775.
# ℹ 304 more rows
# ℹ 1 more variable: ess_tail <dbl>par(mfrow = c(3, 3), mar = c(4, 4, 1, 1))
for (i in 1:9) trace_plot(mod, out_par_pool2$variable[i], 3)
Tenemos un parámetro por cada nido, uno por cada parche y su parámetro de dispersión y la probabilidad de supervivencia promedio. Antes de discutir las implicaciones de los parámetros, grafiquemos las distribuciones predictivas posteriores y estimemos la precisión del modelo para representar las observaciones.
post <- mod$draws(format = 'df')
post <- post[, 2:112]
post <- as.data.frame(apply(post, 2, inv_logit))
ppcheck <- sapply(1:100, simplify = 'array', FUN =
function(x) rbinom(1e3,
size = dat$n,
prob = post$muA[i]))
plot(NULL, xlim = c(-1, 10), ylim = c(0, 0.6))
for (i in 1:100) lines(density(ppcheck[i, ]), lwd = 0.15, col = 'blue')
lines(density(dat$superv), col = 'red', lwd = 3)
abline(v = c(mean(apply(ppcheck, 2, mean)), mean(dat$superv)),
col = c('blue', 'red'), lty = 3)
El modelo se ajusta adecuadamente a los datos, incluso el promedio estimado de huevos supervivientes es básicamente igual al observado.
Veamos las implicaciones de los parámetros
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1))
plot(NULL, xlim = c(0, 1), ylim = c(0, 6), xlab = expression(italic('P')['nido']),
ylab = 'Densidad', main = 'Nido')
for (i in 1:100) lines(density(post[, grep('alpha', colnames(post))][[i]]), lwd = 0.2)
lines(density(apply(post[, grep('alpha', colnames(post))], 2, mean)), lwd = 3, col = 'red')
lines(density(p + var_parche), lwd = 3, col = 'blue')
post2 <- mod$draws(format = 'df')
post2 <- post2[, grepl('sigma', colnames(post2))]
plot(density(post2$sigmaA), main = "Parche", xlim = c(0, 1.5),
xlab = expression(sigma[italic(P)]), ylab = '',
ylim = c(0, 4.5))
lines(density(post2$sigmaT), lty = 2)
Izquierda. Cada línea de densidad negra corresponde a la distribución posterior de \(p\) para cada nido, la línea roja denota el promedio entre nidos —igual que el parámetro \(\mu_i\) en el modelo, y la línea azul corresponde a la probabilidad real (i.e. probabilidad de supervivencia, \(Beta(3, 5)\), más variación dado el parche, \(Normal(0, \tau)\)). Derecha. Distribución posterior de los parámetros de dispersión de los nidos (línea sólida) y los parches (línea discontinua).
La probabilidad de supervivencia promedio es del ~40%.
De hecho, el modelo fue capaz de recuperar la distribución de probabilidad generadora de los datos (línea azul). También, con seguridad, podemos afirmar que la variación en la probabilidad de supervivencia entre nidos es mayor a la producida por los diferentes parches. Podríamos generar un nuevo modelo sin el factor parche, compararlos con looic, y evaluar el aporte del parche en la capacidad predictiva del modelo.
Podríamos concluir (ponele…), que los factores que determinan la probabilidad de supervivencia de huevos en nidos del ave X operan principalmente a escala local —la principal variación proviene del nido, mientras que el aporte del parche es tangencial. Supongamos que la investigadora alcanzó la misma conclusión y, en un intento de identificar dichos factores locales moduladores de la supervivencia, decidió medir la proporción de cobertura vegetal en un radio de 1 m alrededor de cada nido. Veamos:
3 Transiciones divergentes y reparametrización
De nuevo simulemos los datos y definamos el modelo.
n_nidos <- 100
n_parches <- 10
set.seed(555)
p <- rbeta(n_nidos, 3, 5)
set.seed(555)
huevos <- rpois(n_nidos, 10)
set.seed(555)
var_parche <- rnorm(1e3, 0, 0.1)
set.seed(555)
var_parche <- sample(var_parche[var_parche >= 0], n_parches, replace = F)
var_parche <- sapply(var_parche, simplify = 'array', FUN =
function(x) rnorm(n_parches, 0, x))
var_parche <- as.vector(var_parche)
tot_p <- p + var_parche
set.seed(555)
pro_veg <- rbeta(1e3, 2, 1.5)
temp2 <- quantile(pro_veg, probs = seq(0, 1, by = 0.2))
temp <- quantile(tot_p, probs = seq(0, 1, by = 0.2))
veg <- vector('double', 100)
set.seed(555)
for (i in 1:(length(temp)-1)) {
veg[which(tot_p >= temp[i] & tot_p <= temp[i+1])] <-
sample(pro_veg[pro_veg > temp2[i] & pro_veg <= temp2[i+1]],
size = length(which(tot_p >= temp[i] & tot_p <= temp[i+1])),
replace = F) +
rnorm(length(which(tot_p >= temp[i] & tot_p <= temp[i+1])),
0, 0.05)
}
set.seed(555)
beta_veg <- 0.1
set.seed(555)
huevos_obs <- rbinom(n_nidos,
prob = p + var_parche + beta_veg*veg,
size = huevos)
dat <-
list(
n = huevos,
superv = huevos_obs,
parche = rep(1:10, each = 10),
prop_veg = veg,
nido = 1:100,
N = length(huevos),
N_parches = 10
)En principio, ajustemos el modelo anterior, pero incluyamos un parámetro \(\beta\) que controle la relación entre la supervivencia de huevos y la proporción de cobertura vegetal. Esto es:
\[ \begin{align} &Supervivencia_i \sim Binomial(huevos_i,~p_i) \\ &logit(p_i) \sim \alpha[nido_i] + \psi[parche] + \beta\times{prop. vegetación}\\ &\beta \sim Normal(0.5, 1)\\ &\psi[parche] \sim Normal(0, \tau) \\ &\tau \sim exponential(1) \\ &\alpha[nido] \sim Normal(\mu_j, \sigma_j) \\ &\mu_j \sim Normal(5, 2.5) \\ &\sigma_j \sim exponential(1) \\ \end{align} \]
cat(file = 'par_pooling3.stan',
'
data{
int N;
int N_parches;
array[N] int n;
array[N] int superv;
array[N] int parche;
vector[N] prop_veg;
array[N] int nido;
}
parameters{
vector[N] alpha;
vector[N_parches] psi;
real mu;
real beta;
real<lower = 0> sigma;
real<lower = 0> tau;
}
model {
vector[N] p;
for (i in 1:N) {
p[i] = alpha[nido[i]] + psi[parche[i]] + beta*prop_veg[i];
p[i] = inv_logit(p[i]);
}
alpha ~ normal(mu, sigma);
psi ~ normal(0, tau);
mu ~ normal(10, 5);
beta ~ normal(0.5, 1);
sigma ~ exponential(1);
tau ~ exponential(1);
superv ~ binomial(n, p);
}
generated quantities{
vector[N] log_lik;
vector[N] p;
for (i in 1:N){
p[i] = alpha[nido[i]] + psi[parche[i]] + beta*prop_veg[i];
p[i] = inv_logit(p[i]);
}
for (i in 1:N) log_lik[i] = binomial_lpmf(superv[i] | n[i], p[i]);
}
')file <- paste(getwd(), '/par_pooling3.stan', sep = '')
fit_par_pooling3 <- cmdstan_model(file, compile = T)Warning in readLines(stan_file): incomplete final line found on
'D:/github_repos/RPubs/multilevel_mods/par_pooling3.stan'mod2.1 <-
fit_par_pooling3$sample(
data = dat,
iter_sampling = 2e3,
iter_warmup = 500,
chains = 3,
parallel_chains = 3,
thin = 3,
refresh = 500,
seed = 123
)
(out_par_pool3 <- mod2.1$summary())
par(mfrow = c(3, 3), mar = c(4, 4, 1, 1))
for (i in 1:9) trace_plot(mod2.1, out_par_pool3$variable[i], 3)
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
El modelo tiene un un número tal de transiciones divergente que el effective sampling size de algunos parámetros es inferior a 1000 y por lo tanto no deberíamos estar muy confiados de esta estimación. Esto es común en modelos jerárquicos, y sucede porque el algoritmo explora incorrectamente algunas regiones del espacio de los parámetros. Podríamos intentar solucionarlo de varias formas: (i) incrementando el número de iteraciones (no es lo ideal porque no lidiaríamos con el problema de eficiencia en la exploración del algoritmo); (ii) incrementando la tasa de aceptación de valores que explora el algoritmo controlando los “saltos” entre uno y otro valor estimado (leapfrog step); (iii) reparametrización del modelo. Probemos la alternativa ii incrementando el adapt_delta = 0.99:
mod2.2 <-
fit_par_pooling3$sample(
data = dat,
iter_sampling = 2e3,
iter_warmup = 500,
chains = 3,
parallel_chains = 3,
thin = 3,
refresh = 500,
seed = 123,
adapt_delta = 0.99
)
(out_par_pool4 <- mod2.2$summary())
par(mfrow = c(3, 3), mar = c(4, 4, 1, 1))
for (i in 1:9) trace_plot(mod2.2, out_par_pool4$variable[i], 3)
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Bien, tenemos menos transiciones divergentes, pero el effective sampling size continua bajo para algunos coeficientes. Probemos entonces reparametrizando el modelo, específicamente una reparametrización no centrada (non-centered reparametrization). Esta consiste en una estandarización que permite que el algoritmo explore el espacio de los parámetros de manera más eficiente. Veamos:
cat(file = 'par_pooling4.stan',
'
data{
int N;
int N_parches;
array[N] int n;
array[N] int superv;
array[N] int parche;
vector[N] prop_veg;
array[N] int nido;
}
parameters{
vector[N] z_alpha;
vector[N_parches] z_psi;
real mu1;
real mu2;
real beta;
real<lower = 0> sigma;
real<lower = 0> tau;
}
model{
vector[N] p;
vector[N] alpha;
vector[N_parches] psi;
mu1 ~ normal(10, 5);
mu1 ~ normal(0, 5);
z_alpha ~ normal(0, 1);
z_psi ~ normal(0, 1);
sigma ~ exponential(1);
tau ~ exponential(1);
beta ~ normal(0.5, 1);
alpha = mu1 + z_alpha*sigma;
psi = mu2 + z_psi*tau;
for (i in 1:N) {
p[i] = alpha[nido[i]] + psi[parche[i]] + beta*prop_veg[i];
p[i] = inv_logit(p[i]);
}
superv ~ binomial(n, p);
}
generated quantities{
vector[N] log_lik;
vector[N] p;
vector[N] alpha;
vector[N_parches] psi;
alpha = mu1 + z_alpha*sigma;
psi = mu2 + z_psi*tau;
for (i in 1:N) {
p[i] = alpha[nido[i]] + psi[parche[i]] + beta*prop_veg[i];
p[i] = inv_logit(p[i]);
}
for (i in 1:N) log_lik[i] = binomial_lpmf(superv[i] | n[i], p[i]);
}
')file <- paste(getwd(), '/par_pooling4.stan', sep = '')
fit_par_pooling4 <- cmdstan_model(file, compile = T)Warning in readLines(stan_file): incomplete final line found on
'D:/github_repos/RPubs/multilevel_mods/par_pooling4.stan'mod2.3 <-
fit_par_pooling4$sample(
data = dat,
iter_warmup = 500,
iter_sampling = 2e3,
chains = 3,
parallel_chains = 3,
thin = 3,
refresh = 500,
seed = 123
)
(out_par_pool5 <- mod2.3$summary())
par(mfrow = c(3, 3), mar = c(4, 4, 1, 1))
for (i in 1:9) trace_plot(mod2.3, out_par_pool5$variable[i], 3)
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Las transiciones divergentes persisten, pero ahora tenemos suficientes estimados de los parámetros estadísticamente independientes.
Veamos ahora el ajuste a los datos del modelo reparametrizado (i.e. distribución predictiva posterior). Además, verifiquemos que podemos recobrar la distribución de probabilidad generadora de los datos (\(Beta(3, 5)~+~\sigma,~\sigma \sim Normal(0, \phi)\))
post <- mod2.3$draws(format = 'df')
est_p <- sapply(1:100, simplify = 'array', FUN =
function(i) {
# equivalente a usar link(mod2.3, dat)
i_alpha <- dat$nido[i]
i_psi <- dat$parche[i]
p <- post[, grepl('^alpha', colnames(post))][[i_alpha]] +
post[, grepl('^psi', colnames(post))][[i_psi]] +
(post$beta*dat$prop_veg[i])
inv_logit(p)
})
est <- sapply(1:100, simplify = 'array', FUN =
function(x) {
i_alpha <- dat$nido[x]
i_psi <- dat$parche[x]
p <-
post[, grepl('^alpha', colnames(post))][[i_alpha]] +
post[, grepl('^psi', colnames(post))][[i_psi]] +
(post$beta*dat$prop_veg[i])
rbinom(1e3, size = dat$n[x], prob = inv_logit(p))
# equivalente a sim(mod2.3, dat)
})En efecto el modelo está adecuadamente definido.
Comparemos ahora las distribuciones posteriores por nido y parche
par(mfrow = c(1, 3), mar = c(4, 4, 1, 1))
plot(density(dat$superv), ylim = c(0, 0.25), main = 'ppchecks',
xlab = 'N huevos')
for (i in 1:100) lines(density(est[i, ]), lwd = 0.1)
lines(density(dat$superv), col = 'red', lwd = 3)
plot(NULL, xlim = c(0, 1), ylim = c(0, 4), ylab = 'Densidad',
xlab = expression(italic(P)['nido']))
for (i in 1:100) {
lines(density(inv_logit(post[, grepl('alpha', colnames(post))][[i]])),
lwd = 0.3)
}
plot(NULL, xlim = c(0, 1), ylim = c(0, 4), ylab = 'Densidad',
xlab = expression(italic(P)['parche']))
for (i in 1:10) {
lines(density(inv_logit(post[, grepl('psi', colnames(post))][[i]])),
lwd = 0.3)
}
La probabilidad de supervivencia de los huevos es notablemente más variable entre nidos.
Ya sabemos que la cobertura vegetal tiene un efecto positivo sobre la supervivencia de los huevos, nosotros la simulamos, pero veamos sus implicaciones para el modelo.
v <- seq(0, 1, length.out = 100)
b_veg <- sapply(1:100, simplify = 'array', FUN =
function(i) {
inv_logit(post$beta*v[i])
})
ci_veg <- apply(b_veg, 2, quantile, c(0.025, 0.975))
plot(v, apply(b_veg, 2, mean), type = 'l', xlab = 'Cobertura de vegetación (prop)',
ylab = 'Probabilidad de supervivencia')
shade(ci_veg, v, col = col.alpha('cyan4'))
Graficamos el efecto neto de la cobertura vegetación, también podríamos estratificarlo por cada nido o parche de bosque.
4 Ejemplo
Las simulaciones estocásticas nos permite plantear el modelo generador de los datos y estudiar el desempeño de un modelo estadístico, sus supuestos y alcances, para recuperar dichos parámetros iniciales. Es sin duda ilustrativo. Ahora, apliquemos lo ya expuesto a un set de datos empírico. En este caso usaremos datos de un censo de fertilidad, específicamente el uso de anticonceptivos por mujeres bengalíes (Bangladesh) en 1988. Carguemos rethinking::bangladesh y exploremos las variables:
[1] "woman" "district" "use.contraception"
[4] "living.children" "age.centered" "urban"
[7] "district_id" d contiene 6 variables:
woman: identificador de la mujer.distric: distrito al que pertenece.use.contraception: uso (1) o no (0) de anticonceptivos.living.children: número de hijos vivos.age.centered: edad centrada.urban: contexto urbano (1) o no urbano (0). En nuestro caso lo llamaremos rural.
Nuestro objetivo será evaluar cómo el distrito, el contexto y la edad, modulan la probabilidad de que una mujer bengalí, en 1988, usara un método anticonceptivo. Definamos primero el modelo en notación matemática:
\[ \begin{align} &Uso~anticonceptivos \sim Binomial(1, P) \\ &logit(P_i) = \alpha_{urbano_i} + \theta_{distrito_i} + \beta\times edad \\ &\alpha_{urbano_i} \sim Normal(\mu, \sigma) \\ &\theta_{distrito_i} \sim Normal(0, \phi) \\ &\mu \sim Normal(0, 3) \\ &\sigma \sim exponential(1) \\ &\phi \sim exponential(1)\\ &\beta \sim Normal(0.5, 0.2) \end{align} \] Detengámonos un momento en las probabilidades previas, y discutamos mi razonamiento para su elección. Los centros urbanos en general tienen un mayor acceso a establecimientos educativos y de salud. Sería razonable esperar que la probabilidad de usar anticonceptivos sea más baja en contextos rurales comparada con los urbanos. La previa \(\mu \sim Normal(0, 3)\) implica justamente esto en la escala logit y de probabilidad, veamos:
par(mfrow = c(1, 2))
curve(dnorm(x, 0, 1.5), xlim = c(-5, 5), xlab = 'Log-odds(P)', ylab = 'Density')
plot(density(inv_logit(rnorm(1e3, 0, 1.5))), xlab = 'P', ylab = 'Density',
main = '')
Esta previa permite que el algoritmo explore, en escala logit, tanto valores negativos como positivos restringidos a \(-3\) y \(+3\). En probabilidades, implica todo el rango de probabilidades entre 0 y 1.
El contexto social, embebido en el espacial, probablemente modula la elección que hacen las mujeres. Esto es, más allá de vivir en un ambiente rural o urbano, variaciones en el sistema de creencias entre las comunidades en las que viven las mujeres, nos llevarían a esperar oscilaciones en la probabilidad de uso de anticonceptivos de acuerdo al distrito. La previa que hace justamente incorpora dicha variación eso \(\theta_{distrito_i} \sim Normal(0, \phi)\).
Finalmente, con el tiempo llegan la experiencia, madurez, la introspección (PONELE). Sea por razones éticas, ambientales, filosóficas o simplemente cansancio, es razonable suponer que la probabilidad de usar métodos anticonceptivos es mayor en mujeres más adultas. Por lo tanto, empleamos \(\beta \sim Normal(0.15, 0.2)\) una previa que permite al algoritmo explorar principalmente valores positivos de \(\beta\), pero permaneciendo relativamente conservador. Veamos:
betas_i <- rnorm(100, 0.15, 0.2)
plot(NULL, xlim = c(-10, 10), ylim = c(0, 1), xlab = 'Edad', ylab = 'Probabilidad')
for (i in 1:100) curve(inv_logit(0 + betas_i[i]*x), add = T, lwd = 0.1)
curve(inv_logit(0 + mean(betas_i)*x), add = T, lwd = 2, col = 'red')
Una vez definidas las previas, agrupemos las variables en una lista y ajustemos el modelo:
cat(file = 'multilevel_mods/par_pooling5.stan',
'
data{
int N;
int N_distrito;
int N_urbano;
int N_anticonc;
array[N] int anticonc;
array[N] int distrito;
array[N] int urbano;
//vector[N] edad;
}
parameters{
vector[N_distrito] z_alpha;
vector[N_distrito] z_beta;
real beta1;
real alpha1;
real<lower = 0> sigma;
real<lower = 0> phi;
}
model{
vector[N] p;
vector[N_distrito] alpha;
vector[N_distrito] beta;
alpha1 ~ normal(0, 1);
z_alpha ~ normal(0, 1);
sigma ~ exponential(1);
beta1 ~ normal(0, 1);
z_beta ~ normal(0, 1);
phi ~ exponential(1);
alpha = alpha1 + z_alpha*sigma;
beta = beta1 + z_beta*phi;
for (i in 1:N) {
p[i] = alpha[distrito[i]] + beta[distrito[i]]*urbano[i];
p[i] = inv_logit(p[i]);
}
anticonc ~ binomial(1, p);
}
generated quantities{
vector[N] log_lik;
vector[N] p;
vector[N_distrito] alpha;
vector[N_distrito] beta;
alpha = alpha1 + z_alpha*sigma;
beta = beta1 + z_beta*phi;
for (i in 1:N) {
p[i] = alpha[distrito[i]] + beta[distrito[i]]*urbano[i];
p[i] = inv_logit(p[i]);
}
for (i in 1:N) log_lik[i] = binomial_lpmf(anticonc[i] | 1, p[i]);
}
')dat <-
list(
anticonc = d$use.contraception,
urbano = as.integer(ifelse(d$urban == 1, 1, 2)),
distrito = d$district_id,
edad = d$age.centered,
N = nrow(d),
N_distrito = length(unique(d$district_id)),
N_urbano = 2,
N_anticonc = 2
)
file <- paste(getwd(), '/par_pooling5.stan', sep = '')
fit_par_pooling5 <- cmdstan_model(file, compile = T)
me2 <-
fit_par_pooling5$sample(
data = dat,
iter_sampling = 4e3,
iter_warmup = 500,
parallel_chains = 3,
chains = 3,
thin = 3,
seed = 123,
refresh = 500
)
(out_par_pool6 <- me2$summary())
par(mfrow = c(3, 3), mar = c(4, 4, 1, 1))
for (i in 1:9) trace_plot(me2, out_par_pool6$variable[i], 3)
Todos los parámetros tienen más de mil estimados estadísticamente independientes y Rhat ~ 1. Podríamos —más bien deberíamos, graficar las cadenas para corroborar que convergieron a la misma distribución posterior y exploraron el mismo espacio del parámetro. Pero, en este caso, confiaremos en la salida del modelo.
Verifiquemos el ajuste del modelo a los datos:
post2 <- me2$draws(format = 'df')
ppcheck <- sapply(1:length(dat$anticonc), simplify = 'array', FUN =
function(x) {
urbano <- dat$urbano[x]
distrito <- dat$distrito[x]
p <-
post2[, grepl('^alpha', colnames(post2))][[distrito]] +
post2[, grepl('^beta', colnames(post2))][[distrito]]*urbano
rbinom(1e3, 1, inv_logit(p))
})
plot(NULL, xlim = c(-0.3, 1.3), ylim = c(0, 3.2),
xlab = 'Uso de anticonceptivo', ylab = 'Densidad')
for (i in 1:200) lines(density(ppcheck[i, ]), lwd = 0.1)
lines(density(dat$anticonc), col = 'red', lwd = 2)
El modelo es confiable. Grafiquemos ahora los efectos condicionales, y verifiquemos las implicaciones del distrito y el contexto espacial en la probabilidad de que las mujeres usen algún método anticonceptivo. Primero realicemos cálculos y un poco de carpintería:
bar_urb <- unique(tibble(dist = dat$distrito,
urb = dat$urbano))
post_urb_dist <- sapply(1:nrow(bar_urb), simplify = 'array', FUN =
function(x) {
dist <- bar_urb$dist[x]
urb <- bar_urb$urb[x]
p <- post2[, grepl('^alpha', colnames(post2))][[dist]] +
post2[, grepl('^beta', colnames(post2))][[dist]]*urb
inv_logit(p)
})
bar_urb$level <- ifelse(bar_urb$urb == 1, 'Urbano', 'Rural')
post_urb_dist <- as_tibble(post_urb_dist)
colnames(post_urb_dist) <- paste(bar_urb$level, bar_urb$dist, sep = '')
post_urb_dist <- gather(post_urb_dist)
alpha_urb <- post_urb_dist
post_urb_dist <-
post_urb_dist |>
group_by(key) |>
transmute(li = quantile(value, 0.025),
ls = quantile(value, 0.975),
mu = mean(value)) |>
ungroup() |>
unique()
post_urb_dist$dist <- gsub('^([aA-zZ]*)([0-9]{1,})', '\\2', post_urb_dist$key)
post_urb_dist$factor <- gsub('^([aA-zZ]*)([0-9]{1,})', '\\1', post_urb_dist$key)
alpha_urb$key <- ifelse(grepl('^Urban', alpha_urb$key), 'Urbano', 'Rural')
contraste <- split(alpha_urb, alpha_urb$key)
contraste <- tibble(x = sample(contraste$Urbano$value, 1e4) - sample(contraste$Rural$value, 1e4))Ahora el gráfico
plot_contraste <-
ggplot(contraste, aes(x)) +
geom_density(fill = 'seagreen', color = 'seagreen', alpha = 0.5) +
labs(x = quote(italic(P['(Diferencia urbano-rural)'])), y = 'Densidad') +
geom_vline(xintercept = mean(contraste$x), linetype = 3) +
theme_bw() +
theme(panel.grid = element_blank(),
text = element_text(family = 'Times New Roman'))
plot_cond <-
post_urb_dist |>
ggplot(aes(xmin = li, xmax = ls, x = mu,
y = fct_reorder(dist, mu, .desc = T), color = factor)) +
geom_point(position = position_dodge(width = 0.3)) +
geom_linerange(position = position_dodge(width = 0.3)) +
labs(x = quote(italic(P['(usar anticonceptivos)'])), y = 'Distrito') +
theme_bw() +
theme(panel.grid = element_blank(),
legend.title = element_blank(),
legend.position = 'top',
text = element_text(family = 'Times New Roman'))
plot_contraste +
plot_cond +
plot_layout(ncol = 2)
Concluimos que las mujeres en zonas urbanas tienen, en general, un 15% más de probabilidad de usar anticonceptivos comparadas con mujeres en ambientes rurales (figura izquierda). Pero también observamos que, además del contexto espacial, dicha probabilidad está fuertemente influenciada por el distrito (figura derecha).
Podría dedicar una sección a discutir las implicaciones de los resultados para una intervención, distrito específica, para fomentar el uso de anticonceptivos en mujeres bengalíes allá en 1988, pero no es el objetivo del ejercicio. Eso se lo dejo a quién sea que haya llegado hasta acá.
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=Spanish_Argentina.utf8 LC_CTYPE=Spanish_Argentina.utf8
[3] LC_MONETARY=Spanish_Argentina.utf8 LC_NUMERIC=C
[5] LC_TIME=Spanish_Argentina.utf8
time zone: America/Buenos_Aires
tzcode source: internal
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] loo_2.6.0 rstan_2.32.5 StanHeaders_2.32.5 patchwork_1.2.0
[5] magrittr_2.0.3 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
[9] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 tidyr_1.3.0
[13] tibble_3.2.1 ggplot2_3.4.4 tidyverse_2.0.0 rethinking_2.40
[17] posterior_1.5.0 cmdstanr_0.7.1
loaded via a namespace (and not attached):
[1] gtable_0.3.4 shape_1.4.6 tensorA_0.36.2.1
[4] QuickJSR_1.0.9 xfun_0.41 htmlwidgets_1.6.4
[7] processx_3.8.3 inline_0.3.19 lattice_0.21-8
[10] tzdb_0.4.0 vctrs_0.6.5 tools_4.3.1
[13] ps_1.7.5 generics_0.1.3 curl_5.2.0
[16] stats4_4.3.1 fansi_1.0.6 pkgconfig_2.0.3
[19] Matrix_1.6-5 data.table_1.14.10 checkmate_2.3.1
[22] distributional_0.3.2 RcppParallel_5.1.7 lifecycle_1.0.4
[25] compiler_4.3.1 farver_2.1.1 munsell_0.5.0
[28] codetools_0.2-19 htmltools_0.5.7 yaml_2.3.8
[31] pillar_1.9.0 MASS_7.3-60 abind_1.4-5
[34] nlme_3.1-162 tidyselect_1.2.0 digest_0.6.34
[37] mvtnorm_1.2-4 stringi_1.8.3 splines_4.3.1
[40] labeling_0.4.3 fastmap_1.1.1 grid_4.3.1
[43] colorspace_2.1-0 cli_3.6.2 pkgbuild_1.4.3
[46] utf8_1.2.4 withr_3.0.0 scales_1.3.0
[49] backports_1.4.1 timechange_0.2.0 rmarkdown_2.25
[52] matrixStats_1.2.0 gridExtra_2.3 hms_1.1.3
[55] coda_0.19-4 evaluate_0.23 knitr_1.45
[58] V8_4.4.1 mgcv_1.8-42 rlang_1.1.3
[61] Rcpp_1.0.12 glue_1.7.0 rstudioapi_0.15.0
[64] jsonlite_1.8.8 R6_2.5.1