MEEM(object, conLin, control$niterEM) 中的错误/固定效应模型矩阵秩不足

Question

我正在尝试进行多级中介分析（如 here 和 here 所做的那样）。

library(data.table)
library(lme4)
library(nlme)
library(magrittr)
library(dplyr)

set.seed(1)

# Simulated data ------------------------------------------------------------------
dt_1 <- data.table(id = rep(1:10, each=4),
                   time = as.factor(rep(1:4, 10)),
                   x = rnorm(40),
                   m = rnorm(40),
                   y = rnorm(40))

# Melt m and y into z ------------------------------------------------------------------
dt_2 <- dt_1 %>%
  mutate(mm = m, .after=x) %>%
  melt(id.vars = c("id","time","x","mm"),
       na.rm = F,
       variable.name = "dv",
       value.name = "z") %>%
  within({
    dy <- as.integer(dv == "y")
    dm <- as.integer(dv == "m")
    }) %>%
  arrange(id,time)

> head(dt_2,4)
   id time          x         mm dv          z dm dy
1:  1    1 -0.6264538 -0.1645236  m -0.1645236  1  0
2:  1    1 -0.6264538 -0.1645236  y -0.5686687  0  1
3:  1    2  0.1836433 -0.2533617  m -0.2533617  1  0
4:  1    2  0.1836433 -0.2533617  y -0.1351786  0  1

# lme mediation model ------------------------------------------------------------------
model_lme <- lme(fixed = z ~ 0 + 
                             dm + dm:x + dm:time + #m as outcome
                             dy + dy:mm + dy:x + dy:time, #y as outcome
                 random = ~ 0  +  dm:x + dy:mm + dy:x | id, 
                 weights = varIdent(form = ~ 1 | dm), #separate sigma^{2}_{e} for each outcome
                 data = dt_2,
                 na.action = na.exclude)

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

# lmer mediation model ------------------------------------------------------------------
model_lmer <- lmer(z ~ 0 + dm + dm:x + dm:time + dy + dy:mm + dy:x + dy:time +
                     (0  +  dm:x + dy:mm + dy:x | id) + (0 + dm | time), 
                  data = dt_2,
                  na.action = na.exclude)

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

我看过一些关于这些错误 (nlme) / 警告 (lme4) 的帖子（例如，this 和 this），但我没有弄清楚这里有什么问题。

我检查过

X <- model.matrix(~0 + dm + dm:x + dm:time + dy + dy:mm + dy:x + dy:time, data=dt_2)

> caret::findLinearCombos(X)
$linearCombos
$linearCombos[[1]]
[1] 7 1 4 5 6

$remove
[1] 7

但我不太明白输出。

从model_lmer 的摘要中，我验证了dm:time4 和time1:dy 系数丢失，但为什么呢？数据集中有所有可能的组合（0/0、0/1、1/0、1/1）。

Fixed effects:
         Estimate Std. Error t value
dm        0.30898    0.92355   0.335
dy        0.03397    0.27480   0.124
dm:x      0.21267    0.19138   1.111
dm:time1 -0.19713    1.30583  -0.151
dm:time2 -0.30206    1.30544  -0.231
dm:time3 -0.20828    1.30620  -0.159
dy:mm     0.22625    0.18728   1.208
x:dy     -0.37747    0.17130  -2.204
time2:dy  0.29894    0.39123   0.764
time3:dy  0.22640    0.39609   0.572
time4:dy -0.16758    0.39457  -0.425

另一方面，使用time 作为数字不会给出错误/警告：

# lmer mediation model - time as numeric
model_lmer2 <- lmer(z ~ 0 + dm + dm:x + dm:time + dy + dy:mm + dy:x + dy:time +
                     (0  +  dm:x + dy:mm + dy:x | id) + (0 + dm | time), 
                  data = within(dt_2, time <- as.numeric(time)),
                  na.action = na.exclude)

确实可以从dy 知道dm（如果一个是1，另一个是0），但如果这是问题所在，我猜最后一个模型（model_lmer2）仍然会发出警告.

在我的真实数据集中，我最终可以将time 用作数字（虽然不是我的第一种方法），但我想了解将其用作分类有什么问题。

谢谢！

Answer 1

这确实是一个关于 R 中的线性模型构建/公式的通用问题：它不是特定于混合模型的。

让我们看看变量的多重共线性组合所涉及的列的名称（即第 7、1、4、5、6 列）：

cc <- caret::findLinearCombos(X)
colnames(X)[cc$linearCombos[[1]]]
## [1] "dm:time4" "dm"       "dm:time1" "dm:time2" "dm:time3"

这告诉我们dm 的主要作用与dm:time 交互混淆了；一旦我们知道dm:time[i] 对i 的所有级别，dm 的主要作用是多余的。

太糟糕了lme 没有像lmer 那样自动删除列来调整多重共线性，而且lmer 没有超级方便的方法来模拟异方差性 à la varIdent();有可能but it's a nuisance。 可能在nlme 或glmmTMB 中构建自动投放（这也可以轻松地对异方差建模），但目前还没有人使用它）。

...如果您可以指定 dm:time 并将 dm 排除在模型之外，那么这可能是最简单的！

---

您可以试验不同型号规格会发生什么：

lc <- function(f) {
    X <- model.matrix(f, dt_2)
    cc <- caret::findLinearCombos(X)
    lapply(cc$linearCombos, function(z) colnames(X)[z])
}

lc(~0 + dm + dm:time)
lc(~0 + dy + dy:time)
lc(~0 + dm + dm:time + dy + dy:time)
lc(~0 + dy + dy:time + dm + dm:time)

或类似的东西查看模型矩阵的（头部）、模型矩阵的列名等。