||
多层线性模型(Hierarchical Linear Model,HLM),也叫多水平模型(Multilevel Model,MLM),是社会科学常用的高级统计方法之一,它在不同领域也有一些近义词或衍生模型:
线性混合模型(Linear Mixed Model)
混合效应模型(Mixed Effects Model)
随机效应模型(Random Effects Model)
随机系数模型(Random Coefficients Model)
方差成分模型(Variance Components Model)
增长曲线模型(Growth Curve Model)
……
本文分为两个部分:
1. 什么是HLM(入门科普)
2. HLM的自由度计算问题(进阶探讨)
子曾经曰过:“HLM是众多统计方法的「幕后主谋」。”什么意思呢?一般我们收集完数据之后会做什么分析?做实验的人会说:方差分析(ANOVA);做调查的人会说:多元回归(multivariate regression);做文献综述的人会说:元分析(meta-analysis)……但其实,这一切都可以视为HLM的特例:
t检验是ANOVA的一个特例(自变量只有两水平)
ANOVA是回归分析的一个特例(自变量为分类变量)
回归分析的实质是一般线性模型GLM(其推广则是广义线性模型)
GLM是HLM的一个特例(只有Level 1)
元分析可以视为只有组间模型的HLM(只有Level 2)
简而言之:回归分析是几乎所有统计模型的基础,而回归分析的最一般形式则可以归为多层线性模型HLM。
当然,除去上面这些HLM的特例,HLM本身关注的焦点其实是“多层嵌套数据”。
举个栗子,假设我们想知道“智力水平(IQ)能否影响学业成绩(GPA)”,然后我们收集了很多很多学校的数据,迫不及待地画了一个散点图:
一般线性模型(GLM),即最普通的回归分析
这样一看,IQ和GPA并无什么关系呀,甚至还有一点负相关?!于是我们就下结论“IQ与GPA无关”吗?emmm等等,我们漏了什么?
事实上,不同学校之间在很多方面可能都存在固有的差异,比如学生整体的智力情况、教学条件、师资力量等等。现在考虑一个情形:收集的学校里面既有“精英学校”也有“普通学校”,学生入学的标准就是他们的IQ。
然后我们重新对不同学校分别画一下散点图:
多层线性模型(HLM)
你会发现,如果考虑了学校间的固有差异,或者说“基线水平”差异,IQ和GPA之间的关系就不会受到这些因素的混淆了(咳咳,这里的数据都是我瞎编的啊,别当真)。以上就是经典的“组内同质、组间异质”的情况,我们刚刚做的其实就可以视为HLM的其中一种子模型——随机截距-固定斜率模型,也就是假定不同学校的基线水平不同(随机截距),但IQ与GPA之间的变量关系在不同学校中保持相同(固定斜率)。
HLM会把多层嵌套结构数据在因变量上的总方差进行分解:
总方差 = 组内方差(Level 1)+ 组间方差(Level 2)
——咦,你是不是联想到了方差分析?没错!比如在上面的例子中,学生是个体水平(Level 1)的分析单元,IQ和GPA都是在个体水平收集的变量,而学校是群体水平(Level 2)的分析单元,不过我们暂时并没有收集学校水平的任何自变量,只是把学校本身当做一个分组变量(clustering/grouping variable)。换句话说,上面这个例子也可用被称作“随机效应单因素协方差分析(ANCOVA with random effects)”。
现在用公式来表示上面这个例子:
我们既可以用R语言lme4、lmerTest包的lmer函数,也可以用SPSS的MIXED语句进行建模:
# R语言 - 固定斜率:model = lmer(GPA ~ IQ + (1|School), data=data)# SPSS - 固定斜率:MIXED GPA WITH IQ /METHOD=REML /PRINT=DESCRIPTIVES SOLUTION TESTCOV /FIXED=IQ | SSTYPE(3) /RANDOM=INTERCEPT | SUBJECT(School) COVTYPE(VC).
或者可以设置IQ为随机斜率(即理论上假定IQ与GPA之间的关系在不同学校是不一样的):
# R语言 - 随机斜率:model = lmer(GPA ~ IQ + (IQ|School), data=data)# SPSS - 随机斜率:MIXED GPA WITH IQ /METHOD=REML /PRINT=DESCRIPTIVES SOLUTION TESTCOV /FIXED=IQ | SSTYPE(3) /RANDOM=INTERCEPT IQ | SUBJECT(School) COVTYPE(UN).
当然,我们还可以引入学校水平的自变量来对学校间的GPA均值差异进行解释,比如教师数量、教学经费……这些变量由于只在学校层面变化,对于每个学校内的每一个学生而言都只有一种可能的取值,因此必须放在Level 2的方程中作为群体水平自变量,而不能简单地处理为个体水平自变量——这也就是HLM的另一个存在的意义:可以同时纳入分析个体与群体水平的自变量。
我们来看HLM的一般式:
Level 1(组内/个体水平,或重复测量/追踪设计中的时间水平;p个自变量,i个样本量):
Level 2(组间/群体水平,或重复测量/追踪设计中的个体水平;q个自变量,j个样本量):
以上就是关于HLM的入门级科普。我不确定是不是所有读者都对HLM有足够的了解,因为很多同学其实是没有接触过HLM的,我自己也是2017年暑假才开始学习HLM。不仅如此,一些实验取向、认知取向的同学甚至会怀疑回归分析的意义,因为在他们的研究中,t检验和ANOVA就足够用了——而这些分析的本质其实都是回归分析,回归分析的更一般形式则是HLM。
另外,也有不少同学把“多层线性模型”与“分层/逐步多元回归”搞混淆——请注意:
多层线性模型HLM解决的是多层嵌套结构数据(落脚点是数据结构)
分层/逐步多元回归本身是普通的回归分析,解决的是不同自变量的重要性的优先程度(落脚点是变量重要性)
前面这些都只是铺垫,因为下面的内容才是重点。我们来谈谈HLM的自由度计算问题。
自由度(degree of freedom, df),简单来说就是“能够独立变化的数据量”。比如一组有N个数据的样本,其总体平均值是确定的,那么在参数估计中,我们说它的自由度是N – 1,这里的1代表已经确定了的均值,也就是说无论你的数据怎么变,只要给我N – 1个数,我就能知道剩下的那个数是几,因为均值已经确定了。
而在普通的回归分析中,同样道理,截距相当于一个已经确定了的均值,是我们要估计的一个参数,并且每增加一个自变量,都会相应增加一个回归系数,这些回归系数也是我们要估计的参数,也视为确定值,因此剩下的那些不确定的、能够独立变化的数据的个数就是自由度,主要用于假设检验(对回归系数的检验:t = b/SE,服从df = N – k的t分布,其中N为总样本量,k为所有的参数个数,截距也占一个参数)。
一言以蔽之:
df (能够独立变化的数据量)= N(总数据量)– k(待估计的参数个数 [包括截距])
事实上任何一个回归方程都有且仅有一个截距,就好比任何一组数据都有且仅有一个平均值,所以很多时候上面这个表达式也可以用“自变量/预测变量”的个数来表示:
df = N – k = N – p – 1(p表示自变量/预测变量的个数,1表示截距这个参数)
然而,这么一目了然的事情在HLM里面却是一件比较复杂的议题,因为前面我们已经看到,HLM的回归方程不止一个,不仅有Level 1的一个方程,还有Level 2的一系列以Level 1的每个参数(包括截距和斜率)为因变量建立的方程,因此不同的回归系数因其对应的变量有不同的性质(或者说在HLM中所处的不同层级),其自由度也不尽相同。
你可能要问了:
HLM的自由度真的是一个问题吗?
HLM发展了这么多年,难道连一个小小的自由度问题都还没解决?
我相信不仅是我自己有这样的疑问,各位读者也会产生类似的疑问。我遇到HLM自由度的问题也是在前不久,在使用R语言的lme4包和lmerTest包做HLM的时候,lmerTest输出的某个Level 2自变量的自由度(20多)远远小于它本身的Level 2样本量(400多),非常离谱,并且我做的所有模型的自由度都是小数而不是整数。无独有偶,室友
也遇到了HLM自由度的问题(混合线性模型(LMM)中,固定效应(Fixed Effects)中的自由度是如何计算的?)。
碰巧,这几天我查了很多参考书和文献,终于把HLM的自由度问题给搞清楚了!
是的,不要小看它!我们先来看一篇2009年发表在《Trends in Ecology and Evolution》上面的综述文章「Generalized linear mixed models: A practical guide for ecology and evolution」。
文章提到:不同软件在计算HLM自由度的时候,差异非常之大(vary enormously),HLM自由度的计算目前仍然是一个challenge和difficulty。
Bolker et al., 2009
大体分两类:“简单计算法”和“近似估计法”。
在综合了很多中英文教材之后,我把它们归纳如下:
HLM自由度算法总结
HLM自由度算法总结(一览表)
常见自由度类型举例
总结:
HLM软件(软件本身叫做HLM)作为最经典的HLM统计工具,输出的都是由“简单计算法”估计的自由度,并且都是整数,不存在小数。这种方法理解起来也很简单,就是我们上面提到的原则:df (能够独立变化的数据量)= N(总数据量)– k(待估计的参数个数 [包括截距])。只不过N和k都需要具体看是哪一层的模型,对于Level 1的模型,与普通回归基本无异,而对于Level 2的模型,N不再是个体水平的总样本量,而是群体水平的“组数量”,k同样也变成了群体水平的参数个数(注意:当level-1的某变量被设定为随机斜率,df计算略复杂,可参考上面的Table 1)。由于HLM优于其他软件的一点在于它可以准确规定哪个变量属于哪个水平,所以虽然是简单计算法,但其实也是一种理论驱动的方法。
SPSS和R语言的lmerTest包均默认输出Satterthwaite(1946)的自由度近似估计值(approximation),这种算法是基于每一个变量在不同level的残差方差(residual variances)做的校正估计,严格来说是一种数据驱动的方法,主要用于校正方差不齐/异方差、每组内样本量不均衡、组数量较少等情况。这种近似估计往往得不到整数,而是会出现小数,但正常情况下,估计出来的自由度与简单计算法也是处于同一个数量级的,不会偏差太多。而当每组内样本量相等时,Satterthwaite的近似估计自由度与简单计算法得到的自由度就是一致的了,并且都是整数。
R语言的lmerTest包,除了Satterthwaite近似估计,还提供了Kenward-Roger近似估计,但后者不如前者好用(比如百万数量级的数据如果用KR估计就会直接崩溃)。
R语言的lme4包(最根源的HLM程序包,lmerTest只是调用了它并增添了假设检验结果)早期也会报告自由度和假设检验结果,但后来取消了这一功能,我推测可能是因为HLM自由度的估计尚存争议,所以希望用户自己来选择使用哪种自由度。
在明白了上面的理论问题之后,我就自编了一个R函数,从而可以实现像HLM软件那样自己根据理论假设去设定每个变量/参数的类型(如intercept、L1 fixed、L1 random、L2、cross-level interaction),分别计算不同类型的df,并输出假设检验结果和效应量。在实际分析了自己的数据之后,我发现虽然“简单计算法”和“近似估计法”对df的估计有时确实相差挺多,但显著性检验的结果(p值)基本是一致的,这可能是因为我的数据里Level 2的group数量已经很多(几百或几千),即使df有校正,影响也微乎其微。
# 前半部分是lmerTest的输出(t-tests use Satterthwaite's method)# 最后四列是我自编函数的输出(简单计算法)# p值其实是基本一致的,但df差别很大## Fixed effects:# Gamma S.E. t df.apprx p.apprx df.HLM p.HLM r.prtl R2.prtl sig# (Intercept) -0.783 0.008 -99.51 969.5 < 2e-16 NA NA NA NA NA# NU 0.008 0.001 6.90 1016.1 9e-12 1697 7.2e-12 0.165 0.027 ****# CCU 0.003 0.001 3.10 717.9 0.00204 1697 0.00200 0.075 0.006 **# NV 0.003 0.001 2.21 849.1 0.02718 1697 0.02705 0.054 0.003 *# NG 0.002 0.001 1.31 655.3 0.19107 1697 0.19079 0.032 0.001 # SNU 0.001 0.001 0.67 124.1 0.50285 476 0.50193 0.031 0.001 # SNI -0.000 0.000 -0.35 106.8 0.72421 476 0.72367 -0.016 0.000 # sex -0.095 0.001 -87.79 1544855.3 < 2e-16 2509762 < 2e-16 -0.055 0.003 ****# age 0.003 0.000 46.86 1120937.4 < 2e-16 2509762 < 2e-16 0.030 0.001 ****# edu 0.108 0.001 213.28 2492317.9 < 2e-16 2509762 < 2e-16 0.133 0.018 ****# ---# Signif. codes: 0 ‘****’ .0001 ‘***’ .001 ‘**’ .01 ‘*’ .05 ‘.’ .10 ‘ ’ 1# Note #1: df.apprx is estimated by Satterthwaite's (1946) approximation.# Note #2: r.partial is calculated by t-to-r transformation.
然而,如果使用“t-to-r”转换来计算multilevel效应量,那么df的差异就很要命了,因为r = sqrt(t2/(t2+df)),df会直接影响最后算出来的效应量(Satterthwaite常常会低估df,使得效应量被高估)。
所以!自由度的估计方法,对显著性检验(p值)的影响其实并不大,但是对效应量计算(如t-to-r转换)的影响很大。
我个人的选择是:用R作为统计工具(因为SPSS跑大数据模型实在太吃力了),依然使用lme4和lmerTest包,但在自由度方面还是遵循HLM软件所使用的简单计算法,同时也看一眼lmerTest输出的Satterthwaite校正结果,双重保障。
作者:包含吴霜
链接:https://zhuanlan.zhihu.com/p/50048784
来源:知乎 著作权归作者所有。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 09:25
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社