利用主成分分析(PCA)来优化线性混合模型(brief tutorial)

利用主成分分析(PCA)来优化线性混合模型(brief tutorial)

线性混合模型(Linear Mixed Model,LMM)中非常重要的步骤是确定模型应包含哪些随机斜率,特别是在因子设计的实验中。Barr 等人(2013) 认为,应从最大模型(maximum model)开始拟合,如果不能拟合或者出现畸形协方差矩阵再将随机斜率逐个减去:

Recently, Barr, Levy, Scheepers, and Tily (2013) recommended fitting ‘maximal’ models with all possible random effect components included.

但是有两方面的问题:一是模型本身是否能支撑起如此多随机斜率,是不确定的;二是随着实验设计的复杂程度提升,随机斜率上可能的组合数量会急剧提升,因此筛选的工作量也会急剧增多(图 1):

图 1. 因素数量及其内部水平数量与不同随机斜率组合数量的关系

这也是Barr等人的建议提出后,大量研究者报告其建立的全模型难以收敛的原因。之后,Bates等人(2015)提出,要解决这个问题,应该要了解混合模型方法本身的若干限制。

在LMM中随机效应的估计比固定效应的估计复杂得多。比如在一个2*2*2的实验设计中,假如三个因素都是被试内且是项目内设计的,那么模型中一共有8个固定效应和73个随机效应。从现实角度上讲,为了获取主要的数据结果,没有必要拟合这么多高度抽象的参数,尤其是考虑到被试数和项目数时。但是LMM本身是一定要根据你的设计或者你给的参数进行拟合的,这样就导致畸形的协方差矩阵,具体表现就是某个随机效应的方差格外小,或者某两个随机效应之间的相关性格外大。

从原理上讲,出现这样的问题,是因为在估计随机效应时,出现了系数矩阵秩亏,主要是由于实际观测数小于必要观测数引起的,并导致矩阵降维。这样矩阵的维度数量和随机效应的数量出现偏差,或者说存在某个或某些随机效应,其不能代表一个单独的维度或者解释某一部分独特的方差。

好比是你做了一桌子饭菜,结果要宴请两桌子客人,那么有些客人是没饭吃的,没饭吃的这些人无所事事,就会在你家乱搞,此时你要把这些无所事事的人轰出去,因此问题就演变成如何确定这些人呢?

作者提出可利用主成分分析的方法确定有用的随机效应成分的个数,根据模型中具体随机效应的标准差的大小(能解释的变异量)来确定要删除哪个成分。其具体的原则/步骤是:

  1. 前提假设 1:研究假设只与固定效应有关,与随机效应无关;
  2. 前提假设 2:不考虑潜在的协方差
  3. 即使全模型可以收敛,也要检查是否可以削减方差-协方差的维度,维度的个数通过主成分分析中累积解释变异刚达到100%为标准为定;
  4. 检查采用零相关参数的模型是否显著减少模型的拟合程度(可采用AIC或BIC等指标);
  5. 确定主成分个数后,着手减少成分;对削减后的模型再次检查主成分,确保随机效应的结构;
  6. 确定简化后的模型中,每个随机效应成分是否是必需的,方法为每次只去掉一个成分,看是否对模型的拟合程度产生显著影响;
  7. 将零相关的模型扩展到相关模型后,是否会影响拟合程度(稳定的随机效应成分是考察其相关系数的前提);
  8. 一个特例,当实验设计非常复杂时,全模型可能不能收敛,此时首要的步骤是检查其对应的零相关参数模型的主成分维度。因为从全模型到零相关模型本身就是最大的简化。

其中最后一条虽然说是特例,但是其实是非常常见的情形。

这里借助 Schad 等人(2020)编写的函数MixedDesign()(available from OSF)来模拟数据,为大家演示一下:

DF = c()
for(ss in 1:30){
  DF = rbind(DF,mixedDesign(W = 6,B = NULL, n = 30,
                            M = matrix(c(310,300,290,280,280,280),nrow = 1),
                            SD = 20,long = T,)
  )
}
DF = DF %>% as_tibble()
head(DF,3)
# A tibble: 3 x 3
#   id    W_a      DV
#   <fct> <fct> <dbl>
# 1 1     a1     292.
# 2 1     a2     303.
# 3 1     a3     267.

这里模拟了一个单因素6水平设计的实验数据,实验中共30名被试,每名被试在每个条件下共有30个观测值,即共有30 * 30 * 6 = 5400个观测值。六个水平下的均值依次为310, 300, 290, 280, 280, 280。标准差为20。

下面对数据进行一些整理,将变量id的名称改为subj(被试),把变量W_a的名称改为cond(条件):

