Los modelos lineales son una de las herramientas estadísticas más utilizadas actualmente,los datos que contienen variables de agrupación o relación aparecen en muchas áreas, en este caso modelar este tipo de datos mediante modelos lineales puede ser adecuado aunque en ocasiones no se llega a muy buenos resultados, generalmente una buena y mejor alternativa es utilizar modelos mixtos.

#Cargar paquetes
library(nlme)
library(tidyverse)
library(ggpubr)

1 Descripción base de datos

datos <- read.csv2("Base.csv")
datos <- datos %>%
          mutate(Estatus = factor(Estatus, 
                           levels = c('En desarrollo','Desarrollado'))) %>%
          rename(genero = sexo) %>%
          mutate(year = as.factor(year))
head(datos)

1.1 Acerca de la base de datos

Los datos fueron obtenidos de la página del Banco Mundial, la base de datos cuenta con indicadores de desarrollo de 93 países (436 observaciones) sobre una variedad de temas que juntos construyen un retrato del progreso social en todo el mundo. Cubren medidas económicas de salud, diversas medidas de degradación ambiental.

1.2 Variables

  • País: País donde se tienen los indicadores

  • Estatus: Estatus del país
    • Desarrollado
    • En desarrollo
  • Género: Género para la esperanza de vida media.
    • Masculino
    • Femenino

  • Año: Año en el que se tienen los indicadores

  • Población: población total medida en millones

  • Area Forestal: porcentaje de area forestal del total de área del país.

  • Gasto en salud: porcentaje de gasto en salud con respecto al PIB.

  • Incidencia de VIH: Número de nuevos infectados con VIH entre habitantes no infectados de 15 a 49 años, expresadas por cada 100 personas no infectadas en el año anterior.

  • Emisiones Pm2.5:emisiones de material particulado 2.5 medido en microgramos por metro cúbico.

  • Enfermedades: es el porcentaje de personas de 30 años que moriría antes de su 70 cumpleaños de cualquier problema cardiovascular, cáncer, diabetes, o problemas respiratorios crónicos suponiendo que experimentaría las tasas de mortalidad actuales en todas las edades y que no morirá por ninguna otra causa (VIH/SIDA).

  • Esperanza de vida: Años que un recién nacido puede esperar vivir si los patrones de mortalidad por edades imperantes en el momento de su nacimiento siguieran siendo los mismos a lo largo de toda su vida.

2 Objetivo

  • Predecir la experanza de vida mediante modelos mixtos tomando como variable de agrupación el país, además observar la significancia de los efectos fijos y aleatorios mediante pruebas de hipótesis para finalmente escoger el mejor modelo.

  • Observar mediante análisis descrptivo posibles covariables para incluir en los modelos propuestos

3 Analisis descriptivo

3.1 Variable respuesta

datos %>% ggplot(aes(x=EsperanzaDeVida)) + 
  geom_histogram(color="black", fill="salmon",binwidth=2) +
  theme_bw()

Las mayores frecuencias para la esperanza de vida se encuentran entre los 70 y 80 años, sin despreciar que hay esperanza de vida muchísimo más bajos y logran ser menores a 40, la mayor esperanza de vida está entre 85 y 90 años. La distribución parece ser bimodal, esto puede ser explicado por otras variables.

3.2 Exploración de posibles covariables para el modelo

3.2.1 Pais

Para ilustrar el comportamiento de la esperanza de vida por país se seleccionan los paises de Angola, Argentina, Dinamarca, Ecuador, Francia, India, Italia, Sur África y Uruguay.

datos %>% filter(Pais %in% c("France","Ecuador","Uruguay",
                             "Italy","South Africa","Angola",
                             "Denmark","Argentina","India")) %>%
  ggplot(aes(x=year,y=EsperanzaDeVida,color=genero)) +geom_point() +
  theme_bw() +
  facet_wrap(~ Pais)

Es claro que el comportamiento por país es bastante diferente, además en cada país se presentan diferencias para hombres y mujeres, la variable país puede ser considerada como variable de agrupación.

3.2.2 Estatus

datos %>% ggplot(aes(x = Estatus,y = EsperanzaDeVida,fill=Estatus)) + geom_boxplot(alpha=0.55) + theme_bw()

