Ридж-регрессия — это метод, который мы можем использовать для подбора модели регрессии, когда в данных присутствует мультиколлинеарность .
В двух словах, регрессия методом наименьших квадратов пытается найти оценки коэффициентов, которые минимизируют сумму квадратов остатков (RSS):
RSS = Σ(y i – ŷ i )2
куда:
- Σ : греческий символ, означающий сумму
- y i : Фактическое значение отклика для i -го наблюдения
- ŷ i : прогнозируемое значение отклика на основе модели множественной линейной регрессии.
И наоборот, гребневая регрессия стремится минимизировать следующее:
RSS + λΣβ j 2
где j находится в диапазоне от 1 до p переменных-предикторов и λ ≥ 0.
Этот второй член уравнения известен как штраф за усадку.В гребневой регрессии мы выбираем значение λ, которое дает наименьшую возможную тестовую MSE (среднеквадратическую ошибку).
В этом руководстве представлен пошаговый пример выполнения гребневой регрессии в R.
Шаг 1: Загрузите данные
В этом примере мы будем использовать встроенный набор данных R под названием mtcars.Мы будем использовать hp в качестве переменной ответа и следующие переменные в качестве предикторов:
- миль на галлон
- вес
- дрянь
- qсек
Для выполнения гребневой регрессии мы будем использовать функции из пакета glmnet.Этот пакет требует, чтобы переменная ответа была вектором, а набор переменных-предикторов имел класс data.matrix .
Следующий код показывает, как определить наши данные:
#define response variable
y <- mtcars$hp
#define matrix of predictor variables
x <- data.matrix(mtcars[, c('mpg', 'wt', 'drat', 'qsec')])
Шаг 2: Подберите модель регрессии хребта
Далее мы воспользуемся функцией glmnet() , чтобы подогнать модель регрессии хребта, и укажем alpha=0 .
Обратите внимание, что установка альфы, равной 1, эквивалентна использованию регрессии Лассо, а установка альфы на некоторое значение от 0 до 1 эквивалентна использованию эластичной сети.
Также обратите внимание, что гребневая регрессия требует стандартизации данных таким образом, чтобы каждая предикторная переменная имела среднее значение 0 и стандартное отклонение 1.
К счастью , glmnet() автоматически выполняет эту стандартизацию за вас. Если вы уже стандартизировали переменные, вы можете указать standardize=False .
library (glmnet)
#fit ridge regression model
model <- glmnet(x, y, alpha = 0 )
#view summary of model
summary(model)
Length Class Mode
a0 100 -none- numeric
beta 400 dgCMatrix S4
df 100 -none- numeric
dim 2 -none- numeric
lambda 100 -none- numeric
dev.ratio 100 -none- numeric
nulldev 1 -none- numeric
npasses 1 -none- numeric
jerr 1 -none- numeric
offset 1 -none- logical
call 4 -none- call
nobs 1 -none- numeric
Шаг 3: Выберите оптимальное значение для лямбда
Далее мы определим значение лямбда, которое дает наименьшую среднеквадратичную ошибку теста (MSE), используя k-кратную перекрестную проверку .
К счастью, в glmnet есть функция cv.glmnet() , которая автоматически выполняет k-кратную перекрестную проверку, используя k = 10 кратностей.
#perform k-fold cross-validation to find optimal lambda value
cv_model <- cv. glmnet (x, y, alpha = 0 )
#find optimal lambda value that minimizes test MSE
best_lambda <- cv_model$ lambda.min
best_lambda
[1] 10.04567
#produce plot of test MSE by lambda value
plot(cv_model)
Значение лямбда, которое минимизирует тест MSE, оказывается равным 10,04567 .
Шаг 4: Анализ окончательной модели
Наконец, мы можем проанализировать окончательную модель, полученную с помощью оптимального значения лямбда.
Мы можем использовать следующий код, чтобы получить оценки коэффициентов для этой модели:
#find coefficients of best model
best_model <- glmnet(x, y, alpha = 0 , lambda = best_lambda)
coef(best_model)
5 x 1 sparse Matrix of class "dgCMatrix"
s0
(Intercept) 475.242646
mpg -3.299732
wt 19.431238
drat -1.222429
qsec -17.949721
Мы также можем создать график трассировки, чтобы визуализировать, как оценки коэффициентов изменились в результате увеличения лямбда:
#produce Ridge trace plot
plot(model, xvar = " lambda ")
Наконец, мы можем рассчитатьR-квадрат модели на данных обучения:
#use fitted best model to make predictions
y_predicted <- predict (model, s = best_lambda, newx = x)
#find SST and SSE
sst <- sum ((y - mean (y))^2)
sse <- sum ((y_predicted - y)^2)
#find R-Squared
rsq <- 1 - sse/sst
rsq
[1] 0.7999513
R-квадрат оказывается равным 0,7999513.То есть лучшая модель смогла объяснить 79,99% вариации значений отклика обучающих данных.
Полный код R, использованный в этом примере, вы можете найти здесь .