回归分析 | R语言 -- 多元线性回归

2023-02-24 09:51:01Python028

回归分析 | R语言 -- 多元线性回归,第1张

多元线性回归 是简单线性回归的扩展，用于基于多个不同的预测变量（x）预测结果变量（y）。

例如，对于三个预测变量（x），y的预测由以下等式表示： y = b0 + b1*x1 + b2*x2 + b3*x3

回归贝塔系数测量每个预测变量与结果之间的关联。“ b_j”可以解释为“ x_j”每增加一个单位对y的平均影响，同时保持所有其他预测变量不变。

在本节中，依然使用 datarium 包中的 marketing 数据集，我们将建立一个多元回归模型，根据在三种广告媒体（youtube，facebook和报纸）上投入的预算来预测销售。计算公式如下： sales = b0 + b1*youtube + b2*facebook + b3*newspaper

您可以如下计算R中的多个回归模型系数：

请注意，如果您的数据中包含许多预测变量，则可以使用 ~. 以下命令将模型中的所有可用变量简单地包括在内：

从上面的输出中，系数表显示β系数估计值及其显着性水平。列为：

如前所述，您可以使用R函数轻松进行预测 predict() ：

在使用模型进行预测之前，您需要评估模型的统计显着性。通过显示模型的统计摘要，可以轻松地进行检查。

显示模型的统计摘要，如下所示：

摘要输出显示6个组件，包括：

解释多元回归分析的第一步是在模型摘要的底部检查F统计量和关联的p值。

在我们的示例中，可以看出F统计量的p值<2.2e-16，这是非常重要的。这意味着 至少一个预测变量与结果变量显着相关 。

要查看哪些预测变量很重要，您可以检查系数表，该表显示了回归beta系数和相关的t统计p值的估计。

对于给定的预测变量，t统计量评估预测变量和结果变量之间是否存在显着关联，即，预测变量的beta系数是否显着不同于零。

可以看出，youtube和facebook广告预算的变化与销售的变化显着相关，而报纸预算的变化与销售却没有显着相关。

对于给定的预测变量，系数（b）可以解释为预测变量增加一个单位，同时保持所有其他预测变量固定的对y的平均影响。

例如，对于固定数量的youtube和报纸广告预算，在Facebook广告上花费额外的1000美元，平均可以使销售额增加大约0.1885 * 1000 = 189个销售单位。

youtube系数表明，在所有其他预测变量保持不变的情况下，youtube广告预算每增加1000美元，我们平均可以预期增加0.045 * 1000 = 45个销售单位。

我们发现报纸在多元回归模型中并不重要。这意味着，对于固定数量的youtube和报纸广告预算，报纸广告预算的变化不会显着影响销售单位。

由于报纸变量不重要，因此可以 将其从模型中删除 ，以提高模型精度：

最后，我们的模型公式可以写成如下：。 sales = 3.43+ 0.045*youtube + 0.187*facebook

一旦确定至少一个预测变量与结果显着相关，就应该通过检查模型对数据的拟合程度来继续诊断。此过程也称为拟合优度

可以使用以下三个数量来评估线性回归拟合的整体质量，这些数量显示在模型摘要中：

与预测误差相对应的RSE（或模型 sigma ）大致代表模型观察到的结果值和预测值之间的平均差。RSE越低，模型就越适合我们的数据。

将RSE除以结果变量的平均值将为您提供预测误差率，该误差率应尽可能小。

在我们的示例中，仅使用youtube和facebook预测变量，RSE = 2.11，这意味着观察到的销售值与预测值的平均偏差约为2.11个单位。

这对应于2.11 / mean（train.data $ sales）= 2.11 / 16.77 = 13％的错误率，这很低。

R平方（R2）的范围是0到1，代表结果变量中的变化比例，可以用模型预测变量来解释。

对于简单的线性回归，R2是结果与预测变量之间的皮尔森相关系数的平方。在多元线性回归中，R2表示观察到的结果值与预测值之间的相关系数。

R2衡量模型拟合数据的程度。R2越高，模型越好。然而，R2的一个问题是，即使将更多变量添加到模型中，R2总是会增加，即使这些变量与结果之间的关联性很小（James等，2014）。解决方案是通过考虑预测变量的数量来调整R2。

