Introducción: modelos lineales generalizados (GLM)
1 Introducción
En el modelo lineal utilizamos una función de likelihood normal para modelar el parámetro \(\mu\); esto es, asumimos que los valores de la variable respuesta son continuos no acotados a un máximo o límites particulares. En notación matemática:
\[ \begin{aligned} &Y_{i} \sim Normal(\mu_i, \sigma) \\ &\mu_i = \alpha + \beta{x} \end{aligned} \]
En efecto, tamaño, área, temperatura, etcétera, son variables que sin menor dificultad podríamos modelar con una distribución normal. Sin embargo, frecuentemente en ecología analizamos conteos (e.g. abundancias), probabilidades (e.g detección de una especie), o proporciones (e.g. remoción de frutos), etcétera. Modelar dichas variables con una función de likelihood normal resultaría en una imprecisa estimación de los parámetros (Unboxing 1).
En aquellos casos en los que la distribución normal no es la mejor elección, necesitamos una generalización del modelo en términos de la función de likelihood que empleamos para denotar los valores posibles de los datos (i.e. la probabilidad de los datos) y una función de enlace (link function). Entonces, el modelo lineal generalizado (GLM, por sus siglas en inglés) consiste en:
\[ \begin{align} &Y_i \sim likelihood(\phi_i,~\theta) \\ &f(\phi_i) = \alpha + \beta{x} \end{align} \]
Aquí reemplazamos la distribución normal, \(Normal(\mu_i, ~ \sigma)\), por una nueva función de likelihood, en este caso \(likelihood(\phi_i, ~ \theta)\). Elegimos la nueva función basados en nuestro entendimiento de la variable que planeamos modelar (e.g. se trata de conteos, probabilidades, proporciones, variables positivas con valores extremos, variables ordinales) —el paquete fitdistrplus puede facilitar la tarea, ya que permite comparar la distribución de nuestra variable con el esperado teórico por distintas distribuciones de probabilidad. Seguidamente, en lugar de modelar el parámetro \(\phi\) directamente, empleamos la función de enlace \(f\) para estimarlo en la escala del modelo lineal \(\alpha + \beta{x}\).
En GLMs las funciones de enlace más comunes son la función log y logit. En casos en los que la variable respuesta está restringida a valores positivos sin un límite superior definido, empleamos el logaritmo como función de enlace. Esto es, modelamos \(y\) en escala lineal \(log\) y la interpretación en la escala de la variable respuesta corresponde a su inverso matemático \(e^y\). Esto implica que un cambio en una unidad de \(x\) no produce un cambio constante en \(y\), como en el modelo lineal. En este caso, emplear una función de enlace log supone una relación exponencial \(y\sim x\). Veamos:
La función de enlace logit se emplea para modelar probabilidades y proporciones. Justamente, su propósito es evitar que el modelo \(\alpha + \beta{x}\) estime valores de \(y\) fuera del intervalo 0 y 1. El logit o log-odds, se estima como el logaritmo del cociente entre la probabilidad de que el evento ocurra (\(p\)) y no ocurra (\(1-p\)): \[ \begin{align} logit(p_i) = log(\frac{p_i}{1-p_i}) \end{align} \] Entonces,
\[ log(\frac{p_i}{1-p_i}) = \alpha + \beta{x_i} \] Y la función logística del logit, su inverso, corresponde a:
\[ p_i = \frac{exp(\alpha + \beta{x_i})}{1 + exp(\alpha + \beta{x_i})} \]
El log-odds no está restringido a la escala de probabilidad, por el contrario se extiende hacia \(-\infty\) y \(+\infty\). Valores de cero indican una probabilidad de 0.5 e, igual que con la función de enlace log —en realidad en cualquier GLM— los cambios en la variable predictora no generan cambios constantes en la respuesta. Veamos un ejemplo:
2 Ejemplo: modelo logístico
Emplearemos datos de capturas de roedores con trampas Sherman. El estudio se desarrolló en diferentes parches de bosque, y el objetivo es modelar la probabilidad de captura en función a la cobertura de bosque y los años transcurridos desde su fragmentación.
Primero cargamos los datos en R y almacenamos las variables en una lista.
Definimos la previa para el parámetro \(\alpha\) (intercepto) del modelo lineal, considerando los valores posibles de log-odds que podría adquirir y sus implicaciones en términos de probabilidad.
par(mfrow = c(1, 2))
plot(density((rnorm(1e3, 0, 1.5))), xlab = 'log-odds', main = 'Logit')
plot(density(inv_logit(rnorm(1e3, 0, 1.5))), xlab = 'Probabilidad', main = 'Inv. Logit')
Tampoco tenemos muchas certezas sobre la relación entre la relación \(P(capturas) \sim vegetación + años\), pero podríamos suponer que una vegetación más abundante favorece la actividad de roedores y por lo tanto sus capturas. No estamos muy confiados, así que utilicemos \(Normal(\mu=0.5, \alpha = 1)\), y veamos los supuestos que implican nuestras previas en términos del modelo lineal.
alpha <- rnorm(1e3, 0, 1.5)
beta <- rnorm(1e3, 0.5, 1)
par(mfrow = c(1, 2))
plot(NULL, xlim = c(-20, 20), ylim = c(-20, 20),
ylab = 'Log-odds', xlab = 'X')
for (i in 1:100) curve(alpha[i] + beta[i]*x, add = T, lwd = 0.2)
plot(NULL, xlim = c(-20, 20), ylim = c(0, 1),
ylab = 'Probabilidad', xlab = 'X')
for (i in 1:100) curve(inv_logit(alpha[i] + beta[i]*x), add = T, lwd = 0.2)
En efecto, esperamos una relación positiva, pero las previas son lo suficientemente amplias como para esperar relaciones inversas. Ajustemos el modelo:
cat(file = 'binom1.stan',
'
data{
int N;
array[N] int obs;
array[N] int edad;
array[N] int vege;
}
parameters {
real alpha;
real beta;
real tau;
}
model{
vector[N] p;
for (i in 1:N) {
p[i] = alpha + beta*edad[i] + tau*vege[i];
p[i] = inv_logit(p[i]);
}
obs ~ binomial(1, p);
alpha ~ normal(0, 1.5);
beta ~ normal(0.5, 1);
tau ~ normal(0.5, 1);
}
')[1] "D:/github_repos/RPubs/intro_glm/binom1.stan"[1] "D:/github_repos/RPubs/intro_glm/binom1.stan"fit_m <- cmdstan_model(file, compile = T)
m <-
fit_m$sample(
data = dat,
iter_sampling = 2e3,
iter_warmup = 500,
thin = 3,
chains = 3,
parallel_chains = 3,
refresh = 500,
seed = 123
)Running MCMC with 3 parallel chains...
Chain 1 Iteration: 1 / 2500 [ 0%] (Warmup)
Chain 1 Iteration: 500 / 2500 [ 20%] (Warmup)
Chain 1 Iteration: 501 / 2500 [ 20%] (Sampling)
Chain 1 Iteration: 1000 / 2500 [ 40%] (Sampling)
Chain 2 Iteration: 1 / 2500 [ 0%] (Warmup)
Chain 2 Iteration: 500 / 2500 [ 20%] (Warmup)
Chain 2 Iteration: 501 / 2500 [ 20%] (Sampling)
Chain 3 Iteration: 1 / 2500 [ 0%] (Warmup)
Chain 3 Iteration: 500 / 2500 [ 20%] (Warmup)
Chain 3 Iteration: 501 / 2500 [ 20%] (Sampling)
Chain 1 Iteration: 1500 / 2500 [ 60%] (Sampling)
Chain 1 Iteration: 2000 / 2500 [ 80%] (Sampling)
Chain 1 Iteration: 2500 / 2500 [100%] (Sampling)
Chain 1 finished in 0.2 seconds.
Chain 2 Iteration: 1000 / 2500 [ 40%] (Sampling)
Chain 2 Iteration: 1500 / 2500 [ 60%] (Sampling)
Chain 2 Iteration: 2000 / 2500 [ 80%] (Sampling)
Chain 2 Iteration: 2500 / 2500 [100%] (Sampling)
Chain 2 finished in 0.2 seconds.
Chain 3 Iteration: 1000 / 2500 [ 40%] (Sampling)
Chain 3 Iteration: 1500 / 2500 [ 60%] (Sampling)
Chain 3 Iteration: 2000 / 2500 [ 80%] (Sampling)
Chain 3 Iteration: 2500 / 2500 [100%] (Sampling)
Chain 3 finished in 0.2 seconds.
All 3 chains finished successfully.
Mean chain execution time: 0.2 seconds.
Total execution time: 0.5 seconds.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')
}
}
(m_out <- m$summary())# A tibble: 4 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ -13.3 -12.9 1.27 0.994 -15.7 -11.9 1.00 1706.
2 alpha -1.53 -1.53 1.25 1.30 -3.55 0.529 1.00 1782.
3 beta -0.0280 -0.0263 0.0256 0.0247 -0.0729 0.0104 1.00 1649.
4 tau 0.0506 0.0501 0.0195 0.0190 0.0193 0.0848 1.00 1697.
# ℹ 1 more variable: ess_tail <dbl>
No soy fanático de interpretar modelos con tablas de coeficientes. Por ahora, corroboremos que, para cada parámetro, el número de estimados estadísticamente independientes supera los 1000 y las cadenas convergieron a la “misma” distribución posterior (Rhat ~ 1). Todo en orden
Grafiquemos las distribuciones predictivas posteriores y veamos la calidad del ajuste del modelo a los datos.
post <- m$draws(format = 'df')
ppcheck <- sapply(1:length(dat$obs), simplify = 'array', FUN =
function(x) {
log_odds <-
with(post,
{
alpha + beta*dat$edad[x] + tau*dat$vege[x]
})
rbinom(1e3, 1, inv_logit(log_odds))
})
plot(NULL, xlim = c(-0.5, 1.5), ylim = c(0, 1.7))
for (i in 1:500) lines(density(ppcheck[i, ]), lwd = 0.1)
lines(density(dat$obs), col = 'red', lwd = 3)
Todo en orden.
veamos ahora los efectos condicionales de cada predictor sobre la probabilidad de captura de los roedores
edad <- seq(min(dat$edad), max(dat$edad), length.out = 100)
est_beta <-
sapply(edad, simplify = 'array', FUN =
function(x) {
log_odds <-
with(post,
{
alpha + beta*x
})
inv_logit(log_odds)
})
est_beta <-
do.call('rbind',
apply(est_beta, 2, simplify = 'list', FUN =
function(x) {
tibble(li = quantile(x, 0.025),
ls = quantile(x, 0.975),
mu = mean(x))
}))
est_beta$x <- edad
plot_edad <-
est_beta |>
ggplot(aes(x, mu, ymin = li, ymax = ls)) +
geom_ribbon(alpha = 0.5, fill = 'cyan4') +
geom_line(color = 'tan1', linewidth = 1.5) +
theme_bw() +
lims(y = c(0, 1)) +
labs(y = 'Probabilidad', x = 'Años desde la fragmentación') +
theme(panel.grid = element_blank(),
text = element_text(family = 'Times New Roman'))
veg <- seq(min(dat$vege), max(dat$vege), length.out = 100)
est_tau <-
sapply(veg, simplify = 'array', FUN =
function(x) {
log_odds <-
with(post,
{
alpha + tau*x
})
inv_logit(log_odds)
})
est_tau <-
do.call('rbind',
apply(est_tau, 2, simplify = 'list', FUN =
function(x) {
tibble(li = quantile(x, 0.025),
ls = quantile(x, 0.975),
mu = mean(x))
}))
est_tau$x <- veg
plot_vege <-
est_tau |>
ggplot(aes(x, mu, ymin = li, ymax = ls)) +
geom_ribbon(alpha = 0.4, fill = 'cyan4') +
geom_line(color = 'tan1', linewidth = 1.5) +
theme_bw() +
labs(y = ' ', x = 'Porcentaje de vegetación') +
theme(panel.grid = element_blank(),
text = element_text(family = 'Times New Roman'))
plot_edad +
plot_vege +
plot_layout(ncol = 2)
El incremento del porcentaje de vegetación aumenta la probabilidad de capturar roedores, mientras que parches con fragmentación más antigua generan el efecto contrario. Ambos predictores probablemente no son independientes; por ejemplo, el porcentaje de vegetación podría depender de cuán recientemente fue deforestado el parche. Si este fuera el caso, el modelo estaría adecuadamente planteado, ya que estaríamos condicionando la relación \(P(capturas) \sim vegetación\) por la edad de la fragmentación.
Veamos las distribuciones posteriores de los coeficientes de regresión.
plot(density(post$beta), xlim = c(-0.15, 0.16),
ylim = c(0, 20), main = "", lwd = 3,
xlab = 'Coeficiente de regresión', ylab = 'Densidad')
lines(density(post$tau), col = 'tan1', lwd = 3)
abline(v = 0, lty = 3)
text(expression(tau), x = mean(post$tau), y = 10, col = 'tan2')
text(expression(beta), x = mean(post$beta), y = 10)
Sería buena idea comparar nuestra apuesta inicial, i.e. las previas, con lo estimado por el modelo. Veamos cuánto aprendimos:
plot(density(rnorm(1e3, 0.5, 1)), lwd = 2, lty = 3, col = 'red',
main = ' ', xlab = expression('previas'[tau~'-'~beta]))
lines(density(post$beta), xlim = c(-0.15, 0.16), lwd = 1)
lines(density(post$tau), col = 'tan1', lwd = 1)
Las distribuciones posteriores ocupan una región muchísimo menor a lo esperado por la previa.
3 Unboxing 1
Simulemos conteos de gatos (\(15~min~conteo^{-1}\)) en cuatro barrios de Villavicencio (Colombia) — probablemente una de las ciudades más atestadas de gatos en el país, y veamos las implicaciones de usar una función de likelihood normal para modelar conteos.
Simularemos los conteos empleando una distribución de probabilidad \(Poisson(\lambda)\). En la próxima sección entraremos en detalle, por ahora, solo necesitas saber que el parámetro \(\lambda\) puede entenderse como el número esperado de conteos en una unidad muestreal y una tasa de cambio constante de un conteo a otro. Dicho esto, preparemos los datos. Utilizaremos una distribución uniforme para generar \(\lambda\) aleatorios entre 0 y 4. Luego generamos los conteos de gatos.
n_obs <- 50
set.seed(123)
lambdas <- sample(runif(1e3, 0, 4), 4)
barrios <- sapply(lambdas, simplify = 'array', FUN =
function(lambda) {
set.seed(123)
rpois(n_obs, lambda)
})
barrios <- as_tibble(barrios)
colnames(barrios) <- paste(1:4)
barrios <- gather(barrios)
barrios$key <- as.numeric(barrios$key)Conocemos la ciudad, y suponemos que es probable observar ~7 gatos en cada una de las observaciones. Definimos entonces las probabilidades previas.
plot(density(rnorm(1e3, 7, 5)), main = expression(Normal(mu~"=7", sigma~'=5')),
xlab = 'Gatos observados')
abline(v = 0, col = 'red')
Notemos que la previa implicaría que cerca del 8% de de nuestras observaciones (conteos de gatos) puedan adquirir valores negativos, no una muy buena elección. Pero sigamos. Ajustemos los modelos:
cat(file = 'binom2.stan',
'
data{
int N;
int N_barrio;
array[N] int barrio;
array[N] int gatos_norm;
array[N] int gatos_pois;
}
parameters{
vector[N_barrio] alpha_norm;
real<lower = 0> sigma;
vector[N_barrio] alpha_pois;
}
model{
target += normal_lpdf(gatos_norm | alpha_norm[barrio], sigma);
target += normal_lpdf(alpha_norm | 7, 5);
target += exponential_lpdf(sigma | 1);
target += poisson_lpmf(gatos_pois | exp(alpha_pois[barrio]));
target += normal_lpdf(alpha_pois | 7, 5);
}
')[1] "D:/github_repos/RPubs/intro_glm/binom2.stan"fit_norm <- cmdstan_model(file, compile = T)
mod <-
fit_norm$sample(
data = list(gatos_norm = barrios$value,
gatos_pois = barrios$value,
barrio = barrios$key,
N = nrow(barrios),
N_barrio = length(unique(barrios$key))),
iter_sampling = 2e3,
iter_warmup = 500,
thin = 3,
chains = 3,
parallel_chains = 3,
refresh = 500,
seed = 123
)Running MCMC with 3 parallel chains...
Chain 1 Iteration: 1 / 2500 [ 0%] (Warmup)
Chain 2 Iteration: 1 / 2500 [ 0%] (Warmup)
Chain 3 Iteration: 1 / 2500 [ 0%] (Warmup)
Chain 1 Iteration: 500 / 2500 [ 20%] (Warmup)
Chain 1 Iteration: 501 / 2500 [ 20%] (Sampling)
Chain 1 Iteration: 1000 / 2500 [ 40%] (Sampling)
Chain 1 Iteration: 1500 / 2500 [ 60%] (Sampling)
Chain 1 Iteration: 2000 / 2500 [ 80%] (Sampling)
Chain 2 Iteration: 500 / 2500 [ 20%] (Warmup)
Chain 2 Iteration: 501 / 2500 [ 20%] (Sampling)
Chain 2 Iteration: 1000 / 2500 [ 40%] (Sampling)
Chain 2 Iteration: 1500 / 2500 [ 60%] (Sampling)
Chain 2 Iteration: 2000 / 2500 [ 80%] (Sampling)
Chain 3 Iteration: 500 / 2500 [ 20%] (Warmup)
Chain 3 Iteration: 501 / 2500 [ 20%] (Sampling)
Chain 3 Iteration: 1000 / 2500 [ 40%] (Sampling)
Chain 3 Iteration: 1500 / 2500 [ 60%] (Sampling)
Chain 3 Iteration: 2000 / 2500 [ 80%] (Sampling)
Chain 1 Iteration: 2500 / 2500 [100%] (Sampling)
Chain 2 Iteration: 2500 / 2500 [100%] (Sampling)
Chain 3 Iteration: 2500 / 2500 [100%] (Sampling)
Chain 1 finished in 0.6 seconds.
Chain 2 finished in 0.6 seconds.
Chain 3 finished in 0.6 seconds.
All 3 chains finished successfully.
Mean chain execution time: 0.6 seconds.
Total execution time: 0.7 seconds.# A tibble: 10 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ -702. -701. 2.15 2.01 -706. -699. 0.999 1751.
2 alpha_norm[… 2.27 2.27 0.202 0.198 1.94 2.61 1.00 2132.
3 alpha_norm[… 3.73 3.73 0.213 0.216 3.39 4.08 1.00 1973.
4 alpha_norm[… 1.75 1.75 0.201 0.196 1.40 2.09 1.00 2004.
5 alpha_norm[… 0.489 0.496 0.211 0.216 0.141 0.839 0.999 2074.
6 sigma 1.44 1.44 0.0727 0.0719 1.32 1.56 1.00 2095.
7 alpha_pois[… 0.812 0.813 0.0937 0.0942 0.656 0.964 1.00 2178.
8 alpha_pois[… 1.32 1.32 0.0730 0.0699 1.20 1.43 1.00 2159.
9 alpha_pois[… 0.552 0.554 0.108 0.110 0.368 0.724 1.00 2068.
10 alpha_pois[… -0.742 -0.734 0.202 0.194 -1.09 -0.423 1.00 1925.
# ℹ 1 more variable: ess_tail <dbl>
Salida de los modelos:
# A tibble: 10 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ -702. -701. 2.15 2.01 -706. -699. 0.999 1751.
2 alpha_norm[… 2.27 2.27 0.202 0.198 1.94 2.61 1.00 2132.
3 alpha_norm[… 3.73 3.73 0.213 0.216 3.39 4.08 1.00 1973.
4 alpha_norm[… 1.75 1.75 0.201 0.196 1.40 2.09 1.00 2004.
5 alpha_norm[… 0.489 0.496 0.211 0.216 0.141 0.839 0.999 2074.
6 sigma 1.44 1.44 0.0727 0.0719 1.32 1.56 1.00 2095.
7 alpha_pois[… 0.812 0.813 0.0937 0.0942 0.656 0.964 1.00 2178.
8 alpha_pois[… 1.32 1.32 0.0730 0.0699 1.20 1.43 1.00 2159.
9 alpha_pois[… 0.552 0.554 0.108 0.110 0.368 0.724 1.00 2068.
10 alpha_pois[… -0.742 -0.734 0.202 0.194 -1.09 -0.423 1.00 1925.
# ℹ 1 more variable: ess_tail <dbl>En principio parece que ambos modelos estiman de manera similar los gatos observados promedio. Sin embargo, grafiquemos las distribuciones predictivas posteriores.
post_pois <- mod$draws(format = 'df')[, 7:10]
post_pois <- apply(post_pois, 2, FUN =
function(x) {
set.seed(555)
rpois(1e3,
exp(x))
})
post_norm <- mod$draws(format = 'df')[, 2:6]
post_norm <- apply(post_norm[, -5], 2, FUN =
function(x) {
set.seed(555)
rnorm(1e3, x, post_norm$sigma)
})
post_norm <- gather(as_tibble(post_norm))
post_norm$key <- gsub('(V)(.)', '\\2', post_norm$key)
post_norm$mod <- rep('norm', nrow(post_norm))
barrios$key <- as.factor(as.character(barrios$key))
barrios$mod <- rep('obs', nrow(barrios))
gatos_obs <- rbind(barrios, post_norm)
post_pois <- gather(as_tibble(post_pois))
post_pois$key <- gsub('(V)(.)', '\\2', post_pois$key)
post_pois$mod <- rep('pois', nrow(post_pois))
gatos_obs <- rbind(gatos_obs, post_pois)
ggplot(gatos_obs, aes(key, value, fill = mod)) +
geom_boxplot() +
guides(fill = guide_legend(title = 'Modelos')) +
labs(x = 'Barrio', y = 'Gatos observados') +
geom_hline(yintercept = 0, linetype = 2, color = 'red')
En efecto, el modelo Poisson no solo es ligeramente más preciso, sino que además predice los conteos de gatos considerando el límite inferior de nuestra variable, 0!. Mientras que el modelo normal nos lleva a estimaciones imposibles.
Por ejemplo, en el barrio 4 el modelo normal predice que más del 30% de los conteos de gatos fueron inferiores a cero.
key mod value
1 alpha_norm[1] norm 0.060
2 alpha_norm[2] norm 0.006
3 alpha_norm[3] norm 0.112
4 alpha_norm[4] norm 0.363
5 1 obs 0.000
6 2 obs 0.000
7 3 obs 0.000
8 4 obs 0.000
9 alpha_pois[1] pois 0.000
10 alpha_pois[2] pois 0.000
11 alpha_pois[3] pois 0.000
12 alpha_pois[4] pois 0.000En este caso la elección incorrecta de una función de likelihood no afectó la estimación promedio de los parámetros, pero en modelos más complejos podría tener implicaciones más serias.
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] patchwork_1.2.0 magrittr_2.0.3 lubridate_1.9.3 forcats_1.0.0
[5] stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
[9] tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.4 tidyverse_2.0.0
[13] rethinking_2.40 posterior_1.5.0 cmdstanr_0.7.1
loaded via a namespace (and not attached):
[1] tensorA_0.36.2.1 utf8_1.2.4 generics_0.1.3
[4] shape_1.4.6 stringi_1.8.3 lattice_0.21-8
[7] hms_1.1.3 digest_0.6.34 timechange_0.2.0
[10] evaluate_0.23 grid_4.3.1 mvtnorm_1.2-4
[13] fastmap_1.1.1 jsonlite_1.8.8 processx_3.8.3
[16] backports_1.4.1 ps_1.7.5 fansi_1.0.6
[19] scales_1.3.0 abind_1.4-5 cli_3.6.2
[22] rlang_1.1.3 munsell_0.5.0 withr_3.0.0
[25] yaml_2.3.8 tools_4.3.1 tzdb_0.4.0
[28] checkmate_2.3.1 coda_0.19-4 colorspace_2.1-0
[31] vctrs_0.6.5 R6_2.5.1 matrixStats_1.2.0
[34] lifecycle_1.0.4 htmlwidgets_1.6.4 MASS_7.3-60
[37] pkgconfig_2.0.3 pillar_1.9.0 gtable_0.3.4
[40] data.table_1.14.10 loo_2.6.0 glue_1.7.0
[43] xfun_0.41 tidyselect_1.2.0 rstudioapi_0.15.0
[46] knitr_1.45 farver_2.1.1 htmltools_0.5.7
[49] labeling_0.4.3 rmarkdown_2.25 compiler_4.3.1
[52] distributional_0.3.2