How can I speed up calculation of the fixed effects in a GLMM?

Question

I'm doing a simulation study which requires bootstrapping estimates obtained from a generalized linear mixed model (actually, the product of two estimates for fixed effects, one from a GLMM and one from an LMM). To do the study well would require about 1000 simulations with 1000 or 1500 bootstrap replications each time. This takes a significant amount of time on my computer (many days).

How can I speed up the computation of these fixed effects?

To be more specific, I have subjects who are measured repeatedly in three ways, giving rise to variables X, M, and Y, where X and M are continuous and Y is binary. We have two regression equations $$M=\alpha_0+\alpha_1X+\epsilon_1$$ $$Y^*=\beta_0+\beta_1X+\beta_2M+\epsilon_2$$ where Y$^*$ is the underlying latent continuous variable for $Y$ and the errors are not iid.
The statistic we want to bootstrap is $\alpha_1\beta_2$. Thus, each bootstrap replication requires fitting an LMM and a GLMM. My R code is (using lme4)

    stat=function(dat){
        a=fixef(lmer(M~X+(X|person),data=dat))["X"]
        b=fixef(glmer(Y~X+M+(X+M|person),data=dat,family="binomial"))["M"]
        return(a*b)
    }

I realize that I get the same estimate for $\alpha_1$ if I just fit it as a linear model, so that saves some time, but the same trick doesn't work for $\beta_2$.

Do I just need to buy a faster computer? :)

Answer 1

It should help to specify starting values, though it's hard to know how much. As you're doing simulation and bootstrapping, you should know the 'true' values or the un-bootstrapped estimates or both. Try using those in the start = option of glmer.

You could also consider looking into the whether the tolerances for declaring convergence are stricter than you need. I'm not clear how to alter them from the lme4 documentation though.

Answer 2

Two other possibilities too consider, before buying a new computer.

1. Parallel computing - bootstrapping is easy to run in parallel. If your computer is reasonably new, you probably have four cores. Take a look the the multicore library in R.
2. Cloud computing is also a possibility and reasonably cheap. I have colleagues that have used the amazon cloud for running R scripts. They found that it was quite cost effective.

Answer 3

It could possibly be a faster computer. But here is one trick which may work.

Generate a simultation of $Y^*$, but only conditional on $Y$, then just do OLS or LMM on the simulated $Y^*$ values.

Supposing your link function is $g(.)$. this says how you get from the probability of $Y=1$ to the $Y^*$ value, and is most likely the logistic function $g(z)=log \Big(\frac{z}{1-z}\Big)$.

So if you assume a bernouli sampling distribution for $Y\rightarrow Y\sim Bernoulli(p)$, and then use the jeffreys prior for the probability, you get a beta posterior for $p\sim Beta(Y_{obs}+\frac{1}{2},1-Y_{obs}+\frac{1}{2})$. Simulating from this should be like lighting, and if it isn't, then you need a faster computer. Further, the samples are independent, so no need to check any "convergence" diagnostics such as in MCMC, and you probably don't need as many samples - 100 may work fine for your case. If you have binomial $Y's$, then just replace the $1$ in the above posterior with $n_i$, the number of trials of the binomial for each $Y_i$.

So you have a set of simulated values $p_{sim}$. You then apply the link function to each of these values, to get $Y_{sim}=g(p_{sim})$. Fit a LMM to $Y_{sim}$, which is probably quicker than the GLMM program. You can basically ignore the original binary values (but don't delete them!), and just work with the "simulation matrix" ($N\times S$, where $N$ is the sample size, and $S$ is the number of simulations).

So in your program, I would replace the $gmler()$ function with the $lmer()$ function, and $Y$ with a single simultation, You would then create some sort of loop which applies the $lmer()$ function to each simulation, and then takes the average as the estimate of $b$. Something like $$a=\dots$$ $$b=0$$ $$do \ s=1,\dots,S$$ $$b_{est}=lmer(Y_s\dots)$$ $$b=b+\frac{1}{s}(b_{est}-b)$$ $$end$$ $$return(a*b)$$

Let me know if I need to explain anything a bit clearer