统计与R入门 II 回归分析

内容:

  • 回归
  • 零假说显著性检验
  • 中央极限定理
  • 抽样分布
  • 一般线性模型
  • 方差分析
  • 调节
  • 中介
  • 路径模型

回归(regression)

回归:用一个或多个预测变量(predictor)来预测结果变量(outcome variable)值的统计分析

回归分析分为:

  1. 简单回归:使用一个预测变量
  2. 多元回归:使用多个预测变量

简单回归的公式为:Y=B0+B1X1+e,其中:

  • Y是X1的线性函数,是真实值
  • m是回归常数(regression constant)(或称截距intercept)
  • b是回归系数(regression coefficient)(或称斜率slope)
  • e是误差(error)(或称残差residual)

实际使用时的回归模型是:Ŷ=B0+B1X1,其中:

  • Ŷ是预测值

评价回归模型,有如下两个较为重要的量:

  1. R多元回归系数(multiple correlation coefficient):即预测值与观测值之间的相关系数rŶY
  2. R2:Y变量的偏差能被回归模型所能解释的程度

R和R2都是用来评价模型的总体表现的

多元回归的公式为:

  • Y=B0+B1X1+B2X2+...+BKXK+e

R示例

使用longley数据集,数据集包含1947年至1962年16年的包含7个变量的经济数据。

unnamed-chunk-1

尝试构建一个用就业人口预测国民生产总值的简单回归模型

用公式计算非标准化的回归系数、标准化的回归系数,以及运用coef函数和lm()函数计算回归系数

回归常数为-1430.48

回归系数为27.84

计算多元回归系数R与R2

R中可以用summary()来方便地获得这些量

回归系数的计算

选择回归系数的主要思想是:

  • 选择能够使回归模型预测结果最佳的回归系数。

回归模型预测结果最佳,意味着残差(预测误差)最小。

残差e=Ŷ-Y,可能为正也可能为负,求平方和处理。

计算残差的平方和SSRESIDUAL=Σ(Ŷ-Y)2,选择能够使这项最小的回归系数

另一种考虑方式是,残差是模型所不能解释的偏差情况,用韦恩图加以解释:

QQ截图20140101143840

其中:

  • SSX是变量X偏差的平方和SSX=Σ(X-MX)2
  • SSY是变量Y偏差的平方和SSY=Σ(Y-MY)2
  • SSMODEL=SP是外积和即模型所能解释的偏差的平方和,SP= Σ[(X-MX)×(Ŷ-MŶ)]
  • SSRESIDUAL是模型所不能解释的残差的平方和SSRESIDUAL=Σ(Ŷ-Y)2

选择能够使模型所不能解释的残差平方和最小的回归系数

非标准化的回归系数计算公式为:

  • B1=r×(SDY/SDX)

r是皮尔森积差相关系数,是变量Y随着变量X变化的程度,将其乘以Y与X的标准差之比,是考虑到Y与X的不同取值范围。

标准化的回归系数计算公式为:β=r

其原因是:经过标准化处理(均转换为Z值)的X与Y的标准差都为1。

线性回归的前提假设

线性回归的前提假设与相关性分析的前提假设基本一致:

  1. Y是正态分布
  2. X与Y之间是线性关系
  3. 方差齐性
  4. 变量X与Y的可靠性
  5. 变量X与Y的有效性
  6. 变量X与Y是否来自于随机抽样生成的具有代表性的样本

唯一的区别是,X不一定是正态分布的。

为了测试这些假设前提,通常可以绘制预测变量与残差的图表

回顾Anscombe’s quartet

unnamed-chunk-5

计算残差,可以用公式计算,也可以调用predict()函数或者residuals()函数。

绘制图表

unnamed-chunk-7

只有左上方的满足方差齐性,残差与X无关,是随机的。

其他三个数据集都不满足线性回归分析的前提假设。

检验一下longley数据集中构建的线性回归模型

unnamed-chunk-8

应该是满足方差齐性的。

零假说显著性检验(Null Hypothesis Significance Testing,NHST)

设定两个不同的假说:

  • H0:零假说(Null Hypothesis)
  • HA:对立假说(Alternative Hypothesis)

例如,相关性分析中,想证明两个变量之间不是无关的,运用零假说显著性检验,则:

  • 零假说H0是:r=0
  • 对立假说HA是:r>0

