<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/fitting-count-and-zero-inflated-count-glmms-with-mgcv/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.fromthebottomoftheheap.net/2017/05/04/compare-mgcv-with-glmmTMB/"> From the Bottom of the Heap - R</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>
A couple of days ago, Molly Brooks and coauthors posted a <a target="_blank" rel="nofollow" href="http://doi.org/10.1101/132753">preprint</a> on <a target="_blank" rel="nofollow" href="http://biorxiv.org/">BioRχiv</a> illustrating the use of the <strong>glmmTMB</strong> R package for fitting zero-inflated GLMMs <span data-cites="Brooks2017-so">(Brooks et al., 2017)</span>. In the paper, <strong>glmmTMB</strong> is compared with several other GLMM-fitting packages. <strong>mgcv</strong> has recently gained the ability to fit a wider range of families beyond the exponential family of distributions, including zero-inflated Poisson models. <strong>mgcv</strong> can also fit simple GLMMs through a spline equivalent of a Gaussian random effect. So, whilst I was waiting on some Bayesian GAMs to finish sampling, I decided to see how <strong>mgcv</strong> compared against <strong>glmmTMB</strong> on the two examples used in the paper.
</p>
<p>
For this post I’ll be using a couple of packages beyond <strong>glmmTMB</strong> and <strong>mgcv</strong>; make sure you have <strong>ggplot2</strong> and <strong>ggstance</strong> installed if you wish to run through the code below.
</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>&quot;glmmTMB&quot;</span><span>)</span>
<span>library</span><span>(</span><span>&quot;mgcv&quot;</span><span>)</span>
<span>library</span><span>(</span><span>&quot;ggplot2&quot;</span><span>)</span>
<span>theme_set</span><span>(</span><span>theme_bw</span><span>())</span>
<span>library</span><span>(</span><span>&quot;ggstance&quot;</span><span>)</span></code></pre>
</figure>
<p>
There are several ways in which <strong>mgcv</strong> allows GLMMs to be fitted, but the way that interests me here is via <code>gam()</code> and the <em>random effect</em> spline basis. Penalised splines of the type provided in <strong>mgcv</strong> can also be represented in mixed model form, such that GAMs can also be fitted using mixed effect modelling software. The general idea is that the spline is decomposed into two parts:
</p>
<ol type="1">
<li>
the perfectly smooth parts of the basis, namely those functions, including constant and linear functions, in the penalty null space of the spline. These are added to the fixed effects model matrix, whilst,
</li>
<li>
the remaining wiggly parts of the basis are treated as random effects.
</li>
</ol>
<p>
Given this duality between splines and random effects, you can reverse the idea and create a spline basis that is the equivalent of a simple Gaussian i.i.d random effect, such that you can fit a GLMM or GAMM using GAM software like <strong>mgcv</strong>. <strong>mgcv</strong> has the <code>re</code> basis for this, and I’ll exploit that to fit the zero-inflated GLMMs to the two examples.
</p>
<p>
In <span data-cites="Brooks2017-so">Brooks et al. (2017)</span>, two example data sets are used;
</p>
<ol type="1">
<li>
<code>Salamanders</code> — Seven combinations of different salamander species and lif-stagess were sampled repeatedly sampled four times at 23 sites in Applachian streams <span data-cites="Price2016-no">(Price et al., 2016)</span>. Some of the streams were impacted by mountaintop removal and valley filling from coal mining. The data are available from <span data-cites="Price2015-se">Price et al. (2015)</span>, as well as the <strong>glmmTMB</strong> package.
</li>
<li>
<code>Owls</code> — the second example is a well-studied one in mixed modelling papers and textbooks <span data-cites="Zuur2009-vg">(Zuur et al., 2009, <span data-cites="Bolker2013-vl">Bolker et al. (2013)</span>)</span>, and relates to the begging behaviour of owl nestlings. The data were originally reported in <span data-cites="Roulin2007-rq">Roulin and Bersier (2007)</span>.
</li>
</ol>
<h3>
Salamanders<br>
</h3>
<p>
<span data-cites="Brooks2017-so">Brooks et al. (2017)</span> fit several count models to the <code>Salamander</code> data set, including standard Poisson GLMMs, negative binomial GLMMs, with <span>()</span> estimated and modelled via a linear predictor, as well as zero-inflated Poisson (ZIP) and zero-inflated negative binomial (ZINB) models. Of these, <code>gam()</code> can currently fit all but the negative binomial with <span>()</span> modelled via a linear predictor and the ZINB models.
</p>
<p>
The best fitting model of those presented was a negative binomial model, whilst <span data-cites="Brooks2017-so">Brooks et al. (2017)</span> also illustrate how to generate fitted values from the ZIP. Rather than go through fitting all of the <span data-cites="Brooks2017-so">Brooks et al. (2017)</span> models, I restrict fitting here to these two models. A <a target="_blank" rel="nofollow" href="https://gist.github.com/gavinsimpson/8a0f0e072b095295cf5f7af2762e05a7">gist</a> with code to fit all the models that <code>gam()</code> is capable of is available on Github. I have named the models similarly to <span data-cites="Brooks2017-so">Brooks et al. (2017)</span> to facilitate comparison.
</p>
<figure>
<pre><code data-lang="r"><span>nbgam2</span> <span>&lt;-</span> <span>gam</span><span>(</span><span>count</span> <span>~</span> <span>spp</span> <span>*</span> <span>mined</span> <span>+</span> <span>s</span><span>(</span><span>site</span><span>,</span> <span>bs</span> <span>=</span> <span>&quot;re&quot;</span><span>),</span> <span>data</span> <span>=</span> <span>Salamanders</span><span>,</span>
              <span>family</span> <span>=</span> <span>nb</span><span>,</span> <span>method</span> <span>=</span> <span>&quot;ML&quot;</span><span>)</span>
