更改随机效应分组变量的标签会更改 lme4 中的结果

如何解决更改随机效应分组变量的标签会更改 lme4 中的结果

标题说明了一切：更改随机效应分组变量的（据称是任意的）标签（例如，重复测量实验中的对象名称）可以更改 lme4 中的结果输出。最小示例：

require(dplyr)
require(lme4)
require(digest)
df = faithful %>% mutate(subject = rep(as.character(1:8),each = 34),subject2 = rep(as.character(9:16),each = 34))
summary(lmer(eruptions ~ waiting + (waiting | subject),data = df))$coefficients[2,1] # = 0.07564181
summary(lmer(eruptions ~ waiting + (waiting | subject2),1] # = 0.07567655

我认为这是因为 lme4 将它们转换为因子，并且不同的名称会产生不同的因子级别排序。例如。这会产生问题：

df2 = faithful %>% mutate(subject = factor(rep(as.character(1:8),each = 34)),subject2 = factor(rep(as.character(9:16),each = 34)))
summary(lmer(eruptions ~ waiting + (waiting | subject),data = df2))$coefficients[2,1] # = 0.07567655

但这不会：

df3 = faithful %>% mutate(subject = factor(rep(as.character(1:8),subject2 = factor(rep(as.character(1:8),levels = as.character(1:8),labels = as.character(9:16)))
summary(lmer(eruptions ~ waiting + (waiting | subject),data = df3))$coefficients[2,1] # = 0.07564181

这似乎是 lme4 中的一个问题。不同的任意变量标签不应该产生不同的输出，对吗？我错过了什么吗？为什么 lme4 会这样做？

（我知道输出的差异很小，但在其他情况下我得到了更大的差异，足以例如将 ap 值从 .055 更改为 .045。此外，如果这是正确的，我认为这可能会导致轻微的可重复性问题——例如，如果在完成分析后，实验者将他们的人类受试者数据匿名化（通过更改名称），然后将其发布到公共存储库中。）

解决方法

当我拟合这个模型时，我会在拟合时收到一个单一的拟合警告。这不是一个好兆头，因为仅由随机截距解释的方差实际上为 0，而且您还有一个随机斜率。这里的随机效应可能对模型没有任何意义。

其次，我质疑这是否适合这种情况，以下是主动提供的建议，如果您认为不合适，我深表歉意。其次，我会对此发表评论，但不确定如何添加图片。

首先，我做了一些探索性绘图，发现您的因变量和固定效应都具有双峰分布。如果我们绘制如下散点图，我们肯定可以看到它可能不是线性趋势。

然后当我们查看模型残差时，我们会看到异方差性，这是次优的。我不是统计学家，但一些顾问告诉我，这是线性模型中最糟糕的假设之一。

我认为，由于奇异拟合，您可能会看到估计值的不稳定性，但希望其他人能够提供更多信息来解决这个问题。

,

序列 1:8 的第一部分以数字或字符格式给出相同的 order，而第二部分没有：

identical(order(1:8),order(as.character(1:8)))
# [1] TRUE
identical(order(9:16),order(as.character(9:16)))
# [1] FALSE

那是因为数字作为字符按第一个数字排序：

sort(9:16)
# [1]  9 10 11 12 13 14 15 16
sort(as.character(9:16))
# [1] "10" "11" "12" "13" "14" "15" "16" "9" 

因此，如果您使用两个不同但只有一位的字符序列，则似乎没有问题：

library(lme4)
fo1 <- eruptions ~ waiting + (waiting | sub)
fo2 <- eruptions ~ waiting + (waiting | sub2)

df1 <- transform(faithful,sub=rep(as.character(1:8),each=34),sub2=rep(as.character(2:9),each=34))

summary(lmer(fo1,data=df1))$coe[2,1]
# boundary (singular) fit: see ?isSingular
# [1] 0.07564181
summary(lmer(fo2,1]
# boundary (singular) fit: see ?isSingular
# [1] 0.07564181

然而，您的分组变量的顺序在 lmer() 中确实很重要。这可以通过给 subject 和 subject2 相同的级别但不同的顺序来显示：

set.seed(840947)
df2 <- transform(faithful,sub=rep(sample(1:8),sub2=rep(sample(1:8),each=34))

summary(fit2a <- lmer(fo1,data=df2))$coe[2,1]
# boundary (singular) fit: see ?isSingular
# [1] 0.07564179
summary(fit2b <- lmer(fo2,1]
# boundary (singular) fit: see ?isSingular
# [1] 0.07567537

这会再次产生完全不同的系数。可以像这样检查级别和级别顺序：

fit2a@flist$sub
# [1] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
# [33] 4 4 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
# [65] 8 8 8 8 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
# [97] 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
# [129] 3 3 3 3 3 3 3 3 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
# [161] 6 6 6 6 6 6 6 6 6 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
# [193] 1 1 1 1 1 1 1 1 1 1 1 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
# [225] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
# [257] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
# Levels: 1 2 3 4 5 6 7 8

fit2b@flist$sub2
# [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
# [33] 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
# [65] 2 2 2 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
# [97] 6 6 6 6 6 6 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
# [129] 8 8 8 8 8 8 8 8 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
# [161] 7 7 7 7 7 7 7 7 7 7 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
# [193] 3 3 3 3 3 3 3 3 3 3 3 3 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
# [225] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
# [257] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
# Levels: 1 2 3 4 5 6 7 8

已经有一个 ticket filed at github 您应该加入。也许尝试事先找到一个类似的案例，其中存在排序问题，但不是单一拟合。