例如,线性回归分析:

  • 零假说H0是:B=0
  • 对立假说HA是:B!=0

如果对立假说预测了X与Y之间的方向性,则称为:定向检验(directional test)或单尾检验(single tail test)=

否则称为无方向性检验(non-directional test),或双尾检验(two tail test)

考虑如下线性回归分析的显著性检验设定:

  • 零假说H0是:B=0
  • 对立假说HA是:B!=0

假设H0为真,计算获得我们所拥有的数据的条件概率:

  • p = P(D|H0)

如果发现该概率p值非常小,则驳回零假说,否则保留H0

QQ截图20140101165234

运用零假说显著性检验测试的可能结果如上图所示:

  • 如果H0为真,显著性检验测试结果保留H0,则检验结果是正确的
  • 如果H0为假,显著性检验测试结果拒绝H0,则检验结果是正确的
  • 如果H0为真,显著性检验测试结果拒绝H0,则是1型错误或称假警报(false alarm)
  • 如果H0为假,显著性检验测试结果保留H0,则是2型错误或称遗漏(miss)

对p值的解读:

  • 正确的解读是:如果零假说是正确的,那么我们获得手头上数据或者更加极端数据的概率是p,即P(D|H0)
  • 错误的解读是:零假说正确的概率是p,即即P(H0|D)

为了获得p值,需要首先计算出t值,公式为:

t=B/SE,其中:

  • B为非标准化的回归系数
  • SE是的标准误差(standard error)是抽样统计量的标准差
    对于回归系数而言,计算公式为:SE=SQRT[SSRESIDUAL/(DF)/SSX]

t值是一个比例!

回归系数B是我们观测到的两个变量之间的线性相关程度

标准误差SE,是随机误差的情况

t值这个比例是:观测到的情况/随机情况,

如果t值为1,则说明观测到的情况与完全随机的情况是一样的

明显的t值越大,p值越小

中央极限定理部分再仔细解释t值。

R演示

计算以就业人数Employed预测国民生产总值GNP的回归系数的t值

此时SE计算公式为:SE=SQRT[(SSRESIDUAL/(N-2))/SSX]

零假说显著性检验的问题和补救方法

零假说显著性检验的若干问题:

  1. 受到样本量的偏倚:样本量越大,t值越大,p值越小,越有可能拒绝零假说
  2. 随机选择的阀值α:即使是“标准”的阀值0.05的选择也是随机的
  3. 只知道NHST:有些情况下有比NHST更合适的显著性检验
  4. 容易产生误差:如果对同一个数据集做多个NHST容易产生1型错误,很多领域获得的数据的抽样误差较大,NHST容易产生2型错误
  5. 有问题的逻辑:NHST的逻辑是,如果假说成立,则不太可能获得当前数据(p->~q)。现在我们有这样的数据,因此假说不成立(q->~p)。

零假说显著性检验的补救方法:

  1. 受到样本量的偏倚:在NHST以外提供效应量(effect size)作为补充
  2. 随机选择的阀值α:提供效应量(effect size)作为补充,并且解读p值时,不按p值与α的距离汇报"高"显著性或"低"显著性
  3. 只知道NHST:学习其他假说检验方法,考虑模型比较(model comparison)
  4. 容易产生误差:重复试验以避免1型错误,获得随机的有代表性的样本以避免2型错误
  5. 有问题的逻辑:不要错误地解读p值

抽样分布(sampling distributions)

如果知道样本数据的分布情况,我们可以做概率上的推理。

例如:知道人的体温是正态分布的,平均体温为36.5摄氏度。

那么随机选一个人测体温,其体温大于等于36.5摄氏度的概率是50%。

将其转换成Z值,即Z>0的概率为50%

体温大于38摄氏度,对应于Z值>2的概率为2%

抽样分布的概念:

  • 通过多个数据量相同的样本所获得的统计量的分布

例如:

  • 平均值的抽样分布
  • 相关性系数的抽样分布
  • 回归系数的抽样分布

通常我们并没有多个样本,而只是估计抽样分布的情况。

  • 假设我们有一个随机从总体中抽样出的样本量为N的样本
  • 对这个样本,我们计算出平均值
  • 假设我们现今有多个这样的随机样本,数据量均为N
  • 这些样本的平均值一起构成了平均值的抽样分布

有了这样的平均值的抽样分布,我们可以回答这样的问题:

