glmmTMB 上的计划对比

Question

如果这是一个重复的问题，我们深表歉意。许多人 post 正在寻找一种方法来对 glmmTMB 中的条件模型（固定因子）进行 post-hoc 分析。我想在某些组之间进行计划对比，而不是测试每个成对比较（例如 Tukey）。

下面的代码在 nlme:lme 上运行良好。但是，它 returns 以下代码的错误。

Error in modelparm.default(model, ...) : 
  dimensions of coefficients and covariance matrix don't match

有没有办法在 glmmTMB 上进行计划对比？

#filtdens is a dataframe and TRT,DATE,BURN,VEG are factors
filtdens <- merged %>% filter(!BLOCK %in% c("JB2","JB4","JB5") & MEAS =="DENS" & 
                      group == "TOT" & BURN == "N" & VEG == "C")
filtdens$TD <- interaction(filtdens$TRT, filtdens$DATE)
mod2 <- glmmTMB(count~(TD)+(1|BLOCK),
                 data=filtdens,
        zi=~1,
        family=nbinom1(link = "log"))

k1 <- matrix(c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,

       0, 0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0, 0,

       0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0,

       0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1), byrow = T, ncol = 12)

summary(glht(mod2, linfct=k1),test=adjusted("bonferroni"))

Answer 1

一个可重现的例子会有所帮助，但是：this vignette 在开发版本中提供了应该启用 multcomp::linfct 的代码，即：

glht_glmmTMB <- function (model, ..., component="cond") {
    glht(model, ...,
         coef. = function(x) fixef(x)[[component]],
         vcov. = function(x) vcov(x)[[component]],
         df = NULL)
}
modelparm.glmmTMB <- function (model, 
                               coef. = function(x) fixef(x)[[component]],
                               vcov. = function(x) vcov(x)[[component]],
                               df = NULL, component="cond", ...) {
    multcomp:::modelparm.default(model, coef. = coef., vcov. = vcov.,
                        df = df, ...)
}

测试（此示例与 Tukey 一起使用，但我不明白为什么它不能更普遍地工作...）

library(glmmTMB)
data("cbpp",package="lme4")
cbpp_b1 <- glmmTMB(incidence/size~period+(1|herd),
               weights=size,family=binomial,
               data=cbpp)
g1 <- glht(cbpp_b1, linfct = mcp(period = "Tukey"))
summary(g1)

这适用于当前的 CRAN 版本，但是 glmmTMB 的当前 development 版本提供了更多选项（例如 emmeans()；请参阅上面的链接小插图）。您需要通过 devtools::install_github("glmmTMB/glmmTMB/glmmTMB") 安装（您还需要安装编译工具）。

glmmTMB 上的计划对比

Planned Contrasts on glmmTMB

comparison

r

mixed-models