摘要输出中“已调整的R平方”值中的调整是对预测模型中包含的x变量数量的校正。

因此，您应该主要考虑调整后的R平方，对于更多数量的预测变量，它是受罚的R2。

在我们的示例中，调整后的R2为0.88，这很好。

回想一下，F统计量给出了模型的整体重要性。它评估至少一个预测变量是否具有非零系数。

在简单的线性回归中，此检验并不是真正有趣的事情，因为它只是复制了系数表中可用的t检验给出的信息。

一旦我们开始在多元线性回归中使用多个预测变量，F统计量就变得更加重要。

大的F统计量将对应于统计上显着的p值（p <0.05）。在我们的示例中，F统计量644产生的p值为1.46e-42，这是非常重要的。

我们将使用测试数据进行预测，以评估回归模型的性能。

步骤如下：

从上面的输出中，R2为 0.9281111 ，这意味着观察到的结果值与预测的结果值高度相关，这非常好。

预测误差RMSE为 1.612069 ，表示误差率为 1.612069 / mean(testData $ sales) = 1.612069/ 15.567 = 10.35 ％，这很好。

本章介绍了线性回归的基础，并提供了R中用于计算简单和多个线性回归模型的实例。我们还描述了如何评估模型的性能以进行预测。

首先，你可以思考哪些算法可能适用于你的数据集。

其次，我建议尽可能地尝试混合算法并观察哪个方法最适用于你的数据集。

尝试混合算法(如事件模型和树模型)

尝试混合不同的学习算法(如处理相同类型数据的不同算法)

尝试混合不同类型的模型(如线性和非线性函数或者参数和非参数模型)！

今天被粉丝发的文章给难住了，又偷偷去学习了一下竞争风险模型，想起之前写的关于竞争风险模型的做法，真的都是皮毛哟，大家见笑了。想着就顺便把所有的生存分析的知识和R语言的做法和论文报告方法都给大家梳理一遍。

什么时候用生存分析

当你关心结局和结局发生时间的时候，就要考虑生存分析了，这种既有结局又有时间的数据叫做生存数据，英文叫做Time-to-event data. 只不过因为这个方法医学上用来分析存活情况用的多，所以得名生存分析，反正你就记住一个例子，我要研究汽车发生故障，我也应该用生存分析，因为我既关心是不是有故障，我还关心用了多久（跑了多远）才出故障，就是既有time，又有event，Time-to-event data就用生存分析。

基本概念

首先是删失，对象失访了，脱落了，出现结局之前随访结束了，都叫做删失：

删失又分为左删失，区间删失和右删失，图示如下：

比如我想研究得了A病的人的生存情况，存在的所有可能情形为：

第一种，研究的开始的时候有人已经有A病，这个时候人家已经活了一段时间了，具体多久我不知道，叫做左删失；

第二种，入组随访的时候没病，中途得了A病死了，什么时候得的，没记录下来，叫区间删失；

第三种，得了A病，一直活到了研究结束还没死，叫做右删失。

你看，所有的删失情况造成的后果都是我们没法准确估计发生结局的时间，这也是其名字删失的由来，对于这类数据就需记录为删失数据。

生存分析的种类有哪些

具体的种类是为了回答具体的问题，我们做生存分析常常要回答的问题如下：

一个是描述生存情况，一个是比较，再一个就是探究影响因素。

比如我随访了很多病人，我就想知道随着时间变化这群人的生存概率是如何变化的（描述）？我就可以用简单粗暴的Kaplan-Meier method，又叫乘积极限法，基本思想就是此刻的生存概率等于上一时刻的生存概率乘以此刻的存活率。

比如我手上有如下数据：

time是随访时间，status是结局，我就可以写出如下代码拟合出总体人群的生存曲线：

fit1 <- survfit(Surv(time, status) ~ 1, data = mydata)summary（fit1）

并且得到每个时间点，病人生存概率的估计值和标准误：

如果我还恰好收集了病人的性别，我又想看看男病人和女病人生存情况是不是有不同（比较），我可以对生存曲线分组比较，代码如下：

fit2<- survfit(Surv(time, status) ~ sex, data = mydata)summary（fit2）

输出两组生存曲线比较的log-rank test的统计量：

survdiff(Surv(time, status) ~ sex, data = mydata)

并且我们还可以进行复杂分组的组间比较，大家可以翻看我之前的文章。

