預測分析：R語言實作2.4 評估線性回歸模型

<b>2.4　評估線性回歸模型</b>

<b></b>

再次利用lm()函數，用線性回歸模型來拟合資料。我們的兩套資料集會用到上述資料框裡剩下的所有輸入特征。r提供了一種編寫公式的簡寫方式，它可以把某個資料框裡的所有列作為特征，除了被選為輸出的列之外。這是利用一個句号符來完成的，如下列代碼片段所示：

一旦我們準備好了所有的資料，訓練一個線性回歸模型就是一行代碼的事情，但是如果要研究這個模型以便判定其效果，後面就會有重要的工作。幸運的是，我們可以利用summary()函數立刻獲得某些關于這個模型的重要資訊。該函數對于cpu資料集的輸出如下所示：

在重複了一次之前對lm()函數本身的調用之後，summary()函數給出的資訊形成了三個獨立的部分。第一部分是模型殘差（residual）的摘要，它是模型對訓練它的資料中的觀測資料産生的誤差。第二部分是一個表格，包含了模型各系數的預測值及其顯著性檢驗的結果。最後幾行則顯示了該模型的一些總體性能衡量名額。如果在二手車資料集上重複同樣的過程，我們會在該模型的摘要裡注意到如下的一行：

産生這個結果的原因是由于潛在的依賴關系，導緻我們還有1個特征對輸出的作用無法和其他特征厘清楚。這種現象被稱為混疊（aliasing）。alias()指令會顯示我們需要從模型中去除的特征：

正如我們所見，有問題的特征就是saturn特征，是以我們會去除這個特征并重新訓練模型。要從一個線性回歸模型去除某個特征，可以把它加到公式中的句點之後，并在它之前加上一個減号：

2.4.1　殘差分析

殘差就是模型對特定觀測資料産生的誤差。換言之，它是輸出的實際值和預測值之間的差異：

ei＝yi＝

在建構一個良好的回歸模型時，對殘差的分析是非常重要的，因為它能展現這個模型的很多方面，從不能成立的假設和拟合的品質到其他問題（比如離群值）。要了解殘差摘要裡的衡量名額，可以想象把殘差從小到大排序。除了在這個序列裡兩個極端位置上的最小值和最大值之外，這個摘要還會顯示第一個和第三個四分位數（quartile），它們分别是這個序列中處于四分之一和四分之三位置的值。中位數（median）則是序列正中位置的值。四分位間距（interquartile range）是序列處于第一和第三個四分位數之間的部分，根據定義可知，它包含了一半資料。首先來看一下cpu模型裡的殘差摘要，其中有趣的一點是，第一個和第三個四分位數的值與最小和最大值相比是很小的。這是存在一些具有很大殘差誤差的點的首要迹象。在理想情況下，殘差應該具有0中位數和較小的四分位數值。可以通過關注lm()函數所産生模型具有的residuals屬性，來重制summary()函數産生的殘差摘要：

注意，在前面關于二手車模型的示例中，我們需要把殘差的值和輸出變量的平均值進行比較，以獲得殘差是否過大的印象。這樣，前面的結果顯示在訓練資料中二手車的平均售價為$21k，預測結果的 50% 基本是在正确值的±$1.6k範圍内，這個情況看起來是相當合理的。顯然，cpu模型的殘差的絕對值都小得多，這是因為該模型的輸出變量（即釋出的cpu相對性能）比二手車模型裡的price變量值要小得多。

線上性回歸裡，我們假設模型裡不能化簡的誤差是以正态分布的形式随機分布的。一種稱為分位圖（quantile-quantile plot，q-q plot）的診斷圖有助于從視覺上評判該假設的符合程度。這種圖背後的關鍵思想是，我們可以通過比較它們分位數（quantile）的值來比較兩種分布。一個分布的分位數基本上是一個随機變量的距離均等的區間，這樣的每個區間具有相同的機率。例如，四分位數把一個分布劃分為4個機率相等的部分。如果兩個分布相同，那麼分位圖就應該是直線y＝x的圖形。要檢查殘差是否符合正态分布，可以用它們的分布和一個正态分布進行比較，檢視得到的圖形和y＝x的接近程度。

