El siguiente pots tiene como objetivo principal el uso del paquete brms para la modelización bayesiana multinivel avanzada, en este caso utilizaremos varias bases de datos para mostrar su potencial al momento de realizar estimaciones y ajustes de los modelos.
En este caso utilizaremos los datos de Markus Gesmann el cual predice el crecimiento de los pagos acumulados por siniestros de seguros a lo largo del tiempo. A modo de demostración utilizaremos una versión ligeramente simplificada de su modelo el cual es el siguiente:
\[ \begin{equation} cumAY_{, dev} \sim N(\mu _{AY, dev},\sigma) \end{equation} \]
\[ \begin{equation*} \mu _{AY, dev}=ult_{AY}\left(1-exp\left(-\left(\dfrac{dev}{\theta}\right)^\omega\right)\right) \end{equation*} \]
Aquí los pagos acumulados del seguro crecerán con el tiempo y se modela esta dependencia utilizando la variable dev , \(ult_{AY}\) es la perdida ultima de siniestralidad (por estimar) de cada año, constituye un parámetro no lineal junto con los parámetro \(\theta, \omega\) que son los responsables del crecimiento de la perdida acumulada y que se supone que son los mismos todos los años.
url <- paste0("https://raw.githubusercontent.com/mages/",
"diesunddas/master/Data/ClarkTriangle.csv")
Seguros <- read.csv(url)
head(Seguros) AY dev cum premium
1 1991 6 357.848 10000
2 1991 18 1124.788 10000
3 1991 30 1735.330 10000
4 1991 42 2182.708 10000
5 1991 54 2745.596 10000
6 1991 66 3319.994 10000Ahora utilizamos brms para el modelo no lineal propuesto
library(brms)
library(lme4)
nlform <- bf(cum ~ ult * (1 - exp(-(dev / theta)^omega)),
ult ~ 1 + (1|AY), omega ~ 1, theta ~ 1, nl = TRUE)
lprior <- c(prior(normal(5000, 1000), nlpar = "ult"),
prior(normal(1, 2), nlpar = "omega"),
prior(normal(45, 10), nlpar = "theta"))
Modelo1 <- brm(formula = nlform, data = Seguros, family = gaussian(),
prior = lprior, control = list(adapt_delta = 0.9))summary(Modelo1) Family: gaussian
Links: mu = identity; sigma = identity
Formula: cum ~ ult * (1 - exp(-(dev/theta)^omega))
ult ~ 1 + (1 | AY)
omega ~ 1
theta ~ 1
Data: Seguros (Number of observations: 55)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~AY (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(ult_Intercept) 741.21 232.65 422.48 1307.30 1.00 1255 1905
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
ult_Intercept 5293.08 287.76 4693.22 5871.90 1.00 988 1425
omega_Intercept 1.34 0.05 1.24 1.44 1.00 2443 2808
theta_Intercept 46.19 2.13 42.42 50.81 1.00 2524 2365
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 139.61 15.43 113.13 173.90 1.00 2656 2224
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).conditional_effects(Modelo1)
Veamos que en la figura 1 se puede visualizar la pérdida de seguro acumulada a lo largo del tiempo por separado para cada año
conditions <- data.frame(AY = unique(loss$AY))
rownames(conditions) <- unique(loss$AY)
me_year <- conditional_effects(Modelo1, conditions = conditions,
re_formula = NULL, method = "predict")
plot(me_year, ncol = 5, points = TRUE)
La figura 2 muestra mejor esa variación en las perdidas acumuladas entro los años de accidente , esto debido a que los desastres no suelen ocurrir en determinado años.
En el modelo anterior se considero que omega y delta son constantes a o largo de los años lo que no puede ser necesariamente cierto, esto se puede verificar ajustando los intercepto variables de los parámetros no lineales y estimando la correlación a nivel de grupo utilizando la sintaxis ID
library(brms)
nlform2 <- bf(cum ~ ult * (1 - exp(-(dev / theta)^omega)),
ult ~ 1 + (1|ID1|AY), omega ~ 1 + (1|ID1|AY),
theta ~ 1 + (1|ID1|AY), nl = TRUE)
Modelo2 <- update(Modelo1, formula = nlform2,
control = list(adapt_delta = 0.90))Realicemos validación cruzada para comparar el ajuste de los modelos
# Mostrar los resultados
print(results) LOOIC SE
Modelo1 715.960350 19.5546384
Modelo2 717.820450 18.6019357
Modelo1-Modelo2 -1.860101 0.9527027