如果我们从总体中抽取一个样本,这个样本的平均值小于Z=0的概率是多少。

如果抽样分布是正态分布的,答案是0.5。

R演示

总体为1至50,每次抽取样本量为20的样本,4次实验分别抽取20次,100次,500次和2000次,分别绘制抽样平均值的柱图。

unnamed-chunk-10

可以看出,无论总体分布情况如何,抽样次数越多,点统计量的抽样分布越接近于正态分布。

中央极限定律(central limit theorem)

三条规则:

  • 抽样分布的平均值与总体的平均值相同
  • 抽样分布的标准差是抽样分布方差的平方根,σ22/N
  • 如果N>=30,或总体满足正态分布,则抽样分布的形状近似于正态分布

第三条告诉我们样本量足够大时,抽样分布的形状近似于正态分布。

t值得计算公式为t=B/SE

t值也有一个抽样分布,称为t分布

t分布是一个分布族,不同的样本量大小对应不同的t分布,样本量越大t分布越接近正态分布

样本量越小,t分布越宽,达到同样大小的p值所需要的t值越大。

lecture_slides-Stats1.13.L09

因此,在计算出t值后,根据样本量大小找出对应的t分布,在该分布上根据t值计算出p值。

将p值与0.05比较解读为:t值是否属于t分布中5%的极端情况中。

置信区间(confidence intervals)

样本均值的置信区间

任何一个样本统计量,如均值、标准差,都是点估计量(point estimates)。

即,一个样本平均值,代表的是平均值的抽样分布中的一个点。

置信区间是汇报一个可能值的区间估计(interval estimate),而非一个点估计量。置信区间是:根据随机样本,对总体参数所作出的一个区间估计。

例如:95%置信度意味着有95%的概率,区间中包含总体参数的真实值。

抽样误差将会导致不同的样本会有不同的点估计值

置信区间的优点是,将抽样误差考虑进来了,汇报的是区间估计

置信区间受到两方面的影响:

  • 样本大小
  • 总体和样本的方差大小

回顾平均值的标准误差为SE=SD/SQRT(N)很好解释了,置信区间与标准误差之间的关系。

置信区间的计算方法:

  • 上限M+t*SE
  • 下限M+t*SE

t值取决于自由度(样本量决定)和置信度。

R示例

从国民生产总值中抽取一个样本量为10的样本,计算平均值,并汇报置信区间

回归系数的置信区间

回归系数B也是点估计量,从单一样本中获得的B值,是回归系数的抽样分布中的一个点。

回归系数的置信区间也是考虑进标准误差之后的结果。

置信度95解读为:有%95的概率,总体的回归系数在置信区间中。

R示例

用就业人数预测国民生产总值时回归系数的置信区间,confint()是简便算法。

R示例,在图表上绘制回归系数的置信区间

unnamed-chunk-12

多元回归(multiple regression)

简单回归是只使用一个预测变量,多元回归是使用多个预测变量

公式为:Ŷ=B0+B1X1+B2X2+...+BKXK=Σ(BnXn)

其中:

  • Ŷ是结果变量Y的预测值
  • B0是在所有X均为0时的预测值
  • XK是预测变量
  • BK是非标准化的回归系数
  • Y-Ŷ是残差(预测误差)
  • K是预测变量的数量

同简单回归模型一样,评价多元回归模型,有如下两个较为重要的量:

  1. R多元回归系数(multiple correlation coefficient):即预测值与观测值之间的相关系数rŶY
  2. R2:Y变量的偏差能被回归模型所能解释的程度
  • R和R2都是用来评价模型的总体表现的

R示例

在longley数据集中,利用就业人数和总人口两个变量预测国民生产总值;

通过p值发现,这两个变量均是显著的。比较标准化的回归系数,在这个多元回归模型中,相对于就业人数而言,人口是更强的预测变量。

多元回归的回归系数的估计

同简单回归一样,回归系数的值也是要是的模型的预测误差最小,即使残差的平方和最小。

标准化的多元回归模型公式(矩阵形式)为:Ŷ=B(X)

其中:

  • Ŷ是N×1的向量
  • B是K×1的向量
  • X是N×K的矩阵

为了求解B,用Y替代Ŷ

由交换律有:Y=X(B)

等式两边同时乘以XT有:XT(Y)=XT(XB)

两边再同时乘以(XTX)-1得到:B=(XTX)-1(XTY)

