生物统计(6)-多元线性回归

直线回归是研究一个应变量与单个自变量之间呈直线关系的一种统计方法。但是在医学研究中，许多疾病都有多种原因，而且预后也是多种因素决定的，如糖尿病患者的血糖变化可能受胰岛素、糖化血红蛋白、血清总胆固醇、甘油三酯等多种生化指标的影响，这个时候就需要研究多个因素对一个结果变量数量依存的线性关系以及多个因素之间的线性关系，这种关系就称为多元线性回归。例如如人体的体表现积随身高和体重而改变，则体表面积为应变量，身高和体重为自变量，设体表现为Y，身高和体重为X1和X2，则所建立的含有两个自变量二元回归方程为：

其中左侧的y（带帽子的字母）是y均数的估计值或预测值（predicted value，y hat），a为截距，表示各自变量取0时，y的估计值，b1和b2称为偏回归系数（partial regression coefficient），其中b1表示在体重的取值固定不变时，身高每增加一个单位，y估计值平均变化b1个单位；b2表示在身高的取值不变时，体重每增加一个单位，y估计值平均变化b2个单位。一般来讲，对于拟合一个应变量y和m个自变量x1，x2,..xm的多元回归方程，其形式为：

其中，a为截距，意义同上，bi为当其他自变量取值固定不变时，xi每增加一个单位，y的估计值平均变化bi个单位。

多元线性回归资料应满足以下条件：

1. 自变量和应变量之间的关系是线性关系；

2. 各观测单位相互独立；

3. 残差服从正态分布；

4. 残差满足方差齐性。

多元线性回归的原理

对于多元回归方程的建立，也就是求出各常数项和各偏回归系数，原理还是使用最小二乘法来求各待估计系数，即根据观察到的n例数据，使残差平方和最小，即如下所示：

由此可以得到下面的方程组，求解得到各个偏回归系数和截矩，如下所示：

b0的计算如下所示：

这两个公式分别表示自变量的离均差平方和（i=j），两个自变量的离均差积和及自变量Xj与应变量Y的离均差积和。

多元线性回归的计算过程

现在看一个案例：

27名糖尿病患者的血清总胆固醇、甘油三酯、空腹胰岛素、糖化血红蛋白、空腹血糖的测量值位于下表中，试建立血糖与其他几项指标的多元线性回归方程（《医学统计学》第四版.孙振球，第15章，计算过程采用了汤银才的《R语言与统计分析》中的过程）。

第一步，由上表的数据，通过公式（二十一）和公式（二十二）计算出包括应变量在内的各变量离差矩阵，如下所示：

根据公式（十九）列出方程组求解后，可以得到：

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

算得均数分别为：

根据公式（二十）可以计算出截矩，如下所示：

有了这几个系数与截矩，我们就可以得到线性回归方程，如下所示：

在R中计算多元线性回归方程与简单线性回归基本上类似，也是使用lm()函数，如下所示：

mlr <>'https://raw./20170505a/raw_data/master/data_szq_1501.csv',header=T,sep=',') # 原始数据集
mlr_result <>
summary(mlr_result)


计算结果如下所示：

> summary(mlr_result)

Call:
lm(formula = mlr$Y ~ mlr$X1 + mlr$X2 + mlr$X3 + mlr$X4)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6268 -1.2004 -0.2276  1.5389  4.4467 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   5.9433     2.8286   2.101   0.0473 *
mlr$X1        0.1424     0.3657   0.390   0.7006  
mlr$X2        0.3515     0.2042   1.721   0.0993 .
mlr$X3       -0.2706     0.1214  -2.229   0.0363 *
mlr$X4        0.6382     0.2433   2.623   0.0155 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.01 on 22 degrees of freedom
Multiple R-squared:  0.6008,    Adjusted R-squared:  0.5282 
F-statistic: 8.278 on 4 and 22 DF,  p-value: 0.0003121