> names(DF)[[1:2]] = c('subj','cond')
> head(DF,3)
# A tibble: 3 x 3
#   subj  cond    DV
#   <fct> <fct> <dbl>
# 1 1     a1     292.
# 2 1     a2     303.
# 3 1     a3     267.

下面以第一个水平a1为基线,考察各水平与基线之间的差异:

> levels(DF$cond) = paste0('a',1:6)
> levels(DF$cond)
# [1] "a1" "a2" "a3" "a4" "a5" "a6"

建立全模型:

> ModelAll = lmer(data = DF, DV~cond+(1+cond|subj), control = lmerControl(optimizer = 'bobyqa'))
# boundary (singular) fit: see ?isSingular

结果显示出现了畸形协方差矩阵。接下来用rePCA()函数对其随机效应部分进行主成分分析:

> summary(rePCA(ModelAll))
# $subj
# Importance of components:
#                          [,1]   [,2]   [,3]    [,4]      [,5] [,6]
# Standard deviation     0.2931 0.1850 0.1354 0.02149 7.642e-07    0
# Proportion of Variance 0.6184 0.2462 0.1320 0.00332 0.000e+00    0
# Cumulative Proportion  0.6184 0.8647 0.9967 1.00000 1.000e+00    1

看到,整个随机效应共有6个成分,具体是哪六个可以看一下模型随机效应的结果:

> VarCorr(ModelAll)
# Groups   Name        Std.Dev. Corr                              
# subj     (Intercept)  2.5795                                    
#          conda2       3.2674  -0.665                            
#          conda3       2.2161  -0.657 -0.110                     
#          conda4       3.3645  -0.833  0.554  0.647              
#          conda5       3.4440  -0.677  0.123  0.672  0.282       
#          conda6       2.7756  -0.837  0.938  0.147  0.611  0.458
# Residual             19.5547 

6个随机效应为subj上的随机截距,和5个随机斜率。其中随机斜率即为cond产生的5个虚拟变量。

主成分分析的结果(第二行为每个因素能够解释的方差比例)表明,只存在4个成分,其解释的方差比例大于0。换句话说,有2个成分是多余的,因此需要对模型进行优化。

将全模型简化为零相关系数(zero-correlation-parameter, ZCP)的模型:

> ModelZCP = lmer(data = DF, DV~cond+(1+cond||subj), control = lmerControl(optimizer = 'bobyqa')) # 用 '||' 代替 ‘|’
# boundary (singular) fit: see ?isSingular
# Warning message:

# Model failed to converge with 2 negative eigenvalues: -2.0e-03 -2.1e+01 

在零相关系数的模型下,模型依然无法收敛,依然存在畸形协方差矩阵,需要继续优化,删减随机斜率。

这时遇到一个难题:在建立以上两个模型时,我们直接把原始的因子变量放入模型,在删减随机斜率时,只能一次性删掉所有的斜率,但并不是所有的斜率都导致模型出现畸形协方差,我们要删除的只是部分随机斜率。

为了解决这个问题,我们需要在建模前首先生成该因子的虚拟变量,并用虚拟变量建模,而不是用原始因子建模:

相关内容请见:

张光耀:线性(混合)模型中如何对因子变量事先生成虚拟变量zhuanlan.zhihu.com图标
> Dummy = model.matrix(~cond,data = DF)
> head(Dummy)
#   (Intercept) conda2 conda3 conda4 conda5 conda6
# 1           1      0      0      0      0      0
# 2           1      1      0      0      0      0
# 3           1      0      1      0      0      0
# 4           1      0      0      1      0      0
# 5           1      0      0      0      1      0
# 6           1      0      0      0      0      1
> DF[paste0('a',2:6)] = Dummy[,2:6]
> head(DF)
#   A tibble: 6 x 8
#   subj  cond     DV    a2    a3    a4    a5    a6
#   <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 1     a1     292.     0     0     0     0     0
# 2 1     a2     303.     1     0     0     0     0
# 3 1     a3     267.     0     1     0     0     0
# 4 1     a4     292.     0     0     1     0     0
# 5 1     a5     305.     0     0     0     1     0
# 6 1     a6     247.     0     0     0     0     1

接下来用虚拟变量建模:

> ModelAll = lmer(data = DF, DV ~ a2 + a3 + a4 + a5 + a6 + (1 + a2 + a3 + a4 + a5 + a6|subj), control = lmerControl(optimizer = 'bobyqa'))
# boundary (singular) fit: see ?isSingular

全模型随机效应主成分分析的结果:

> summary(rePCA(ModelAll))
# $subj
# Importance of components:
#                          [,1]   [,2]   [,3]    [,4]      [,5] [,6]
# Standard deviation     0.2931 0.1850 0.1354 0.02149 7.642e-07    0
# Proportion of Variance 0.6184 0.2462 0.1320 0.00332 0.000e+00    0
# Cumulative Proportion  0.6184 0.8647 0.9967 1.00000 1.000e+00    1

可以看到,用原始因子建模和用虚拟变量建模,随机效应主成分分析结果也是一致的。

建立零相关系数模型:

> ModelZCP = lmer(data = DF, DV ~ a2 + a3 + a4 + a5 + a6 + (1 + a2 + a3 + a4 + a5 + a6||subj), control = lmerControl(optimizer = 'bobyqa'))
# boundary (singular) fit: see ?isSingular

模型依然出现畸形协方差,需要削减随机斜率。

由主成分分析结果可知,如果采用每个成分解释方差的比例需要大于0的标准,需要削减2个随机斜率。零相关模型中,随机效应为:

> VarCorr(ModelZCP)
 Groups   Name        Std.Dev.
 subj     (Intercept)  0.57799
 subj.1   a2           1.99048
 subj.2   a3           0.00000
 subj.3   a4           1.72675
 subj.4   a5           2.46727
 subj.5   a6           0.00000
 Residual             19.61359

可以看到,存在两个随机效应a3a6,其标准差为0,需要删去:

> ModelZCP2 = lmer(data = DF, DV ~ a2 + a3 + a4 + a5 + a6 + (1 + a2 + a4 + a5||subj), control = lmerControl(optimizer = 'bobyqa'))

删减后,模型可以收敛,且不存在畸形协方差。下面检验优化后的模型和全模型(注意,不是零相关模型)之间的差异:

> anova(ModelZCP2, ModelAll)
# refitting model(s) with ML (instead of REML)
# Data: DF
# Models:
# ModelZCP2: DV ~ a2 + a3 + a4 + a5 + a6 + (1 + a2 + a4 + a5 || subj)
# ModelAll: DV ~ a2 + a3 + a4 + a5 + a6 + (1 + a2 + a3 + a4 + a5 + a6 | subj)
#           Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)
# ModelZCP2 11 47514 47586 -23746    47492                        
# ModelAll  28 47536 47721 -23740    47480 11.46     17     0.8317

差异不显著。此时如果再把优化后的模型扩充为相关系数模型,模型依然会出问题:

> ModelZCP3 = lmer(data = DF, DV ~ a2 + a3 + a4 + a5 + a6 + (1 + a2 + a4 + a5|subj), control = lmerControl(optimizer = 'bobyqa'))
boundary (singular) fit: see ?isSingular
Warning message:
Model failed to converge with 1 negative eigenvalue: -9.2e+00 

因此采用ModelZCP2即可。

关于此方法, 有四点需要说明,第一,关于对零相关模型优化后,是否需要再考察含有相关系数的模型的问题,我还没有确切的想法和结论。按照Bates等人 (2015) 的看法,对零相关模型优化后,需要再考察相关系数模型,但实际操作中,由优化后的零相关模型扩展到含有相关系数模型时,后者往往还会出问题,我个人在数据分析时一般不再考察含有相关系数的模型,但承认考察相关系数模型的重要性。

第二,关于删除主成分的标准的问题,我个人倾向采用“先松后严”的方法,先以较宽松的标准(比如对方差的解释比例 > 0)削减,若削减后仍有问题,再采用较严格的标准(比如对方差的解释比例 > 1%)削减。

第三,关于删除随机斜率时应该逐个删除,还是一次性删除。我个人倾向于采用逐个删除的方法,但如果存在若干标准差过于小的随机斜率(比如上述例子中,两个斜率的标准差都是0),我会一次性地删除。

第四,这里的例子中我是按照全模型的主成分分析结果来确定要删几个成分,但在实际删除时,是从零相关模型中删除的。考虑到全模型和零相关模型的主成分分析结果存在差别,我个人在实际操作中会综合考虑这两个模型的主成分分析结果。

  • 参考文献

Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal.Journal of memory and language,68(3), 255 – 278.

Bates, D., Kliegl, R., Vasishth, S., & Baayen, H. (2015). Parsimonious mixed models.arXiv preprint arXiv:1506.04967.

Schad, D. J., Vasishth, S., Hohenstein, S., & Kliegl, R. (2020). How to capitalize on a priori contrasts in linear (mixed) models: A tutorial.Journal of Memory and Language,110, 104038.

编辑于 01-20

文章被以下专栏收录