探索預測模式
本節介紹兩種迴歸預測模型：多變量線性迴歸(multiple linear regression)及迴歸樹(regression trees)。這兩種方法是完全不同形式的迴歸作法，而且很容易說明及執行。在真實的資料採礦項目中還是有許多其他的方法可供採用。這兩種方法對於缺漏值的處理方式不同，多線性迴歸不允許有缺漏值，所以必須先處理缺漏值的問題。迴歸樹則不在乎缺漏值，所以可以直接使用原始資料集。本節會利用前述200筆資料來探索預測模式，並用這些模式預測另外140筆的測試資料。
多變量線性迴歸
多變量線性迴歸是種常用資料分析的技術。這種模式是由多個變量Xi乘上特定係數βi加總後的結果。
================================
//首先處理缺漏值：
> data(algae)
> algae <- algae[-manyNAs(algae), ]
> clean.algae <- knnImputation(algae,k=10)
//使用lm()取得迴歸模型
> lm.a1 <- lm(a1 ~ ., data = clean.algae[, 1:12])
================================
lm()第一個參數a1 ~ .表示a1依據第二個參數指定的欄位統計迴歸模型，第二個參數指定前12欄位最為迴歸統計之來源。（注意第12個欄位是a1必須包含）如果你想依據某幾個欄位來分析第一個參數可以改寫成a1∼mxPH + NH4。'.'則代表第二參數的所有欄位。
如果想看此迴歸模型的內容可使用summary():
================================
================================
summary()會給出一些診斷的資訊，首先是線性模型是否符合的殘差(residuals)，這個殘差必須是平均0（或愈接近0）愈表示模型可用，而且是一個常態分配。
對於每一個變量，R顯示其值及標準差。
要驗證每一變量的重要度，先假設其參數值為H0:βi=0(null)，然後運用t-test來測試這個假設。t=βi/(標準差)，最後算出Pr(>|t|)這一欄，其值如果是0.0001表示我們有99.99％的信心這個參數不是0(null)。如果結果不到90％的信心度（一個星號以上），βi=0(null)便為真。
另一個診斷指標R²(multiple 及 adjusted)表示模型的適合度，亦即模型資料的變異比率(proportion of variance)。其值愈接近1愈佳（幾乎100％呈現其變動），值愈小，適合度愈差。adjusted係數更是嚴苛，因為它參考數個迴歸模型的參數。
最後我們也測試null假設：所有變數皆無關，即H0:β1=
β2=...=βm=0。F-statistic即代表這個假設。若p-level值為0.0001即表示有99.99％的信心度null假設為非。即可以排除null假設。如果p-level。0.1表示無法通過測試，上述的t-test便無意義了。
有些診斷方式可以採用繪圖方式。如plot(lm.a1):
 |
每一個適合變數值相對殘差值
標出編號者表示殘差太大 |
 |
殘差是否符合常態分配
分配是否貼近直線 |
這個模型資料的變異比率(proportion of variance)不是很顯著（R² adjusted約32.0％）。因為F-statistic的p-level值非常小，故可以推翻目標變數與預測不相依的假設。
由於部分變數不是很顯著相關，必須將其從模型中排除。有許多方法可供選擇，下面是使用後退淘汰法(backward elimination)，我們使用stats套件中的anova()來處理：
執行結果由於season的Pr(>F)接近1，代表其對於降低適合度殘差貢獻度最差。所以必需移除：
================================
> lm2.a1 <- update(lm.a1,.~.-season)
================================
update()是stats套件的方法。可以變動現有的線性模型。lms.a1即是從lm.a1中移除season變數。我們可以看看執行結果：
變異比率稍有改善但還不是很顯著（R² adjusted約32.8％），我們可以再次使用anova()來比對前後兩次的差距：
anova()使用F-test分析兩者的變異是否顯著，在這個案例中，雖然方差和(sum of the squared errors)降低(-447.62)，但比對結果變異不是恨顯著（Pr(>F)=0.6971約只有30％的信心度說兩者有差異）。我們還可以再次執行直到移除所有的變數。但我們可以使用stats套件中的step()方法，它會幫我們自動反覆執行：
我們可以檢視執行結果：
最後的結果顯示變異比率(33%)仍然不理想，表示此領域問題並不適合線性假設。
迴歸樹
這是另一種迴歸模型，執行下列指令：
================================
> library(rpart)
> data(algae)
> algae <- algae[-manyNAs(algae), ]
> rt.a1 <- rpart(a1 ~ ., data = algae[, 1:12])
================================
執行rpart()產生一迴歸樹，其參數與前述的lm()一樣。檢視rt.a1的內容：
================================
 |
split：分支條件
n：樣本數
deviance：與平均值偏差的平方和
yval：a1的平均值 |
================================
迴歸樹是解釋性(explanatory)變量邏輯推論的層級架構(hierarchy)。它會自動選取相關的變數，因此並不是所有的變數都會用到。閱讀方式，從root節點開始，然後依數字由小而大。每個數字代表一個節點。除了末節點，每一個節點都有兩個分支，每個分支都有其進入的條件。只要從root開始，依條件進行，直到末節點其最終結果（預測值）便是yval。我們可以使用graphics套件的plot(rt.a1)及text(rt.a1)或DMwR套件的prettyTree(rt.a1)來繪製樹狀圖：
 |
| 用Plot()繪製的樹狀圖 |
 |
| 用text()繪製的樹狀圖 |
 |
