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

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

• @BR whats the bottle neck here? Basically what's taking time in Rprof. Commented Feb 5, 2011 at 22:14
• One way is to just ignore the "mixed" part of the GLMM. Just fit an ordinary GLM, the estimates won't change much, but their standard errors probably will Commented Feb 5, 2011 at 22:40
• @probabilityislogic. In addition to your remark, I also think, whether answer would differ much depends on the group size and individual behavior in group. As Gelman and Hill say: mixed effects model results would be between pooling and no pooling. (Obv. this is for Bayesian hierarchical models, but mixed models are a classical way of doing the same.) Commented Feb 5, 2011 at 23:29
• @probabilityislogic: That works for LMM, but it seems to fail for GLMM (meaning that I ran models with and without the extra M on the same data and ended up with significantly different results). Unless, of course, there is an error in the implementation of glmer.
– B R
Commented Feb 5, 2011 at 23:47
• @suncoolsu: the bottleneck is the estimation of the GLMM, which may a couple of seconds (especially with several random effects). But do that 1000*1000 times, and that's 280 hours of computation. Fitting a GLM takes about 1/100 of the time.
– B R
Commented Feb 5, 2011 at 23:55

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.

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.
• Thanks for the answer. Somehow, I've overlooked the fact that I have two cores (my computer isn't very new). I should've looked at multicore a long time ago.
– B R
Commented Feb 6, 2011 at 23:16

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

• Thanks for the answer, it'll take me a bit to digest it (and I already have plans for my Saturday night). It is different enough that it is not clear to me if it gives the same answer as the GLMM approach, but I need to think about it more.
– B R
Commented Feb 6, 2011 at 0:16