datos %>% ggplot(aes(x =EsperanzaDeVida,fill=Estatus)) + geom_density(alpha=0.55) + theme_bw()

Es claro que el estatus es una variable que influye en la esperanza de vida, como se ve en el boxplot y en el gráfico de densidad, la media varia para cada categoría del estatus, es evidente que para los países que son desarrollados la esperanza de vida es mucho más alta pues su media se encuentra alrededor de 80 años, mientras que para los países que están en desarrollo los datos están más dispersos y la media puede estar entre 60 y 70 años.

3.2.3 Genero

datos %>% ggplot(aes(x = genero,y = EsperanzaDeVida,fill=genero)) + geom_boxplot(alpha=0.55) + theme_bw()

datos %>% ggplot(aes(x =EsperanzaDeVida,fill=genero)) + geom_density(alpha=0.55) + theme_bw()

Aunque la esperanza de vida para hombres y mujeres no se diferencia por mucho, se puede concluir que las mujeres esperan vivir más años que los hombres,en ambas densidades se ve bimodalidad y puede ser debida a otros factores.

3.2.4 Variables continuas

area <- datos %>% ggplot(aes(x=AreaForestalP, y=EsperanzaDeVida)) + geom_point(alpha = 0.6, size = 1) + theme_bw()
gasto <- datos %>% ggplot(aes(x=GastoSaludPIB, y=EsperanzaDeVida)) + geom_point(alpha = 0.6, size = 1)+ theme_bw()
vih <-  datos %>% ggplot(aes(x=IncidenciaVIHP, y=EsperanzaDeVida)) + geom_point(alpha = 0.6, size = 1)+ theme_bw()
pm2.5 <- datos %>% ggplot(aes(x=Pm2.5, y=EsperanzaDeVida)) + geom_point(alpha = 0.6, size = 1)+ theme_bw()
enfermedades <- datos %>% ggplot(aes(x=Enfermedades, y=EsperanzaDeVida)) + geom_point(alpha = 0.6, size = 1)+ theme_bw()
poblacion <- datos %>% ggplot(aes(x=pob, y=EsperanzaDeVida)) + geom_point(alpha = 0.6, size = 1)+ theme_bw()
ggarrange(area,gasto,vih,pm2.5,enfermedades,poblacion,
          ncol = 3, nrow = 2)

Para observar el comportamiento de cada variable continua se realiza un gráfico de dispersión contra la variable respuesta.

Como se puede evidenciar para los gráficos de dispersión en los cuales se grafican las variables continuas vs la esperanza de vida, algunas como el porcentaje de área forestal y el gasto en salud no parecen presentar un patrón muy definido y por la tanto no explican la variable respuesta, en cuanto a la emisión de Pm2.5, y las enfermedades sí parecen estar relacionadas con la esperanza de vida y presentan un patrón no aleatorio, por otra parte las variables ncidencia del VIH y la población no presentan un patrón muy claro pero pueden ser incluidas en los modelos.

4 Modelos a considerar

Se considerarán modelos con efectos fijos y efectos aleatorios considerando que \(y_{ij}\sim N(\mu_{ij},\sigma^2_y)\) donde Y:Esperanza de vida.