<span>nbm2</span> <span>&lt;-</span> <span>glmmTMB</span><span>(</span><span>count</span> <span>~</span> <span>spp</span> <span>*</span> <span>mined</span> <span>+</span> <span>(</span><span>1</span> <span>|</span> <span>site</span><span>),</span> <span>data</span> <span>=</span> <span>Salamanders</span><span>,</span>
                <span>family</span> <span>=</span> <span>nbinom2</span><span>)</span></code></pre>
</figure>
<p>
As <code>glmmTMB()</code> is currently only capable of fitting models using maximum likelihood, not REML, I use the Laplace approximate maximum likelihood estimation method for <code>gam()</code>. The new <code>nb</code> family in <strong>mgcv</strong> is for the negative binomial distribution with the (fixed) dispersion parameter <span>()</span> estimated as a model parameter, in the same way that <code>MASS::glm.nb()</code> and <code>lme4::glmer.nb()</code> models do.
</p>
<p>
In the <code>gam()</code> model, the random effect is specified using the standard <code>s()</code> smooth function with the <code>“re”</code> basis selected. The named variable, here <code>site</code>, should be stored as a factor in the data object to avoid problems.
</p>
<p>
The figure below compares the coefficient estimates returned by <code>glmmTMB()</code> and <code>gam()</code>; they are very similar, which is encouraging.
</p>
<figure>
<pre><code data-lang="r"><span>nb2.coefs</span> <span>&lt;-</span> <span>data.frame</span><span>(</span><span>estimate</span> <span>=</span> <span>c</span><span>(</span><span>coef</span><span>(</span><span>summary</span><span>(</span><span>nbm2</span><span>))</span><span>$</span><span>cond</span><span>[,</span> <span>&quot;Estimate&quot;</span><span>],</span> <span>coef</span><span>(</span><span>nbgam2</span><span>)[</span><span>c</span><span>(</span><span>1</span><span>:</span><span>14</span><span>)]),</span>
                        <span>model</span>    <span>=</span> <span>rep</span><span>(</span><span>c</span><span>(</span><span>&quot;glmmTMB&quot;</span><span>,</span> <span>&quot;mgcv::gam&quot;</span><span>),</span> <span>each</span> <span>=</span> <span>14</span><span>),</span>
                        <span>term</span>     <span>=</span> <span>rep</span><span>(</span><span>names</span><span>(</span><span>coef</span><span>(</span><span>nbgam2</span><span>)[</span><span>c</span><span>(</span><span>1</span><span>:</span><span>14</span><span>)]),</span> <span>2</span><span>))</span>

<span>ggplot</span><span>(</span><span>nb2.coefs</span><span>,</span> <span>aes</span><span>(</span><span>x</span> <span>=</span> <span>estimate</span><span>,</span> <span>y</span> <span>=</span> <span>term</span><span>,</span> <span>colour</span> <span>=</span> <span>model</span><span>,</span> <span>shape</span> <span>=</span> <span>model</span><span>))</span> <span>+</span>
    <span>geom_point</span><span>(</span><span>position</span> <span>=</span> <span>position_dodgev</span><span>(</span><span>height</span> <span>=</span> <span>0.3</span><span>))</span> <span>+</span>
    <span>labs</span><span>(</span><span>y</span> <span>=</span> <span>NULL</span><span>,</span>
         <span>x</span> <span>=</span> <span>&quot;Regression estimate&quot;</span><span>,</span>
         <span>title</span> <span>=</span> <span>&quot;Comparing mgcv with glmmTMB&quot;</span><span>,</span>
         <span>subtitle</span> <span>=</span> <span>&quot;Salamander: Negative Binomial&quot;</span><span>)</span></code></pre>
</figure>
<figure>
<img src="https://i2.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/compare-mgcv-with-glmmtmb-salamander-nb2-coefs-1.png?w=456" data-recalc-dims="1"><br>
</figure>
<p>
The values (posterior modes, or means) for the <code>site</code> random effect can also be compared
</p>
<figure>
<pre><code data-lang="r"><span>nbgam2.r</span> <span>&lt;-</span> <span>coef</span><span>(</span><span>nbgam2</span><span>)[</span><span>-</span><span>c</span><span>(</span><span>1</span><span>:</span><span>14</span><span>)]</span>
<span>nbm2.r</span>   <span>&lt;-</span> <span>ranef</span><span>(</span><span>nbm2</span><span>)</span><span>$</span><span>cond</span><span>$</span><span>site</span><span>[,</span><span>1</span><span>]</span>
<span>nms</span> <span>&lt;-</span> <span>sub</span><span>(</span><span>&quot;s\\(site\\)\\.&quot;</span><span>,</span> <span>&quot;Site &quot;</span><span>,</span> <span>names</span><span>(</span><span>nbgam2.r</span><span>))</span>
<span>ranefs</span> <span>&lt;-</span> <span>data.frame</span><span>(</span><span>ranef</span> <span>=</span> <span>c</span><span>(</span><span>unname</span><span>(</span><span>nbgam2.r</span><span>),</span> <span>nbm2.r</span><span>),</span>
                     <span>model</span> <span>=</span> <span>rep</span><span>(</span><span>c</span><span>(</span><span>&quot;glmmTMB&quot;</span><span>,</span> <span>&quot;mgcv::gam&quot;</span><span>),</span> <span>each</span> <span>=</span> <span>length</span><span>(</span><span>nbgam2.r</span><span>)),</span>
                     <span>site</span>  <span>=</span> <span>rep</span><span>(</span><span>nms</span><span>,</span> <span>2</span><span>))</span>