但是大家明白，KM法始终是在做单因素分析，而且都是做的分类变量的单因素分析，我们经常还会有的需求是考虑各种混杂的情况下去探讨影响生存时间的因素，这个时候我们就要用到The Cox regression model了。

模型形式如下：

上面的式子把h0移项到左边，等号两边同时取对数就成了一个线性模型，和广义线性模型的理解一样一样的，b就是我们的影响因素的系数，解释为lnHR的改变量，其为正就是危险因素，为负就是保护因素：

The quantities exp(bi) are called hazard ratios (HR). A value of bi greater than zero, or equivalently a hazard ratio greater than one, indicates that as the value of the ith covariate increases, the event hazard increases and thus the length of survival decreases.

再观测一下上面的式子，我们拟合了预测因素对风险比例（t时刻风险比上基础风险）的模型，这个时候暗含的假设就是就是对于每个人在任何时刻风险只有一个常数，就是在所有预测因素的作用下，各个时间的风险是不变的，这个就叫做Proportional Hazards Assumption, 比例风险假设。

做COX比例风险模型的时候都有必要验证这个假设是不是满足，具体方法如下：

我们需要先做一个cox模型，拟合cox模型需要用到coxph函数，代码如下：

res.cox <- coxph(Surv(time, status) ~ age + sex + wt.loss, data = lung)res.cox

模型输出结果里面会有预测变量的β值（coef）和其标准误，有HR（exp(coef)）值，还有预测变量的显著性检验结果：

我们将这个模型对象直接喂给cox.zph，就可以得到风险比例假设的检验结果了：

test.ph <- cox.zph(res.cox)test.ph

可以看到，我们的p值都大于0.05，即都满足ph假设。

我们还可以从图上看，看scaled Schoenfeld residuals是不是和时间独立，如果独立则满足ph假设：

ggcoxzph(test.ph)

有竞争事件时

上面写了只有一个截距事件时生存分析的不同目的，描述，比较，和影响因素探讨，接着再来看存在竞争事件的时候该如何描述，比较，和探讨影响因素。

生存分析的另外情况就是竞争风险模型，就是time to event的event有多种，一种发生了另外一种就不可能发生了，这个就是竞争，比如好多文献中都会列举的经典例子：就是在造血干细胞移植人群中，我们关心其疾病复发风险，但是有些人因为移植死了（TRM），死了之后肯定也谈不上复发，如果你把因为移植死了的人都作为删失数据，肯定也是不对的，这个里面就会有两个竞争结局相关性造成的问题，此时我们应该用竞争风险模型来分析。

For example, disease relapse is an event of interest in studies of allogeneic hematopoietic stem cell transplantation (HSCT), as is mortality related to complications of transplantation (transplant-related mortality or TRM). Relapse and TRM are not independent in this setting because these two events are likely related to immunologic effector mechanisms following HSCT, whereby efforts to reduce TRM may adversely affect the risk of relapsemoreover, patients who die from TRM cannot be at further risk of relapse.

在比例风险中我们的结局事件只有一个，我们关心的是哪些因素对事件发生的风险有影响。在竞争风险模型中我们关注的地方又变了，因为我们有好几个结局事件，这个时候我们会关心在竞争事件存在的情况下各个结局事件的累计发病函数是如何随着时间变化的，以及如何来比较不同的累计发病函数，以及如何探讨影响因素（competing risks regression analysis）。

之前写的在非竞争风险模型中累计发病率的组间比较可以用KM法，将纵坐标换一换，用log-rank检验，在竞争风险模型中我们需要用Fine and Gray提到的方法来做（Gray’s test），如果比如说我发现两组（治疗组和对照组）的累计发病风险不一样，我肯定还想探讨哪些因素影响累计发病风险，之前是用COX比例风险模型做的，在竞争结局存在的情况下我们依然是得用Fine and Gray提出的模型（Fine and Gray Model），叫做竞争风险回归模型：

Fine and Gray (6) proposed a method for direct regression modeling of the effect of covariates on the cumulative incidence function for competing risks data. As in any other regression analysis, modeling cumulative incidence functions for competing risks can be used to identify potential prognostic factors for a particular failure in the presence of competing risks, or to assess a prognostic factor of interest after adjusting for other potential risk factors in the model.