Del análisis descriptivo se observaron variables que podrían ser significativas, las enfermedades, las emisiones de Pm2.5, la población, la incidencia del VIH, la población, el género, el estatus. En todos los modelos se tendrá en cuenta la variable de agrupación país.

  • Modelo 1: Modelo lineal con covariables enfermedades,incidencia del VIH, género, estatus e intercepto aleatorio. \[y_{ij} \sim N(\mu_{ij},\sigma^2_y)\] \[\mu_{ij} = \beta_0 + \beta_1 Enfermedades_{ij} + \beta_2 IncVIH_{ij} + \beta_3Genero_{masculino,i} + \beta_4 Estatus_{desarrollado,ij}+b_{0i}\] \[b_{0} \sim N(0,\sigma_{b0}^2)\]

  • Modelo 2: Modelo lineal con covariables enfermedades,Incidencia del VIH, Pm2.5,género, estatus e intercepto aleatorio. \[y_{ij} \sim N(\mu_{ij},\sigma^2_y)\] \[\mu_{ij} = \beta_0 + \beta_1 Enfermedades_{ij} + \beta_2 IncVIH_{ij} + \beta_3 Pm2.5{ij} + \beta_4Genero_{masculino,i} + \beta_5 Estatus_{desarrollado,ij}+b_{0i}\] \[b_{0} \sim N(0,\sigma_{b0}^2)\]

  • Modelo 3: Modelo lineal con covariables enfermedades,Incidencia del VIH, pm2.5, poblacion, género, estatus e intercepto aleatorio. \[y_{ij} \sim N(\mu_{ij},\sigma^2_y)\] \[\mu_{ij} = \beta_0 + \beta_1 Enfermedades_{ij} + \beta_2 IncVIH_{ij} + \beta_3 Pm2.5_{ij} + \beta_4 Poblacion_{ij} \beta_5Genero_{masculino,i} + \] \[\beta_6 Estatus_{desarrollado,ij}+b_{0i}\]

  • Modelo 4: Modelo lineal con covariables enfermedades,Incidencia del VIH, pm2.5, género, estatus, pendiente e intercepto aleatorio \[y_{ij} \sim N(\mu_{ij},\sigma^2_y)\] \[\mu_{ij} = \beta_0 + \beta_1 Enfermedades_{ij} + \beta_2 IncVIH_{ij} + \beta_3 Pm2.5_{ij} + \beta_4 Genero_{masculino,i} +\] \[ \beta_5 Estatus_{desarrollado,ij}+b_{0i}+b_{1i}Enfermedades_{ij}\] \[\begin{pmatrix} b_0 \\ b_1 \end{pmatrix} \sim N \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix},\begin{pmatrix} \sigma^2_{b0} & \sigma_{b01} \\ \sigma_{b01} & \sigma^2_{b1} \end{pmatrix} \right)\]

  • Modelo 5: Modelo lineal con covariables enfermedades,Incidencia del VIH, pm2.5, poblacion, género, estatus, pendiente e intercepto aleatorio.b_{0i}+b_{1i} \[y_{ij} \sim N(\mu_{ij},\sigma^2_y)\] \[\mu_{ij} = \beta_0 + \beta_1 Enfermedades_{ij} + \beta_2 IncVIH_{ij} + \beta_3 Pm2.5_{ij} + \beta_4 Poblacion{ij} +\] \[ \beta_5Genero_{masculino,i} + \beta_6 Estatus_{desarrollado,ij}+ b_{0i}+b_{1i}Enfermedades_{ij}\] \[\begin{pmatrix} b_0 \\ b_1 \end{pmatrix} \sim N \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix},\begin{pmatrix} \sigma^2_{b0} & \sigma_{b01} \\ \sigma_{b01} & \sigma^2_{b1} \end{pmatrix} \right)\]

5 Construcción de los modelos

Para los ajustes de los modelos se utilizará la función lme del paquete nlme.

*Modelo 1

mod1 <- lme(EsperanzaDeVida ~ Enfermedades + IncidenciaVIHP + genero +
               Estatus, random = ~1 | Pais, data=datos, method = "ML")
summary(mod1)
## Linear mixed-effects model fit by maximum likelihood
##  Data: datos 
##        AIC      BIC    logLik
##   2240.422 2268.966 -1113.211
## 
## Random effects:
##  Formula: ~1 | Pais
##         (Intercept) Residual
## StdDev:    7.480325 1.999512
## 
## Fixed effects: EsperanzaDeVida ~ Enfermedades + IncidenciaVIHP + genero + Estatus 
##                        Value Std.Error  DF   t-value p-value
## (Intercept)         89.31761 1.6928418 340  52.76194   0.000
## Enfermedades        -0.94589 0.0661064 340 -14.30862   0.000
## IncidenciaVIHP      -4.72514 0.8916999 340  -5.29903   0.000
## generoMasculino     -4.97683 0.1926262 340 -25.83673   0.000
## EstatusDesarrollado  6.47878 2.1239679  91   3.05032   0.003
##  Correlation: 
##                     (Intr) Enfrmd InVIHP gnrMsc
## Enfermedades        -0.847                     
## IncidenciaVIHP       0.054 -0.218              
## generoMasculino     -0.057  0.000  0.000       
## EstatusDesarrollado -0.367  0.171  0.062  0.000
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.90875315 -0.46904905  0.02489179  0.47795932  2.56480588 
## 
## Number of Observations: 436
## Number of Groups: 93