从结果中可以看出来，截矩，X1，X2，X3，X4的偏回归系数分别为5.9433，0.1424，0.3515，-0.2706，0.6382，对回归的结果进行t检验可以看出来，回归方程的系数的显著性不高, 有的没有通过检验，这说明如果选择全部变量构造方程, 效果并不好，这就涉及到变量选择的问题，以建立“最优”的回归方程。如果先不考虑变量的筛选问题，那么我们计算得到的多元线性回归方程就是下面的这个样子，与前面手工计算的结果一样：

多元线性回归中变量的筛选

前面提到，有的变量没有通过t检验，例如X1，因此我们要对变量进行筛选，在R中，进行“最优”回归方程的方法是“逐步回归法”的计算函数step( ),它是以Akaike信息统计量为准则(简称AIC准则), 通过选择最小的AIC信息统计量, 来达到删除或增加变量的目的。AIC的数值越小，就说明对变量的筛选效果越好，它的用法如下所示：

step(object, scope, scale = 0,
     direction = c('both', 'backward', 'forward'),
     trace = 1, keep = NULL, steps = 1000, k = 2, ...)


其中，object是线性模型或广义线性模型的结果，scope是确定逐步回归的搜索区域，direction确定逐步搜索的方向，both是一切子集回归法，backward是向后法，forward是向前法，默认值为both，现在使用setp()对上述的多元线性回归结果做逐步回归，如下所示：

lm.step <>


计算结果如下所示：

> lm.step <>
Start:  AIC=42.16
mlr$Y ~ mlr$X1 + mlr$X2 + mlr$X3 + mlr$X4

         Df Sum of Sq     RSS    AIC
- mlr$X1  1    0.6129  89.454 40.343
                 88.841 42.157
- mlr$X2  1   11.9627 100.804 43.568
- mlr$X3  1   20.0635 108.905 45.655
- mlr$X4  1   27.7939 116.635 47.507

Step:  AIC=40.34
mlr$Y ~ mlr$X2 + mlr$X3 + mlr$X4

         Df Sum of Sq     RSS    AIC
                 89.454 40.343
- mlr$X3  1    25.690 115.144 45.159
- mlr$X2  1    26.530 115.984 45.356
- mlr$X4  1    32.269 121.723 46.660


结论: 用全部变量作回归方程时, AIC统计量的值为42.16,如果去掉变量X1，AIC统计量的值为40.34，如下去掉变量X2，AIC统计量的值为43.568，依次类推。由于去掉X1命名AIC统计量达到最小，因此R会自动掉变量X1，进入下一轮的计算，在下一轮中，无论去掉哪个变量，AIC的统计量的值均会升高，因此R会自动终止计算，得到“最优”回归方程，现在使用summary()提取相关回归信息，如下所示：

summary(lm.step)
# 
# Call:
#   lm(formula = mlr$Y ~ mlr$X2 + mlr$X3 + mlr$X4)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -3.2692 -1.2305 -0.2023  1.4886  4.6570 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)   
# (Intercept)   6.4996     2.3962   2.713  0.01242 * 
#   mlr$X2        0.4023     0.1541   2.612  0.01559 * 
#   mlr$X3       -0.2870     0.1117  -2.570  0.01712 * 
#   mlr$X4        0.6632     0.2303   2.880  0.00845 **
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.972 on 23 degrees of freedom
# Multiple R-squared:  0.5981,    Adjusted R-squared:  0.5456 
# F-statistic: 11.41 on 3 and 23 DF,  p-value: 8.793e-05


结论：从上面的结果我们可以看出来，回归系数的显著性水平有很大提高，每个自变量的检验均是显著的，𢅟得到了“最优”的回归方程，如下所示：

回归诊断

由样本数据得到回归方程后，为了确定回归方程及引入的自变量是否有统计学意义，必须进一步做假设检验。方差分析法可以将回归方程中所有自变量作为一个整体来检验它们与应变量Y与之间是否具有线性关系，并对回归方程的预测或解释能力作出综合评价。在此基础上进一步对各变量的重要性作出评价。

方差分析

方差分析法原理

