EJERCICIO: Modelos de Regresión Múltiple
D. S. Fernandez del Viso
2022-10-18
Datos
Usaremos los datos bmi en el archivo
reg-lin-mul-BMI.csv. Estos son datos de individuos adultos
entre 21 y 79 años, con las siguientes variables: BMI, índice
de masa corporal (\(kg/m^{2}\));
Age, edad (años); Cholesterol, niveles de
colesterol en sangre, (\(mg/dL\));
Glucose, niveles de glucosa en la sangre, (\(mg/dL\)).
chol <- read.csv("data/reg-lin-mul-BMI.csv")
head(chol)
## BMI Age Cholesterol Glucose
## 1 19.283 21 178 95
## 2 24.542 57 250 98
## 3 24.738 46 176 102
## 4 47.868 47 171 105
## 5 44.220 61 222 101
## 6 29.881 74 156 72
Estadísticas descriptivas
Exploración con estadísticos descriptivos. Nos ayuda a detectar
posibles errores en los datos y tener una idea de las características de
la población muestreada.
## BMI Age Cholesterol Glucose
## Min. :19.28 Min. :21.00 Min. :119.0 Min. : 72.0
## 1st Qu.:25.10 1st Qu.:35.50 1st Qu.:163.0 1st Qu.: 91.5
## Median :29.59 Median :48.50 Median :183.0 Median :100.0
## Mean :30.11 Mean :48.41 Mean :186.5 Mean :117.6
## 3rd Qu.:34.34 3rd Qu.:60.75 3rd Qu.:204.0 3rd Qu.:107.8
## Max. :47.87 Max. :76.00 Max. :283.0 Max. :345.0
hist(chol$Glucose, breaks = 20, xlim = c(50,200))
Matriz de Correlaciones
Para explorar posible multicolinealidad, y eliminar variables muy
correlacionadas entre sí o considerarlas para interacción.
library(Hmisc)
#matriz de correlación
rm_matrix_corr <- as.matrix(chol)
rcorr(rm_matrix_corr)
## BMI Age Cholesterol Glucose
## BMI 1.00 0.19 0.28 0.22
## Age 0.19 1.00 0.17 0.27
## Cholesterol 0.28 0.17 1.00 0.08
## Glucose 0.22 0.27 0.08 1.00
##
## n= 58
##
##
## P
## BMI Age Cholesterol Glucose
## BMI 0.1614 0.0326 0.0952
## Age 0.1614 0.2029 0.0427
## Cholesterol 0.0326 0.2029 0.5373
## Glucose 0.0952 0.0427 0.5373
Gráfica de matriz de
regresión
#gráfica de regresiones en parejas, con línea de regresión
library(car)
## Loading required package: carData
scatterplotMatrix(chol, ~ BMI + Cholesterol + Glucose + Age,
smooth = list(lty = 2), id = TRUE,
regLine = list(lty = 1, col = "red"),
col = "blue")
## Warning in applyDefaults(legend, defaults = list(coords = NULL, pt.cex = cex, :
## unnamed legend arguments, will be ignored
Modelo de Regresión Múltiple Básico
Para escoger la variable dependiente debemos tener alguna evidencia
de la relación causa-efecto entre las variables. Conocer el sistema
estudiado o usar la literatura es lo sugerido en esta etapa. Vamos a
usar la investigación de Laclaustra y col. (2018) para formular el
modelo; según este trabajo el colesterol (LDL) está relacionado con el
BMI, siendo el colestero la la variable dependiente. Añadiremos las
otras variables y empezaremos con un modelo básico completo, sin
interacciones.
rm1 <- lm(Cholesterol ~ ., data = chol)
summary(rm1)
##
## Call:
## lm(formula = Cholesterol ~ ., data = chol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.107 -22.062 -3.733 17.464 76.962
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 133.315499 23.183220 5.751 4.27e-07 ***
## BMI 1.348919 0.696178 1.938 0.0579 .
## Age 0.269970 0.297106 0.909 0.3676
## Glucose -0.004541 0.080281 -0.057 0.9551
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.54 on 54 degrees of freedom
## Multiple R-squared: 0.09333, Adjusted R-squared: 0.04295
## F-statistic: 1.853 on 3 and 54 DF, p-value: 0.1486
## [1] 577.9346
Examen de supuestos
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
##
## recode
## The following objects are masked from 'package:Hmisc':
##
## src, summarize
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# para multicolinealidad
car::vif(rm1) %>%
knitr::kable()
| BMI |
1.071069 |
| Age |
1.096912 |
| Glucose |
1.113279 |
# otros
par(mfrow=c(2,2))
plot(rm1)
Multicolinealidad
Evaluada utilizando la función vif (“Variance Inflation
Factor”), que indica variables que tienen multicolinealidad cuando su
vif es mucho mayor que 1.
Multicollinearity
in Regression Analysis: Problems, Detection, and Solutions
Linealidad
Gráfica Residuals vs Fitted
Normalidad
Gráfica Q-Q
Homocedasticidad
Gráfica Scale-Location
Valores extremos
Gráfica Residuals vs Leverage
What
is a Residuals vs. Leverage Plot?
Autocorrelación
(Independencia)
Prueba Durbin-Watson La prueba de Durbin-Watson nos
permite evaluar si ocurre autocorrelación entre los valores residuales
de la variable dependiente. Por ejemplo, una variable que depende del
tiempo, presenta autocorrelación. La \(H_0\) es que la autocorrelación en los
residuales del modelo es 0.
library(car)
dbt1 <- durbinWatsonTest(rm1, simulate = TRUE)
dbt1
## lag Autocorrelation D-W Statistic p-value
## 1 -0.1965377 2.379202 0.152
## Alternative hypothesis: rho != 0
Variación del Modelo Básico, rm2
Eliminando Glucose
# modelo rm2
rm2 <- lm(Cholesterol ~ BMI + Age, data = chol)
summary(rm2)
##
## Call:
## lm(formula = Cholesterol ~ BMI + Age, data = chol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -59.977 -22.301 -3.873 17.606 77.149
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 133.1879 22.8631 5.825 3.07e-07 ***
## BMI 1.3418 0.6784 1.978 0.053 .
## Age 0.2660 0.2861 0.930 0.357
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.24 on 55 degrees of freedom
## Multiple R-squared: 0.09327, Adjusted R-squared: 0.0603
## F-statistic: 2.829 on 2 and 55 DF, p-value: 0.06771
## [1] 575.9381
# supuestos
# multicolinealidad
car::vif(rm2) %>%
knitr::kable()
| BMI |
1.035969 |
| Age |
1.035969 |
# linealidad, homocedasticidad, normalidad, extremos
par(mfrow=c(2,2))
plot(rm2)
# autocorrelación
dbt2 <- durbinWatsonTest(rm2, simulate = TRUE)
dbt2
## lag Autocorrelation D-W Statistic p-value
## 1 -0.1978905 2.382569 0.122
## Alternative hypothesis: rho != 0
Variación del Modelo Básico, rm3
Interacción Glucose * Age
# modelo rm2
rm3 <- lm(Cholesterol ~ BMI + Age:Glucose, data = chol)
summary(rm3)
##
## Call:
## lm(formula = Cholesterol ~ BMI + Age:Glucose, data = chol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -58.710 -23.227 -3.879 17.417 76.748
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.420e+02 2.072e+01 6.856 6.5e-09 ***
## BMI 1.395e+00 6.940e-01 2.010 0.0493 *
## Age:Glucose 4.087e-04 1.127e-03 0.363 0.7182
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.46 on 55 degrees of freedom
## Multiple R-squared: 0.08122, Adjusted R-squared: 0.04781
## F-statistic: 2.431 on 2 and 55 DF, p-value: 0.09735
## [1] 576.704
# supuestos
# multicolinealidad
car::vif(rm3) %>%
knitr::kable()
## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif
| BMI |
1.069743 |
| Age:Glucose |
1.069743 |
# linealidad, homocedasticidad, normalidad, extremos
par(mfrow=c(2,2))
plot(rm3)
# autocorrelación
dbt3 <- durbinWatsonTest(rm3, simulate = TRUE)
dbt3
## lag Autocorrelation D-W Statistic p-value
## 1 -0.2181471 2.426426 0.098
## Alternative hypothesis: rho != 0
Variación del Modelo Básico, rm4
Interacción BMI * Age
# modelo rm4
rm4 <- lm(Cholesterol ~ BMI + Age + Age:BMI, data = chol)
summary(rm4)
##
## Call:
## lm(formula = Cholesterol ~ BMI + Age + Age:BMI, data = chol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.336 -19.759 -3.487 24.347 74.836
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.49401 69.44551 0.468 0.6417
## BMI 4.89352 2.41136 2.029 0.0474 *
## Age 2.35477 1.39125 1.693 0.0963 .
## BMI:Age -0.07271 0.04742 -1.533 0.1310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32.84 on 54 degrees of freedom
## Multiple R-squared: 0.1311, Adjusted R-squared: 0.08283
## F-statistic: 2.716 on 3 and 54 DF, p-value: 0.05365
## [1] 575.4662
# supuestos
# multicolinealidad
car::vif(rm4) %>%
knitr::kable()
## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif
| BMI |
13.40860 |
| Age |
25.09798 |
| BMI:Age |
43.86474 |
# linealidad, homocedasticidad, normalidad, extremos
par(mfrow=c(2,2))
plot(rm4)
# autocorrelación
dbt4 <- durbinWatsonTest(rm4, simulate = TRUE)
dbt4
## lag Autocorrelation D-W Statistic p-value
## 1 -0.1980423 2.36812 0.15
## Alternative hypothesis: rho != 0
PROCEDIMIENTO “STEPWISE”
#formulación de un modelo nulo y un modelo completo
modNulo <- lm(Cholesterol ~ 1, data = chol)
modFull <- lm(Cholesterol ~ Age + BMI + Glucose + BMI:Age + Age:Glucose + BMI:Glucose,
data = chol)
#procedimiento stepwise
cholstep <- step(modNulo,
scope = list(lower=modNulo, upper=modFull,
direction="both"))
## Start: AIC=411.02
## Cholesterol ~ 1
##
## Df Sum of Sq RSS AIC
## + BMI 1 5294.7 61710 408.25
## <none> 67004 411.02
## + Age 1 1928.9 65076 411.33
## + Glucose 1 457.9 66547 412.62
##
## Step: AIC=408.25
## Cholesterol ~ BMI
##
## Df Sum of Sq RSS AIC
## <none> 61710 408.25
## + Age 1 954.9 60755 409.34
## + Glucose 1 29.6 61680 410.22
## - BMI 1 5294.7 67004 411.02
##
## Call:
## lm(formula = Cholesterol ~ BMI, data = chol)
##
## Coefficients:
## (Intercept) BMI
## 142.528 1.459
Modelo seleccionado con
stepwise
# modelo step
step <- lm(Cholesterol ~ BMI, data = chol)
summary(step)
##
## Call:
## lm(formula = Cholesterol ~ BMI, data = chol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.211 -24.396 -2.752 18.375 75.331
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 142.5275 20.5132 6.948 4.22e-09 ***
## BMI 1.4593 0.6658 2.192 0.0326 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.2 on 56 degrees of freedom
## Multiple R-squared: 0.07902, Adjusted R-squared: 0.06257
## F-statistic: 4.805 on 1 and 56 DF, p-value: 0.03255
AICstep <- AIC(step)
AICstep
## [1] 574.8426
# supuestos
# multicolinealidad no se puede correr con menos de 2 términos
# linealidad, homocedasticidad, normalidad, extremos
par(mfrow=c(2,2))
plot(step)
# autocorrelación
dbtstep <- durbinWatsonTest(step, simulate = TRUE)
dbtstep
## lag Autocorrelation D-W Statistic p-value
## 1 -0.2150966 2.41667 0.112
## Alternative hypothesis: rho != 0
Métodos gráficos de selección de modelos
Comparación de coeficientes
library(ggplot2)
library(coefplot)
#cálculo para todos los modelos
modbmi1 <- lm(Cholesterol ~ ., data=chol)
modbmi2 <- lm(Cholesterol ~ BMI + Age, data=chol)
modbmi3 <- lm(Cholesterol ~ BMI + Glucose, data=chol)
modbmi4 <- lm(Cholesterol ~ Age + Glucose, data=chol)
modbmi5 <- lm(Cholesterol ~ BMI, data = chol)
#comparando coeficientes de todos los modelos
multiplot(modbmi1, modbmi2, modbmi3, modbmi4, modbmi5, pointSize = 2, intercept=FALSE)
Selección
usando R-cuadrado ajustado y Mallow’s Cp para mejores modelos
Otro método visual de selección, basado en el \(R^{2}-ajustado\) y en Cp (Mallow’s Cp):
library(leaps)
modR2Cp <- regsubsets(Cholesterol ~ Age + BMI + Glucose + BMI:Glucose + Age:Glucose + Age:BMI, data = chol, nbest = 2)
plot(modR2Cp, scale="adjr2")
plot(modR2Cp, scale="Cp")
Interpretación
El mejor modelo con R-cuadrado ajustado está en la parte superior, y
contiene Age, BMI y su interacción.
Para el caso del estadístico Cp, al igual que con el AIC, el valor
más bajo representa el mejor modelo, en este caso uno de BMI solamente.
En general el valor de Cp se acerca al modelo con un número similar de
parámetros (incluyendo el intercepto), en nuestro caso dos.