Regresión de crestas en R (paso a paso)

La regresión de crestas es un método que podemos utilizar para ajustar un modelo de regresión cuando la multicolinealidad está presente en los datos.

En pocas palabras, la regresión de mínimos cuadrados intenta encontrar estimaciones de coeficientes que minimicen la suma de los residuos al cuadrado (RSS):

RSS = Σ (y _i – ŷ _i ) 2

dónde:

• Σ : símbolo griego que significa suma
• y _i : El valor de respuesta real para la i- ^ésima observación
• ŷ _i : el valor de respuesta predicho basado en el modelo de regresión lineal múltiple

Por el contrario, la regresión de crestas busca minimizar lo siguiente:

RSS + λΣβ _j ²

donde j varía de 1 a p variables predictoras y lambda ≥ 0.

Este segundo término de la ecuación se conoce como penalización por contracción . En la regresión de crestas, seleccionamos un valor para λ que produce el MSE de prueba más bajo posible (error cuadrático medio).

Este tutorial proporciona un ejemplo paso a paso de cómo realizar la regresión de crestas en R.

Paso 1: cargue los datos

Para este ejemplo, usaremos el conjunto de datos integrado de R llamado mtcars . Usaremos hp como variable de respuesta y las siguientes variables como predictores:

• mpg
• peso
• maldita sea
• qsec

Para realizar la regresión de crestas, usaremos funciones del paquete glmnet . Este paquete requiere que la variable de respuesta sea ​​un vector y que el conjunto de variables predictoras sea de la clase data.matrix .

El siguiente código muestra cómo definir nuestros datos:

#define variable de respuesta
y <- mtcars $ hp

#definir matriz de variables predictoras
 x <- data.matrix (mtcars [, c ('mpg', 'wt', 'drat', 'qsec')])

Paso 2: ajuste el modelo de regresión de crestas

A continuación, usaremos la función glmnet () para ajustar el modelo de regresión de la cresta y especificaremos alpha = 0 .

Tenga en cuenta que establecer alfa igual a 1 es equivalente a usar Regresión de lazo y establecer alfa en algún valor entre 0 y 1 es equivalente a usar una red elástica.

También tenga en cuenta que la regresión de crestas requiere que los datos estén estandarizados de manera que cada variable predictora tenga una media de 0 y una desviación estándar de 1.

Afortunadamente glmnet () realiza automáticamente esta estandarización por usted. Si ya estandarizó las variables, puede especificar estandarizar = Falso .

biblioteca (glmnet)

#fit modelo de regresión de la cresta
modelo <- glmnet (x, y, alpha = 0 )

#Ver resumen del resumen del modelo
 (modelo)

          Modo de clase de longitud   
a0 100 -ninguno- numérico
beta 400 dgCMatrix S4     
df 100 -ninguno- numérico
dim 2 -ninguno- numérico
lambda 100 -ninguno- numérico
dev.ratio 100 -ninguno- numérico
nulldev 1 -ninguno- numérico
n pasa 1 -ninguno- numérico
jerr 1 -ninguno- numérico
offset 1 -ninguno- lógico
llamar 4 -ninguno- llamar   
nobs 1 -ninguno- numérico

Paso 3: elija un valor óptimo para Lambda

A continuación, identificaremos el valor lambda que produce el error cuadrático medio de prueba (MSE) más bajo mediante el uso de una validación cruzada de k veces .

Afortunadamente, glmnet tiene la función cv.glmnet () que realiza automáticamente una validación cruzada de k veces utilizando k = 10 veces.

# Realice una validación cruzada de k-fold para encontrar el valor lambda óptimo
 cv_model <- cv. glmnet (x, y, alpha = 0 )

# encontrar un valor lambda óptimo que minimice la prueba MSE
 best_lambda <- cv_model $ lambda . min
best_lambda

[1] 10.04567

#producir gráfico de prueba MSE por valor lambda
parcela (cv_model) 

El valor lambda que minimiza la prueba MSE resulta ser 10.04567 .

Paso 4: Analizar el modelo final

Por último, podemos analizar el modelo final producido por el valor lambda óptimo.

Podemos usar el siguiente código para obtener las estimaciones de coeficientes para este modelo:

# encontrar coeficientes del mejor modelo
 best_model <- glmnet (x, y, alpha = 0 , lambda = best_lambda)
coef (mejor_modelo)

Matriz dispersa 5 x 1 de clase "dgCMatrix"
                    s0
(Intercepción) 475.242646
mpg -3.299732
peso 19.431238
maldición -1.222429
qsec -17.949721

También podemos producir un diagrama de seguimiento para visualizar cómo cambiaron las estimaciones de coeficientes como resultado del aumento de lambda:

#produce diagrama de 
trazado de Ridge (modelo, xvar = " lambda ")

Por último, podemos calcular el R cuadrado del modelo sobre los datos de entrenamiento:

#utilice el mejor modelo ajustado para hacer predicciones
 y_predicted <- predict (model, s = best_lambda, newx = x)

#find SST y SSE 
sst <- sum ((y - mean (y)) ^ 2)
sse <- sum ((y_predicted - y) ^ 2)

#find R-Squared
 rsq <- 1 - sse / sst
rsq

[1] 0,7999513

El R cuadrado resulta ser 0,7999513 . Es decir, el mejor modelo fue capaz de explicar el 79,99% de la variación en los valores de respuesta de los datos de entrenamiento.

Puede encontrar el código R completo utilizado en este ejemplo aquí .

• https://r-project.org
• https://www.python.org/
• https://www.stata.com/