6.常用统计分析方法——多重线性回归

多重线性回归

多重线性回归模型(Multiple Linear Regression Model) 即描述因变量Y如何随多个自变量X的改变而改变。假定对 n 例观察对象逐一测定了因变量Y与 m 个自变量 X 1 、 X 2 、…、 X m 的数值，多重线性回归模型的一般形式为：

式中 β 0 为常数项，又称截距。 β 1 、 β 2 ，…， β m 称为偏回归系数(partial regression coefficient)或简称回归系数。

e 则是去除 m 个自变量对Y的影响后的 随机误差(即残差) 。

一、多重线性回归的条件：

线性
独立
正态
等方差

二、多重线性回归分析一般可分为两个步骤：

1. 拟合多重线性回归方程

使用lm(因变量~自变量1+自变量2+自变量3...+自变量n，data=mydata)的形成拟合线性回归方程。

2. 模型检验及评价

复相关系数

R= $\sqrt{R^{2}}$ ：称为复相关系数(multiple correlation coefficient)，可用来度量因变量Y与多个自变量间的线性相关程度，亦即观察值Y与估计值之间的相关程度

如果只有一个自变量时， $R=\left | r \right |$ ， r 是简单相关系数。

案例：ISwR包中自带数据集cystfibr是25例囊性纤维化患者的基本信息和肺功能信息。其中包含的变量有性别、年龄、身高、体重，BMI，用力呼气量（fev1），残气量（rv），功能性残气量（frc），最大肺活量（tlc），最大呼气压力（pemax）。

library(ISwR)
data(cystfibr)
str(cystfibr)#查看数据结构类型
#分类变量因子化
cystfibr$sex <- factor(cystfibr$sex,labels = c("male","female"))

因为结局指标肺功能的参数较多，现选择其中fev1作为结局的分析指标。

cor(cystfibr[,6:10])

可以发现肺功能指标之间存在较高的相关性，因此选择其中一个作为分析。

fit1 <- lm(fev1 ~ age + sex + height + weight + bmp, data = cystfibr)
#线性拟合
summary(fit1)

多重线性的分析结果可以F统计量=3.992，p=0.012<0.05，说明回归方程的假设是有统计学差异的。在回归系数的假设检验中，只有性别的回归系数有统计学差异。female表示女性与男性相比，fev1低10.04。

自变量的筛选和最优模型：

一、自变量筛选的统计学标准：

1.残差平方和（ $SS_{残差}$ ）缩小或决定系数 $R^{2}$ 增大 若某一自变量被引入模型后使残差平方和缩小很多或决定系数增大很多，说明该变量对Y的作用很大，可以被引入；反之，说明其中对y的作用很小，不被引入。

2. 残差的均方（MS残差）缩小或校正决定系数 $R^2_{ad}$ 增大 。虽然自变量的增加可使分子位置的残差平方和减小，但同时分母位置上的自由度（n-p-1）也将减小，如果增加的自变量无统计学意义，分子中残差平方和的减少不多，而分母中自由度却减小很多，从而残差的均方并不减少。要使残差的均方缩小，必须是引入模型中的自变量对y的影响较大，使得残差平方和减少很多。

3. $C_{p}$ 统计量缩小

$SS_{残差,p}$ 表示p个自变量做回归时的残差平方和， $MS_{残差,全}$ 表示用全部自变量作回归时的残差均方。标准：选择 $C_{p}$ 较小，并且 $C_{p}$ 的值接近未知参数个数的模型作为较优模型。使 $C_{p}$ 统计量较小的标准符合了统计建模的原则：一方面入选的自变量不能太多，以使方程易于理解和解释，并且降低工作量及研究经费；另一方面自变量不能太少，以保证一定的估计和预测精度。

4. 赤池信息量准则AIC

AIC越小，说明模型用较少的参数获得了足够的拟合优度。

二、多重共线性

多重共线性是指：自变量之间由于存在相关性，而使得模型估计失真或难以估计准确。

方差膨胀因子（VIF） ：表征自变量之间的共线性程度。

容忍度：是方差膨胀因子的倒数，如果容忍度小于0.25，则可能存在共线性问题；如果容忍度小于0.10，则提示共线性问题严重。

案例中的height（身高）、体重和BMI之间是相关的。查看相关系数

cor(cystfibr[,3:5])

变量age对应的方差膨胀因子VIF

#比如age对应的VIF
lm.age <- summary(lm(age~sex+height+weight+bmp,data=cystfibr))
#计算VIF
VIF.age <- 1/(1-lm.age$r.squared);VIF.age

当VIF>5时，此时 $R_j^{2}$ >80%时可以认为存在多重共线性，此时用最小二乘法计算的回归系数是不准确的。

car包中的vif（）函数可以快速计算所有变量对应的VIF值