Así, se obtiene la siguiente ecuación ajustada y los parámetros estimados:

\[y_{ij} \sim N(\mu_{ij},\hat\sigma^2_y = 1.9995^2)\] \[\hat{\mu_{ij}} = 89.31761 -0.94589 Enfermedades_{ij} -4.72514 IncVIH_{ij} - 4.97683 Genero_{masculino,i} \] \[+ 6.47878Estatus_{desarrollado,ij}+b_{0i}\]

\[b_{0} \sim N(0,\widehat{\sigma_{b0}}^2 = 7.4803)\]

*Modelo 2

mod2 <- lme(EsperanzaDeVida ~ Enfermedades + IncidenciaVIHP+Pm2.5 + 
              genero +Estatus, random = ~1 | Pais, data=datos, method = "ML")
summary(mod2)
## Linear mixed-effects model fit by maximum likelihood
##  Data: datos 
##        AIC      BIC    logLik
##   2212.528 2245.149 -1098.264
## 
## Random effects:
##  Formula: ~1 | Pais
##         (Intercept) Residual
## StdDev:    6.785562 1.964724
## 
## Fixed effects: EsperanzaDeVida ~ Enfermedades + IncidenciaVIHP + Pm2.5 + genero +      Estatus 
##                        Value Std.Error  DF   t-value p-value
## (Intercept)         92.08846 1.7044241 339  54.02908   0.000
## Enfermedades        -0.87766 0.0644378 339 -13.62027   0.000
## IncidenciaVIHP      -4.21815 0.8675437 339  -4.86218   0.000
## Pm2.5               -0.13151 0.0234802 339  -5.60076   0.000
## generoMasculino     -4.97683 0.1894948 339 -26.26368   0.000
## EstatusDesarrollado  4.58488 1.9672217  91   2.33064   0.022
##  Correlation: 
##                     (Intr) Enfrmd InVIHP Pm2.5  gnrMsc
## Enfermedades        -0.760                            
## IncidenciaVIHP       0.098 -0.188                     
## Pm2.5               -0.326 -0.139 -0.154              
## generoMasculino     -0.056  0.000  0.000  0.000       
## EstatusDesarrollado -0.401  0.153  0.037  0.177  0.000
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.70238055 -0.49998841 -0.00223797  0.48423840  2.44454160 
## 
## Number of Observations: 436
## Number of Groups: 93

Así, se obtiene la siguiente ecuación ajustada y los parámetros estimados:

\[y_{ij} \sim N(\mu_{ij},\hat\sigma^2_y = 1.964724^2)\] \[\hat{\mu_{ij}} = 92.08846 -0.87766Enfermedades_{ij} -4.21815 IncVIH_{ij}-0.13151 Pm2.5_{ij}\] \[-4.97683 Genero_{masculino,i} + 4.58488 Estatus_{desarrollado,ij}+b_{0i}\] \[b_{0} \sim N(0,\widehat{\sigma_{b0}}^2 = 6.785562^2)\] * Modelo 3

mod3 <- lme(EsperanzaDeVida ~ Enfermedades + IncidenciaVIHP+Pm2.5 + pob 
            +genero + Estatus, random = ~1 | Pais, data=datos, method = "ML")
summary(mod3)
## Linear mixed-effects model fit by maximum likelihood
##  Data: datos 
##        AIC      BIC    logLik
##   2209.703 2246.402 -1095.852
## 
## Random effects:
##  Formula: ~1 | Pais
##         (Intercept) Residual
## StdDev:    6.749198   1.9538
## 
## Fixed effects: EsperanzaDeVida ~ Enfermedades + IncidenciaVIHP + Pm2.5 + pob +      genero + Estatus 
##                        Value Std.Error  DF   t-value p-value
## (Intercept)         91.60181 1.7116548 338  53.51652  0.0000
## Enfermedades        -0.87189 0.0642127 338 -13.57816  0.0000
## IncidenciaVIHP      -4.08562 0.8658864 338  -4.71842  0.0000
## Pm2.5               -0.13522 0.0234408 338  -5.76874  0.0000
## pob                  0.01117 0.0051132 338   2.18487  0.0296
## generoMasculino     -4.97683 0.1886608 338 -26.37979  0.0000
## EstatusDesarrollado  4.79170 1.9612172  91   2.44323  0.0165
##  Correlation: 
##                     (Intr) Enfrmd InVIHP Pm2.5  pob    gnrMsc
## Enfermedades        -0.758                                   
## IncidenciaVIHP       0.088 -0.184                            
## Pm2.5               -0.313 -0.142 -0.158                     
## pob                 -0.130  0.041  0.070 -0.073              
## generoMasculino     -0.055  0.000  0.000  0.000  0.000       
## EstatusDesarrollado -0.404  0.154  0.040  0.173  0.048  0.000
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.706441666 -0.500426888 -0.002412211  0.495289244  2.457794158 
## 
## Number of Observations: 436
## Number of Groups: 93

