bootmer 和 glmer 的限制出错

Question

当我尝试引导 glmer 对象时，我遇到了一个错误，我找不到任何文档：

glm2 <- glmer(RT~valence+location+first_location+Trial_num + 
                      (1+Trial_num|id)+(1|Trial_num), 
                      family=inverse.gaussian(log), 
                      control = glmerControl(optimizer = "nloptwrap",
                                      calc.derivs = FALSE), data=df_long)

错误是：

lme4::.simulateFun(object = , : 找不到函数“sfun”

这与我尝试bootMer 或confint 无关：

bootMer_out <- bootMer(glm2,FUN=fixef, nsim=300)

confint_out <- confint(glm2, method="boot")

当我作为 lmer 对象运行时，我没有引导问题。即

lm2 <- glmer(RT~valence+location+first_location+Trial_num + (1+Trial_num|id)+(1|Trial_num), family=inverse.gaussian(log), control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data=df_long))  

是否与链接功能有关？有解决方法吗？我在simulateFun 文档中也找不到函数'sfun'。我总是可以单独对数据进行转换并使用 lmer 而不是 glmer，但是如果有人有一些见解会很棒（因为我现在很好奇）。

Answer 1

正如@user20650 所指出的，您需要为逆高斯族添加一个模拟方法。

例如，我将这些添加到我的 lme4 fork 的 predict.R 下的分支中：

rinvgauss <- function(n, mu, lambda) {
    # transcribed from https://en.wikipedia.org/wiki/Inverse_Gaussian_distribution
    nu <- rnorm(n)
    y <- nu^2
    x <- mu + (mu^2 * y)/(2*lambda) - (mu/(2*lambda)) * sqrt(4*mu*lambda*y + mu^2*y^2)
    z <- runif(n)
    ifelse(z <= mu/(mu + x), x, mu^2/x)
}

inverse.gaussian_simfun <- function(object, nsim, ftd = fitted(object),
                                    wts = weights(object)) {
    if (any(wts != 1)) message("using weights as inverse variances")
    dispersion <- sum((weights(object, 'working') * 
        resid(object, 'working')^2)[weights(object, 'working')>0])/df.residual(object)
    rinvgauss(nsim * length(ftd), mu = ftd, 
                     lambda = wts/dispersion)
}
# ... skip a few
simfunList <- list(gaussian = gaussian_simfun,
       binomial = binomial_simfun,
       poisson  = poisson_simfun,
       Gamma    = Gamma_simfun,
       negative.binomial = negative.binomial_simfun,
       inverse.gaussian = inverse.gaussian_simfun)

这是一个例子：

# devtools::install_github('aforren1/lme4', ref = 'add_invgauss_simulate')
library(lme4)
set.seed(1)
dat <- data.frame(y = lme4:::rinvgauss(1000, 3, 4),
                  x = runif(1000), 
                  subj = factor(rep(1:10, 100)))

mod <- glmer(y ~ x + (1|subj), 
             data = dat, 
             family = inverse.gaussian(link='log'))

# ~60 secs on my laptop
(boots <- confint(mod, method = 'boot', nsim = 100, parm = 'beta_')) 
                 2.5 %   97.5 %
(Intercept)  1.0044813 1.248774
x           -0.2158155 0.161213

(walds <- confint(mod, method = 'Wald', parm = 'beta_'))
                2.5 %    97.5 %
(Intercept)  1.000688 1.2289971
x           -0.205546 0.1644621

您可以看到 bootstrap 方法给出的结果（大致）与 Wald 方法相同。