R示例

solve()是矩阵求逆(inverse),需要MASS包裹,%*%是矩阵乘法,t()是矩阵转置(transpose)


 

 

 

我将常数项加在了X的最后一列,将其与lm()的结果比较,结果是一样的,:

一般线性模型(general linear model,GLM)

一般线性模型是在许多常见的统计分析,例如多元回归和方差分析(ANOVA)中采用的数学框架(mathematical framework)

GLM的特点如下:

  • 线性的(linear):变量对之间假设是呈线性关系的
  • 累加的(addictive):如果是用多个变量来预测一个结果变量,则每一个预测变量的效果都被认为是累加的

可以利用GLM来做一系列的检验,例如测试变量之间的非累加性质等。

下面举例说明一般线性模型的几种实例:

简单回归:Y=B0+B1X1+e,其中

  • Y为薪水
  • X1为工作年限

多元回归:Y=B0+B1X1+B2X2+B3X3+e,其中

  • Y为薪水
  • X1为工作年限
  • X2为获奖次数
  • X3为(工作年限*获奖次数)

X3是非累加的,加上这个变量可以用来测试获奖次数是否是调节变量,获奖次数是否能调节工作年限对薪水的影响。

单因素方差分析(one way ANOVA): Y=B0+B1X1+e,其中

  • Y为薪水
  • X1为性别

在这个方差分析中,性别是类别变量,而非连续型变量。

因素方差分析(factorial ANOVA):Y=B0+B1X1+B2X2+B3X3+e,其中

  • Y为薪水
  • X1为性别
  • X2为民族
  • X3为性别*民族

这个方差分析可以分析性别和民族之间的交互作用。

方差分析(Analysis of Variance,ANOVA)

是在因变量为类别变量而结果变量为连续变量时适用的分析方法。

方差分析最常用于有超过两个实验组时的随机实验所获得的数据

如果实验组只有2个,可以用非独立或独立t检验(indenpendent t-test, dependent t-test)

R示例

glm()函数,以npk数据集为例,做因素方差分析

虚拟编码(dummy coding)

虚拟编码是在回归分析中,将类别型的预测变量编码为数值的方法。

例如:

  • 因变量为:学科类别{计算机,统计,心理学}
  • 结果变量为:发文量

虚拟编码时,通常选择一组作为参考组,例如选择计算机作为参考组,上例的虚拟编码结果可能如下所示:

 D1  D2
 计算机  0  0
 统计  1 0
 心理学  0 1

虚拟编码(D1,D2)数量为2,永远是组数-1。

R示例

以iris数据集为例,对Species做虚拟编码处理,用以做回归分析,预测Sepal.Length。演示两种dummy coding方法,C()和factor()

setosa被处理成参考组,截距的值实际为setosa的回归系数,这个值与种类为setosa的Sepal.Length平均值是一样的。

对结果可以这样解读:每从setosa向versicolor有1个单位的改变,Sepal.Length的值将增加0.93。

这种组间差异的显著性可以看对应的p值,可以发现,上例中都是显著的。

另外还有编码方式效果编码(effects coding),分为无权重的效果编码和带权重的效果编码。

调节(moderation)

也称交互作用(interaction)

一个调节变量(moderator variable)Z,Z变量的取值能够影响模型中其他变量之间的关系。

例如,如果Z是X与Y之间的调节变量,则可以发现X与Y之间的关系可以表达呈Z的函数。

以实验研究为例:研究中我们通过调整自变量X,研究因变量Y的变化。

  • 如果存在一个变量Z,我们发现X随着Y的变化情况根据Z值的不同呈现出不一致的情况,则此时发现Z为调节变量。

以相关性分析为例:假设X与Y之间存在相关性

  • 存在调节变量Z意味着:X与Y之间的相关性对应于Z的分布呈现出不一致的情形。
  • 即:在Z值不同的时候,X与Y之间的相关性不同

调节模型示例:

假设X与Z均为连续型变量:

  • Y=B0+B1X+B2Z+B3(X×Z)+e

如果检验B3发现其影响是显著的,则证明Z是起到调节作用的。

假设X为类型变量,Z为连续变量,X有三种不同类型:

  • Y=B0+B1(D1)+B2(D2)+B3(Z)+B4(D1×Z)+B5(D2×Z)+e

调节的检验:

