r - 混合 glm 零膨胀模型的 Bootstrap 方法

Question

我想使用 glmmTMB 包引导混合 glm 零膨胀模型 (m_F)，但尽管使用了 coef 或 fixef 系数规范，我总是输出错误:

Error in bres[i, ] <- coef(bfit) : 
  incorrect number of subscripts on matrix

我的例子:

library(glmmTMB)
library(boot)
my.ds <- read.csv("https://raw.githubusercontent.com/Leprechault/trash/main/ds.desenvol.csv")
str(my.ds)
# 'data.frame': 400 obs. of  4 variables:
#  $ temp       : num  0 0 0 0 0 0 0 0 0 0 ...
#  $ storage    : int  5 5 5 5 5 5 5 5 5 5 ...
#  $ rep        : chr  "r1" "r2" "r3" "r4" ...
#  $ development: int  0 23 22 27 24 25 24 22 0 22 ...

# Fit a GLM mixed Hurdle (zero-inflated) log-link Gamma model
m_F <- glmmTMB(development ~ poly(temp,2) + (1 | storage), data = my.ds,
               family = ziGamma(link = "log"),
               ziformula = ~ 1)
summary(m_F)

# Create a bootstrap aproach
nboot <- 1000
bres <- matrix(NA,nrow=nboot,
                  ncol=length(coef(m_F)),
                  dimnames=list(rep=seq(nboot),
                                coef=names(coef(m_F))))
set.seed(1000)
bootsize <- 100
for (i in seq(nboot)) {
  bdat <- my.ds[sample(nrow(my.ds),size=bootsize,replace=TRUE),]
  bfit <- update(m_F, data=bdat)  ## refit with new data
  bres[i,] <- coef(bfit)
}

请问有什么帮助吗？

Answer 1

我的回答与@RuiBarradas 的有点相似，但更接近您的原始代码。要点是 coef() 并没有按照您的想法进行； (1) 约定(最初由 nlme 包设置)是对于混合模型 coef() 返回组级系数的矩阵(或矩阵列表)，而 fixef() 返回固定效应(人口水平)系数； (2) 对于 glmmTMB，fixef() 返回条件、零膨胀和分散模型的固定效应向量的列表 ( unlist() 将其折叠回具有串联名称的向量)。

要记住的另一点是，对于具有分组结构的数据集，在个体观察级别进行引导可能并不明智(您可以在组级别或组内级别或两者进行引导；您可以引导残差(如果你有一个线性模型 - 这不适用于带有计数数据的 GLMM)；你也可以使用 lme4::bootMer 进行参数引导，这几乎是唯一的选择当您的 GLMM 具有交叉随机效应时)。

PS bootsize 在这里做什么？ Bootstrap 的标准方法是对与原始大小相同的数据集重新采样并进行替换。仅对数据集的四分之一进行重采样(nrow(my.ds) == 400，bootsize == 100)定义明确但非常不寻常——您是否正在做一些有意使用特定的非标准 Bootstrap ......？

sum_fun <- function(fit) {
    unlist(fixef(fit))
}

bres <- matrix(NA,
               nrow=nboot,
               ncol=length(sum_fun(m_F)),
               dimnames=list(rep=seq(nboot),
                             coef=names(sum_fun(m_F))))
set.seed(1000)
bootsize <- 100
pb <- txtProgressBar(max = bootsize, style = 3)
for (i in seq(nboot)) {
    setTxtProgressBar(pb, i)
    bdat <- my.ds[sample(nrow(my.ds), size=bootsize,replace=TRUE),]
    bfit <- update(m_F, data=bdat)  ## refit with new data
    bres[i,] <- sum_fun(bfit)
}