直線回歸是研究一個應(yīng)變量與單個自變量之間呈直線關(guān)系的一種統(tǒng)計方法。但是在醫(yī)學(xué)研究中，許多疾病都有多種原因，而且預(yù)后也是多種因素決定的，如糖尿病患者的血糖變化可能受胰島素、糖化血紅蛋白、血清總膽固醇、甘油三酯等多種生化指標(biāo)的影響，這個時候就需要研究多個因素對一個結(jié)果變量數(shù)量依存的線性關(guān)系以及多個因素之間的線性關(guān)系，這種關(guān)系就稱為多元線性回歸。例如如人體的體表現(xiàn)積隨身高和體重而改變，則體表面積為應(yīng)變量，身高和體重為自變量，設(shè)體表現(xiàn)為Y，身高和體重為X1和X2，則所建立的含有兩個自變量二元回歸方程為：

其中左側(cè)的y（帶帽子的字母）是y均數(shù)的估計值或預(yù)測值（predicted value，y hat），a為截距，表示各自變量取0時，y的估計值，b1和b2稱為偏回歸系數(shù)（partial regression coefficient），其中b1表示在體重的取值固定不變時，身高每增加一個單位，y估計值平均變化b1個單位；b2表示在身高的取值不變時，體重每增加一個單位，y估計值平均變化b2個單位。一般來講，對于擬合一個應(yīng)變量y和m個自變量x1，x2,..xm的多元回歸方程，其形式為：

其中，a為截距，意義同上，bi為當(dāng)其他自變量取值固定不變時，xi每增加一個單位，y的估計值平均變化bi個單位。

多元線性回歸資料應(yīng)滿足以下條件：

1. 自變量和應(yīng)變量之間的關(guān)系是線性關(guān)系；

2. 各觀測單位相互獨立；

3. 殘差服從正態(tài)分布；

4. 殘差滿足方差齊性。

多元線性回歸的原理

對于多元回歸方程的建立，也就是求出各常數(shù)項和各偏回歸系數(shù)，原理還是使用最小二乘法來求各待估計系數(shù)，即根據(jù)觀察到的n例數(shù)據(jù)，使殘差平方和最小，即如下所示：

由此可以得到下面的方程組，求解得到各個偏回歸系數(shù)和截矩，如下所示：

b0的計算如下所示：

這兩個公式分別表示自變量的離均差平方和（i=j），兩個自變量的離均差積和及自變量Xj與應(yīng)變量Y的離均差積和。

多元線性回歸的計算過程

現(xiàn)在看一個案例：

27名糖尿病患者的血清總膽固醇、甘油三酯、空腹胰島素、糖化血紅蛋白、空腹血糖的測量值位于下表中，試建立血糖與其他幾項指標(biāo)的多元線性回歸方程（《醫(yī)學(xué)統(tǒng)計學(xué)》第四版.孫振球，第15章，計算過程采用了湯銀才的《R語言與統(tǒng)計分析》中的過程）。

第一步，由上表的數(shù)據(jù)，通過公式（二十一）和公式（二十二）計算出包括應(yīng)變量在內(nèi)的各變量離差矩陣，如下所示：

根據(jù)公式（十九）列出方程組求解后，可以得到：

b1=0.1424, b2=0.3515,  b3=-0.2706, b4=0.6382

算得均數(shù)分別為：

根據(jù)公式（二十）可以計算出截矩，如下所示：

有了這幾個系數(shù)與截矩，我們就可以得到線性回歸方程，如下所示：

在R中計算多元線性回歸方程與簡單線性回歸基本上類似，也是使用lm()函數(shù)，如下所示：

計算結(jié)果如下所示：

從結(jié)果中可以看出來，截矩，X1，X2，X3，X4的偏回歸系數(shù)分別為5.9433，0.1424，0.3515，-0.2706，0.6382，對回歸的結(jié)果進(jìn)行t檢驗可以看出來，回歸方程的系數(shù)的顯著性不高, 有的沒有通過檢驗，這說明如果選擇全部變量構(gòu)造方程, 效果并不好，這就涉及到變量選擇的問題，以建立“最優(yōu)”的回歸方程。如果先不考慮變量的篩選問題，那么我們計算得到的多元線性回歸方程就是下面的這個樣子，與前面手工計算的結(jié)果一樣：

多元線性回歸中變量的篩選

前面提到，有的變量沒有通過t檢驗，例如X1，因此我們要對變量進(jìn)行篩選，在R中，進(jìn)行“最優(yōu)”回歸方程的方法是“逐步回歸法”的計算函數(shù)step( ),它是以Akaike信息統(tǒng)計量為準(zhǔn)則(簡稱AIC準(zhǔn)則), 通過選擇最小的AIC信息統(tǒng)計量, 來達(dá)到刪除或增加變量的目的。AIC的數(shù)值越小，就說明對變量的篩選效果越好，它的用法如下所示：