在多元线性回归中，为什么要做方差分析呢，这是因为在线性回归中，把因变量的总变异，也就是因变量的所有观测值到其均值的距离平方和（通常用SST）表示，分解成对应于各个自变量的若干平方和(SS,sum of squares)，每一个这样的平方和代表了一个自变量对总变异平方的贡献，剩下不能被自变量解释的，就是误差平方和或残差平方和(SSE,error sum of squares)，自变量及残差平方和所对应的平方和加起来等于因变量的总变异。

用简单的话来讲，就是方差分析表把由自变量所导致的因㸄的变化和误差所导致的因变量的变化找出来，并进行比较（通过F检验），如果这两个变化很不一样，那么则认为自变量在回归中的确显著影响了因变量，表示回归有道理，否则，没有理由认为因变量和误差有所不同。

具体来讲，如果我们有两个自变量，分别为A和B，则把SSt分解为对应于这两个自变量的两个平方和SSA，SSB以及残差平方和SSE，即SST=SSA+SSB+SSE。那么如果SSA与SSE相比较很显著，即变量A的贡献显著大于误差的贡献，则称变量A显著，否则，如果SSA和SSE的贡献差不多，那么变量A就不显著。这种显著需要正式的检验，如果前面关于误差的四个假定成立（回归模型的误差项独立，且服从正态分布），那么这些平方和（比如这里的SSA、SSB、SSE）都有卡方分布（因为它们是正态分布变量的平方和），而除以各自的自由度之后，得到它们的平均（叫均方，mean square，记为MSA，MSB，MSE），根据F分布的定义，MSA/MSE和MSB/MSE都有F分布，如果MSA/MSE很大，那么A就是和残差相比显著，则相应的p值就小。

方差分析法过程

方程的假设检验如下所示：

将应变量Y的总变异分解成两部分，即如下所示：

其中左侧部分是回归平方和，最右侧部分是残差平方和，上式可记作：

SS总=SS回+SS残

其中，回归平方和可用下式计算：

残差平方和为：SS残=SS总-SS回

用统计量F检验假设H0是否成立，如下所示：

方差分析见下表：

变异来源	自由度	SS	MS	F
总变异	n-1	SS总
回归	m	SS回	SS回/m	MS回/MS残
残差	n-m-1	SS残	SS残/(n-m-1)

如果F>F(α,(回归,残差))，则在α水准上拒绝H0，接受H1，认为应变量Y与m个自变量之x1,x2,xm间存在线性回归关系。现在我们从刚才的多元线性回归方程中计算各部分的变异：

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

方差分析的结果如下表所示：

变异来源	自由度	SS	MS	F	P
总变异	26	222.5519
回归	3	133.0978	44.36593	11.4	<>
残差	23	89.4541	3.89

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

注，在R中查找相应F值的函数为qf()，用法如下所示：

qf(p, df1, df2, ncp, lower.tail = TRUE, log.p = FALSE)


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

qf(0.99,3,23)
# [1] 4.764877


在R中，可以通过anova()来查看多元线性回归的方差分析结果，如下所示：

mlr_result2 <>
# 使用x2,x3,x4变量做多元线性回归

summary(mlr_result2)
# 查看计算结果
# Call:
#   lm(formula = mlr$Y ~ mlr$X2 + mlr$X3 + mlr$X4)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -3.2692 -1.2305 -0.2023  1.4886  4.6570 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)   
# (Intercept)   6.4996     2.3962   2.713  0.01242 * 
#   mlr$X2        0.4023     0.1541   2.612  0.01559 * 
#   mlr$X3       -0.2870     0.1117  -2.570  0.01712 * 
#   mlr$X4        0.6632     0.2303   2.880  0.00845 **
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.972 on 23 degrees of freedom
# Multiple R-squared:  0.5981,    Adjusted R-squared:  0.5456 
# F-statistic: 11.41 on 3 and 23 DF,  p-value: 8.793e-05

aov_result <>
aov_result
# Analysis of Variance Table
# 
# Response: mlr$Y
# Df Sum Sq Mean Sq F value   Pr(>F)   
# mlr$X2     1 46.787  46.787 12.0297 0.002082 **
#   mlr$X3     1 54.042  54.042 13.8950 0.001104 **
#   mlr$X4     1 32.269  32.269  8.2968 0.008447 **
#   Residuals 23 89.454   3.889                    
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
sum(aov_result$`Sum Sq`)
# SS总
# [1] 222.5519
sum(aov_result$`Sum Sq`)-aov_result$`Sum Sq`[4]
# SS回归
# [1] 133.0978
222.5519-133.0978
# SS残
89.4541


