线性混合模型置信区间问题
Linear mixed model confidence intervals question
希望你能解开我脑海中的一些困惑。
线性混合模型构造为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, mean 76.4, SE:8.96）。
我遗漏了什么，我应该报告什么数字？
总而言之，对于 Ca 的含量，我每次处理总共有 6 个测量值（每年三个，每个来自不同的复制）。我将使用我的语言在代码中保留名称。想法是测试某些生产实践是否会影响谷物中特定矿物质的含量。此示例的模型中保留了没有残差的随机效应。
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 校正，这或多或少是最准确的可用选项 — 不[=48= 的唯一原因] 使用 K-R 的情况是你有一个大数据集，它变得非常慢。
加载包，模拟数据，拟合模型
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 没有明确告诉我们 what/whether 它正在使用有限样本校正：我们可以在文档或代码中四处挖掘以尝试找出答案。或者，我们可以告诉 emmeans not 使用有限样本校正：
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）。这对我来说似乎有问题，除非你的数据集非常小......
根据您更新后的 post，Kenward-Roger 似乎确实只估计了 1.5 分母 df。
一般来说，dicey/not 建议在分组变量具有少量水平的情况下尝试拟合随机效应（尽管参见 here 反驳）。我会尝试将 LETO （只有两个级别）视为固定效果，即
CA ~ SISTEM + LETO + (1 | LETO:ponovitev) + (1 | LETO:SISTEM)
看看是否有帮助。 （我预计你会得到 7 df 的数量级，这将使你的 CIs ± 2.4 SE 而不是 ± 6 SE ...）
希望你能解开我脑海中的一些困惑。
线性混合模型构造为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, mean 76.4, SE:8.96）。
我遗漏了什么，我应该报告什么数字？
总而言之，对于 Ca 的含量，我每次处理总共有 6 个测量值（每年三个，每个来自不同的复制）。我将使用我的语言在代码中保留名称。想法是测试某些生产实践是否会影响谷物中特定矿物质的含量。此示例的模型中保留了没有残差的随机效应。
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 校正，这或多或少是最准确的可用选项 — 不[=48= 的唯一原因] 使用 K-R 的情况是你有一个大数据集，它变得非常慢。
加载包，模拟数据，拟合模型
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 没有明确告诉我们 what/whether 它正在使用有限样本校正：我们可以在文档或代码中四处挖掘以尝试找出答案。或者，我们可以告诉 emmeans not 使用有限样本校正：
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）。这对我来说似乎有问题，除非你的数据集非常小......
根据您更新后的 post，Kenward-Roger 似乎确实只估计了 1.5 分母 df。
一般来说，dicey/not 建议在分组变量具有少量水平的情况下尝试拟合随机效应（尽管参见 here 反驳）。我会尝试将 LETO （只有两个级别）视为固定效果，即
CA ~ SISTEM + LETO + (1 | LETO:ponovitev) + (1 | LETO:SISTEM)
看看是否有帮助。 （我预计你会得到 7 df 的数量级，这将使你的 CIs ± 2.4 SE 而不是 ± 6 SE ...）