<span>ranefs</span> <span>&lt;-</span> <span>transform</span><span>(</span><span>ranefs</span><span>,</span> <span>site</span> <span>=</span> <span>factor</span><span>(</span><span>site</span><span>,</span> <span>nms</span><span>[</span><span>order</span><span>(</span><span>nbgam2.r</span><span>)]))</span>

<span>ggplot</span><span>(</span><span>ranefs</span><span>,</span> <span>aes</span><span>(</span><span>x</span> <span>=</span> <span>ranef</span><span>,</span> <span>y</span> <span>=</span> <span>site</span><span>,</span> <span>colour</span> <span>=</span> <span>model</span><span>,</span> <span>shape</span> <span>=</span> <span>model</span><span>))</span> <span>+</span>
    <span>geom_point</span><span>(</span><span>position</span> <span>=</span> <span>position_dodgev</span><span>(</span><span>height</span> <span>=</span> <span>0.5</span><span>))</span> <span>+</span>
    <span>labs</span><span>(</span><span>y</span> <span>=</span> <span>NULL</span><span>,</span>
         <span>x</span> <span>=</span> <span>&quot;Random effect&quot;</span><span>,</span>
         <span>title</span> <span>=</span> <span>&quot;Comparing mgcv with glmmTMB&quot;</span><span>,</span>
         <span>subtitle</span> <span>=</span> <span>&quot;Salamanders: Negative Binomial&quot;</span><span>)</span></code></pre>
</figure>
<figure>
<img src="https://i0.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/compare-mgcv-with-glmmtmb-slamander-nb2-ranefs-1.png?w=456" data-recalc-dims="1"><br>
</figure>
<p>
As the figure above shows, these too are essentially equivalent for the two fits.
</p>
<p>
The <code>summary()</code> output for the <code>glmmTMB()</code> model conveniently provides some additional useful information, in the context of GLMMs most notably the estimated variances (or standard deviations) of the random effect terms. As <code>gam()</code> wasn’t designed with GLMMs specifically in mind, the same information is not provided in the the <code>summary()</code> method for <code>gam()</code> model fits. However, Simon Wood has provided the <code>gam.vcomp()</code> function, which can be used to return the variance components of the model in a way that allows comparison with other mixed-models specific software.
</p>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>nbm2</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text"> Family: nbinom2  ( log )
Formula:          count ~ spp * mined + (1 | site)
Data: Salamanders

     AIC      BIC   logLik deviance df.resid 
  1663.4   1734.8   -815.7   1631.4      628 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 site   (Intercept) 0.2842   0.5331  
Number of obs: 644, groups:  site, 23

Overdispersion parameter for nbinom2 family ():    1 