有很多其他的方法可以檢查模型殘差是否是正态分布。有一個很好的工具是r語言的nortest包，它實作了很多知名的正态性檢驗方法，包括anderson-darling檢驗和lilliefors檢驗。此外，stats包包含了shapiro.test()函數，可以用于進行shapiro-wilk正态性檢驗。

下列代碼會對兩個資料集建立分位圖：

兩個資料集的分位圖如下所示：

兩個模型的殘差看起來和理論上的正态分布分位數相當貼近，雖然拟合程度并非完美，但是這對于大部分實際環境資料來說是典型的情況。第二種對線性回歸非常有用的診斷圖被稱為殘差圖（residual plot）。它是反映訓練資料中觀測資料的殘差和對應拟合值之間關系的一種圖表。換言之，它是配對（, ei）的一個圖。殘差圖有兩個重要屬性是我們特别感興趣的。首先，我們希望通過檢查殘差是否在不同的拟合值範圍中保持大小一緻，來确認有關常數方差的假設是成立的。其次，我們必須檢驗在殘差中不存在某種模式。但是，如果我們觀察到了某種模式，它就有可能是一個迹象，表明對應的模型對于涉及的特征是非線性的，或者在模型裡有缺失的特征沒有被納入。實際上，一種發現新特征的有用方法就是尋找和模型殘差相關的新特征。

這兩張圖在圖的左側部分都顯示了細微的殘差遞減模式。稍微更令人擔心的是，對于兩個輸出變量的較大值，殘差的方差（the variance of the residual）看起來有點大，這有可能表明誤差并不是同方差的（homoscedastic）。這個問題在第二個有關二手車資訊的圖中更為顯著。在前面的兩個殘差圖中，我們也标記了某些較大的殘差（絕對值大小）。我們很快還會看到可能有一些潛在的離群值。另一種獲得殘差圖的方法是對lm()函數本身産生的模型調用plot()函數。這樣會建立4個診斷圖，包括殘差圖和分位圖。

2.4.2　線性回歸的顯著性檢驗

在檢查了殘差的摘要之後，我們應該關注的下一個東西是模型所産生的系數表。在這裡，每個估算出來的系數都配套了額外的一組數值，在最後還有一組星号或點。首先，這樣一連串的數值看起來可能有點混亂，但是把所有這些資訊包括進來是有充分理由的。當對某些資料收集測量值并指定一組特征來建構一個線性回歸模型時經常會遇到的情況是，裡面有一個或多個特征和我們要預測的輸出實際上并不相關。當然，這是我們在之前收集資料的時候并沒有在意的問題。理想情況下，我們會希望這個模型不但找到和輸出變量确實存在依賴關系的那些特征所對應的最佳系數值，而且能告訴我們哪些特征是不需要的。

判定模型是否需要某個特征的一種可參考的方法是訓練兩個模型而不是一個。第二個模型會包含第一個模型的所有特征，再去掉需要确認其顯著性的那個特征。然後就可以通過它們殘差的分布情況來檢驗兩個模型是否有差别。實際上這就是r語言會對我們在每個模型裡指定的所有特征進行的處理。對于每個系數，都會針對它對應的特征和輸出變量無關的無效假設（null hypothesis）建構一個置信區間（confidence interval）。具體而言，對于每個系數，我們都會考慮一個包括除了該系數對應的特征之外的所有其他特征的線性模型。然後，我們會檢驗把這個特征加入模型是否能顯著改變殘差誤差的分布，這種改變會是該特征和輸出存線上性關系的證據，如果成立，則該系數就應該是一個非0值。r語言的lm()函數會自動為我們進行這些檢驗。

在統計學裡，置信區間會把一個點估計（point estimate）和該估計的精度合并在一起。具體做法是指定一個區間，在某種程度的可信度下，可以期望要估計的參數的實際值會處于該區間内。某個參數的95%置信區間實際上告訴我們的是，如果要從同一個實驗收集 100 個資料樣本，并在每個樣本裡為這個估計的參數建構一個95%的置信區間，那麼這些資料樣本裡有 95 個會滿足目标參數的實際值，處于它所對應的置信區間内。對于同樣程度的可信度，為方差較大的點估計（例如對隻有很少的資料點作出估計的情況）建構的置信區間，和對較小方差作出的估計相比，往往會定義出範圍更大的區間。

