R與生物專題 | 第三十六講 R-多元線性回歸中的交互作用及更優模型選擇

在「R與生物統計專題」中，我們會從介紹R的基本知識展開到生物統計原理及其在R中的實現。以從淺入深，層層遞進的形式在投必得學術公眾號更新。

在前一講中，我們介紹了如何用多個預測變量（x）建立用於預測連續型數值的結果變量（y）的多元線性回歸模型（第三十五講 R-多元線性回歸）。

例如，要預測血糖，根據在胰島素水平和懷孕次數，模型方程式為Glucose = b0 + b1*Insulin + b2*Pregnancies，其中b0是截距；b1 和b2是分別與預測變量血壓、懷孕次數和懷孕次數相關聯的回歸係數。

上面的方程，也稱為加法模型，僅研究預測變量的主要影響。它的前提條件是，假設給定的預測變量與結果之間的關係獨立於其他預測變量。

但是，各個預測變量之間是不是相互獨立的呢？如果懷孕次數的增加可以提高胰島素水平對血糖的影響，則認為懷孕次數與胰島素水平之間有協同效應，在統計學上，稱為交互作用。

具有兩個預測變量（x1和x2）之間的相互作用的多元線性回歸方程可寫為：

y = b0 + b1*x1 + b2*x2 + b3*(x1*x2)

考慮我們的示例，它變為：

Glucose = b0 + b1*Insulin + b2* Pregnancies + b3*(Insulin* Pregnancies)

也可以寫成：

Glucose = b0 + (b1 + b3* Pregnancies)*Insulin + b2* Pregnancies

或為：

Glucose = b0 + b1*Insulin + (b2 +b3*Insulin)* Pregnancies

b3 可以解釋為胰島素提高，而懷孕次數額外增加一個單位（反之亦然）。

tidyverse 用於數據可視化

caret 用於簡化機器學習的工作流程

library(tidyverse)library(caret)

我們將使用一組糖尿病的數據,它包含768人糖尿病相關的數據，包括懷孕情況，血糖，血壓，皮膚厚度，胰島素水平，懷孕次數，糖尿病譜系功能，年齡和糖尿病診斷結果（Outcome）。

如需獲取數據diabetes.csv，請關注投必得醫學公眾號，後臺回復「diabetes.csv」獲取數據。

導入數據：

my_data<-read.csv('diabetes.csv')
檢查數據：
[1] 768   9
輸出結果
Pregnancies Glucose BloodPressure SkinThickness Insulin  PREGNANCIES DiabetesPedigreeFunction Age Outcome1           6     148            72            35       0 33.6                    0.627  50       12           1      85            66            29       0 26.6                    0.351  31       03           8     183            64             0       0 23.3                    0.672  32       14           1      89            66            23      94 28.1                    0.167  21       05           0     137            40            35     168 43.1                    2.288  33       16           5     116            74             0       0 25.6                    0.201  30       0
數據清理
new_data<-my_data[my_data$Glucose>0 & my_data$Insulin >0,]dim(new_data)
#[1] 393   9

研究問題：胰島素水平和懷孕次數、及其交互作用預測血糖水平。
set.seed(123)training.samples <- new_data$Glucose %>%createDataPartition(p = 0.8, list = FALSE)train.data  <- new_data[training.samples, ]test.data <- new_data[-training.samples, ]

標準線性回歸模型可以計算如下：
model1 <- lm(Glucose ~ Insulin + Pregnancies, data = train.data)summary(model1)
輸出結果
Call:lm(formula = Glucose ~ Insulin + Pregnancies, data = train.data)

Residuals:Min      1Q  Median      3Q     Max-59.368 -17.110  -3.895  13.805  83.981

Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 93.91665    2.61539  35.909  < 2e-16 ***Insulin      0.15017    0.01182  12.700  < 2e-16 ***Pregnancies  1.55894    0.42649   3.655 0.000301 ***Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1