Conditional model:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -3.3750     0.7576  -4.455 8.40e-06 ***
sppPR              0.9306     0.8773   1.061 0.288829    
sppDM              2.2485     0.7878   2.854 0.004314 ** 
sppEC-A            0.7143     0.9052   0.789 0.430029    
sppEC-L            1.8130     0.8130   2.230 0.025741 *  
sppDES-L           2.5111     0.7795   3.221 0.001275 ** 
sppDF              2.5765     0.7801   3.303 0.000957 ***
minedno            4.1619     0.7932   5.247 1.55e-07 ***
sppPR:minedno     -2.5831     0.9328  -2.769 0.005617 ** 
sppDM:minedno     -2.1495     0.8258  -2.603 0.009245 ** 
sppEC-A:minedno   -1.5828     0.9461  -1.673 0.094339 .  
sppEC-L:minedno   -1.3383     0.8493  -1.576 0.115100    
sppDES-L:minedno  -1.9358     0.8164  -2.371 0.017729 *  
sppDF:minedno     -2.7426     0.8217  -3.338 0.000844 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
</figure>
<p>
Now the <code>gam()</code> version, conveniently with a confidence interval
</p>
<figure>
<pre><code data-lang="r"><span>gam.vcomp</span><span>(</span><span>nbgam2</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Standard deviations and 0.95 confidence intervals:

          std.dev    lower     upper
s(site) 0.5325309 0.327768 0.8652132

Rank: 1/1</code></pre>
</figure>
<p>
One further analysis that <span data-cites="Brooks2017-so">Brooks et al. (2017)</span> do with the <code>Salamander</code> data (in their Appendix B) is to demonstrate how to generate and plot fitted values from the model. To do this, the analyst needs to consider whether to and how to marginalise over or condition on the random effects. The Appendix has some details on this more generally (via a linked reference) and more-specific pointers on how to go about doing this with <code>glmmTMB()</code> models. In the next few code chunks I will illustrate how to achieve the result from their section <em>Alternative prediction method</em>, where the aim is to predict at the population mode by setting the random effect component to 0. To illustrate this, <span data-cites="Brooks2017-so">Brooks et al. (2017)</span> use the more complex ZIP model with linear predictors for both the mean and the zero-inflation components of the model. I fit those models first
</p>
<figure>
<pre><code data-lang="r"><span>## glmmTMB()
</span><span>zipm3</span> <span>&lt;-</span> <span>glmmTMB</span><span>(</span><span>count</span> <span>~</span> <span>spp</span> <span>*</span> <span>mined</span> <span>+</span> <span>(</span><span>1</span> <span>|</span> <span>site</span><span>),</span> <span>zi</span> <span>=</span> <span>~</span> <span>spp</span> <span>*</span> <span>mined</span><span>,</span>
                 <span>data</span> <span>=</span> <span>Salamanders</span><span>,</span> <span>family</span> <span>=</span> <span>poisson</span><span>)</span>

<span>## gam()
</span><span>zipgam3</span> <span>&lt;-</span> <span>gam</span><span>(</span><span>list</span><span>(</span><span>count</span> <span>~</span> <span>spp</span> <span>*</span> <span>mined</span> <span>+</span> <span>s</span><span>(</span><span>site</span><span>,</span> <span>bs</span> <span>=</span> <span>&quot;re&quot;</span><span>),</span>
                    <span>~</span> <span>spp</span> <span>*</span> <span>mined</span><span>),</span>
               <span>data</span> <span>=</span> <span>Salamanders</span><span>,</span> <span>family</span> <span>=</span> <span>ziplss</span><span>,</span> <span>method</span> <span>=</span> <span>&quot;REML&quot;</span><span>)</span></code></pre>
</figure>
<p>
The <code>glmmTMB()</code> model has the zero-inflation linear predictor specified via the <code>ziformula</code> argument (abbreviated to <code>zi</code> above). With <code>gam()</code> however, multiple linear predictors are specified via a list of formula objects, only the first of which has a response (left-hand-side). The first formula, with the response, is for the Poisson mean, whilst the second is for the zero-inflation component. Note also that we use the special <code>ziplss()</code> family and that now the model is being estimated using REML, because that is the only option available for these models, which Simon Wood calls <strong>general smooth models</strong> <span data-cites="Wood2016-fx">(Wood et al., 2016)</span>. Do note that there is (as of writing) no <code>link</code> argument(s) for the <code>ziplss()</code> family. This is due to the way the model is parameterised internally in the software. This will require us to pay particular attention to the implementation shortly.
</p>
<p>
To recreate part of Figure B.3 in Appendix B <span data-cites="Brooks2017-so">(Brooks et al., 2017)</span>, the code below predicts from the fitted <code>gam()</code> model for all combinations of the factors <code>mined</code> and <code>spp</code>. Notice how we have to specify a <code>site</code> in the prediction data, otherwise <code>predict()</code> will throw a tantrum. To set the random effect for <code>site</code> to zero, use the <code>exclude</code> argument. To exclude (i.e. set to zero) any model term, you supply a character vector or list of terms to <code>exclude</code>. For smooth terms, these must be named as they appear in <code>summary(model)</code>, hence the use of <code>“s(site)”</code>. The final step is to call <code>predict()</code> with <code>type = “link”</code>. This will return a two column matrix (or a list of two-column matrices if <code>se.fit = TRUE</code> is also used).
</p>
<figure>
<pre><code data-lang="r"><span>## Newdata
</span><span>newd0</span> <span>&lt;-</span> <span>newd</span> <span>&lt;-</span> <span>as.data.frame</span><span>(</span><span>cbind</span><span>(</span><span>unique</span><span>(</span><span>Salamanders</span><span>[,</span> <span>c</span><span>(</span><span>&quot;mined&quot;</span><span>,</span><span>&quot;spp&quot;</span><span>)]),</span> <span>site</span> <span>=</span> <span>&quot;R -1&quot;</span><span>))</span>
<span>rownames</span><span>(</span><span>newd0</span><span>)</span> <span>&lt;-</span> <span>rownames</span><span>(</span><span>newd</span><span>)</span> <span>&lt;-</span> <span>NULL</span>
<span>pred</span> <span>&lt;-</span> <span>predict</span><span>(</span><span>zipgam3</span><span>,</span> <span>newd</span><span>,</span> <span>exclude</span> <span>=</span> <span>&quot;s(site)&quot;</span><span>,</span> <span>type</span> <span>=</span> <span>&quot;link&quot;</span><span>)</span>
<span>head</span><span>(</span><span>pred</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">        [,1]       [,2]
1 0.36061171 -3.7727169
2 0.94857203  0.3875087
3 0.07601834 -2.6504763
4 0.35637762 -1.3459854
5 0.55867674 -1.4747253
6 1.14836206  0.3266343</code></pre>
</figure>
<p>
The first column is the predicted value of the response from the Poisson part of the model <em>on the scale of the linear predictor</em> (the log scale). The second column is the predicted value from the zero-inflation component and is on the complementary log-log scale. Both of these need to be back transformed to the respective response scales and then multiplied together. To do this for the zero-inflation part, I copied the code from the base R <code>binomial()</code> family with the appropriate link specified. The second line of code below adds the predicted values for each combination of <code>mined</code> and <code>spp</code> to the prediction data object. Note that each component is back-transformed using the appropriate link, and then multiplied together.
</p>
<figure>
<pre><code data-lang="r"><span>ilink</span> <span>&lt;-</span> <span>binomial</span><span>(</span><span>link</span> <span>=</span> <span>&quot;cloglog&quot;</span><span>)</span><span>$</span><span>linkinv</span>
<span>newd</span> <span>&lt;-</span> <span>transform</span><span>(</span><span>newd</span><span>,</span> <span>fitted</span> <span>=</span> <span>exp</span><span>(</span><span>pred</span><span>[,</span><span>1</span><span>])</span> <span>*</span> <span>ilink</span><span>(</span><span>pred</span><span>[,</span><span>2</span><span>]))</span></code></pre>
</figure>
<p>
A plot of the predicted values is then easily produced
</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>newd</span><span>,</span> <span>aes</span><span>(</span><span>x</span> <span>=</span> <span>spp</span><span>,</span> <span>y</span> <span>=</span> <span>fitted</span><span>,</span> <span>colour</span> <span>=</span> <span>mined</span><span>))</span> <span>+</span>
    <span>geom_point</span><span>()</span></code></pre>
</figure>
<figure>
<img src="https://i2.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/compare-mgcv-with-glmmtmb-salamander-population-mode-plot-3-1.png?w=456" data-recalc-dims="1"><br>
</figure>
<p>
Because of the way the <code>gam()</code> model is implemented, I could also have computed the Bayesian credible intervals using the Bayesian covariance matrix of the model parameters via the <code>se.fit</code> argument to <code>predict()</code>. I’ll perhaps save that for another day…
</p>
<h3>
Owls<br>
</h3>
<p>
The <code>Owls</code> data are also available in the <strong>glmmTMB</strong> package, which I load and then do a little processing of the data to simplify the name of the response variable and to mean centre the <code>ArrivalTime</code> covariate.
</p>
<figure>
<pre><code data-lang="r"><span>data</span><span>(</span><span>Owls</span><span>,</span> <span>package</span> <span>=</span> <span>&quot;glmmTMB&quot;</span><span>)</span>
<span>names</span><span>(</span><span>Owls</span><span>)</span> <span>&lt;-</span> <span>sub</span><span>(</span><span>&quot;SiblingNegotiation&quot;</span><span>,</span> <span>&quot;NCalls&quot;</span><span>,</span> <span>names</span><span>(</span><span>Owls</span><span>))</span>
<span>Owls</span> <span>&lt;-</span> <span>transform</span><span>(</span><span>Owls</span><span>,</span> <span>cArrivalTime</span> <span>=</span> <span>ArrivalTime</span> <span>-</span> <span>mean</span><span>(</span><span>ArrivalTime</span><span>))</span></code></pre>
</figure>
<p>
Two ZIP models are considered
</p>
<ol type="1">
<li>
a ZIP with constant zero-inflation (an intercept-only model for the zero-inflation), and
</li>
<li>
a ZIP with complex zero-inflation, where one covariate and a random effect for <code>Nest</code> are included in the linear predictor of the zero-inflation part of the model.
</li>
</ol>
<p>
The constant zero-inflation models are fitted using the <code>ziformula</code> argument for <code>glmmTMB</code> with <code>family = poisson</code>, whilst for <code>gam()</code> we use a list of two formula objects, the second for the ZI linear predictor, and the <code>ziplss</code> family. Note that this model could also be fitted using the <code>Zip()</code> family in <strong>mgcv</strong> but that employs a different, simpler fitting algorithm so to facilitate comparison with the more complex model I use <code>ziplss()</code> instead.
</p>
<figure>
<pre><code data-lang="r"><span>m1.tmb</span> <span>&lt;-</span> <span>glmmTMB</span><span>(</span><span>NCalls</span> <span>~</span> <span>(</span><span>FoodTreatment</span> <span>+</span> <span>cArrivalTime</span><span>)</span> <span>*</span> <span>SexParent</span> <span>+</span> <span>offset</span><span>(</span><span>logBroodSize</span><span>)</span> <span>+</span> <span>(</span><span>1</span> <span>|</span> <span>Nest</span><span>),</span>
                  <span>ziformula</span> <span>=</span> <span>~</span> <span>1</span><span>,</span> <span>data</span> <span>=</span> <span>Owls</span><span>,</span> <span>family</span> <span>=</span> <span>poisson</span><span>)</span>
<span>m1.gam</span> <span>&lt;-</span> <span>gam</span><span>(</span><span>list</span><span>(</span><span>NCalls</span> <span>~</span> <span>(</span><span>FoodTreatment</span> <span>+</span> <span>cArrivalTime</span><span>)</span> <span>*</span> <span>SexParent</span> <span>+</span> <span>offset</span><span>(</span><span>logBroodSize</span><span>)</span> <span>+</span> <span>s</span><span>(</span><span>Nest</span><span>,</span> <span>bs</span> <span>=</span> <span>&quot;re&quot;</span><span>),</span>
                   <span>~</span> <span>1</span><span>),</span>
              <span>data</span> <span>=</span> <span>Owls</span><span>,</span> <span>family</span> <span>=</span> <span>ziplss</span><span>())</span></code></pre>
</figure>
<p>
Again note that these models are not estimated in the same way; <code>glmmTMB()</code> estimates the model parameters using maximum likelihood, whilst only REML estimation is available for the <code>ziplss()</code> family with <code>gam()</code>. In <code>gam()</code>, the intercept-only ZI linear predictor is specified with the formula <code>~ 1</code>.
</p>
<p>
To compare the estimates of the model coefficients I wrote a little function to extract the estimated values and their standard errors from the two model objects
</p>
<figure>
<pre><code data-lang="r"><span>createCoeftab</span> <span>&lt;-</span> <span>function</span><span>(</span><span>TMB</span><span>,</span> <span>GAM</span><span>,</span> <span>GAMrange</span><span>)</span> <span>{</span>
    <span>bTMB</span> <span>&lt;-</span> <span>fixef</span><span>(</span><span>TMB</span><span>)</span><span>$</span><span>cond</span><span>[</span><span>-1</span><span>]</span>
    <span>bGAM</span> <span>&lt;-</span> <span>coef</span><span>(</span><span>GAM</span><span>)[</span><span>GAMrange</span><span>]</span>
    <span>seTMB</span> <span>&lt;-</span> <span>diag</span><span>(</span><span>vcov</span><span>(</span><span>TMB</span><span>)</span><span>$</span><span>cond</span><span>)[</span><span>-1</span><span>]</span>
    <span>seGAM</span> <span>&lt;-</span> <span>diag</span><span>(</span><span>vcov</span><span>(</span><span>GAM</span><span>))[</span><span>GAMrange</span><span>]</span>
    <span>nms</span> <span>&lt;-</span> <span>names</span><span>(</span><span>bTMB</span><span>)</span>
    <span>nms</span> <span>&lt;-</span> <span>sub</span><span>(</span><span>&quot;FoodTreatment&quot;</span><span>,</span> <span>&quot;FT&quot;</span><span>,</span> <span>nms</span><span>)</span>
    <span>nms</span> <span>&lt;-</span> <span>sub</span><span>(</span><span>&quot;cArrivalTime&quot;</span><span>,</span> <span>&quot;ArrivalTime&quot;</span><span>,</span> <span>nms</span><span>)</span>
    <span>df</span> <span>&lt;-</span> <span>data.frame</span><span>(</span><span>model</span>    <span>=</span> <span>rep</span><span>(</span><span>c</span><span>(</span><span>&quot;glmmTMB&quot;</span><span>,</span> <span>&quot;mgcv::gam&quot;</span><span>),</span> <span>each</span> <span>=</span> <span>length</span><span>(</span><span>bGAM</span><span>)),</span>
                     <span>term</span>     <span>=</span> <span>rep</span><span>(</span><span>nms</span><span>,</span> <span>2</span><span>),</span>
                     <span>estimate</span> <span>=</span> <span>unname</span><span>(</span><span>c</span><span>(</span><span>bTMB</span><span>,</span> <span>bGAM</span><span>)))</span>
    <span>df</span> <span>&lt;-</span> <span>transform</span><span>(</span><span>df</span><span>,</span>
                    <span>upper</span> <span>=</span> <span>estimate</span> <span>+</span> <span>sqrt</span><span>(</span><span>c</span><span>(</span><span>seTMB</span><span>,</span> <span>seGAM</span><span>)),</span>
                    <span>lower</span> <span>=</span> <span>estimate</span> <span>-</span> <span>sqrt</span><span>(</span><span>c</span><span>(</span><span>seTMB</span><span>,</span> <span>seGAM</span><span>)))</span>
    <span>df</span>
<span>}</span></code></pre>
</figure>
<p>
Passing each of the models to <code>createCoeftab()</code>
</p>
<figure>
<pre><code data-lang="r"><span>m1.coefs</span> <span>&lt;-</span> <span>createCoeftab</span><span>(</span><span>m1.tmb</span><span>,</span> <span>m1.gam</span><span>,</span> <span>GAMrange</span> <span>=</span> <span>2</span><span>:</span><span>6</span><span>)</span></code></pre>
</figure>
<p>
results in a tidy data frame suitable for plotting with <code>ggplot()</code>.
</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>m1.coefs</span><span>,</span> <span>aes</span><span>(</span><span>x</span> <span>=</span> <span>estimate</span><span>,</span> <span>y</span> <span>=</span> <span>term</span><span>,</span> <span>colour</span> <span>=</span> <span>model</span><span>,</span> <span>shape</span> <span>=</span> <span>model</span><span>,</span> <span>xmax</span> <span>=</span> <span>upper</span><span>,</span> <span>xmin</span> <span>=</span> <span>lower</span><span>))</span> <span>+</span>
    <span>geom_pointrangeh</span><span>(</span><span>position</span> <span>=</span> <span>position_dodgev</span><span>(</span><span>height</span> <span>=</span> <span>0.3</span><span>))</span> <span>+</span>
    <span>labs</span><span>(</span><span>y</span> <span>=</span> <span>NULL</span><span>,</span>
         <span>x</span> <span>=</span> <span>&quot;Regression estimate&quot;</span><span>,</span>
         <span>title</span> <span>=</span> <span>&quot;Comparing mgcv with glmmTMB&quot;</span><span>,</span>
         <span>subtitle</span> <span>=</span> <span>&quot;Owls: ZIP with constant zero-inflation&quot;</span><span>,</span>
         <span>caption</span> <span>=</span> <span>&quot;Bars are ±1 SE&quot;</span><span>)</span></code></pre>
