奇异随机效应协方差矩阵

获得+1或-1的随机效应相关性估计意味着优化算法达到了“边界”：相关性不能高于+1或低于-1。即使没有明确的收敛错误或警告，这也可能表明收敛存在一些问题，因为我们不期望真正的相关性位于边界上。正如您所说，这通常意味着没有足够的数据来可靠地估计所有参数。马图舍克等人。2017说，在这种情况下，权力可能会受到影响。

达到边界的另一种方法是获得 0 的方差估计：为什么我的混合模型中随机效应的方差为零，尽管数据有一些变化？

这两种情况都可以看作是获得随机效应的退化协方差矩阵（在您的示例中，输出协方差矩阵是$4\times 4$); 零方差或完全相关意味着协方差矩阵不是满秩的，并且[至少]其中一个特征值为零。这一观察立即表明，还有其他更复杂的方法来获得退化协方差矩阵：可以有一个$4\times 4$协方差矩阵没有任何零点或完全相关，但仍然秩不足（奇异）。贝茨等人。2015 Parsimonious Mixed Models（未发表的预印本）推荐使用主成分分析（PCA）来检查获得的协方差矩阵是否是奇异的。如果是，他们建议以与上述单一情况相同的方式处理这种情况。

那么该怎么办？

如果没有足够的数据来可靠地估计模型的所有参数，那么我们应该考虑简化模型。以您的示例模型为例X*Cond + (X*Cond|subj)，有多种可能的方法可以简化它：

1. 去除其中一种随机效应，通常是最高阶相关：

   X*Cond + (X+Cond|subj)
   

2. 去掉所有相关参数：

   X*Cond + (X*Cond||subj)
   

   更新：正如@Henrik 所指出的，||只有在其左侧的所有变量都是数字的情况下，该语法才会删除相关性。如果涉及分类变量（例如Cond），则应该使用他方便的afex包（或繁琐的手动解决方法）。有关更多详细信息，请参阅他的答案。

3. 通过将术语分成几个来摆脱一些相关参数，例如：

   X*Cond + (X+Cond|subj) + (0+X:Cond|subj)
   

4. 以某种特定方式约束协方差矩阵，例如，按照您的建议，将一个特定的相关性（触及边界的相关性）设置为零。没有内置的方法lme4可以实现这一点。请参阅@BenBolker 对 SO 的回答，了解如何通过一些智能黑客来实现这一目标。

与你所说的相反，我不认为Matuschek 等人。2017特别推荐#4。Matuschek 等人的要点。2017 年和贝茨等人。2015 年似乎是从最大模型 a la Barr 等人开始的。2013然后降低复杂度，直到协方差矩阵满秩。（此外，他们通常会建议进一步降低复杂性，以增加功率。）更新：相比之下，Barr 等人。建议仅在模型未收敛时降低复杂度；他们愿意容忍奇异的协方差矩阵。请参阅@Henrik 的回答。

如果有人同意 Bates/Matuschek，那么我认为可以尝试不同的方法来降低复杂性，以便找到能够在“损害最小”的同时完成工作的方法。看我上面的列表，原来的协方差矩阵有 10 个参数；#1 有 6 个参数，#2 有 4 个参数，#3 有 7 个参数。如果不拟合它们，就不可能说哪种模型会摆脱完美的相关性。

但是如果你对这个参数感兴趣怎么办？

上述讨论将随机效应协方差矩阵视为令人讨厌的参数。如果您对必须“放弃”以获得有意义的全等级解决方案的相关参数特别感兴趣，您会提出一个有趣的问题。

请注意，将相关参数固定为零不一定会产生ranef不相关的 BLUP ( )；实际上，它们甚至可能根本不会受到太大影响（请参阅@Placidia 的演示答案）。因此，一种选择是查看 BLUP 的相关性并进行报告。

另一个可能不太有吸引力的选择是使用treat subjectas a fixed effect Y~X*cond*subj，获取每个主题的估计值并计算它们之间的相关性。这相当于Y~X*cond对每个主题分别运行单独的回归并从中获得相关估计。

另请参阅Ben Bolker 的混合模型常见问题解答中关于奇异模型的部分：

过度拟合的混合模型导致奇异拟合是很常见的。从技术上讲，奇点意味着一些$\theta$（方差-协方差Cholesky分解）Cholesky因子的对角元素对应的参数正好为零，是可行空间的边缘，或者等价于方差-协方差矩阵有一些零特征值（即是半正定而不是正定)，或（几乎等价地）一些方差估计为零或一些相关性估计为 +/-1。

我同意阿米巴的回答中所说的一切，它很好地总结了当前关于这个问题的讨论。我将尝试添加一些额外的点，否则请参阅我最近的混合模型课程的讲义，其中也总结了这些点。