Así, se obtiene la siguiente ecuación ajustada y los parámetros estimados:

\[y_{ij} \sim N(\mu_{ij},\sigma^2_y = 1.9538^2)\] \[hat{\mu_{ij}} = 91.60181 -0.87189Enfermedades_{ij} -4.08562IncVIH_{ij} -0.13522Pm2.5_{ij}\] \[0.01117Poblacion_{ij}-4.97683 Genero_{masculino,i} + 4.79170Estatus_{desarrollado,ij}+b_{0i}\] \[b_{0} \sim N(0,\widehat{\sigma_{b0}}^2 = 6.749198^2)\] *Modelo 4

mod4 <- lme(EsperanzaDeVida ~ Enfermedades + IncidenciaVIHP+Pm2.5 
            +genero + Estatus, random = ~1 +Enfermedades| Pais, data=datos, method = "ML")
summary(mod4)
## Linear mixed-effects model fit by maximum likelihood
##  Data: datos 
##        AIC      BIC    logLik
##   2164.996 2205.772 -1072.498
## 
## Random effects:
##  Formula: ~1 + Enfermedades | Pais
##  Structure: General positive-definite, Log-Cholesky parametrization
##              StdDev       Corr  
## (Intercept)  4.087201e-06 (Intr)
## Enfermedades 3.227998e-01 -0.012
## Residual     1.833226e+00       
## 
## Fixed effects: EsperanzaDeVida ~ Enfermedades + IncidenciaVIHP + Pm2.5 + genero +      Estatus 
##                        Value Std.Error  DF   t-value p-value
## (Intercept)         92.72845 1.4427164 339  64.27351  0.0000
## Enfermedades        -0.96570 0.0747915 339 -12.91189  0.0000
## IncidenciaVIHP      -3.00855 0.9136665 339  -3.29283  0.0011
## Pm2.5               -0.09806 0.0219629 339  -4.46481  0.0000
## generoMasculino     -4.97683 0.1768120 339 -28.14758  0.0000
## EstatusDesarrollado  2.80687 1.3419628  91   2.09162  0.0393
##  Correlation: 
##                     (Intr) Enfrmd InVIHP Pm2.5  gnrMsc
## Enfermedades        -0.760                            
## IncidenciaVIHP       0.171 -0.219                     
## Pm2.5               -0.257 -0.201 -0.133              
## generoMasculino     -0.061  0.000  0.000  0.000       
## EstatusDesarrollado -0.580  0.290 -0.036  0.173  0.000
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.0997021 -0.5136637  0.0149966  0.5155984  2.6789700 
## 
## Number of Observations: 436
## Number of Groups: 93

Así, la ecuación ajustada y los parámetros estimados quedan de la siguiente forma: \[y_{ij} \sim N(\mu_{ij},\hat\sigma^2_y = 1.833226^2 )\] \[\hat{\mu_{ij}} = 92.72845-0.96570Enfermedades_{ij}-3.00855IncVIH_{ij}-0.09806 Pm2.5_{ij}\] \[-4.97683Genero_{masculino,i}+ 2.80687Estatus_{desarrollado,ij}+b_{0i}\] Para obtener \(\hat{\sigma_{b01}}\) se despeja de la correlación de la correlación \(\rho\), por lo tanto $ = {b1} {b0} = (-0.012)(4.08720110^{-6})(0.3227998) = -1.58321710^{-7} $ \[\begin{pmatrix} b_0 \\ b_1 \end{pmatrix} \sim N \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix},\begin{pmatrix} \widehat{\sigma^2_{b0}} = 4.087201\times10^{-12} & \widehat{\sigma_{b01}} = -1.583217\times 10^{-07} \\ \widehat{\sigma_{b01}} = -1.583217\times 10^{-07} & \widehat{\sigma^2_{b1}} =0.3227998^2 \end{pmatrix} \right)\]