</figure>
<figure>
<img alt="Comparison of estimated model fixed effect parameters for the constant zer-inflation model fitted to the owl nestling behaviour data." src="https://i2.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/compare-mgcv-with-glmmtmb-plot-m1-coefs-1.png?w=456" data-recalc-dims="1">
<figcaption>
Comparison of estimated model fixed effect parameters for the constant zer-inflation model fitted to the owl nestling behaviour data.<br>
</figcaption>
</figure>
<p>
As can be seen in the figure, the estimates from the two functions are quite similar.
</p>
<p>
The more-complex models with covariates in the ZI linear predictor are fitted next
</p>
<figure>
<pre><code data-lang="r"><span>m2.tmb</span> <span>&lt;-</span> <span>glmmTMB</span><span>(</span><span>NCalls</span> <span>~</span> <span>(</span><span>FoodTreatment</span> <span>+</span> <span>cArrivalTime</span><span>)</span> <span>*</span> <span>SexParent</span> <span>+</span>
                      <span>offset</span><span>(</span><span>logBroodSize</span><span>)</span> <span>+</span> <span>(</span><span>1</span> <span>|</span> <span>Nest</span><span>),</span>
                  <span>ziformula</span> <span>=</span> <span>~</span> <span>FoodTreatment</span> <span>+</span> <span>(</span><span>1</span> <span>|</span> <span>Nest</span><span>),</span> <span>data</span> <span>=</span> <span>Owls</span><span>,</span> <span>family</span> <span>=</span> <span>poisson</span><span>)</span>
