尝试在R中运行lmer()时出现错误

Question

我的问题是这样的。我有一个在R语言中的数据集，需要在其中运行一个混合效应模型。以下是代码：

data <- read.csv("D:/blahblah.csv")
analysis.data <- lmer(intdiff ~ stress_limit * word_position * follows + (1|speaker), data)
summary(analysis.data)

当我尝试运行脚本时，它返回以下错误：

 Error in mer_finalize(ans) : Downdated X'X is not positive definite, 15.

我已经追踪到错误是由于"follows"参数引起的，因为当我只使用stress_limit和word_position时，它可以正常运行。如果有帮助的话，"follows"中的数据只有三个字符串：n或l、辅音、元音。我尝试用下划线替换空格，但没有成功。在这种情况下，lmer()函数的内部工作是否有阻止使用"follows"的原因？任何帮助都将是极好的！

更多信息：intdiff包含数值，stress_limit是字符串（Stressed或Unstressed），单词位置也是字符串（Word Medial或Word Initial）。

编辑：这里是一个可重现错误的数据样本：

structure(list(intdiff = c(11.45007951, 12.40144758, 13.47898367, 
6.279497762, 18.19461897, 16.15539707), word_position = structure(c(2L, 
2L, 2L, 1L, 1L, 1L), .Label = c("Word Initial", "Word Medial"
), class = "factor"), follows = structure(c(4L, 4L, 4L, 1L, 2L, 
4L), .Label = c("Consonant", "n or l", "Pause", "Vowel"), class = "factor"), 
stress_limit = structure(c(2L, 1L, 1L, 2L, 2L, 2L), .Label = c("Stressed", 
"Unstressed"), class = "factor"), speaker = structure(c(2L, 
2L, 2L, 2L, 2L, 2L), .Label = c("f11r", "f13r", "f15a", "f16a", 
"m09a", "m10a", "m12r", "m14r"), class = "factor")), .Names = c("intdiff", 
"word_position", "follows", "stress_limit", "speaker"), row.names = c(NA, 
6L), class = "data.frame")

我也尝试了lme()函数，但是返回了以下错误：

Error in MEEM(object, conLin, control$niterEM) : 
Singularity in backsolve at level 0, block 1

我原始帖子中的代码是我正在使用的确切代码，不包含库(lme4)调用，因此我想不出遗漏任何信息。

我使用的R版本是2.15.2

- Shakesbeery

---

你的实际数据框架 data 有多少行？ - Sven Hohenstein

数据框大约有1110行。但是，如果给出一个样本，数据是可以预测的。 - Shakesbeery

你的预测变量是否包含应力限制、单词位置和后续所有可能的组合，还是有一些缺失（因为它们不可行或者你没有测量到）？with(data,all(table(stress_limit,word_position,follows)>=1)) 是否为真？（这似乎已经变成了一个统计问题而不是一个编程问题...） - Ben Bolker

结果是假的，所以我猜确实并非所有可能的组合都存在。这是否意味着完全无法运行混合效应模型？或者当我运行函数时可以考虑到这一点吗？ - Shakesbeery

Answer 1

没有可重复的例子很难确定： 如何创建一个完美的 R 可重复例子？

但是，猜测：这些问题通常是由于设计矩阵中的共线性导致的。对连续预测变量（intdiff）进行居中可能有所帮助。您还可以直接探索设计矩阵。

X <- model.matrix( ~ stress_limit * word_position * follows, data)

两个变量之间的共线性： cor（X）。不幸的是，我暂时没有关于检测多重共线性（即不仅限于对之间的共线性，而是多于两个预测变量之间的共线性）的建议，不过你可以查看用于计算方差膨胀因子的工具（例如library("sos"); findFn("VIF")）。

作为交叉验证，lme也应该能够处理你的模型：

library(nlme)
lme(intdiff ~ stress_limit * word_position * follows, 
   random=~1|speaker, data=data)

当我在lme4的开发版本中运行你的测试数据（在github上可用）时，我得到了Error in lmer(intdiff ~ stress_limit * word_position * follows + (1 | : rank of X = 5 < ncol(X) = 12。然而，对于这么小的输入数据集（6个观测值），不可能拟合12个参数。很难确定你的问题出在哪里。你的3个变量的所有12种组合实际上都出现在你的数据中吗？如果有缺失的话，那么你需要遵循开发版本的帮助中给出的建议：

与一些更简单的建模框架（例如“lm”和“glm”）自动检测完全共线的预测变量不同，“[gn]lmer”不能处理不满秩的设计矩阵。例如，在具有未观察到的交互水平组合的模型的情况下，由用户定义一个新的变量是用户的责任（例如，在从“droplevels( interaction(a,b) )”结果创建的数据中创建“ab”）。

特别地，你可以按以下方式拟合该模型：

data <- transform(data,
       allcomb=interaction(stress_limit,word_position,follow,drop=TRUE))
lme(intdiff ~ allcomb, random=~1|speaker, data=data)

这将为您提供一种单向方差分析（ANOVA），将实际存在于数据中的唯一级别组合视为类别。 您需要自行弄清它们的含义。

另一种方法是减少模型中的交互作用，直到您达到没有任何缺失组合的一组为止；如果您很幸运，(stress_limit+word_position+follow)^2（所有二元交互作用）将会起作用，但您可能还需要进一步缩小模型（例如，stress_limit + word_position*follow）。

另一种测试方法是在您的拟议模型上使用lm()，并检查估计系数中是否有NA值。

这些方法失去的主要是方便性/易解释性，因为无法从数据中估计缺失组合的参数...

- Ben Bolker

---

首先，非常感谢您的回答！我不太理解问题中共线性方面，也许您可以解释一下吗？我在谷歌上搜索了这个术语，但并没有起到很好的澄清作用（因为我是一名初学统计学的人）。除此之外，我已经改进了上面提供的信息，也许可以帮助我们解决这个问题？ - Shakesbeery

@Shakesbeery 是的，这就是难点。工具只能做数据允许的事情：-(。这让我想起了我第一次接触最小二乘分析。我欣然地想到（减少本科物理实验数据），即使底层模型是y ~x^2，我也可以得到一个更“酷”的8阶多项式拟合。结果失败了。这种事情发生在我们所有人身上。 - Carl Witthoft

太棒了，这是一条非常有帮助和信息丰富的回复！我想我可以从这里开始自己解决问题了。再次感谢您的帮助。 - Shakesbeery