# Librerías:
library(rstan)
library(dplyr)
library(gt)
library(loo)
options(mc.cores = parallel::detectCores())
datos <- data.frame("temp" = seq(1999, 2006),
"juegos" = c(66,68,80,82,65,66,80,42),
"canastas" = c(554,701,749,868,516,573,978,399),
"intentos" = c(1183,1510,1597,1924,1178,1324,2173,845))
datos <- datos |> mutate(porcentaje = canastas/intentos)Es necesario discutir las posibles implicancias de las elecciones de A y B dada la especificación del modelo. Con un valor de \(\theta_t = 0\), se tiene que \(\pi_t = 0.5\), esto es un \(50\%\) de probabilidad de acertar canasta. En base a un desconocimiento sobre las proporciones de acierto de los jugadores de basketball, es razonble asumir una posición conservadora, y proponer que la probabilidad es cercana al \(50\%\), esto es \(theta_t\) cercano a \(0\), y que estos valores son estables para los distintos valores de \(\theta_t\). Dado que \(\theta_t \sim N(\mu_{\theta}, \sigma^2_{\theta})\) y \(\mu_{\theta} \sim N(0, A)\), un alto valor de A indicaría que todos los \(\theta_t\) podrían estar centrados en valores muy distantes de 0, lo que no parece razonable.
Dado que \(\theta_t \sim N(\mu_{\theta}, \sigma^2_{\theta})\), con \(\sigma \sim U(0, B)\), para lograr variabilidad entre los distintos \(\theta_t\), parece razonble tener un alto valor de B.
Por estas razones se propone la utilización de un bajo valor para A y un alto valor para B. Sin embargo, a continuación se presentan los resultados de una simulación MCMC en Stan para el modelo, con distintos valores de A y B. Se tiene AB 1 (A = 1, B = 3), AB 2 (A = 1, B = 10000), AB 3 (A = 10000, B = 3) y AB 4 (A = 10000, B = 10000). Para cada ajuste se utilizaron 3 cadenas, de tamaño 30000, con 15000 de burning. Es de notar que AB 2 corresponde al modelo a priori más razonable según los argumentos expuestos anteriormente.
A continuación se muestran las densidades posteriori simuladas para cada valor de \(\theta_t\), comparando los 4 modelos.
modelo_stan_1 <- '
data {
int<lower=0> t; // numero de temporadas
int y[t]; // vector de canastas
int N[t]; // vector de intentos
real<lower=0> A;
real<lower=0> B;
}
parameters {
real theta[t];
real mu;
real<lower=0> sigma;
}
transformed parameters {
real<lower=0, upper=1> pi[t];
for (j in 1:t)
pi[j] = exp(theta[j]) / (1 + exp(theta[j]));
}
model {
y ~ binomial(N, pi);
theta ~ normal(mu, sigma);
mu ~ normal(0, A);
sigma ~ uniform(0, B);
}
generated quantities {
real log_lik[t];
int<lower=0> y_new;
real theta_new = normal_rng(mu, sigma);
real pi_new = exp(theta_new) / (1 + exp(theta_new));
y_new = binomial_rng(1450, pi_new);
for (i in 1:t) {
log_lik[i] = binomial_lpmf(y[i] | N[i], pi[i]);
}
}
'
datos_stan_1.1 <- list(t = length(datos$temp), y = datos$canastas, N = datos$intentos,
A = 1, B = 3)
datos_stan_1.2 <- list(t = length(datos$temp), y = datos$canastas, N = datos$intentos,
A = 1, B = 10000)
datos_stan_1.3 <- list(t = length(datos$temp), y = datos$canastas, N = datos$intentos,
A = 10000, B = 3)
datos_stan_1.4 <- list(t = length(datos$temp), y = datos$canastas, N = datos$intentos,
A = 10000, B = 10000)
modelo_compilado_1 <- stan_model(model_code = modelo_stan_1)
set.seed(2024)
ajuste_modelo_1.1 <- sampling(modelo_compilado_1, data = datos_stan_1.1, chains = 3,
iter = 30000, pars = c("theta", "pi", "lp__", "log_lik", "y_new", "pi_new"), warmup = 15000)
ajuste_modelo_1.2 <- sampling(modelo_compilado_1, data = datos_stan_1.2, chains = 3,
iter = 30000, pars = c("theta", "pi", "lp__", "log_lik", "y_new", "pi_new"), warmup = 15000)
ajuste_modelo_1.3 <- sampling(modelo_compilado_1, data = datos_stan_1.3, chains = 3,
iter = 30000, pars = c("theta", "pi", "lp__", "log_lik", "y_new", "pi_new"), warmup = 15000)
ajuste_modelo_1.4 <- sampling(modelo_compilado_1, data = datos_stan_1.4, chains = 3,
iter = 30000, pars = c("theta", "pi", "lp__", "log_lik", "y_new", "pi_new"), warmup = 15000)
theta.1 <- extract(ajuste_modelo_1.1, "theta")$theta
theta.2 <- extract(ajuste_modelo_1.2, "theta")$theta
theta.3 <- extract(ajuste_modelo_1.3, "theta")$theta
theta.4 <- extract(ajuste_modelo_1.4, "theta")$thetafor (i in 1:8) {
plot(density(theta.2[,i]), col = "blue3", lwd = 3, xlim = c(-0.31, 0), main= paste("Posteriori de theta ",i))
lines(density(theta.1[,i]), col = "red", lwd = 3)
lines(density(theta.3[,i]), col = "green4", lwd = 3)
lines(density(theta.4[,i]), col = "purple4", lwd = 3)
legend("topright", legend = c(paste("AB 1, mediana: ", round(median(theta.1[,i]),3)),
paste("AB 2, mediana: ", round(median(theta.2[,i]),3)),
paste("AB 3, mediana: ", round(median(theta.3[,i]),3)),
paste("AB 4, mediana: ", round(median(theta.4[,i]),3))),
col = c("red", "blue3", "green4", "purple4"), lwd = 2)
}
Es de destacar que, a pesar de tener valores de A y B muy distintos, no se observan grandes diferencias en las distribuciones a posteriori, ya sea por la forma de la densidad, por las estimaciones puntuales de la mediana. Es por esto, y por los argumentos expuestos al comienzo, que se presenta más razonable utilizar valores a priori para A y B de 1 y 10000, respectivamente. Además, se debe tener en cuenta que los valores de A no afectan en gran medida el agrupamiento entre los diferentes \(\theta_t\), en cambio, un alto valor de B sí implicaría que el nivel de agrupamiento de los \(\theta_t\) sea menor, por lo tanto se tendría mayor flexibilidad para el modelamiento.
En base a los valores de A y B definidos, a continuación se muestran las estimaciones de los parámetros \(\pi_t\) y sus intervalos de credibilidad al \(95\%\), para la distribución a posteriori en base al MCMC.
ajuste_modelo_1 <- ajuste_modelo_1.2
thetas_1 <- theta.2
plot(ajuste_modelo_1, pars = "pi")ci_level: 0.8 (80% intervals)outer_level: 0.95 (95% intervals)
Se observa una cierta tendencia de cambio para los valores de \(\pi_t\) entre las diferentes temporadas, aunque el cambio no es demasiado grande.
A continuación se muestran gráficos de traceplot y autocorrelación para mostrar que los parámetros convergen bien en el MCMC.
traceplot(ajuste_modelo_1, pars = c("theta"))
acf(thetas_1[,1],main=expression(theta[1999]))
acf(thetas_1[,2],main=expression(theta[2000]))
acf(thetas_1[,3],main=expression(theta[2001]))
acf(thetas_1[,4],main=expression(theta[2002]))
acf(thetas_1[,5],main=expression(theta[2003]))
acf(thetas_1[,6],main=expression(theta[2004]))
acf(thetas_1[,7],main=expression(theta[2005]))
acf(thetas_1[,8],main=expression(theta[2006]))
A continuación, se construye el modelo de agrupamiento total, es decir, que todos los \(\theta_t\) son idénticos. Al igual que para la pregunta 1, el modelamiento se realiza a partir del programa Stan, específicamente con el paquete stanr, a partir de 3 cadenas de tamaño 30000, con 15000 de burning.
modelo_stan_2 <- '
data {
int<lower=0> t; // numero de temporadas
int y[t]; // vector de canastas
int N[t]; // vector de intentos
real<lower=0> A;
real<lower=0> B;
}
parameters {
real theta;
real mu;
real<lower=0> sigma;
}
transformed parameters {
real<lower=0, upper=1> pi;
pi = exp(theta) / (1 + exp(theta));
}
model {
y ~ binomial(N, pi);
theta ~ normal(mu, sigma);
mu ~ normal(0, A);
sigma ~ uniform(0, B);
}
generated quantities {
real log_lik[t];
for (i in 1:t) {
log_lik[i] = binomial_lpmf(y[i] | N[i], pi);
}
}
'
datos_stan_2 <- list(t = length(datos$temp), y = datos$canastas, N = datos$intentos,
A = 1, B = 10000)
modelo_compilado_2 <- stan_model(model_code = modelo_stan_2)
set.seed(2024)
ajuste_modelo_2 <- sampling(modelo_compilado_2, data = datos_stan_2, chains = 3, iter = 30000,
pars = c("theta", "pi", "lp__", "log_lik"), warmup = 15000)
theta_2 <- extract(ajuste_modelo_2, "theta")$thetaA continuación se presenta en un gráfico, una comparación entre las distribuciones a posteriori del \(\theta\) del nuevo modelo, y los distintos \(\theta_t\) para el primer modelo, de la pregunta 1.
colores <- c("red3", "green", "orange", "purple", "yellow", "pink", "cyan", "magenta")
plot(density(theta_2), col = "blue", lwd = 3, xlim = c(-0.31, 0),
main= paste("Comparación entre theta agrupado y desagrupados"))
for (i in 1:8) lines(density(theta.2[,i]), lty = 2, lwd = 2, col = colores[i])
legend("topright", legend = c(paste("theta, median: ", round(median(theta_2),3)),
paste("theta_1, median: ", round(median(theta.2[,1]),3)),
paste("theta_2, median: ", round(median(theta.2[,2]),3)),
paste("theta_3, median: ", round(median(theta.2[,3]),3)),
paste("theta_4, median: ", round(median(theta.2[,4]),3)),
paste("theta_5, median: ", round(median(theta.2[,5]),3)),
paste("theta_6, median: ", round(median(theta.2[,6]),3)),
paste("theta_7, median: ", round(median(theta.2[,7]),3)),
paste("theta_8, median: ", round(median(theta.2[,8]),3))),
col = c("blue", colores), lwd = 2, lty = c(1, rep(2, 8)))
Se observa que las distintas distribuciones no difieren demasiado en cuanto a la centralidad, dado que la distribución del \(\theta\) agrupado logra abarcar a las distribuciones de los \(\theta_t\) individuales. Sin embargo, si es de notar que el modelo desagrupado presenta una mayor variabilidad entre sus distintas estimaciones de parámetros.
A continuación se presentan las diferencias entre la mediana de \(\theta_t\) y la mediana de \(theta\) agrupado.
thetas_2 <- theta_2
res <- data_frame("i" = seq(1, 8),
"theta_i - theta (median)" = round(apply(theta.2, 2, median) - median(theta_2), 3))
res |> gt()| i | theta_i - theta (median) |
|---|---|
| 1 | 0.014 |
| 2 | 0.011 |
| 3 | 0.017 |
| 4 | -0.005 |
| 5 | -0.017 |
| 6 | -0.024 |
| 7 | -0.007 |
| 8 | 0.014 |
Como se observa, las diferencias entre estimaciones no son demasiado grandes. En base a estos dos gráficos anteriores, no parece haber evidencia contundente para escoger un modelo agrupado o desagrupado, sino que será necesario realizar algún tipo de comparación de modelos más preciso, como serían los criterios de información, que se realizarán en la pregunta 4.
A continuación se muestran gráficos de traceplot y autocorrelación para mostrar que los parámetros convergen bien en el MCMC.
traceplot(ajuste_modelo_2, pars = c("theta"))
acf(thetas_2,main=expression(theta[1999]))
Para el nuevo modelo, al igual que en la pregunta 1, es necesario discutir acerca de cuáles serían los mejores valores par C y D. Dado que \(\theta_t \sim N(\theta_{t-1}, \sigma^2_{\theta})\) y \(\theta_{1999} \sim N(0, C)\) el valor de \(\pi_{1999}\) está centrado en una probabilidad del \(50\%\), por lo que sería una elección conservadora escoger un bajo valor de C. En cambio, dado que \(\sigma^2_{\theta} \sim U(0,D)\), sería lo indicado escoger un alto valor de D para que no exista demasiada agrupación entre los distintos valores de \(\theta_t\).
Por estas razones se propone la utilización de un bajo valor para C y un alto valor para D. Sin embargo, a continuación se presentan los resultados de una simulación MCMC en Stan para el modelo, con distintos valores de C y D. Se tiene CD 1 (A = 3, B = 3), CD 2 (A = 3, B = 10000), CD 3 (A = 10000, B = 3) y CD 4 (A = 10000, B = 10000). Para cada ajuste se utilizaron 3 cadenas, de tamaño 30000, con 15000 de burning. Es de notar que CD 2 corresponde al modelo a priori más razonable según los argumentos expuestos anteriormente.
A continuación se muestran las densidades posteriori simuladas para cada valor de \(\theta_t\), comparando los 4 modelos.
modelo_stan_3 <- '
data {
int<lower=0> t; // numero de temporadas
int y[t]; // vector de canastas
int N[t]; // vector de intentos
real<lower=0> C;
real<lower=0> D;
}
parameters {
real theta[t];
real mu;
real<lower=0> sigma;
}
transformed parameters {
real<lower=0, upper=1> pi[t];
for (j in 1:t)
pi[j] = exp(theta[j]) / (1 + exp(theta[j]));
}
model {
y ~ binomial(N, pi);
sigma ~ uniform(0, D);
theta[1] ~ normal(0, C);
for (j in 2:t)
theta[j] ~ normal(theta[j-1], sigma);
}
generated quantities {
real log_lik[t];
int<lower=0> y_new;
real theta_new = normal_rng(theta[t], sigma);
real pi_new = exp(theta_new) / (1 + exp(theta_new));
y_new = binomial_rng(1450, pi_new);
for (i in 1:t) {
log_lik[i] = binomial_lpmf(y[i] | N[i], pi[i]);
}
}
'
datos_stan_3.1 <- list(t = length(datos$temp), y = datos$canastas, N = datos$intentos,
C = 3, D = 3)
datos_stan_3.2 <- list(t = length(datos$temp), y = datos$canastas, N = datos$intentos,
C = 3, D = 10000)
datos_stan_3.3 <- list(t = length(datos$temp), y = datos$canastas, N = datos$intentos,
C = 10000, D = 3)
datos_stan_3.4 <- list(t = length(datos$temp), y = datos$canastas, N = datos$intentos,
C = 10000, D = 10000)
modelo_compilado_3 <- stan_model(model_code = modelo_stan_3)
set.seed(2024)
ajuste_modelo_3.1 <- sampling(modelo_compilado_3, data = datos_stan_3.1, chains = 3, iter = 30000,
pars = c("theta", "pi", "lp__", "log_lik", "y_new", "pi_new"), warmup = 15000)
ajuste_modelo_3.2 <- sampling(modelo_compilado_3, data = datos_stan_3.2, chains = 3, iter = 30000,
pars = c("theta", "pi", "lp__", "log_lik", "y_new", "pi_new"), warmup = 15000)
ajuste_modelo_3.3 <- sampling(modelo_compilado_3, data = datos_stan_3.3, chains = 3, iter = 30000,
pars = c("theta", "pi", "lp__", "log_lik", "y_new", "pi_new"), warmup = 15000)
ajuste_modelo_3.4 <- sampling(modelo_compilado_3, data = datos_stan_3.4, chains = 3, iter = 30000,
pars = c("theta", "pi", "lp__", "log_lik", "y_new", "pi_new"), warmup = 15000)
theta.1 <- extract(ajuste_modelo_3.1, "theta")$theta
theta.2 <- extract(ajuste_modelo_3.2, "theta")$theta
theta.3 <- extract(ajuste_modelo_3.3, "theta")$theta
theta.4 <- extract(ajuste_modelo_3.4, "theta")$thetafor (i in 1:8) {
plot(density(theta.4[,i]), col = "purple4", lwd = 3, xlim = c(-0.31, 0), main= paste("Posteriori de theta ",i), ylim = c(0, 15))
lines(density(theta.1[,i]), col = "red", lwd = 3)
lines(density(theta.2[,i]), col = "blue3", lwd = 3)
lines(density(theta.3[,i]), col = "green4", lwd = 3)
legend("topright", legend = c(paste("CD 1, mediana: ", round(median(theta.1[,i]),3)),
paste("CD 2, mediana: ", round(median(theta.2[,i]),3)),
paste("CD 3, mediana: ", round(median(theta.3[,i]),3)),
paste("CD 4, mediana: ", round(median(theta.4[,i]),3))),
col = c("red", "blue3", "green4", "purple4"), lwd = 2)
}
Es de destacar que, a pesar de tener valores de C y D muy distintos, no se observan grandes diferencias en las distribuciones a posteriori, ya sea por la forma de la densidad, por las estimaciones puntuales de la mediana. Es por esto, y por los argumentos expuestos al comienzo, que se presenta más razonable utilizar valores a priori para C y D de 3 y 10000, respectivamente.
En base a los valores de C y D definidos, a continuación se muestran las estimaciones de los parámetros \(\pi_t\) y sus intervalos de credibilidad al \(95\%\), para la distribución a posteriori en base al MCMC.
ajuste_modelo_3 <- ajuste_modelo_3.2
thetas_3 <- theta.2
plot(ajuste_modelo_3, pars = "pi")ci_level: 0.8 (80% intervals)outer_level: 0.95 (95% intervals)
Se observa una cierta tendencia de cambio para los valores de \(\pi_t\) entre las diferentes temporadas. Sin embargo, los cambios no son tan marcados al comparar con el modelo 1 de la pregunta 1, como habría de esperar dado la construcción de este modelo en que existe una dependencia temporal entre los distintos \(\theta_t\).
A continuación se muestran gráficos de traceplot y autocorrelación para mostrar que los parámetros convergen bien en el MCMC.
traceplot(ajuste_modelo_3, pars = c("theta"))
acf(thetas_3[,1],main=expression(theta[1999]))
acf(thetas_3[,2],main=expression(theta[2000]))
acf(thetas_3[,3],main=expression(theta[2001]))
acf(thetas_3[,4],main=expression(theta[2002]))
acf(thetas_3[,5],main=expression(theta[2003]))
acf(thetas_3[,6],main=expression(theta[2004]))
acf(thetas_3[,7],main=expression(theta[2005]))
acf(thetas_3[,8],main=expression(theta[2006]))
Para comparar las estimaciones puntuales de \(\theta_t\), se presentan a continuación gráficos que incluyen las densidades de las simulaciones a posteriori de los parámetros de \(\theta_t\) para cada modelo.
for (i in 1:8) {
plot(density(thetas_2), col = "blue3",lty = 2, lwd = 3, xlim = c(-0.31, 0), main= paste("Posteriori de theta ",i))
lines(density(thetas_1[,i]), col = "red", lwd = 3)
lines(density(thetas_3[,i]), col = "green4", lwd = 3)
legend("topright", legend = c(paste("model 1, mediana: ", round(median(thetas_1[,i]),3)),
paste("model 2, mediana: ", round(median(thetas_2),3)),
paste("model 3, mediana: ", round(median(thetas_3[,i]),3))),
col = c("red", "blue3", "green4"), lwd = 2, lty = c(1,2,1))
}
Al menos para los modelos 1 y 3, no se observan grandes diferencias entre las estimaciones puntuales de \(\theta_t\), ni tampoco en presentan grandes diferencias de variabilidad. Y, al comparar con el modelo 2 (que corresponde al modelo agrupado), tampoco existen tan grandes diferencias, teniendo en cuenta la simplicidad del modelo 2.
A continuación se presenta una tabla en que se muestran las diferencias entre medianas de los modelos 1 y 2, 2 y 3, y 1 y 3, respectivamente.
temporada <- seq(1999, 2006)
diff12 <- round(apply(thetas_1, MARGIN = 2, median) - median(thetas_2), 4)
diff23 <- round(median(thetas_2) - apply(thetas_3, MARGIN = 2, median), 4)
diff13 <- round(apply(thetas_1, MARGIN = 2, median) - apply(thetas_3, MARGIN = 2, median), 4)
cbind(temporada, diff12, diff23, diff13) |> as.data.frame() |> gt()| temporada | diff12 | diff23 | diff13 |
|---|---|---|---|
| 1999 | 0.0140 | -0.0314 | -0.0174 |
| 2000 | 0.0109 | -0.0261 | -0.0152 |
| 2001 | 0.0175 | -0.0206 | -0.0032 |
| 2002 | -0.0046 | 0.0050 | 0.0003 |
| 2003 | -0.0167 | 0.0253 | 0.0086 |
| 2004 | -0.0237 | 0.0323 | 0.0086 |
| 2005 | -0.0067 | 0.0136 | 0.0069 |
| 2006 | 0.0139 | -0.0066 | 0.0073 |
En general, se observa que las estimaciones puntuales difieren bastante poco, en especial al comparar el modelo 1 con el 3. Por lo cual no parece existir gran diferencia en este aspecto entre escoger un modelo u otro.
A continuación se presentan dos gráficos que muestran la distribución de las predicciones de los valores de \(y_{2007}\) y \(\pi_{2007}\) para la temporada siguiente a la última de la que se tiene registro en este set de datos. Las estimaciones se realizaron desde el mismo modelo Stan al realizar MCMC, y con un valor asumido de intentos de canasta \(N_{2007} = 1450\). Se muestran, además, la mediana de la predicción y los intervalos de credibilidad centrados al \(95\%\).
y_new1 <- extract(ajuste_modelo_1, "y_new")$y_new
y_new3 <- extract(ajuste_modelo_3, "y_new")$y_new
quantiles_y_new1 <- quantile(y_new1, c(0.025, 0.975))
quantiles_y_new3 <- quantile(y_new3, c(0.025, 0.975))
plot(density(y_new1), col = "blue3", lwd = 3, main= paste("Predicción de y_1999"))
lines(density(y_new3), col = "red", lwd = 3)
abline(v = quantiles_y_new1, col = "blue3", lty = 2)
abline(v = quantiles_y_new3, col = "red", lty = 2)
legend("topright", legend = c(paste("model 1, mediana: ", round(median(y_new1),3)),
paste("model 3, mediana: ", round(median(y_new3),3))),
col = c("blue3", "red"), lwd = 2)
pi_new1 <- extract(ajuste_modelo_1, "pi_new")$pi_new
pi_new3 <- extract(ajuste_modelo_3, "pi_new")$pi_new
quantiles_pi_new1 <- quantile(pi_new1, c(0.025, 0.975))
quantiles_pi_new3 <- quantile(pi_new3, c(0.025, 0.975))
plot(density(pi_new1), col = "blue3", lwd = 3, main= paste("Predicción de pi_1999"))
lines(density(pi_new3), col = "red", lwd = 3)
abline(v = quantiles_pi_new1, col = "blue3", lty = 2)
abline(v = quantiles_pi_new3, col = "red", lty = 2)
legend("topright", legend = c(paste("model 1, mediana: ", round(median(pi_new1),3)),
paste("model 3, mediana: ", round(median(pi_new3),3))),
col = c("blue3", "red"), lwd = 2)
Se observa que las predicciones puntuales, a partir de la mediana, son muy similares entre ambos modelos. Sin embargo, es de notar que los intervalos de credibilidad para el modelo 1 son claramente más acotados, por lo que su predicción parece ser más precisa que para el modelo 3.
Para finalizar el análisis de los distintos modelo, a continuación se presenta una tabla que muestra los criterios de información calculados BIC, AIC, DIC y WAIC (calculado a partir del paquete loo, que utiliza leave-one-out crossvalidation para el cálculo).
# BIC
pi_1 <- colMeans(extract(ajuste_modelo_1, "pi")$pi)
BIC_1 <- 2 * sum(dbinom(datos$canastas, datos$intentos, pi_1, log = T)) - (8+2) * log(length(datos$canastas))
pi_2 <- mean(extract(ajuste_modelo_2, "pi")$pi)
BIC_2 <- 2 * sum(dbinom(datos$canastas, datos$intentos, pi_2, log = T)) - (1+2) * log(length(datos$canastas))
pi_3 <- colMeans(extract(ajuste_modelo_3, "pi")$pi)
BIC_3 <- 2 * sum(dbinom(datos$canastas, datos$intentos, pi_3, log = T)) - (8+1) * log(length(datos$canastas))
BIC <- c(BIC_1, BIC_2, BIC_3)
# AIC
pis_1 <- extract(ajuste_modelo_1, "pi")$pi
pi_1 <- colMeans(pis_1)
AIC_1 <- 2 * (8+2) - 2 * sum(dbinom(datos$canastas, datos$intentos, pi_1, log = T))
pis_2 <- extract(ajuste_modelo_2, "pi")$pi
pi_2 <- mean(pis_2)
AIC_2 <- 2 * (1+2) - 2 * sum(dbinom(datos$canastas, datos$intentos, pi_2, log = T))
pis_3 <- extract(ajuste_modelo_3, "pi")$pi
pi_3 <- colMeans(pis_3)
AIC_3 <- 2 * (8+2) - 2 * sum(dbinom(datos$canastas, datos$intentos, pi_3, log = T))
AIC <- c(AIC_1, AIC_2, AIC_3)
# DIC
log_p_bay_1 <- sum(dbinom(datos$canastas, datos$intentos, pi_1, log = T))
log_p_bay_2 <- sum(dbinom(datos$canastas, datos$intentos, pi_2, log = T))
log_p_bay_3 <- sum(dbinom(datos$canastas, datos$intentos, pi_3, log = T))
N <- length(pis_2)
suma <- matrix(0,nrow=3,ncol=N)
for (i in 1:N) {
suma[1,i] <- sum(dbinom(datos$canastas, datos$intentos, pis_1[i,], log = T))
suma[2,i] <- sum(dbinom(datos$canastas, datos$intentos, pis_2[i], log = T))
suma[3,i] <- sum(dbinom(datos$canastas, datos$intentos, pis_3[i,], log = T))
}
suma <- rowMeans(suma)
p_dic_1 <- 2 * (log_p_bay_1 - suma[1])
p_dic_2 <- 2 * (log_p_bay_2 - suma[2])
p_dic_3 <- 2 * (log_p_bay_3 - suma[3])
DIC_1 <- -2 * log_p_bay_1 + 2 * p_dic_1
DIC_2 <- -2 * log_p_bay_2 + 2 * p_dic_2
DIC_3 <- -2 * log_p_bay_3 + 2 * p_dic_3
DIC <- c(DIC_1, DIC_2, DIC_3)
# WAIC
waic.1 <- loo(ajuste_modelo_1)
waic.2 <- loo(ajuste_modelo_2)
waic.3 <- loo(ajuste_modelo_3)
WAIC <- c(waic.1$estimates[3,1], waic.2$estimates[3,1], waic.3$estimates[3,1])comparaciones <- data.frame("Modelo" = c(1, 2, 3),
"BIC" = BIC,
"AIC" = AIC,
"DIC" = DIC,
"WAIC" = WAIC)
comparaciones |> gt() |> tab_header(title = "Comparación de modelos en base a criterios de información")| Comparación de modelos en base a criterios de información | ||||
| Modelo | BIC | AIC | DIC | WAIC |
|---|---|---|---|---|
| 1 | -85.95495 | 85.16053 | 72.7870 | 72.54265 |
| 2 | -75.77307 | 75.53474 | 71.5028 | 71.58076 |
| 3 | -83.00260 | 84.28762 | 71.9103 | 70.93656 |
Para los criterios BIC, AIC y DIC se tiene que el modelo 2 es el que mejor se ajusta a los datos. Y según el criterio WAIC correspone el modelo 3. Sin embargo, se deben realizar aclaraciones. BIC y AIC presentan al modelo 2 como mejor modelo con bastante diferencia; esto se debe a su fuerte manera de penalizar la cantidad de parámetros utilizadas en el modelo, dado que el modelo 2 utiliza un solo valor de \(\theta_t\), en vez de los 8 que utilizan los demás modelos. Además, estos dos criterios no corresponden a criterios estrictamente bayesianos. Luego, las diferencias de DIC y WAIC para modelos 2 y 3 corresponden a 0.4075 y 0.6442, respectivamente. La diferencia entre los DIC es tan baja que no es suficiente para escoger uno de los dos modelos. Y la diferencias de WAIC para ambos modelos también es baja, pero debe considerarse con mayor atención dado que el WAIC es un criterio completamente bayesiano y, por lo tanto, más adecuado para comparar modelos. Entonces, se escoge al modelo 3 como el más adecuado, aunque con precaución.