| 用prettyTree()繪製的樹狀圖 |
我們也可以用summary()來檢視其內容：
================================
> summary(rt.a1)
----------------------------------------------------------
Call:
rpart(formula = a1 ~ ., data = algae[, 1:12])
n= 198
CP nsplit rel error xerror xstd
1 0.40573990 0 1.0000000 1.0089717 0.1308742
2 0.07188523 1 0.5942601 0.6757552 0.1183296
3 0.03088731 2 0.5223749 0.6591388 0.1164799
4 0.03040753 3 0.4914876 0.7006920 0.1219898
5 0.02787181 4 0.4610800 0.7020964 0.1220327
6 0.02775354 5 0.4332082 0.7020964 0.1220327
7 0.01812406 6 0.4054547 0.6953859 0.1139465
8 0.01634372 7 0.3873306 0.6959028 0.1125755
9 0.01000000 9 0.3546432 0.6957636 0.1127715
Variable importance
PO4 oPO4 NH4 Cl mxPH Chla NO3 mnO2 size season speed
25 20 15 15 9 7 3 2 1 1 1
Node number 1: 198 observations, complexity param=0.4057399
mean=16.99646, MSE=456.5722
left son=2 (147 obs) right son=3 (51 obs)
Primary splits:
PO4 < 43.818 to the right, improve=0.4048567, (1 missing)
oPO4 < 18.889 to the right, improve=0.3793450, (0 missing)
NH4 < 51.27 to the right, improve=0.3625269, (0 missing)
Cl < 7.2915 to the right, improve=0.3583409, (8 missing)
Chla < 1.15 to the right, improve=0.2533869, (10 missing)
Surrogate splits:
oPO4 < 17.5415 to the right, agree=0.944, adj=0.78, (1 split)
NH4 < 37.639 to the right, agree=0.893, adj=0.58, (0 split)
Cl < 9.0275 to the right, agree=0.858, adj=0.44, (0 split)
Chla < 1.05 to the right, agree=0.822, adj=0.30, (0 split)
mxPH < 7.295 to the right, agree=0.817, adj=0.28, (0 split)
Node number 2: 147 observations, complexity param=0.07188523
mean=8.979592, MSE=212.7831
left son=4 (140 obs) right son=5 (7 obs)
Primary splits:
Cl < 7.8065 to the right, improve=0.2071337, (1 missing)
Chla < 1.15 to the right, improve=0.1959676, (1 missing)
oPO4 < 51.118 to the right, improve=0.1651094, (0 missing)
NH4 < 49.25 to the right, improve=0.1494842, (0 missing)
PO4 < 125 to the right, improve=0.1393822, (0 missing)
Surrogate splits:
Chla < 0.6 to the right, agree=0.959, adj=0.143, (1 split)
Node number 3: 51 observations, complexity param=0.03040753
mean=40.10392, MSE=440.0541
left son=6 (28 obs) right son=7 (23 obs)
Primary splits:
mxPH < 7.87 to the left, improve=0.12171490, (1 missing)
PO4 < 6.35 to the right, improve=0.10576260, (1 missing)
Cl < 7.544 to the right, improve=0.10428070, (7 missing)
NH4 < 18.381 to the right, improve=0.10356000, (0 missing)
oPO4 < 10.625 to the right, improve=0.09644168, (0 missing)
Surrogate splits:
size splits as RRL, agree=0.78, adj=0.522, (1 split)
NO3 < 1.1875 to the right, agree=0.74, adj=0.435, (0 split)
oPO4 < 3.111 to the left, agree=0.70, adj=0.348, (0 split)
season splits as LLRR, agree=0.60, adj=0.130, (0 split)
NH4 < 22.0355 to the left, agree=0.60, adj=0.130, (0 split)
Node number 4: 140 observations, complexity param=0.03088731
mean=7.492857, MSE=154.4488
left son=8 (84 obs) right son=9 (56 obs)
Primary splits:
oPO4 < 51.118 to the right, improve=0.12913450, (0 missing)
PO4 < 125 to the right, improve=0.09908251, (0 missing)
NH4 < 41.875 to the right, improve=0.05847356, (0 missing)
NO3 < 3.2725 to the right, improve=0.05343570, (0 missing)
Chla < 3.65 to the right, improve=0.04761161, (1 missing)
Surrogate splits:
PO4 < 125 to the right, agree=0.857, adj=0.643, (0 split)
Cl < 27.8665 to the right, agree=0.721, adj=0.304, (0 split)
NO3 < 3.313 to the right, agree=0.679, adj=0.196, (0 split)
mnO2 < 9.5 to the left, agree=0.664, adj=0.161, (0 split)
season splits as RLLL, agree=0.657, adj=0.143, (0 split)
Node number 5: 7 observations
mean=38.71429, MSE=451.1098
Node number 6: 28 observations, complexity param=0.02775354
mean=33.45, MSE=409.0275
left son=12 (18 obs) right son=13 (10 obs)
Primary splits:
mxPH < 7.045 to the right, improve=0.2296931, (1 missing)
PO4 < 6.25 to the right, improve=0.2174386, (1 missing)
oPO4 < 12.375 to the right, improve=0.1721865, (0 missing)
NH4 < 17.1 to the right, improve=0.1098949, (0 missing)
season splits as LRRR, improve=0.0944271, (0 missing)
Surrogate splits:
NH4 < 11.25 to the right, agree=0.852, adj=0.556, (1 split)
oPO4 < 1.125 to the right, agree=0.852, adj=0.556, (0 split)
PO4 < 6.5835 to the right, agree=0.852, adj=0.556, (0 split)
speed splits as L-R, agree=0.778, adj=0.333, (0 split)
NO3 < 1.9675 to the left, agree=0.778, adj=0.333, (0 split)
Node number 7: 23 observations, complexity param=0.02787181
mean=48.20435, MSE=358.3091
left son=14 (12 obs) right son=15 (11 obs)
Primary splits:
PO4 < 15.177 to the right, improve=0.3057413, (0 missing)
NH4 < 20.4165 to the right, improve=0.2692864, (0 missing)
Cl < 7.544 to the right, improve=0.2055829, (0 missing)
Chla < 0.85 to the right, improve=0.1534699, (1 missing)
oPO4 < 6.25 to the right, improve=0.1013330, (0 missing)
Surrogate splits:
NH4 < 20.4165 to the right, agree=0.913, adj=0.818, (0 split)
Cl < 5.8595 to the right, agree=0.826, adj=0.636, (0 split)
NO3 < 1.353 to the right, agree=0.826, adj=0.636, (0 split)
oPO4 < 5 to the right, agree=0.783, adj=0.545, (0 split)
Chla < 0.85 to the right, agree=0.739, adj=0.455, (0 split)
Node number 8: 84 observations
mean=3.846429, MSE=40.96606
Node number 9: 56 observations, complexity param=0.01812406
mean=12.9625, MSE=274.8113
left son=18 (24 obs) right son=19 (32 obs)
Primary splits:
mnO2 < 10.05 to the right, improve=0.10646520, (0 missing)
PO4 < 101.894 to the left, improve=0.08815216, (0 missing)
oPO4 < 24.3335 to the left, improve=0.07637520, (0 missing)
size splits as LLR, improve=0.06017653, (0 missing)
mxPH < 8.35 to the right, improve=0.05440345, (0 missing)
Surrogate splits:
PO4 < 101.894 to the left, agree=0.750, adj=0.417, (0 split)
size splits as LRR, agree=0.696, adj=0.292, (0 split)
season splits as LRRR, agree=0.679, adj=0.250, (0 split)
NH4 < 89.8 to the left, agree=0.661, adj=0.208, (0 split)
mxPH < 8.025 to the right, agree=0.643, adj=0.167, (0 split)
Node number 12: 18 observations
mean=26.39444, MSE=285.8983
Node number 13: 10 observations
mean=46.15, MSE=379.7645
Node number 14: 12 observations
mean=38.18333, MSE=253.9597
Node number 15: 11 observations
mean=59.13636, MSE=243.086
Node number 18: 24 observations
mean=6.716667, MSE=52.02806
Node number 19: 32 observations, complexity param=0.01634372
mean=17.64687, MSE=390.6975
left son=38 (9 obs) right son=39 (23 obs)
Primary splits:
NO3 < 3.1875 to the right, improve=0.09580105, (0 missing)
Chla < 2.55 to the left, improve=0.08399898, (0 missing)
oPO4 < 24.917 to the left, improve=0.07524892, (0 missing)
mnO2 < 9.4 to the left, improve=0.06578127, (0 missing)
Cl < 43.7085 to the right, improve=0.04807023, (0 missing)
Surrogate splits:
mxPH < 7.55 to the left, agree=0.844, adj=0.444, (0 split)
NH4 < 224.643 to the right, agree=0.812, adj=0.333, (0 split)
PO4 < 206.7225 to the right, agree=0.812, adj=0.333, (0 split)
Cl < 84.0465 to the right, agree=0.750, adj=0.111, (0 split)
Node number 38: 9 observations
mean=7.866667, MSE=28.56444
Node number 39: 23 observations, complexity param=0.01634372
mean=21.47391, MSE=480.3263
left son=78 (13 obs) right son=79 (10 obs)
Primary splits:
mnO2 < 8 to the left, improve=0.15906320, (0 missing)
PO4 < 118.6 to the left, improve=0.10091960, (0 missing)
NH4 < 168.75 to the left, improve=0.07651249, (0 missing)
NO3 < 1.2495 to the left, improve=0.07260629, (0 missing)
mxPH < 8.26 to the right, improve=0.06930695, (0 missing)
Surrogate splits:
season splits as RLLL, agree=0.696, adj=0.3, (0 split)
size splits as LLR, agree=0.696, adj=0.3, (0 split)
speed splits as RLL, agree=0.696, adj=0.3, (0 split)
Cl < 46.9725 to the right, agree=0.696, adj=0.3, (0 split)
NH4 < 216.653 to the left, agree=0.696, adj=0.3, (0 split)
Node number 78: 13 observations
mean=13.80769, MSE=224.5807
Node number 79: 10 observations
mean=31.44, MSE=637.0704
================================
樹的建構有兩個步驟：一開始先建一個大樹，然後依據統計評估從底部修剪。因為大樹可以完美的符合訓練資料，但是可能擷取到假性的關連性，因而應用在新資料時，效能可能不佳 。這種現象稱為『過度擬合(overfitting)』。過度擬合在建模時常發生，特別是假設太過寬鬆時。
rpart()建立的樹是依據下列三原則決定其規模：1.偏差低於某門檻。2.節點中的樣本數低於某門檻。3.樹的深度高於某值。這三者是依據下列參數：cp、minsplit、maxdepth，其預設值分別為0.01、20、30。
printcp()透過cp值產生一群子樹，並評估這些子樹的效能：
================================
>printcp(rt.a1)
----------------------------------------------------------
Regression tree:
rpart(formula = a1 ~ ., data = algae[, 1:12])
Variables actually used in tree construction:
[1] Cl mnO2 mxPH NO3 oPO4 PO4
Root node error: 90401/198 = 456.57
n= 198
CP nsplit rel error xerror xstd
1 0.405740 0 1.00000 1.01123 0.130994
2 0.071885 1 0.59426 0.76443 0.115210
3 0.030887 2 0.52237 0.71454 0.114962
4 0.030408 3 0.49149 0.69195 0.106685
5 0.027872 4 0.46108 0.68825 0.106771
6 0.027754 5 0.43321 0.68398 0.107168
7 0.018124 6 0.40545 0.67585 0.104675
8 0.016344 7 0.38733 0.69122 0.105072
9 0.010000 9 0.35464 0.69876 0.099859
================================
註：由於是採隨機運算，因此上列xerror及xstd數字每次執行可能不同。
rpart()產生的是其中第9回，其cp=0.01。但第7回的相對誤差最低，是比較佳的選擇。
另一種選擇的規則是1-SE，選擇最小的樹（第2回）。誤差0.76443小於0.67585+0.104675=0.780525。
當決定所採用的樹，便可以透過cp值來修剪：
================================
> rt2.a1 <- prune(rt.a1, cp = 0.08)
> rt2.a1
----------------------------------------------------------
n= 198
node), split, n, deviance, yval
* denotes terminal node
1) root 198 90401.29 16.996460
2) PO4>=43.818 147 31279.12 8.979592 *
3) PO4< 43.818 51 22442.76 40.103920 *
================================
上述建立樹及修剪樹的過程可以用rpartXse()一次完成：
================================
> (rt.a1 <- rpartXse(a1 ~ ., data = algae[, 1:12]))
----------------------------------------------------------
n= 198
node), split, n, deviance, yval
* denotes terminal node
1) root 198 90401.29 16.996460
2) PO4>=43.818 147 31279.12 8.979592 *
3) PO4< 43.818 51 22442.76 40.103920 *
================================
R也可以讓你自行決定修剪樹的哪些節點：
================================
> first.tree <- rpart(a1 ~ ., data = algae[, 1:12])
> snip.rpart(first.tree, c(4, 7))
----------------------------------------------------------
n= 198
node), split, n, deviance, yval
* denotes terminal node
1) root 198 90401.290 16.996460
2) PO4>=43.818 147 31279.120 8.979592
4) Cl>=7.8065 140 21622.830 7.492857 *
5) Cl< 7.8065 7 3157.769 38.714290 *
3) PO4< 43.818 51 22442.760 40.103920
6) mxPH< 7.87 28 11452.770 33.450000 *
7) mxPH>=7.87 23 8241.110 48.204350
14) PO4>=15.177 12 3047.517 38.183330 *
================================
snip.rpart()第二個參數是欲修剪的節點編號。
你也可以在圖形上直接點選欲修剪的節點：
================================
> prettyTree(first.tree)
> snip.rpart(first.tree)
----------------------------------------------------------
node number: 2 n= 147
response= 8.979592
Error (dev) = 31279.12
node number: 6 n= 28
response= 33.45
Error (dev) = 11452.77
n= 198
node), split, n, deviance, yval
1) root 198 90401.290 16.996460
2) PO4>=43.818 147 31279.120 8.979592 *
3) PO4< 43.818 51 22442.760 40.103920
6) mxPH< 7.87 28 11452.770 33.450000 *
7) mxPH>=7.87 23 8241.110 48.204350
14) PO4>=15.177 12 3047.517 38.183330 *
15) PO4< 15.177 11 2673.945 59.136360 *
================================
以上是點擊2、6兩個節點的結果。
模型的評估與選擇
評估模型最常用的標準是模型的預測能力(predictive performance of the models)，其評估方式是使用真實資料預測並比對目標值，計算出其平均誤差。平均絕對誤差(mean absolute error (MAE))即是其中之一。在R中使用predict()來獲得預測結果：
================================
> lm.predictions.a1 <- predict(final.lm, clean.algae)
> rt.predictions.a1 <- predict(rt.a1, algae)
================================
predict(模型,資料集)
接著計算平均絕對誤差：
================================
> (mae.a1.lm <- mean(abs(lm.predictions.a1 - algae[, "a1"])))
----------------------------------------------------------
[1] 13.10681
================================
> (mae.a1.rt <- mean(abs(rt.predictions.a1 - algae[, "a1"])))
----------------------------------------------------------
[1] 11.61717
================================
另一個常用的誤差評估值是均方誤差(mean squared error (MSE))
================================
> (mse.a1.lm <- mean((lm.predictions.a1 - algae[, "a1"])^2))
----------------------------------------------------------
[1] 295.5407
================================
> (mse.a1.rt <- mean((rt.predictions.a1 - algae[, "a1"])^2))
----------------------------------------------------------
[1] 271.3226
================================
MSE有一個缺點是評估的目標值並非同一單位。因此較不常用，所以一般採用歸一化均方誤差(normalized mean squared error(NMSE))。這個統計值計算模型效能與預測基準（目標值的平均值）的比值：
================================
> (nmse.a1.lm <- mean((lm.predictions.a1-algae[,'a1'])^2)/mean((mean(algae[,'a1'])-algae[,'a1'])^2))
----------------------------------------------------------
[1] 0.6473034
================================
> (nmse.a1.rt <- mean((rt.predictions.a1-algae[,'a1'])^2)/mean((mean(algae[,'a1'])-algae[,'a1'])^2))
----------------------------------------------------------
[1] 0.5942601
================================
NMSE是無單位的，一般從0到1。如果模型效能佳，即愈趨近預測基準，則NMSE小於1，且值愈小愈好。如果大於1表示比用平均值作為預測值的還差。
我們可以使用regr.eval()各種的迴歸評估指標：
================================
> regr.eval(algae[, "a1"], rt.predictions.a1, train.y = algae[, "a1"])
----------------------------------------------------------
mae mse rmse nmse nmae
11.6171709 271.3226161 16.4718735 0.5942601 0.6953711
================================
我們也可以視覺化模型的效能，如誤差分布圖：
================================
> old.par <- par(mfrow = c(1, 2))
> plot(lm.predictions.a1, algae[, "a1"], main = "Linear Model",xlab = "Predictions", ylab = "True Values")
> abline(0, 1, lty = 2)
> plot(rt.predictions.a1, algae[, "a1"], main = "Regression Tree",xlab = "Predictions", ylab = "True Values")
> abline(0, 1, lty = 2)
> par(old.par)
================================
 |
