Regresión lineal múltiple
20/abril/2022
Ejemplo 1
En este ejemplo realizaremos un análisis básico de un modelo de regresión lineal múltiple. Los datos corresponden al cuadro 3.1 de las notas. El problema planteado es el siguiente:
Un analista hace un estudio químico y espera que el rendimiento de cierta sustancia se vea afectado por dos factores. Por experimentos similares, se sabe que los factores x1 y x2 no están relacionados, por ello el analista decide utilizar un modelo de regresión lineal múltiple. Calcule el modelo de regresión y grafiquelo sobre las observaciones.
library("scatterplot3d")
#Datos del cuadro 3.1
y<-c(251.3,251.3, 248.3, 267.5, 273.0, 276.5, 270.3, 274.9 , 285.0 , 290.0 , 297.0 , 302.5 , 304.5 , 309.3 , 321.7 , 330.7 , 349.0 )
x1<-c(41.9 , 43.4 , 43.9 , 44.5 , 47.3 , 47.5 , 47.9 , 50.2 , 52.8 , 53.2 ,56.7 , 57.0 , 63.5 , 64.3, 71.1 ,77.0 , 77.8)
x2<-c( 29.1 ,29.3 , 29.5,29.7 ,29.9 , 30.3 , 30.5 , 30.7, 30.8,30.9 ,31.5,31.7,31.9 , 32.0 ,32.1 , 32.5 , 32.9)Una vez que tenemos los datos, procedemos a graficar y ajustar el modelo
#Graficar
scatterplot3d(x1,x2,y,type = "h", color = "blue", angle=55, pch = 16)
#Modelo
model<-lm(y~x1+x2)
#Graficar el plano
scatterplot3d(x1,x2,y, color = "blue", angle=40, pch = 16)$plane3d(model)
¿Qué puede decir sobre el modelo a partir de la gráfica?, ¿será buena idea que el investigador se guie por las experiencias similares para asumir el cumplimiento de los supuestos?
Veamos un resumen:
summary(model)##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.0119 -4.0282 -0.3633 4.2258 8.3784
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -151.7580 98.6751 -1.538 0.14635
## x1 1.2539 0.3878 3.233 0.00601 **
## x2 12.0008 3.8480 3.119 0.00755 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.431 on 14 degrees of freedom
## Multiple R-squared: 0.9688, Adjusted R-squared: 0.9643
## F-statistic: 217.4 on 2 and 14 DF, p-value: 2.875e-11¿Que información podemos extraer de la consulta anterior?, ¿es posible ver indicios de algún problema a partir de los resultados anteriores?, si la respuesta es si, ¿qué podría hacer para comprobarlo?.
Ejemplo 2
Para el siguiente ejemplo trabajaremos con datos provenientes de un paquete de R.
library("datarium")
#Datos
data("marketing")
?marketing
head(marketing)## youtube facebook newspaper sales
## 1 276.12 45.36 83.04 26.52
## 2 53.40 47.16 54.12 12.48
## 3 20.64 55.08 83.16 11.16
## 4 181.80 49.56 70.20 22.20
## 5 216.96 12.96 70.08 15.48
## 6 10.44 58.68 90.00 8.64Observe que si queremos usar todas las variables disponibles para describir el comportameinto de ‘sales’, no es posible graficar. ¿Qué podría hacer en este caso para verificar los supuestos básicos sobre los datos?
Ajustemos el modelo y veamos los resultados principales
modelo1 <- lm(sales ~ youtube + facebook + newspaper, data = marketing)
summary(modelo1)##
## Call:
## lm(formula = sales ~ youtube + facebook + newspaper, data = marketing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.5932 -1.0690 0.2902 1.4272 3.3951
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.526667 0.374290 9.422 <2e-16 ***
## youtube 0.045765 0.001395 32.809 <2e-16 ***
## facebook 0.188530 0.008611 21.893 <2e-16 ***
## newspaper -0.001037 0.005871 -0.177 0.86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.023 on 196 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8956
## F-statistic: 570.3 on 3 and 196 DF, p-value: < 2.2e-16¿Puede identificar algún problema en el modelo?, ¿los cambios en ‘newspaper’ son significativos?, ¿qué podria hacer para solucionar este problema?
Una posible solución es simplemente eliminar la variable que no es significativa
modelo2 <- lm(sales ~ youtube + facebook, data = marketing)
summary(modelo2)##
## Call:
## lm(formula = sales ~ youtube + facebook, data = marketing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.5572 -1.0502 0.2906 1.4049 3.3994
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.50532 0.35339 9.919 <2e-16 ***
## youtube 0.04575 0.00139 32.909 <2e-16 ***
## facebook 0.18799 0.00804 23.382 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.018 on 197 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8962
## F-statistic: 859.6 on 2 and 197 DF, p-value: < 2.2e-16Ahora que tenemos solamente dos variables explicativas, es posible graficar
scatterplot3d(marketing$youtube,marketing$facebook, marketing$sales, color = "darkblue", angle=60, pch = 16)$plane3d(modelo2)
¿Qué puede decir sobre este nuevo modelo?. Si tuviera que elegir, ¿cual de los dos modelos sería mejor y por qué?, ¿qué otro modelo podría intentar?.
Extras
Finalmente, para encontrar intervalos de confianza, predicciones, pruebas de hipotesis, etc., es posible usar las funciones que ya conocemos, por ejemplo:
#Intervalos de confianza
confint(model)## 2.5 % 97.5 %
## (Intercept) -363.3950078 59.879008
## x1 0.4220968 2.085689
## x2 3.7476056 20.254018También hay que mencionar que en la mayoría de los casos tendremos un número grande de variables, por lo que escribir la suma de todas ella no es óptimo. Para esos casos podemos usar un punto que representa a todas las variables de la base:
model <- lm(sales ~., data = marketing)Eliminar variables de forma fácil:
#Modelo usando todos menos 1 predictor
model1 <- lm(sales ~. -newspaper, data = marketing)O actualizar un modelo que ya teniamos
model2 <- update(model, ~. -newspaper)Para estos ultimos modelos, verifique que utilizan las variables correctas comparándolos con los del segundo ejemplo.