這個課程將介紹R在廣義線性模型上的實作。
傳統的線性迴歸分析中，目標變數通常是沒有範圍限制的數值型變數。當目標變數是非負整數（例：來客數）、或是類別型變數（例：有無點擊）時，我們可以運用廣義線性模型來做學習。
請同學先輸入：?glm打開廣義非線性模型的說明頁面。
?glmglm的參數和lm很接近，主要都是使用formula與data。formula接受formula物件，而data接受data.frame物件。比較不同的是family參數。
請問glm的family參數預設是？
gaussian
當family是gaussian時，glm的行為和lm一樣。同學可以輸入：glm(dist~speed,data=cars)看看。
glm(dist ~ speed, data = cars)但是當我們要作分類問題時，family就要做更改。
請同學輸入skip()檢查mlbench套件
check_then_install("mlbench", "2.1.1")請同學載入mlbench套件
library(mlbench)請同學輸入：data(Ionosphere)載入一個物理實驗的資料集
data(Ionosphere)上一個單元，我們使用AIC來比較模型的好壞。這個單元，我們把資料分成training和testing資料集合，並用模型在testingdataset上的表現做比較。
同學可以用sample函數來自1:351中抽取testingdataset的列序號。舉例來說，test.i中就是一組用sample產生的列序號。train.i則是運用setdiff產生的列序號。
請同學輸入：df.train<-Ionosphere[train.i,]建立trainingdataset的data.frame
df.train <- Ionosphere[train.i,]請同學輸入：df.test<-Ionosphere[test.i,]建立testingdataset的data.frame
df.test <- Ionosphere[test.i,]我們先觀察資料集。請同學輸入：summary(df.train)
summary(df.train)同學有沒有注意到，有一個欄位的值，沒有變化。請同學以字串形式輸入該欄位的名稱。
"V2"
因此，我們模型中就不考慮“V2”。我們只考慮“V1”,“V3”~“V34”當成解釋變數。請同學輸入：x.name<-setdiff(paste0("V",1:34),"V2")來產生一個代表所有要放到模型中的解釋變數
x.name <- setdiff(paste0("V", 1:34), "V2")接著，我們輸入：f<-reformulate(x.name,"Class")來建立我們需要的formula物件。
f <- reformulate(x.name, "Class")接著，請同學輸入：m1<-glm(formula=f,family="binomial",data=df.train)建立一個logisticregression的物件。
m1 <- glm(formula = f, family = "binomial", data = df.train)請有興趣的同學記住：如果要作兩個類別的分類問題，使用glm時，family=“binomial”就等價於logisticregression。
螢幕上出現的：“glm.fit:fittedprobabilitiesnumerically0or1occurred”是因為logisticregression在某些類別組合中目標變數全部是相同類別時，會有數值問題。這邊我們先略過，等到下一個單元再給一個解決辦法。
運用p1.train<-predict(m1,df.train,type="response")就可以算出模型m1在df.train上的預測結果。
p1.train <- predict(m1, df.train, type = "response")我們真正使用的函數，其實是predict.glm，所以有興趣的同學可以打開predict.glm去深入理解這個函數的威能。
請同學輸入：plot(Class~p1.train,data=df.train)看一看預測結果與實際結果在trainingdataset上的比較。
plot(Class ~ p1.train, data = df.train)我們可以看到，隨著p1.train的值越來越大，good所佔有的比率越來越高。而當p1.train<0.1時，所有的資料都是bad；當p1.train超過某個超過0.7的值時，所有的資料都是good。我們判斷，m1是有學到trainingdataset上的趨勢。
接著我們看看testingdataset上的表現。請同學輸入：p1<-predict(m1,df.test,type="response")
p1 <- predict(m1, df.test, type = "response")請同學輸入：plot(Class~p1,df.test)看一看預測結果在df.test上的表現。
plot(Class ~ p1, df.test)我們可以看到，當預測機率為0~0.2時，結果全部都是“bad”。但是當預測結果超過0.8時，“good”的比率大致上是0.1左右。由圖判斷，我們學出來的模型在testingdataset上是有道理，真的有學到判斷bad/good的方法。
運用上一個單元學到的知識，我們可以利用交互作用將模型複雜化。formula的處理，我是將http://stackoverflow.com/a/29691154/1182304中寫的函數，匯入到課程環境中給同學使用。他們並不是內建的。請同學輸入：f2<-interact_rhs(f)
f2 <- interact_rhs(f)我們看看f2的結果是不是列出解釋變數的所有二次交互作用
f2接著，請輸入：m2<-glm(formula=f2,family="binomial",data=df.train)用更複雜的模型來作學習。
m2 <- glm(formula = f2, family = "binomial", data = df.train)事實上這個模型已經過度複雜，參數超過了資料的個數，所以會完美的詮釋trainingdataset。
請同學輸入：p2.train<-predict(m2,df.train,type="response")
p2.train <- predict(m2, df.train, type = "response")接著我們輸入：plot(Class~p2.train,df.train)
plot(Class ~ p2.train, df.train)我們可以看到，如果切點在0.1的話，m2完美的將bad/good給切割開來了。
我們再輸入p2<-predict(m2,df.test,type="response")
p2 <- predict(m2, df.test, type = "response")我們再輸入plot(Class~p2,df.test)以視覺化的方式看看在df.test上的表現
plot(Class ~ p2, df.test)我們可以看到，預測的不同機率值下，bad/good的比率沒有明顯變化。顯示p2的分類結果很不好。
我們也可以用一些數值指標來比較p1與p2的表現。一種常用的指標是LogarithmicLoss：-(ylog(p)+(1-y)log(1-p))這裡的y是以0或1來代表分類結果。p則是預測發生1的機率。當完美預測（y=1時p=1,y=0時p=0），LogarithmicLoss的結果會趨近於0。當結果越不準確，LogarithmicLoss會越大
透過R的向量化運算，我們可以很快速的算出這個指標。請同學先輸入：y<-df.test$Class=="good"代表解答
y <- df.test$Class == "good"我們計算p1的LogarithmicLoss。請同學輸入：-sum(y*log(p1)+(1-y)*log(1-p1))
-sum(y * log(p1) + (1 - y) * log(1-p1))我們計算p2的LogarithmicLoss。請同學輸入：-sum(y*log(p2)+(1-y)*log(1-p2))
-sum(y * log(p2) + (1 - y) * log(1-p2))p1的LogarithmicLoss遠小於p2，代表它的預測準確率勝過p2。
這是一個標準的overfitting的範例。m2在trainingdataset上完美的詮釋了類別的變化，但是在testingdataset上表現卻不如m1。這表示說，並不是描述trainingdataset越好，在其他資料上的表現就會越好。這個現象，讓機器學習的問題更具挑戰。
事實上，我們可以選到更好的模型組合。請同學挑戰看看。
# 方便起見，同學可以使用這個函數計算 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))
}
# 請找出一個在df.test上的logloss小於24.5的模型
answer_02 <- local({
NULL
# 請填寫你的程式碼
step(m1, trace = 0)
})
stopifnot(class(answer_02) == c("glm", "lm"))
stopifnot(logloss(df.test$Class == "good", predict(answer_02, df.test, type = "response")) < 24.5)