setwd("~/ea1130")
library(pacman)
p_load("prettydoc", "DT", "xfun", "dplyr", "psych", "GGally", "ggplot2","readr", "gridExtra")Caso de estudio: Esperanza de vida en Estados Unidos.
Esquema del lenguaje natural
En el siguiente caso de estudio se estudia un aspecto que es muy importante para la sociedad hoy en día, que es la esperanza de vida en Estados unidos, basándose en diferentes variables. Para poder explicar este fenómeno se usaron distintas variables que pudieran ser relevantes, entre las cuales destacan:
- habitantes.
- analfabetismo.
- ingresos.
- esperanza de vida.
- asesinatos.
- universitarios.
- heladas.
- área.
- densidad poblacional.
Para realizar esta investigación usamos datos del dataset de Rstudio State.x77 que muestra los datos de los 50 estados de Estados Unidos.
Objetivo
En este caso de estudio se estará realizando paso a paso un modelo que nos permita predecir la esperanza de vida media de los habitantes.
Nuestra hipotesis es que la esperanza de vida se puede predecir tomando en cuenta las helas, los asesinatos y los ingresos. Se tomarán en cuenta debido a que las heladas pueden afectar tanto a los cultivos como a la salud de las personas, cierta población puede depender de los cultivos y si éstos mueren por la helada no tendrán de donde sacar ingresos o de que comer, por otro lado las heladas pueden afectar a la salud haciendo a las personas más propensos a enfermarse, los asesinatos afectan directamente a la esperanza de vida, pues si son muvhos los que ocurren hay riesgo de que las personas no avancen de cierta edad y tener un buen ingreso es fundamental para acceder a las necesidades basicas.
¿Qué es la esperanza de vida?
La esperanza de vida se define como el número de años de vida que restan, término medio, a una persona de no variar la tendencia de la mortalidad (Haupt y Kane, 2003: 58). En otras palabras, es una síntesis de las tasas específicas de mortalidad. El indicador resume la intensidad de la mortalidad que experimenta una población, ya sea al nacimiento o en una edad alcanzada específica. Con frecuencia la esperanza de vida se utiliza como indicador de desarrollo humano o bienestar, como en el Programa de las Naciones Unidas para el Desarrollo (PNUD, 2012); también se suele emplear para fines comparativos entre poblaciones residentes en distintas unidades geográficas (países, entidades, zonas metropolitanas, etc.). Resulta importante su monitoreo de manera periódica con el fin de estudiar su evolución respecto a sí misma, garantizar la comparabilidad con otras poblaciones, y en ocasiones, propiciar la toma de decisiones informada en materia de salud pública y desarrollo social.
Fuente del concepto de esperanza de vida: http://www.scielo.org.mx/pdf/educm/v32n1/2448-6515-educm-32-01-00097.pdf
Esquema del lenguaje natural
Datos
datos <- as.data.frame(state.x77)
datatable(datos)Renombran y modificar el conjunto de datos
Para que sean más entendibles las variables se renombran, cambiandolas del inglés al español. *****
datos <- rename(habitantes = Population,
analfabetismo = Illiteracy,
ingresos = Income,
esp_vida = `Life Exp`,
asesinatos = Murder,
universitarios = `HS Grad`,
heladas = Frost,
area = Area,.data = datos)
datos <- mutate(.data=datos, densidad_pobl = habitantes * 1000 / area)Tabla con las varianbles cambiadas a idioma español
datatable(datos)1.Analizar la relación entre variables
Lo primero que tenemos que hacer para poder realizar el modelo lineal multiple es estudiar la relacion que existen entre las variables, esto nos ayudara a saber cuales son las mejores predicciones para el modelo
round(cor(x = datos, method = "pearson"), 3)## habitantes ingresos analfabetismo esp_vida asesinatos
## habitantes 1.000 0.208 0.108 -0.068 0.344
## ingresos 0.208 1.000 -0.437 0.340 -0.230
## analfabetismo 0.108 -0.437 1.000 -0.588 0.703
## esp_vida -0.068 0.340 -0.588 1.000 -0.781
## asesinatos 0.344 -0.230 0.703 -0.781 1.000
## universitarios -0.098 0.620 -0.657 0.582 -0.488
## heladas -0.332 0.226 -0.672 0.262 -0.539
## area 0.023 0.363 0.077 -0.107 0.228
## densidad_pobl 0.246 0.330 0.009 0.091 -0.185
## universitarios heladas area densidad_pobl
## habitantes -0.098 -0.332 0.023 0.246
## ingresos 0.620 0.226 0.363 0.330
## analfabetismo -0.657 -0.672 0.077 0.009
## esp_vida 0.582 0.262 -0.107 0.091
## asesinatos -0.488 -0.539 0.228 -0.185
## universitarios 1.000 0.367 0.334 -0.088
## heladas 0.367 1.000 0.059 0.002
## area 0.334 0.059 1.000 -0.341
## densidad_pobl -0.088 0.002 -0.341 1.000Análisis con histogramas
multi.hist( x = datos, dcol = c("blue","red"), dlty = c("dotted", "solid"),
main = ""
)
Podemos aprecias los historigramas de las 9 variables, en la mayoria de los casos podremos ver un cambio casi nulo e imperceptible, pero los datos ahí estan, Ahora veremos lo anterios pero ahora representando utilizando ggplot y ggally incluyendo una recta. *****
ggpairs(datos, lower = list(continuous ="smooth"),
diag = list (continuos = "barDiag"), axisLabels = "none"
)
Entonces, con los resultados de los análisis realizados hasta el momento, podemos obtener las siguientes conclusiones preliminares:
Por logica podiamos decir que las variables relacionadas a la esperanza de vida son los asesinatos, el alfabetismo y universitarios.
El punltimo y antepenultimo est muy relacionadas.
Generar el primer modelo de regresión lineal simple
Emplearemos el método mixto, para esto utilizaremos todas las variables (para este caso, los predictores) y despues seleccionaremos los mejores predictores empleando la medición Akaitke(AIC).
modelo <- lm(esp_vida ~ asesinatos, data=datos)
summary(modelo)##
## Call:
## lm(formula = esp_vida ~ asesinatos, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.81690 -0.48139 0.09591 0.39769 2.38691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 72.97356 0.26997 270.30 < 2e-16 ***
## asesinatos -0.28395 0.03279 -8.66 2.26e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8473 on 48 degrees of freedom
## Multiple R-squared: 0.6097, Adjusted R-squared: 0.6016
## F-statistic: 74.99 on 1 and 48 DF, p-value: 2.26e-11Tomando en cuenta lo anterior, usaremos esta formula para representar el modelo.
\[ y = 72.97356 - 0.28395x \]
Procedemos a encontrar Y a traves del residuo y producto.
y = 72.97356 - (0.28395*11.6)
y## [1] 69.67974Obetenemos los resultados.
y - 67.96## [1] 1.71974Procedemos a graficar los resultados anteriores.
plot(datos$asesinatos,col="#104E8B", datos$esp_vida, xlab="asesinatos", ylab="esperanza de vida")
abline(modelo)
Ahora entraremos a detalles, esta ves con múltiples predictores
modelo2 <- lm(esp_vida ~ habitantes + ingresos + analfabetismo + asesinatos + universitarios + heladas + area + densidad_pobl, data = datos )
summary(modelo2) ##
## Call:
## lm(formula = esp_vida ~ habitantes + ingresos + analfabetismo +
## asesinatos + universitarios + heladas + area + densidad_pobl,
## data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47514 -0.45887 -0.06352 0.59362 1.21823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.995e+01 1.843e+00 37.956 < 2e-16 ***
## habitantes 6.480e-05 3.001e-05 2.159 0.0367 *
## ingresos 2.701e-04 3.087e-04 0.875 0.3867
## analfabetismo 3.029e-01 4.024e-01 0.753 0.4559
## asesinatos -3.286e-01 4.941e-02 -6.652 5.12e-08 ***
## universitarios 4.291e-02 2.332e-02 1.840 0.0730 .
## heladas -4.580e-03 3.189e-03 -1.436 0.1585
## area -1.558e-06 1.914e-06 -0.814 0.4205
## densidad_pobl -1.105e-03 7.312e-04 -1.511 0.1385
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7337 on 41 degrees of freedom
## Multiple R-squared: 0.7501, Adjusted R-squared: 0.7013
## F-statistic: 15.38 on 8 and 41 DF, p-value: 3.787e-10El modelo una vez con todas las variables introducidas tienen un R2 de 0.7501 lo que significa que el 75.01% de las variables se relacionan en la esperanza de vida.el p-value del modelo tiene 3.787-10 lo que significa que el modelo no es por azar.
Seleccionando los mejores predictores
Se van a emplear la estrategia de stepwise mixto y el valor matemático empleado para determinar la calidad del modelo va a ser el criterio de información de Akaike (AIC), éste es una medida de la calidad relativa de un modelo estadístico, nos proporciona un medio para la selección del modelo.
AIC maneja un trade-off entre la bondad de ajuste del modelo y la complejidad del modelo. Nos Ofrece una estimación relativa de la información perdida cuando se utiliza un modelo determinado para representar el proceso que genera los datos.
IC no rd capaz de proporcionar una prueba de un modelo en el sentido de probar una hipótesis nula, esto quiere decir que AIC no puede decir nada acerca de la calidad del modelo en un sentido absoluto. AIC no dará ningún aviso si todos los modelos candidatos encajan mal.
step(object = modelo2, direction="both", trace = 1)## Start: AIC=-22.89
## esp_vida ~ habitantes + ingresos + analfabetismo + asesinatos +
## universitarios + heladas + area + densidad_pobl
##
## Df Sum of Sq RSS AIC
## - analfabetismo 1 0.3050 22.373 -24.208
## - area 1 0.3564 22.425 -24.093
## - ingresos 1 0.4120 22.480 -23.969
## <none> 22.068 -22.894
## - heladas 1 1.1102 23.178 -22.440
## - densidad_pobl 1 1.2288 23.297 -22.185
## - universitarios 1 1.8225 23.891 -20.926
## - habitantes 1 2.5095 24.578 -19.509
## - asesinatos 1 23.8173 45.886 11.707
##
## Step: AIC=-24.21
## esp_vida ~ habitantes + ingresos + asesinatos + universitarios +
## heladas + area + densidad_pobl
##
## Df Sum of Sq RSS AIC
## - area 1 0.1427 22.516 -25.890
## - ingresos 1 0.2316 22.605 -25.693
## <none> 22.373 -24.208
## - densidad_pobl 1 0.9286 23.302 -24.174
## - universitarios 1 1.5218 23.895 -22.918
## + analfabetismo 1 0.3050 22.068 -22.894
## - habitantes 1 2.2047 24.578 -21.509
## - heladas 1 3.1324 25.506 -19.656
## - asesinatos 1 26.7071 49.080 13.072
##
## Step: AIC=-25.89
## esp_vida ~ habitantes + ingresos + asesinatos + universitarios +
## heladas + densidad_pobl
##
## Df Sum of Sq RSS AIC
## - ingresos 1 0.132 22.648 -27.598
## - densidad_pobl 1 0.786 23.302 -26.174
## <none> 22.516 -25.890
## - universitarios 1 1.424 23.940 -24.824
## + area 1 0.143 22.373 -24.208
## + analfabetismo 1 0.091 22.425 -24.093
## - habitantes 1 2.332 24.848 -22.962
## - heladas 1 3.304 25.820 -21.043
## - asesinatos 1 32.779 55.295 17.033
##
## Step: AIC=-27.6
## esp_vida ~ habitantes + asesinatos + universitarios + heladas +
## densidad_pobl
##
## Df Sum of Sq RSS AIC
## - densidad_pobl 1 0.660 23.308 -28.161
## <none> 22.648 -27.598
## + ingresos 1 0.132 22.516 -25.890
## + analfabetismo 1 0.061 22.587 -25.732
## + area 1 0.043 22.605 -25.693
## - habitantes 1 2.659 25.307 -24.046
## - heladas 1 3.179 25.827 -23.030
## - universitarios 1 3.966 26.614 -21.529
## - asesinatos 1 33.626 56.274 15.910
##
## Step: AIC=-28.16
## esp_vida ~ habitantes + asesinatos + universitarios + heladas
##
## Df Sum of Sq RSS AIC
## <none> 23.308 -28.161
## + densidad_pobl 1 0.660 22.648 -27.598
## + ingresos 1 0.006 23.302 -26.174
## + analfabetismo 1 0.004 23.304 -26.170
## + area 1 0.001 23.307 -26.163
## - habitantes 1 2.064 25.372 -25.920
## - heladas 1 3.122 26.430 -23.877
## - universitarios 1 5.112 28.420 -20.246
## - asesinatos 1 34.816 58.124 15.528##
## Call:
## lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
## heladas, data = datos)
##
## Coefficients:
## (Intercept) habitantes asesinatos universitarios heladas
## 7.103e+01 5.014e-05 -3.001e-01 4.658e-02 -5.943e-03Según la llamada final de este proceso (AIC) tenemos que el mejor modelo es: *****
modelo3 <- (lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
heladas, data = datos))
summary(modelo3)##
## Call:
## lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
## heladas, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47095 -0.53464 -0.03701 0.57621 1.50683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.103e+01 9.529e-01 74.542 < 2e-16 ***
## habitantes 5.014e-05 2.512e-05 1.996 0.05201 .
## asesinatos -3.001e-01 3.661e-02 -8.199 1.77e-10 ***
## universitarios 4.658e-02 1.483e-02 3.142 0.00297 **
## heladas -5.943e-03 2.421e-03 -2.455 0.01802 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7197 on 45 degrees of freedom
## Multiple R-squared: 0.736, Adjusted R-squared: 0.7126
## F-statistic: 31.37 on 4 and 45 DF, p-value: 1.696e-12intervalo de confianza de los coeficientes
En la siguiente tabla se puede apreciar el intervalo de confianza para cada uno de los coeficientes parciales de regresión:
confint(lm(formula = esp_vida ~ habitantes + asesinatos + universitarios + heladas, data = datos))## 2.5 % 97.5 %
## (Intercept) 6.910798e+01 72.9462729104
## habitantes -4.543308e-07 0.0001007343
## asesinatos -3.738840e-01 -0.2264135705
## universitarios 1.671901e-02 0.0764454870
## heladas -1.081918e-02 -0.0010673977Se puede observar que por cada unidad que aumenta el predictor universitarios, la esperanza de vida aumenta en promedio 0.04658 unidades, manteniéndose constantes el resto de predictores
Validación de condiciones
plot1 <- ggplot(data = datos, aes(habitantes, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = datos, aes(asesinatos, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = datos, aes(universitarios, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot4 <- ggplot(data = datos, aes(heladas, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2, plot3, plot4)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Grafica Habitantes: Se puede apreciar que hay una gran concentración de los datos entre 0 y 6000, es por esto que para tener un margen de error minimo no se tomarán los datos posteriores a 6000.
Grafica Asesinatos: En esta gráfica podemos darnos cuentas que los datos están mayormente concentrados entre -1 y 1, esto nos dice que la tendencia de los asesinatos eran mayormente eventos aislados.
Grafica Universitarios: Se puede observar que la cantidad de universitarios fue incrementando.
Grafica Heladas: En la gráfica podemos observar que, a pesar de la dispersión de los datos, claramente las heladas fueron disminuyendo y esto se puede deber a multiples factores.
Distribución normal de los residuos:
qqnorm(modelo3$residuals, col="#8B3A3A")
qqline(modelo3$residuals)
Podemos observar que la gráfica tiende a seguir la línea recta entre el -1 y el 1, por lo que se podría tomar como un buen indicador de que los residuos están distribuidos de fomra normal.
Prueba de shapiro wilk para determinar normalidad
shapiro.test(modelo$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.98431, p-value = 0.7416Al realizar el analisis de Shapiro-wilk se puede concluir que que los datos si están tocando la línea de tendencia y por lo tanto el análisis gráfico y el test de hipótesis confirman la normalidad.
Comparación de valores ajustados y valores residuales
ggplot(data = datos, aes(modelo$fitted.values, modelo$residuals)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Sin tomar en cuenta los valores de los extremos, se puede decir que el modelo está bastante acertado, pues los datos no están tan dispersos.
Conclusión
En conclusión, son varios los factores que pueden afectar a la esperanza de vida, pero se puede predecir tomando en cuenta la cantidad de habitantes, porcentaje de universitarios, el número de asesinatos y el promedio de heladas en una ciudad.