我們來看看cpu模型的摘要輸出的一個快照，它顯示了cpu模型中截距和 myct 特征的系數：

我們目前關注的是myct特征，它所在行裡的第一個數值是其系數的估算值，這個值大約是0.05 (5.210×10-2)。标準誤（standard error）是這個估算值的标準差（standard deviation），它等于下一個數值給出的0.01885。我們可以通過對在0和該系數估算值之間的标準誤個數進行計數，判定該系數是否實際為0（代表該特征與輸出線性無關）的置信度。為此，我們可以用系數估算值除以标準誤，這就是t值（t-value）的定義，也就是該行裡的第三個值：

這樣，我們的myct系數值和0值之間幾乎隔開了3個标準誤的數值，這個名額很好地表明該系數不太可能為0。t值越大，就越能說明我們應該給模型中的這個特征賦予非0的系數值。我們可以把這個絕對數值轉變為一個機率，它可以告訴我們該系數的真實值為0的可能性。這個機率是從學生t分布（student抯 t-distribution）得到的，被稱為p值（p-value）。對于myct特征，這個機率是0.006335，是很小的。可以利用pt()函數得到這個值：

pt()函數是t分布的分布函數，而t分布是對稱的。為了了解p值以上面這種方式進行計算的原因，要注意我們感興趣的是t值的絕對值比我們計算的值更大的機率。要得到這個結果，我們首先要獲得該t分布上尾或右尾（the upper / right tail）部分的機率，然後将其乘以2，這樣就把下尾（the lower tail）部分也包含進來了。運用基本的分布函數在r語言裡是一項非常重要的技能，如果前面這個示例看起來難度太大，我們在關于r語言的線上教程裡還有一些示例可供參考。t分布是以自由度作為參數的。

自由度的數量實際上是在計算某個統計量（例如某個系數的估計值）時能夠自由變化的變量數。線上性回歸主題下，它就相當于訓練資料中觀測資料的個數減去模型中參數的個數（即回歸系數的個數）。對于我們的cpu模型，這個數字就是179－7＝172。對于具有更多資料點的二手車模型，這個數字是664。自由度這個名字來源于它和獨立次元（或者說，作為某個系統的輸入的資訊片段）的數量之間的關系，進而反映系統能夠（在不違反任何對于輸入的限制條件的條件下）自由配置的程度。

作為一般的經驗規則，我們會希望p值小于0.05，這相當于我們希望系數不為0值的估計具有95%的置信區間。每個系數旁邊的星号數量給我們提供了一種關于置信度水準的快速視覺輔助，其中一顆星對應95%的經驗規則，而兩顆星代表了99%的置信區間。是以，在模型摘要裡，每個不帶任何星号的系數對應的就是一個按照經驗規則我們沒有把握将其納入模型的特征。在cpu模型裡，chmin特征是唯一可疑的特征，而其他特征的p值都非常小。二手車模型的情況就不同了。這裡，有4個特征都是可疑的，截距也是如此。

線上性回歸模型主題下正确了解p值的含義是很重要的。首先，我們不能也不應該對p值進行互相比較來評判哪個特征是最重要的。其次，較高的p值并不一定代表某個特征和輸出之間不存線上性關系；它隻表明在所有其他模型特征存在時，這個特征對于輸出變量不能提供任何新的資訊。最後，我們必須永遠記住，95%的經驗規則并非絕對真理，它隻是在特征（以及對應的系數）的數量不是特别大的時候确實有用。在95%的置信度下，如果在模型裡有1000個特征，可以預期得到錯誤系數的結果是50個。是以，線性回歸系數的顯著性檢驗對于高維問題并不是那麼有用。

最後一個顯著性檢驗值實際上出現在lm()輸出摘要的最下部，并且是在最後一行。這一行給我們提供了f統計量（f statistic），其命名來源于檢查兩個（理想狀況下是正态的）分布的方差之間是否存在統計顯著性（statistical significance）的f檢驗（f test）。這個示例中的f統計量嘗試評價某個所有系數都為0的模型中的殘差方差（the variance of the residuals）是否和我們訓練的模型的殘差方差有顯著性的差異。

