我正在尝试进行多级中介分析(如这里和这里所做的(。
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(的帖子(例如，这个和这个(，但我不知道这里的问题是什么。
我检查了
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用作数字(尽管这不是我的第一种方法(，但我想了解将其用作分类的问题是什么。
谢谢！
这实际上是一个关于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的相互作用相混淆；一旦我们知道了所有级别的i的dm:time[i]，则dm的主要作用是多余的。
令人遗憾的是，lme没有像lmer那样自动删除列以调整多重共线性，lmer也没有一种超级方便的方法来对varIdent()的异方差进行建模；这是可能的，但令人讨厌。将自动丢弃构建到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)
或者查看模型矩阵的(标题(、模型矩阵的列名等类似内容