這個課程跟同學介紹如何用R自資料中建立線性模型、學習參數與做預測和選模。
在R中要建立線性模型，首先要把資料給整理好。舉例來說，先把資料整理成data.frame是很常見的作法。
舉例來說，請同學輸入：cars看一看這個車速與煞車距離的資料集。
cars我們想要研究車速（speed）與煞車後的滑行距離（dist）之間的關係。而這些資料，都已經放到cars資料集合之中了。
cars的資料中，每一欄都代表一個面向的數據。舉例來說，第一欄位都是車速、第二欄位都是煞車的滑行距離。
cars的每一列，則代表一筆觀測資料。舉例來說，第一列的資料顯示，第一筆觀測數據，車速是4mph，而煞車後的滑行距離則是2ft。這些配對關係是很重要，第一列的車速就是配對到第一列的煞車滑行距離。
在建立模型之前，建議先盡可能的看看資料的長相。請同學先輸入skip()自動檢查有無安裝ggplot2套件。
check_then_install("ggplot2", "2.0.0")接著，請載入ggplot2套件
library(ggplot2)請輸入：g<-ggplot(cars,aes(x=speed,y=dist))+geom_point()來建立ggplot2的繪圖物件，並且把結果寫入到變數g
g <- ggplot(cars, aes(x = speed, y = dist)) + geom_point()請輸入print(g)看看上一題答案的結果
print(g)由於兩筆數據間有呈現線性關係的趨勢，所以用線性模型應該是合理的選擇。請同學輸入：?lm看看建立線性模型的說明文件
?lmlm的第一個參數是formula，是R中常常用來描述數據間關係的語法。data代表的是R要取得資料的data.frame。請同學輸入：m<-lm(dist~speed,cars)，將lm的輸出寫入m
m <- lm(dist ~ speed, cars)上一個答案中的dist~speed是R的formula物件。formula物件的目的是描述數據間的關係。~左方的名稱，是應變變數，也是我們感興趣的變數名稱；~右方的名稱，是獨立變數，我們將拿獨立變數來解釋應變變數的變化。
請同學輸入m看一看結果
m請同學輸入class(m)看看m的型態
class(m)R中計算而得的數據模型，通常都會有自己的型態。例如lm產生的線性模型物件，型態就是：“lm”(linearmodel)。
請同學輸入：mode(m)，看看m真正的型態
mode(m)lm輸出的物件，本質上只是一個帶有名字的list。在說明文件中，副標題為Value段落，解釋了該list中每個元素的名字與意義。
請同學輸入：coef(m)
coef(m)線性模型所學到的參數，都紀錄在coef(m)的輸出中。我們可以看到一個帶有名字的數值向量（numericvector）。這裡的Intercept，代表的就是迴歸線的截距，也就是當車速為0時的煞車滑動距離。而speed，則代表speed變數的參數。也就是說，如果我們知道車速為x時，線性模型預測滑行的距離為：-17.58+3.93*x
coef(m)我們可以使用：g+stat_smooth(method="lm",formula=y~x)來畫出計算出的迴歸線。
g + stat_smooth(method = "lm", formula = y ~ x)在geom_smooth中，method參數代表我們要建立模型的指令，也就是lm。只是這裡要給字串（要用雙引號包覆lm）。而formula=y~x代表這條線所描述的關係。圖形中的灰色區域，代表預測值的95%信賴區間。
summary(m)則鉅細靡遺的描述線性模型的性質。請同學試試看。
summary(m)Residuals區以敘述統計描述模型誤差的分佈。所謂的Residuals，就是在學習的資料（也就是lm的參數data所對應到的data.frame，例：cars）中，觀測值與預測值得差距。
Coefficients則顯示的各參數檢定結果(虛無假設為對應的參數值為0)，星號表示顯著，即代表有明顯正顯示該參數不為0。
最下面則顯示各種重要的統計指標：ResidualStandarderror代表的是模型誤差的變異數；
R-squared則類似相關係數的平方，描述資料接不接近線性模型。如果R-squared很靠近一，表示線性模型很接近資料；如果靠近0，則代表模型和資料不像。
最後的F-statistic則用來檢定此模型和只有平均數的模型相比，是不是有明顯的差異。
接著我們講解如何使用模型來作預測。請同學輸入?predict看看predict函數的說明。
?predictR的predict是可用於各種模型的函數（在軟體工程中，又稱作「介面」）。除了線性模型以外的模型，也可能可以用predict函數來做預測。在Usage中，我們看到：predict(object,...)表示predict的第一個參數（名稱為object），一定是模型物件。舉例來說，我們手上的m物件就是代表線性模型，在使用predict函數時是第一個參數或名稱要為object。
如何查詢所有可以和predict搭配使用的模型呢？請輸入：methods(predict)
methods(predict)同學有沒有看到predict.lm這個名稱顯示在螢幕上呢？還記得class(m)的結果是lm嗎？R在處理predict(m,...)這樣的語法時，因為predict是一個介面函數（在R中稱為GenericFunction），所以R會自動轉換成呼叫：predict.lm(m,...)，把第一個參數的型態（這裡是:lm）接到predict之後。這種把型態接到GenericFunction的名稱之後的作法，在R中稱為S3物件導向。
既然我們知道predict(m,...)實際上是predict.lm(m,...)。要查詢我們目前可用的參數就需要輸入：?predict.lm。請同學試試看。
?predict.lm同學這次可以看到許多可以設定的參數了。我們只針對最重要的newdata做說明。這裡的newdata如果不給的話，predict.lm會拿m裡面剛剛學習用的資料（也就是cars，之後我們稱之為「trainingdataset」，因為這個資料是用於做訓練的）
如果我們要預測車速為20時滑行距離為多少時，就需要利用data.frame指令來建立暫時的data.frame物件，告訴R我們想預測的情境是車速20。
接著，同學如果要預測車速為20時，煞車滑行距離為多少時，可以輸入：predict(m,data.frame(speed=20))
predict(m, data.frame(speed = 20))predict.lm也可以輸出其他資訊，例如預測的標準差。有興趣的同學可以細讀Arguments段落與Value段落。
R也可以輸出預測的信賴區間。請同學輸入：predict(m,data.frame(speed=20),interval="predict")標準差的信心水準則是透過參數level控制。
predict(m, data.frame(speed = 20), interval = "predict")請問同學，不設定level參數時，上一題答案所計算的信心水準是多少？
0.95
我們上面介紹了如何用R的lm來建立線性模型，與運用predict來計算預測結果。接著我們來介紹調校模型與選擇模型的方法。
我們要練習的資料集是iris，請同學輸入適當的指令從套件載入該資料集。
data(iris)請同學輸入?iris看看資料的背景
?iris我們要利用iris中的其他資料，來預測Sepal.Length欄位的值
iris是一個data.frame。請同學輸入：sapply(iris,class)看一看每個欄位的型態
sapply(iris, class)請問唯一屬於類別型變數的欄位名稱是？請輸入字串
Species
因為Sepal.Length是數值型變數，所以我們可以使用lm做學習。
使用formula時，如果輸入：xxx~.就代表xxx是目標變數，而其他欄位通通都當成解釋變數。請同學輸入：m.iris1<-lm(Sepal.Length~.,iris)
m.iris1 <- lm(Sepal.Length ~ ., iris)請用coef(m.iris1)列出所有學到的模型參數。
coef(m.iris1)我們想要比較coef(m.iris1)是不是大部分都屬於iris的欄位名稱。請同學先輸入：tmp1<-names(coef(m.iris1))
tmp1 <- names(coef(m.iris1))接著請同學輸入：tmp2<-colnames(iris)取得iris的欄位名稱。
tmp2 <- colnames(iris)最後我們運用setdiff（集合的減法）來找出tmp1中不屬於tmp2的元素。請同學輸入：setdiff(tmp1,tmp2)
setdiff(tmp1, tmp2)除了(Intercept)之外，同學應該還會看到Speciesversicolor與Speciesvirginica。這是因為Species是類別型變數，所以R把欄位名稱與類別名稱接在一起當成參數名稱。另外考量到解釋性，Speciessetosa（欄位名稱+第一個類別名稱）會被移除，只留下剩下的兩個類別。
線性模型處理類別型變數的概念，很接近「分組取平均」。
我們請同學複習一下dplyr的group_by與summarise。
check_then_install("dplyr", "0.4.3")
library(dplyr)
# 請計算出iris的三種花的種類，所帶有的平均Sepal.Length
answer_01 <- local({
data(iris)
group_by(iris, Species) %>%
summarise(mean(Sepal.Length))
})
stopifnot(is.data.frame(answer_01))
stopifnot(colnames(answer_01) == c("Species", "mean(Sepal.Length)"))
stopifnot(sum(answer_01[[2]]) == 17.53)
# 完成後，請同學存檔並回到console輸入：`submit()`為了方便比較，請同學輸入：answer_01把上一題的答案輸出到console上
answer_01接著請同學輸入：lm(Sepal.Length~Species,iris)，看看學習出來的結果。
lm(Sepal.Length ~ Species, iris)模型中的參數(Intercept)，其實就是setosa中Sepal.Length的平均值；模型中的參數Speciesversicolor則是versicolor與setosa的Sepal.Length的平均值得差；Speciesvirginica同理。
我們也可以透過適當的設定formula來讓模型更複雜。請同學輸入：m.iris2<-lm(Sepal.Length~.^2,iris)這種作法除了考量到各種欄位對Sepal.Length的個別影響之外，欄位的影響力也會受到其它欄位的影響。
m.iris2 <- lm(Sepal.Length ~ .^2, iris)請同學輸入m.iris2看看這個模型的參數
m.iris2除了原本在m.iris1看到的參數之外，還多了許多中間有:的參數。舉例來說，Sepal.Width:Petal.Width代表的是Sepal.Width與Petal.Width對Sepal.Length（目標變數）的影響力，會受到彼此的影響。Sepal.Width:Petal.Width這個參數，在統計上稱為：Sepal.Width與Petal.Width的「交互作用」。或是「二次交互作用項」。因為這類的交互作用，是可以擴充到三、甚至四個變數的狀況。
在formula中的Sepal.Length~.^2中，.就代表除了Sepal.Length之外所有在iris中的欄位。而^2則代表所有這些欄位之間的兩兩組合（二次交互作用項），都要放到模型之中。同理，^3就代表所有的三次交互作用項。我們也可以手動指定要放入模型的交互作用，例如：Sepal.Length~.+Sepal.Width:Petal.Width就代表除了m.iris1的參數之外，我們只增加了Sepal.Width:Petal.Width這一個二次交互作用項。
我們可以一直放高次的交互作用到模型之中，讓模型變得非常複雜。結果導致我們能夠計算出很多個模型，但是哪一個才是「最好用」的模型呢？
我們可以用指令residuals(m.iris1)來看到「實際的Sepal.Length」與「m.iris1預測的Sepal.Length」的差距。這些在trainingdataset上看到的模型與實際觀測值之間的差距，稱為residuals。請同學試試看這個指令
residuals(m.iris1)透過比較residuals的大小，我們可以知道模型在trainingdataset描述Sepal.Length的變化上，做的好不好。請輸入：sum(residuals(m.iris1)^2)，計算residuals的平方和
sum(residuals(m.iris1)^2)請輸入：sum(residuals(m.iris2)^2)，計算第二個模型的residuals的平方和
sum(residuals(m.iris2)^2)請問哪一個模型的residuals的平方和比較小？描述trainingdataset的變化描述的比較好？
m.iris2
事實上，模型越複雜，trainingdataset上的Residuals的平方和就越小。
請同學計算一個有所有三次交互作用項的模型，residuals的平方和為多少。同學可以藉此機會確認，這個模型的residual的平方和，會小於m.iris2的residuals的平方和。
# 請同學計算一個有所有三次交互作用項的模型
# Residuals的平方和為多少。
answer_02 <- local({
m.iris3 <- lm(Sepal.Length ~ .^3, iris) # 請修改這段程式碼
sum(residuals(m.iris3)^2)
})
stopifnot(class(answer_02) == "numeric")
stopifnot(length(answer_02) == 1)
# 完成後，請同學存檔並回到console輸入：`submit()`
事實上，透過residuals的平方和，我們可以計算出RSquared。RSquared是很多統計學課程中會提到的「描述線性模型表現」的指標（我們略過嚴謹的數學定義）。數學上可以證明，RSquared的值會介於0與1之間。0代表模型完全沒有解釋到目標變數在trainingdataset上的變化；1代表模型100%解釋了目標變數在trainingdataset的行為（同時也代表residuals都是0）。RSquared等價於「相關係數的平方」。請同學輸入summary(m.iris2)，這個詳細的表格會回報模型的RSquared。
summary(m.iris2)在RBasic中有介紹如何取出RSquared的值，所以我們就不再贅述了。
實務上，RSquared並不代表模型比較「好用」（例如：預測的比較準）。模型在trainingdataset上的表現，「不代表」模型在trainingdataset以外的資料上會表現的比較好。透過方法挑選出在trainingdataset以外的資料上比較好的模型，稱為「選模」。
統計學提供了一些統計指標，在只使用trainingdataset的狀況下，去預測模型在在trainingdataset以外的資料上的表現。R則有提供相關的實作，甚至連選模的過程都實作了。同學可以試試看：m.iris<-step(m.iris2,steps=1)
m.iris <- step(m.iris2, steps = 1)R會去計算AIC(TheAkaikeinformationcriterion)來評估模型在testingdataset的好壞。AIC越小，模型就越好。接著R會嘗試去拿掉一個欄位、交互作用項或是新增一個欄位、交互作用項，看看模型的AIC會不會提昇。我們可以看到，當拿掉：Sepal.Width:Species時，AIC會變低（原本的AIC，即是m.iris2的解釋變數組合，從「所有的二次交互作用」(.^2)變成「所有的二次交互作用扣除Sepal.Width:Species」(.^2-Sepal.Width:Species)。
m.iris<-step(m.iris2)則會一直重複上述的動作，直到模型的AIC不會變得更低。並把結果的模型輸出（寫入變數m.iris）。請同學試試看。
m.iris <- step(m.iris2)最後我們請同學利用另一個資料做練習。
check_then_install("mlbench", "2.1.1") # 請同學安裝套件mlbench，版本為2.1.1或更新
library(mlbench)
data(BostonHousing)
# 同學應該先看看資料。請輸入：
# ?BostonHousing # 了解資料的背景
# sapply(BostonHousing, class) # 看看每欄的型態
# summary(BostonHousing) # 看看數字分佈
# nrow(BostonHousing) # 看看資料個數
# 請問這份資料的目標是哪一個變數？請閱讀說明文件後作答
answer_03_01 <- "medv"
# 第一回合：請同學做出一個R Squared超過0.9的模型
answer_03_02 <- local({
lm(medv ~ .^2, BostonHousing)
})
stopifnot(class(answer_03_02) == "lm")
stopifnot(summary(answer_03_02)$r.squared > 0.9)
# 第二回合：請同學做出一個AIC低於1089的模型
answer_03_03 <- local({
step(answer_03_02, trace = 0)
})
stopifnot(class(answer_03_03) == "lm")
stopifnot(extractAIC(answer_03_03)[2] < 1089)
# 第三回合：請同學做出一個AIC低於-520模型
# (本題為挑戰題，submit後不檢查)
answer_03_04 <- local({
b.breaks <- c(-Inf, quantile(BostonHousing$b))
df <-
BostonHousing %>%
mutate(b.cut = cut(b, breaks = b.breaks))
answer_03_04 <-
lm(formula = medv ~ .^3, data = df)
})