| 誤差分布圖 |
從上圖觀察，顯示模型效能不佳，因為效能佳的時，誤差值的點應接近虛線[abline(0,1,lty=2)]部分。虛線是代表X軸等於Y軸部分。同樣的也可以使用identify()手動標示個別的點：
================================
> plot(lm.predictions.a1,algae[,'a1'],main="Linear Model",xlab="Predictions",ylab="True Values")
> abline(0,1,lty=2)
> algae[identify(lm.predictions.a1,algae[,'a1']),]
================================
以滑鼠點選以選出差異大的點。按Esc結束後，R會顯示所選點的資料。
從上圖中可以發現有些點是負值，在此案例中的領域規則裏，負值是無意義的，而大部分情況下，應該是0。因此我們可以利用領域規則，加上ifelse()判斷<0的情況來改善這個模型：（請注意評估指標，是否對於模型的效能有所提昇）
================================
> sensible.lm.predictions.a1 <- ifelse(lm.predictions.a1 < 0, 0, lm.predictions.a1)
> regr.eval(algae[, "a1"], lm.predictions.a1, stats = c("mae", "mse"))
-----------------------------------------------------------
mae mse
13.10681 295.54069
================================
> regr.eval(algae[, "a1"], sensible.lm.predictions.a1, stats = c("mae", "mse"))
-----------------------------------------------------------
mae mse
12.48276 286.28541
================================
ifelse()有三個參數，第一個是邏輯判斷，如果此邏輯判斷為True則傳回第二參數，否則傳回第三參數。
從上述的評估，迴歸樹因為NMSE比較低似乎是比較適合用來預估此140個測試樣本。但這樣的推理有個陷阱。因為我們不知道這些樣本的目標便數值，而要去評估哪些模型的效能是比較佳的。重點是從不知道真正目標值的資料中，獲得模型效能可信賴的評估。使用訓練資料計算效能值標不是可信賴的，因為其結果會有偏差。事實上，有些模型可以從訓練資料得到零預測誤差的。但是這樣的效能並不能同樣在目標值未知的新樣本中也有同樣的效能。這個現象就如同前述的過度擬合。
因此需要一些模型可以評估運用於含未知值資料的模型效能。
k-fold交叉驗證(k-fold cross-validation)便是針對小型資料集一種最常用的方法。從訓練資料中隨機取得k個相同數量的子集，對每一個k子集使用其餘k-1個子集建構模型，病理用此k子集評估其效能並記錄模型的效能，然後針對所有其他的子集進行相同的步驟。最後得到k個效能評估。所有這些模型效能評估所用的資料都不是用以建構模型的資料。一般k值都是使用10。有時我們還會執行多次，以獲得更有信賴的評估。
一般面對預測的任務，需要下列的決策：
- 選擇各式模型（也可能是同模型但參數不同）。
- 選擇用於比對模型的評估指標。
- 選擇試驗的方法論，以獲得這些指標可信賴的評估。
我們使用experimentalComparison()來選擇及評估。這個方法有三個參數：1.用以比對的資料集。2.各式模型。3.試驗程序的參數。
使用者提供實作模型的函數，這樣的函數必需針對所提供的訓練資料及測試資料，實作完整的「訓練」+「測試」+「評估」的流程。在函數中，逐一呼叫這些函數，這些函數傳回使用者想要的評估指標（此例使用NMSE）。如下所示：
================================
> cv.rpart <- function(form,train,test,...) {
+ m <- rpartXse(form,train,...)
+ p <- predict(m,test)
+ mse <- mean((p-resp(form,test))^2)
+ c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2))
+}
> cv.lm <- function(form,train,test,...) {
+ m <- lm(form,train,...)
+ p <- predict(m,test)
+ p <- ifelse(p < 0,0,p)
+ mse <- mean((p-resp(form,test))^2)
+ c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2))
+}
================================
所有使用者實作的函數需有前三個參數：1.(form)公式(formula)2.訓練資料3.測試資料。其餘的參數(...)是experimentalComparison用來學習評估之用。其中resp()是依據公式取得資料集中的目標便數值。
接下來便是進入學習及測試流程：
================================
> res <- experimentalComparison(
+ c(dataset(a1 ~ .,clean.algae[,1:12],'a1')),
+ c(variants('cv.lm'),
+ variants('cv.rpart',se=c(0,0.5,1))),
+ cvSettings(3,10,1234))
-----------------------------------------------------------
##### CROSS VALIDATION EXPERIMENTAL COMPARISON #####
** DATASET :: a1
++ LEARNER :: cv.lm variant -> cv.lm.v1
3 x 10 - Fold Cross Validation run with seed = 1234
Repetition 1
Fold: 1 2 3 4 5 6 7 8 9 10
Repetition 2
Fold: 1 2 3 4 5 6 7 8 9 10
Repetition 3
Fold: 1 2 3 4 5 6 7 8 9 10
++ LEARNER :: cv.rpart variant -> cv.rpart.v1
3 x 10 - Fold Cross Validation run with seed = 1234
Repetition 1
Fold: 1 2 3 4 5 6 7 8 9 10
Repetition 2
Fold: 1 2 3 4 5 6 7 8 9 10
Repetition 3
Fold: 1 2 3 4 5 6 7 8 9 10
++ LEARNER :: cv.rpart variant -> cv.rpart.v2
3 x 10 - Fold Cross Validation run with seed = 1234
Repetition 1
Fold: 1 2 3 4 5 6 7 8 9 10
Repetition 2
Fold: 1 2 3 4 5 6 7 8 9 10
Repetition 3
Fold: 1 2 3 4 5 6 7 8 9 10
++ LEARNER :: cv.rpart variant -> cv.rpart.v3
3 x 10 - Fold Cross Validation run with seed = 1234
Repetition 1
Fold: 1 2 3 4 5 6 7 8 9 10
Repetition 2
Fold: 1 2 3 4 5 6 7 8 9 10
Repetition 3
Fold: 1 2 3 4 5 6 7 8 9 10
================================
experimentalComparison()方法的第 一個參數，是由dataset(<formula>,<data frame>,<label>)定義的資料集。第二個參數是由variants()定義的使用者定義包含訓練+測試+評估過程的函數。第三個參數是函數的設定值，第一個參數表示k-flods執行次數，第二個參數指定k值，第三個參數指定隨機種子。
variants()第一個參數是函數名稱，第二個參數是該函數所需之參數。執行時依據第二參數指定的參數數量，執行第一參數指定函數相同的次數，並將第二參數的每一個參數傳入函數中。
我們使用summary()來檢視實際執行結果：
================================
> summary(res)
-----------------------------------------------------------
== Summary of a Cross Validation Experiment ==
3 x 10 - Fold Cross Validation run with seed = 1234
* Data sets :: a1
* Learners :: cv.lm.v1, cv.rpart.v1, cv.rpart.v2, cv.rpart.v3
* Summary of Experiment Results:
-> Datataset: a1
*Learner: cv.lm.v1
nmse
avg 0.7196105
std 0.1833064
min 0.4678248
max 1.2218455
invalid 0.0000000
*Learner: cv.rpart.v1
nmse
avg 0.6440843
std 0.2521952
min 0.2146359
max 1.1712674
invalid 0.0000000
*Learner: cv.rpart.v2
nmse
avg 0.6873747
std 0.2669942
min 0.2146359
max 1.3356744
invalid 0.0000000
*Learner: cv.rpart.v3
nmse
avg 0.7167122
std 0.2579089
min 0.3476446
max 1.3356744
invalid 0.0000000
================================
也可以使用圖形來觀察：
================================
> plot(res)
================================
欲查看某一模型的結果：
================================
> getVariant("cv.rpart.v1", res)
----------------------------------------------------------
Learner:: "cv.rpart"
Parameter values
se = 0
================================
我們也可以同時執行7個水藻變數的預測評估：
================================
> DSs <- sapply(names(clean.algae)[12:18],
+ function(x,names.attrs) {
+ f <- as.formula(paste(x,"~ ."))
+ dataset(f,clean.algae[,c(names.attrs,x)],x)
+}, + names(clean.algae)[1:11])
>DSs
----------------------------------------------------------
$a1
Task Name :: a1
Formula :: a1 ~ .
<environment: 0x000000001071e1b8>
Task Data ::
'data.frame': 198 obs. of 12 variables:
$ a1 : num 0 1.4 3.3 3.1 9.2 15.1 2.4 18.2 25.4 17 ...
$ season: Factor w/ 4 levels "autumn","spring",..: 4 2 1 2 1 4 3 1 4 4 ...
$ size : Factor w/ 3 levels "large","medium",..: 3 3 3 3 3 3 3 3 3 3 ...
$ speed : Factor w/ 3 levels "high","low","medium": 3 3 3 3 3 1 1 1 3 1 ...
$ mxPH : num 8 8.35 8.1 8.07 8.06 8.25 8.15 8.05 8.7 7.93 ...
$ mnO2 : num 9.8 8 11.4 4.8 9 13.1 10.3 10.6 3.4 9.9 ...
$ Cl : num 60.8 57.8 40 77.4 55.4 ...
$ NO3 : num 6.24 1.29 5.33 2.3 10.42 ...
$ NH4 : num 578 370 346.7 98.2 233.7 ...
$ oPO4 : num 105 428.8 125.7 61.2 58.2 ...
$ PO4 : num 170 558.8 187.1 138.7 97.6 ...
$ Chla : num 50 1.3 15.6 1.4 10.5 ...
$a2
Task Name :: a2
Formula :: a2 ~ .
<environment: 0x0000000010317b08>
Task Data ::
'data.frame': 198 obs. of 12 variables:
$ a2 : num 0 7.6 53.6 41 2.9 14.6 1.2 1.6 5.4 0 ...
$ season: Factor w/ 4 levels "autumn","spring",..: 4 2 1 2 1 4 3 1 4 4 ...
$ size : Factor w/ 3 levels "large","medium",..: 3 3 3 3 3 3 3 3 3 3 ...
$ speed : Factor w/ 3 levels "high","low","medium": 3 3 3 3 3 1 1 1 3 1 ...
$ mxPH : num 8 8.35 8.1 8.07 8.06 8.25 8.15 8.05 8.7 7.93 ...
$ mnO2 : num 9.8 8 11.4 4.8 9 13.1 10.3 10.6 3.4 9.9 ...
$ Cl : num 60.8 57.8 40 77.4 55.4 ...
$ NO3 : num 6.24 1.29 5.33 2.3 10.42 ...
$ NH4 : num 578 370 346.7 98.2 233.7 ...
$ oPO4 : num 105 428.8 125.7 61.2 58.2 ...
$ PO4 : num 170 558.8 187.1 138.7 97.6 ...
$ Chla : num 50 1.3 15.6 1.4 10.5 ...
$a3
Task Name :: a3
Formula :: a3 ~ .
<environment: 0x0000000010898510>
Task Data ::
'data.frame': 198 obs. of 12 variables:
$ a3 : num 0 4.8 1.9 18.9 7.5 1.4 3.2 0 2.5 0 ...
$ season: Factor w/ 4 levels "autumn","spring",..: 4 2 1 2 1 4 3 1 4 4 ...
$ size : Factor w/ 3 levels "large","medium",..: 3 3 3 3 3 3 3 3 3 3 ...
$ speed : Factor w/ 3 levels "high","low","medium": 3 3 3 3 3 1 1 1 3 1 ...
$ mxPH : num 8 8.35 8.1 8.07 8.06 8.25 8.15 8.05 8.7 7.93 ...
$ mnO2 : num 9.8 8 11.4 4.8 9 13.1 10.3 10.6 3.4 9.9 ...
$ Cl : num 60.8 57.8 40 77.4 55.4 ...
$ NO3 : num 6.24 1.29 5.33 2.3 10.42 ...
$ NH4 : num 578 370 346.7 98.2 233.7 ...
$ oPO4 : num 105 428.8 125.7 61.2 58.2 ...
$ PO4 : num 170 558.8 187.1 138.7 97.6 ...
$ Chla : num 50 1.3 15.6 1.4 10.5 ...
$a4
Task Name :: a4
Formula :: a4 ~ .
<environment: 0x0000000010364358>
Task Data ::
'data.frame': 198 obs. of 12 variables:
$ a4 : num 0 1.9 0 0 0 0 3.9 0 0 2.9 ...
$ season: Factor w/ 4 levels "autumn","spring",..: 4 2 1 2 1 4 3 1 4 4 ...
$ size : Factor w/ 3 levels "large","medium",..: 3 3 3 3 3 3 3 3 3 3 ...
$ speed : Factor w/ 3 levels "high","low","medium": 3 3 3 3 3 1 1 1 3 1 ...
$ mxPH : num 8 8.35 8.1 8.07 8.06 8.25 8.15 8.05 8.7 7.93 ...
$ mnO2 : num 9.8 8 11.4 4.8 9 13.1 10.3 10.6 3.4 9.9 ...
$ Cl : num 60.8 57.8 40 77.4 55.4 ...
$ NO3 : num 6.24 1.29 5.33 2.3 10.42 ...
$ NH4 : num 578 370 346.7 98.2 233.7 ...
$ oPO4 : num 105 428.8 125.7 61.2 58.2 ...
$ PO4 : num 170 558.8 187.1 138.7 97.6 ...
$ Chla : num 50 1.3 15.6 1.4 10.5 ...
$a5
Task Name :: a5
Formula :: a5 ~ .
<environment: 0x0000000008812940>
Task Data ::
'data.frame': 198 obs. of 12 variables:
$ a5 : num 34.2 6.7 0 1.4 7.5 22.5 5.8 5.5 0 0 ...
$ season: Factor w/ 4 levels "autumn","spring",..: 4 2 1 2 1 4 3 1 4 4 ...
$ size : Factor w/ 3 levels "large","medium",..: 3 3 3 3 3 3 3 3 3 3 ...
$ speed : Factor w/ 3 levels "high","low","medium": 3 3 3 3 3 1 1 1 3 1 ...
$ mxPH : num 8 8.35 8.1 8.07 8.06 8.25 8.15 8.05 8.7 7.93 ...
$ mnO2 : num 9.8 8 11.4 4.8 9 13.1 10.3 10.6 3.4 9.9 ...
$ Cl : num 60.8 57.8 40 77.4 55.4 ...
$ NO3 : num 6.24 1.29 5.33 2.3 10.42 ...
$ NH4 : num 578 370 346.7 98.2 233.7 ...
$ oPO4 : num 105 428.8 125.7 61.2 58.2 ...
$ PO4 : num 170 558.8 187.1 138.7 97.6 ...
$ Chla : num 50 1.3 15.6 1.4 10.5 ...
$a6
Task Name :: a6
Formula :: a6 ~ .
<environment: 0x0000000008277740>
Task Data ::
'data.frame': 198 obs. of 12 variables:
$ a6 : num 8.3 0 0 0 4.1 12.6 6.8 8.7 0 0 ...
$ season: Factor w/ 4 levels "autumn","spring",..: 4 2 1 2 1 4 3 1 4 4 ...
$ size : Factor w/ 3 levels "large","medium",..: 3 3 3 3 3 3 3 3 3 3 ...
$ speed : Factor w/ 3 levels "high","low","medium": 3 3 3 3 3 1 1 1 3 1 ...
$ mxPH : num 8 8.35 8.1 8.07 8.06 8.25 8.15 8.05 8.7 7.93 ...
$ mnO2 : num 9.8 8 11.4 4.8 9 13.1 10.3 10.6 3.4 9.9 ...
$ Cl : num 60.8 57.8 40 77.4 55.4 ...
$ NO3 : num 6.24 1.29 5.33 2.3 10.42 ...
$ NH4 : num 578 370 346.7 98.2 233.7 ...
$ oPO4 : num 105 428.8 125.7 61.2 58.2 ...
$ PO4 : num 170 558.8 187.1 138.7 97.6 ...
$ Chla : num 50 1.3 15.6 1.4 10.5 ...
$a7
Task Name :: a7
Formula :: a7 ~ .
<environment: 0x000000000fe02028>
Task Data ::
'data.frame': 198 obs. of 12 variables:
$ a7 : num 0 2.1 9.7 1.4 1 2.9 0 0 0 1.7 ...
$ season: Factor w/ 4 levels "autumn","spring",..: 4 2 1 2 1 4 3 1 4 4 ...
$ size : Factor w/ 3 levels "large","medium",..: 3 3 3 3 3 3 3 3 3 3 ...
$ speed : Factor w/ 3 levels "high","low","medium": 3 3 3 3 3 1 1 1 3 1 ...
$ mxPH : num 8 8.35 8.1 8.07 8.06 8.25 8.15 8.05 8.7 7.93 ...
$ mnO2 : num 9.8 8 11.4 4.8 9 13.1 10.3 10.6 3.4 9.9 ...
$ Cl : num 60.8 57.8 40 77.4 55.4 ...
$ NO3 : num 6.24 1.29 5.33 2.3 10.42 ...
$ NH4 : num 578 370 346.7 98.2 233.7 ...
$ oPO4 : num 105 428.8 125.7 61.2 58.2 ...
$ PO4 : num 170 558.8 187.1 138.7 97.6 ...
$ Chla : num 50 1.3 15.6 1.4 10.5 ...
================================
> res.all <- experimentalComparison(
+ DSs,
+ c(variants('cv.lm'),
+ variants('cv.rpart',se=c(0,0.5,1))
+), + cvSettings(5,10,1234))
>summary(res.all)
-----------------------------------------------------------
== Summary of a Cross Validation Experiment ==
5 x 10 - Fold Cross Validation run with seed = 1234
* Data sets :: a1, a2, a3, a4, a5, a6, a7
* Learners :: cv.lm.v1, cv.rpart.v1, cv.rpart.v2, cv.rpart.v3
* Summary of Experiment Results:
-> Datataset: a1
*Learner: cv.lm.v1
nmse
avg 0.7077282
std 0.1639373
min 0.4661104
max 1.2218455
invalid 0.0000000
*Learner: cv.rpart.v1
nmse
avg 0.6423100
std 0.2399321
min 0.2146359
max 1.1712674
invalid 0.0000000
*Learner: cv.rpart.v2
nmse
avg 0.6569726
std 0.2397636
min 0.2146359
max 1.3356744
invalid 0.0000000
*Learner: cv.rpart.v3
nmse
avg 0.6875212
std 0.2348946
min 0.2886238
max 1.3356744
invalid 0.0000000
-> Datataset: a2
*Learner: cv.lm.v1
nmse
avg 1.0449317
std 0.6276144
min 0.4098820
max 3.6764115
invalid 0.0000000
*Learner: cv.rpart.v1
nmse
avg 1.0426327
std 0.2005522
min 0.7584102
max 1.9152968
invalid 0.0000000
*Learner: cv.rpart.v2
nmse
avg 1.01626123
std 0.07435826
min 0.99287105
max 1.49933643
invalid 0.00000000
*Learner: cv.rpart.v3
nmse
avg 1.000000e+00
std 2.389599e-16
min 1.000000e+00
max 1.000000e+00
invalid 0.000000e+00
-> Datataset: a3
*Learner: cv.lm.v1
nmse
avg 1.0216439
std 0.3522588
min 0.6818271
max 2.3354516
invalid 0.0000000
*Learner: cv.rpart.v1
nmse
avg 1.0406844
std 0.1823151
min 0.9234394
max 2.1525733
invalid 0.0000000
*Learner: cv.rpart.v2
nmse
avg 1.000000e+00
std 1.980954e-16
min 1.000000e+00
max 1.000000e+00
invalid 0.000000e+00
*Learner: cv.rpart.v3
nmse
avg 1.000000e+00
std 1.980954e-16
min 1.000000e+00
max 1.000000e+00
invalid 0.000000e+00
-> Datataset: a4
*Learner: cv.lm.v1
nmse
avg 2.1119756
std 3.1181959
min 0.3990221
max 15.4079825
invalid 0.0000000
*Learner: cv.rpart.v1
nmse
avg 1.0073953
std 0.1065607
min 0.4651041
max 1.4588424
invalid 0.0000000
*Learner: cv.rpart.v2
nmse
avg 1.000000e+00
std 2.774424e-16
min 1.000000e+00
max 1.000000e+00
invalid 0.000000e+00
*Learner: cv.rpart.v3
nmse
avg 1.000000e+00
std 2.774424e-16
min 1.000000e+00
max 1.000000e+00
invalid 0.000000e+00
-> Datataset: a5
*Learner: cv.lm.v1
nmse
avg 0.9316803
std 0.3478931
min 0.4980249
max 2.2706313
invalid 0.0000000
*Learner: cv.rpart.v1
nmse
avg 1.1793298
std 0.5659329
min 0.6366705
max 3.5927541
invalid 0.0000000
*Learner: cv.rpart.v2
nmse
avg 1.01137139
std 0.07049532
min 1.00000000
max 1.49445097
invalid 0.00000000
*Learner: cv.rpart.v3
nmse
avg 1.000000e+00
std 2.841612e-16
min 1.000000e+00
max 1.000000e+00
invalid 0.000000e+00
-> Datataset: a6
*Learner: cv.lm.v1
nmse
avg 0.9359697
std 0.6045963
min 0.1211761
max 3.9247377
invalid 0.0000000
*Learner: cv.rpart.v1
nmse
avg 1.0191041
std 0.1991436
min 0.8225254
max 2.3774810
invalid 0.0000000
*Learner: cv.rpart.v2
nmse
avg 1.000000e+00
std 2.451947e-16
min 1.000000e+00
max 1.000000e+00
invalid 0.000000e+00
*Learner: cv.rpart.v3
nmse
avg 1.000000e+00
std 2.451947e-16
min 1.000000e+00
max 1.000000e+00
invalid 0.000000e+00
-> Datataset: a7
*Learner: cv.lm.v1
nmse
avg 1.5238076
std 1.3762973
min 0.8670194
max 8.2205224
invalid 0.0000000
*Learner: cv.rpart.v1
nmse
avg 1.889286
std 3.528353
min 1.000000
max 20.829858
invalid 0.000000
*Learner: cv.rpart.v2
nmse
avg 1.463952
std 2.229237
min 1.000000
max 15.763627
invalid 0.000000
*Learner: cv.rpart.v3
nmse
avg 1.0295055
std 0.2086351
min 1.0000000
max 2.4752731
invalid 0.0000000
================================
其中DSs是7個資料集(a1～a7)，用來放入experimentalComparison()中逐一執行。sapply()第一及第三個參數，傳入第二參數（function(x,names.attrs)）作為參數。而且第一參數如果是向量，則每一個向量中的每一元素皆會逐一傳入。
paste(s1,s2)是將s1及s2字串連接成新字串。
as.formula(s)是將s字串轉換成R的公式(formula)物件。
同樣也可以圖形顯示：
================================
> plot(res.all)
================================
 |
