Modelización Bayesiana Multinivel Avanzada con R con el paquete brms

Author

Eduardo Canales

Introducción

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.

Pesca de Peces

Los siguientes datos trata del número de peces capturados por varios grupos de personas, los datos se obtuvieron del sitio web (https: //stats.idre.ucla.edu/stata/dae/zero-inflated-poisson-regression) de la UCLA. Una descripción breve de los datos: "Los biólogos de vida salvaje desean modelar cuantos peces son capturados por los pescadores en un parque estatal y se realizan las siguientes preguntas a los visistantes; Cuánto tiempo se quedaron, Cuántas personas había en el grupo, si había niños en el grupo y cuántos peces se pescaron. Algunos visitantes no pescan, pero no hay datos sobre si una persona pescó o no, algunos visitantes si pescaron pero no capturan ningún pez, por lo que hay un exceso de ceros en los datos debido a las personas que nos pescaba"

La base de datos Peces contiene la siguientes variables:

nofish: Esta variable representa el número de peces que no fueron capturados en un día de pesca en particular.
livebait: Es una variable binaria que indica si se utilizó cebo vivo durante la pesca (1) o no (0).
camper: Esta variable muestra si el pescador era un campista (yes) o no (no).
persons: Representa el número de personas presentes en la pesca.
child: Es una variable binaria que indica si había niños presentes durante la pesca (1) o no (0).
xb: Esta variable corresponde a un predictor continuo o covariable en un modelo estadístico.
zg: Esta variable es un predictor continuo o covariable en un modelo estadístico.
count: Representa el número de peces capturados en un día de pesca en particular.

Peces<- read.csv("http://stats.idre.ucla.edu/stat/data/fish.csv")
Peces$camper<-factor(Peces$camper,labels = c("no","yes"))
head(Peces,n=5)

  nofish livebait camper persons child         xb         zg count
1      1        0     no       1     0 -0.8963146  3.0504048     0
2      0        1    yes       1     0 -0.5583450  1.7461489     0
3      0        1     no       1     0 -0.4017310  0.2799389     0
4      0        1    yes       2     1 -0.9562981 -0.6015257     0
5      0        1     no       1     0  0.4368910  0.5277091     1

Para este caso se elige como predictores el numero de personas por grupo, el numero de niños, así como si el grupo esta formado por campistas o no. Dado que muchos grupos puede que no pesquen nada simplemente no lo intenten, por lo cual se ajusta un Modelo de Poisson de inflación cero.

library("brms")
Mod1<-brm(count~persons+child+camper,
          data = Peces,
          family = zero_inflated_poisson("log"))

summary(Mod1)

 Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: count ~ persons + child + camper 
   Data: Peces (Number of observations: 250) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -1.01      0.18    -1.37    -0.66 1.00     2332     2504
persons       0.87      0.05     0.79     0.96 1.00     2459     2723
child        -1.36      0.09    -1.56    -1.18 1.00     2467     2432
camperyes     0.80      0.09     0.62     0.99 1.00     2934     2357

Family Specific Parameters: 
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi     0.41      0.05     0.31     0.49 1.00     2740     2627

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).

Figure 1: Gráficos de efectos condicionales del modelo Mod1.

Figure 2: Gráficos de efectos condicionales del modelo Mod1.

Figure 3: Gráficos de efectos condicionales del modelo Mod1.

Ahora, se intenta predecir adicionalmente la probabilidad de inflación cero en función del número de hijos

library("brms")
Mod2<-brm(bf(count~persons+child+camper, zi~child),
          data = Peces,
          family = zero_inflated_poisson())

summary(Mod2)

 Family: zero_inflated_poisson 
  Links: mu = log; zi = logit 
Formula: count ~ persons + child + camper 
         zi ~ child
   Data: Peces (Number of observations: 250) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept       -1.08      0.18    -1.44    -0.72 1.00     3075     2843
zi_Intercept    -0.96      0.26    -1.48    -0.47 1.00     3815     2509
persons          0.89      0.05     0.80     0.99 1.00     3137     2892
child           -1.18      0.10    -1.37    -0.99 1.00     2782     2878
camperyes        0.77      0.09     0.60     0.96 1.00     4050     3185
zi_child         1.22      0.28     0.69     1.77 1.00     3536     2874

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).

Realizamos una comparación del ajuste del modelo mediante validación cruzada leavel-one-out-al

# Mostrar los resultados
print(results)

              LOOIC          SE
Mod1      1632.2478 357.9252602
Mod2      1614.3993 358.0806199
Mod1-Mod2   17.8485  -0.1553597

Alquiler de viviendas

Los siguientes datos para analizar son los alquileres de viviendas en Múnich a partir de 1999, los cuales contienen información sobre 3000 apartamentos aproximadamente.

La base de datos rent99 se enucentra dentro del paquete gamlss.data donde contiene la siguientes variables:

rent: Alquiler absoluto
rentsqm Alquiler por metro cuadrado
area El tamaño del apartamento
yearc Año de construcción .
district El distrito de Múnich
location Ubicación
bath Baños del apartamentos.

kitchen Cocina.

data("rent99",package = "gamlss.data")
head(rent99,n=5)

      rent  rentsqm area yearc location bath kitchen cheating district
1 109.9487 4.228797   26  1918        2    0       0        0      916
2 243.2820 8.688646   28  1918        2    0       0        1      813
3 261.6410 8.721369   30  1918        1    0       0        1      611
4 106.4103 3.547009   30  1918        2    0       0        0     2025
5 133.3846 4.446154   30  1918        2    0       0        1      561

El objetivo en este caso es predecir el alquiler por metro cuadrado con el tamaño del apartamento y el año de construcción , teniendo en cuenta el distrito de Múnich

library("brms")
Mod3<- brm(rentsqm ~ t2(area, yearc) + (1|district), data = rent99,
                 chains = 2, cores = 2)

summary(Mod3)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rentsqm ~ t2(area, yearc) + (1 | district) 
   Data: rent99 (Number of observations: 3082) 
  Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 2000

Smooth Terms: 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sds(t2areayearc_1)     3.91      1.61     1.46     7.64 1.00     1610     1202
sds(t2areayearc_2)     4.13      1.87     1.36     8.69 1.00     1318     1316
sds(t2areayearc_3)     6.35      2.16     3.22    11.71 1.00     1561     1325

Group-Level Effects: 
~district (Number of levels: 336) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.59      0.06     0.47     0.71 1.00      679      982

Population-Level Effects: 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept         7.80      0.11     7.59     8.01 1.00     2051     1604
t2areayearc_1    -1.01      0.08    -1.16    -0.85 1.00     2557     1701
t2areayearc_2     0.74      0.16     0.42     1.06 1.00     1887     1429
t2areayearc_3    -0.11      0.15    -0.39     0.20 1.00     2088     1415

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.96      0.03     1.91     2.01 1.00     2847     1421

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(Mod3)

En el ejemplo anterior solo se considero que media de la distribución solo varía según la zona y el año, en consecuencia se ajusta splines y efectos de distrito tanto para la localización como el parámetro de escala, que denominamos sigma en los modelos gausiano

library("brms")
bform <- bf(rentsqm ~ t2(area, yearc) + (1|ID1|district),
            sigma ~ t2(area, yearc) + (1|ID1|district))
Mod4 <- brm(bform, data = rent99, chains = 2, cores = 2)

summary(Mod4)

 Family: gaussian 
  Links: mu = identity; sigma = log 
Formula: rentsqm ~ t2(area, yearc) + (1 | ID1 | district) 
         sigma ~ t2(area, yearc) + (1 | ID1 | district)
   Data: rent99 (Number of observations: 3082) 
  Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 2000

Smooth Terms: 
                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sds(t2areayearc_1)           4.08      1.60     1.60     7.73 1.00     1485
sds(t2areayearc_2)           4.36      2.00     1.52     8.98 1.00     1351
sds(t2areayearc_3)           5.83      2.10     2.74    10.64 1.00     1416
sds(sigma_t2areayearc_1)     0.54      0.49     0.01     1.81 1.00      925
sds(sigma_t2areayearc_2)     1.70      0.83     0.56     3.70 1.00      890
sds(sigma_t2areayearc_3)     0.75      0.60     0.03     2.26 1.00      787
                         Tail_ESS
sds(t2areayearc_1)           1281
sds(t2areayearc_2)           1218
sds(t2areayearc_3)           1350
sds(sigma_t2areayearc_1)     1230
sds(sigma_t2areayearc_2)      646
sds(sigma_t2areayearc_3)      886

Group-Level Effects: 
~district (Number of levels: 336) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.60      0.06     0.49     0.73 1.00
sd(sigma_Intercept)                0.11      0.02     0.06     0.15 1.00
cor(Intercept,sigma_Intercept)     0.73      0.17     0.36     0.99 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                       828     1108
sd(sigma_Intercept)                 758     1428
cor(Intercept,sigma_Intercept)      579      709

Population-Level Effects: 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept               7.75      0.11     7.53     7.97 1.00     1522     1355
sigma_Intercept         0.71      0.03     0.66     0.78 1.00     1261     1495
t2areayearc_1          -1.00      0.08    -1.17    -0.84 1.00     1840     1631
t2areayearc_2           0.77      0.15     0.47     1.07 1.00     1618     1524
t2areayearc_3          -0.09      0.16    -0.40     0.21 1.00     1701     1529
sigma_t2areayearc_1     0.03      0.02    -0.02     0.08 1.00     1786     1507
sigma_t2areayearc_2     0.06      0.04    -0.03     0.13 1.00     1507     1508
sigma_t2areayearc_3     0.01      0.04    -0.07     0.09 1.00     1780     1424

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).