其中，object是線性模型或廣義線性模型的結(jié)果，scope是確定逐步回歸的搜索區(qū)域，direction確定逐步搜索的方向，both是一切子集回歸法，backward是向后法，forward是向前法，默認(rèn)值為both，現(xiàn)在使用setp()對上述的多元線性回歸結(jié)果做逐步回歸，如下所示：

計算結(jié)果如下所示：

結(jié)論: 用全部變量作回歸方程時, AIC統(tǒng)計量的值為42.16,如果去掉變量X1，AIC統(tǒng)計量的值為40.34，如下去掉變量X2，AIC統(tǒng)計量的值為43.568，依次類推。由于去掉X1命名AIC統(tǒng)計量達(dá)到最小，因此R會自動掉變量X1，進(jìn)入下一輪的計算，在下一輪中，無論去掉哪個變量，AIC的統(tǒng)計量的值均會升高，因此R會自動終止計算，得到“最優(yōu)”回歸方程，現(xiàn)在使用summary()提取相關(guān)回歸信息，如下所示：

結(jié)論：從上面的結(jié)果我們可以看出來，回歸系數(shù)的顯著性水平有很大提高，每個自變量的檢驗均是顯著的，??得到了“最優(yōu)”的回歸方程，如下所示：

回歸診斷

由樣本數(shù)據(jù)得到回歸方程后，為了確定回歸方程及引入的自變量是否有統(tǒng)計學(xué)意義，必須進(jìn)一步做假設(shè)檢驗。方差分析法可以將回歸方程中所有自變量作為一個整體來檢驗它們與應(yīng)變量Y與之間是否具有線性關(guān)系，并對回歸方程的預(yù)測或解釋能力作出綜合評價。在此基礎(chǔ)上進(jìn)一步對各變量的重要性作出評價。

方差分析

方差分析法原理

在多元線性回歸中，為什么要做方差分析呢，這是因為在線性回歸中，把因變量的總變異，也就是因變量的所有觀測值到其均值的距離平方和（通常用SST）表示，分解成對應(yīng)于各個自變量的若干平方和(SS,sum of squares)，每一個這樣的平方和代表了一個自變量對總變異平方的貢獻(xiàn)，剩下不能被自變量解釋的，就是誤差平方和或殘差平方和(SSE,error sum of squares)，自變量及殘差平方和所對應(yīng)的平方和加起來等于因變量的總變異。

用簡單的話來講，就是方差分析表把由自變量所導(dǎo)致的因?的變化和誤差所導(dǎo)致的因變量的變化找出來，并進(jìn)行比較（通過F檢驗），如果這兩個變化很不一樣，那么則認(rèn)為自變量在回歸中的確顯著影響了因變量，表示回歸有道理，否則，沒有理由認(rèn)為因變量和誤差有所不同。

具體來講，如果我們有兩個自變量，分別為A和B，則把SSt分解為對應(yīng)于這兩個自變量的兩個平方和SSA，SSB以及殘差平方和SSE，即SST=SSA+SSB+SSE。那么如果SSA與SSE相比較很顯著，即變量A的貢獻(xiàn)顯著大于誤差的貢獻(xiàn)，則稱變量A顯著，否則，如果SSA和SSE的貢獻(xiàn)差不多，那么變量A就不顯著。這種顯著需要正式的檢驗，如果前面關(guān)于誤差的四個假定成立（回歸模型的誤差項獨立，且服從正態(tài)分布），那么這些平方和（比如這里的SSA、SSB、SSE）都有卡方分布（因為它們是正態(tài)分布變量的平方和），而除以各自的自由度之后，得到它們的平均（叫均方，mean square，記為MSA，MSB，MSE），根據(jù)F分布的定義，MSA/MSE和MSB/MSE都有F分布，如果MSA/MSE很大，那么A就是和殘差相比顯著，則相應(yīng)的p值就小。

方差分析法過程

方程的假設(shè)檢驗如下所示：

將應(yīng)變量Y的總變異分解成兩部分，即如下所示：

其中左側(cè)部分是回歸平方和，最右側(cè)部分是殘差平方和，上式可記作：

SS總=SS回+SS殘

其中，回歸平方和可用下式計算：

殘差平方和為：SS殘=SS總-SS回

用統(tǒng)計量F檢驗假設(shè)H0是否成立，如下所示：

方差分析見下表：

如果F>F(α,(回歸,殘差))，則在α水準(zhǔn)上拒絕H0，接受H1，認(rèn)為應(yīng)變量Y與m個自變量之x1,x2,xm間存在線性回歸關(guān)系?，F(xiàn)在我們從剛才的多元線性回歸方程中計算各部分的變異：

