Regresión Lineal Multiple
2025-04-05
Angel Felipe Ospina Beltrán
Jenny Katalina Beltrán Garzón
Juan Esteban Vargas Velandia
Juan Pablo Rodriguez Ariza
Es necesario activar los paquetes dplyr, MASS, ISLR, corrplot, psych, GGally.
library(dplyr)
library(MASS)
library(ISLR)
library(corrplot)
library(psych)
library(GGally)Descripción del dataset:
A continuación, se presenta una base de datos generada por la Inteligencia Artificial ChatGPT, la cual incluye variables económicas como el Producto Interno Bruto (PIB), la tasa de desempleo, la inversión y el Índice de Desarrollo Humano (IDH). El objetivo de este análisis es examinar las relaciones y correlaciones entre estas variables, y posteriormente construir un modelo de regresión lineal que permita predecir el comportamiento de una o varias de ellas a partir de las demás.
1. Exploración de los datos
datos_pib <- data.frame(
PIB = c(25.3, 28.1, 30.4, 22.5, 27.9, 24.2, 31.0, 26.5, 29.3, 23.8,
28.7, 30.1, 27.5, 25.0, 29.8, 24.9, 26.8, 31.5, 23.6, 27.0),
Desempleo = c(7.2, 6.8, 5.5, 8.1, 6.3, 7.5, 5.0, 6.9, 5.8, 7.8,
6.4, 5.6, 6.7, 7.3, 5.9, 7.6, 6.6, 4.9, 7.9, 6.5),
Inflacion = c(3.1, 2.8, 2.4, 3.5, 2.9, 3.2, 2.1, 3.0, 2.5, 3.3,
2.7, 2.3, 2.9, 3.4, 2.6, 3.2, 2.8, 2.0, 3.6, 2.7),
Inversion = c(520, 610, 750, 480, 690, 530, 800, 580, 710, 500,
660, 730, 600, 540, 740, 520, 590, 810, 490, 620),
IDH = c(0.820, 0.835, 0.850, 0.800, 0.830, 0.815, 0.860, 0.825, 0.845, 0.805,
0.840, 0.855, 0.830, 0.810, 0.850, 0.812, 0.828, 0.865, 0.802, 0.832))
head(datos_pib)## PIB Desempleo Inflacion Inversion IDH
## 1 25.3 7.2 3.1 520 0.820
## 2 28.1 6.8 2.8 610 0.835
## 3 30.4 5.5 2.4 750 0.850
## 4 22.5 8.1 3.5 480 0.800
## 5 27.9 6.3 2.9 690 0.830
## 6 24.2 7.5 3.2 530 0.815pairs(datos_pib,main="IDH",pch = 21,
bg = c("red", "green3", "blue","darkblue"))
2. Analizar relación entre variables
• cor() -> Calcula los coeficientes de correlación entre variables (solo acepta vectores numéricos)
• corrplot() -> Representación gráfica de los coeficientes de correlación
• ggpairs() -> Combina en un único gráfico diagramas de dispersión, distribución de las variables y los valores de correlación.
A continuación estudiamos la relación entre las variables para identificar cuáles pueden ser los mejores predictores o si hay alguna con una relación de tipo no lineal o detectar indicios de colinealidad.
mat_cor <- round(cor(subset(datos_pib, select = ), method = "pearson"), digits = 3)
mat_cor## PIB Desempleo Inflacion Inversion IDH
## PIB 1.000 -0.979 -0.961 0.974 0.987
## Desempleo -0.979 1.000 0.974 -0.984 -0.981
## Inflacion -0.961 0.974 1.000 -0.949 -0.984
## Inversion 0.974 -0.984 -0.949 1.000 0.969
## IDH 0.987 -0.981 -0.984 0.969 1.000Valores de correlación r próximos a 1 o -1 indican una alta correlación entre variables. También podemos representarlos gráficamente con la siguiente función:
corrplot(round(cor(subset(datos_pib, select =)), digits = 3), type = "lower")
require(psych)
ggplot(datos_pib, aes(x = PIB)) +
geom_histogram(aes(y = ..density..), binwidth = 1.5, fill = "grey70", color = "black", alpha = 0.8) +
stat_function(fun = dnorm,
args = list(mean = mean(datos_pib$PIB), sd = sd(datos_pib$PIB)),
color = "red", size = 1.2, linetype = "dashed") +
geom_density(color = "blue", size = 1.2) +
labs(
title = "Distribución del PIB",
x = "PIB",
y = "Densidad"
) +
theme_minimal(base_size = 14)
x <- seq(0, max(datos_pib$PIB), length = 40)
x## [1] 0.0000000 0.8076923 1.6153846 2.4230769 3.2307692 4.0384615
## [7] 4.8461538 5.6538462 6.4615385 7.2692308 8.0769231 8.8846154
## [13] 9.6923077 10.5000000 11.3076923 12.1153846 12.9230769 13.7307692
## [19] 14.5384615 15.3461538 16.1538462 16.9615385 17.7692308 18.5769231
## [25] 19.3846154 20.1923077 21.0000000 21.8076923 22.6153846 23.4230769
## [31] 24.2307692 25.0384615 25.8461538 26.6538462 27.4615385 28.2692308
## [37] 29.0769231 29.8846154 30.6923077 31.5000000f <- dnorm(x, mean = mean(datos_pib$PIB), sd = sd(datos_pib$PIB))
f## [1] 6.417662e-24 1.306646e-22 2.429305e-21 4.124294e-20 6.393818e-19
## [6] 9.051369e-18 1.170069e-16 1.381185e-15 1.488797e-14 1.465421e-13
## [11] 1.317142e-12 1.081050e-11 8.102189e-11 5.545008e-10 3.465335e-09
## [16] 1.977568e-08 1.030531e-07 4.903810e-07 2.130834e-06 8.454906e-06
## [21] 3.063453e-05 1.013577e-04 3.062286e-04 8.448466e-04 2.128400e-03
## [26] 4.896342e-03 1.028569e-02 1.973052e-02 3.456104e-02 5.528131e-02
## [31] 8.074453e-02 1.076939e-01 1.311633e-01 1.458736e-01 1.481441e-01
## [36] 1.373837e-01 1.163401e-01 8.996356e-02 6.352537e-02 4.096104e-02El siguiente paquete permite combinar los diagramas de dispersión, la distribución y los valores de correlación:
require(GGally)
ggpairs(datos_pib, lower = list(continuous = "smooth"),
diag = list(continuous = "bar"), axisLabels = "none")
De lo analizado hasta ahora podemos concluir que:
Los Valores de correlación r próximos a 1 o -1 indican una alta correlación entre las variables, en este caso absolutamente todas las variables económicas tienen una correlacion por encima del 0,96 en valor absoluto, si la correlacion es positiva la relacion es directa y si la correlacion es negativa la relacion es inversa.
Las variables mas fuertemente correlacionadas son IDH-PIB-(r:0,987),IDH-Desempleo-(r:-0,981),IDH-Inflacion-(r:-0,984) e Inflacion-Inversion-(r:-0,984).
Notandose que el Indice de desarrollo humano esta más fuertemente relacionado con la mayoria de las variables, sea de forma directa o inversa.
3. Se genera el modelo
modelo.lineal <- lm(IDH ~ . , data = datos_pib)
summary(modelo.lineal)##
## Call:
## lm(formula = IDH ~ ., data = datos_pib)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0032159 -0.0008964 -0.0002676 0.0008314 0.0037049
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.589e-01 4.779e-02 15.878 8.68e-11 ***
## PIB 3.916e-03 9.708e-04 4.033 0.001083 **
## Desempleo 2.378e-03 4.342e-03 0.548 0.591941
## Inflacion -2.147e-02 5.158e-03 -4.163 0.000833 ***
## Inversion 1.697e-05 2.886e-05 0.588 0.565252
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002143 on 15 degrees of freedom
## Multiple R-squared: 0.9905, Adjusted R-squared: 0.988
## F-statistic: 393.1 on 4 and 15 DF, p-value: 5.517e-15De lo analizado hasta ahora podemos concluir que:
El modelo con todas las variables introducidas como predictores es capaz de explicar el 99,05% de la varianza observada en el índice de desarrollo humano. Esto se refleja en un R² ajustado de 0,988, lo cual indica un excelente nivel de ajuste del modelo.
El p-value del modelo es significativo (5.517x10-15), lo que significa que existe una relación entre al menos una de las variables predictoras y el IDH, por lo tanto, el modelo es útil para explicar o predecir el comportamiento del IDH.
Los predictores que parecen tener una relación estadísticamente significativa con la variable respuesta son: PIB e inflación, dado que, ambos tienen valores p menores a 0.01, lo que indica que su relación con el IDH es significativa al 99% de confianza.
4. Selección de predictores
. step(modelo) -> Selección de predictores por stepwise y AIC
. update() -> Reajustar y actualizar el modelo
. confint(modelo) -> Intervalos de confianza para los coeficientes del modelo
En este ejemplo vamos a utilizar el método del stepwise mixto y el AIC con la función step( ) para determinar la calidad del modelo:
step(modelo.lineal, direction = "both", trace = 1)## Start: AIC=-241.58
## IDH ~ PIB + Desempleo + Inflacion + Inversion
##
## Df Sum of Sq RSS AIC
## - Desempleo 1 1.3770e-06 7.0243e-05 -243.19
## - Inversion 1 1.5880e-06 7.0454e-05 -243.13
## <none> 6.8866e-05 -241.58
## - PIB 1 7.4692e-05 1.4356e-04 -228.89
## - Inflacion 1 7.9553e-05 1.4842e-04 -228.22
##
## Step: AIC=-243.19
## IDH ~ PIB + Inflacion + Inversion
##
## Df Sum of Sq RSS AIC
## - Inversion 1 3.8600e-07 7.0630e-05 -245.08
## <none> 7.0243e-05 -243.19
## + Desempleo 1 1.3770e-06 6.8866e-05 -241.58
## - PIB 1 7.3901e-05 1.4414e-04 -230.81
## - Inflacion 1 1.0967e-04 1.7991e-04 -226.38
##
## Step: AIC=-245.08
## IDH ~ PIB + Inflacion
##
## Df Sum of Sq RSS AIC
## <none> 7.0630e-05 -245.08
## + Inversion 1 3.8600e-07 7.0243e-05 -243.19
## + Desempleo 1 1.7600e-07 7.0454e-05 -243.13
## - Inflacion 1 1.1711e-04 1.8774e-04 -227.52
## - PIB 1 1.6630e-04 2.3693e-04 -222.87##
## Call:
## lm(formula = IDH ~ PIB + Inflacion, data = datos_pib)
##
## Coefficients:
## (Intercept) PIB Inflacion
## 0.778706 0.003993 -0.019948El PIB y desempleo (la variable con mayor p-value) ha sido la única variable excluida en el proceso de selección. Reajustamos el modelo excluyendo dichas variables:
modelo.lineal1 <- update(modelo.lineal, formula = ~ . -PIB )
modelo.lineal2 <- update(modelo.lineal1, formula = ~ . -Desempleo )
summary(modelo.lineal2)##
## Call:
## lm(formula = IDH ~ Inflacion + Inversion, data = datos_pib)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0047007 -0.0019732 -0.0002435 0.0022102 0.0048716
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.694e-01 2.546e-02 34.143 < 2e-16 ***
## Inflacion -2.805e-02 4.697e-03 -5.973 1.51e-05 ***
## Inversion 6.580e-05 1.989e-05 3.308 0.00416 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002912 on 17 degrees of freedom
## Multiple R-squared: 0.9802, Adjusted R-squared: 0.9779
## F-statistic: 421.2 on 2 and 17 DF, p-value: 3.297e-15Los intervalos de confianza para cada uno de los coeficientes serían:
confint(modelo.lineal2)## 2.5 % 97.5 %
## (Intercept) 8.156504e-01 0.9230946050
## Inflacion -3.796056e-02 -0.0181429552
## Inversion 2.383353e-05 0.0001077624Cada una de las pendientes de un modelo de regresión lineal múltiple (coeficientes parciales de regresión de los predictores) se define del siguiente modo: Si el resto de variables se mantienen constantes, por cada unidad que aumenta el predictor en cuestión, la variable (Y) varía en promedio tantas unidades como indica la pendiente. Para este ejemplo, por cada unidad que aumenta el predictor inflación, las IDH disminuyen en promedio 0.019948 unidades, manteniéndose constantes el resto de predictores.
5. Validación de condiciones
Relación lineal entre los predictores numéricos y la variable respuesta:
Uno de los supuestos fundamentales en la regresión lineal múltiple es que exista una relación lineal entre cada predictor numérico y la variable dependiente. Este supuesto puede evaluarse mediante gráficos de dispersión entre los residuos del modelo y cada uno de los predictores.
Si la relación entre las variables es efectivamente lineal, los residuos deberían distribuirse aleatoriamente alrededor de cero, sin mostrar patrones sistemáticos y manteniendo una variabilidad constante a lo largo del eje X. Esto también permite identificar posibles datos atípicos o no linealidades.
Aunque los diagramas de dispersión entre la variable dependiente y los predictores pueden ser útiles en el análisis exploratorio inicial, los gráficos de residuos vs predictores son una herramienta más precisa para validar este supuesto en el contexto del modelo ajustado.
library(ggplot2)
library(gridExtra)
names(datos_pib)## [1] "PIB" "Desempleo" "Inflacion" "Inversion" "IDH"plot1 <- ggplot(data = datos_pib, aes(Inflacion, modelo.lineal2$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = datos_pib, aes(Inversion, modelo.lineal2$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2)
qqnorm(modelo.lineal2$residuals,col="black")
qqline(modelo.lineal2$residuals,col="blue")
Al representar los residuos frente a los valores ajustados por el modelo, se espera que los residuos se distribuyan de forma aleatoria alrededor de cero, sin seguir ningún patrón definido. Además, deben mantener una variabilidad aproximadamente constante a lo largo del eje X. Si, en cambio, se observa algún patrón sistemático (como una forma cónica, mayor dispersión en los extremos o agrupamientos), esto sugiere que la varianza de los errores no es constante, lo que implica una violación del supuesto de homocedasticidad.
ggplot(data = datos_pib, aes(modelo.lineal2$fitted.values, modelo.lineal2$residuals)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()
par(mfrow=c(1,1))
plot(modelo.lineal2)
# Detección y visualización de observaciones influyentes.
require(car)
par(mfrow=c(1,1))
influencePlot(modelo.lineal2)
## StudRes Hat CookD
## 5 -1.4892916 0.3323460 0.34342017
## 10 -1.8455526 0.1266150 0.14418579
## 15 2.0981902 0.2369767 0.37975535
## 19 0.5334249 0.2517012 0.03330496library(ggplot2)
# Gráfico residuos estudentizados frente a valores ajustados por el modelo
ggplot(data = datos_pib, aes(x = predict(modelo.lineal2),
y = abs(rstudent(modelo.lineal2))))+
geom_hline(yintercept = 3, color = "grey", linetype = "dashed")+
# Se identifican en rojo las observaciones con residuos estandarizados absolutos > 3
geom_point(aes(color = ifelse(abs(rstudent(modelo.lineal2)) > 3, "red", "black")))+
scale_color_identity()+
labs(title = "Distribución de los residuos estudentizados",
x = "Predicción modelo",
y = "Residuos estudentizados")+
theme_bw() + theme(plot.title = element_text(hjust = 0.5))
# Detección de los residuos estudentizados > 3 considerados como outliers
which(rstudent(modelo.lineal2) > 3)## named integer(0)outlierTest(modelo.lineal2)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 15 2.09819 0.05212 NA# Test de hipótesis para el análisis de normalidad de los residuos
shapiro.test(modelo.lineal2$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo.lineal2$residuals
## W = 0.96201, p-value = 0.5847# Test de contraste de homocedasticidad Breusch-Pagan
library(lmtest)
bptest(modelo.lineal2)##
## studentized Breusch-Pagan test
##
## data: modelo.lineal2
## BP = 0.30319, df = 2, p-value = 0.8593corrplot(cor(dplyr::select(datos_pib, Desempleo, Inflacion,
IDH)), method = "number", type = "lower")
A partir del análisis de correlación entre las variables predictoras, se observa una posible presencia de colinealidad. Esta situación ocurre cuando dos o más variables independientes están altamente correlacionadas entre sí, lo que puede afectar la estabilidad de los coeficientes del modelo.
Para confirmar esta sospecha, se recomienda calcular el Factor de Inflación de la Varianza (VIF), que permite cuantificar el grado de colinealidad existente en el modelo.
# Factores de inflación de la varianza
vif(modelo.lineal2)## Inflacion Inversion
## 10.01535 10.01535Hasta el momento podemos concluir que:
Luego de revisar los graficos de dispersion del modelo, se evidencia que el supuesto de linealidad se cumple ya que los residuos se distribuyen aleatoriamente alrededor de 0 y su variabilidad de los residuos parece constante por lo tanto la relacion de los predictores y la variable dependiente es aproximadamente lineal.
Los gráficos de residuos frente a los predictores (Inflación e Inversión) no muestran patrones sistemáticos evidentes, lo que sugiere que la relación entre los predictores y la variable dependiente es aproximadamente lineal. Esto valida el supuesto fundamental de que el modelo puede describirse mediante una combinación lineal de los predictores.
El test de Shapiro ayuda a confirmar si los residuos siguen una distribución normal. Un valor-p mayor a 0.05 respalda la hipótesis de normalidad; en caso contrario, se debe considerar aplicar transformaciones o usar modelos robustos.
El gráfico de residuos ajustados muestra una dispersión bastante uniforme en torno a cero, sin forma cónica ni patrones visibles. No obstante, el test de Breusch-Pagan permite confirmar estadísticamente si la varianza de los errores es constante.
Para que el modelo cumpla las condiciones, no deben haber outliers en los residuos estudentializados. Para el modelo lineal multiple que explica el IDH, es posible apreciar que el outlier test con valor P de Bonferroni no es lo suficientemente significativo para considerar la existencia de outliers en el modelo, haciendo que este cumpla la condicion.
Otro supuesto fundamental para el funcionamiento del modelo de regresión lineal multiple es la no colinealidad. El test VIF permite analizar la colinealidad entre las variables, si el resultado del test es mayor a 10 indica una alta colinealidad. Al ejecutar el VIF test para el modelo lineal es posible ver que el supuesto de colinealidad se rompe para los datos, pues para todas las variables existe esta caracteristica.
6. Conclusión
Si bien el modelo de regresión lineal multiple encontrado pasa todas las pruebas luego de eliminar las variables que no son relevantes para el modelo, queda en duda la aplicabilidad en la magnitud de los coeficientes predictores debido a la existencia de colinealidad. No obstante, es posible destacar la capacidad del modelo para explicar la variabilidad de los datos prorpocionados, pues el modelo solo falla con uno de los test aplicados.