Bayesian bootstrap in `R` with `update()`

bboot_vab.R

library(fixest)
library(data.table)
library(insight)

bboot_vab <- function(model, reps = 1e3, conf.level = .95, cluster = NULL) {
    data <- insight::get_data(model)
    setDT(data)
    if (!is.null(cluster)) {
        if (anyNA(data[[cluster]])) {
            stop("The cluster variable cannot include missing values.")
        }
        cluster_idx <- unique(data[[cluster]])
        wts <- rexp(reps * length(cluster_idx))
        wts <- matrix(wts, ncol = reps)
        wts <- data.table(cluster_idx, wts)
        colnames(wts)[1] <- cluster
        wts <- merge(data, wts, all.x = TRUE, sort = FALSE)
        wts <- as.matrix(wts[, (ncol(data)+1):ncol(wts)])    } else {
        wts <- rexp(reps * nrow(data))
        wts <- matrix(wts, ncol = reps)
    }
    fit <- function(i) {
        coef(update(model, weights = wts[, i]))
    }
    boot_mat <- t(sapply(seq_len(reps), fit))
    alpha <- 1 - conf.level
    ci <- c(alpha / 2, 1 - alpha / 2)
    out <- t(apply(boot_mat, 2, \(x) c(mean(x), quantile(x, ci))))
    return(out)
}

model <- fepois(Sepal.Length ~ Sepal.Width + Petal.Length | Species,
                data = iris)
bboot_vab(model, cluster = "Species")

Not sure whether that does the trick. I'll have a look sometime, if I find something, I'll post it here!