library(car)
vif(fit1)

解决多重共线性问题的方法 ：

（1）剔除造成共线性的自变量，重新建立模型；

（2）采取主成分分析法，将一组具有共线性的自变量合并成少数不相关的变量；

（3）采用逐步回归法，限制有较强相关性的自变量同时进入方程

自变量筛选常用的算法：

前进法：自变量由少到多，一个一个引入方程。将部分回归平方和最大且能使F检验拒绝H0者入选为第一个自变量；规定一个界值 $F_{引入}$ ，接着将余下的变量中部分回归平方和最大，并使F检验拒绝H0者入选为第二个自变量，如此不断引入新的自变量，直到不能拒绝H0为止。
后退法：自变量由多到少，一个一个从方程中剔除。首先对全部候选变量进行总的回归，每次剔除一个部分回归平方和最小，而使F检验不能拒绝H0者，直到不能再剔除时为止。
逐步法：将上述两种法结合起来，在向前引入的每一步之后，都要考虑从已引入方程的变量中剔除相形见绌者。事先规定两个界值 $F_{引入}$ 和 $F_{剔除}$ （ $F_{引入}$ > $F_{剔除}$ )，当方程外侯选变量中部分回归平方和最大者F值大于等于 $F_{引入}$ 时，引入相应的变量；当方程内变量中部分回归平方和最小者F值小于等于 $F_{剔除}$ 时，剔除相应的变量。如此引入和剔除交替进行，直到方程外无可引入的变量，而方程也无可剔除的变量时为止。
最优子集回归法：对于p个自变量的线性回归问题，所有可能的自变量子集作回归方程，共有（ $2^{p}-1$ ）个。根据某种变量的选择准则，通过比较各子集符合准则的程度，从中选择出一个或几个最优的回归子集。仅适合于自变量个数不多的情况。

函数drop1( ) 可以从包含所有自变量的模型中手动一次剔除一个变量以进行逐步回归。。依据AIC准则，在较少的参数情况下实现足够的拟合度

drop1(fit1)

结果说明：在不踢除的任何自变量，模型AIC值为113.81，剔除age自变量，模型AIC为111.83，以此类推。比较所有的AIC值可知，需要剔除变量age。

另外可以通过设置drop1函数中参数test得到似然比检验或F检验的结果，

drop1(fit1,test = "F")

根据结果可以看到提出sex后，p有统计学差异，其余变量剔除没有统计学差异，age对应的p值最大，所以可以考虑先剔除age，再次使用drop1()判断哪一个变量需要从模型中剔除

函数step（） 可以根据AIC值进行逐步回归自动选择最优模型，其中参数direction可以设置为默认的both（逐步法）、forward（向前法）、backward（向后法）

fit2 <- step(fit1)

最后生成只包含sex和bmp的指标时模型AIC最小。

summary(fit2)

没有剔除之前的决定系数为0.5123，只保留sex和bmp的决定系数为0.507，总体的决定系数相差不大，但调整后的决定系数明显提高，回归系数的检验p都具有统计学差异。

函数confint（） 计算回归系数的置信区间：

confint(fit2)

回归模型的假设检验：

上述的t检验只表明回归系数与0的差异有无统计学意义，要想检验变量对于模型有无统计学意义，需要进行F检验的结论，这在 自变量是多个水平的分类变量时尤为必要 。

anova(fit2)

三、哑变量设置

多重线性回归分析中的自变量可以是连续型的变量，也可以是二分类变量，但不能把无序多分类变量直接纳入模型。必须先将无序多分类变量变换成为多个二分类变量之后，才能将她们引入回归模型。若某个无序多分类变量有k各可能的类别，则可用（k-1）个二分类变量而取而代之，这个过程称为“哑元化”，得到的（k-1）个二分类变量称为“哑变量”。

四、交互效应

具体判断有无交互效应可以参考简单线性回归中的分层线性回归。

五、回归诊断

使用gvlma包中的 gvlma（）函数 对模型进行综合检测

#install.packages("gvlma")
library(gvlma)
gvlma(fit2)

global stat及其他检验都大于0.05，说明满足回归条件。如果有p>0.05，则用简单线性回归中提到的残差分析几种方式来判断哪些假设没有被满足。

第二种方法：残差图

res <- fit2$residuals
hist(res)#样本量较小难以评价
qqnorm(res)#QQ图
qqline(res)
shapiro.test(res)#正态性检验
plot(fitted(fit2),res,xlab = "Fitted values",pch=15,col="blue") #绘制残差图
abline(h=0,col="red")

最终的模型可以用epiDisplay包中的regress.display()汇总输出
regress.display(fit2)
write.csv(regress.display(fit2)$table,file = "线性回归练习模型输出版.csv")
				什么是回归分析？