抑制相关参数（变形虫答案中的选项 2 和 3）||仅适用于数值协变量lmer，不适用于因子。Reinhold Kliegl 的代码对此进行了详细讨论。

但是，如果调用中的 参数（另请参见函数） ，我的afex包提供了抑制因素之间相关性的功能。它本质上是通过实施 Reinhold Kliegl 讨论的方法来实现的（即，将因素转换为数值协变量并指定这些因素的随机效应结构）。expand_re = TRUEmixed()lmer_alt()

一个简单的例子：

library("afex")
data("Machines", package = "MEMSS") # same data as in Kliegl code

# with correlation:
summary(lmer(score ~ Machine + (Machine  | Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr       
#  Worker   (Intercept) 16.6405  4.0793              
#           MachineB    34.5467  5.8776    0.48      
#           MachineC    13.6150  3.6899   -0.37  0.30
#  Residual              0.9246  0.9616              
# Number of obs: 54, groups:  Worker, 6

## crazy results:
summary(lmer(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr     
#  Worker   (Intercept)  0.2576  0.5076            
#  Worker.1 MachineA    16.3829  4.0476            
#           MachineB    74.1381  8.6103   0.80     
#           MachineC    19.0099  4.3600   0.62 0.77
#  Residual              0.9246  0.9616            
# Number of obs: 54, groups:  Worker, 6

## as expected:
summary(lmer_alt(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  16.600   4.0743  
#  Worker.1 re1.MachineB 34.684   5.8894  
#  Worker.2 re1.MachineC 13.301   3.6471  
#  Residual               0.926   0.9623  
# Number of obs: 54, groups:  Worker, 6

对于那些不知道的人afex，混合模型的主要功能是为固定效应提供 p 值，例如：

(m1 <- mixed(score ~ Machine + (Machine  || Worker), data=Machines, expand_re = TRUE))
# Mixed Model Anova Table (Type 3 tests, KR-method)
# 
# Model: score ~ Machine + (Machine || Worker)
# Data: Machines
#    Effect      df        F p.value
# 1 Machine 2, 5.98 20.96 **    .002
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

summary(m1)  
# [...]
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  27.4947  5.2435  
#  Worker.1 re1.Machine1  6.6794  2.5845  
#  Worker.2 re1.Machine2 13.8015  3.7150  
#  Residual               0.9265  0.9626  
# Number of obs: 54, groups:  Worker, 6
# [...]

来自 Barr 等人的 Dale Barr。(2013) 论文在建议减少随机效应结构方面比在变形虫的答案中提出的更为谨慎。在最近的一次推特交流中，他写道：

• “减少模型会引入未知的反保守风险，如果有的话，应该谨慎行事。” 和
• “我主要担心的是人们了解与模型简化相关的风险，并且最小化这种风险需要一种比通常采用的更保守的方法（例如每个斜率测试为 0.05）。”

所以建议谨慎。

作为审稿人之一，我还可以提供一些关于为什么我们Bates 等人的见解。(2015)论文仍未发表。我和另外两位审稿人（已签名，但在此未透露姓名）对 PCA 方法提出了一些批评（似乎没有原则，并且没有证据表明它在权力方面更胜一筹）。此外，我认为这三个人都批评该论文没有关注如何指定随机效应结构的问题，而是试图引入 GAMM。因此，Bates 等人 (2015) 的论文演变为Matuschek 等人的论文。（2017 年）论文通过模拟和Baayen 等人解决了随机效应结构的问题。(2017)介绍 GAMM 的论文。

我对贝茨等人的完整评论。可以在这里找到草稿。IIRC，其他评论有类似的要点。

在使用最大似然估计时，我也遇到过这个问题——只有我使用通过 MLwiN 软件实现的 Goldstein IGLS 算法，而不是 R 中的 LME4。但是，在每种情况下，当我使用相同的方法切换到 MCMC 估计时，问题都已解决软件。我什至有超过 3 的相关性，这在我改变估计时解决了。使用 IGLS，相关性在估计后计算为协方差除以相关方差乘积的平方根的乘积 - 这没有考虑每个组成估计中的不确定性。

IGLS 软件不“知道”协方差意味着相关性，而只是计算常数、线性、二次等方差函数的估计值。相比之下，MCMC 方法建立在多元正态分布的样本假设之上，该样本对应于具有良好特性和完全误差传播的方差和协方差，因此在估计方差时考虑了协方差估计中的不确定性反之亦然。

MLwin 是具有 IGLS 估计和非负定方差协方差矩阵的 MCMC 估计链，可能需要通过在开始采样之前将协方差更改为零来进行更改。

有关工作示例，请参见

使用 MLwiN 3 第 1 卷（2017 年 9 月更新）开发用于分析上下文、异质性和变化的多层次模型；第 2 卷也在 RGate 上

https://www.researchgate.net/publication/320197425_Vol1Training_manualRevisedSept2017

第 10 章附录