Error with bootmer and confint for glmer

I'm running into an error that I can't find any documentation on when I try to bootstrap a glmer object:

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)

The error is:

Error in lme4::.simulateFun(object = , : could not find function "sfun

This is regardless of whether I try bootMer or confint:

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

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

When I run as an lmer object I don't have the issue with bootstrapping. i.e.

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

Does it have to do with the link function? Is there a workaround? I couldn't find function 'sfun' in the simulateFun documentation either. I could always just do the transformation on the data separately and use lmer instead of glmer, but if anyone has some insight that would be great (since I'm curious now).

Solution

As pointed out by @user20650, you'll need to add a simulation method for the inverse gaussian family.

For example, I added these to a branch on my lme4 fork under 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)

Here's an example:

# 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

You can see that the bootstrap method gives (roughly) the same results as the Wald method.