<span>m2.gam</span> <span>&lt;-</span> <span>gam</span><span>(</span><span>list</span><span>(</span><span>NCalls</span> <span>~</span> <span>(</span><span>FoodTreatment</span> <span>+</span> <span>cArrivalTime</span><span>)</span> <span>*</span> <span>SexParent</span> <span>+</span>
                       <span>offset</span><span>(</span><span>logBroodSize</span><span>)</span> <span>+</span> <span>s</span><span>(</span><span>Nest</span><span>,</span> <span>bs</span> <span>=</span> <span>&quot;re&quot;</span><span>),</span>
                   <span>~</span> <span>FoodTreatment</span> <span>+</span> <span>s</span><span>(</span><span>Nest</span><span>,</span> <span>bs</span> <span>=</span> <span>&quot;re&quot;</span><span>)),</span>
              <span>data</span> <span>=</span> <span>Owls</span><span>,</span> <span>family</span> <span>=</span> <span>ziplss</span><span>())</span></code></pre>
</figure>
<p>
As before, we gather the model coefficients
</p>
<figure>
<pre><code data-lang="r"><span>m2.coefs</span> <span>&lt;-</span> <span>createCoeftab</span><span>(</span><span>m2.tmb</span><span>,</span> <span>m2.gam</span><span>,</span> <span>GAMrange</span> <span>=</span> <span>2</span><span>:</span><span>6</span><span>)</span></code></pre>
</figure>
<p>
and plot them
</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>m2.coefs</span><span>,</span> <span>aes</span><span>(</span><span>x</span> <span>=</span> <span>estimate</span><span>,</span> <span>y</span> <span>=</span> <span>term</span><span>,</span> <span>colour</span> <span>=</span> <span>model</span><span>,</span> <span>shape</span> <span>=</span> <span>model</span><span>,</span>
                     <span>xmax</span> <span>=</span> <span>upper</span><span>,</span> <span>xmin</span> <span>=</span> <span>lower</span><span>))</span> <span>+</span>
    <span>geom_pointrangeh</span><span>(</span><span>position</span> <span>=</span> <span>position_dodgev</span><span>(</span><span>height</span> <span>=</span> <span>0.3</span><span>))</span> <span>+</span>
    <span>labs</span><span>(</span><span>y</span> <span>=</span> <span>NULL</span><span>,</span>
         <span>x</span> <span>=</span> <span>&quot;Regression estimate&quot;</span><span>,</span>
         <span>title</span> <span>=</span> <span>&quot;Comparing mgcv with glmmTMB&quot;</span><span>,</span>
         <span>subtitle</span> <span>=</span> <span>&quot;Owls: ZIP with complex zero-inflation&quot;</span><span>,</span>
         <span>caption</span> <span>=</span> <span>&quot;Bars are ±1 SE&quot;</span><span>)</span></code></pre>
