<h2>Choosing R packages for mixed effects modelling based on the car you drive</h2>
<p>There are many roads you can take to fit a mixed effects model<br>
(sometimes termed hierarchical models) in <code>R</code>. There are numerous<br>
packages that each deploy different engines to fit mixed effects models.</p>
<p>In this driveability review, I look at some of the dominant packages for<br>
mixed effects models The focus of this review is on the practicality of<br>
each package, rather than their numerical accuracy. Information on<br>
numerical accuracy can be found in the academic literature, see the<br>
bibliography below. I will lookat each package’s handling, speed and<br>
adaptability for tackling different types of questions.</p>
<p><em>Summary of the packages strengths and weaknesses. Speeds are relative<br>
to <code>lme4</code> in a test case problem, see below for details.</em></p>
<table>
<colgroup>
<col width="15%">
<col width="12%">
<col width="13%">
<col width="12%">
<col width="12%">
</colgroup>
<thead>
<tr>
<th>Parameter</th>
<th><code>lme4</code></th>
<th><code>rjags</code></th>
<th><code>INLA</code></th>
<th><code>RStan</code></th>
</tr>
</thead>
<tbody>
<tr>
<td> Relative speed </td>
<td>1</td>
<td>100</td>
<td>10</td>
<td>60</td>
</tr>
<tr>
<td> Handling- ease of coding<sup>1</sup> </td>
<td>Straightforward</td>
<td>Hard</td>
<td>Moderate</td>
<td>Hard</td>
</tr>
<tr>
<td> User manual </td>
<td>Accessible</td>
<td>Brief and technical</td>
<td>Very technical<sup>2</sup></td>
<td>Highly accessible</td>
</tr>
<tr>
<td> Online documentation </td>
<td>Extensive</td>
<td>Extensive</td>
<td>Nascent</td>
<td>Nascent</td>
</tr>
<tr>
<td> Support from developers<sup>3</sup> </td>
<td>Good</td>
<td>Good</td>
<td>Outstanding</td>
<td>Good</td>
</tr>
<tr>
<td> Testing for accuracy </td>
<td>Extensive</td>
<td>Extensive for some types of models</td>
<td>Nascent</td>
<td>Nascent</td>
</tr>
<tr>
<td> Adaptability – types of models you can fit </td>
<td>Limited</td>
<td>Your imagination is the limit</td>
<td>Extensive</td>
<td>Almost as much as JAGS</td>
</tr>
<tr>
<td> Best for </td>
<td>Simple problems, quick implementation</td>
<td>Customising complex models</td>
<td>Spatio-temporal models</td>
<td>Developing new types of models</td>
</tr>
<tr>
<td> How you will look driving it </td>
<td>Traditional, but elegant</td>
<td>You like to take the scenic route</td>
<td>Edgy and slightly mysterious</td>
<td>Hanging with the cool kids</td>
</tr>
</tbody>
</table>
<ol>
<li>but see below for some add-on packages which make life easier for<br>
STAN and JAGS</li>
<li>= lots of equations</li>
<li>In my experience. But, heck these packages are all free and all give better user support than many<br>
paid services!</li>
</ol>
<p>The first is<br>
<a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/lme4/index.html"><code>lme4</code></a>,<br>
meaning linear mixed effects models with <code>S4</code> classes. <code>lme4</code> is like an<br>
older model sports-car – fast, respectable, well known and able to<br>
handle common types of questions. <code>lme4</code> uses maximum likelihood methods<br>
for estimation, so interpret its statistics as a<br>
<a target="_blank" rel="nofollow" href="https://stats.stackexchange.com/questions/22/bayesian-and-frequentist-reasoning-in-plain-english">frequentist</a>.</p>
<p>The next three packages are Bayesian techniques.</p>
<p>The second package is<br>
<a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/rjags/index.html"><code>rjags</code></a>,<br>
which is an interface to the <a target="_blank" rel="nofollow" href="http://mcmc-jags.sourceforge.net/">“Just Another Gibbs<br>
Sampler”</a> program. <code>rjags</code> is the<br>
Land-Rover of packages for mixed effects models: well tested, able to<br>
tackle just about any terrian, but slow and prone to breaking down<br>
(usually a user error, rather than the package itself).</p>
<div>
<img alt="" src="https://i0.wp.com/www.seascapemodels.org/Images/cars-glmm.png?w=456" data-recalc-dims="1">
<p> My take on reviewing R packages for Generalised Linear Mixed Effects. Clockwise from top left: The `lme4` package is fast but a bit older, like an early noughties Ferrari; the JAGS program goes anywhere like a Range Rover, but breaks often and is slow; the Stan language is modern and efficient like a hybrid BMW and INLA is fast and modern but requires technical expertise like a Tesla. </p>
</div>
<p>The third package is <a target="_blank" rel="nofollow" href="http://www.r-inla.org/"><code>INLA</code></a>, standing for<br>
“Integrated Nested Laplace Approximation”. <code>INLA</code> is like an electric<br>
sports car: it is a modern technique that is lightning fast, but there<br>
has been limited testing to date and you better understand what is under<br>
the hood before you attempt to drive it.</p>
<p>The final package is<br>
<a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/rstan/index.html"><code>RStan</code></a>,<br>
which is the <code>R</code> interface to the <a target="_blank" rel="nofollow" href="http://mc-stan.org/">Stan modelling<br>
language</a>. <code>RStan</code> is like like BMW’s hybrid,<br>
combining computational speed with flexibility and with exceptional<br>
documentation and support.</p>
<p>Now, let’s look at each packages computational speed, handling (ease of<br>
use), driver support (documentation, existing knowledge and forums) and<br>
flexibility.</p>
<p>Skip to the end for a biblography of examples and some cool packages<br>
that build on these tools.</p>
<h3>A note about philosophy</h3>
<p>The frequentist and Bayesian approaches I review here come from<br>
<a target="_blank" rel="nofollow" href="https://stats.stackexchange.com/questions/22/bayesian-and-frequentist-reasoning-in-plain-english">fundamentally different statistical<br>
philosophies</a>.<br>
Here I will tackle the approaches pragmatically, looking at estimation<br>
speed and so on, however you might exclude one or the other a-prior<br>
based on your bent towards Bayesian or frequentist thinking.</p>
<h3>Speed</h3>
<p>I simulated a simple hierarchical data-set to test each of the models.<br>
The script is available <a target="_blank" rel="nofollow" href="http://www.seascapemodels.org/data/glmm-comparison-script.R">here</a>. The<br>
data-set has 100 binary measurements. There is one fixed covariate<br>
(continuous ) and one random effect with five groups. The linear<br>
predictor was transformed to binomial probabilities using the logit<br>
function. For the Bayesian approaches, slightly different priors were<br>
used for each package, depending on what was available. See the script<br>
for more details on priors.</p>
<p>For each package I fit the same linear model with binomial errors and a<br>
logit link. To fit the models I used the functions <code>glmer()</code> from<br>
<code>lme4</code>, <code>inla()</code> from <code>INLA</code>, custom code for <code>rjags</code> and the<br>
<code>stan_glmer()</code> from the package<br>
<a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/rstanarm/index.html"><code>rstanarm</code></a><br>
to fit the Stan model (you can also write your own custom code for<br>
<code>RStan</code>, but for convenience I used the canned version of mixed effects<br>
models).</p>
<p><em>Computational speeds for fitting a binomial Generalised Linear Model to<br>
100 samples with one random and one fixed effect. Relative times are<br>
relative to <code>lme4</code></em></p>
<table>
<thead>
<tr>
<th>Measurement</th>
<th><code>lme4</code></th>
<th><code>rjags</code></th>
<th><code>INLA</code></th>
<th><code>RStan</code></th>
</tr>
</thead>
<tbody>
<tr>
<td><strong>Speed (Seconds)</strong></td>
<td>0.063</td>
<td>6.2</td>
<td>0.59</td>
<td>3.8</td>
</tr>
<tr>
<td><strong>Speed (relative)</strong></td>
<td>1</td>
<td>99</td>
<td>9.4</td>
<td>60</td>
</tr>
</tbody>
</table>
<p>So <code>lme4</code> is by far the fastest, clocking in at around 6 tenths of a<br>
second (on my 3GHz Intel i7 Macbook Pro). <code>INLA</code> comes in second, ten<br>
times slower than <code>lme4</code>, then followed an order of magnitude later by<br>
<code>RStan</code> and finally <code>rjags</code> came plodding along behind about 100 times<br>
slower than <code>lme4</code>.</p>
<p>If you run the code I have provided you will see that all packages come<br>
up with similar estimates for the model parameters.</p>
<h3>Handling – ease of use</h3>
<p><code>lme4</code> leads the pack when it comes to ease of use. If you have your<br>
data sorted, then you can go from zero to fitted model in one line of<br>
code. <code>INLA</code> is similar, however you should generally not rely on<br>
default priors but specify your own, which requires additional code and<br>
thinking.</p>
<p>Use of <code>RStan</code> and <code>rjags</code> depends on how you use them. You can write<br>
your own model code, which can lead to quite lengthy scripts, but also<br>
gives you greater number of choices for how you design your model.<br>
However, there are also many packages that provide wrapper’s to <code>RStan</code><br>
and <code>rjags</code> that will fit common types of models (like mixed effects<br>
models). These tend to be as simple as <code>INLA</code>’s functions, remembering<br>
that you should think carefully about <a target="_blank" rel="nofollow" href="https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations">prior<br>
choice</a>.</p>
<p>If you are going to write your own models, you will need to get good at<br>
debugging, so there is an extra overhead there.</p>
<p><strong>Model checking</strong> is another important consideration.</p>
<p>For users familiar with GLMs, <code>lme4</code> probably provides the most familiar<br>
support for model checking. <code>INLA</code> can be more difficult to use, because<br>
it requires the user to hand-code many model checking routines. <code>INLA</code>,<br>
<code>rjags</code> and <code>RStan</code> also require the user to carefully check that the<br>
fitting algorithms have performed properly.</p>
<p>Each of these packages use different types of algorithms, so the<br>
statistics you will use to check algorithm performance are different.<br>
See each package’s documentation for further advice.</p>
<h3>Driver support – documentation and support</h3>
<p><code>lme4</code> and <code>RStan</code> have the highest quality user-manuals. Both are well<br>
written and provide numerous examples that users will find helpful.</p>
<p>The user manual for <code>rjags</code> and JAGS is somewhat brief, but the user can<br>
easily find many helpful examples on the web, and these packages are<br>
covered by numerous books. So if you are willing to broaden your search<br>
you should have no trouble finding help.</p>
<p><code>RStan</code> and <code>INLA</code> are both relatively new, so fewer examples can be<br>
found on the web. While <code>RStan</code> has an excellent manual for both the<br>
mathsy and non-mathsy types, documentation of <code>INLA</code> can be difficult to<br>
follow for novice users, because much of it is equations.</p>
<p><strong>Web support and chat groups</strong> are also an important aspect of model<br>
development. Help for all packages can be found on sites like<br>
<a target="_blank" rel="nofollow" href="http://stackoverflow.com/questions/tagged/lme4">Stackoverflow</a>, or you<br>
can post your own new questions there.</p>
<p>Personally, I have found that support for <code>INLA</code>, <a target="_blank" rel="nofollow" href="https://groups.google.com/forum/#!forum/r-inla-discussion-group">through their<br>
forum</a>,<br>
is exceptional. On several occaisons I have posted questions for which I<br>
could not find documentation and have generally recieved a response from<br>
<code>INLA</code>’s developers before the end of the working day. This is amazing,<br>
because it is better than you would expect from many paid services.</p>
<p>I have not participated in forums for the other packages, so I can’t add<br>
comments here.</p>
<p>One final note on support is that <code>lme4</code> and JAGS are both pretty well<br>
studied approaches to modelling and there are numerous academic papers<br>
dealing with their biases and accuracy. <code>INLA</code> and <code>RStan</code> are<br>
relatively new on the scene, so for new and different types of problems<br>
you might want to run your own simulation study to check that your model<br>
is performing as expected.</p>
<h3>Adaptability – ability to tackle different types of problems</h3>
<p>Your imagination is the limit with <code>rjags</code> – take it on well driven<br>
routes, like a mixed effects model, or off-road on new adventures with<br>
Bayesian structural equation modelling – it can do it all.</p>
<p><code>rjags</code> is followed closely by <code>RStan</code>. <code>RStan</code> has a few limitations,<br>
but can basically do anything JAGS can do, and often faster (e.g. <a target="_blank" rel="nofollow" href="http://onlinelibrary.wiley.com/doi/10.1111/2041-210X.12681/full">this<br>
ecological<br>
study</a>).<br>
<code>RStan</code> will probably eventually replace <code>rjags</code>, but to date <code>rjags</code><br>
persists because of the extensive range of well documented examples<br>
users can build on.</p>
<p><code>INLA</code> comes in fourth with a diverse range of built in likelihoods and<br>
model functions (you can even ask the developers to add more and they<br>
might do it!). <code>INLA</code> is becoming particularly popular for modelling<br>
spatial and temporal auto-correlation, partly due to its speed. You can<br>
do these types of models with <code>RStan</code> and <code>rjags</code> but computations might<br>
be impossibly slow.</p>
<p><code>lme4</code> is the least flexible of packages – though there are some options<br>
to customise it’s models. The similar <code>nlme</code> package also provides a<br>
range of tools for fitting random effects for spatial and temporal<br>
autocorrelation.</p>
<p>There are however, some clever tricks you can do with <code>lme4</code> to fit a<br>
broader range of models. For instance, you can include a random effect<br>
where every sample is a separate group in a Poisson GLM to get a Poisson<br>
with extra variation.</p>
<h3>Conclusion</h3>
<p>That’s all. I hope this blog helps you make a good choice before<br>
investing in learning a new tool for mixed effects models. <a target="_blank" rel="nofollow" href="https://twitter.com/bluecology?lang=en">Let me<br>
know</a> how you go and if you<br>
found this useful.</p>
<h3>Cool applications of mixed effects models in R</h3>
<p><a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/rstanarm/vignettes/rstanarm.html">rstanarm</a> for mixed effects models coded like <code>glm()</code> but using Stan for estimation.</p>
<p><a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/glmmBUGS/index.html">glmmBUGS</a> for mixed effects models coded like <code>glm()</code> but using JAGS for estimation.</p>
<p><a target="_blank" rel="nofollow" href="https://www.jstatsoft.org/article/view/v075i07/v75i07.pdf">A new function in the mgcv package (links to<br>
pdf)</a> for<br>
auto-writing code for Bayesian Generalised Additive Models.</p>
<p><a target="_blank" rel="nofollow" href="http://www.r-inla.org/barrier-model">A new study</a> showing how you can<br>
fit spatial models with barriers.</p>
<p>Fit <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/boral/index.html">Bayesian Latent Variable<br>
models</a> for<br>
analysing multivariate data (e.g. ecological communities).</p>
<p><a target="_blank" rel="nofollow" href="http://glmm.wikidot.com/faq">The GLMM in<br>
R FAQ</a></p>
<h3>Bibliography</h3>
<p>Unless otherwise indicated, these resources are open access.</p>
<p><a target="_blank" rel="nofollow" href="http://mc-stan.org/documentation/">Stan</a> user manual.</p>
<p><a target="_blank" rel="nofollow" href="https://www.google.com.au/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwjn2ZSp9KPTAhVDTrwKHRueAbMQFggjMAA&url=http%3A%2F%2Fwww.stats.ox.ac.uk%2F~nicholls%2FMScMCMC15%2Fjags_user_manual.pdf&usg=AFQjCNH7r4e1fZFtwIh9qlUGLcuHXXXZdQ">JAGS</a><br>
user manual.</p>
<p><a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf">Introduction to lme4<br>
(pdf)</a>.</p>
<p><a target="_blank" rel="nofollow" href="http://www.r-inla.org/home">INLA</a> project page.</p>
<p><a target="_blank" rel="nofollow" href="https://groups.google.com/forum/#!forum/r-inla-discussion-group">INLA</a><br>
forum.</p>
<p><a target="_blank" rel="nofollow" href="http://www.springer.com/gp/book/9780387874579">Mixed Effects Models and Extensions in Ecology with<br>
R</a> my go-to book for<br>
theory and code for fitting likelihood based mixed effects models (e.g.<br>
with <code>lme4</code>). You will have to buy it.</p>
<p><a target="_blank" rel="nofollow" href="http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0075957">A good example for how to design a simulation study to test your<br>
models</a></p>
<p><a target="_blank" rel="nofollow" href="http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2012.00211.x/full">Example for INLA for spatial<br>
models</a><br>
in ecology.</p>
<p><a target="_blank" rel="nofollow" href="https://academic.oup.com/icesjms/article/72/5/1245/758333/Mixed-effects-a-unifying-framework-for-statistical">A review of mixed effects models in fisheries<br>
science</a><br>
(good for other disciplines too).</p>
<p><a target="_blank" rel="nofollow" href="https://r-nimble.org/">NIMBLE</a> another R package for Bayesian<br>
modelling.</p>
<p><a target="_blank" rel="nofollow" href="http://glmmadmb.r-forge.r-project.org/">glmmADMB</a> another R<br>
package giving greater flexibility for fitting using maximum likelihood.</p>
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/zaPZO1_vG98">