*Modelo 5

mod5 <- lme(EsperanzaDeVida ~ Enfermedades + IncidenciaVIHP+Pm2.5 + pob 
            +genero + Estatus, random = ~1 +Enfermedades| Pais, data=datos, method = "ML")
summary(mod5)
## Linear mixed-effects model fit by maximum likelihood
##  Data: datos 
##        AIC      BIC    logLik
##   2162.629 2207.483 -1070.315
## 
## Random effects:
##  Formula: ~1 + Enfermedades | Pais
##  Structure: General positive-definite, Log-Cholesky parametrization
##              StdDev       Corr  
## (Intercept)  2.207429e-06 (Intr)
## Enfermedades 3.216917e-01 -0.014
## Residual     1.823314e+00       
## 
## Fixed effects: EsperanzaDeVida ~ Enfermedades + IncidenciaVIHP + Pm2.5 + pob +      genero + Estatus 
##                        Value Std.Error  DF   t-value p-value
## (Intercept)         92.35144 1.4490939 338  63.73047  0.0000
## Enfermedades        -0.96012 0.0745860 338 -12.87265  0.0000
## IncidenciaVIHP      -2.91150 0.9112808 338  -3.19495  0.0015
## Pm2.5               -0.10139 0.0219465 338  -4.61999  0.0000
## pob                  0.00954 0.0045897 338   2.07805  0.0385
## generoMasculino     -4.97683 0.1760609 338 -28.26766  0.0000
## EstatusDesarrollado  2.89700 1.3388456  91   2.16380  0.0331
##  Correlation: 
##                     (Intr) Enfrmd InVIHP Pm2.5  pob    gnrMsc
## Enfermedades        -0.757                                   
## IncidenciaVIHP       0.164 -0.217                            
## Pm2.5               -0.245 -0.203 -0.137                     
## pob                 -0.127  0.038  0.047 -0.074              
## generoMasculino     -0.061  0.000  0.000  0.000  0.000       
## EstatusDesarrollado -0.580  0.290 -0.035  0.169  0.034  0.000
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.103397258 -0.523661011  0.005928518  0.523347778  2.691527979 
## 
## Number of Observations: 436
## Number of Groups: 93

Así, la ecuación ajustada y los parámetros estimados quedan de la siguiente forma: \[y_{ij} \sim N(\mu_{ij},\sigma^2_y = 1.823314^2 )\] \[\mu_{ij} = 92.35144-0.96012 Enfermedades_{ij}-2.91150IncVIH_{ij} -0.10139 Pm2.5_{ij}\] \[+0.00954Población_{ij}-4.97683 Genero_{masculino,i}+ 2.89700Estatus_{desarrollado,ij}+b_{0i}\] Para obtener \(\sigma_{b01}\) se despeja de la correlación de la correlación \(\rho\), por lo tanto \(\hat\sigma_{b01} = \rho \times \sigma_{b1} \times \sigma_{b0} = (-0.014)(2.20749\times10^{-6})(0.3216917) = -9.941562e-08\times 10^{-8}\)

\[\begin{pmatrix} b_0 \\ b_1 \end{pmatrix} \sim N \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix},\begin{pmatrix} \hat\sigma^2_{b0} = 2.207429\times10^{-12} & \hat\sigma_{b01} = 0\\ \hat\sigma_{b01} = 0 & \hat\sigma^2_{b1} =0.3.216917^2 \end{pmatrix} \right)\]

6 Comparación de los modelos

Para la comparación de los modelos ajustados se utilizarán pruebas de hipótesis, la correlación entre los valores ajustados y los valores observados, para verificar si los efectos fijos o aleatorios aportan al modelo, inicialmente se considera el modelo 1 como el modelo base, se procede a verificar si la variable Pm2.5 incluida en el modelo 2 es significativa:

6.1 Pruebas de hipótesis

  • \(H_0\): la variable Pm2.5 no aporta al modelo vs \(H_1\): la variable Pm2.5 aporta al modelo