</figure>
<figure>
<img alt="Comparison of estimated model fixed effect parameters for the complex zer-inflation model fitted to the owl nestling behaviour data." src="https://i0.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/compare-mgcv-with-glmmtmb-plot-m2-coefs-1.png?w=456" data-recalc-dims="1">
<figcaption>
Comparison of estimated model fixed effect parameters for the complex zer-inflation model fitted to the owl nestling behaviour data.<br>
</figcaption>
</figure>
<p>
and likewise as before, the estimates of the fixed effect terms are very similar indeed.
</p>
<h3>
Conclusions<br>
</h3>
<p>
The comparisons shown above show that <code>mgcv::gam()</code> and <code>glmmTMB()</code> produce very similar estimates for the two models. And some crude timings showed that <code>gam()</code> was 20–40% faster than <code>glmmTMB()</code> at fitting the examples discussed in the paper. So all is roses, right!? Who needs <code>glmmTMB()</code>?
</p>
<p>
That would however, be totally the wrong message to take from this comparison. Most notably, and something that isn’t surfaced in these simple examples is that <code>gam()</code> is limited in the complexity of the random effects it can efficiently represent in models:
</p>
<ul>
<li>
it can’t do correlated random effects for random slopes and intercepts models (as far as I can tell anyway), and, and this is probably the deal breaker,
</li>
<li>
model fitting with <code>gam()</code> gets bogged down quickly if the number of levels in a random effect gets large. <a target="_blank" rel="nofollow" href="https://twitter.com/jaimedash">Jamie Ashander</a> did some quick tests with a larger version of the Salamander with 100s of <code>site</code>s and <code>glmmTMB()</code> totally dominated <code>gam()</code>.
</li>
</ul>
<p>
And that’s fine; <code>gam()</code> was not designed to fit GLMMs — there are no less than <strong>three</strong> implementations <em>by Simon Wood alone</em> of functions to fit GAMs with complex random effects in mixed model software (<code>gamm()</code> to fit with <code>lme()</code>, <code>gamm4()</code> to fit using <code>lmer()</code> and <code>glmer()</code>, and <code>jagam()</code> in <strong>mgcv</strong> to fit via JAGS). Furthermore, <code>glmmTMB()</code> is currently more flexible in the range of models that it can fit than any these implementations, except for JAGS, because the <code>nb</code>, <code>Zip</code>, and <code>ziplss</code> families only work with <code>gam()</code>.
</p>
<p>
What the above comparison illustrates, however, is that if you either don’t have complex or many random effects or that you don’t mind running models overnight, <code>gam()</code> is a good option for fitting GLMMs. Plus you have the advantage of estimating smooth functions of covariates, which is one area where <code>glmmTMB()</code> is currently very lacking compared to <code>gam()</code>.
</p>
<p>
That said, it should be possible to emulate what Paul-Christian Bürkner has done in his <a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=brms"><strong>brms</strong> package</a> (and similar implementations by Simon Wood in <code>gamm4()</code>) to use <strong>mgcv</strong> to set up the correct model matrices for the random effect representation of splines which can then be fitted using <code>glmmTMB()</code>.
</p>
<p>
Finally, this was a fun exercise to replicate the analyses in <span data-cites="Brooks2017-so">Brooks et al. (2017)</span>, motivated by a desire to understand what <strong>mgcv</strong> and <code>gam()</code> are doing with these random effect splines. It wasn’t intended as a prize-fight between two title contenders — hopefully this write-up didn’t come across that way. I also learned a lot more about <strong>glmmTMB</strong>, which is shaping up nicely and looks like it’ll have a place in my modelling toolbox.
</p>
<div>
<h3>
References<br>
</h3>
<div>
<p>
Bolker, B. M., Gardner, B., Maunder, M., Berg, C. W., Brooks, M., Comita, L., et al. (2013). Strategies for fitting nonlinear ecological models in r, AD model builder, and BUGS. <em>Methods in ecology and evolution / British Ecological Society</em> 4, 501–512. doi:<a target="_blank" rel="nofollow" href="http://doi.org/10.1111/2041-210X.12044">http://doi.org/10.1111/2041-210X.12044</a>.
</p>
</div>
<div>
<p>
Brooks, M. E., Kristensen, K., Benthem, K. J. van, Magnusson, A., Berg, C. W., Nielsen, A., et al. (2017). Modeling Zero-Inflated count data with glmmTMB. <em>bioRxiv</em>, 132753. doi:<a target="_blank" rel="nofollow" href="http://doi.org/10.1101/132753">http://doi.org/10.1101/132753</a>.
</p>
</div>
<div>
<p>
Price, S. J., Muncy, B. L., Bonner, S. J., Drayer, A. N., and Barton, C. D. (2015). Data from: Effects of mountaintop removal mining and valley filling on the occupancy and abundance of stream salamanders. doi:<a target="_blank" rel="nofollow" href="http://doi.org/10.5061/dryad.5m8f6">http://doi.org/10.5061/dryad.5m8f6</a>.
</p>
</div>
<div>
<p>
Price, S. J., Muncy, B. L., Bonner, S. J., Drayer, A. N., and Barton, C. D. (2016). Effects of mountaintop removal mining and valley filling on the occupancy and abundance of stream salamanders. <em>The Journal of applied ecology</em> 53, 459–468. doi:<a target="_blank" rel="nofollow" href="http://doi.org/10.1111/1365-2664.12585">http://doi.org/10.1111/1365-2664.12585</a>.
</p>
</div>
<div>
<p>
Roulin, A., and Bersier, L.-F. (2007). Nestling barn owls beg more intensely in the presence of their mother than in the presence of their father. <em>Animal behaviour</em> 74, 1099–1106.
</p>
</div>
<div>
<p>
Wood, S. N., Pya, N., and Säfken, B. (2016). Smoothing parameter and model selection for general smooth models. <em>Journal of the American Statistical Association</em> 111, 1548–1563. doi:<a target="_blank" rel="nofollow" href="http://doi.org/10.1080/01621459.2016.1180986">http://doi.org/10.1080/01621459.2016.1180986</a>.
</p>
</div>
<div>
<p>
Zuur, A. F., Ieno, E. N., Walker, N., Saveliev, A. A., and Smith, G. M. (2009). <em>Mixed effects models and extensions in ecology with r:</em> Springer New York doi:<a target="_blank" rel="nofollow" href="http://doi.org/10.1007/978-0-387-87458-6">http://doi.org/10.1007/978-0-387-87458-6</a>.
</p>
</div>
</div>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.fromthebottomoftheheap.net/2017/05/04/compare-mgcv-with-glmmTMB/"> From the Bottom of the Heap - R</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/pEGT1oNd1xc"></a></div>