解释95%CI在线性混合效应模型之间的两个电平因子模型估计差异

问题描述 投票:0回答:1

这是我的数据帧(请复制并粘贴复制):

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) 

我然后运行FA值与NLME封装的线性混合模型:

library(nlme)
lmm <- lme(FA ~ Volume + Group, ~ 1|Subject, data = data)
summary(lmm)

现在我想以确定用于“组”因子(控制与病人)的两个电平之间在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,这是没有意义的,因为英足总应该是从0比值为1的事实是,我的因变量是一个比率值的问题?如果是这样,我需要改变以某种方式我的数据(即数变换)来纠正呢?任何反馈将不胜感激!

r statistics mixed-models confidence-interval nlme
1个回答
0
投票

正如@ 42-指出,这里的问题是模型本身。由于FD被限制为[0,1],我们不能使用lme其中假定正常的错误。

Model definition

我不知道你的数据/实验,但也许是一个测试模型可以工作的任何细节;具体来说,我们可以使用形式的不同截距层次模型

enter image description here

在这里我们通过分对数链路FDμ的均值连接到特定Subject截距和所述预测VolumeGroup

Implementation

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

在估计一些评论

请注意,估计是在给定的罗吉特(对数几率)的规模;对于Group = Control估计是那么0.503 - 0.369 * 0 = 0.503,并Group = Patient0.503 - 0.369 * 1 = 0.134Group = PatientGroup = Control(再次在分对数刻度)之间的差别只是对GroupPatient系数是-0.369

Comparison of marginal means

那么我建议使用emmeans任何后续的分析;在这种情况下,我们可以使用emmeans::pairs的估计边际均值(EMM)的两个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

请注意,结果Logit变换的分摊比额表(而不是在响应规模)。要获得FD响应为Group = Patient的比例和Group = Control您需要手动将这些估计。

说明:这里emmeans返回的EMM对Group,和pairs执行不同级别的Group的成对比较。然后,我们使用confint返回(默认95%)置信区间。

这种做法的好处是,你不会有,如果你有> 2个水平Group任何改变; 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
© www.soinside.com 2019 - 2024. All rights reserved.