SS總=222.5519

SS回=0.1424×67.6962+0.3515×89.8025+0.2706×142.4337+0.6382×84.5570=133.7107

SS殘=222.5519-133.7107

方差分析的結(jié)果如下表所示：

查F界值表可知，F(xiàn)(0.01,(3,23))=4.76，F(xiàn)>4.76，P<>

注，在R中查找相應(yīng)F值的函數(shù)為qf()，用法如下所示：

使用我要查找在0.01水平上，兩個自由度為3和23的F值，如下所示：

在R中，可以通過anova()來查看多元線性回歸的方差分析結(jié)果，如下所示：

決定系數(shù)R平方

根據(jù)方差分析的結(jié)果，還可以得到多元線性回歸的決定系數(shù)R平方，它的計算公式如下所示：

從上面的結(jié)果可以看出來R的平方=0.5981，這個數(shù)字從前面的方差分析結(jié)果中也能看出來，即Multiple R-squared:  0.5981。這個數(shù)字表明，血糖含量變量的58%可由甘油三脂(X2)、胰島素(X3)、糖化血紅蛋白(X4)的變化來解釋。

復(fù)相關(guān)系數(shù)

R稱為復(fù)相關(guān)系數(shù)（multiple correlation coefficient），可用來度量應(yīng)變量Y與多個自變量的線性相關(guān)程度，亦即觀察值Y與估計值之間的相關(guān)程度，在這個案例中，復(fù)相關(guān)系數(shù)R=0.7734（0.5981開根號），如果只有一個自變量時，R=|r|，r為簡單相關(guān)系數(shù)。

t檢驗

方差分析和決定系數(shù)是將所有自變量作為一個整體來檢驗和說明它們與Y的相關(guān)程度及解釋能力的，并未指明方程中的每一個自變量對Y的影響如何，但是在醫(yī)學(xué)研究中，我們往往更加關(guān)心是對各自變量的解釋。對每一自變量的作用進(jìn)行檢驗和稱量它們對Y的作用大小，可以采用偏回歸平方和法，t檢驗法，標(biāo)準(zhǔn)化回歸系數(shù)法。

而在R中的，通常采用的是t檢驗法，在前面的多元線性回歸中，我們就可以看到這個t檢驗法的結(jié)果，如下所示：

根據(jù)t界值表，我們可以發(fā)現(xiàn)，剔除了X1變量后，剩下的變量的p值都小于0.05（在R中，一個星號表示大于0.01，小于0.05）,說明這些變量有統(tǒng)計學(xué)意義，對于同一資料而言，不同自變量的t值間可以相互比較，t的絕對值越大，說明該自變量對Y的回歸所起的作用越大。

殘差分析

殘差分析是檢查資料是否符合模型條件的一個有用的工具。如果資料與模型的假設(shè)偏離較大，最小二乘估計、假設(shè)檢驗及區(qū)間估計則有可能出現(xiàn)問題。繪制多元線性回歸的殘差圖，如下所示：

共線性診斷

共線性問題是指擬合多元線性回歸時，自變量之間存在線性關(guān)系或氛這線性關(guān)系，自變量之間的線性關(guān)系將會隱藏變量的顯著性，增加參數(shù)估計的誤差，還會產(chǎn)生一個很不穩(wěn)定的模型，因此共線性診斷就是要找出哪些變量之間存在共線性關(guān)系，共線性診斷的方法主要有特征值法，條件指數(shù)法，方差膨脹因子法。具體的原理可以參考相應(yīng)的統(tǒng)計學(xué)書籍，例如《R語言與統(tǒng)計分析》（湯銀才），這里只說一種方法，即方差膨脹因子法，使用的函數(shù)為car包中的vif()，VIF的全稱是Variance Inflation Factor，即方差膨脹因子，一般原則下，如果VIF值大于4就表明存在多重共線性問題，它的使用如下所示：

從結(jié)果我們可以看出，這幾個變量的VIF值都沒有大于4，因此不存在多重共線性問題。

參考資料

1. 孫振球. 醫(yī)學(xué)統(tǒng)計學(xué).第4版[M]. 人民衛(wèi)生出版社, 2014.

2. DawnGriffiths. 深入淺出統(tǒng)計學(xué)[M]. 電子工業(yè)出版社, 2012.

3. 王炳順. 醫(yī)學(xué)統(tǒng)計學(xué)及SAS應(yīng)用[J]. 2007.

4. 王斌會. 多元統(tǒng)計分析及R語言建模[M]. 暨南大學(xué)出版社, 2010.

5. 湯銀才. R語言與統(tǒng)計分析[M]. 高等教育出版社, 2008.

6. 卡巴科弗. R語言實戰(zhàn)[M]. 人民郵電出版社, 2016.