| 7種測試之NMSE |
從上圖可以發現，大部分結果不是很好，因為NMSE明顯>1。{問題：是指平均值嗎？}
我們也可以使用bestScores()找出最佳的模型：
================================
> bestScores(res.all)
----------------------------------------------------------
$a1
system score
nmse cv.rpart.v1 0.64231
$a2
system score
nmse cv.rpart.v3 1
$a3
system score
nmse cv.rpart.v2 1
$a4
system score
nmse cv.rpart.v2 1
$a5
system score
nmse cv.lm.v1 0.9316803
$a6
system score
nmse cv.lm.v1 0.9359697
$a7
system score
nmse cv.rpart.v3 1.029505
================================
除了a1以外，其結果是令人失望的。從圖『7種測試之NMSE』顯示結果的變動性，可以提供一個觀點，使用『綜效(ensemble)』試驗方法可能是一個好的替代方案。綜效是一種模型建構方法，主要是解決個別模型的限制，作法是建立一大堆各式模型，然後綜和這些模型的預測。有許多途徑可以做到綜效，除了多樣的模型獲得方式（無不同的訓練樣本、不同的變數、不同的模型技術），也可以不同預測方式（如投票、平均等）。『隨機森林(Random forests)』是由一大群樹狀模型（迴歸樹或分類樹）所組成，每棵樹都是完整的（未經修剪），而且每一步成長節點的最佳切割是從隨機屬性子集中選擇。最終結果則是平均所有樹的結果。randomForest套件randomForest()即實作這個概念。下列程式重複前述交叉驗證試驗，其中包含第三種隨機森林：
================================
> library(randomForest)
> cv.rf <- function(form,train,test,...) {
+ m <- randomForest(form,train,...)
+ p <- predict(m,test)
+ mse <- mean((p-resp(form,test))^2)
+ c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2))
+}
> res.all <- experimentalComparison(
+ DSs,
+ c(variants('cv.lm'),
+ variants('cv.rpart',se=c(0,0.5,1)),
+ variants('cv.rf',ntree=c(200,500,700))
+), + cvSettings(5,10,1234))
================================
執行過程會花一點時間，最後使用bestScores()檢視執行結果：
================================
> bestScores(res.all)
----------------------------------------------------------
$a1
system score
nmse cv.rf.v3 0.5447361
$a2
system score
nmse cv.rf.v3 0.7777851
$a3
system score
nmse cv.rf.v2 0.9946093
$a4
system score
nmse cv.rf.v3 0.9591182
$a5
system score
nmse cv.rf.v1 0.7907947
$a6
system score
nmse cv.rf.v3 0.9126477
$a7
system score
nmse cv.rpart.v3 1.029505
================================
結果顯示，最佳的分數來自隨機森林模型，同樣的結果並不總是好的，尤其是a7。bestScores()顯示的資料，並未告訴我們這些最佳模型與其他模型之差異是否是統計上的顯著，亦即使用其他隨機資料可以得到類似結果的信賴度。使用compAnalysis()便可以得到這樣的資訊。cv.rf.v3在a1、a2、a4、a6表現最佳，下列驗證統計的顯著程度：
================================
> compAnalysis(res.all,against='cv.rf.v3',datasets=c('a1','a2','a4','a6'))
-----------------------------------------------------------
================================
其中sign.x提供了我們想要的資訊，無標示即表示cv.rf.v3於其他的模型差異的統計性顯著度低於95％，加號表示模型的平均評估指標顯著高於cv.rf.v3，而NMSE應該是愈低愈好。負號代表相對的意義。
你可以發現，不同隨機森林之間的差異，及與其他模型之間的差異，並不是統計上顯著的。相對其他模型，在大多數情況下隨機森林擁有較明顯的優勢。
*******************************************
問題：與其他模型有顯著差異，代表何種意義？
*******************************************
================================
