glmmTMB零膨胀公式指定AR-1协方差结构参数估计异常咨询

问题描述 #

在glmmTMB的条件公式中指定AR-1协方差结构时，能按预期估计2个参数（sigma和rho），但在零膨胀公式中指定该结构却会估计大量参数。以下是复现示例代码：

n <- 25                                              ## Number of time points
x <- MASS::mvrnorm(mu = rep(0,n),
                   Sigma = .7 ^ as.matrix(dist(1:n)) )    ## Simulate the process using the MASS package
y <- x + rnorm(n)                                   ## Add measurement noise

times <- factor(1:n, levels=1:n)
head(levels(times))

group <- factor(rep(1,n))

dat0 <- data.frame(y, times, group)

library(glmmTMB)
model <- glmmTMB(y ~ ar1(factor(times) + 0 | group), ziformula = ~ ar1(factor(times) + 0 | group), data=dat0, family = gaussian())
summary(model)

原因解析 #

核心问题是glmmTMB的零膨胀公式（ziformula）不支持指定AR-1这类特殊协方差结构，具体细节如下：

• 条件模型中的ar1(factor(times)+0|group)是合规用法：这里ar1()用于定义随机效应的协方差结构，表示同一group内不同time对应的随机效应服从AR-1相关结构，因此仅需估计2个参数（随机效应的标准差sigma和自相关系数rho）。
• 零膨胀公式的语法限制：glmmTMB目前不允许在ziformula的随机效应部分使用ar1()这类协方差结构语法。当你在ziformula中写入ar1(...)时，程序不会将其识别为协方差结构定义，而是错误地将ar1()当作生成预测变量的函数，对每个time水平生成一个独立参数，最终参数数量等于time的水平数（示例中为25个），导致参数远超预期。

可行替代方案 #

如果需要为零膨胀部分加入时间相关的结构，可以尝试以下方式：

• 将time作为连续变量，拟合自回归项（例如加入滞后1期的time变量作为固定效应）；
• 在ziformula中使用普通随机效应（如~(1|times)），但这种方式无法实现AR-1的相关结构；
• 更换支持零膨胀模型+特殊协方差结构的统计包（如brms）。

内容的提问来源于stack exchange，提问作者StatisticsFanBoy