MEEM(object，conLin，control$niterEM) /固定效果模型矩阵的误差是秩不足的
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)的确，我们可以从dm中了解dy (如果其中一个是1，另一个是0)，但如果这是问题所在，我猜最后一个模型(model_lmer2)仍然会发出警告。
在我的真实数据集中，我最终可以使用time作为数字(虽然不是我的第一种方法)，但是我想了解使用它作为绝对的方法有什么问题。
谢谢!
回答 1
Stack Overflow用户
发布于 2021-12-29 23:27:54
这确实是一个关于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()；这是可能的但这是件讨厌的事。可以在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)或者类似的东西，看看模型矩阵的(头)，模型矩阵的列名，等等。
https://stackoverflow.com/questions/70524441
复制相似问题
腾讯云开发者
