Regresión Ridge
Ricardo Rodriguez Ovilla.
3/9/2022
Regresión Lineal
Consideremos \(p\) características y \(n\) observaciones, la matriz de dimensión \(n \cdot p\). Sea \(Y\) la variable respuesta. El objetivo de toda regresión es explicar \(Y\) en términos de \(X\) a traves de una relación funcional \[ Y_i = f(X_i) \]
Partiendo del supuesto que \(X\) se relaciona de forma lineal con la variable \(Y\) surge el modelo de regresión lineal:
\[ Y_i = X_i \cdot B = \beta_1x_1 + \beta_2x_2 + ...+ \beta_px_p + \epsilon_i \]
En este caso los coeficientes \(\beta_j ,\ \ \ \ j = 1,...,p\) representan el tamaño del efecto de la covariable \(j\) sobre \(Y\), es decir por cada unida de cambio de la variable \(j\) (manteniendo las demás variables constantes) se observa un cambio de \(\beta_j\) en la variable respuesta.
En el modelo anterior se necesita especificar la adicion del valor aleatorio \(\epsilon_i\) tradicionalmente si \(\epsilon_i \backsim N(0,\sigma^{2})\) y \(\epsilon_i \perp \epsilon_j\) para \(j \neq i\) se tiene que:
\[cov(\epsilon_i,\epsilon_i) = \left\{ \begin{array}{lcc}\sigma^2 & si & i = j \\\\ 0&si & i \neq j \\\end{array}\right.\]
Considerar que \(\epsilon_i\) es una variable aleatoria también implica que \(Y_i\) lo es, es decir \(Y_i\) se distribuye de forma normal por lo que para calcular los parámetros basta con encontrar los primeros dos momentos.
\[ E(Y_i) = E( X_i \cdot B + \epsilon_i) = E( X_i \cdot B ) + E(\epsilon_i) = X_i \cdot B\]
\[Var(Y_i) = E[(Y_i-E(Y_i)^{2})] = E[Y_i^{2} -2Y_iE(Y_i)+E(Y_i)^{2}] = E(Y_i^{2})-[E(Y_i)]^{2} \] \[= E[(X_iB)^{2} + 2 \epsilon_i X_iB + (\epsilon_i)^{2}] - [E(Y_i)]^{2} \] \[ = E(\epsilon_i^{2}) = var(\epsilon_i)\] Por lo que se concluye que \(Y_i \backsim N( X_i \cdot B , var(\epsilon_i))\).
Podemos escribir la relacion lineal en forma matricial de la siguiente forma:
\[ Y = XB + \epsilon \] Donde \(\epsilon = (\epsilon_1,...,\epsilon_n)^{T}\).
Los parámetros desconocidos en la regresión lineal son \(B \ \ \ \ \sigma^{2}\) los cuales necesitan ser aprendidos de los datos. Partiendo de que \(Y_i \backsim N( X_i \cdot B , \sigma^{2})\) se plantea lo siguiente:
\[f_{Y_i}(y_i) = \frac{1}{\sqrt{2\pi\sigma^{2}}}e^{[-\frac{(Y_i-X_i \cdot B)^{2}}{2\sigma^2}]} \]
Tomando la función likelihood tenemos lo siguiente:
\[L(*|Y,X)=\prod_{i=1}^n\left(\sqrt{2 \pi} \sigma^2\right)^{-1} \exp \left[-\left(y_i-x_i \beta\right)^2 / 2 \sigma^2\right]\] \[= \sqrt{2 \pi \sigma^2}^{-n} \prod_{i=1}^n\exp \left[-\left(y_i-x_i \beta\right)^2 / 2 \sigma^2\right]\] Si a esta expresion tomamos el logaritmo natural obtenemos que:
\[ = log[\sqrt{2 \pi \sigma^2}^{-n} \prod_{i=1}^n\exp \left[-\left(y_i-x_i \beta\right)^2 / 2 \sigma^2\right]] \] \[ = -nlog(\sqrt{2 \pi \sigma^2}) + log[\prod_{i=1}^n\exp \left[-\left(y_i-x_i \beta\right)^2 / 2 \sigma^2\right]] \] \[ = -nlog(\sqrt{2 \pi \sigma^2}) + \sum_{i=1}^{n}log(\exp \left[-\left(y_i-x_i \beta\right)^2 / 2 \sigma^2\right]) = -nlog(\sqrt{2 \pi \sigma^2}) + \sum_{i=1}^{n}-(y_i-x_i \beta)^2 / 2 \sigma^2 \] \[ log[L(*|Y,X)] = -nlog(\sqrt{2 \pi \sigma^2}) -[2\sigma^{2}]^{-1} \sum_{i=1}^{n}(y_i-x_i \beta)^2 \] Para obtener la expresión de \(\beta\):
\[\frac{\partial}{\partial\beta}(-nlog(\sqrt{2 \pi \sigma^2}) -[2\sigma^{2}]^{-1} \sum_{i=1}^{n}(y_i-x_i \beta)^2) \] \[ = -[2\sigma^{2}]^{-1}\frac{\partial}{\partial\beta} 2 \sum_{i=1}^{n}(y_i-x_i \beta)\sum_{i=1}^{n}(-x_i) \] \[0 = X^{T}(Y-XB) = X^{T}Y - X^{T}XB \Rightarrow X^{T}Y = X^{T}XB \]
Pre multiplicando por \((X^{T}X)^{-1}\):
\[B = (X^{T}X)^{-1} X^{T}Y \]
Poder multiplicar por \((X^{T}X)^{-1}\) exige que dicha matriz inversa exista.
Para \(\sigma^2\) :
\[log[L(*|Y,X)] = -nlog(\sqrt{2 \pi \sigma^2}) -[2\sigma^{2}]^{-1} \sum_{i=1}^{n}(y_i-x_i \beta)^2\] \[\Rightarrow \frac{\partial}{\partial\sigma} (-nlog(\sqrt{2 \pi \sigma^2}) -[2\sigma^{2}]^{-1} \sum_{i=1}^{n}(y_i-x_i \beta)^2) \]
\[ -n(\frac{\sqrt(2\pi)}{\sqrt(2\pi\sigma)})-\frac{1}{2} \sum_{i=1}^{n}(y_i-x_i \beta)^2(-2\sigma^{3}) = \frac{-n}{\sigma} + \sum_{i=1}^{n}(y_i-x_i \beta)^2(\sigma^{-3}) \] igualando a cero se obtiene que:
\[ 0 = -n\sigma^{-1} + \sum_{i=1}^{n}(y_i-x_i \beta)^2(\sigma^{-3}) = -n\sigma^{2}+ \sum_{i=1}^{n}(y_i-x_i \beta)^2\] \[ \Rightarrow \sigma^{2} = \frac{\sum_{i=1}^{n}(y_i-x_i \beta)^2}{n} \]
Las siguientes propiedades se incluyen sin demostración pero son sencillas de hacer:
\[ E(B) = B\] \[ Var(B) = \sigma^2(X^TX)^{-1}\]
Si la matriz \(X\) posee covariables altamente correlacionadas linealmente entonces es imposible separar la contribución individual de las variables.Esto se refleja en la regresion lineal con incertidumbre en el calculo de los parámetros.
Por ejemplo consideremos la siguiente base en la cual existen 100 variables predictoras altamente correlacionadas y la variable objetivo que indica el contenido en grasa medido por técnicas químicas.
Matriz de correlaciones
# Predictores
pred = datos %>% select( V11,V22,V33,V44,V87,V12,V17,V35)
cor = cor(pred)
ggcorrplot(cor, hc.order = TRUE, type = "lower",
lab = TRUE)
Para comprobar metricas se dividiran los datos en dos conjuntos de la siguiente forma:
set.seed(777)
pred = datos %>% select( V11,V22,V33,V44,V87,V12,V17,V35,fat)
id_train <- sample(1:nrow(pred), size = 0.7*nrow(pred), replace = FALSE)
datos_train <- pred[id_train, ]
datos_test <- pred[-id_train, ]
modelo <- lm(fat ~ ., data = datos_train)
summary(modelo)##
## Call:
## lm(formula = fat ~ ., data = datos_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.5478 -1.4735 -0.0275 1.8595 6.1188
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.451 3.183 2.027 0.044578 *
## V11 -12007.838 2392.001 -5.020 1.53e-06 ***
## V22 505.540 442.984 1.141 0.255715
## V33 -1682.722 450.604 -3.734 0.000272 ***
## V44 -158.843 21.549 -7.371 1.31e-11 ***
## V87 3.549 4.184 0.848 0.397748
## V12 15631.968 3093.214 5.054 1.32e-06 ***
## V17 -4209.664 1065.143 -3.952 0.000122 ***
## V35 1920.771 352.167 5.454 2.15e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.607 on 141 degrees of freedom
## Multiple R-squared: 0.9575, Adjusted R-squared: 0.9551
## F-statistic: 397.2 on 8 and 141 DF, p-value: < 2.2e-16# Coeficientes del modelo
# ==============================================================================
df_coeficientes <- modelo$coefficients %>%
enframe(name = "predictor", value = "coeficiente")
t = summary(modelo)$coefficients[,4]
t = ifelse(t>0.05,"No significante","")
df_coeficientes = cbind(df_coeficientes , t)
df_coeficientes %>%
filter(predictor != "(Intercept)") %>%
ggplot(aes(x = predictor, y = coeficiente, fill = predictor)) +
geom_col() +
labs(title = "Coeficientes") +
geom_text(aes(label = t), vjust = -1, colour = "black") +
theme_bw() +
theme(axis.text.x = element_text(size = 5, angle = 45))
Regresión Ridge
Para matrices \((X^{T}X)\) con problemas de singularidad se propuso el siguiente “fix” que se añade a la ecuación normal:
\[B(\lambda) = (X^{T}X + \lambda I)^{-1} X^{T}Y \]
donde \(\lambda \ \ \ \epsilon \ \ \ [0,\infty)\) lo que esto producirá distintas regresiones lineales. Este escalar es también conocido como parámetro de penalización
Si definimos el operador \(W(\lambda) = (X^{T}X + \lambda I)^{-1} X^{T}X\) entonces podemos poner el estimador ridge de la siguiente forma:
\[W(\lambda)B = W(\lambda) (X^{T}X )^{-1} X^{T}Y=(X^{T}X + \lambda I)^{-1} X^{T}X(X^{T}X)^{-1}X^{T}Y\] \[= (X^{T}X + \lambda I)^{-1}X^{T}Y = B(\lambda) \]
Es decir el operador \(W(\lambda)\) transforma así el estimador de máxima verosimilitud del parámetro de regresión en su contraparte regularizada.
Momentos del estimador
De la forma anteriormente encontrada se pueden obtener varios datos que caracterizan a los coeficientes de la regresion ridge:
\[ E[B(\lambda)] = E[(X^{T}X + \lambda I)^{-1}X^{T}Y] =(X^{T}X + \lambda I)^{-1}X^{T}E(Y) \] \[=(X^{T}X + \lambda I)^{-1}X^{T}(X^TX)B = B-\lambda(X^{T}X + \lambda I)B \]
Por lo que se puede concluir que el estimador ridge esta sesgado.
\[Var[B(\lambda)]=Var[W(\lambda)B] = W(\lambda)Var(B)W(\lambda) = \sigma^2W(\lambda)(X^TX)^{-1}W(\lambda)^{-1}\]
De aqui que cuando tomamos \(\lambda>>\):
\[\lim_{\lambda \to \infty}Var[B(\lambda)] =\lim_{\lambda \to \infty} \sigma^2W(\lambda)(X^TX)^{-1}W(\lambda)^{-1} = 0_p \] Continuando con el ejemplo numérico y tomando distintos valores de \(\lambda\) podemos notar la convergencia a cero tal como la demostracion lo indica.
Para \(\lambda = 0\) :
- Podemos descartar el término de regularización por completo y volver a nuestra función de error al cuadrado. Este valor de lambda es cuando no usamos la regularización en absoluto y no estamos tratando de resolver el problema del sobreajuste.
Para \(0<\lambda < \infty\) :
- Aquí es cuando se resuelve el problema del sobreajuste. Sin embargo, el valor debe elegirse con cuidado. Si es demasiado menor, todavía habrá un sobreajuste, y si es demasiado grande, el modelo puede ajustarse por debajo de los datos.
Para \(\lambda>>\) :
- Los pesos del modelo serán todos cero, porque tenemos mucho peso en el término de regularización, es decir, en los cuadrados de los pesos.
regularizacion <- modelo$beta %>%
as.matrix() %>%
t() %>%
as_tibble() %>%
mutate(lambda = modelo$lambda)
regularizacion <- regularizacion %>%
pivot_longer(
cols = !lambda,
names_to = "predictor",
values_to = "coeficientes"
)
regularizacion %>%
ggplot(aes(x = lambda, y = coeficientes, color = predictor)) +
geom_line() +
scale_x_log10(
breaks = trans_breaks("log10", function(x) 10^x),
labels = trans_format("log10", math_format(10^.x))
) +
labs(title = "Coeficientes del modelo en función de la regularización") +
theme_bw() +
theme(legend.position = "none")
Eleccion de \(\lambda\)
La forma más común de encontrar lo mejor es mediante la llamada validación cruzada de exclusión:
- Se elige \(p\) posibles valores para el parámetro de penalización.
- Para \(i = 1,...,n\) se excluye la \(i-ésima\) observación \((y_i,x_i)\).
- Se ajusta el estimador rigde para las \(n-1\) observaciones considerando los \(p\) valores de \(\lambda\).
- Se calcula el MSE
- Se elige el valor de \(\lambda\) que minimiza el MSE.
En otras palabras, igualamos el valor del parámetro de penalización que genera el MSE más bajo en el ejercicio de validación cruzada dejando uno fuera.
set.seed(123)
cv_error <- cv.glmnet(
x = x_train,
y = y_train,
alpha = 0,
nfolds = 10,
type.measure = "mse",
standardize = TRUE
)
plot(cv_error)
Una vez realizado el proceso anterior descrito podemos extraer los coeficientes de la regresión final:
Recordatorio
Los supuestos son los mismos que los utilizados en la regresión múltiple: linealidad, varianza constante (sin valores atípicos) e independencia. Dado que la regresión rigde no proporciona intervalos de confianza, no es necesario asumir la normalidad.
# Coeficientes del modelo
# ==============================================================================
df_coeficientes <- coef(modelo) %>%
as.matrix() %>%
as_tibble(rownames = "predictor") %>%
rename(coeficiente = s0)
df_coeficientes %>%
filter(predictor != "(Intercept)") %>%
ggplot(aes(x = predictor, y = coeficiente)) +
geom_col() +
labs(title = "Coeficientes del modelo Ridge") +
theme_bw() +
theme(axis.text.x = element_text(size = 6, angle = 45))
Observaciones Finales
Ventajas:
Resuelve el problema del sobreajuste ya que una regresión múltiple no reconoce las características menos importantes y las usa todas, lo que lleva a un sobreajuste.
En conjuntos de datos en los que tenemos una cantidad de características (n) mayor que la cantidad de ejemplos de entrenamiento (m) , la regresión de cresta se vuelve crucial, ya que funciona significativamente mejor que el método de suma de cuadrados normal.
Desventajas:
Aunque mejora la precisión de la prueba, utiliza todas las características de entrada en el conjunto de datos , a diferencia de los métodos paso a paso que solo seleccionan algunas características importantes para la regresión
La regresión de cresta reduce los coeficientes a valores muy bajos si la característica no es importante, pero no los hará completamente cero , por lo tanto, todavía usamos la característica en nuestro modelo. La regresión de Lasso supera este inconveniente.