Para obtener el resultado de esta prueba y de las pruebas posteriores se utilizará la función anova

anova(mod1,mod2)

Como el valor p es pequeño y menor que 0.0001, se rechaza \(H_0\), la variable Pm2.5 es significativa en el modelo.

En el modelo 3 se incluye la variable Población, por tanto se realiza el siguiente test (en este caso se compara con el modelo 2)

  • \(H_0\): la variable Población no aporta al modelo vs \(H_1\): la variable Población aporta al modelo
anova(mod2,mod3)

En este caso el valor p también es algo pequeño, podría decirse que aporta al modelo pero en una medida muy pequeña además de las estimaciones obtenidas para dicho modelo el \(\beta\) para dicho parámetro fue muy cercano a cero, por lo tanto se concluye que no es significativo.

Ahora, se incluirá otro efecto aleatorio que es la pendiente, en este caso se deben hacer pruebas de hipótesis sobre la varianza, la pendiente se incluyó para la variable Enfermedades y se tuvieron como base los modelos 2 y 3

Con modelo 2 como base * \(H_0\): la pendiente aleatoria no aporta al modelo vs \(H_1\): la pendiente aleatoria no aporta al modelo.

anova(mod2,mod4)

Con modelo 3 como base

  • \(H_0\): la pendiente aleatoria no aporta al modelo vs \(H_1\): la pendiente aleatoria no aporta al modelo.
anova(mod3,mod5)

En ambos modelos la pendiente aleatoria resulta significativa.

6.2 Correlación entre valores observados y valores ajustados

correlaciones <- data.frame(cor_mod1 = cor(fitted(mod1),datos$EsperanzaDeVida),
                            cor_mod2 = cor(fitted(mod2),datos$EsperanzaDeVida),
                            cor_mod3 = cor(fitted(mod3),datos$EsperanzaDeVida),
                            cor_mod4 = cor(fitted(mod4),datos$EsperanzaDeVida),
                            cor_mod5 = cor(fitted(mod5),datos$EsperanzaDeVida))
correlaciones

Para esta medida todos los modelos obtuvieron un buen desempeño.

6.3 Criterios AIC y BIC

anova(mod1,mod2,mod3,mod4,mod5)

En cuanto a estas medidas los modelos que obtuvieron mejor desempeño fueron los modelos 4 y 5.

6.4 Conclusión

Todos los modelos obtuvieron una correlación muy alta en cuanto a sus valores ajustados y los valores observados, con las pruebas de hipótesis se concluye que el mejor modelo es el modelo 4, en caso de presentar violación en los supuestos se tendrá en observación el modelo 5. Los efectos aleatorios mejoran significativamente los modelos, pues la inclusión de dichos efectos mejoraron el modelo.

7 Análisis de residuales

residuales <- residuals(mod4)
ajustados <- fitted(mod4)
resid <- data.frame(residuales = residuales, fitt = ajustados)

7.1 Gráficos de residuales

res_vs_fitt <- resid %>% ggplot(aes(x=fitt, y=residuales)) + geom_point(alpha = 0.6, size = 1)+ 
  theme_bw() +
  geom_abline(intercept = 0,slope = 0, color = "red", size = 1.5, alpha = 0.8) +
  xlab("Ajustados")

res_vs_nrow <- resid %>% ggplot(aes(x=1:nrow(resid), y=residuales)) + geom_point(alpha = 0.6, size = 1)+ 
  theme_bw() +
  geom_abline(intercept = 0,slope = 0, color = "red", size = 1.5, alpha = 0.8) +
  xlab("nrow")

dens_res <- ggplot(resid,aes(x=residuales)) +
  geom_density(alpha=0.55) +
  theme_bw()

qqplot_Res <- ggplot(resid, aes(sample = residuales)) +
  stat_qq() +
  stat_qq_line(color = "red") +
  theme_bw()

ggarrange(res_vs_fitt,res_vs_nrow,dens_res,qqplot_Res)

Al ver el gráfico de residuales vs ajustados puede observarse un patrón aleatorio, aunque para algunos puntos parece haber una concentración de datos para valores altos de los valores ajustados no parece presentar un problema grande en cuanto a varianza constante, por otra parte los residuales se encuentran alrededor del cero por tanto cumple el supuesto de media cero.