假设X与Z均为连续型变量,构建两个模型:模型1为不含调节

  • Y=B0+B1X+B2Z+e

模型2为含调节:

  • Y=B0+B1X+B2Z+B3(X×Z)+e

假设X为类型变量,Z为连续变量,X有三种不同类型,构建两个模型:

模型1为不含调节:

  • Y = B0+B1(D1)+B2(D2)+B3Z+e

模型2为含调节:

  • Y=B0+B1(D1)+B2(D2)+B3(Z)+B4(D1×Z)+B5(D2×Z)+e

构建完两个模型后:

  • 比较R2的值
  • 评价与调节效应有关的预测变量的回归系数,如(X×Z),(D1×Z),(D1×Z)对应的回归系数

R示例:

检验结果的解读:

  • 发现p=0.3080.表明(失业人数)对(就业人数与国民生产总值之间关系)的影响是不显著的。
  • R2从0.984增加到0.985,增加不显著。
  • ANOVA的结果P=0.31也同样表明两个模型在统计上没有显著的改变。

预测变量的中心化(centralization)

中心化处理是对变量值进行转变,转变为以0为均值

公式为:XC=X-MX

如下图所示,中心化处理后改变的只是回归常数和中心点的位置:

QQ截图20140102221909

QQ截图20140102221916

这么做原因:概念上:

  • 如果所有的预测变量都不可能为0,则截距/回归常数B0的解读是无意义的
  • 如果不存在调节效应,无论Z的取值,B1的取值是稳定的
  • 如果存在调节效应,Z的取值不同,B1的取值是变化的

经过中心化处理后,回归系数,R,R2都不会,但可以使得对分析结果的解读更加容易。

统计上:

  • 避免多元共线性(multicolinearity),如果一般线性模型中的两个变量之间相关性很高,则两者是冗余的,预测两者分别对应的回归系数会很困难

中介(mediation)

中介分析(mediation analysis)被用来更好地理解观测到的自变量对因变量的影响机制,即:是否有一个中间变量M,X通过M进而影响Y。

如果X与Y是相关的,但是有中介变量M在其中起作用(X->M->Y),意味着:

  • Y=B0+B1M+e
  • M=B0+B1X+e

因此合起来:

  • Y=B0+B1M+B2X+e

我们可以看回归系数B2是仍然否显著,如果仍然显著,则M可能并不起作用,但如果B2不再显著,则M可能为中介变量。

根据中介变量M对X与Y之间的关系起到作用的程度,则可分为:

  • 起到部分作用:部分中介
  • 起到全部作用:完全中介

要检验中介效应构建三个模型:

  • lm(Y~X)
  • lm(M~X)
  • lm(Y~X+M)

其中:

  • lm(Y~X)中X的回归系数应该显著
  • lm(M~X)中X的回归系数应该显著
  • lm(Y~X+M)中M的回归系数显著,看X的回归系数如何变化

R示例

结果发现,显著性由8.4e-12 降至5.0e-05,回归系数由27降至11,表明人口应该是起到了部分中介的作用。

路径模型(path model)

中介分析通常使用路径模型来分析,路径模型是结构化方程的一种形式,是一种图形模型。

路径模型中:

  • 矩形代表观测到的变量(X,Y,M)
  • 圆圈代表未观测到的变量(误差e)
  • 三角代表常量
  • 箭头代表关系

没有中介效应的模型:

QQ截图20140102170054

  • Y受到X的影响,影响的程度为B1。
  • Y受到常数1和误差e的影响,影响程度均为1。

有中介效应的模型:

QQ截图20140102225306

  • Y受到X的影响,影响的程度为c'。
  • Y受到常数1和误差e的影响,影响程度均为1。
  • Y受到M的影响,影响程度为b
  • M受到X的影响,影响程度为a
  • M受到误差e的影响,程度为1

要检验的是c'和B1之间的差距。

另外有a*b是X对Y的间接影响的程度。

可以用Sobel检验做中介效应的检验,其检验的是这个X->M->Y的间接影响是否是显著的。

零假说为:间接影响为0

检验中z值的计算公式:z=(Ba*Bb)/SQRT[(Ba2*SEb2)+(Bb2*SEa2)]

注意这个z值本质仍然是一个比例,是观测到的情况/随机的情况。

R示例:

sobel检验

上例中这个z 值的大小表明,影响是显著的。

Leave a Reply

Your email address will not be published. Required fields are marked *