Chapter 27 Linear regression complex cases
27.1 Cars
<- mtcars
df
# Subset of necessary data from mtcars
<- df[,c(1,3:7)]
df.sub
# Linear regression model
<- lm(mpg ~ hp, df)
fit
fitsummary(fit)
# Plot regression models using ggplot
# parameter geom_smooth builds regression model
library(ggplot2)
# auto model
ggplot(df, aes(hp, mpg))+geom_point(size=2)+geom_smooth()
# for linear model: method = 'lm'
ggplot(df, aes(hp, mpg))+geom_point(size=2)+geom_smooth(method = "lm")
# Split data into two groups by am - Transmission (0 = automatic, 1 = manual)
ggplot(df, aes(hp, mpg, col=factor(am)))+
geom_point(size=2)+
geom_smooth(method = "lm")
# Prediction of values using linear regression model
<- data.frame(mpg=df$mpg, fitted=fit$fitted.values)
fitted_values_mpg
fitted_values_mpgView(fitted_values_mpg)
# Lets predict galons of petrol for given horse powers
<- data.frame(hp=c(100, 150, 129, 300))
new_hp
predict(fit, new_hp)
$mpg <- predict(fit, new_hp)
new_hp
new_hp
# Lets make regression model for cylinders as numeric (not factor)
<- lm(mpg ~ cyl, df)
fit
fit
# Regression line for cyclinders
ggplot(df, aes(cyl, mpg))+
geom_point()+
geom_smooth(method="lm")+
theme(axis.text=element_text(size=25),
axis.title=element_text(size=25, face="bold"))
27.2 Linear regression modeling, compair with kNN
library('GGally')
library('lmtest')
library('FNN')
# константы
<- 12345
my.seed <- 0.85
train.percent
# загрузка данных
<- 'https://sites.google.com/a/kiber-guu.ru/msep/mag-econ/salary_data.csv?attredirects=0&d=1'
fileURL
# преобразуем категориальные переменные в факторы
<- read.csv(fileURL, row.names = 1, sep = ';', as.is = T)
wages.ru $male <- as.factor(wages.ru$male)
wages.ru$educ <- as.factor(wages.ru$educ)
wages.ru$forlang <- as.factor(wages.ru$forlang)
wages.ru
# обучающая выборка
set.seed(my.seed)
<- sample(seq_along(wages.ru$salary),
inTrain nrow(wages.ru) * train.percent)
<- wages.ru[inTrain, c(colnames(wages.ru)[-1], colnames(wages.ru)[1])]
df.train <- wages.ru[-inTrain, -1]
df.test
# Variable description
# salary – среднемесячная зарплата после вычета налогов за последние 12 месяцев (рублей);
# male – пол: 1 – мужчина, 0 – женщина;
# educ – уровень образования:
# 1 – 0-6 классов,
# 2 – незаконченное среднее (7-8 классов),
# 3 - незаконченное среднее плюс что-то ещё,
# 4 – законченное среднее,
# 5 – законченное среднее специальное, 6 – законченное высшее образование и выше;
# forlang - иност. язык: 1 – владеет, 0 – нет;
# exper – официальный стаж c 1.01.2002 (лет).
summary(df.train)
<- ggpairs(df.train)
ggp print(ggp, progress = F)
# цвета по фактору male
<- ggpairs(df.train[, c('exper', 'male', 'salary')],
ggp mapping = ggplot2::aes(color = male))
print(ggp, progress = F)
# цвета по фактору educ
<- ggpairs(df.train[, c('exper', 'educ', 'salary')],
ggp mapping = ggplot2::aes(color = educ))
print(ggp, progress = F)
# цвета по фактору forlang
<- ggpairs(df.train[, c('exper', 'forlang', 'salary')],
ggp mapping = ggplot2::aes(color = forlang))
print(ggp, progress = F)
# Linear regression model
.1 <- lm(salary ~ . + exper:educ + exper:forlang + exper:male, data = df.train)
modelsummary(model.1)
## Exclude uninfluencial parameters
# Exclude eper:educ as paramaeters are not important
.2 <- lm(salary ~ . + exper:forlang + exper:male, data = df.train)
modelsummary(model.2)
# Exclude male1:exper
.3 <- lm(salary ~ . + exper:forlang, data = df.train)
modelsummary(model.3)
# forlang1 is less important, has no sence
.4 <- lm(salary ~ male + educ + exper, data = df.train)
modelsummary(model.4)
$educ <- as.numeric(df.train$educ)
df.train$educ <- as.numeric(df.test$educ)
df.test
.6 <- lm(salary ~ ., data = df.train)
modelsummary(model.6)
# Model 6 is week, let's add exper:male interactions
$educ <- as.numeric(df.train$educ)
df.train
.7 <- lm(salary ~ . + exper:male, data = df.train)
modelsummary(model.7)
# Obviously the best decision is not to use interactions for modeling
# Test remainers
# тест Бройша-Пагана
bptest(model.6)
# статистика Дарбина-Уотсона
dwtest(model.6)
# графики остатков
par(mar = c(4.5, 4.5, 2, 1))
par(mfrow = c(1, 3))
plot(model.7, 1)
plot(model.7, 4)
plot(model.7, 5)
### Comparison with kNN-method
par(mfrow = c(1, 1))
# фактические значения y на тестовой выборке
<- wages.ru[-inTrain, 1]
y.fact <- predict(model.6, df.test)
y.model.lm <- sum((y.model.lm - y.fact)^2) / length(y.model.lm)
MSE.lm
# kNN требует на вход только числовые переменные
<- as.data.frame(apply(df.train, 2, as.numeric))
df.train.num <- as.data.frame(apply(df.test, 2, as.numeric))
df.test.num
for (i in 2:50){
<- knn.reg(train = df.train.num[, !(colnames(df.train.num) %in% 'salary')],
model.knn y = df.train.num[, 'salary'],
test = df.test.num, k = i)
<- model.knn$pred
y.model.knn if (i == 2){
<- sum((y.model.knn - y.fact)^2) / length(y.model.knn)
MSE.knn else {
} <- c(MSE.knn,
MSE.knn sum((y.model.knn - y.fact)^2) / length(y.model.knn))
}
}
# график
par(mar = c(4.5, 4.5, 1, 1))
plot(2:50, MSE.knn, type = 'b', col = 'darkgreen',
xlab = 'значение k', ylab = 'MSE на тестовой выборке')
lines(2:50, rep(MSE.lm, 49), lwd = 2, col = grey(0.2), lty = 2)
legend('bottomright', lty = c(1, 2), pch = c(1, NA),
col = c('darkgreen', grey(0.2)),
legend = c('k ближайших соседа', 'регрессия (все факторы)'),
lwd = rep(2, 2))
Source
Course ‘Math modeling’ practical work, State University of Management, Moscow
27.3 More complex example
# Linear regression modeling, compair with kNN
# Source: Course 'Math modeling' practical work, State University of Management, Moscow
# link: https://sites.google.com/a/kiber-guu.ru/r-practice/home
library('GGally')
library('lmtest')
library('FNN')
# variable
<- 12345
my.seed <- 0.85
train.percent
<- read.csv("~/Projects/data_analysis/DATA/salary.csv", row.names = 1, sep = '\t',
wages.ru as.is = T) # transform data into factors when possible
$male <- as.factor(wages.ru$male)
wages.ru$educ <- as.factor(wages.ru$educ)
wages.ru$forlang <- as.factor(wages.ru$forlang)
wages.ru
# test data
set.seed(my.seed)
<- sample(seq_along(wages.ru$salary),
inTrain nrow(wages.ru) * train.percent)
<- wages.ru[inTrain, c(colnames(wages.ru)[-1], colnames(wages.ru)[1])]
df.train <- wages.ru[-inTrain, -1]
df.test
summary(df.train)
<- ggpairs(df.train)
ggp print(ggp, progress = F)
# цвета по фактору male
<- ggpairs(df.train[, c('exper', 'male', 'salary')],
ggp mapping = ggplot2::aes(color = male))
print(ggp, progress = F)
# цвета по фактору educ
<- ggpairs(df.train[, c('exper', 'educ', 'salary')],
ggp mapping = ggplot2::aes(color = educ))
print(ggp, progress = F)
# цвета по фактору forlang
<- ggpairs(df.train[, c('exper', 'forlang', 'salary')],
ggp mapping = ggplot2::aes(color = forlang))
print(ggp, progress = F)
# Linear regression model
.1 <- lm(salary ~ . + exper:educ + exper:forlang + exper:male, data = df.train)
modelsummary(model.1)
## Exclude uninfluencial parameters
# Exclude eper:educ as paramaeters are not important
.2 <- lm(salary ~ . + exper:forlang + exper:male, data = df.train)
modelsummary(model.2)
# Exclude male1:exper
.3 <- lm(salary ~ . + exper:forlang, data = df.train)
modelsummary(model.3)
# forlang1 is less important, has no sence
.4 <- lm(salary ~ male + educ + exper, data = df.train)
modelsummary(model.4)
$educ <- as.numeric(df.train$educ)
df.train$educ <- as.numeric(df.test$educ)
df.test
.6 <- lm(salary ~ ., data = df.train)
modelsummary(model.6)
# Model 6 is week, let's add exper:male interactions
$educ <- as.numeric(df.train$educ)
df.train
.7 <- lm(salary ~ . + exper:male, data = df.train)
modelsummary(model.7)
# Obviously the best decision is not to use interactions for modeling
# Test remainers
# тест Бройша-Пагана
bptest(model.6)
# статистика Дарбина-Уотсона
dwtest(model.6)
# графики остатков
par(mar = c(4.5, 4.5, 2, 1))
par(mfrow = c(1, 3))
plot(model.7, 1)
plot(model.7, 4)
plot(model.7, 5)
### Comparison with kNN-method
par(mfrow = c(1, 1))
# фактические значения y на тестовой выборке
<- wages.ru[-inTrain, 1]
y.fact <- predict(model.6, df.test)
y.model.lm <- sum((y.model.lm - y.fact)^2) / length(y.model.lm)
MSE.lm
# kNN требует на вход только числовые переменные
<- as.data.frame(apply(df.train, 2, as.numeric))
df.train.num <- as.data.frame(apply(df.test, 2, as.numeric))
df.test.num
for (i in 2:50){
<- knn.reg(train = df.train.num[, !(colnames(df.train.num) %in% 'salary')],
model.knn y = df.train.num[, 'salary'],
test = df.test.num, k = i)
<- model.knn$pred
y.model.knn if (i == 2){
<- sum((y.model.knn - y.fact)^2) / length(y.model.knn)
MSE.knn else {
} <- c(MSE.knn,
MSE.knn sum((y.model.knn - y.fact)^2) / length(y.model.knn))
}
}
# график
par(mar = c(4.5, 4.5, 1, 1))
plot(2:50, MSE.knn, type = 'b', col = 'darkgreen',
xlab = 'значение k', ylab = 'MSE на тестовой выборке')
lines(2:50, rep(MSE.lm, 49), lwd = 2, col = grey(0.2), lty = 2)
legend('bottomright', lty = c(1, 2), pch = c(1, NA),
col = c('darkgreen', grey(0.2)),
legend = c('k ближайших соседа', 'регрессия (все факторы)'),
lwd = rep(2, 2))
27.4 NEXT part
# Linear regression (part 1)
# Linear Regression (part 2)
<- read.table("~/DataAnalysis/R_data_analysis/DATA/Diamond.dat", header=T)
x
# Check 3 models
# 1. ax+b
# 2. ax^2+b
# 3. ax^2+bx+c
.1 <- lm(x[,2]~x[,1])
ressummary(res.1)
<- x[,1]*x[,1]
ves2 .2 <- lm(x[,2]~ves2)
ressummary(res.2)
.3 <- lm(x[,2]~x[,1]+ves2)
ressummary(res.3)
# weight of diamands ~ weight^2
plot(x[,1]~ves2)
# Conclusion: from 3 models, the 2d is optimal
### Временные ряды
.01 <- read.table("~/DataAnalysis/R_data_analysis/DATA/series_g.csv", header=T, sep=";")
ser.g
# visual inspection
.01
ser.gdim(ser.g.01)
names(ser.g.01)
# plot vaiable
plot(ser.g.01$series_g, type="l")
<- log(ser.g.01$series_g)
log.ser.g
plot(log.ser.g, type="l")
<- 1:(144+12)
time.
.01 <- rep(c(1,0,0,0,0,0,0,0,0,0,0,0), 12+1)
month.02 <- rep(c(0,1,0,0,0,0,0,0,0,0,0,0), 12+1)
month.03 <- rep(c(0,0,1,0,0,0,0,0,0,0,0,0), 12+1)
month.04 <- rep(c(0,0,0,1,0,0,0,0,0,0,0,0), 12+1)
month.05 <- rep(c(0,0,0,0,1,0,0,0,0,0,0,0), 12+1)
month.06 <- rep(c(0,0,0,0,0,1,0,0,0,0,0,0), 12+1)
month.07 <- rep(c(0,0,0,0,0,0,1,0,0,0,0,0), 12+1)
month.08 <- rep(c(0,0,0,0,0,0,0,1,0,0,0,0), 12+1)
month.09 <- rep(c(0,0,0,0,0,0,0,0,1,0,0,0), 12+1)
month.10 <- rep(c(0,0,0,0,0,0,0,0,0,1,0,0), 12+1)
month.11 <- rep(c(0,0,0,0,0,0,0,0,0,0,1,0), 12+1)
month.12 <- rep(c(0,0,0,0,0,0,0,0,0,0,0,1), 12+1)
month
145:(144+12)] <-NA
log.ser.g[
.02 <- data.frame(log.ser.g, time., month.01, month.02, month.03,
ser.g.04, month.05, month.06, month.07, month.08, month.09, month.10,
month.11, month.12)
month
.02
ser.g
.01 <- lm(log.ser.g ~ time. + month.02 + month.03 + month.04 +
res.05 + month.06 + month.07 + month.08 + month.09 + month.10 +
month.11 + month.12, ser.g.02)
month
summary(res.01)
.01$fitted.values
res
plot(ser.g.02$log.ser.g, type="l")
lines(res.01$fitted.values)
plot(ser.g.02$log.ser.g, type="l", col="green")
lines(res.01$fitted.values, col="red")
= predict.lm(res.01, ser.g.02)
x.lg
x.lg
plot(x.lg, type="l", col="red")
lines(ser.g.02$log.ser.g, col="green")
<- exp(x.lg)
y
plot(y, type="l", col="red")
lines(ser.g.01$series_g, col="green")
27.5 NEXT Part
# провести отбор оптимального подмножества переменных;
# отобрать предикторы методами пошагового включения и исключения;
# как построить ридж- и лассо-регрессию;
# как использовать снижение размерности: PCR и PLS;
# как применять эти методы в сочетании с кросс-валидацией.
library('ISLR') # набор данных Hitters
library('leaps') # функция regsubset() -- отбор оптимального
# подмножества переменных
library('glmnet') # функция glmnet() -- лассо
library('pls') # регрессия на главные компоненты -- pcr()
# и частный МНК -- plsr()
<- 1
my.seed
?Hittersfix(Hitters)
names(Hitters)
dim(Hitters)
sum(is.na(Hitters$Salary))
<- na.omit(Hitters)
Hitters dim(Hitters)
sum(is.na(Hitters$Salary))
## Отбор оптимального подмножества
# подгоняем модели с сочетаниями предикторов до 8 включительно
<- regsubsets(Salary ~ ., Hitters)
regfit.full summary(regfit.full)
# подгоняем модели с сочетаниями предикторов до 19 (максимум в данных)
<- regsubsets(Salary ~ ., Hitters, nvmax = 19)
regfit.full <- summary(regfit.full)
reg.summary
reg.summary
# структура отчёта по модели (ищем характеристики качества)
names(reg.summary)
# R^2 и скорректированный R^2
round(reg.summary$rsq, 3)
# на графике
plot(1:19, reg.summary$rsq, type = 'b',
xlab = 'Количество предикторов', ylab = 'R-квадрат')
# сода же добавим скорректированный R-квадрат
points(1:19, reg.summary$adjr2, col = 'red')
# модель с максимальным скорректированным R-квадратом
which.max(reg.summary$adjr2)
### 11
points(which.max(reg.summary$adjr2),
$adjr2[which.max(reg.summary$adjr2)],
reg.summarycol = 'red', cex = 2, pch = 20)
legend('bottomright', legend = c('R^2', 'R^2_adg'),
col = c('black', 'red'), lty = c(1, NA),
pch = c(1, 1))
# C_p
$cp
reg.summary
# число предикторов у оптимального значения критерия
which.min(reg.summary$cp)
### 10
# график
plot(reg.summary$cp, xlab = 'Число предикторов',
ylab = 'C_p', type = 'b')
points(which.min(reg.summary$cp),
$cp[which.min(reg.summary$cp)],
reg.summarycol = 'red', cex = 2, pch = 20)
# BIC
$bic
reg.summary
# число предикторов у оптимального значения критерия
which.min(reg.summary$bic)
### 6
# график
plot(reg.summary$bic, xlab = 'Число предикторов',
ylab = 'BIC', type = 'b')
points(which.min(reg.summary$bic),
$bic[which.min(reg.summary$bic)],
reg.summarycol = 'red', cex = 2, pch = 20)
# метод plot для визуализации результатов
?plot.regsubsetsplot(regfit.full, scale = 'r2')
plot(regfit.full, scale = 'adjr2')
plot(regfit.full, scale = 'Cp')
plot(regfit.full, scale = 'bic')
# коэффициенты модели с наименьшим BIC
round(coef(regfit.full, 6), 3)
## Отбор путём пошагового включения и исключения переменных
# Пошаговое включение
<- regsubsets(Salary ~ ., data = Hitters, nvmax = 19, method = 'forward')
regfit.fwd summary(regfit.fwd)
<- regsubsets(Salary ~ ., data = Hitters,
regfit.bwd nvmax = 19, method = 'backward')
summary(regfit.bwd)
round(coef(regfit.full, 7), 3)
round(coef(regfit.fwd, 7), 3)
round(coef(regfit.bwd, 7), 3)
### Нахождение оптимальной модели при помощи методов проверочной выборки и перекрёстной проверки
set.seed(my.seed)
<- sample(c(T, F), nrow(Hitters), rep = T)
train <- !train
test
# обучаем модели
<- regsubsets(Salary ~ ., data = Hitters[train, ],
regfit.best nvmax = 19)
# матрица объясняющих переменных модели для тестовой выборки
<- model.matrix(Salary ~ ., data = Hitters[test, ])
test.mat
# вектор ошибок
<- rep(NA, 19)
val.errors # цикл по количеству предикторов
for (i in 1:19){
<- coef(regfit.best, id = i)
coefi <- test.mat[, names(coefi)] %*% coefi
pred # записываем значение MSE на тестовой выборке в вектор
<- mean((Hitters$Salary[test] - pred)^2)
val.errors[i]
}round(val.errors, 0)
# находим число предикторов у оптимальной модели
which.min(val.errors)
### 10
# коэффициенты оптимальной модели
round(coef(regfit.best, 10), 3)
# функция для прогноза для функции regsubset()
<- function(object, newdata, id, ...){
predict.regsubsets <- as.formula(object$call[[2]])
form <- model.matrix(form, newdata)
mat <- coef(object, id = id)
coefi <- names(coefi)
xvars %*% coefi
mat[, xvars]
}
# набор с оптимальным количеством переменных на полном наборе данных
<- regsubsets(Salary ~ ., data = Hitters,
regfit.best nvmax = 19)
round(coef(regfit.best, 10), 3)
# k-кратная кросс-валидация
# отбираем 10 блоков наблюдений
<- 10
k set.seed(my.seed)
<- sample(1:k, nrow(Hitters), replace = T)
folds
# заготовка под матрицу с ошибками
<- matrix(NA, k, 19, dimnames = list(NULL, paste(1:19)))
cv.errors
# заполняем матрицу в цикле по блокам данных
for (j in 1:k){
<- regsubsets(Salary ~ ., data = Hitters[folds != j, ],
best.fit nvmax = 19)
# теперь цикл по количеству объясняющих переменных
for (i in 1:19){
# модельные значения Salary
<- predict(best.fit, Hitters[folds == j, ], id = i)
pred # вписываем ошибку в матрицу
<- mean((Hitters$Salary[folds == j] - pred)^2)
cv.errors[j, i]
}
}
# усредняем матрицу по каждому столбцу (т.е. по блокам наблюдений),
# чтобы получить оценку MSE для каждой модели с фиксированным
# количеством объясняющих переменных
<- apply(cv.errors, 2, mean)
mean.cv.errors round(mean.cv.errors, 0)
# на графике
plot(mean.cv.errors, type = 'b')
points(which.min(mean.cv.errors), mean.cv.errors[which.min(mean.cv.errors)],
col = 'red', pch = 20, cex = 2)
# перестраиваем модель с 11 объясняющими переменными на всём наборе данных
<- regsubsets(Salary ~ ., data = Hitters, nvmax = 19)
reg.best round(coef(reg.best, 11), 3)
# из-за синтаксиса glmnet() формируем явно матрицу объясняющих...
<- model.matrix(Salary ~ ., Hitters)[, -1]
x
# и вектор значений зависимой переменной
<- Hitters$Salary
y
### Гребневая регрессия
# вектор значений гиперпараметра лямбда
<- 10^seq(10, -2, length = 100)
grid
# подгоняем серию моделей ридж-регрессии
<- glmnet(x, y, alpha = 0, lambda = grid)
ridge.mod
# размерность матрицы коэффициентов моделей
dim(coef(ridge.mod))
## [1] 20 100
# значение лямбда под номером 50
round(ridge.mod$lambda[50], 0)
## [1] 11498
# коэффициенты соответствующей модели
round(coef(ridge.mod)[, 50], 3)
# норма эль-два
round(sqrt(sum(coef(ridge.mod)[-1, 50]^2)), 2)
# всё то же для лямбды под номером 60
# значение лямбда под номером 50
round(ridge.mod$lambda[60], 0)
# коэффициенты соответствующей модели
round(coef(ridge.mod)[, 60], 3)
# норма эль-два
round(sqrt(sum(coef(ridge.mod)[-1, 60]^2)), 1)
# мы можем получить значения коэффициентов для новой лямбды
round(predict(ridge.mod, s = 50, type = 'coefficients')[1:20, ], 3)
## Метод проверочной выборки
set.seed(my.seed)
<- sample(1:nrow(x), nrow(x)/2)
train <- -train
test <- y[test]
y.test
# подгоняем ридж-модели с большей точностью (thresh ниже значения по умолчанию)
<- glmnet(x[train, ], y[train], alpha = 0, lambda = grid,
ridge.mod thresh = 1e-12)
plot(ridge.mod)
# прогнозы для модели с лямбда = 4
<- predict(ridge.mod, s = 4, newx = x[test, ])
ridge.pred round(mean((ridge.pred - y.test)^2), 0)
# сравним с MSE для нулевой модели (прогноз = среднее)
round(mean((mean(y[train]) - y.test)^2), 0)
# насколько модель с лямбда = 4 отличается от обычной ПЛР
<- predict(ridge.mod, s = 0, newx = x[test, ], exact = T,
ridge.pred x = x[train, ], y = y[train])
round(mean((ridge.pred - y.test)^2), 0)
# predict с лямбдой (s) = 0 даёт модель ПЛР
lm(y ~ x, subset = train)
round(predict(ridge.mod, s = 0, exact = T, type = 'coefficients',
x = x[train, ], y = y[train])[1:20, ], 3)
## Подбор оптимального значения лямбда с помощью перекрёстной проверки
# k-кратная кросс-валидация
set.seed(my.seed)
# оценка ошибки
<- cv.glmnet(x[train, ], y[train], alpha = 0)
cv.out plot(cv.out)
# значение лямбда, обеспечивающее минимальную ошибку перекрёстной проверки
<- cv.out$lambda.min
bestlam round(bestlam, 0)
## [1] 212
# MSE на тестовой для этого значения лямбды
<- predict(ridge.mod, s = bestlam, newx = x[test, ])
ridge.pred round(mean((ridge.pred - y.test)^2), 0)
## [1] 96016
# наконец, подгоняем модель для оптимальной лямбды,
# найденной по перекрёстной проверке
<- glmnet(x, y, alpha = 0)
out round(predict(out, type = 'coefficients', s = bestlam)[1:20, ], 3)
## Лассо
<- glmnet(x[train, ], y[train], alpha = 1, lambda = grid)
lasso.mod plot(lasso.mod)
# Подбор оптимального значения лямбда с помощью перекрёстной проверки
set.seed(my.seed)
<- cv.glmnet(x[train, ], y[train], alpha = 1)
cv.out plot(cv.out)
<- cv.out$lambda.min
bestlam <- predict(lasso.mod, s = bestlam, newx = x[test, ])
lasso.pred round(mean((lasso.pred - y.test)^2), 0)
# коэффициенты лучшей модели
<- glmnet(x, y, alpha = 1, lambda = grid)
out <- predict(out, type = 'coefficients', s = bestlam)[1:20, ]
lasso.coef round(lasso.coef, 3)
round(lasso.coef[lasso.coef != 0], 3)
### Лабораторная работа 3: регрессия при помощи методов PCR и PLS
### 6.7.1 Регрессия на главные компоненты
# кросс-валидация
set.seed(2) # непонятно почему они сменили зерно; похоже, опечатка
<- pcr(Salary ~ ., data = Hitters, scale = T, validation = 'CV')
pcr.fit summary(pcr.fit)
# график ошибок
validationplot(pcr.fit, val.type = 'MSEP')
# Подбор оптиального M: кросс-валидация на обучающей выборке
set.seed(my.seed)
<- pcr(Salary ~ ., data = Hitters, subset = train, scale = T,
pcr.fit validation = 'CV')
validationplot(pcr.fit, val.type = 'MSEP')
# MSE на тестовой выборке
<- predict(pcr.fit, x[test, ], ncomp = 7)
pcr.pred round(mean((pcr.pred - y.test)^2), 0)
# подгоняем модель на всей выборке для M = 7
# (оптимально по методу перекрёстной проверки)
<- pcr(y ~ x, scale = T, ncomp = 7)
pcr.fit summary(pcr.fit)
# Регрессия по методу частных наименьших квадратов
set.seed(my.seed)
<- plsr(Salary ~ ., data = Hitters, subset = train, scale = T,
pls.fit validation = 'CV')
summary(pls.fit)
# теперь подгоняем модель для найденного оптимального M = 2
# и оцениваем MSE на тестовой
<- predict(pls.fit, x[test, ], ncomp = 2)
pls.pred round(mean((pls.pred - y.test)^2), 0)
## [1] 101417
# подгоняем модель на всей выборке
<- plsr(Salary ~ ., data = Hitters, scale = T, ncomp = 2)
pls.fit summary(pls.fit)