解释线性混合效应模型中两个水平因子之间模型估计差异的 95% CI
Interpreting 95% CI for model estimated difference between two level factor in linear mixed-effects model
这是我的数据框（请复制粘贴复制）：
Control <- replicate(2, c("112", "113", "116", "118", "127", "131", "134", "135", "136", "138", "143", "148", "149", "152", "153", "155", "162", "163"))
EPD <- replicate(2, c("101", "102", "103", "104", "105", "106", "107", "108", "109", "110", "114", "115", "117", "119", "120", "122", "124", "125", "126", "128", "130", "133", "137", "139", "140", "141", "142", "144", "145", "147"))
Subject <- c(Control, EPD)
Control_FA_L <- c(0.43, 0.39, 0.38, 0.58, 0.37, 0.5, 0.35, 0.36, 0.72, 0.38, 0.45, 0.30, 0.47, 0.30, 0.67, 0.34, 0.42, 0.29)
Control_FA_R <- c(0.36, 0.49, 0.55, 0.59, 0.33, 0.41, 0.32, 0.50, 0.59, 0.52, 0.32, 0.40, 0.49, 0.33, 0.46, 0.39, 0.37, 0.33)
EPD_FA_L <- c(0.25, 0.39, 0.36, 0.42, 0.21, 0.40, 0.43, 0.16, 0.31, 0.41, 0.39, 0.40, 0.35, 0.29, 0.31, 0.24, 0.39, 0.36, 0.54, 0.38, 0.34, 0.28, 0.42, 0.33, 0.40, 0.36, 0.42, 0.28, 0.40, 0.41)
EPD_FA_R <- c(0.26, 0.36, 0.36, 0.61, 0.22, 0.33, 0.36, 0.34, 0.35, 0.37, 0.39, 0.45, 0.30, 0.31, 0.50, 0.31, 0.29, 0.43, 0.41, 0.21, 0.38, 0.28, 0.66, 0.33, 0.50, 0.27, 0.46, 0.37, 0.26, 0.39)
FA <- c(Control_FA_L, Control_FA_R, EPD_FA_L, EPD_FA_R)
Control_Volume_L <- c(99, 119, 119, 146, 127, 96, 100, 132, 103, 103, 107, 142, 140, 134, 117, 117, 133, 143)
Control_Volume_R <- c(93, 123, 114, 152, 122, 105, 98, 138, 111, 110, 115, 137, 142, 140, 124, 102, 153, 143)
EPD_Volume_L <- c(132, 115, 140, 102, 130, 131, 110, 124, 102, 111, 93, 92, 94, 104, 92, 115, 144, 118, 104, 132, 90, 102, 94, 112, 106, 105, 79, 114, 104, 108)
EPD_Volume_R <- c(136, 116, 143, 105, 136, 137, 103, 121, 105, 115, 97, 97, 93, 108, 91, 117, 147, 111, 97, 129, 85, 107, 91, 116, 113, 101, 75, 108, 95, 98)
Volume <- c(Control_Volume_L, Control_Volume_R, EPD_Volume_L, EPD_Volume_R)
Group <- c(replicate(36, "Control"), replicate(60, "Patient"))
data <- data.frame(Subject, FA, Volume, Group)
然后我 运行 使用 nlme 包的 FA 值的线性混合模型：
library(nlme)
lmm <- lme(FA ~ Volume + Group, ~ 1|Subject, data = data)
summary(lmm)
我现在想确定 "Group" 因子（对照和患者）两个水平之间 FA 模型估计差异的 95% 置信区间。我通常会通过执行以下代码来继续：
# Compute 95% Confidence Interval for Group factor
# True difference in STN FA between Control and EPD subjects
0.0857851 # Value from mixed model
# Multiply 97.5 percentile point of normal distribution by std error from mixed model
1.96 * 0.02555076 # 95% CI: 0.086 ± 0.050 mm^3 (p = .0016) - !!CI includes values > 1!!
我很难理解这意味着什么。我计算的置信区间包括大于 1 的值，这没有意义，因为 FA 应该是从 0 到 1 的比率值。我的因变量是比率值这一事实是问题所在吗？如果是这样，我是否需要以某种方式转换我的数据（即日志转换）来纠正这个问题？任何反馈将不胜感激！
正如@42-所指出的，这里的问题在于模型本身。由于 FD 被限制为 [0, 1] 我们不能使用假定正常错误的 lme。
模型定义
我不知道关于您的 data/experiment 的任何细节，但 Beta 模型可能会起作用；具体来说，我们可以使用
形式的变截距层次模型
我们通过 logit-link 将 FD μ 的均值连接到特定于 Subject 的截距和预测变量 Volume 和 Group .
实施
库glmmTMB允许实现这样一个混合效应模型
library(glmmTMB)
lmm <- glmmTMB(
FA ~ Volume + Group + (1 | Subject),
data = data,
family = "beta_family")
summary(lmm)$coef$cond
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) 0.502858259 0.348506927 1.442893 0.1490505719
#Volume -0.006464251 0.002782781 -2.322947 0.0201820253
#GroupPatient -0.369273205 0.104832100 -3.522520 0.0004274642
对估算的一些评论
请注意，估计值是在 logit（对数几率）尺度上给出的； Group = Control 的估计值为 0.503 - 0.369 * 0 = 0.503，而 Group = Patient 的估计值为 0.503 - 0.369 * 1 = 0.134。 Group = Patient 和 Group = Control 之间的 差异 （同样在对数尺度上）只是 GroupPatient 的系数，即 -0.369。
边际均值比较
然后我建议使用 emmeans 进行任何后续分析；在这种情况下，我们可以使用 emmeans::pairs 来比较两个 Group 水平
的估计边际均值 (EMM)
library(emmeans)
confint(pairs(emmeans(lmm, "Group")))
# contrast estimate SE df lower.CL upper.CL
# Control - Patient 0.3692732 0.1048321 91 0.1610371 0.5775093
#
#Results are given on the log odds ratio (not the response) scale.
#Confidence level used: 0.95
请注意，结果是在逻辑尺度上给出的（而不是在响应尺度上）。要获得 Group = Patient 和 Group = Control 的 FD 响应比率，您需要手动转换这些估计值。
解释：这里emmeansreturn是Group的EMM，pairs对Group的不同水平进行了成对比较。然后我们使用 confint 到 return（默认 95%）的置信区间。
好处是，如果 Group 的级别超过 2 个，则无需更改任何内容； pairs 将执行成对比较，并自动更正 p 值以进行多重假设检验。
如需更多阅读，请查看优秀的插图 Comparisons and contrasts in emmeans。
您还可以在比值比尺度上获得估计的边际均值和置信区间（这避免了必须手动从对数比尺度转换为比值比尺度）
confint(pairs(emmeans(lmm, "Group"), type = "response"))
# contrast odds.ratio SE df lower.CL upper.CL
# Control / Patient 1.446683 0.1516588 91 1.174729 1.781595
Confidence level used: 0.95
Intervals are back-transformed from the log odds ratio scale
这是我的数据框（请复制粘贴复制）：
Control <- replicate(2, c("112", "113", "116", "118", "127", "131", "134", "135", "136", "138", "143", "148", "149", "152", "153", "155", "162", "163"))
EPD <- replicate(2, c("101", "102", "103", "104", "105", "106", "107", "108", "109", "110", "114", "115", "117", "119", "120", "122", "124", "125", "126", "128", "130", "133", "137", "139", "140", "141", "142", "144", "145", "147"))
Subject <- c(Control, EPD)
Control_FA_L <- c(0.43, 0.39, 0.38, 0.58, 0.37, 0.5, 0.35, 0.36, 0.72, 0.38, 0.45, 0.30, 0.47, 0.30, 0.67, 0.34, 0.42, 0.29)
Control_FA_R <- c(0.36, 0.49, 0.55, 0.59, 0.33, 0.41, 0.32, 0.50, 0.59, 0.52, 0.32, 0.40, 0.49, 0.33, 0.46, 0.39, 0.37, 0.33)
EPD_FA_L <- c(0.25, 0.39, 0.36, 0.42, 0.21, 0.40, 0.43, 0.16, 0.31, 0.41, 0.39, 0.40, 0.35, 0.29, 0.31, 0.24, 0.39, 0.36, 0.54, 0.38, 0.34, 0.28, 0.42, 0.33, 0.40, 0.36, 0.42, 0.28, 0.40, 0.41)
EPD_FA_R <- c(0.26, 0.36, 0.36, 0.61, 0.22, 0.33, 0.36, 0.34, 0.35, 0.37, 0.39, 0.45, 0.30, 0.31, 0.50, 0.31, 0.29, 0.43, 0.41, 0.21, 0.38, 0.28, 0.66, 0.33, 0.50, 0.27, 0.46, 0.37, 0.26, 0.39)
FA <- c(Control_FA_L, Control_FA_R, EPD_FA_L, EPD_FA_R)
Control_Volume_L <- c(99, 119, 119, 146, 127, 96, 100, 132, 103, 103, 107, 142, 140, 134, 117, 117, 133, 143)
Control_Volume_R <- c(93, 123, 114, 152, 122, 105, 98, 138, 111, 110, 115, 137, 142, 140, 124, 102, 153, 143)
EPD_Volume_L <- c(132, 115, 140, 102, 130, 131, 110, 124, 102, 111, 93, 92, 94, 104, 92, 115, 144, 118, 104, 132, 90, 102, 94, 112, 106, 105, 79, 114, 104, 108)
EPD_Volume_R <- c(136, 116, 143, 105, 136, 137, 103, 121, 105, 115, 97, 97, 93, 108, 91, 117, 147, 111, 97, 129, 85, 107, 91, 116, 113, 101, 75, 108, 95, 98)
Volume <- c(Control_Volume_L, Control_Volume_R, EPD_Volume_L, EPD_Volume_R)
Group <- c(replicate(36, "Control"), replicate(60, "Patient"))
data <- data.frame(Subject, FA, Volume, Group)
然后我 运行 使用 nlme 包的 FA 值的线性混合模型：
library(nlme)
lmm <- lme(FA ~ Volume + Group, ~ 1|Subject, data = data)
summary(lmm)
我现在想确定 "Group" 因子（对照和患者）两个水平之间 FA 模型估计差异的 95% 置信区间。我通常会通过执行以下代码来继续：
# Compute 95% Confidence Interval for Group factor
# True difference in STN FA between Control and EPD subjects
0.0857851 # Value from mixed model
# Multiply 97.5 percentile point of normal distribution by std error from mixed model
1.96 * 0.02555076 # 95% CI: 0.086 ± 0.050 mm^3 (p = .0016) - !!CI includes values > 1!!
我很难理解这意味着什么。我计算的置信区间包括大于 1 的值，这没有意义，因为 FA 应该是从 0 到 1 的比率值。我的因变量是比率值这一事实是问题所在吗？如果是这样，我是否需要以某种方式转换我的数据（即日志转换）来纠正这个问题？任何反馈将不胜感激！
正如@42-所指出的，这里的问题在于模型本身。由于 FD 被限制为 [0, 1] 我们不能使用假定正常错误的 lme。
模型定义
我不知道关于您的 data/experiment 的任何细节，但 Beta 模型可能会起作用；具体来说，我们可以使用
形式的变截距层次模型我们通过 logit-link 将 FD μ 的均值连接到特定于 Subject 的截距和预测变量 Volume 和 Group .
实施
库glmmTMB允许实现这样一个混合效应模型
library(glmmTMB)
lmm <- glmmTMB(
FA ~ Volume + Group + (1 | Subject),
data = data,
family = "beta_family")
summary(lmm)$coef$cond
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) 0.502858259 0.348506927 1.442893 0.1490505719
#Volume -0.006464251 0.002782781 -2.322947 0.0201820253
#GroupPatient -0.369273205 0.104832100 -3.522520 0.0004274642
对估算的一些评论
请注意，估计值是在 logit（对数几率）尺度上给出的； Group = Control 的估计值为 0.503 - 0.369 * 0 = 0.503，而 Group = Patient 的估计值为 0.503 - 0.369 * 1 = 0.134。 Group = Patient 和 Group = Control 之间的 差异 （同样在对数尺度上）只是 GroupPatient 的系数，即 -0.369。
边际均值比较
然后我建议使用 emmeans 进行任何后续分析；在这种情况下，我们可以使用 emmeans::pairs 来比较两个 Group 水平
library(emmeans)
confint(pairs(emmeans(lmm, "Group")))
# contrast estimate SE df lower.CL upper.CL
# Control - Patient 0.3692732 0.1048321 91 0.1610371 0.5775093
#
#Results are given on the log odds ratio (not the response) scale.
#Confidence level used: 0.95
请注意，结果是在逻辑尺度上给出的（而不是在响应尺度上）。要获得 Group = Patient 和 Group = Control 的 FD 响应比率，您需要手动转换这些估计值。
解释：这里emmeansreturn是Group的EMM，pairs对Group的不同水平进行了成对比较。然后我们使用 confint 到 return（默认 95%）的置信区间。
好处是，如果 Group 的级别超过 2 个，则无需更改任何内容； pairs 将执行成对比较，并自动更正 p 值以进行多重假设检验。
如需更多阅读，请查看优秀的插图 Comparisons and contrasts in emmeans。
您还可以在比值比尺度上获得估计的边际均值和置信区间（这避免了必须手动从对数比尺度转换为比值比尺度）
confint(pairs(emmeans(lmm, "Group"), type = "response"))
# contrast odds.ratio SE df lower.CL upper.CL
# Control / Patient 1.446683 0.1516588 91 1.174729 1.781595
Confidence level used: 0.95
Intervals are back-transformed from the log odds ratio scale