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

Question

我想使用m_F包引导混合glm零膨胀模型( glmmTMB )，但是尽管系数规范使用了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()将其折叠回一个带有缩名的向量)。

要记住的另一点是，对于具有分组结构的数据集(您可以在组级或组内级别进行引导，或者两者兼而有之；如果您有一个线性模型--对于带有计数数据的GLMMs，这是行不通的)；您也可以使用lme4::bootMer进行参数自举，这几乎是具有交叉随机效应的GLMMs的唯一选择)。

bootsize在这里做什么？引导的标准方法是重采样与原始数据集相同的数据集，并进行替换。仅重采样四分之一的数据集(nrow(my.ds) == 400，bootsize == 100)定义得很好，但非常不寻常--您是在故意进行某种非标准的引导吗？

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)
}