式中β1-βn代表了預測變量X1-Xn的回歸係數（斜率），β0是線性方程的截距。回歸的目的是找到一組最優的模型參數（斜率和截距值的組合），使響應變量的觀測值和預測值之間的殘差平方和最小化。
從技術上來說，簡單多項式回歸也可以視為多元線性回歸的特例，例如二次回歸有兩個預測變量（X1=X和X2=X2），三次回歸有三個預測變量（X1=X、X2=X2和X3=X3）。
本篇繼續通過一個示例，展示在R語言中擬合多元線性回歸的方法過程。
下文的示例數據、R代碼等，可見網盤附件（提取碼，54ua）：
https://pan.baidu.com/s/1ixB_7VQ5dXHu7ckyv1r-GQ
節選自馬裡蘭州生物資源調查（https://dnr.maryland.gov/streams/Pages/mbss.aspx）。
就節選部分數據為例，記錄了所調查的馬裡蘭州河流中每75米長的區段水域內，魚類物種Rhinichthys cataractae的豐度，並測量了每段水域中相應的環境特徵。
其中第一列代表了調查河流區段的位置信息，其餘各列依次為：
fish，水域中魚類物種Rhinichthys cataractae的個體數量，代表了物種豐度；
acre，水域流域面積（英畝，acre）；
do2，水域溶解氧含量（毫克/升，mg/L）；
depth，水域最大深度（釐米，cm）；
no3，水域硝酸鹽濃度（毫克/升，mg/L）；
so4，水域硫酸鹽濃度（毫克/升，mg/L）；
temp，水域溫度（攝氏度，℃）。
一個生物學目的是探索可能影響魚類物種Rhinichthys cataractae豐度的環境因素，並對物種豐度變化的原因作出假設。
線性關係無疑是最簡單的類型，儘管實際中變量間關係並非線性這樣簡單，但由於線性關係比較直觀且易於理解，大多數實際分析中我們更期望首選將問題簡化為線性模型去解釋。
因此，不妨首先站在單個預測變量與響應變量角度上，通過一元線性回歸初步觀察每個環境因素和物種豐度的線性關係。
#讀取魚類物種豐度和水體環境數據
dat <- read.delim('fish_data.txt', sep = '\t', row.names = 1)
#繪製二維散點圖觀測各環境變量與魚類物種豐度的關係
library(ggplot2)
dat_plot <- reshape2::melt(dat, id = 'fish')
p <- ggplot(dat_plot, aes(value, fish)) +
geom_point() +
facet_wrap(~variable, ncol = 3, scale = 'free') +
geom_smooth(method = 'lm')
p
#擬合各環境變量與魚類物種豐度的一元回歸，並提取各自 R2 和 p 值
env <- c('acre', 'do2', 'depth', 'no3', 'so4', 'temp')
R2_adj <- c()
p_value <- c()
for (i in env) {
fit_stat <- summary(lm(dat[['fish']]~dat[[i]])) #一元線性回歸
R2_adj <- c(R2_adj, fit_stat$adj.r.squared) #提取校正後 R2
p_value <- c(p_value, fit_stat$coefficients[2,4]) #提取顯著性 p 值
}
env_stat <- data.frame(row.names = env, R2_adj, p_value)
env_stat #數據框中存儲了各環境變量與魚類物種豐度的一元回歸的 R2 和 p 值
#在散點圖中添加各環境變量與魚類物種豐度的一元回歸的 R2 和 p 值作為標識
#註：該 R2 和 p 值僅為單個變量一元回歸的 R2 和 p 值，和下文即將提到的多元回歸的 R2 和 p 值存在區別
env_stat$fish <- max(dat$fish) * 0.8
for (i in env) env_stat[i,'value'] <- max(dat[[i]]) * 0.8 #這句和上一句，定義文字在圖中的展示坐標，這裡僅粗略設置下
env_stat$variable <- rownames(env_stat)
env_stat$label <- paste('R2.adj =', round(env_stat$R2_adj, 3), '\np =', round(env_stat$p_value, 3)) #文字標籤
env_stat <- env_stat[c('fish', 'variable', 'value', 'label')]
dat_plot$label <- NA
dat_plot <- rbind(dat_plot, env_stat) #和先前的數據框合併，便於作圖
p + geom_text(data = dat_plot, aes(label = label), size = 3)
通過分別查看單個環境變量和物種豐度的關係，可以看到acre（流域面積）、no3（硝酸鹽濃度）以及depth（水域深度）與物種豐度間的線性關係較明顯，其餘三個環境因素則很弱。
接下來擬合環境因素與物種豐度的多元線性回歸，以綜合看待環境因素的效應，為可能導致物種豐度變化的原因提供見解。
類似一元線性回歸，多元線性回歸在R中也是通過函數lm()執行。
#首先不妨使用全部環境變量擬合與魚類物種豐度的多元線性回歸
#自變量之間用「+」相連
fit1 <- lm(fish~acre+do2+depth+no3+so4+temp, data = dat)
summary(fit1) #展示擬合方程的簡單統計
#如果數據集中所有的自變量都使用到，且不考慮交互作用，則上式可以簡化如下
#使用「.」代替所有自變量
fit1 <- lm(fish~., data = dat)
summary(fit1) #展示擬合方程的簡單統計
整體而言，全模型p<0.001是顯著的，表明多元線性回歸中存在部分環境因素可以解釋物種豐度變化。回歸方程校正後的R2=0.25，算不上很高，這是可預期的，正如上所提到的，物種豐度和大部分環境因素間關係實際上並非線性這樣簡單。儘管如此，仍可以對結果作出一些解釋。
關注單個預測變量的p值，p值越小代表該預測變量（環境因素）和響應變量（物種豐度）間線性關係越明顯。例如，acre（流域面積）和no3（硝酸鹽濃度）的p值顯著，且回歸係數大於0，表明在其它環境因素保持不變時，它們的增加有助於物種豐度的增加；depth（水域深度）接近臨界顯著性，暗示水深的增加也有助於提升物種的豐度；temp（水域溫度）的p值不顯著，表明物種的豐度並非隨著溫度的增加而呈線性增長。
通過比較，可以看到環境因素與物種豐度各自的一元回歸和整體的多元線性回歸之間趨勢基本一致，儘管存在微小差異。
如果期望從中提取某些關鍵的信息，或者用作其它統計，這是一些參考用法。
#提取或統計重要數值，例如
coefficients(fit1) #各變量的斜率和截距項
#也可以 names(summary(fit)) 後查看主要的內容項，然後從中提取，例如
summary(fit1)$adj.r.squared #校正後 R2
此外，如期望作圖觀測變量間關係，由於多元線性回歸是多維結構，很難直接在三維空間中表現出來。一種可選替代方案是通過多個一元回歸，單獨觀測每個環境因素和物種豐度的關係，這個在上文已經展示了。但需注意的是，兩種方法還是存在一些差異的，這個在上文比較二者就看出來了。
#作圖，三維空間中最多展示 1 個響應變量和 2-3 個預測變量間的線性關係，例如
library(car)
scatter3d(fish~acre+no3, data = dat)
#要麼就考慮使用多個一元回歸散點圖替代，但是它和多元線性回歸是存在一些區別的
由於未觀察到物種豐度與do2（水域溶解氧含量）、so4（水域硫酸鹽濃度）和temp（水域溫度）呈線性響應，因此不妨：
（1）考慮轉化數據，並繼續觀察轉化後的變量間線性關係是否明顯；
（2）直接把它們從回歸方程中剔除。
轉化數據需要思考很多問題，懶得思考，就直接選擇把它們剔除了，來看新擬合的多元線性回歸方程。
#只考慮使用 3 個和魚類物種豐度線性關係較為明顯的環境變量擬合多元線性回歸
#這裡將處於線性關係臨界值的 depth 也算在內
fit2 <- lm(fish~acre+depth+no3, data = dat)
summary(fit2) #展示擬合方程的簡單統計
整體而言，去除了3個不顯著的預測變量（環境因素）後，方程的R2並未發生明顯變化，也就是精度未降低，且全模型p<0.001仍然是顯著的。
如果和上文使用全部環境因素擬合的最初的回歸方程相比，肯定更期望選擇後者，因為在二者精度相差不大時，後者更簡約、更容易理解且解釋輕鬆。
如果期望通過一些統計檢驗方法比較二者，以選擇更優的一個，如下提供兩個示例參考。
前後兩個多元線性回歸可視為嵌套模型，所謂嵌套模型，即其中之一的一些項完全包含在另一模型中。此時可以通過anova()比較兩個嵌套模型的擬合優度，以檢驗不含do2（水域溶解氧含量）、so4（水域硫酸鹽濃度）和temp（水域溫度）的回歸方程與包含這3個變量的回歸方程的預測效果是否一樣好。
還可以通過如AIC（Akaike Information Criterion，赤池信息準則）等值進行比較。AIC會在回歸方程複雜性與其擬合優度之間進行權衡，AIC值較小的回歸方程優先選擇，表明較少的預測變量已經獲得了足夠的擬合度。
#anova() 比較兩個嵌套模型的擬合優度
anova(fit2, fit1)
#AIC 比較兩個回歸方程
AIC(fit2, fit1)
anova()對是否應該添加變量do2、so4和temp到多元線性回歸方程中進行了檢驗，由於p>0.05，表明添加這3個變量前後兩回歸方程的差別不大，表明可以將這3個變量直接刪除。類似地，AIC值也顯示剔除這3個變量更佳。
本篇未涉及預測變量（環境因素）間的交互作用，帶交互項的多元回歸在放下節展示。需注意，考慮交互項時，就不再是線性方程了。
本篇未探討預測變量（環境因素）間可能存在的共線性問題，它比較複雜，關於共線性以及變量選擇等過程，將在下下節細講。
http://www.biostathandbook.com/multipleregression.htmlRobert I. Kabacoff. R語言實戰（第二版）（王小寧 劉擷芯 黃俊文 等 譯）. 人民郵電出版社, 2016.
生物空間站致力於分享生信技能幹貨、實用軟體教程、科研經驗、發布最新科研熱點資訊、解讀生命科學前沿進展、最新實驗技術。溫馨提示：添加主編微信，邀請你進生物空間站交流群，獲取更多有趣的知識及資料。添加時，請備註一下單位/學校+專業。
轉載、合作等事宜請聯繫BioSpaceStation