首先看竞争风险时候累计发病曲线的比较方法。我手上有数据如下：

其中dis是疾病类型，2分类，ftime是时间，status是结局事件，结局事件有3个水平，多的1个水平是竞争事件。现在我关心不同疾病人群各个事件累积发病曲线有无不同，我可以用cuminc函数结合plot画出各个组的累计发病曲线：

fit=cuminc (ftime, status, dis, cencode = 0) plot(fit)

fit对象的结果中还有每条曲线的比较，从比较结果可以看出两组在事件2的累积发病曲线上是有显著差异的：

上面介绍的方法相当于没有竞争风险的时候的KM法，通过上面的方法我们知道有了不同风险结局累计发病曲线的差异，继续我们会继续看影响因素，要做的就是竞争风险回归模型了，需要用到的函数就是crr。

比如我手上有如下数据，除了时间，结局还包括每个病人的像sex，age等等协变量，我想探讨说这些因素是如何影响病人某个结局的，我就可以写出一个竞争风险回归模型：

mod1 <- crr(ftime,Status,x)summary(mod1)

上面的代码中x是自变量矩阵，运行代码输出竞争风险回归模型的结果如下：

到这儿基本就写完了所有的生存分析的情况，接着再结合一篇论文看看报告方法，论文就是下面这篇，是我随意查的：

S. Chen, H. Sun, M. Heng, X. Tong, P. Geldsetzer, Z. Wang, P. Wu, J. Yang, Y. Hu, C.

Wang, T. Bärnighausen, Factors Predicting Progression to Severe COVID-19: A Competing Risk Survival Analysis of 1753 Patients in Community Isolation in Wuhan, China, Engineering (2021), doi: https://doi.org/10.1016/j.eng.2021.07.021

这个作者用竞争风险模型探讨了重型新冠进程的影响因素，作者报告了每个影响因素的HR和置信区间，p值，这些都很容易做出来。还报告了固定时间点的累积发病的置信区间。在R语言中固定时间点的累计发病置信区间可以用CumIncidence这个函数计算出来：

The CumIncidence() function allows for the pointwise confidence intervals, by simply adding a further argument, level, where we specify the desired confidence level.

比如我想计算第10天各组的累计发病风险的置信区间，我就可以写出如下代码：

CumIncidence (ftime, status, dis, cencode = 0,t=10,level = 0.95)

输出结果如下：

图中就是各个组在时间为10的时候结局累计发病风险的置信区间。

小结

也欢迎大家的意见和建议，大家想了解什么统计方法都可以在文章下留言，说不定我看见了就会给你写教程哦，有疑问欢迎私信。

如果你是一个大学本科生或研究生，如果你正在因为你的统计作业、数据分析、模型构建,科研统计设计等发愁，如果你在使用SPSS, R，Python，Mplus, Excel中遇到任何问题，都可以联系我。因为我可以给您提供最好的，最详细和耐心的数据分析服务。

如果你对Z检验，t检验，方差分析，多元方差分析，回归，卡方检验，相关，多水平模型，结构方程模型，中介调节，量表信效度等等统计技巧有任何问题，请私信我，获取详细和耐心的指导。

If you are a student and you are worried about you statistical #Assignments, #Data #Analysis, #Thesis, #Reports, #Composing, #Quizzes, Exams.. And if you are facing problem in #SPSS, #R-Programming, #Excel, Mplus, then contact me. Because I could provide you the best services for your Data Analysis.

Are you confused with statistical Techniques like z-test, t-test, ANOVA, MANOVA, Regression, Logistic Regression, Chi-Square, Correlation, Association, SEM, multilevel model, mediation and moderation etc. for your Data Analysis...??

Then Contact Me. I will solve your Problem...

加油吧，打工人！

R数据分析：用R语言做潜类别分析LCA

R数据分析：ROC曲线与模型评价实例

R数据分析：用R语言做meta分析

R数据分析：贝叶斯定理的R语言模拟

R数据分析：什么是人群归因分数Population Attributable Fraction

R数据分析：有调节的中介

R数据分析：如何用R做多重插补，实例操练

R数据分析：如何用R做验证性因子分析及画图，实例操练

R数据分析：手把手教你画列线图(Nomogram)及解读结果

R数据分析：倾向性评分匹配完整实例（R实现）