Para el gráfico Normal QQplot para los residuales se puede evidenciar que para las colas de la distribución los valores se alejan de la linea y por tanto hay problemas en el supuesto de normalidad.

7.2 Pruebas para los efectos aleatorios

b0 <- ranef(mod4)[,1]
b1 <- ranef(mod4)[,2]
mix_efec <- data.frame(b0=b0,b1=b1)

qqplot_b0 <- ggplot(mix_efec, aes(sample = b0)) +
  stat_qq(alpha = 0.6) +
  stat_qq_line(color = "red") +
  theme_bw()

qqplot_b1 <- ggplot(mix_efec, aes(sample = b1)) +
  stat_qq(alpha=0.6) +
  stat_qq_line(color = "red") +
  theme_bw()

qqplot_b0

qqplot_b1

colMeans(mix_efec)
##            b0            b1 
## -3.630888e-21 -9.652221e-16

Para la normalidad de los efectos \(b_0\) y \(b_1\) se puede observar que hay problemas en el supuesto de que estos distribuyen normal, pues hay varios puntos que están alejados de la linea(para ambos efectos), este supuesto debe revisarse aunque es bastante flexible, al calcular las medias se puede ver claramente que ambas son cero.

8 Mejor modelo

El mejor modelo es el modelo 4, este modelo considera las covariables enfermedades, incidencia del VIH, Pm2.5, poblacion, género, estatus, pendiente e intercepto aleatorio, se puede concluir que la inclusión de efectos aleatorios en el modelo permiten mejorar de manera significativa datos que tengan variables de agrupación, en este caso país tal y como se vio en el análisis descriptivo el comportamiento de cada país es diferente.

\[y_{ij} \sim N(\mu_{ij},\hat\sigma^2_y = 1.833226^2 )\] \[\hat\mu_{ij} = 92.72845-0.96570Enfermedades_{ij}-3.00855IncVIH_{ij}-0.09806 Pm2.5_{ij}\] \[-4.97683Genero_{masculino,i}+ 2.80687Estatus_{desarrollado,ij}+b_{0i}\] \[\begin{pmatrix} b_0 \\ b_1 \end{pmatrix} \sim N \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix},\begin{pmatrix} \hat\sigma^2_{b0} = 4.087201\times10^{-12} &\hat \sigma_{b01} = -1.583217\times 10^{-07} \\ \hat\sigma_{b01} = -1.583217\times 10^{-07} & \hat\sigma^2_{b1} =0.3227998^2 \end{pmatrix} \right)\]

Supongamos entonces que queremos estimar la esperanza de vida de un hombre de Argentina donde el porcentaje de 30 años que morirá antes de los 70 por __enfermedades_ cardiovasculares, cáncer, diabetes es de 22%%, la incidencia del VIH es de 0.05, el material particulado pm2.5 es de 12 microgramos por metro, cabe aclarar que Argentina es un país en desarrollo

new_data <- data.frame(Pais = "Argentina",
                       genero = "Masculino",
                       Enfermedades = 20,
                       IncidenciaVIHP = 0.04,
                       Pm2.5 = 15,
                       Estatus = "En desarrollo")
predict(mod4,new_data)
## Argentina 
##  71.49381 
## attr(,"label")
## [1] "Predicted values"

Luego, la esperanza de vida para una persona con dichas características es \(71.49381\)

Veamos ahora las mismas características pero para una persona de género femenino

new_data2 <- data.frame(Pais = "Argentina",
                       genero = "Femenino",
                       Enfermedades = 20,
                       IncidenciaVIHP = 0.04,
                       Pm2.5 = 15,
                       Estatus = "En desarrollo")
predict(mod4,new_data2)
## Argentina 
##  76.47064 
## attr(,"label")
## [1] "Predicted values"

El valor estimado fue \(76.47\) años,es decir, la diferencia para hombres y mujeres en Argentina es de aproximadamente 5 años, en este caso las mujeres tienen una esperanza de vida más alta # Conclusiones * El uso de modelos mixtos permitió realizar un buen ajuste, las covariables excogidas lograron explicar en gran medida la variable respuesta. * La inclusión de los efectos aleatorios mejoraron los modelos puesto que en los datos había una variable de agrupación. # Referencias