換言之，f檢驗會告訴我們訓練的模型是否能解釋輸出中的某些方差，如果能，我們就知道至少有一個系數必須為非0值。雖然這個檢驗對于我們有很多系數的情況不是那麼有用，但它還是能檢驗系數的整體顯著性，而不會遇到t檢驗對于單個系數的那種問題。摘要部分對此顯示了一個非常小的p值，是以我們知道我們的系數裡至少有一個是非0值。我們可以調用anova()函數重制運作過的f檢驗，該函數代表方差分析（analysis of variance）。該檢驗會用隻有截距而沒有特征的方式建立的空模型（null model）和訓練模型進行比較。我們要用cpu資料集對這個函數進行如下示範：

注意，空模型的公式是prp～1，其中的1代表截距。

2.4.3　線性回歸的性能衡量名額

在摘要裡，最後的細節是和該模型作為一個整體的性能以及該線性模型和資料的拟合程度相關的。要了解我們如何評估線性回歸的拟合，必須首先指出的是，線性回歸模型的訓練準則是使均方差（mse）在訓練資料上最小化。換言之，把一個線性模型拟合到一組資料點相當于找到一條直線，它的斜率和位置要讓它與這些點的距離的平方的總和（或平均值）最小。因為我們把某個資料點和它在該直線上的預測值之間的誤差稱為殘差，我們可以定義殘差平方和（residual sum of square，rss）為所有殘差的平方之和：

換言之，rss就是誤差平方和（sum of squared error，sse），是以我們可以通過下面的簡單公式把它和我們已經熟悉的均方差（mse）聯系起來：

除了某些曆史原因之外，rss之是以成為值得注意的重要衡量名額，是因為它和另一個被稱為rse的重要衡量名額相關，我們後續會讨論到rse。為此，我們需要先對訓練線性回歸模型的過程有個直覺的了解。如果用人工資料多次進行簡單線性回歸實驗，每次改變随機種子（random seed）以獲得不同的随機樣本，就會得到一組很可能非常接近實際總體回歸線的回歸線，正如進行單次實驗所看到的那樣。這表明了一個事實，即線性模型總體而言是具有低方差的特點的。當然，嘗試近似的未知函數也有可能是非線性的，如果是這樣，那麼即使總體回歸線也不太可能對于非線性函數的資料産生良好的拟合。這是因為線性假設是非常嚴格的，是以，線性回歸是一種帶有較高偏誤的方法。

我們要定義一個叫作殘差标準差（residual standard error，rse）的衡量名額，它會估計我們的模型和目标函數之間的标準差。這就是說，它會大緻地衡量模型和總體回歸線之間的平均距離。它是以輸出變量的機關衡量的，并且是一個絕對值。是以，它需要和y的值進行比較，以便判定它對于特定樣本是高還是低。帶有k個輸入特征的模型的一般rse可以按如下公式計算：

對于簡單線性回歸，其特征數k＝1：

可以利用前面的公式對兩個模型計算rse值，結果如下所示：

要解釋這兩個模型的rse值，就需要把它們和輸出變量的均值進行比較：

注意，在cpu模型中，其取值為61.3的rse值和二手車模型大約為2947的rse值相比是相當小的。但是，就它們和對應輸出變量的均值的近似程度而言，我們可以看到反而是二手車模型的rse表現了更優的拟合。

現在，雖然rse作為一個絕對數值在和輸出變量的均值進行比較方面是有用的，但是我們經常需要一個相對數值用來在不同的訓練場景下進行比較。為此，在評估線性回歸模型的拟合情況時，我們也經常觀察r2統計量（r2 statistic）。在摘要資訊裡，它是用“multiple r-squared”來訓示的。在提供它的計算公式之前，我們首先要講解總平方和（total sum of square，tss）的概念。總平方和是和輸出變量裡的總方差成比例的，它的設計思路是，在進行回歸之前，用它來衡量輸出變量内在的變異性大小。tss的計算公式是：

r2統計量背後的思想是，如果一個線性回歸模型是對實際總體模型比較接近的拟合，它就應該能完整地獲得輸出中的所有方差。實際上，我們經常把r2統計量當做一個相對數值，用它來指明某個回歸模型解釋了輸出變量的方差中多大的比例。當運用回歸模型得到了輸出變量的一個估計值時，可以看到，觀測資料中的誤差被稱為殘差，而rss實際上是與預測值和輸出函數真實值之間留下的差異成比例的。是以，可以把r2統計量（即線性回歸模型所解釋的輸出y中方差的數量）定義為起始方差（tss）和終結方差（rss）之間的內插補點相對于這個起始方差（tss）的比例。用公式表示，這個計算就是：

