引用本文： 李多多, 俞興, 韓晟, 朱賀, 苑藝, 沈捷, 林景峰, 黎霞, 甘葉娜, 劉建平. 多元線性回歸分析數據可視化的 R Studio 軟件實踐. 中國循證醫學雜志, 2021, 21(4): 482-490. doi: 10.7507/1672-2531.202008172 復制
數據可視化是將數據以視覺表現展現的科學研究技術，用繪制圖表的方式來表示數據及其之間關系的方法，讓研究者可準確而快捷地了解數據的趨勢、異常值和關聯性等信息[1]。人類的眼睛對圖案和色彩較為敏感，數據可視化技術是構建人與數據合作的重要橋梁。在大數據時代，數據可視化技術應用廣泛，在分析海量信息后提出決策分析方面表現出明顯的優勢[2]。多元線性回歸分析（multiple linear regression analysis，MLRA）是研究多個自變量與一個因變量間是否存在線性變化關系的研究方法，通過建立回歸模型來衡量不同自變量對因變量的影響能力，是醫學研究中最常用的統計學方法之一[3]。MLRA 首先依據研究目的確定自變量（Xi）與因變量（Y）的函數關系，然后進行建模、參數估計和檢驗，最后進行模型的賦值運用[4]。MLRA 的本質是構建具有實際臨床應用價值的模型。而確定模型能否構建和使用，在建模前需要滿足變量間的線性關系和方差齊性等條件，在建模后還需進行顯著性檢驗和殘差分析等，通過數據可視化技術可直接觀察模型是否構建合理，有無異常值和異常趨勢等，較單純公式計算的結果更直觀。
R 語言軟件是由新西蘭奧克蘭大學 Ross Ihaka 和 Robert Gentleman 開發，R 開發核心團隊負責完善，目前眾多研究領域廣泛使用的一個免費開放源代碼統計軟件（官方網站：https://www.r-project.org/），常用于統計分析、繪圖、數據挖掘等[5]。目前隨著計算機技術的快速發展，基于 R 語言集成開發環境的 R Studio 軟件，是一組可幫助我們直接便捷高效利用 R 來進行軟件開發的專業集成工具，該軟件可在各種 UNIX 平臺，Windows 和 MacOS 上進行編譯和運行，它在傳統 R 語言統計基礎上更提升了數據可視化功能，使曾經可能難以實現的制圖過程輕松求解[6]。本文結合具體實例介紹如何使用 R 語言代碼，基于 R Studio 實現 MLRA 的數學建模、賦值應用及相關數據可視化，以期有效提升 MLRA 在醫學研究工作中的服務能力。
在安裝 R Studio 軟件（下載網址：http://www.R Studio.com/products/R Studio/）前我們須在電腦上安裝 R 語言軟件。本文所選用的數據為孫振球和徐勇勇主編的“十二五”規劃教材《醫學統計學》第十五章多元線性回歸分析中的教學案例：27 例糖尿病患者的血糖及有關變量的測量結果[3]（表 1），數據需要提前錄入 Excel 軟件并保存為擴展名為 csv（逗號分隔）類型的文件，本例我們在 D 盤（亦可根據個人電腦硬盤分區進行選擇）建立名為 GLU.csv 的 Excel 表，對“列”名稱進行修改，采用英文縮寫便于代碼書寫和軟件運行，具體為“膽固醇”定義為“CHOL”，“甘油三酯”定義為“TG”，“胰島素”定義為“Insulin”，“糖化血紅蛋白”定義為“HbA1c”，“空腹血糖”定義為“GLU”，并對應每列錄入數據。以上步驟實現了數據清洗、整理，以便 R 軟件的調用和讀取。對以上數據進行分析數據格式化（表 2）。根據研究目的和內容，我們確定多個自變量為Xi（i=1，2，3…i）和它們對應的因變量為Y。對應法則f的解析式Y=f（Xi）表示[7]。
表1
27 例糖尿病患者血糖及可能相關變量
表2
多元回歸分析數據格式
1 加載數據到 R Studio 中
在構建回歸方程之前，我們首先需要為所有確定變量進行賦值。啟動 R Studio 軟件后，我們進行數據加載。
輸入命令：
> fichier <- "D:/GLU.csv"
> mydata <-read.table（fichier,header=TRUE,sep=","）
選擇 D 盤 GLU.csv 為目標文件，讀取文件數據，header=TRUE 為邏輯值，以指定不同變量名稱均由文件第一行給出，sep=“,”意為每一行的各個值之間都以“,”來進行分隔。
此外還需載入后續工作中用到的多個函數包（packages）。
輸入命令：
> library（ggplot2）
> library（reshape）
> library（corrplot）
函數包是由廣大第三方開發者根據一定的工作主題，編寫一系列函數集成體系化的工具套裝，并發布在 CRAN（The Comprehensive R Archive Network）工具倉庫中。各研究者可根據當前工作需要，自由組合選用不同函數包中的函數。因此，函數包是 R 語言靈活性和活躍生命力的體現，為科研工作帶來極大的便利。
正式工作開展前，需要將計劃使用到的函數包用 library（）函數加載到工作環境中。本次演示中用到的第三方工具庫包括：ggplot2、reshape 和 corrplot。
2 確認數據是否滿足多元回歸分析條件
2.1 全域關系
我們通過命令求解所有變量之間的關系，并繪制全域散點圖觀察，確認數據是否滿足多元回歸分析條件。
輸入命令：
> attach（mydata）
> pairs（mydata［,2:6］）
> cor（mydata［,2:6］）
輸出結果：見圖 1。
圖1
全域散點圖
這個數據集在第一步中被命名為“mydata”，我們將該數據集通過連接函數 attach（）連接到搜索路徑中，后續使用該數據時，不必引用數據集名稱，可只引用每一列的列名，即字段名。pairs（）函數對每一對變量生成全域散點圖，這類圖形可將變量之間的關系形象化，再連接到統計函數 cor（）求解出上圖對應的代表“斜率”的矩陣，為下文的高級運算提供重要的基礎數據。
2.2 線性擬合
觀察上圖似乎難以發現變量之間存在明顯的線性關系，這就需要我們將散點進一步擬合尋找自變量與因變量是否可能存在線性關系，以確認觀察值的獨立性、因變量的正態性、變量間的線性關系是否滿足進一步構建多元回歸分析模型的條件。
輸入命令：
> add.cor <- function（x,y）{
usr <- par（"usr"）; on.exit（par（usr））
par（usr = c（0, 1, 0, 1））
text（0.5,0.5,round（cor（x,y）,2））}
> add.abline <- function（x, y）{
points（x, y）
abline（lm（y～x）, col = 'red'）}
> pairs（mydata［,2:6］,lower.panel = add.abline,upper.panel=add.cor）
輸出結果：見圖 2。
圖2
全域散點線性擬合圖
執行上述命令后生成全域散點線性擬合圖，觀察圖中最下一行除 Insulin 外其余三項均與 GLU 呈現正相關，且散點具有一定的線性回歸趨勢，故本例數據可建立線性模型。
3 建立多元線性回歸模型
在實例中我們設X1為膽固醇值（變量 CHOL），X2為甘油三酯值（變量 TG），X3為胰島素值（變量 Insulin），X4為糖化血紅蛋白百分比（變量 HbA1c），則解析式Y=f（Xi）可具體為多元線性回歸方程[3]：
=b0+b1X1+b2X2+b3X3+b4X4，亦可改寫為GLU（Y/X1，X2，X3，X4）=b0+b1CHOL+b2TG+b3Insulin+b4HbA1c，b0，b1，b2，b3和 b4為下文模型中偏回歸系數β0，β1，β2，β3和β4參數的估計值；
表示在X1，X2，X3，X4取值時Y的平均值。
3.1 參數估計
公式中的回歸參數b0，b1，b2，b3，b4都是我們不知道的，所以參數估計就是要通過自變量的不同數據來估計這些參數從而確定X和Y之間的關系。參數估計的核心目標是要盡可能計算出一條直線，并使這條直線上的每個點的Y值和實際數據Y值之差的平方和最小，即最小二乘法（least sum of squares，LS）：令（Y1實際?Y1預測）2+（Y2實際?Y2預測）2+···+（Yn實際?Yn預測）2值最小。本例中我們可通過 R 軟件的 lm（）函數來求得模型中參數β0，β1，β2，β3，β4進行估計，從而得到多元線性回歸方程b0，b1，b2，b3，b4參數值。
輸入命令：
> model1 <- lm（GLU～CHOL+TG+Insulin+HbA1c）
> model1
輸出結果：
R 軟件的 lm（Y～X）函數為線性模型的估計命令，model1 為模型構建命令，相應數據為模型參數的估計值，具體為b0=5.9433，b1=0.1424，b2=0.3515，b3=?0.2706，b4=0.6382，帶入數據后擬合的 MLRA 方程為：
=5.9433+0.1424X1+0.3515X2-0.2706X3+0.6382X4，如果我們想進一步得出圖示上述模型中 27 位患者X1-4對Y的函數映射關系，則執行如下操作：
輸入命令：
> datamel <- melt（mydata,id.vars = c（"ID"））
> ggplot（datamel, aes（x = ID））+
geom_line（aes（y = value, col = variable）,size = 1）
輸出結果：見圖 3。
圖3
模型函數 X1-4 與 Y 映射關系
3.2 多重共線性的識別與處理
如一些自變量之間存在較強線性關系時，我們稱它們之間存在多重共線性（multicollinearity），其共線性程度越強越容易引起建模的不良后果，如參數估計值的標準誤過大、t檢驗結果過小、回歸參數正負符號與客觀實際不一致等。因此我們首先要進行共線性的識別，通過矩陣運算對比自變量之間的相關系數，如果其絕對值越接近 1 則認為二者共線性越強。
輸入命令：
> corr <- cor（mydata［,2:6］）
> corr
輸出結果：
輸入命令：
> par（mfrow = c（1, 2））
> corrplot（corr,method="number"）+
corrplot（corr,method="circle"）
輸出結果：見圖 4。
圖4
相關系數圖
圖 4 為矩陣運算的數值，其顏色參照右側標尺，越接近深紅色和深藍色表明共線性越強，圖 4（左）中除左上至右下的對角線表示自身比較外，其余數值均小于 0.65，說明自變量間共線性較弱，不需要對方程進行調整；圖 4（右）為數值的可視化呈現使讀者一目了然。
3.3 顯著性檢驗
MLRA 的顯著性檢驗是關于回歸模型中因變量與自變量間函數關系的假設是否顯著的檢驗，也就是探討二者之間的線性關系是否密切的檢驗。同樣是需要經過t檢驗，F檢驗，和決定系數R2相關檢驗[3]。在 R 軟件中這三種檢驗的計算可同時采用函數 summary（）來完成。
輸入命令：
> summary（model1）
輸出結果：
從輸出的結果可觀察到標準誤差系數表（coefficients）、殘差的五數概括（residuals）、決定系數 R2（multiple R-squared，adjusted R-squared）和殘余標準誤差（residual standard error）等，回歸方程參數的估計值在 Estimate 中也有列出，對應于H0：βi =0 和H1：βi≠0 采用t檢驗統計量的輸出數值在t value 中列出。Std. Error 為回歸系數估計值b（Estimate）的標準誤，故t檢驗的公式與上述結果的對應為t value=Estimate/ Std. Error，依據t界值表[1]得 t0.05/2，22=2.074，t4>|t3|>2.074，P值小于 0.05，表明b3和b4顯著，對同一資料不同自變量的t值是可進行相互比較的，t的絕對值越大就說明自變量對因變量Y的回歸起的作用越大。本例說明t4對應的糖化血紅蛋白值對空腹血糖的作用較其他變量更大，綜合結果我們可進一步明確，對空腹血糖影響大小的順序依次為糖化血紅蛋白（HbA1c）、胰島素（Insulin）、甘油三酯（CHOL）和總膽固醇（TG）；結果中我們還還可得到多元線性回歸決定系數R2[3]，即 R2=SS回/SS總=133.7100/222.5510=0.6008，即輸出數據中的 multiple R-squared，調整后的 R2為 adjusted R-squared=0.5282，這表明表示空腹血糖含量（Y）變化的 52.82% 可通過預測變量“CHOL”、“TG”、“Insulin”和“HbA1c”之間的線性關系來解釋，故提示模型擬合度不良。
3.4 殘差分析
殘差即觀察值與估計值（擬合值）之間的差，其蘊含著預設模型基本假設的重要信息。如果多元回歸模型正確，則我們可將殘差看作誤差的觀測值，且應符合模型的假設條件。我們利用殘差提供的信息，進一步檢測對應模型假設的合理性和其數據可靠性的方法，我們稱之為殘差分析，它是檢查資料是否符合模型條件的有效工具，如果資料與模型的假設偏離較大，則采用最小二乘估計、假設檢驗及區間估計有可能出現問題。我們可對函數 lm（）結合函數 plot（）方法得到殘差圖，殘差圖包括兩個子圖形，即殘差圖—擬合值圖和殘差—正態 QQ 圖。
輸入命令：
> par（mfrow=c（1,2））
> plot（model1）
輸出結果：見圖 5。
圖5
標準化殘差圖和殘差正態 QQ 圖
圖 5（左）理想情況為數據點均勻分布在 y=0 兩側，呈現出隨機的分布，紅色線呈現出一條平穩的曲線并沒有明顯的形狀特征；圖 5（右）理想情況為數據點按對角直線排列，趨于一條直線，并被對角直接穿過，直觀上符合正態分布；求解運算后輸出圖形具體情況為圖 5（左）中散點分布不理想，紅色線呈一個近似的“U”形，提示該模型中可能缺失了一個重要的未知變量；圖 5（右）中，除殘差最大和最小的 26 與 13，其余基本符合正態分布。
3.5 模型確立
多元線性回歸模型的構建形式為：Y=β0+β1X1+β2X2+…+βmXm+e[1]，β0為截距，β1，β2，…βm即回歸系數（partial regression coefficient），e即去除m個自變量對Y影響后的殘差，因為在模型和函數關系構建之初，殘差e是不可觀測的，故在構建回歸方程時不需要考慮殘差。根據上文的顯著性檢驗和殘差分析，MLRA 方程的各自變量參數均無修正，則模型即等于上述 MLRA 方程
=5.9433+0.1424X1+0.3515X2?0.2706X3+0.6382X4。
4 根據實際情況運用模型進行預測
多元線性回歸模型可應用于統計區間的預測，即采用多個自變量來計算因變量變化的統計區間。我們可把預設的自變量X代入回歸方程，通過 R Studio 軟件的函數 predict（）得到Y的個體預測值、Y個體值 95% 預測區間和Y總體均數的 95% 置信區間。如果我們僅針對某一患者已有自變量的觀測值對Y進行預測，我們可直接將數據代入方程對應的 R 軟件命令進行求解，如我們假設有一位患者，其化驗單顯示總膽固醇 10.50 mmol/L，甘油三酯 9.50 mmol/L，胰島素 1.82 μU/mL，糖化血紅蛋 9.85%，我們就可依據此進行估計空腹血糖的上述區間。
輸入命令：
> newdata <- data.frame（CHOL=10.50,TG=9.50,Insulin=1.82,HbA1c=9.85）
> predict（model1,newdata）
> predict（model1,newdata,interval = "pred",level =0.95）
> predict（model1,newdata,interval = "conf",level =0.95）
輸出結果：
我們輸入 predict（model1，newdata）命令可計算出Y的個體預測值為 16.5717，為避免誤差，我們可在“predict”方法中指定“interval”的參數，通過“level”參數指定置信水平為 0.95，對預測區間可使用“interval="pred"”，求解得（11.56282，21.58057）；Y總體均數的 95% 置信區間使用“interval="conf"”，求解得（13.79310，19.35029）。
如果我們欲對多位患者的自變量進行Y值預測，則應建立數據庫以便 R 軟件整體導入進行運算和分析。我們仍以上文 27 例患者觀測指標的數據庫為例，為避免數據重合和圖像重合難以閱讀觀察，我們首先計算 Y 值的預測值和預測區間。數據導入及求解步驟如下。
輸入命令：
> library（ggplot2）
> library（reshape2）
> myprediction<-predict（model1,interval="pred",level =0.95）
> head（myprediction,27）
輸出結果：
上述結果是 myprediction 的賦值矩陣，我們需要把它變成數據框，然后和原始數據的Y進行合并，即 mydata 中的 GLU 值。
輸入命令：
> myprediction <- data.frame（myprediction）
> myprediction <- cbind（y = mydata$GLU, myprediction）
> myprediction
輸出結果：
上述表格是寬表格形式，如果進行可視化制圖，我們需要將表格轉換成長表格形式才能滿足函數的讀取要求。
輸入命令：
> myprediction_long <- melt（myprediction,id.vars = 'x'）
> myprediction_long$x <- as.numeric（myprediction_long$x）
> ggplot（myprediction_long）+
geom_line（aes（x = x, y = value, color = variable）,size=1）
輸出結果：見圖 6。
圖6
27 例患者空腹血糖（Y）實測值，預測值和預測區間上下限
如果我們將上述代碼中的“myprediction”替換為“myconfidence”，“interval="pred"”改為“interval="conf"”，“預測區間上下限”改為“置信區間上下限”，則再進一步計算出置信區間上下限及繪圖（圖 7）。
圖7
27 例患者空腹血糖（Y）實際值，預測值和置信區間上下限
5 討論
多元線性回歸分析（MLRA）的本質是為了估計或預測，如文中在建模后代入新患者的糖化血紅蛋白觀測值以推算空腹血糖預測值的統計區間，需要應用矩陣運算，且在建模前構建預測回歸方程時應選擇具有較高 R2值的方程以確保方程在引入自變量計算后的結果具有統計學意義（P<0.05）[3]。
MLRA 在醫學領域有著廣泛的應用，原因在于多因素混雜分析是醫學研究中經常遇到的問題。我們已知的絕大多數疾病都是多種原因共同作用的結果，其預后也往往是由多因素所決定的[8]，如本文示例中影響空腹血糖的因素不止上述幾種，還可能有年齡、飲食習慣、工作環境和家族史等，在此基礎上，我們在這些因素中還要明確不同因素的影響權重和可能存在的邏輯關系。這些因素可能影響著臨床試驗的各個方面，造成我們可能需要面臨在難以保證基線相同的混雜情況下對不同干預方法進行比較這樣的難題，以上這些問題很多都可利用多元線性回歸分析來處理。我們控制其混雜因素的一個簡單而有效的方法就是將具體觀測值通過函數的映射關系代入方程進行分析，但這就決定了 MLRA 計算出的預測區間不能說明實際的效果，僅可用于佐證之需。
目前，多元線性回歸方程假設檢驗應用較為普遍的軟件是 SPSS 軟件、SAS 軟件和 Stata 軟件等。SPSS 軟件雖然簡便，但是受其模塊化界面操作的限制，難以完成多種統計，SAS、Stata 軟件雖然具有編程等強大功能，但其高昂的租賃年費也常令我們的工作面臨困境。R 語言因其邏輯運算功能強大，統計制圖多樣，編程語言簡明科學等優勢，始終受到科研工作者的青睞，但 R 語言傳統自帶的環境操作起來較為不便，R Studio 作為 R 語言的高性能操作平臺。它兼容了 R 語言數據處理、計算、制圖以及通過全球不同地區 Comprehensive R Archive Network（CRAN）鏡像站下載升級數據包以隨時不斷豐富數學計算函數模塊等幾乎全部功能外，簡化了編程的難度，提升了操作的便捷性，且它還具有更加優良的數據可視化功能，更具有支持純 R 腳本、腳本文檔混排（rmarkdown）、腳本文檔混排成書（bookdown）、交互式網絡應用（shiny）等優勢。其唯一的障礙即要求我們掌握 R 語言的編寫，才能完成所需要的統計分析和圖形繪制。故本文結合實例提供讀者代碼的具體編寫方法，希望 R Studio 軟件今后能為更多研究者所采用，為醫學研究的發展做出貢獻。
數據可視化是將數據以視覺表現展現的科學研究技術，用繪制圖表的方式來表示數據及其之間關系的方法，讓研究者可準確而快捷地了解數據的趨勢、異常值和關聯性等信息[1]。人類的眼睛對圖案和色彩較為敏感，數據可視化技術是構建人與數據合作的重要橋梁。在大數據時代，數據可視化技術應用廣泛，在分析海量信息后提出決策分析方面表現出明顯的優勢[2]。多元線性回歸分析（multiple linear regression analysis，MLRA）是研究多個自變量與一個因變量間是否存在線性變化關系的研究方法，通過建立回歸模型來衡量不同自變量對因變量的影響能力，是醫學研究中最常用的統計學方法之一[3]。MLRA 首先依據研究目的確定自變量（Xi）與因變量（Y）的函數關系，然后進行建模、參數估計和檢驗，最后進行模型的賦值運用[4]。MLRA 的本質是構建具有實際臨床應用價值的模型。而確定模型能否構建和使用，在建模前需要滿足變量間的線性關系和方差齊性等條件，在建模后還需進行顯著性檢驗和殘差分析等，通過數據可視化技術可直接觀察模型是否構建合理，有無異常值和異常趨勢等，較單純公式計算的結果更直觀。
R 語言軟件是由新西蘭奧克蘭大學 Ross Ihaka 和 Robert Gentleman 開發，R 開發核心團隊負責完善，目前眾多研究領域廣泛使用的一個免費開放源代碼統計軟件（官方網站：https://www.r-project.org/），常用于統計分析、繪圖、數據挖掘等[5]。目前隨著計算機技術的快速發展，基于 R 語言集成開發環境的 R Studio 軟件，是一組可幫助我們直接便捷高效利用 R 來進行軟件開發的專業集成工具，該軟件可在各種 UNIX 平臺，Windows 和 MacOS 上進行編譯和運行，它在傳統 R 語言統計基礎上更提升了數據可視化功能，使曾經可能難以實現的制圖過程輕松求解[6]。本文結合具體實例介紹如何使用 R 語言代碼，基于 R Studio 實現 MLRA 的數學建模、賦值應用及相關數據可視化，以期有效提升 MLRA 在醫學研究工作中的服務能力。
在安裝 R Studio 軟件（下載網址：http://www.R Studio.com/products/R Studio/）前我們須在電腦上安裝 R 語言軟件。本文所選用的數據為孫振球和徐勇勇主編的“十二五”規劃教材《醫學統計學》第十五章多元線性回歸分析中的教學案例：27 例糖尿病患者的血糖及有關變量的測量結果[3]（表 1），數據需要提前錄入 Excel 軟件并保存為擴展名為 csv（逗號分隔）類型的文件，本例我們在 D 盤（亦可根據個人電腦硬盤分區進行選擇）建立名為 GLU.csv 的 Excel 表，對“列”名稱進行修改，采用英文縮寫便于代碼書寫和軟件運行，具體為“膽固醇”定義為“CHOL”，“甘油三酯”定義為“TG”，“胰島素”定義為“Insulin”，“糖化血紅蛋白”定義為“HbA1c”，“空腹血糖”定義為“GLU”，并對應每列錄入數據。以上步驟實現了數據清洗、整理，以便 R 軟件的調用和讀取。對以上數據進行分析數據格式化（表 2）。根據研究目的和內容，我們確定多個自變量為Xi（i=1，2，3…i）和它們對應的因變量為Y。對應法則f的解析式Y=f（Xi）表示[7]。
表1
27 例糖尿病患者血糖及可能相關變量
表2
多元回歸分析數據格式
1 加載數據到 R Studio 中
在構建回歸方程之前，我們首先需要為所有確定變量進行賦值。啟動 R Studio 軟件后，我們進行數據加載。
輸入命令：
> fichier <- "D:/GLU.csv"
> mydata <-read.table（fichier,header=TRUE,sep=","）
選擇 D 盤 GLU.csv 為目標文件，讀取文件數據，header=TRUE 為邏輯值，以指定不同變量名稱均由文件第一行給出，sep=“,”意為每一行的各個值之間都以“,”來進行分隔。
此外還需載入后續工作中用到的多個函數包（packages）。
輸入命令：
> library（ggplot2）
> library（reshape）
> library（corrplot）
函數包是由廣大第三方開發者根據一定的工作主題，編寫一系列函數集成體系化的工具套裝，并發布在 CRAN（The Comprehensive R Archive Network）工具倉庫中。各研究者可根據當前工作需要，自由組合選用不同函數包中的函數。因此，函數包是 R 語言靈活性和活躍生命力的體現，為科研工作帶來極大的便利。
正式工作開展前，需要將計劃使用到的函數包用 library（）函數加載到工作環境中。本次演示中用到的第三方工具庫包括：ggplot2、reshape 和 corrplot。
2 確認數據是否滿足多元回歸分析條件
2.1 全域關系
我們通過命令求解所有變量之間的關系，并繪制全域散點圖觀察，確認數據是否滿足多元回歸分析條件。
輸入命令：
> attach（mydata）
> pairs（mydata［,2:6］）
> cor（mydata［,2:6］）
輸出結果：見圖 1。
圖1
全域散點圖
這個數據集在第一步中被命名為“mydata”，我們將該數據集通過連接函數 attach（）連接到搜索路徑中，后續使用該數據時，不必引用數據集名稱，可只引用每一列的列名，即字段名。pairs（）函數對每一對變量生成全域散點圖，這類圖形可將變量之間的關系形象化，再連接到統計函數 cor（）求解出上圖對應的代表“斜率”的矩陣，為下文的高級運算提供重要的基礎數據。
2.2 線性擬合
觀察上圖似乎難以發現變量之間存在明顯的線性關系，這就需要我們將散點進一步擬合尋找自變量與因變量是否可能存在線性關系，以確認觀察值的獨立性、因變量的正態性、變量間的線性關系是否滿足進一步構建多元回歸分析模型的條件。
輸入命令：
> add.cor <- function（x,y）{
usr <- par（"usr"）; on.exit（par（usr））
par（usr = c（0, 1, 0, 1））
text（0.5,0.5,round（cor（x,y）,2））}
> add.abline <- function（x, y）{
points（x, y）
abline（lm（y～x）, col = 'red'）}
> pairs（mydata［,2:6］,lower.panel = add.abline,upper.panel=add.cor）
輸出結果：見圖 2。
圖2
全域散點線性擬合圖
執行上述命令后生成全域散點線性擬合圖，觀察圖中最下一行除 Insulin 外其余三項均與 GLU 呈現正相關，且散點具有一定的線性回歸趨勢，故本例數據可建立線性模型。
3 建立多元線性回歸模型
在實例中我們設X1為膽固醇值（變量 CHOL），X2為甘油三酯值（變量 TG），X3為胰島素值（變量 Insulin），X4為糖化血紅蛋白百分比（變量 HbA1c），則解析式Y=f（Xi）可具體為多元線性回歸方程[3]：=b0+b1X1+b2X2+b3X3+b4X4，亦可改寫為GLU（Y/X1，X2，X3，X4）=b0+b1CHOL+b2TG+b3Insulin+b4HbA1c，b0，b1，b2，b3和 b4為下文模型中偏回歸系數β0，β1，β2，β3和β4參數的估計值；表示在X1，X2，X3，X4取值時Y的平均值。
3.1 參數估計
公式中的回歸參數b0，b1，b2，b3，b4都是我們不知道的，所以參數估計就是要通過自變量的不同數據來估計這些參數從而確定X和Y之間的關系。參數估計的核心目標是要盡可能計算出一條直線，并使這條直線上的每個點的Y值和實際數據Y值之差的平方和最小，即最小二乘法（least sum of squares，LS）：令（Y1實際?Y1預測）2+（Y2實際?Y2預測）2+···+（Yn實際?Yn預測）2值最小。本例中我們可通過 R 軟件的 lm（）函數來求得模型中參數β0，β1，β2，β3，β4進行估計，從而得到多元線性回歸方程b0，b1，b2，b3，b4參數值。
輸入命令：
> model1 <- lm（GLU～CHOL+TG+Insulin+HbA1c）
> model1
輸出結果：
R 軟件的 lm（Y～X）函數為線性模型的估計命令，model1 為模型構建命令，相應數據為模型參數的估計值，具體為b0=5.9433，b1=0.1424，b2=0.3515，b3=?0.2706，b4=0.6382，帶入數據后擬合的 MLRA 方程為：=5.9433+0.1424X1+0.3515X2-0.2706X3+0.6382X4，如果我們想進一步得出圖示上述模型中 27 位患者X1-4對Y的函數映射關系，則執行如下操作：
輸入命令：
> datamel <- melt（mydata,id.vars = c（"ID"））
> ggplot（datamel, aes（x = ID））+
geom_line（aes（y = value, col = variable）,size = 1）
輸出結果：見圖 3。
圖3
模型函數 X1-4 與 Y 映射關系
3.2 多重共線性的識別與處理
如一些自變量之間存在較強線性關系時，我們稱它們之間存在多重共線性（multicollinearity），其共線性程度越強越容易引起建模的不良后果，如參數估計值的標準誤過大、t檢驗結果過小、回歸參數正負符號與客觀實際不一致等。因此我們首先要進行共線性的識別，通過矩陣運算對比自變量之間的相關系數，如果其絕對值越接近 1 則認為二者共線性越強。
輸入命令：
> corr <- cor（mydata［,2:6］）
> corr
輸出結果：
輸入命令：
> par（mfrow = c（1, 2））
> corrplot（corr,method="number"）+
corrplot（corr,method="circle"）
輸出結果：見圖 4。
圖4
相關系數圖
圖 4 為矩陣運算的數值，其顏色參照右側標尺，越接近深紅色和深藍色表明共線性越強，圖 4（左）中除左上至右下的對角線表示自身比較外，其余數值均小于 0.65，說明自變量間共線性較弱，不需要對方程進行調整；圖 4（右）為數值的可視化呈現使讀者一目了然。
3.3 顯著性檢驗
MLRA 的顯著性檢驗是關于回歸模型中因變量與自變量間函數關系的假設是否顯著的檢驗，也就是探討二者之間的線性關系是否密切的檢驗。同樣是需要經過t檢驗，F檢驗，和決定系數R2相關檢驗[3]。在 R 軟件中這三種檢驗的計算可同時采用函數 summary（）來完成。
輸入命令：
> summary（model1）
輸出結果：
從輸出的結果可觀察到標準誤差系數表（coefficients）、殘差的五數概括（residuals）、決定系數 R2（multiple R-squared，adjusted R-squared）和殘余標準誤差（residual standard error）等，回歸方程參數的估計值在 Estimate 中也有列出，對應于H0：βi =0 和H1：βi≠0 采用t檢驗統計量的輸出數值在t value 中列出。Std. Error 為回歸系數估計值b（Estimate）的標準誤，故t檢驗的公式與上述結果的對應為t value=Estimate/ Std. Error，依據t界值表[1]得 t0.05/2，22=2.074，t4>|t3|>2.074，P值小于 0.05，表明b3和b4顯著，對同一資料不同自變量的t值是可進行相互比較的，t的絕對值越大就說明自變量對因變量Y的回歸起的作用越大。本例說明t4對應的糖化血紅蛋白值對空腹血糖的作用較其他變量更大，綜合結果我們可進一步明確，對空腹血糖影響大小的順序依次為糖化血紅蛋白（HbA1c）、胰島素（Insulin）、甘油三酯（CHOL）和總膽固醇（TG）；結果中我們還還可得到多元線性回歸決定系數R2[3]，即 R2=SS回/SS總=133.7100/222.5510=0.6008，即輸出數據中的 multiple R-squared，調整后的 R2為 adjusted R-squared=0.5282，這表明表示空腹血糖含量（Y）變化的 52.82% 可通過預測變量“CHOL”、“TG”、“Insulin”和“HbA1c”之間的線性關系來解釋，故提示模型擬合度不良。
3.4 殘差分析
殘差即觀察值與估計值（擬合值）之間的差，其蘊含著預設模型基本假設的重要信息。如果多元回歸模型正確，則我們可將殘差看作誤差的觀測值，且應符合模型的假設條件。我們利用殘差提供的信息，進一步檢測對應模型假設的合理性和其數據可靠性的方法，我們稱之為殘差分析，它是檢查資料是否符合模型條件的有效工具，如果資料與模型的假設偏離較大，則采用最小二乘估計、假設檢驗及區間估計有可能出現問題。我們可對函數 lm（）結合函數 plot（）方法得到殘差圖，殘差圖包括兩個子圖形，即殘差圖—擬合值圖和殘差—正態 QQ 圖。
輸入命令：
> par（mfrow=c（1,2））
> plot（model1）
輸出結果：見圖 5。
圖5
標準化殘差圖和殘差正態 QQ 圖
圖 5（左）理想情況為數據點均勻分布在 y=0 兩側，呈現出隨機的分布，紅色線呈現出一條平穩的曲線并沒有明顯的形狀特征；圖 5（右）理想情況為數據點按對角直線排列，趨于一條直線，并被對角直接穿過，直觀上符合正態分布；求解運算后輸出圖形具體情況為圖 5（左）中散點分布不理想，紅色線呈一個近似的“U”形，提示該模型中可能缺失了一個重要的未知變量；圖 5（右）中，除殘差最大和最小的 26 與 13，其余基本符合正態分布。
3.5 模型確立
多元線性回歸模型的構建形式為：Y=β0+β1X1+β2X2+…+βmXm+e[1]，β0為截距，β1，β2，…βm即回歸系數（partial regression coefficient），e即去除m個自變量對Y影響后的殘差，因為在模型和函數關系構建之初，殘差e是不可觀測的，故在構建回歸方程時不需要考慮殘差。根據上文的顯著性檢驗和殘差分析，MLRA 方程的各自變量參數均無修正，則模型即等于上述 MLRA 方程=5.9433+0.1424X1+0.3515X2?0.2706X3+0.6382X4。
4 根據實際情況運用模型進行預測
多元線性回歸模型可應用于統計區間的預測，即采用多個自變量來計算因變量變化的統計區間。我們可把預設的自變量X代入回歸方程，通過 R Studio 軟件的函數 predict（）得到Y的個體預測值、Y個體值 95% 預測區間和Y總體均數的 95% 置信區間。如果我們僅針對某一患者已有自變量的觀測值對Y進行預測，我們可直接將數據代入方程對應的 R 軟件命令進行求解，如我們假設有一位患者，其化驗單顯示總膽固醇 10.50 mmol/L，甘油三酯 9.50 mmol/L，胰島素 1.82 μU/mL，糖化血紅蛋 9.85%，我們就可依據此進行估計空腹血糖的上述區間。
輸入命令：
> newdata <- data.frame（CHOL=10.50,TG=9.50,Insulin=1.82,HbA1c=9.85）
> predict（model1,newdata）
> predict（model1,newdata,interval = "pred",level =0.95）
> predict（model1,newdata,interval = "conf",level =0.95）
輸出結果：
我們輸入 predict（model1，newdata）命令可計算出Y的個體預測值為 16.5717，為避免誤差，我們可在“predict”方法中指定“interval”的參數，通過“level”參數指定置信水平為 0.95，對預測區間可使用“interval="pred"”，求解得（11.56282，21.58057）；Y總體均數的 95% 置信區間使用“interval="conf"”，求解得（13.79310，19.35029）。
如果我們欲對多位患者的自變量進行Y值預測，則應建立數據庫以便 R 軟件整體導入進行運算和分析。我們仍以上文 27 例患者觀測指標的數據庫為例，為避免數據重合和圖像重合難以閱讀觀察，我們首先計算 Y 值的預測值和預測區間。數據導入及求解步驟如下。
輸入命令：
> library（ggplot2）
> library（reshape2）
> myprediction<-predict（model1,interval="pred",level =0.95）
> head（myprediction,27）
輸出結果：
上述結果是 myprediction 的賦值矩陣，我們需要把它變成數據框，然后和原始數據的Y進行合并，即 mydata 中的 GLU 值。
輸入命令：
> myprediction <- data.frame（myprediction）
> myprediction <- cbind（y = mydata$GLU, myprediction）
> myprediction
輸出結果：
上述表格是寬表格形式，如果進行可視化制圖，我們需要將表格轉換成長表格形式才能滿足函數的讀取要求。
輸入命令：
> myprediction_long <- melt（myprediction,id.vars = 'x'）
> myprediction_long$x <- as.numeric（myprediction_long$x）
> ggplot（myprediction_long）+
geom_line（aes（x = x, y = value, color = variable）,size=1）
輸出結果：見圖 6。
圖6
27 例患者空腹血糖（Y）實測值，預測值和預測區間上下限
如果我們將上述代碼中的“myprediction”替換為“myconfidence”，“interval="pred"”改為“interval="conf"”，“預測區間上下限”改為“置信區間上下限”，則再進一步計算出置信區間上下限及繪圖（圖 7）。
圖7
27 例患者空腹血糖（Y）實際值，預測值和置信區間上下限
5 討論
多元線性回歸分析（MLRA）的本質是為了估計或預測，如文中在建模后代入新患者的糖化血紅蛋白觀測值以推算空腹血糖預測值的統計區間，需要應用矩陣運算，且在建模前構建預測回歸方程時應選擇具有較高 R2值的方程以確保方程在引入自變量計算后的結果具有統計學意義（P<0.05）[3]。
MLRA 在醫學領域有著廣泛的應用，原因在于多因素混雜分析是醫學研究中經常遇到的問題。我們已知的絕大多數疾病都是多種原因共同作用的結果，其預后也往往是由多因素所決定的[8]，如本文示例中影響空腹血糖的因素不止上述幾種，還可能有年齡、飲食習慣、工作環境和家族史等，在此基礎上，我們在這些因素中還要明確不同因素的影響權重和可能存在的邏輯關系。這些因素可能影響著臨床試驗的各個方面，造成我們可能需要面臨在難以保證基線相同的混雜情況下對不同干預方法進行比較這樣的難題，以上這些問題很多都可利用多元線性回歸分析來處理。我們控制其混雜因素的一個簡單而有效的方法就是將具體觀測值通過函數的映射關系代入方程進行分析，但這就決定了 MLRA 計算出的預測區間不能說明實際的效果，僅可用于佐證之需。
目前，多元線性回歸方程假設檢驗應用較為普遍的軟件是 SPSS 軟件、SAS 軟件和 Stata 軟件等。SPSS 軟件雖然簡便，但是受其模塊化界面操作的限制，難以完成多種統計，SAS、Stata 軟件雖然具有編程等強大功能，但其高昂的租賃年費也常令我們的工作面臨困境。R 語言因其邏輯運算功能強大，統計制圖多樣，編程語言簡明科學等優勢，始終受到科研工作者的青睞，但 R 語言傳統自帶的環境操作起來較為不便，R Studio 作為 R 語言的高性能操作平臺。它兼容了 R 語言數據處理、計算、制圖以及通過全球不同地區 Comprehensive R Archive Network（CRAN）鏡像站下載升級數據包以隨時不斷豐富數學計算函數模塊等幾乎全部功能外，簡化了編程的難度，提升了操作的便捷性，且它還具有更加優良的數據可視化功能，更具有支持純 R 腳本、腳本文檔混排（rmarkdown）、腳本文檔混排成書（bookdown）、交互式網絡應用（shiny）等優勢。其唯一的障礙即要求我們掌握 R 語言的編寫，才能完成所需要的統計分析和圖形繪制。故本文結合實例提供讀者代碼的具體編寫方法，希望 R Studio 軟件今后能為更多研究者所采用，為醫學研究的發展做出貢獻。
