這個課程將跟同學介紹近代機器學習常用的一個技術。
在之前的課程，我們是選擇解釋變數的組合、交互作用等方式，來挑選出好的模型。
這堂課程介紹的作法，是一次放入所有可能的變數，然後透過修改「評分」項目，讓演算法來自動找出應該要留在模型裡的變數。這種方法，專業術語叫做：Regularization。
R的glmnet套件在線性模型上實作了Regularization的相關功能。這個套件的作者群，也是機器學習領域中鼎鼎大名的學者。
請同學先安裝glmnet套件。
check_then_install("glmnet", "2.0.3")接著請載入glmnet套件
library(glmnet)摸索套件的第一步：找尋vignettes。請同學輸入：vignette(package="glmnet")
vignette(package = "glmnet")請同學輸入vignette("glmnet_beta",package="glmnet")打開套件的簡介文件。
vignette("glmnet_beta", package = "glmnet")由於閱讀這份說明文件，需要統計學的專業知識，所以我們就先暫停在這邊。有興趣的同學可以再自行鑽研。
這裡我們先用iris這個資料集來介紹glmnet與regularization。請同學先輸入：m0<-lm(Sepal.Length~.,iris)這是第一個單元介紹的線性模型，透過這個指令，我們可以取得最小方差直線。我們等等會拿這個模型與glmnet的結果做比較。
m0 <- lm(Sepal.Length ~ ., iris)glmnet的參數並不接受formula，而是接受x(矩陣)與y(向量)。請同學先輸入：X<-model.matrix(Sepal.Length~.,iris)建立一個矩陣好填入glmnet的x參數。model.matrix是lm使用formula時，在背後產生線性代數矩陣的函數。
X <- model.matrix(Sepal.Length ~ ., iris)接著，請同學輸入：y<-iris$Sepal.Length建立一個向量好填入glmnet的y參數。
y <- iris$Sepal.Length我們先輸入：m<-glmnet(X[,-1],y,lambda=0)這裡使用X[,-1]是因為glmnet會自動加入Intercept，所以我們要把Intercept從X中移除（Intercept是第一欄）。關於lambda，我們先賣個關子。
m <- glmnet(X[,-1], y, lambda = 0)請同學輸入：cbind(coef(m,s=0),coef(m0))比較一下glmnet學出來的模型與lm學出來的模型。coef的s=0參數等等會補充說明。
cbind(coef(m, s = 0), coef(m0))我們應該看到非常接近的數字。這是因為兩者要解的目標函數是一樣的，但是glmnet是用數值解，而lm是用理論直接計算，所以lm是比較準確的。
接著，我們慢慢加大lambda參數。請同學出入：m<-glmnet(X[,-1],y,lambda=seq(1,0,by=-0.1))
m <- glmnet(X[,-1], y, lambda = seq(1, 0, by = -0.1))請同學看看m的輸出結果。
m事實上，glmnet剛剛在背後一口氣算了時一個模型，分別對應到11個lambda（即seq(1.0,0.0,by=-0.1)）
我們可以運用coef(m,s=0.5)來查詢當lambda=0.5時，學出來的參數。請同學試試看
coef(m, s = 0.5)同學應該可以注意到，除了Petal.Length之外，其他的參數學出來都是“.”。這就是glmnet透過regularization的技術，去挑選變數。挑選的結果，在lambda=0.5之下，只有挑出Petal.Length而已，其他的變數都是0。
接著請同學輸入：coef(m,s=0.1)
coef(m, s = 0.1)這次同學可以看到，當lambda是0.1時，glmnet挑出的變數是Sepal.Width與Petal.Length。
所以lambda到底是什麼呢？機器學習的原理，就是在訂定評分機制後，找到能夠最優化評分的參數。而regularization，就是在原本的評分之外，加上對參數數值的懲罰。具體作法是當參數的值越大，評分就會越低。而lambda就是控制懲罰強度的參數。
舉例來說，當lambda=0的時候，代表線性模型的參數值的懲罰加權為0（也就是不懲罰），所以學習出來的模型就等價於原本的線性模型。但是當lambda=1時，參數值得懲罰讓機器學習最後把大部分的參數刪除，只留下Petal.Length。這就是regularization。
glmnet所實作的regularization有兩種。一種是運用模型參數的絕對值的和做懲罰，這種又稱為L1-Regularization。統計上的別名是Lasso。另一種是以模型參數的平方的和作為懲罰，這種稱為L2-Regularization，統計上又稱為RidgeRegression。
glmnet中的參數alpha控制了L1和L2Regularization的權重。當alpha=0，就只剩下L2-Regularization；當alpha=1，就只剩下L1-Regularization。當alpha介於0、1之間，兩者並存時，統計上稱作elasticnet。
回到選變數的問題。當我們使用預設的alpha=1時，隨著lambda的加大或縮小，模型參數變成0的順序，常常被用來當成評估變數重要性。舉例來說，在m之中，我們可以看到Petal.Length可能是最重要的參數，接著才是Sepal.Width。
我們回過頭解釋m的顯示結果。請同學輸入m
mm會顯示一個有三個欄位的矩陣。第一欄位（Df，DegreeofFreedom的簡稱）就代表有多少個不為0的模型參數；第二欄位（%Dev）就代表這個模型在trainingdataset上解釋了多少比率的變化；第三個欄位（Lambda）則代表lambda的值，也就是Regularization的強度（懲罰模型參數的值的強度）。
要注意的是，類似的選模只有在alpha>0的時候才會有效果。當alpha=0時，參數通常都不會是0。請同學試試看：glmnet(X[,-1],y,alpha=0,lambda=seq(1,0,by=-0.1))
glmnet(X[,-1], y, alpha = 0, lambda = seq(1, 0, by = -0.1))透過Df，同學可以發現所有的變數都不為0。
類似glm，glmnet除了傳統的線性模型之外，也將類似的概念實作於廣義線性模型上。glmnet會偵測y的型態，或是指定family後，就會以對應的廣義線性模型做處理。我們直接在最後的關卡練習glmnet提供的logisticregression。
實務上，glmnet這套作法雖然解決選模的問題，但是帶來挑選lambda的問題。像lambda這類參數，很難透過傳統的統計方法找到好的值。所以一般都會用trainingdatasetv.s.testingdataset的概念來解決。也就是說，看看不同lambda選出來的模型，在testingdataset上的表現優劣，來挑選好的lambda。
最後我們請同學用上一個單元的使用過得資料集來作練習。
check_then_install("mlbench", "2.1.1")
library(mlbench)
# 方便起見，同學可以使用這個函數計算 Logarithmic Loss
logloss <- function(y, p, tol = 1e-4) {
# tol 的用途是避免對0取log所導致的數值問題
p[p < tol] <- tol
p[p > 1 - tol] <- 1 - tol
-sum(y * log(p) + (1 - y) * log(1-p))
}
data(Ionosphere)
test.i <- c(4L, 6L, 9L, 13L, 14L, 22L, 31L, 33L, 50L, 52L, 61L, 63L, 68L,
79L, 91L, 99L, 119L, 135L, 154L, 155L, 160L, 162L, 166L, 194L,
200L, 219L, 233L, 236L, 237L, 242L, 244L, 248L, 250L, 257L, 261L,
276L, 278L, 283L, 292L, 310L, 312L, 315L, 319L, 323L, 325L, 327L,
335L, 337L, 338L, 344L)
df.test <- Ionosphere[test.i,-2] # remove V2
train.i <- setdiff(seq_len(nrow(Ionosphere)), test.i)
df.train <- Ionosphere[train.i,-2]
# 上一個單元，我們請同學找到df.test上的logloss小於26的模型
# 我們想要讓同學透過glmnet找到更好的模型
# X.train 是等等要放到glmnet的x參數的值
X.train <- model.matrix(Class ~ ., df.train)[,-1]
# y.train 是等等要放到glmnet的y參數的值
y.train <- df.train$Class
# 請同學透過找到適當的lambda、alplha組合
# 利用glmnet在df.train上學出一個模型，
# 它在df.test上logloss小於8
answer_03 <- local({
NULL
# 你的最終物件應該要由以下的程式碼產生：
# glmnet(x = X.train, y = y.train, lambda = ?, alpha = ?, family = "binomial")
glmnet(x = X.train, y = y.train, lambda = 0.00514326, alpha = 1, family = "binomial")
})
stopifnot(class(answer_03) == c("lognet", "glmnet"))
stopifnot(length(answer_03$lambda) == 1) # 你的答案應該要指定lambda
if (interactive()) { # 讓自動測試略過以下程式碼
stopifnot(logloss(
df.test$Class == "good",
predict(answer_03, model.matrix(Class ~ ., df.test)[,-1], type = "response")) <
8
)
}
# 完成以後請存檔並回到console輸入`submit()`