【譯文】R語言線性回歸入門

作者 Troy Walters

譯者 錢亦欣

基於最小二乘法的線性回歸是你在學習數據科學和機器學習時最先遇到的模型之一，它不僅簡單易懂，還能在很多問題中發揮作用，並且已經集成在了很多種編程語言之中。大部分用戶對R語言中的lm()函數肯定不陌生，它讓你能簡易而快速地擬合一個線性回歸模型。然而，這個函數並不現實參數估計和很多檢驗統計量的計算過程，所以本文就打算手把手地計算實現線性回歸模型的參數估計。

本文中，我只會使用矩陣、向量和矩陣操作符來獲得參數估計（基於最小二乘法的線性回歸模型的參數估計有解析表達式）。在讀完本文之後，你會發現自己動手也很簡單，並且你可以把這個過程套用在任何數據集之上（儘管使用lm()函數會更快更穩健）。

我將使用MASS包中大名鼎鼎的Boston數據集作為例子，它包含了14個關於經濟，地理和人口方面的變數，一共有506個觀測值。每個觀測代表1990年代波士頓某個地區的特徵。我們的被解釋變數將會是每個地區的房價的中位數，對應列名是"medv"，讓我們先看看數據集長什麼樣：

library(MASS)str(Boston)## "data.frame": 506 obs. of 14 variables:## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...## $ rm : num 6.58 6.42 7.18 7 7.15 ...## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...## $ black : num 397 397 393 395 397 ...## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

你可以運行help(Boston)來看看每個變數的定義。

線性回歸的表達式一般如下：

其中，y是響應變數（被解釋變數），X是設計矩陣（特徵變數），是需要估計的參數向量，是殘差項，它反映了哪些沒被納入模型的因素對y的影響。

萬事俱備，第一步，我們要聲明響應變數

y <- Boston$medv

然後，我們需要確定輸入特徵的設計矩陣，也就是公式中的X。我們可以使用as.matrix()函數完成這一操作，吧『medv』以外的所有變數納入。我們還需要在X中添加一列元素全為1的列來估計截距項。

# 設計矩陣X <- as.matrix(Boston[-ncol(Boston)])# 生成用來估計截距項的列int <- rep(1, length(y))# 添加用來估計截距項的列X <- cbind(int, X)

現在被解釋變數和設計矩陣都有了，我們可以通過如下的計算式來得到係數的估計值：

這個表達式的推導超出了本文的範疇，如果你有興趣，可以參看Wikipedia 或者任何一本相關的教科書。

我們可以使用剛剛賦值好的X和y，結合%*%操作符來實現矩陣乘法。t()函數能將矩陣轉置，solve()函數能夠得到任意可逆矩陣的逆矩陣。

# 實現解析表達式betas <- solve(t(X) %*% X) %*% t(X) %*% y# 結果保留兩位小數，更美觀betas <- round(betas, 2)

現在我們有係數向量了，它代表了y和解釋變數之間的線性關係。比如，房間數每增長1單位，房價就增長3810刀。我們可以遵循這個方式來解釋估計結果，但需要注意每個變數的量綱。注意，變數『chas』是個只取0，1兩個值的虛擬變數。1代表該區毗鄰查理斯河，0代表相反情況。而對應的係數則表明，河邊的房子的房價中位數比內陸的搞2690刀（在其他變數一定時）。

譯者註：作者沒有考慮變數的顯著性就直接解釋結果了，這是不對的。如果變數的係數並不顯著，我們可以在統計意義上認為這個係數是0，也就是該變數和y之間沒有明顯的相關性。

print(betas)## [,1]## int 36.46## crim -0.11## zn 0.05## indus 0.02## chas 2.69## nox -17.77## rm 3.81## age 0.00## dis -1.48## rad 0.31## tax -0.01## ptratio -0.95## black 0.01## lstat -0.52

現在，讓我們來看看我們得到的結果是否正確，調用lm()建模並比較結果，我知道你們很想知道這個。

# 線性回歸模型lm.mod <- lm(medv ~ ., data=Boston)# 保留兩位小數lm.betas <- round(lm.mod$coefficients, 2)# 結果放到一個數據框內results <- data.frame(our.results=betas, lm.results=lm.betas)print(results)## our.results lm.results## int 36.46 36.46## crim -0.11 -0.11## zn 0.05 0.05## indus 0.02 0.02## chas 2.69 2.69## nox -17.77 -17.77## rm 3.81 3.81## age 0.00 0.00## dis -1.48 -1.48## rad 0.31 0.31## tax -0.01 -0.01## ptratio -0.95 -0.95## black 0.01 0.01## lstat -0.52 -0.52

係數估計的結果一模一樣，現在你知道怎麼手把手得到係數估計值了。但是使用lm()函數無疑更快，並且它還提供了諸如擬合優度，t值，p值等檢驗統計量，還是用QR分解等技術使得計算更加穩健。可無論如何，我希望展示最為原始的解法。

線性回歸的係數估計還可以通過其他方法得到，比如在特徵數量很大的時候使用梯度下降法就是不錯的選擇。這個時候再利用解析表達式結算可能就很費時間。而隨機梯度下降則在特徵數和樣本數都很大的時候適用，我會在後續文章中介紹它們。

你可以登陸我的github獲取所有代碼.

祝好！

註：原文刊載於Datascience+網站

原文鏈接：Linear Regression from Scratch in R

推薦閱讀：