决定系数R平方

根据方差分析的结果，还可以得到多元线性回归的决定系数R平方，它的计算公式如下所示：

从上面的结果可以看出来R的平方=0.5981，这个数字从前面的方差分析结果中也能看出来，即Multiple R-squared:  0.5981。这个数字表明，血糖含量变量的58%可由甘油三脂(X2)、胰岛素(X3)、糖化血红蛋白(X4)的变化来解释。

复相关系数

R称为复相关系数（multiple correlation coefficient），可用来度量应变量Y与多个自变量的线性相关程度，亦即观察值Y与估计值之间的相关程度，在这个案例中，复相关系数R=0.7734（0.5981开根号），如果只有一个自变量时，R=|r|，r为简单相关系数。

t检验

方差分析和决定系数是将所有自变量作为一个整体来检验和说明它们与Y的相关程度及解释能力的，并未指明方程中的每一个自变量对Y的影响如何，但是在医学研究中，我们往往更加关心是对各自变量的解释。对每一自变量的作用进行检验和称量它们对Y的作用大小，可以采用偏回归平方和法，t检验法，标准化回归系数法。

而在R中的，通常采用的是t检验法，在前面的多元线性回归中，我们就可以看到这个t检验法的结果，如下所示：

# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)   
# (Intercept)   6.4996     2.3962   2.713  0.01242 * 
#   mlr$X2        0.4023     0.1541   2.612  0.01559 * 
#   mlr$X3       -0.2870     0.1117  -2.570  0.01712 * 
#   mlr$X4        0.6632     0.2303   2.880  0.00845 **


根据t界值表，我们可以发现，剔除了X1变量后，剩下的变量的p值都小于0.05（在R中，一个星号表示大于0.01，小于0.05）,说明这些变量有统计学意义，对于同一资料而言，不同自变量的t值间可以相互比较，t的绝对值越大，说明该自变量对Y的回归所起的作用越大。

残差分析

残差分析是检查资料是否符合模型条件的一个有用的工具。如果资料与模型的假设偏离较大，最小二乘估计、假设检验及区间估计则有可能出现问题。绘制多元线性回归的残差图，如下所示：

par(mfrow=c(2,2))
plot(mlr_result2)


共线性诊断

共线性问题是指拟合多元线性回归时，自变量之间存在线性关系或氛这线性关系，自变量之间的线性关系将会隐藏变量的显著性，增加参数估计的误差，还会产生一个很不稳定的模型，因此共线性诊断就是要找出哪些变量之间存在共线性关系，共线性诊断的方法主要有特征值法，条件指数法，方差膨胀因子法。具体的原理可以参考相应的统计学书籍，例如《R语言与统计分析》（汤银才），这里只说一种方法，即方差膨胀因子法，使用的函数为car包中的vif()，VIF的全称是Variance Inflation Factor，即方差膨胀因子，一般原则下，如果VIF值大于4就表明存在多重共线性问题，它的使用如下所示：

library(car)
vif(mlr_result2)
# mlr$X2   mlr$X3   mlr$X4 
# 1.051763 1.123518 1.178342 


从结果我们可以看出，这几个变量的VIF值都没有大于4，因此不存在多重共线性问题。

参考资料

1. 孙振球. 医学统计学.第4版[M]. 人民卫生出版社, 2014.

2. DawnGriffiths. 深入浅出统计学[M]. 电子工业出版社, 2012.

3. 王炳顺. 医学统计学及SAS应用[J]. 2007.

4. 王斌会. 多元统计分析及R语言建模[M]. 暨南大学出版社, 2010.

5. 汤银才. R语言与统计分析[M]. 高等教育出版社, 2008.

6. 卡巴科弗. R语言实战[M]. 人民邮电出版社, 2016.