如何解决线性混合模型置信区间问题
希望你能消除我脑海中的一些困惑。
线性混合模型由lmertest构建：
MODEL <- lmer(Ca content ~ SYstem +(1 | YEAR/replicate) +
(1 | YEAR:SYstem),data = IOSDV1)
当我试图获得主效应特定水平的置信区间时，乐趣开始出现。
命令 emmeans 和 lsmeans 产生相同的间隔（示例；SYstem A3: 23.9-128.9,mean 76.4,SE:8.96）。
但是，命令 as.data.frame(effect("SYstem",MODEL)) 会产生不同的、更窄的置信区间（示例；SYstem A3: 58.0-94.9,SE:8.96）。
我遗漏了什么，我应该报告什么数字？
总而言之，对于 Ca 的含量，我每次处理有 6 次总测量值（每年 3 次，每次来自不同的重复）。我将使用我的语言在代码中保留名称，如使用。想法是测试某些生产实践是否会影响谷物中特定矿物质的含量。本例的模型中保留了无残差的随机效应。
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModlmertest']
Formula: CA ~ SIstem + (1 | Leto/ponovitev) + (1 | Leto:SIstem)
Data: IOSDV1
REML criterion at convergence: 202.1
Scaled residuals:
Min 1Q Median 3Q Max
-1.60767 -0.74339 0.04665 0.73152 1.50519
Random effects:
Groups Name Variance Std.Dev.
Leto:SIstem (Intercept) 0.0 0.0
ponovitev:Leto (Intercept) 0.0 0.0
Leto (Intercept) 120.9 11.0
Residual 118.7 10.9
Number of obs: 30,groups: Leto:SIstem,10; ponovitev:Leto,8; Leto,2
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 76.417 8.959 1.548 8.530 0.0276 *
SIstem[T.C0] -5.183 6.291 24.000 -0.824 0.4181
SIstem[T.C110] -13.433 6.291 24.000 -2.135 0.0431 *
SIstem[T.C165] -7.617 6.291 24.000 -1.211 0.2378
SIstem[T.C55] -10.883 6.291 24.000 -1.730 0.0965 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) SIstem[T.C0 SIstem[T.C11 SIstem[T.C16
SIstem[T.C0 -0.351
SIstem[T.C11 -0.351 0.500
SIstem[T.C16 -0.351 0.500 0.500
SIstem[T.C5 -0.351 0.500 0.500 0.500
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular
> ls_means(MODEL,ddf="Kenward-Roger")
Least Squares Means table:
Estimate Std. Error df t value lower upper Pr(>|t|)
SIstemA3 76.4167 8.9586 1.5 8.5299 23.9091 128.9243 0.02853 *
SIstemC0 71.2333 8.9586 1.5 7.9514 18.7257 123.7409 0.03171 *
SIstemC110 62.9833 8.9586 1.5 7.0305 10.4757 115.4909 0.03813 *
SIstemC165 68.8000 8.9586 1.5 7.6797 16.2924 121.3076 0.03341 *
SIstemC55 65.5333 8.9586 1.5 7.3151 13.0257 118.0409 0.03594 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Confidence level: 95%
degrees of freedom method: Kenward-Roger
> emmeans(MODEL,spec = c("SIstem"))
SIstem emmean SE df lower.CL upper.CL
A3 76.4 8.96 1.53 23.9 129
C0 71.2 8.96 1.53 18.7 124
C110 63.0 8.96 1.53 10.5 115
C165 68.8 8.96 1.53 16.3 121
C55 65.5 8.96 1.53 13.0 118
degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
> as.data.frame(effect("SIstem",MODEL))
SIstem fit se lower upper
1 A3 76.41667 8.958643 57.96600 94.86734
2 C0 71.23333 8.958643 52.78266 89.68400
3 C110 62.98333 8.958643 44.53266 81.43400
4 C165 68.80000 8.958643 50.34933 87.25067
5 C55 65.53333 8.958643 47.08266 83.98400
非常感谢。
解决方法
我很确定这与可怕的“分母自由度”问题有关，即正在使用哪种（如果有）有限样本校正。 tl;dr emmeans 正在使用 Kenward-Roger 校正，这或多或少是最准确的可用选项——不使用 KR 的唯一原因是如果你有一个大数据集，它变得慢得难以忍受。
加载包，模拟数据，拟合模型
library(lmerTest)
library(emmeans)
library(effects)
dd <- expand.grid(f=factor(letters[1:3]),g=factor(1:20),rep=1:10)
set.seed(101)
dd$y <- simulate(~f+(1|g),newdata=dd,newparams=list(beta=rep(1,3),theta=1,sigma=1))[[1]]
m <- lmer(y~f+(1|g),data=dd)
比较默认的emmeans和效果
emmeans(m,~f)
## f emmean SE df lower.CL upper.CL
## a 0.848 0.212 21.9 0.409 1.29
## b 1.853 0.212 21.9 1.414 2.29
## c 1.863 0.212 21.9 1.424 2.30
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
as.data.frame(effect("f",m))
## f fit se lower upper
## 1 a 0.8480161 0.2117093 0.4322306 1.263802
## 2 b 1.8531805 0.2117093 1.4373950 2.268966
## 3 c 1.8632228 0.2117093 1.4474373 2.279008
effects 没有明确告诉我们什么/是否使用有限样本校正：我们可以在文档或代码中挖掘以试图找出答案。或者，我们可以告诉emmeans 不要使用有限样本校正：
emmeans(m,~f,lmer.df="asymptotic")
## f emmean SE df asymp.LCL asymp.UCL
## a 0.848 0.212 Inf 0.433 1.26
## b 1.853 0.212 Inf 1.438 2.27
## c 1.863 0.212 Inf 1.448 2.28
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
测试表明，这些大约相当于 0.001 的容差（可能足够接近）。 原则上我们应该能够指定 KR=TRUE 来让 effects 使用 Kenward-Roger 校正，但我还没有能够让它工作。 >
但是，我还要说您的示例有点奇怪。如果我们以标准误差为单位计算均值和下 CI 之间的距离，对于 emmeans，我们得到 (76.4-23.9)/8.96 = 5.86，这意味着 非常 效应自由度（例如约 1.55）。这对我来说似乎有问题，除非您的数据集非常小......
从您更新的帖子来看，Kenward-Roger 似乎确实只估计了 1.5 分母 df。
一般来说，在分组变量具有少量级别的情况下尝试拟合随机效应是冒险的/不推荐的（尽管请参阅 here 以获得反驳）。我会尝试将 LETO（只有两个级别）视为固定效果，即
CA ~ SISTEM + LETO + (1 | LETO:ponovitev) + (1 | LETO:SISTEM)
看看是否有帮助。 （我预计你会得到 7 df 的数量级，这将使你的 CI ± 2.4 SE 而不是 ± 6 SE ...）
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。