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

Python07

回归分析 | 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软件的meta分析数据包实现相关系数的Meta分析,第二大部分如何用R语言进行多变量的meta分析。

想获取R语言相关系数meta分析的程序模板的同学请在公众号(全哥的学习生涯)内回复“相关系数”即可。

meta数据包提供实现相关系数的Meta分析命令是:metacor(),这个命令通过加权的倒方差法运用相关系数和纳入的样本数来实现相关系数的随机效用模型和固定效用模型的合并,得到合并的相关系数及95%可信区间。具体的命令如下:

metacor(cor, n,studlab, data= NULL, subset=NULL, sm=.settings$smcor)

cor为每一个纳入研究的相关系数, n为样本量, studlab纳入研究的标签向量, data为相应的的数据集,sm选项为合并的方法,包括ZCOR和COR,其中ZCOR是合并之前先做Fisher Z变换,COR是直接合并。具体的步骤如下:

library(meta)

data<-read.csv(“C:/Users/86187/Desktop/data.csv”)

录入的数据见图1。

data<-metacor(r,n,data=m1,sm="ZCOR")

在这里合并的方法用的是Fisher Z变换。对样本的相关系数做Fisher Z变换是因为Fisher Z变换可以使样本的相关系数的分布正态分布,尤其是在样本量较小的时候,这样便于进一步估计。一般来说,不管是随机还是固定效应都会先对相关系数做Fisher Z变换。只有很少的情况下才直接用相关系数直接来做分析,比如样本量很大的时候,如果直接合并相关系数,当相关系数值接近1的时候,小样本量研究得到的权重会非常大。因此在这里推荐合并的方法都用(ZCOR)Fisher Z变换。Meta分析的结果见图2。

结果显示,异质性检验Q=6.16, P=0.0461, I2=67.5,可以认为有统计学意义上的异质性。选用随机效用模型,COR=0.8427, 95%CI: 0.6264-0.9385, z=4.8724, P<0.0001, 有统计学差异。

具体的命令如下:

forest(a)

从森林图中,非常简单和直观地看到Meta分析的统计结果,见图3

关于这两个方法的介绍请看我之前公众号(全哥的学习生涯)的推送文章(如何用R语言进行meta分析,详细教程一)的内容。敏感性分析和剪补法的结果图分别见图4和图5。

通常Meta分析假定效应量来自于独立的研究,因此统计结果也是独立的。然而,许多研究不能满足独立性的假设,比如多个治疗组与一个共同的对照组比较的研究和多个结局变量的研究就可能产生效应量之间的相关。多变量meta 分析(multivariate meta⁃analysis)作为单变量meta分析的一个拓展,可合并估计多个研究的多个相关参数,这些参数可以是多个结局或多组间的比较。当同一总体中的测量结局相关时,分别对每个结局进行Meta 分析,测量结局之间的相关结构就可能被忽略。多变量Meta分析在随机对照研究中有多种应用,最简单的是在临床试验中把每个组的结局分别处理,其他的应用还有同时探索两个临床结局的治疗效应,或同时探索成本效益的治疗效应,比较多个治疗的联合试验,以及在观察性研究中评估暴露量与疾病之间的相关性,还有在诊断试验和网络干预中的应用。

本次数据来源请见文末的参考文献,主要研究肝硬化的非手术治疗方式预防其出血的危险性,以初次出血的例数为指标,其中三个组分别是:β⁃受体阻滞剂(A),硬化疗法(B),对照组(C),目的是评价这三种非手术治疗方式预防肝硬化出血的效果。,Bled表示初次出血的例数,Total表示干预组的总例数。YAC和YBC分别表示A、B两组相对于C组估计的ln(OR),即干预组的肝硬化初次出血的危险性是对照组的倍数的自然对数;SAA、SBB和SAB则表示其对应方差及两者之间的协方差。对于包含0的研究(研究10和研究20),在每个组增加0.5个初次出血的例数。整理后见表1。

随后安装调用程序包,并进行加载:

install.packages(‘mvmeta’)

library(mvmeta)。

随后将肝硬化初次出血整理后的数据集data(至少包含YAC、YBC、SAA、SAB、SBB变量)保存为csv格式,然后利用下面命令将其导入R语言。

mvmeta 的语句:mvmeta(formula,S,data,subset,method=“reml”,bscov=“unstr”,model=TRUE,contrasts=NULL,offset,na.action,control=list())

其中formula 表示结局变量名称(即YAC、YBC);S 表示研究内(协)方差(即SAA、SAB、SBB);data 表示数据集名称;method 表示所用的估计方法:固定效应模型时选择FIXED;随机效应模型时则选择

限制性最大似然估计(REML)、最大似然估计(ML)、矩估计(MM)、方差成分法(VC)的其中之一,默认为REML。由输出结果中Q 检验的P 值和I2 统计量来判断异质性以及选择何种效应模型。

mvmeta包中主要提供了多变量Meta分析与多变量的Meta 回归,另外也提供了单变量的Meta 分析和Meta 回归。但对于后两者,在R 语言中的metafor、meta、rmeta 及metalik 等包提供了更多、更详尽和有效的功能。多变量Meta 程序为library(mvmeta),调用mvmeta软件包。

model<-mvmeta(cbind(Ya,Yb),S=S,data=cirrhosis)

model <- mvmeta(cbind(Ya,Yb)~X,S=S,data=cirrhosis),此处X代表协变量。

model<-mvmeta(Y,S=S,data=cirrhosis),此处Y为单变量的效应量,S为效应量方差。

model<-mvmeta(Y~X,S=S,data=cirrhosis),此处X代表协变量。

运行以上程序后,最后将结果输出。

单变量和多变量Meta分析都是采用ln(OR)值做分析。单变量Meta分析时YAC和YBC的Q检验P 值均小于0.05,I2统计量分别为57.7%和77.8%。多变量Meta分析Q检验P<0.05,I2统计量为73.9%。可知两种Meta 分析均存在异质性,都用随机效应模型。估计方法选择默认的REML法。

表2 是单变量Meta 分析结果,可得:AC 与BC的OR 值及95%可信区间分别为0.5281(0.2802,0.9955)、0.5406(0.3095,0.9443),表明初次出血的危险性由于干预而降低,即β⁃受体阻滞剂、硬化疗法可以预防肝硬化出血,两者为保护因素。

多变量Meta 分析的结果:YAC 为-0.6755(-1.3073,-0.0438),YBC 为-0.5938(-1.1444,-0.043 2),研究间相关系数为0.436 5(见表3),A组与B组的治疗效果呈正相关。OR 值及95%可信区间分别为0.508 9(0.2705,0.9571)、0.5522(0.318 4,0.957 7),多变量Meta 分析的结果说明β⁃受体阻滞剂预防肝硬化出血的效果是最好,其次是硬化疗法。OR 值的95%可信区间不包含1,上下限均小于1,说明两种疗法与对照组比较的初次出血危险性均小于1,差异有统计学意义。

最后,如果屏幕前的你对R语言学习还有什么问题或者看法,可以在我的公众号(全哥的学习生涯)给我留言,公众号里也有我的个人联系方式,我也希望可以结合更多志同道合的伙伴。

感谢你的阅读。