從這個公式我們可以看到，r2的範圍在0和1之間。接近1的值表明良好的拟合，因為它意味着大部分輸出變量裡的方差已經被回歸模型所解釋。另一方面，一個較小的值表明模型裡的誤差中還有顯著的方差，說明模型拟合得不夠好。讓我們來看看，對于這兩個模型，r2統計量是如何手工計算出來的：

我們用到了lm()函數訓練的模型的fitted.values屬性，它是模型對于訓練資料作出的預測。計算出的兩個值都相當大，其中二手車模型又一次顯示出了略微更好一些的拟合。現在我們看到了評價線性回歸模型的兩個重要衡量名額，即rse和r2統計量。現在，我們要考慮兩個變量之間的線性關系是否還有一種更總體性的衡量名額，也可以運用到示例中。從統計學角度，我們可以想起相關系數的概念，它正好是描述這個問題的。

兩個随機變量x和y之間的相關系數（correlation）可以由下面的公式得出：

其實，在簡單回歸的情況下，輸出變量和輸入特征之間相關系數的平方就是r2統計量，這個結果進一步增強了後者作為一種有用的衡量名額的重要性。

2.4.4　比較不同的回歸模型

當需要對用同一組輸入特征訓練的兩個不同回歸模型之間進行比較時，r2統計量就會非常有用。不過，我們還經常需要對輸入特征數量并不相等的兩個模型進行比較。例如，在特征選擇的過程中，我們會需要知道把某個特征納入模型是否是個好的思路。r2統計量的一個局限性是它對具有更多輸入參數的模型往往會産生更大的值。

調整後的r2（adjusted r2）試圖糾正一個問題，即r2統計量對具有更多輸入特征的模型總是會産生更大的值，因而容易受到過拟合的影響。調整後的r2總體上比r2本身更小，這是我們可以通過檢查模型摘要裡的值來檢驗的。調整後的r2的計算公式如下所示：

其中n和k的定義和在r2統計量裡是一樣的。現在，讓我們在r語言裡實作這個函數，并對我們的兩個模型計算調整後的r2：

另外還有幾個常用的性能衡量名額也可以用來比較帶有不同數量特征的模型。赤池資訊準則（akaike information criterion，aic）使用了一種資訊論方法，通過平衡模型的複雜性和精确度來評估模型的相對品質。在我們按方差最小化方法訓練的線性回歸模型上，它和另一個著名統計量—馬洛斯cp（mallow抯 cp）是成比例的，是以這兩個性能衡量名額是可以互換使用的。第三個衡量名額是貝葉斯資訊準則（bayesian information criterion，bic）。它和前兩個衡量名額相比，對帶有更多變量的模型往往會産生更不利的結果。

2.4.5　在測試集上的性能

到目前為止，我們在訓練資料上觀察了模型的性能。這對于判定一個線性模型能否很好地拟合資料有重要意義，但不會讓我們對它處理未知資料的預測精确度有清晰的認識。為此，我們要轉到測試資料集。要利用模型進行預測，可以調用predict()函數。這是r語言裡的一個通用函數，有很多包擴充了它。對于用lm()訓練過的模型，我們隻需要提供模型和帶有我們需要進行預測的觀測資料的資料框即可：

下一步，要定義我們自己的函數，用來計算均方差（mse）：

對于每個模型，我們已經調用了compute_mse()函數來傳回訓練和測試的均方差。在這個示例中，碰巧兩個測試均方差的值都比訓練均方差的值更小。測試均方差是略微大于還是略微小于訓練均方差并不是特别重要的問題。重要的問題是測試均方差并沒有顯著大于訓練均方差，因為測試均方差顯著大于訓練均方差就表明我們的模型對資料是存在過拟合的。注意，尤其是對于cpu模型，在原始資料集裡的觀測資料數量很少，這就導緻測試集也是很小的。是以，就上述模型對于未知資料的預測性能而言，我們會在估計的精确度方面持保留态度，因為利用一個小的測試集作出的預測會具有更大的方差。