Métodos estadísticos I

Tarea 5

Para este caso trabajaremos sobre problemas de colinealidad. Utilizaremos datos vistos en clase, los ordenamos con base en la variable respuesta considerando los pasos que continuaran:

datos_ridge <- read.csv("datos ridge.csv")
datos_ridge <- datos_ridge[order(datos_ridge$y), ]

Para buscar correlaciones podemos utilizar la matriz de corrleaciones:

cor(datos_ridge[, 1:3])

##         x1      x2      x3
## x1  1.0000  0.2236 -0.9582
## x2  0.2236  1.0000 -0.2402
## x3 -0.9582 -0.2402  1.0000

Observemos que la variables 1 y la variable 3 tienen una alta correlación. Sin embargo podemos obtener otras formas para detectar colinealidad, por ejemplo el indicador VIF, primero debemos ajustar el modelo lineal de mínimos cuadrados ordinarios:

library(car)

## Warning: package 'car' was built under R version 3.0.3

ols <- lm(y ~ ., data = datos_ridge)
sqrt(vif(ols))

##    x1    x2    x3 
## 3.496 1.030 3.511

Obsevemos que para la variable 1 y 3 tenemos un valor del VIF mayor a 3, por lo cual concluimos que existe colinealidad entre las varibles.

La primer alternativa es una regresion Ridge, para la cual debemos de selecciónar un \( \lambda \) adecuado, utilizamos el paquete “MASS”, y probamos \( \lambda \)'s desde 0 hasta 0.1 con una diferencia entre cada lambda de 0.001:

library(MASS)
select(lm.ridge(y ~ ., data = datos_ridge, lambda = seq(0, 0.1, 0.001)))

## modified HKB estimator is 0.1391 
## modified L-W estimator is 0.1162 
## smallest value of GCV  at 0.1

A partir de los resultados escogemos un valor de 0.1 para \( \lambda \). Con este valor ajustamos el modelo.

Una segunda alternativa consiste en utilizar los componentes principales de la matriz de variables predictoras, como nuevas variables regresoras, por definición los vectores resultantes son independientes. Una vez ajustado este modelo podemos comparar los estimadores de los coeficientes de regresión de los tres modelos:

ridge <- lm.ridge(y ~ ., data = datos_ridge, lambda = 0.1)
componentes <- lm(datos_ridge$y ~ princomp(datos_ridge[, 1:3])$scores)

# comparamos los coefficientes:

coeficientes <- matrix(nrow = 3, ncol = 4, dimnames = (list(c("OLS", "Ridge", 
    "PC"), c("Intercepto", "v1", "v2", "v3"))))
coeficientes[1, ] <- coef(ols)
coeficientes[2, ] <- coef(ridge)[1:4]
coeficientes[3, ] <- coef(componentes)

coeficientes

##       Intercepto     v1      v2     v3
## OLS      -121.27 0.1269  0.3482 -19.02
## Ridge    -110.77 0.1189  0.3454 -38.41
## PC         36.11 0.1395 -0.3489 -19.02

Vamos a hacer algo más, comparemos los residuales de los modelos:

plot(datos_ridge$y,ols$residuals,type='p',ylab='Residuales',ylab='') abline(h=0,lty=2) lines(datos_ridge$y,componentes$residuals,pch=16) lines(datos_ridge$y,datos_ridge$y-as.matrix(M)%*%as.matrix(coef(ridge)))

M <- cbind(rep(1, nrow(datos_ridge)), datos_ridge[, 1:3])
plot(datos_ridge$y, ols$residuals, type = "p", ylab = "Residuales", xlab = "y")
abline(h = 0, lty = 2)
lines(datos_ridge$y, componentes$residuals, pch = 16, col = 2)
lines(datos_ridge$y, datos_ridge$y - as.matrix(M) %*% as.matrix(coef(ridge)), 
    col = 4)
legend("topleft", col = c(2, 4), bty = "n", legend = c("Ridge", "Componentes\n      Principales"), 
    lty = 1)

En la gráfica anterior los puntos muestran los residuales del modelo lineal de mínimos cuadrados ordinarios, notemos que son prácticamente idénticos a los del modelo de componentes principales; y la diferencia con respecto al modelo de regresión Ridge. Ahora, si tomamos en la raiz del error cuadrado medio tenemos:

Conclusión: No me gustan. Y el de componentes principales arroja resultados muy similares al de OLS, pero la interpretación es más complicada.