En este ejemplo analizaremos el cumplimiento de los supuestos vistos en clase para el modelo de regresión lineal múltiple usando simulaciones

library("scatterplot3d")
library(car)
## Loading required package: carData

Primero creamos un conjunto de datos que sabemos que si cumplen con los supuestos:

set.seed(16)
x<-1:10
y0<-3*x+2
plot(x,y0)

error1<-rnorm(10,0,1)
y<-y0+error1
plot(x,y)
model1<-lm(y~x)
abline(model1)

summary(model1)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5466 -0.5060  0.1112  0.7226  1.0245 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.0185     0.6502   3.105   0.0146 *  
## x             3.0210     0.1048  28.830 2.27e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9518 on 8 degrees of freedom
## Multiple R-squared:  0.9905, Adjusted R-squared:  0.9893 
## F-statistic: 831.2 on 1 and 8 DF,  p-value: 2.267e-09

Veamos qué pasa si los errores no tienen mendia cero:

error2<-rnorm(10,10,1)
y2<-y0+error2
plot(x,y2)
model2<-lm(y2~x)
abline(model2)

summary(model2)
## 
## Call:
## lm(formula = y2 ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9873 -0.6230  0.2964  0.8786  1.2055 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  12.7096     0.7837   16.22 2.10e-07 ***
## x             2.9358     0.1263   23.24 1.25e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.147 on 8 degrees of freedom
## Multiple R-squared:  0.9854, Adjusted R-squared:  0.9836 
## F-statistic: 540.3 on 1 and 8 DF,  p-value: 1.247e-08

¿qué estimador se ve afectado por la media de los errores y en qué cantidad?

Veamos ahora qué ocurre con el supuesto de multicolinealidad:

x1<-1:10
x2<-2*x1+rnorm(10,0,.1)
y<-10+3*x1+7*x2+rnorm(10,0,1)
scatterplot3d(x1,x2,y,type = "h", color = "blue", angle=80, pch = 16)

model<-lm(y~x1+x2)
summary(model)
## 
## Call:
## lm(formula = y ~ x1 + x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9512 -0.1892  0.1392  0.2570  0.5239 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.6052     0.4031  23.829 5.83e-08 ***
## x1           -0.7455     4.1750  -0.179  0.86335    
## x2            8.8884     2.0672   4.300  0.00357 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5439 on 7 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 
## F-statistic: 4.128e+04 on 2 and 7 DF,  p-value: 5.609e-15
scatterplot3d(x1,x2,y, color = "blue", angle=10, pch = 16)$plane3d(model)

Usando el código anterior, realice diferentes experimentos cambiando la varianza al generar x2. ¿Qué pasa cuando la varianza al generar x2 aumenta?

En el siguiente conjunto de datos las variables x1 y x2 no tienen relación:

x1<-c(1,1,1,2,2,2,3,3,3)
x2<-c(1,2,3,1,2,3,1,2,3)
y<-3+9*x1+10*x2+rnorm(9,0,3)
scatterplot3d(x1,x2,y,type = "h", color = "blue", angle=80, pch = 16)

model<-lm(y~x1+x2)
summary(model)
## 
## Call:
## lm(formula = y ~ x1 + x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0289 -0.9866 -0.2687  1.0369  5.8186 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.5913     3.5736   0.165 0.874008    
## x1            9.1375     1.2139   7.528 0.000285 ***
## x2           11.3240     1.2139   9.329  8.6e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.973 on 6 degrees of freedom
## Multiple R-squared:  0.9599, Adjusted R-squared:  0.9466 
## F-statistic: 71.85 on 2 and 6 DF,  p-value: 6.44e-05
scatterplot3d(x1,x2,y, color = "blue", angle=50, pch = 16)$plane3d(model)

Veamos otra forma de eliminar la multicolinealidad: El siguiente conjunto de datos contiene informacion sobre la densidad de hueso del femur

data<-read.csv("MulticollinearityExample.csv", header=TRUE)

Queremos un modelo que incluya el porcentaje de grasa, peso, actividad y el producto de grasa por peso ,el cual llamaremos interacción.

model<-lm(Femoral.Neck~X.Fat+Weight.kg+Activity+X.Fat*Weight.kg,data)
par(mfrow=c(2,2))
plot(model)

summary(model)
## 
## Call:
## lm(formula = Femoral.Neck ~ X.Fat + Weight.kg + Activity + X.Fat * 
##     Weight.kg, data = data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.178453 -0.042754 -0.006129  0.033937  0.186795 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.549e-01  1.317e-01   1.176  0.24274    
## X.Fat            5.571e-03  4.087e-03   1.363  0.17632    
## Weight.kg        1.447e-02  2.852e-03   5.073 2.19e-06 ***
## Activity         2.238e-05  7.276e-06   3.075  0.00281 ** 
## X.Fat:Weight.kg -2.142e-04  7.394e-05  -2.898  0.00476 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07051 on 87 degrees of freedom
## Multiple R-squared:  0.5623, Adjusted R-squared:  0.5422 
## F-statistic: 27.95 on 4 and 87 DF,  p-value: 6.242e-15

¿Qué indican las pruebas de hipótesis simples?

Recuerde que el VIF nos ayuda a detectar la presencia de multicolinealidad:

vif(model)
##           X.Fat       Weight.kg        Activity X.Fat:Weight.kg 
##       14.931555       33.948375        1.053005       75.059251

Considerando lor resultados de las pruebas de hipótesis y el VIF para este modelo, ¿qué puede inferir?

Centrar los valores es otro método que ayuda a reducir la multicolinearidad. Este procedimiento consiste en restar la media a cada una de las variables. En el conjunto de datos, las restas ya estan guardadas como .S, pero de no estarlo, solo habría que tomar cada columna y restarle su media.

model2<-lm(Femoral.Neck~X.Fat.S+Weight.S+Activity.S+X.Fat.S*Weight.S,data)
plot(model2)

summary(model2)
## 
## Call:
## lm(formula = Femoral.Neck ~ X.Fat.S + Weight.S + Activity.S + 
##     X.Fat.S * Weight.S, data = data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.178453 -0.042754 -0.006129  0.033937  0.186795 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       8.216e-01  9.734e-03  84.403  < 2e-16 ***
## X.Fat.S          -5.982e-03  1.928e-03  -3.103  0.00259 ** 
## Weight.S          8.348e-03  1.066e-03   7.829  1.1e-11 ***
## Activity.S        2.238e-05  7.276e-06   3.075  0.00281 ** 
## X.Fat.S:Weight.S -2.142e-04  7.394e-05  -2.898  0.00476 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07051 on 87 degrees of freedom
## Multiple R-squared:  0.5623, Adjusted R-squared:  0.5422 
## F-statistic: 27.95 on 4 and 87 DF,  p-value: 6.242e-15
vif(model2)
##          X.Fat.S         Weight.S       Activity.S X.Fat.S:Weight.S 
##         3.323870         4.745648         1.053005         1.991063