Regresión
24 junio 2021
Contents
1 Introducción
Los supuestos subyacentes de un modelo de regresión lineal son:
- Varianza constante (supuesto de homocedasticidad)
- Los residuales se distribuyen normalmente
- Sin multicolinealidad entre predictores (o muy poca)
- Relación lineal entre la variable de respuesta y los predictores
En este ejemplo se considera la esperanza de vida como variable de respuesta y un modelo con fines de inferencia, i.e., un modelo complicado con términos de interacción solo para obtener una mayor precisión de predicción, si no que interesa un modelos que sea muy interpretable. La paquetería encesaria es:
library(dplyr)
library(gapminder)
library(leaps)
library(glmnet)
library(ggplot2)
library(MASS)
library(usdm)
library(car)
library(corrplot)
library(texreg)
library(MPV)2 Datos
El análisis toma el conjunto de datos gapminder y se va a fusionar con otro conjunto de https://www.kaggle.com/kumarajarshi/life-expectancy-who
data("gapminder")
gap.2002 <- filter(gapminder, year == 2002)
gap.2007 <- filter(gapminder, year == 2007)
gap.2002 <- rename(gap.2002, Country = country)
gap.2007 <- rename(gap.2007, Country = country)lifeExp <- read.csv("Life Expectancy Data.csv")
lifeExp.2002 <- filter(lifeExp, Year == 2002)
lifeExp.2007 <- filter(lifeExp, Year == 2007)Se toman sólo los años de 2002 y 2007 y se fusionan en país para cada año. Después, podemos hacer una rbind durante estos dos años.
gapExp.2002 <- merge(lifeExp.2002, gap.2002, by = "Country")
gapExp.2007 <- merge(lifeExp.2007, gap.2007, by = "Country")
gapExp <- rbind(gapExp.2002, gapExp.2007)Se descartan las variables que aparecen dos veces en la base de datos, la variable Hepatitis.B por su gran cantidad de valores de NA y las variables Status y continent por falta de significado físico.
gapExp <- dplyr::select(gapExp, -c(Year, year, Life.expectancy, Population, GDP, Country, Hepatitis.B, Income.composition.of.resources))
gapExp <- na.omit(gapExp)
gapExp.muu = dplyr::select(gapExp, -c(Status, continent))
(data.frame(Variables = colnames(gapExp.muu)))## Variables
## 1 Adult.Mortality
## 2 infant.deaths
## 3 Alcohol
## 4 percentage.expenditure
## 5 Measles
## 6 BMI
## 7 under.five.deaths
## 8 Polio
## 9 Total.expenditure
## 10 Diphtheria
## 11 HIV.AIDS
## 12 thinness..1.19.years
## 13 thinness.5.9.years
## 14 Schooling
## 15 lifeExp
## 16 pop
## 17 gdpPercap3 Comprobación de Supuestos
Al final, se cuenta con \(236\) observaciones, \(16\) predictores y una variable de respuesta (lifeExp).
dim(gapExp.muu)## [1] 236 17Primero, se ajusta un modelo con todos los predictores incluidos para evaluar el supuesto de varianza constante.
all.reg <- lm(lifeExp ~. ,data = gapExp.muu)
summary(all.reg)##
## Call:
## lm(formula = lifeExp ~ ., data = gapExp.muu)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.9242 -2.1037 0.5952 2.6685 15.4555
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.600e+01 2.288e+00 20.102 < 2e-16 ***
## Adult.Mortality -9.800e-03 2.855e-03 -3.432 0.000715 ***
## infant.deaths 6.868e-02 3.998e-02 1.718 0.087223 .
## Alcohol -2.530e-01 1.143e-01 -2.214 0.027846 *
## percentage.expenditure -1.164e-04 1.933e-04 -0.602 0.547543
## Measles 1.260e-05 5.671e-05 0.222 0.824411
## BMI 6.356e-02 2.327e-02 2.732 0.006814 **
## under.five.deaths -4.963e-02 2.757e-02 -1.800 0.073158 .
## Polio 4.774e-02 2.293e-02 2.082 0.038483 *
## Total.expenditure 4.224e-02 1.473e-01 0.287 0.774576
## Diphtheria 3.442e-02 1.974e-02 1.744 0.082620 .
## HIV.AIDS -6.101e-01 4.800e-02 -12.712 < 2e-16 ***
## thinness..1.19.years -4.188e-01 1.255e-01 -3.336 0.000999 ***
## thinness.5.9.years 2.438e-01 1.279e-01 1.905 0.058057 .
## Schooling 1.280e+00 1.814e-01 7.056 2.23e-11 ***
## pop 2.791e-09 5.893e-09 0.474 0.636276
## gdpPercap 1.688e-04 5.072e-05 3.328 0.001025 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.613 on 219 degrees of freedom
## Multiple R-squared: 0.8694, Adjusted R-squared: 0.8598
## F-statistic: 91.1 on 16 and 219 DF, p-value: < 2.2e-16plot(all.reg, which = 1)
Se observa que no se cumple el supuesto de homocedasticidad pues hay un alza y luego una baja muy notorias. La forma de corregirlo es realizando transformaciones.
plot(all.reg, which = 2)
En la gráfica anterior se ve que, aproximadamente los residuales están distribuidos normalmente. La variación descendente de los residuales a la izquierda y la variación ascendente a la derecha de la gráfica sugiere que la distribución de los residuos tiene una cola más pesada que la distribución teórica.
Para el supuesto de no multicolinealidad se observan los factores de inflación de la varianza (VIF). Cuando el factor de inflación de la varianza está por encima de \(5\), existe multicolinealidad. A continuación se muestra que infant.deaths.1 y under.five.deaths.1 tienen factores de inflación de varianza muy alta. Una forma de lidiar con eso es centrar estas dos variables y ver si los valores de VIF disminuyen.
vif.1 <- usdm::vif(gapExp.muu)
vif.1## Variables VIF
## 1 Adult.Mortality 1.897464
## 2 infant.deaths 452.041800
## 3 Alcohol 2.212311
## 4 percentage.expenditure 2.194532
## 5 Measles 3.891227
## 6 BMI 2.430624
## 7 under.five.deaths 407.652335
## 8 Polio 2.546104
## 9 Total.expenditure 1.254948
## 10 Diphtheria 2.246485
## 11 HIV.AIDS 2.466771
## 12 thinness..1.19.years 4.314330
## 13 thinness.5.9.years 4.573660
## 14 Schooling 5.769637
## 15 lifeExp 7.655801
## 16 pop 9.316206
## 17 gdpPercap 4.342583centre <- function(x) { x - mean(x)}
gapExp.centered <- data.frame(gapExp.muu, infant.deaths = centre(gapExp.muu$infant.deaths), under.five.deaths = centre(gapExp.muu$under.five.deaths))
gapExp.centered <- dplyr::select(gapExp.centered, -c(under.five.deaths, infant.deaths))
usdm::vif(gapExp.centered)## Variables VIF
## 1 Adult.Mortality 1.897464
## 2 Alcohol 2.212311
## 3 percentage.expenditure 2.194532
## 4 Measles 3.891227
## 5 BMI 2.430624
## 6 Polio 2.546104
## 7 Total.expenditure 1.254948
## 8 Diphtheria 2.246485
## 9 HIV.AIDS 2.466771
## 10 thinness..1.19.years 4.314330
## 11 thinness.5.9.years 4.573660
## 12 Schooling 5.769637
## 13 lifeExp 7.655801
## 14 pop 9.316206
## 15 gdpPercap 4.342583
## 16 infant.deaths.1 452.041800
## 17 under.five.deaths.1 407.652335Desafortunadamente, centrar las variables no ayuda a reducir los valores de VIF. Esto obliga a desechar una de estas variables. Eliminando under.five.deaths.1:
gapExp.centered <- dplyr::select(gapExp.centered, -under.five.deaths.1)
usdm::vif(gapExp.centered)## Variables VIF
## 1 Adult.Mortality 1.894701
## 2 Alcohol 2.114592
## 3 percentage.expenditure 2.192347
## 4 Measles 3.189797
## 5 BMI 2.428769
## 6 Polio 2.545950
## 7 Total.expenditure 1.253784
## 8 Diphtheria 2.230882
## 9 HIV.AIDS 2.460571
## 10 thinness..1.19.years 4.286697
## 11 thinness.5.9.years 4.466637
## 12 Schooling 5.729722
## 13 lifeExp 7.544128
## 14 pop 5.524167
## 15 gdpPercap 4.330769
## 16 infant.deaths.1 3.554883Ahora, todos los valores de VIF están por debajo de \(10\), una mejora considerable. Para continuar verificando supuestos, se debe trabajar con los valores de VIF que están por encima de 5.
También se asume una relación lineal entre la variable de respuesta y los predictores. Esta suposición se comprueba con un diagrama de dispersión.
ggplot(gapExp, aes(x = HIV.AIDS , y = lifeExp, col = continent)) +
geom_point() +
geom_smooth(method = "lm", formula= y ~ log(x), se = F, aes(group = 1), col = "black") +
theme_minimal()
ggplot(gapExp, aes(x = log(HIV.AIDS + 1) , y = lifeExp, col = continent)) +
geom_point() +
geom_smooth(method = "lm", formula= y ~ x, se = F, aes(group = 1), col = "black") +
theme_minimal()
ggplot(gapExp, aes(x = gdpPercap , y = lifeExp, col = continent)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ log(x), aes(group = 1), col = "black", se = F) +
theme_minimal()
ggplot(gapExp, aes(x = log(gdpPercap) , y = lifeExp, col = continent)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ x, aes(group = 1), col = "black", se = F) +
theme_minimal()
La curva de color negro representa una curva logarítmica. Los puntos de datos siguen esta curva bastante de cerca, por lo tanto, es altamente recomendable emplear una transformación logarítmica en los predictores HIV.AIDS y gdpPercap. De los gráficos anteriores se observa que la relación lineal es más fuerte después de que estas variables se hayan transformado logarítmicamente.
ggplot(gapExp, aes(x = pop , y = lifeExp, col = continent)) +
geom_point() +
geom_smooth(method = "lm", aes(group = 1), col = "black", se = F) +
theme_minimal() 
ggplot(gapExp, aes(x = log(pop) , y = lifeExp, col = continent)) +
geom_point() +
geom_smooth(method = "lm", aes(group = 1), col = "black", se = F) +
theme_minimal()
ggplot(gapExp, aes(x = log(pop) , y = lifeExp, col = continent)) +
geom_point() +
geom_smooth(method = "lm", aes(group = 1), col = "black", se = F) +
theme_minimal()
ggplot(gapExp, aes(x = infant.deaths , y = lifeExp, col = continent)) +
geom_point() +
geom_smooth(method = "lm", aes(group = 1), col = "black", se = F) +
theme_minimal()
ggplot(gapExp, aes(x = log(infant.deaths + 1) , y = lifeExp, col = continent)) +
geom_point() +
geom_smooth(method = "lm", aes(group = 1), col = "black", se = F) +
theme_minimal()
A través de las visualizaciones, las transormaciones tienen un aspecto muy prometedor y parece mejorar la relación lineal de la variable de respuesta con los predictores. Para asegurar esto, se verifican los coeficientes de correlación antes y después de las transformaciones.
attach(gapExp)
N <- cor(gapExp[, c(3,5,6,12,17:19)])[5,1:7, drop = F]
corrplot(N, method = "number", cl.pos = "n")
gapExp$lHIV.AIDS <- log(gapExp$HIV.AIDS + 1)
gapExp$lgdpPercap <- log(gapExp$gdpPercap)
gapExp$lpop <- log(gapExp$pop)
gapExp$lMeasles <- log(gapExp$Measles + 1)
gapExp$linfant.deaths <- log(gapExp$infant.deaths + 1)
gapExp <- dplyr::select(gapExp, -c(HIV.AIDS, gdpPercap, pop, Measles, infant.deaths))
M<-cor(gapExp[, c(4,14:19)])[2,1:7, drop = F]
corrplot(M, method = "number", cl.pos = "n")
Se ve que el coeficiente de correlación aumentó para cada variable transformada logarítmicamente. Ahora existe una relación lineal más fuerte entre estos predictores y lifeExp. Exactamente lo que se buscaba.
4 Construcción de Modelos
4.1 Primera propuesta
Después de que preparar la información y verificar todas las suposiciones necesarias para construir un modelo de regresión aceptable, se comienza con la construcción y selección del mejor modelo. Con el método de selección hacia atrás se elege el mejor modelo comenzando con el modelo de mínimos cuadrados completo que contiene todos los predictores, y luego se van eliminando iterativamente los predictores menos útiles, uno a la vez, del subconjunto. La función regsubsets del paquete leaps permite realizar la selección paso a paso hacia hacia atrás, usando el argumento method = "backward". Centramos las variables para evitar el problema de multicolinealidad en el modelo y para incluir un término polinomial. Esto se debe a que es probable que los términos lineales y cuadrático tengan altas correlaciones.
centre <- function(x) { x - mean(x)}
gapExp.centered <- data.frame(gapExp, Polio = centre(gapExp$Polio),
Diphtheria = centre(gapExp$Diphtheria),
thinness..1.19.years = centre(gapExp$thinness..1.19.years),
thinness.5.9.years = centre(gapExp$thinness.5.9.years))back <- regsubsets(lifeExp ~ Adult.Mortality + Alcohol + percentage.expenditure + BMI + poly(Polio, 2) + Total.expenditure + poly(Diphtheria, 2) + poly(thinness..1.19.years, 2) + poly(thinness.5.9.years, 2) + Schooling + lHIV.AIDS + lgdpPercap + lpop + lMeasles + linfant.deaths, data = gapExp, method = "backward")
plot(back, scale = "Cp")
Para un modelo de mínimos cuadrados que contiene \(d\) predictores, la estimación de la estadística \(C_p\) de Mallow es un técnicas para ajustar el error cuadrático medio (MSE) de la prueba. Se calcula utilizando la ecuación:
\[C_p=\frac{1}{n}(RSS + 2dσ̂ ^2). \]
Aquí, \(σ̂ ^2\) es una estimación de la varianza del error \(ϵ\). La estadística \(C_p\) agrega una penalización de \(2dσ̂^2\) al RSS, donde la penalización aumenta a medida que aumenta el número de predictores en el modelo, por lo que al determinar cuál subconjunto de predictores es el mejor, se elige el modelo con el valor \(C_p\) más bajo.
En la gráfica anterior, podemos ver que el algoritmo selecciona las variables Alcohol, thinness..1.19.years, thinness.5.9.years, lHIV.AIDS, lgdpPercap, lpop, lMeasles, and linfant.deaths. El valor más bajo de \(C_p\) es \(20\) y en donde hay una cuadrícula negra en la gráfica, se incluyen los valores para el modelo por el método de selección hacia atrás.
Otra forma de visualizar esto, es através de una tibble, empleando la función tidy de la librería parsnip, una parte fundamental de la paquetería en tidymodels. tidy resume la información sobre los componentes de un modelo de acuerdo al argumento que pasemos en method dentro de la función regsubsets.
require(tidymodels)
tidy(back)## # A tibble: 8 x 24
## `(Intercept)` Adult.Mortality Alcohol percentage.expen… BMI `poly(Polio, 2)…
## <lgl> <lgl> <lgl> <lgl> <lgl> <lgl>
## 1 TRUE FALSE FALSE FALSE FALSE FALSE
## 2 TRUE FALSE FALSE FALSE FALSE FALSE
## 3 TRUE FALSE FALSE FALSE FALSE FALSE
## 4 TRUE FALSE FALSE FALSE FALSE FALSE
## 5 TRUE FALSE FALSE FALSE FALSE FALSE
## 6 TRUE FALSE TRUE FALSE FALSE FALSE
## 7 TRUE FALSE TRUE FALSE FALSE FALSE
## 8 TRUE FALSE TRUE FALSE FALSE FALSE
## # … with 18 more variables: poly(Polio, 2)2 <lgl>, Total.expenditure <lgl>,
## # poly(Diphtheria, 2)1 <lgl>, poly(Diphtheria, 2)2 <lgl>,
## # poly(thinness..1.19.years, 2)1 <lgl>, poly(thinness..1.19.years, 2)2 <lgl>,
## # poly(thinness.5.9.years, 2)1 <lgl>, poly(thinness.5.9.years, 2)2 <lgl>,
## # Schooling <lgl>, lHIV.AIDS <lgl>, lgdpPercap <lgl>, lpop <lgl>,
## # lMeasles <lgl>, linfant.deaths <lgl>, r.squared <dbl>, adj.r.squared <dbl>,
## # BIC <dbl>, mallows_cp <dbl>La ventaja de contar con esta tibble es que se pueden llamar a los distintos componentes del modelo, como la \(R²\) o el BIC (Criterio de Información Bayesiana).
Con esta depuración, se propone la nueva modelización para la base de datos. Se emplea el marco de referencia de la meta-paquetería tidymodels, pues es una colección de librerías para modelación y aprendizaje automático utilizando los principios del tidyverse.
require(dotwhisker)
lm_mod <- linear_reg() %>%
set_engine("lm")
lm_fit <- lm_mod %>%
fit(lifeExp ~ Alcohol + thinness..1.19.years + thinness.5.9.years + lHIV.AIDS +
lgdpPercap + lpop + lMeasles + linfant.deaths, data = gapExp)
broom::tidy(lm_fit)
broom::tidy(lm_fit) %>%
dwplot(dot_args = list(size = 2, color = "black"),
whisker_args = list(color = "black"),
vline=geom_vline(xintercept=0,colour='grey50',linetype=2))## # A tibble: 9 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 9.89 4.34 2.28 2.35e- 2
## 2 Alcohol -0.312 0.0807 -3.86 1.45e- 4
## 3 thinness..1.19.years -0.385 0.0925 -4.16 4.48e- 5
## 4 thinness.5.9.years 0.258 0.0884 2.92 3.85e- 3
## 5 lHIV.AIDS -5.52 0.324 -17.0 1.95e-42
## 6 lgdpPercap 2.51 0.329 7.65 5.85e-13
## 7 lpop 3.07 0.336 9.14 3.72e-17
## 8 lMeasles -0.276 0.0882 -3.13 1.97e- 3
## 9 linfant.deaths -3.55 0.396 -8.96 1.22e-16
Para obtener las métricas de rendimiento, se cuenta con la función glance:
glance(lm_fit)## # A tibble: 1 x 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.924 0.921 3.46 344. 2.77e-122 8 -623. 1267. 1301.
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>También se cuenta con la función vip, de la librería homónima, para graficar la importancia de la variable para cada predictor en el modelo. El valor de importancia se determina en función de las estadísticas \(F\) y los coeficientes de estimación en el objeto lm.
require(vip)
vip(lm_fit)
La resultante se ve muy bien y cada variable es estadísticamente significativa a un nivel de \(0.1\%\). La \(R^2\) ajustada es \(0.9211\). Esto significa que todas las variables de nuestro modelo están explicando el \(92.11\%\) de la variación en la variable de respuesta lifeExp.
car::vif(lm(lifeExp ~ Alcohol + thinness..1.19.years + thinness.5.9.years +
lHIV.AIDS + lgdpPercap + lpop + lMeasles + linfant.deaths,
data = gapExp))## Alcohol thinness..1.19.years thinness.5.9.years
## 1.914408 3.957565 3.812058
## lHIV.AIDS lgdpPercap lpop
## 1.962144 3.849195 5.288542
## lMeasles linfant.deaths
## 1.745723 8.907579Al verificar los valores de VIF, se tiene un problema con la multicolinealidad. Esto se debe a que algunas de las variables están por encima de \(5\). Esto es problemático porque los errores estándar de los predictores en la regresión podrían estar inflados, i.e., son más altos de lo que podrían ser sin un problema de multicolinealidad.
Se vuelven a analizar las visualizaciones para el objeto del tipo lm y así contar con un diagnóstico más robusto.
plot(lm(lifeExp ~ Alcohol + thinness..1.19.years + thinness.5.9.years +
lHIV.AIDS + lgdpPercap + lpop + lMeasles + linfant.deaths,
data = gapExp), c(1,2))
La q-q plot parece decente y que los residuales están distribuidos aproximadamente normalmente con una cola izquierda más gruesa que la distribución teórica. Esto es porque los puntos se desvían de la línea de la izquierda. La gráfica de residuales vs la ajustada se ve mejor que la de el modelo inicial, considerando todos los predictores. Sin embargo, la varianza constante no se puede suponer todavía porque podemos ver que la línea roja tiene forma de arco y, por lo tanto, sigue un patrón. Se busca lograr que la línea roja sea una línea recta horizontal en \(y = 0\) para poder asumir homocedasticidad.
La estadística PRESS nos da la suma de cuadrados del error residual predicho y mide el ajuste del modelo. Cuanto más bajo, mejor.
PRESS(lm(lifeExp ~ Alcohol + thinness..1.19.years + thinness.5.9.years +
lHIV.AIDS + lgdpPercap + lpop + lMeasles + linfant.deaths,
data = gapExp))## [1] 3114.46Se concluye no utilizar el modelo construido por selección hacia atrás debido a la alta multicolinealidad entre los predictores. Se observa que linfant.deaths tiene el valor VIF más alto de todos los predictores. Se centró esta variable y el valor de VIF sigue siendo alto. Por lo tanto, se descarta esta variable.
4.2 Segunda propuesta
back <- regsubsets(lifeExp ~ Adult.Mortality + Alcohol + percentage.expenditure + BMI + poly(Polio, 2) + Total.expenditure + poly(Diphtheria, 2) + poly(thinness..1.19.years, 2) + poly(thinness.5.9.years, 2) + Schooling + lHIV.AIDS + lgdpPercap + lpop + lMeasles, data = gapExp, method = "backward")
plot(back, scale = "Cp")
tidy(back)## # A tibble: 8 x 23
## `(Intercept)` Adult.Mortality Alcohol percentage.expen… BMI `poly(Polio, 2)…
## <lgl> <lgl> <lgl> <lgl> <lgl> <lgl>
## 1 TRUE FALSE FALSE FALSE FALSE FALSE
## 2 TRUE FALSE FALSE FALSE FALSE FALSE
## 3 TRUE FALSE FALSE FALSE FALSE FALSE
## 4 TRUE FALSE FALSE FALSE FALSE FALSE
## 5 TRUE FALSE FALSE FALSE FALSE FALSE
## 6 TRUE FALSE FALSE FALSE FALSE FALSE
## 7 TRUE TRUE FALSE FALSE FALSE FALSE
## 8 TRUE TRUE FALSE FALSE FALSE FALSE
## # … with 17 more variables: poly(Polio, 2)2 <lgl>, Total.expenditure <lgl>,
## # poly(Diphtheria, 2)1 <lgl>, poly(Diphtheria, 2)2 <lgl>,
## # poly(thinness..1.19.years, 2)1 <lgl>, poly(thinness..1.19.years, 2)2 <lgl>,
## # poly(thinness.5.9.years, 2)1 <lgl>, poly(thinness.5.9.years, 2)2 <lgl>,
## # Schooling <lgl>, lHIV.AIDS <lgl>, lgdpPercap <lgl>, lpop <lgl>,
## # lMeasles <lgl>, r.squared <dbl>, adj.r.squared <dbl>, BIC <dbl>,
## # mallows_cp <dbl>De la gráfica y la tibble anteriores, se ve que el algoritmo de selección hacia atrás eligió Adult.Mortality, Diftheria (el término lineal y el cuadrático), thinness..1.19.years, thinness.5.9.years, Schooling, lHIV.AIDS y lgdpPercap. Como antes, las cuadrículas en negro se incluyen en el modelo y los campos en blanco se excluyen. El mejor modelo tiene un valor \(C_p\) de \(7.3\), menor que para el modelo anterior.
lm_mod2 <- linear_reg() %>%
set_engine("lm")
lm_fit2 <- lm_mod2 %>%
fit(lifeExp ~ Adult.Mortality + poly(Diphtheria, 2) + thinness..1.19.years +
thinness.5.9.years + Schooling + lHIV.AIDS + lgdpPercap, data = gapExp)
broom::tidy(lm_fit2)
broom::tidy(lm_fit2) %>%
dwplot(dot_args = list(size = 2, color = "black"),
whisker_args = list(color = "black"),
vline=geom_vline(xintercept=0,colour='grey50',linetype=2))## # A tibble: 9 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 45.5 2.58 17.6 1.89e-44
## 2 Adult.Mortality -0.00652 0.00225 -2.89 4.21e- 3
## 3 poly(Diphtheria, 2)1 20.7 4.45 4.66 5.32e- 6
## 4 poly(Diphtheria, 2)2 18.1 4.36 4.15 4.67e- 5
## 5 thinness..1.19.years -0.391 0.0973 -4.02 7.81e- 5
## 6 thinness.5.9.years 0.242 0.0915 2.65 8.61e- 3
## 7 Schooling 0.505 0.138 3.67 3.00e- 4
## 8 lHIV.AIDS -6.04 0.345 -17.5 4.66e-44
## 9 lgdpPercap 2.48 0.372 6.66 2.01e-10
Una vez más observamos las visualizaciones de diagnóstico:
Y, las métricas de rendimiento:
## # A tibble: 1 x 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.915 0.912 3.65 306. 4.28e-117 8 -636. 1292. 1326.
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
La \(R^2\) ajustada es de alrededor de \(0.91\) y todos los predictores son estadísticamente significativos a un nivel del \(0.1\%\). La q-q plot se parece a al modelo anterior. La gráfica de residuales vs ajustada ha mejorado y se ha eliminado la forma del arco en la tendencia. La línea roja es casi horizontal ahora y los residuales se dispersan aleatoriamente alrededor de la línea \(y = 0\).
## GVIF Df GVIF^(1/(2*Df))
## Adult.Mortality 1.793482 1 1.339210
## poly(Diphtheria, 2) 1.925111 2 1.177915
## thinness..1.19.years 3.939165 1 1.984733
## thinness.5.9.years 3.673942 1 1.916753
## Schooling 4.321024 1 2.078707
## lHIV.AIDS 1.995428 1 1.412596
## lgdpPercap 4.427470 1 2.104155Todos los valores de VIF para este modelo están por debajo de \(5\), lo que indica que casi no hay multicolinealidad entre los predictores.
## [1] 3463.621Sin embargo, la estadística PRESS para el segundo modelo es más alta (peor) en comparación con el anterior.
En conclusión, se prefiere el segundo modelo debido al mejor aspecto de la gráfica de residuales vs ajustada y también debido a los valores de VIF más bajos. El valor \(C_p\) para el segundo modelo es \(7.3\), que es un valor más bajo en comparación a \(20\) para el primer modelo. Por lo tanto, considerando todo esto, se va con el segundo modelo como el mejor modelo.