关于glmmTMB中simulate_new()参数设置与功效分析的技术问询

我打算基于预实验数据开展功效分析，为后续实验设计提供依据。目前用glmmTMB拟合了包含固定预测变量和随机效应的负二项模型，之前用simr对glmer模型做功效分析，但simr不支持glmmTMB，所以计划用glmmTMB内置函数实现：

• 用预实验数据拟合模型
• 用估计系数通过simulate_new()模拟新数据
• 对模拟数据拟合模型并评估显著性
• 重复步骤2-3计算功效

我想明确simulate_new()的newparams参数各值含义，发现帮助文档用的是任意值而非拟合得到的系数，想确认我用模型拟合结果设置newparams的代码是否正确。

帮助文档示例代码 #

sim_count <- simulate_new(~ mined + (1|site),
             newdata = Salamanders,
             zi = ~ mined,
             family = nbinom2,
             newparams = list(beta = c(2, 1),
                         betazi = c(-0.5, 0.5), ## logit-linear model for zi
                         betadisp = log(2), ## log(NB dispersion)
                         theta = log(1)) ## log(among-site SD)
)

我的模型拟合结果 #

Family: nbinom2  ( log )
Formula:          count ~ mined + (1 | site)
Zero inflation:         ~mined
Data: Salamanders

     AIC      BIC   logLik deviance df.resid 
  1741.4   1768.2   -864.7   1729.4      638 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 site   (Intercept) 0.1338   0.3658  
Number of obs: 644, groups:  site, 23

Dispersion parameter for nbinom2 family (): 0.965 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.5573     0.3711  -1.501    0.133    
minedno       1.4674     0.3616   4.059 4.94e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   0.2404     0.4863   0.494  0.62103   
minedno      -2.2561     0.7829  -2.882  0.00396 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我的代码 #

sim_count <- simulate_new(~ mined + (1|site),
             newdata = Salamanders,
             zi = ~ mined,
             family = nbinom2,
             newparams = list(beta = c(-0.56, 1.47),
                         betazi = c(0.24, -2.26), ## logit-linear model for zi
                         betadisp = log(0.97), ## log(NB dispersion)
                         theta = log(0.37)) ## log(among-site SD)

回答 #

你的newparams参数设置逻辑完全正确，各参数与模型输出的对应关系准确：

• beta：对应条件模型的系数，顺序为截距项→minedno的系数，你取的-0.56（对应原模型-0.5573）、1.47（对应原模型1.4674）符合要求
• betazi：对应零膨胀模型的系数，顺序为截距项→minedno的系数，0.24（原0.2404）、-2.26（原-2.2561）设置正确
• betadisp：负二项分布离散参数的对数，原模型离散参数为0.965，log(0.97)的转换符合simulate_new()的参数要求
• theta：随机效应标准差的对数，原模型site组随机截距标准差为0.3658，log(0.37)的设置没问题

补充两个细节：

1. 你的代码末尾缺少闭合的括号和大括号，完整代码应补上：

sim_count <- simulate_new(~ mined + (1|site),
             newdata = Salamanders,
             zi = ~ mined,
             family = nbinom2,
             newparams = list(beta = c(-0.56, 1.47),
                         betazi = c(0.24, -2.26), ## logit-linear model for zi
                         betadisp = log(0.97), ## log(NB dispersion)
                         theta = log(0.37)) ## log(among-site SD)
             )

2. 如果追求精度，建议直接从模型对象中提取参数而非手动四舍五入，避免误差且更自动化：

# 假设拟合好的模型对象名为mod
mod <- glmmTMB(count ~ mined + (1|site), zi=~mined, family=nbinom2, data=Salamanders)
newparams <- list(
  beta = fixef(mod)$cond,
  betazi = fixef(mod)$zi,
  betadisp = log(mod$fit$par["betadisp"]),
  theta = log(sigma(mod)$cond)
)
# 传入simulate_new
sim_count <- simulate_new(~ mined + (1|site),
             newdata = Salamanders,
             zi = ~ mined,
             family = nbinom2,
             newparams = newparams
)

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