回归分析是研究自变量与因变量之间数量变化关系的一种分析方法，它主要是通过因变量Y与影响它的自变量Xi（i1,2,3…）之间的回归模型，衡量自变量Xi对因变量Y的影响能力的，进而可以用来预测因变量Y的发展趋势。
????????????相关分析与回归分析的联系与区别：
两者均为研究与测度两个或两个以上变量之间关系的方法。
在实际工作中，一般先进行相关分析，计算相关函数，然后建立回归模型，最后用回归模型进行推算或预测。
✨相关分析研究的都是随机变量，并且不分自变量与因变量；
				现实中大多数的结果都是多个原因引起的，所谓一因一果对应的事情是非常少见的。从医学领域来看，几乎任何疾病都不是由一种原因造成的，而是多钟病因共同作用的结果。因此，多重线性回归在实际中应用更为广泛。
多重线性回归（multiple linear regression），有的教材也称之为多元线性回归，就是指一个因变量，多个自变量的线性回归。多重线性回归比单因素的线性回归复杂之处不在于多了几个变量，更为
				多重线性回归分析简介：
简单线性回归是只考虑单因素影响的预测模型，事实上，影响因变量的因素往往不止一个，可能会有多个影响因素，也就是研究一个因变量与多个自变量的线性回归问题，就用到多重线性回归分析了。
多重线性回归的作用：
1、探索对于因变量具有影响作用的因素；
2、控制混杂因素，评价多个自变量对因变量的独立效应；
3、用已知的自变量来估计和预测因变量的值及其变化。
多重线性回归与多元线性回归
多重线性回归：是指包含两个或两个以上自变量的线性回归模型；
多元线性回归：是指包含两个或两个以上因变量的线性回
				P133，这是第二次作业，考察多重线性回归。这个youtube频道真是精品，用R做统计。这里是R代码的总结。
连续变量和类别型变量总要分开讨论；
多重线性回归可以写成矩阵形式的一元一次回归；相当于把多变量当成列向量压缩一下；矩阵有着非常优美的简单的计算法则，大大简化了计算的复杂度；
在建多重线性回归模型时我们会遇到很多问题：
选哪些变量建模，一元的很简单，可以判断有无显著性，多元就复杂...
				一元线性回归是SPSS中常用的统计分析方法之一。它用于研究自变量与因变量之间的线性关系，并通过回归方程描述这种关系。下面以某医院的病人年龄和其住院天数为例进行一元线性回归分析。
首先，我们收集了100位病人的数据，其中自变量是病人的年龄，因变量是病人的住院天数。我们将这些数据输入SPSS软件进行分析。
在SPSS软件中，首先选择"回归"，然后选择"线性"，将因变量（住院天数）拖放到"因变量"栏中，将自变量（年龄）拖放到"解释变量"栏中。
在回归分析结果中，我们关注回归方程的系数、显著性和决定系数R^2。
回归方程的系数是关键指标之一，它告诉我们自变量（年龄）对因变量（住院天数）的影响程度。如果系数为正数，则表示自变量的增加与因变量的增加有正相关关系；如果系数为负数，则表示自变量的增加与因变量的减少有负相关关系。通过系数的大小，我们可以判断自变量对因变量的影响强弱。
显著性水平是判断系数是否统计上显著的指标。在假设检验中，一般认为当p值小于0.05时，系数是显著的，即存在影响关系；当p值大于0.05时，系数是不显著的，即没有影响关系。
决定系数R^2可以解释回归方程的拟合程度，其取值范围在0到1之间。R^2的值越接近1，说明回归方程可以更好地解释因变量的变异；R^2的值越接近0，说明回归方程对因变量的解释能力较弱。
通过分析回归方程的系数、显著性和决定系数R^2，我们可以得出结论：病人的年龄对住院天数有显著影响，且年龄与住院天数呈正相关关系。这些可以为医院制定合理的住院时间和医疗方案提供科学依据。
                    松子小鱼: 
                    # 建立 logistic 回归模型使用全部数据
> alb = openxlsx::read.xlsx("Validate  data.xlsx", sheet =3, rowNames = T)
> ddist <- datadist(alb)
> options(datadist='ddist')
> f <- lrm(Group ~ PON1 + HPR + THBS1 + Age, data = alb, x = TRUE, y = TRUE)  
Error in eval(predvars, data, callenv) : object 'Group' not found
这个报错什么意思啊
                7. 预测模型评价——DCA曲线
                    2401_84451762: 
                    你好，上传的附件在哪
                实用技巧——绘制森林图
                    Jing 我愿。: 
                    这个se计算有误，多了负号而且忘记除2了？ 赞同，而且这段代码过渡的时候，似乎是修改了列名result2$characteristic <- ifelse(is.na(result2$number),
                                 result2$characteristic,
                                 paste0("    ",result2$characteristic))