Task 3. Chaper 6, ex 11

library(MASS)
library(leaps)
library(glmnet)
library(pls)
data("Boston")

11a

Best Subset Selection

set.seed(1)
predict.regsubsets = function(object, newdata, id, ...) {
    form = as.formula(object$call[[2]])
    mat = model.matrix(form, newdata)
    coefi = coef(object, id = id)
    mat[, names(coefi)] %*% coefi
}

k = 10
p = ncol(Boston) - 1
folds = sample(1:k, nrow(Boston), replace = TRUE)
cv.errors = matrix(NA, k, p)
for (i in 1:k) {
    best.fit = regsubsets(crim ~ ., data = Boston[folds != i, ], nvmax = p)
    for (j in 1:p) {
        pred = predict(best.fit, Boston[folds == i, ], id = j)
        cv.errors[i, j] = mean((Boston$crim[folds == i] - pred)^2)
    }
}

err.best.seq = apply(cv.errors, 2, mean)
plot(err.best.seq ,type = "b")

err.best <- err.best.seq[which.min(err.best.seq)]

Можем видеть, что метод кросс-валидации выбрал модель с 12 переменными. Оценка test MSE - 41.0345657.

reg.best <- regsubsets(crim~., data=Boston, nvmax=ncol(Boston)-1)

reg.summary <- summary(reg.best)
par ( mfrow =c(3,1) )
plot(reg.summary$adjr2 ,xlab =" Number of Variables ",ylab =" Adjusted RSq ", type ="l")
max.agjr2 <- which.max(reg.summary$adjr2)
points (max.agjr2, reg.summary$adjr2[max.agjr2] , col =" red ", cex =2, pch =20)

plot(reg.summary$cp ,xlab =" Number of Variables ", ylab =" Cp", type="l")
min.cp <- which.min(reg.summary$cp )
points (min.cp, reg.summary$cp[min.cp] , col =" red ", cex =2, pch =20)

min.bic <- which.min(reg.summary$bic)
plot(reg.summary$bic ,xlab =" Number of Variables ", ylab =" BIC", type="l")
points (min.bic, reg.summary$bic[min.bic] , col =" red ", cex =2, pch =20)

На графиках представлены завистимости \(c_p\), BIC и Adjusted RSq от количества переменных в модели. Красными точками отмечено оптимальное количество переменных по каждому из критериев. Как можно видеть, ни один не совпал с количеством переменных по кросс-валидации.

coef(reg.best, 12)
##   (Intercept)            zn         indus          chas           nox 
##  16.985713928   0.044673247  -0.063848469  -0.744367726 -10.202169211 
##            rm           dis           rad           tax       ptratio 
##   0.439588002  -0.993556631   0.587660185  -0.003767546  -0.269948860 
##         black         lstat          medv 
##  -0.007518904   0.128120290  -0.198877768

Коэффициенты регрессии по методу Best Subset Selection.

Lasso Regression

x <- model.matrix(crim ~ ., Boston)[, -1]
y <- Boston$crim
set.seed(1)
train = sample(nrow(x), nrow(x)/2)
test = (-train)

fit.lasso <- cv.glmnet(x[train,], y[train], alpha=1)
plot(fit.lasso)

lambda <- fit.lasso$lambda.min
pred.lasso <- predict(fit.lasso, s=lambda, newx=x[test,])
err.lasso <- mean((pred.lasso - y[test])^2)
predict(fit.lasso, s=lambda, type="coefficients")
## 14 x 1 sparse Matrix of class "dgCMatrix"
##                        1
## (Intercept)  7.204613974
## zn           0.037084733
## indus       -0.034037718
## chas        -0.518884600
## nox         -9.663821285
## rm           1.187198535
## age          .          
## dis         -0.971658064
## rad          0.517060059
## tax          .          
## ptratio     -0.220503429
## black       -0.002432899
## lstat        0.176698759
## medv        -0.193172676

По методу Лассо было выбрано 11 переменных в модель. test MSE = 38.3766396.

Ridge Regression

set.seed(1)
fit.ridge <- cv.glmnet(x[train,], y[train], alpha=0)
plot(fit.ridge)

lambda <- fit.ridge$lambda.min
pred.ridge <- predict(fit.ridge, s=lambda, newx=x[test,])
err.ridge <- mean((pred.ridge - y[test])^2)
predict(fit.ridge, s=lambda, type="coefficients")
## 14 x 1 sparse Matrix of class "dgCMatrix"
##                        1
## (Intercept)  2.861544809
## zn           0.034393951
## indus       -0.058853509
## chas        -0.785807275
## nox         -6.750896984
## rm           1.108697728
## age          0.005922744
## dis         -0.809673182
## rad          0.398888728
## tax          0.004266196
## ptratio     -0.131627213
## black       -0.004306092
## lstat        0.190997048
## medv        -0.157552605

Здесь test MSE = 38.3671924

PCR

set.seed(2)
pcr.fit = pcr(crim ~ ., data = Boston, subset = train, scale = TRUE, validation = "CV")
summary(pcr.fit)
## Data:    X dimension: 253 13 
##  Y dimension: 253 1
## Fit method: svdpc
## Number of components considered: 13
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV           8.892    7.426    7.419    7.145    7.135    7.126    7.154
## adjCV        8.892    7.424    7.418    7.139    7.122    7.121    7.146
##        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
## CV       7.157    7.113    7.032     7.050     7.087     7.106     6.960
## adjCV    7.146    7.114    7.014     7.031     7.068     7.083     6.938
## 
## TRAINING: % variance explained
##       1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
## X       49.04    60.72    69.75    76.49    83.02    88.40    91.73
## crim    30.39    30.93    36.63    37.31    37.35    37.98    38.85
##       8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
## X       93.77    95.73     97.36     98.62     99.57    100.00
## crim    39.94    41.89     42.73     42.73     43.55     45.48
validationplot(pcr.fit, val.type ="MSEP")

По PCR было выбрано 13 переменных, то есть фактически размерность не уменьшилась. test MSE = 44.502241

11b

Сравним test MSE всех четырех моделей, для того чтобы выбрать наилучшую.

Best Subset Selection: 41.0345657

Lasso Regression: 38.3766396

Ridge Regression: 38.3671924

PCR: 44.502241

Наименьшая ошибка у Ridge Regression, однко для Lasso Regression ошибка совсем чуть-чуть больше.

11c

Так как нет большой разницы между MSE для ridge и lasso, я предпочту выбрать модель lasso, так как интерпретация коэффициентов регресии проще из-за того, что некоторые коэффициенты (в нашем случае их 3) в точности равны нулю.