Residual standard error: 24.45 on 312 degrees of freedomMultiple R-squared:  0.3731,  Adjusted R-squared:  0.369F-statistic: 92.83 on 2 and 312 DF,  p-value: < 2.2e-16
predictions <- model1 %>% predict(test.data)
# (a) 預測誤差, RMSERMSE(predictions, test.data$Glucose)
[1] 25.86233
R2(predictions, test.data$Glucose)
[1] 0.3167198

在R中，使用*運算符號來表示變量之間的交互作用：
model2 <- lm(Glucose ~ Insulin + Pregnancies + BloodPressure: Pregnancies,data = new_data)
model2 <- lm(Glucose ~ Insulin* Pregnancies, data = train.data)summary(model2)
輸出結果
Call:lm(formula = Glucose ~ Insulin * Pregnancies, data = train.data)

Residuals:Min      1Q  Median      3Q     Max-56.176 -16.922  -4.239  13.203  85.592

Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept)         92.133076   3.301245  27.909   <2e-16 ***Insulin              0.162520   0.018282   8.890   <2e-16 ***Pregnancies          2.142610   0.784902   2.730   0.0067 **Insulin:Pregnancies -0.003746   0.004229  -0.886   0.3763Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1

Residual standard error: 24.45 on 311 degrees of freedomMultiple R-squared:  0.3746,  Adjusted R-squared:  0.3686F-statistic: 62.11 on 3 and 311 DF,  p-value: < 2.2e-16
# 用該模型計算驗證數據集中預測值
# 模型在驗證數據集中的預測評估情況
predictions <- model2 %>% predict(test.data)RMSE(predictions, test.data$Glucose)
[1] 25.61056
# 計算 R平方
R2(predictions, test.data$Glucose)
[1] 0.3299708


可以看出，當我們考慮胰島素水平和懷孕次數的交互作用時，該交互作用並不能預測血糖水平，其係數與0間沒有差異。即最終模型中，我們不應該考慮胰島素水平和懷孕次數間的交互作用，其對血糖沒有交互作用。
我們可以進一步使用anova()函數，對兩個預測模型進行對比比較其R平方。
輸出結果
Analysis of Variance TableModel 1: Glucose ~ Insulin + PregnanciesModel 2: Glucose ~ Insulin * PregnanciesRes.Df    RSS Df Sum of Sq      F Pr(>F)1    312 1864542    311 185984  1    469.35 0.7848 0.3763
從結果我們可以看出，model2並沒有比model1好。殘差平方和（Residual sum of squares，RSS）在model1中為186454，在model2中為185984。兩者進行卡方檢驗，P = 0.2763，即兩者殘差平方和沒有統計學差異。
所以，最後更好的模型為model1。
有時交互作用項很重要，但不是主要影響。如果模型中交互作用沒有統計學意義，最終的模型中，不應該加入交互作用項。

如果模型中包括交互作用，並且交換作用有統計學意義，那麼在模型中，即便原預測因素的係數與0沒有統計學差異，也應該包括在最終的模型中。

1. Alboukadel Kassambara, Machine Learning Essentials: Practical Guide in R

好了，本期講解就先到這裡。小夥伴們趕緊試起來吧。
在之後的更新中，我們會進一步為您介紹R的入門，以及常用生物統計方法和R實現。歡迎關注，投必得學術手把手帶您走入R和生物統計的世界。
提前預告一下，下一講我們將學習多元線性回歸中的多重共線性和方差膨脹因子。
當然啦，R語言的掌握是在長期訓練中慢慢積累的。一個人學習太累，不妨加入「R與統計交流群」，和數百位碩博一起學習。

快掃二維碼撩客服，
帶你進入投必得學術交流群，
讓我們共同進步！
↓↓

來源：投必得醫學R與生物統計
聲明：本文僅做學術分享，版權歸原作者所有，並不代表本平臺的觀點，如有侵權，請先聯繫topedit2021刪除，萬分感謝！
發表SCI論文很迷茫？
來找「投必得」幫忙↓↓↓