Regresión poisson bayesiano
Rojas Rojas Henry
2023-09-24
Práctica calificada N°03
El ejercicio propuesto, La base de datos CASOSCOVID contiene información de los casos de contagio por covid-19 en el transcurso de 15 días, el detalle de la base se presenta a continuación:
Día: Número de día Casos: Casos de COVID por día Miembros hogar: Cantidad promedio del número de miembros por hogar. Animales: Número promedio de mascotas por hogar. *Integrantes vacunados: Cantidad promedio de miembros por hogar que recibieron al menos la primera vacuna contra la COVID.
#Cargar los datos
datos <- read.table("CASOSCOVID.txt", header = T, sep = "\t")
attach(datos)
head(datos,5) #Observar los 5 datos de la base de datos## DIA CASOS MIEMBROS_HOGAR ANIMALES INTEGRANTES_VACUNADOS
## 1 dia1 150 2 3 1
## 2 dia2 180 7 0 3
## 3 dia3 220 8 3 5
## 4 dia4 130 3 0 2
## 5 dia5 145 5 3 3De acuerdo a los datos, Se pide:
Regresión poisson clásico
- Estimar el modelo de regresión poisson clásico e interpretar los parámetros.
# Ajusta el modelo de regresión de Poisson
modelo_poisson <- glm(CASOS ~ MIEMBROS_HOGAR+ANIMALES+INTEGRANTES_VACUNADOS, data = datos, family = "poisson")
# Obtén un resumen de los resultados del modelo
summary(modelo_poisson)##
## Call:
## glm(formula = CASOS ~ MIEMBROS_HOGAR + ANIMALES + INTEGRANTES_VACUNADOS,
## family = "poisson", data = datos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.38811 0.07439 58.987 < 2e-16 ***
## MIEMBROS_HOGAR 0.10212 0.01716 5.952 2.64e-09 ***
## ANIMALES 0.06324 0.01996 3.168 0.00153 **
## INTEGRANTES_VACUNADOS -0.01646 0.02983 -0.552 0.58104
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 178.912 on 14 degrees of freedom
## Residual deviance: 83.578 on 11 degrees of freedom
## AIC: 193.62
##
## Number of Fisher Scoring iterations: 4Nota. Se observa que los coeficiente del modelo de regresión poisson clasico, nos expresan los siguiente: B(Miembros_Hogar)=0.10212 : Nos da ha entender que, si la cantidad promedio de miembros del hogar aumenta, y el número de promedio de de mascotas por hogar y cantidad promedio de miembros por hogar que recibieron al menos la primera vacuna contra la COVID se mantienen constantes, los casos de covid por día aumentaría. B(Animales)=0.06324 : Nos da ha entender que, si el número de promedio de de mascotas por hogar aumenta, y la cantidad promedio de miembros del hogar y la cantidad promedio de miembros por hogar que recibieron al menos la primera vacuna contra la COVID se mantienen constantes, los casos de covid por día aumentaría. *B(Integrantes_vacunados)=-0.01646 : Nos da ha entender que, si la cantidad promedio de miembros por hogar que recibieron al menos la primera vacuna contra la COVID aumenta, y el número de promedio de de mascotas por hogar y la cantidad promedio de miembros del hogar se mantienen constantes, los casos de covid por día disminuiria.
- Evaluar los supuestos del modelo.
#Evaluación de supuestos
# Verifica la bondad de ajuste del modelo
#Bondad de ajuste del modelo
plot(modelo_poisson, which = 1)
Nota. Se observa que los datos se ajustan a la recta, lo que nos da ha entender que el modelo es bueno.
library(AER)## Loading required package: car## Loading required package: carData## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Loading required package: survival#Examina la sobredispersión o subdispersión:
dispersion_test <- dispersiontest(modelo_poisson)
dispersion_test##
## Overdispersion test
##
## data: modelo_poisson
## z = 2.9092, p-value = 0.001812
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion
## 5.409008Nota. De acuerdo a la dispersión de los datos, se observa que, es significativo, p-value menor al 5%, dando a entender que la verdadera dispersión es mayor que 1.
Regresión poisson bayesiano
Considerando los datos del anterior ejercicio:
- Estimar el modelo de regresión bayesiano.
#BAYESIANO
# install.packages("brms")
# Carga el paquete brms
library(brms)## Loading required package: Rcpp## Loading 'brms' package (version 2.20.3). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:survival':
##
## kidney## The following object is masked from 'package:stats':
##
## ar# Utiliza la función brm() para ajustar el modelo
model <- brm(CASOS ~ MIEMBROS_HOGAR+ANIMALES+INTEGRANTES_VACUNADOS, family = poisson(), data = datos)## Compiling Stan program...## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 7.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.73 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.095 seconds (Warm-up)
## Chain 1: 0.075 seconds (Sampling)
## Chain 1: 0.17 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.145 seconds (Warm-up)
## Chain 2: 0.082 seconds (Sampling)
## Chain 2: 0.227 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.121 seconds (Warm-up)
## Chain 3: 0.1 seconds (Sampling)
## Chain 3: 0.221 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.097 seconds (Warm-up)
## Chain 4: 0.124 seconds (Sampling)
## Chain 4: 0.221 seconds (Total)
## Chain 4:# Resumen del modelo
summary(model)## Family: poisson
## Links: mu = log
## Formula: CASOS ~ MIEMBROS_HOGAR + ANIMALES + INTEGRANTES_VACUNADOS
## Data: datos (Number of observations: 15)
## 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
## Intercept 4.39 0.07 4.24 4.53 1.00 3321
## MIEMBROS_HOGAR 0.10 0.02 0.07 0.14 1.01 1407
## ANIMALES 0.06 0.02 0.02 0.10 1.00 2792
## INTEGRANTES_VACUNADOS -0.02 0.03 -0.07 0.04 1.01 1538
## Tail_ESS
## Intercept 2786
## MIEMBROS_HOGAR 1849
## ANIMALES 2578
## INTEGRANTES_VACUNADOS 2360
##
## 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).Nota. Se observa que los coeficiente del modelo de regresión poisson bayesiano.
- Comparar los parámetros del modelo poisson clásico y bayesiano.
library(insight)
#Para el clasico
clas <- get_parameters(modelo_poisson)
clasico<-data.frame(clas)
#Para bayesiana
post <- get_parameters(model)
Bayess<-data.frame(purrr::map_dbl(post,median))
#En una sola tabla
Comparacion<-data.frame(clasico,Bayess)
Comparacion## Parameter Estimate
## b_Intercept (Intercept) 4.38810964
## b_MIEMBROS_HOGAR MIEMBROS_HOGAR 0.10211804
## b_ANIMALES ANIMALES 0.06323648
## b_INTEGRANTES_VACUNADOS INTEGRANTES_VACUNADOS -0.01646320
## purrr..map_dbl.post..median.
## b_Intercept 4.38545222
## b_MIEMBROS_HOGAR 0.10249796
## b_ANIMALES 0.06390041
## b_INTEGRANTES_VACUNADOS -0.01741101Nota. Se observa que los valores se asemejan, entre los dos metódos, tienen la mista tendencia de sus estimadores.
# Gráficos y diagnósticos
plot(model)
Nota. Se observa que, los datos son semejantes, uno ensima el otro, buen ajuste.
- Estimar el número de casos de COVID en un día, para los siguientes valores:
- Miembros del hogar = 10
- Animales = 3
- Integrantes vacunados = 6
# Predicciones
new_data <- data.frame(MIEMBROS_HOGAR=10,ANIMALES=3,INTEGRANTES_VACUNADOS=6)
predictions <- predict(model, newdata = new_data, summary = TRUE)
# Muestra las predicciones
print(predictions)## Estimate Est.Error Q2.5 Q97.5
## [1,] 244.7217 22.13583 203 290Nota. Considerando los valores expuestos, se espera una cantidad de 246 casos de covid 19 aproximadamente.
- A MODO PRÁCTICA