【发布时间】：2021-02-10 18:39:13
【问题描述】：
我正在使用 R 中的 GLMMadaptive 包运行障碍对数正态模型。连续部分和零部分都在固定效应中定义了分类变量。我想对这些分类变量进行方差分析，以检测是否存在主效应。
我已经看到，使用 glmmTMB 包，您可以分别在条件模型和零部分模型上分别运行 ANOVA，如 here 所示。
对于 GLMMadaptive 包是否有类似的策略可用？ （据我所知，glmmTMB 不支持障碍对数正态模型）。也许使用 emmeans 包中的 joint_tests 函数？如果是这样，您如何定义要测试零部分模型？因为emmeans::joint_tests(hurdlemodel) 只给出了模型条件部分的 F 检验。
或者作为一种替代方法，您能否将排除感兴趣变量的模型与完整模型的拟合度进行比较，正如this vignette 中随机效应的相关性所证明的那样？
非常感谢！
Russ Lenth 在 cmets 中的建议在下面实现，使用 GLMMadaptive two-part model vignette 中的数据和模型：
library(GLMMadaptive)
library(emmeans)
# data generating code from the vignette:
{
set.seed(1234)
n <- 100 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 5 # maximum follow-up time
# we construct a data frame with the design:
# everyone has a baseline measurement, and then measurements at random follow-up times
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))
# design matrices for the fixed and random effects non-zero part
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ 1, data = DF)
# design matrices for the fixed and random effects zero part
X_zi <- model.matrix(~ sex, data = DF)
Z_zi <- model.matrix(~ 1, data = DF)
betas <- c(1.5, 0.05, 0.05, -0.03) # fixed effects coefficients non-zero part
shape <- 2 # shape/size parameter of the negative binomial distribution
gammas <- c(-1.5, 0.5) # fixed effects coefficients zero part
D11 <- 0.5 # variance of random intercepts non-zero part
D22 <- 0.4 # variance of random intercepts zero part
# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))
# linear predictor non-zero part
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, 1, drop = FALSE]))
# linear predictor zero part
eta_zi <- as.vector(X_zi %*% gammas + rowSums(Z_zi * b[DF$id, 2, drop = FALSE]))
# we simulate negative binomial longitudinal data
DF$y <- rnbinom(n * K, size = shape, mu = exp(eta_y))
# we set the extra zeros
DF$y[as.logical(rbinom(n * K, size = 1, prob = plogis(eta_zi)))] <- 0
}
#create categorical time variable
DF$time_categorical[DF$time<2.5] <- "early"
DF$time_categorical[DF$time>=2.5] <- "late"
DF$time_categorical <- as.factor(DF$time_categorical)
#model with interaction in fixed effects zero part and adding nesting in zero part as in model above
km3 <- mixed_model(y ~ sex * time_categorical, random = ~ 1 | id, data = DF,
family = hurdle.lognormal(), n_phis = 1,
zi_fixed = ~ sex * time_categorical, zi_random = ~ 1 | id)
#### ATTEMPT at QDRG function in emmeans ####
coef_zero_part <- fixef(km3, sub_model = "zero_part")
vcov_zero_part <- vcov(km3)[9:12,9:12]
qd_km3 <- emmeans::qdrg(formula = ~ sex * time_categorical, data = DF,
coef = coef_zero_part, vcov = vcov_zero_part)
输出：
> joint_tests(qd_km3)
model term df1 df2 F.ratio p.value
sex 1 Inf 11.509 0.0007
time_categorical 1 Inf 0.488 0.4848
sex:time_categorical 1 Inf 1.077 0.2993
> emmeans(qd_km3, pairwise ~ sex|time_categorical)
$emmeans
time_categorical = early:
sex emmean SE df asymp.LCL asymp.UCL
male -1.592 0.201 Inf -1.99 -1.198
female -1.035 0.187 Inf -1.40 -0.669
time_categorical = late:
sex emmean SE df asymp.LCL asymp.UCL
male -1.914 0.247 Inf -2.40 -1.429
female -0.972 0.188 Inf -1.34 -0.605
Confidence level used: 0.95
$contrasts
time_categorical = early:
contrast estimate SE df z.ratio p.value
male - female -0.557 0.270 Inf -2.064 0.0390
time_categorical = late:
contrast estimate SE df z.ratio p.value
male - female -0.942 0.306 Inf -3.079 0.0021
检查对比是否与零部分固定效果相对应：
> fixef(km3, sub_model = "zero_part")
(Intercept) sexfemale time_categoricallate sexfemale:time_categoricallate
-1.5920415 0.5568072 -0.3220390 0.3849780
> (-1.5920) - (-1.5920 + 0.5568)
[1] -0.5568 #matches contrast within "early" level of "time_categorical"
> (-1.5920 + -0.3220) - (-1.5920 + -0.3220 + 0.5568 + 0.3850)
[1] -0.9418 #matches contrast within "late" level of "time_categorical"
【问题讨论】：
-
您也许可以使用
emmeans::qdrg()来创建所需的对象。请参阅其文档。您显然无法使用object参数。您需要为模型的条件部分或零部分指定数据、固定效应公式，以及相关模型部分的相关回归系数和 vcov 矩阵。对于后者，您可能必须选择系数和协方差矩阵的子集。 -
感谢您的建议！我一直在尝试，但它似乎还没有工作（可能是由于我缺乏洞察力......）。我在上面的帖子中添加了我的代码。
-
看起来基本没问题。尝试列出您用于 coef 和 vcov 的内容，并确保其中没有多余的内容。
-
它仍然不起作用，尽管
vcov元素确实对零部分没有选择性。我现在将其编辑为一个 4x4 矩阵，其中包含行/列“zi_(Intercept)”、“zi_sexfemale”、“zi_time_categoricallate”和“zi_sexfemale:time_categoricallate”。coef元素是一个包含 4 列“(Intercept)”、“sexfemale”、“time_categoricallate”和“sexfemale:time_categoricallate”的矩阵，id的每个嵌套都有一行。列的命名是否重要？还是可能有其他问题？谢谢！ -
为什么你的 coef 是矩阵？它应该是一个长度为 4 的向量。
标签： r anova mixed-models