MODELOS DE RESPUESTA CUANTITATIVA Modelos de regresión Ridge y Lasso
Curso: Modelos Lineales II
Alumno : Alvaro Cesar Tauma Salvador (20181160)
1 MODELOS DE RESPUESTA CUANTITATIVA Modelos de regresión Ridge y Lasso
Ejercicio1
Se tiene información de la producción minera del Perú para cinco productos en el periodo de 2001 -2018.
Datos en archivo MLG_PT_Semana_ 05 _ 01 csv
1. Formule el MLG en términos del estudio
2. Ajuste los datos de la producción minera a una regresión lineal múltiple Realice al análisis de la multicolinealidad
3.Ajuste los datos de la producción minera a una regresión lineal múltiple Ridge Use la estimación sesgada por mínimos cuadrados
4.Realice la comparación de ambos estimaciones
1.1 MODELO DE REGRESIÓN RIDGE
1.1.1 1. Formule el MLG en términos del estudio
Variable respuesta: Y= Producción total minera (TM)
Variables predictoras:
X1=Producción cobre (TM)
X2=Producción zinc (TM)
X3=Producción oro (Krgs.)
X4=Producción plata (Krgs.)
X5=Producción plomo (TM)
a .Componente aleatorio:
$Yi∼N(μi,σ^2),E(Y)=μ $
b .Componente sistematico (Predictor Lineal)
\(ni=X'β=β0+β1X1+β2X2+β3X3+β4X4+β5X5\)
c .Funcion de enlace : Distribución:Normal , Función Lineal: identidad
\(E(Y)=g(μ)=1(μi)=ηi=X'β = β0+β1X1+β2X2+β3X3+β4X4+β5X5\)
\(E(Y)=X'β\)
d .Modelo de Regresion Lineal Multiple
\(Y=β0+β1X1+β2X2+β3X3+β4X4+β5X5\)
#Entrada de datos
Datos <- read.csv("MLII_PT_Clase_04_Datos_01.csv",sep = ";")
head(Datos)## X1_Cobre X2_Zinc X3_Oro X4_Plata X5_Plomo Y_Total
## 1 345620 493408 160 44815 187248 4798673
## 2 331407 479285 206 42553 201191 5014445
## 3 363688 543590 142 46904 213448 4953285
## 4 328020 581753 177 55921 223727 4051263
## 5 370004 548705 197 58973 208485 4863254
## 6 401327 612903 224 64647 240950 4694722X = model.matrix (Y_Total~ . , Datos) [, -1 ] # separar la ultima columna
head(X)## X1_Cobre X2_Zinc X3_Oro X4_Plata X5_Plomo
## 1 345620 493408 160 44815 187248
## 2 331407 479285 206 42553 201191
## 3 363688 543590 142 46904 213448
## 4 328020 581753 177 55921 223727
## 5 370004 548705 197 58973 208485
## 6 401327 612903 224 64647 240950y = Datos$Y_Total# Ajuste a la regresión lineal múltiple (mínimos cuadrados)
Modelo=lm(Y_Total~ X1_Cobre+X2_Zinc+X3_Oro+X4_Plata+X5_Plomo, Datos)
summary(Modelo) #obtenemos los coef. estimados##
## Call:
## lm(formula = Y_Total ~ X1_Cobre + X2_Zinc + X3_Oro + X4_Plata +
## X5_Plomo, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1292889 -422945 75776 459001 1628524
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.812e+06 9.702e+05 4.959 2.08e-05 ***
## X1_Cobre 3.712e+00 7.337e-01 5.059 1.55e-05 ***
## X2_Zinc 1.729e+00 1.230e+00 1.406 0.1692
## X3_Oro -5.132e+00 9.590e+00 -0.535 0.5961
## X4_Plata 1.444e+00 6.365e-01 2.268 0.0300 *
## X5_Plomo -1.330e+01 5.597e+00 -2.377 0.0234 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 702200 on 33 degrees of freedom
## Multiple R-squared: 0.9789, Adjusted R-squared: 0.9757
## F-statistic: 306.8 on 5 and 33 DF, p-value: < 2.2e-161.1.2 2. Ajuste los datos de la producción minera a una regresión lineal múltiple Realice al análisis de la multicolinealidad
# Análisis de la multicolinealidad
n=dim(X)[1]
n # cantidad de observaciones ## [1] 39p=dim(X)[2]
p # cantidad de variables predictoras ## [1] 5# Examinar matriz de correlaciones
cor(X) #vemos que as correlaciones son altas## X1_Cobre X2_Zinc X3_Oro X4_Plata X5_Plomo
## X1_Cobre 1.0000000 0.8444156 0.7438040 0.8840808 0.6775012
## X2_Zinc 0.8444156 1.0000000 0.9199822 0.9505124 0.8684625
## X3_Oro 0.7438040 0.9199822 1.0000000 0.9609244 0.8229567
## X4_Plata 0.8840808 0.9505124 0.9609244 1.0000000 0.7987709
## X5_Plomo 0.6775012 0.8684625 0.8229567 0.7987709 1.0000000Interpretación
Se puede observar que existe una alta correlación significativa entre las variables predictoras Se puede concluir, que existe una multicolinealidad en el modelo de regresión .
#Factores de inflación de variancia (VIF)
1/(1-summary(lm(X1_Cobre~X2_Zinc+X3_Oro+X4_Plata+X5_Plomo,Datos))$ r.squared)## [1] 15.305511/(1-summary(lm(X2_Zinc ~X1_Cobre+X3_Oro+X4_Plata+X5_Plomo,Datos))$ r.squared)## [1] 16.131341/(1-summary(lm(X3_Oro~X1_Cobre+X2_Zinc+X4_Plata+X5_Plomo,Datos))$ r.squared)## [1] 49.127131/(1-summary(lm(X4_Plata~X1_Cobre+X2_Zinc+X3_Oro+X5_Plomo,Datos))$ r.squared)## [1] 100.07951/(1-summary(lm(X5_Plomo~X1_Cobre+X2_Zinc+X3_Oro+X4_Plata,Datos))$ r.squared)## [1] 5.195703library(faraway)
vif(Modelo)## X1_Cobre X2_Zinc X3_Oro X4_Plata X5_Plomo
## 15.305512 16.131344 49.127134 100.079505 5.195703si el vif > 10 => hay multicolinealidad, la mayoria supera el calor de 10, por lo tanto podemos decir que hay una multicolinalidad alta ntra las v.v predictoras.
Interpretación
El análisis de los factores de inflación de variancia, resultan que los cuatro primeros superan el valor de 10 y sólo el quinto es menor a 10 Se puede concluir, que existe una alta multicolinealidad entre las variable predictoras
# Examinar valores propios
XX=t(X)%*%X
Lambda=eigen(XX)$values; Lambda## [1] 3.157506e+14 9.414956e+12 2.410785e+12 3.716266e+10 5.406394e+09# Calcular el índice k
k=max(Lambda)/min(Lambda); k## [1] 58403.17#si k>1000 entonces exites multicolinealidadLambda## [1] 3.157506e+14 9.414956e+12 2.410785e+12 3.716266e+10 5.406394e+09# Calcular el índice k por predictor
for (i in 1:p) {cat("kj: ",i, " = ", max(Lambda)/Lambda[i],"\n") }## kj: 1 = 1
## kj: 2 = 33.53713
## kj: 3 = 130.9742
## kj: 4 = 8496.448
## kj: 5 = 58403.17Interpretación
El análisis de los valores propios de la matriz (resultaron con valores altos Se tiene un valor del índice k superior a 1000 lo que indica una severa multicolinealidad Los índices k para cada predictor resultaron los dos últimos mayores a 1000 lo que indica una moderada multicolinealidad.
1.1.3 3.Ajuste los datos de la producción minera a una regresión lineal múltiple Ridge Use la estimación sesgada por mínimos cuadrados
library(glmnet) ## Loading required package: Matrix## Loaded glmnet 4.1-7##### Ajuste a la regresi?n Ridge ####
Modelo_R=glmnet(X, y, alpha=0) #alpha=0 Ridge ; alpha=1 LAsso .................
plot(Modelo_R, xvar="lambda", label=TRUE)
#muestra como al aumenar lambda, los coef de reg tienden a cero# Estimaciones Ridge para varios valores de lambda
M1<-glmnet(X,y,alpha=0,lambda=0.9)
M2<-glmnet(X,y,alpha=0,lambda=10)
M3<-glmnet(X,y,alpha=0,lambda=100)
M4<-glmnet(X,y,alpha=0,lambda=1000)
M5<-glmnet(X,y,alpha=0,lambda=10000)
cbind(coef(M1),coef(M2),coef(M3),coef(M4),coef(M5))## 6 x 5 sparse Matrix of class "dgCMatrix"
## s0 s0 s0 s0
## (Intercept) 4.829986e+06 4.829961e+06 4.829943e+06 4.829031e+06
## X1_Cobre 3.772705e+00 3.772835e+00 3.774899e+00 3.792681e+00
## X2_Zinc 1.761267e+00 1.761359e+00 1.762684e+00 1.774369e+00
## X3_Oro -4.212231e+00 -4.210142e+00 -4.177659e+00 -3.895680e+00
## X4_Plata 1.383774e+00 1.383628e+00 1.381410e+00 1.362060e+00
## X5_Plomo -1.357937e+01 -1.357970e+01 -1.358658e+01 -1.364286e+01
## s0
## (Intercept) 4.813067e+06
## X1_Cobre 3.917936e+00
## X2_Zinc 1.860630e+00
## X3_Oro -1.826502e+00
## X4_Plata 1.220047e+00
## X5_Plomo -1.400434e+01# Nos muestra que a medida que aumenta el lambda los estimadores salen mas pequeños.# Seleccionando el mejor lambda (óptimo)
set.seed(1230)
# Usando la validaci?n cruzada
cv.Ridge=cv.glmnet(X, y, alpha=0, nfolds=10, type.measure="mse")
#nfolds para indicar la separación, type.measure "el tipo de medida "
plot(cv.Ridge)
cv.Ridge$lambda.min# Lambda que consigue el mínimo test-error## [1] 429069.1# Lambda optimo para el test-error
cv.Ridge$lambda.1se## [1] 622505.8# Estimaci?n de coeficientes para el lambda óptimo
Modelo_R_Final=glmnet(X, y, alpha=0, lambda= cv.Ridge$lambda.1se) # aqui esta usando el lambda optimo sacado de
coef(Modelo_R_Final)## 6 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 3.278525e+06
## X1_Cobre 3.463234e+00
## X2_Zinc 2.179510e+00
## X3_Oro 5.694463e+00
## X4_Plata 7.077771e-01
## X5_Plomo -6.613563e+00options(scipen=999)
coef(Modelo_R_Final)## 6 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 3278525.3177302
## X1_Cobre 3.4632342
## X2_Zinc 2.1795101
## X3_Oro 5.6944632
## X4_Plata 0.7077771
## X5_Plomo -6.6135626MODELO ESTIMADO POR REGRESION RIDGE : \[ \hat{y_{i}}= 3278525+ 3.463234*X_{1i}+2.179510*X_{2i}+ 5.694463*X_{3i}+0.7077771*X_{4i}-6.613563*X_{5i}\]
1.1.4 4.Realice la comparación de ambos estimaciones
coef(Modelo_R_Final)## 6 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 3278525.3177302
## X1_Cobre 3.4632342
## X2_Zinc 2.1795101
## X3_Oro 5.6944632
## X4_Plata 0.7077771
## X5_Plomo -6.6135626coef(Modelo)## (Intercept) X1_Cobre X2_Zinc X3_Oro X4_Plata
## 4811788.685430 3.712085 1.728900 -5.132357 1.443709
## X5_Plomo
## -13.304929Modelo_R_Final##
## Call: glmnet(x = X, y = y, alpha = 0, lambda = cv.Ridge$lambda.1se)
##
## Df %Dev Lambda
## 1 5 96.85 622500 1 - cv.Ridge$cvm/var(y)## [1] -0.0189896994 -0.0122603384 -0.0099618985 -0.0090251632 -0.0079987987
## [6] -0.0068744097 -0.0056428509 -0.0042941674 -0.0028175323 -0.0012011807
## [11] 0.0005676586 0.0025028340 0.0046193436 0.0069334118 0.0094625885
## [16] 0.0122257107 0.0152431156 0.0185366424 0.0221296845 0.0260472464
## [21] 0.0303159793 0.0349642002 0.0400218901 0.0455206654 0.0514937181
## [26] 0.0579757166 0.0650026625 0.0726116925 0.0808408207 0.0897286104
## [31] 0.0993137690 0.1096346565 0.1207287029 0.1326317278 0.1453771592
## [36] 0.1589951538 0.1735116217 0.1889471661 0.2053159547 0.2226245518
## [41] 0.2408707286 0.2600423945 0.2801162110 0.3010570127 0.3228168371
## [46] 0.3453344619 0.3685352181 0.3923311816 0.4166217932 0.4412949389
## [51] 0.4662285031 0.4912923835 0.5163509245 0.5412657001 0.5658985487
## [56] 0.5901147370 0.6137861173 0.6367941320 0.6590325228 0.6804096167
## [61] 0.7008500801 0.7202973149 0.7387243017 0.7560873460 0.7723840448
## [66] 0.7876286425 0.8018452199 0.8150697592 0.8273478823 0.8387324808
## [71] 0.8492813472 0.8590549081 0.8681141471 0.8765204434 0.8843278787
## [76] 0.8915817575 0.8983413670 0.9046490302 0.9105427254 0.9160553176
## [81] 0.9212145874 0.9260417176 0.9305554455 0.9347734494 0.9387098932
## [86] 0.9423774184 0.9457822503 0.9489342486 0.9518429986 0.9545076373
## [91] 0.9569449961 0.9591642196 0.9611728013 0.9629811745 0.9646007380
## [96] 0.9660423877 0.9673188030 0.9684419230 0.9694230422 0.97020083391.2 MODELO DE REGRESIÓN LASSO
Ejercicio 2
Se tiene información de la producción minera del Perú para cinco productos en el periodo de 2001 2018 Datos en archivo MLG_PT_Semana_ 05 _ 01 csv.
1.Ajuste los datos de la producción minera a una regresión lineal múltiple Aplique los métodos de selección de variables
2. Ajuste los datos de la producción minera a una regresión lineal múltiple Lasso Encuentre el modelo óptimo de regresión Lasso que se ajuste a la producción minera
3. Realice la comparación de los modelos estimados
1.2.1 1 Ajuste los datos de la producción minera a una regresión lineal múltiple Aplique los métodos de selección de variables
Datos<-read.table("MLII_PT_Clase_04_Datos_01.csv",header=TRUE,sep=";")
head(Datos)## X1_Cobre X2_Zinc X3_Oro X4_Plata X5_Plomo Y_Total
## 1 345620 493408 160 44815 187248 4798673
## 2 331407 479285 206 42553 201191 5014445
## 3 363688 543590 142 46904 213448 4953285
## 4 328020 581753 177 55921 223727 4051263
## 5 370004 548705 197 58973 208485 4863254
## 6 401327 612903 224 64647 240950 4694722str(Datos)## 'data.frame': 39 obs. of 6 variables:
## $ X1_Cobre: int 345620 331407 363688 328020 370004 401327 399302 406325 316355 368168 ...
## $ X2_Zinc : int 493408 479285 543590 581753 548705 612903 604018 619209 498583 620956 ...
## $ X3_Oro : int 160 206 142 177 197 224 299 287 313 318 ...
## $ X4_Plata: int 44815 42553 46904 55921 58973 64647 64492 64334 53065 62120 ...
## $ X5_Plomo: int 187248 201191 213448 223727 208485 240950 201873 203960 161182 203034 ...
## $ Y_Total : int 4798673 5014445 4953285 4051263 4863254 4694722 4698445 4609905 3831735 4172608 ...# Ajuste a la regresión lineal múltiple (mínimos cuadrados)
Modelo=lm(Y_Total~ ., Datos)
summary(Modelo)##
## Call:
## lm(formula = Y_Total ~ ., data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1292889 -422945 75776 459001 1628524
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4811788.6854 970232.0055 4.959 0.0000208 ***
## X1_Cobre 3.7121 0.7337 5.059 0.0000155 ***
## X2_Zinc 1.7289 1.2300 1.406 0.1692
## X3_Oro -5.1324 9.5899 -0.535 0.5961
## X4_Plata 1.4437 0.6365 2.268 0.0300 *
## X5_Plomo -13.3049 5.5974 -2.377 0.0234 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 702200 on 33 degrees of freedom
## Multiple R-squared: 0.9789, Adjusted R-squared: 0.9757
## F-statistic: 306.8 on 5 and 33 DF, p-value: < 0.00000000000000022R2=summary(Modelo)$r.squared; R2## [1] 0.9789392library(leaps)
# Selección de variables: Medidas de bondad de ajuste
Modelo_S<-regsubsets(Y_Total~., data=Datos, nvmax=5) #numero máximo de variables que puedes incorporar en el modelonames(summary(Modelo_S))# Medidas para la selección de modelos## [1] "which" "rsq" "rss" "adjr2" "cp" "bic" "outmat" "obj"summary(Modelo_S)$adjr2 # nos da todos los r2 ajustados (seleccionar el mayor)## [1] 0.9274212 0.9720049 0.9754244 0.9762571 0.9757482summary(Modelo_S)$cp # (seleccionar el menor)## [1] 75.730327 8.556583 4.467214 4.286420 6.000000options(scipen=999)
# Identificar el modelo con mayor R2
which.max(summary(Modelo_S)$adjr2)## [1] 4coef(object=Modelo_S, id=which.max(summary(Modelo_S)$adjr2)) #el mejor modelo es con los 4## (Intercept) X1_Cobre X2_Zinc X4_Plata X5_Plomo
## 4922573.764230 4.039232 1.804452 1.127270 -14.558157max(summary(Modelo_S )$adjr2)## [1] 0.9762571Modelomayorr2=lm(Y_Total~ X1_Cobre+X2_Zinc+X4_Plata+X5_Plomo, Datos)
summary(Modelomayorr2)$adj.r.squared## [1] 0.9762571Interpretación
Seleccionar el mejor modelo usando la medida del R 2 ajustado
MEJOR MODELO USANDO LA MEDIDA DE \(R^2\) AJUSTADO :
\[ \hat{y_{i}}= 4922573.764230+ 4.039232*X_{1i}+1.804452*X_{2i}+ 1.127270*X_{4i}-14.558157*X_{5i}\]
Se identifica el modelo con el índice 4 .
1.2.2 2. Ajuste los datos de la producción minera a una regresión lineal múltiple Lasso Encuentre el modelo óptimo de regresión Lasso que se ajuste a la producción minera**
#Selección de variables: Metodos Forward, Backward y Stepwise
# Identificar el modelo con mayor R2 con el m?todo Backward
Modelo_B<-regsubsets(Y_Total~., data=Datos, nvmax=5, method="backward")
coef(object = Modelo_B, id=which.max(summary(Modelo_B)$adjr2))## (Intercept) X1_Cobre X2_Zinc X4_Plata X5_Plomo
## 4922573.764230 4.039232 1.804452 1.127270 -14.558157####Identificar el modelo con mayor R2 con el método Forward
Modelo_F<-regsubsets(Y_Total~.,data=Datos, nvmax=5, method="forward")
coef(object = Modelo_F, id=which.max(summary(Modelo_F)$adjr2))## (Intercept) X1_Cobre X2_Zinc X4_Plata X5_Plomo
## 4922573.764230 4.039232 1.804452 1.127270 -14.558157##### Ajuste a la regresión Lasso ####
library(glmnet)
X = model.matrix (Y_Total~ . , Datos)[,-1] # separar la última columna
y = Datos$Y_Total
Modelo_L=glmnet(X, y, alpha=1)# Seleccionando el major lambda (?ptimo)
set.seed(100)
cv.Lasso=cv.glmnet(X, y, alpha=1, nfolds=10, type.measure="mse")
cv.Lasso$lambda.1se## [1] 218573.7Modelo_L_Final=glmnet(X, y, alpha=1, lambda= cv.Lasso$lambda.1se)
coef(Modelo_L_Final)## 6 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 3066219.610251
## X1_Cobre 4.018070
## X2_Zinc .
## X3_Oro .
## X4_Plata 1.064968
## X5_Plomo .MODELO ESTIMADO POR LASSO: \[ \hat{y_{i}}= 3066219.610251+ 4.018070*X_{1i}+1.064968*X_{4i}\]