<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/fitting-gams-with-brms-part-1/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://www.fromthebottomoftheheap.net/2018/04/21/fitting-gams-with-brms/"> 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>
Regular readers will know that I have a somewhat unhealthy relationship with GAMs and the <strong>mgcv</strong> package. I use these models all the time in my research but recently we’ve been hitting the limits of the range of models that <strong>mgcv</strong> can fit. So I’ve been looking into alternative ways to fit the GAMs I want to fit but which can handle the kinds of data or distributions that have been cropping up in our work. The <strong>brms</strong> package <span data-cites="brms-2017">(Bürkner, 2017)</span> is an excellent resource for modellers, providing a high-level R front end to a vast array of model types, all fitted using <a rel="nofollow" href="http://mc-stan.org/" target="_blank">Stan</a>. <strong>brms</strong> is the perfect package to go beyond the limits of <strong>mgcv</strong> because <strong>brms</strong> even uses the smooth functions provided by <strong>mgcv</strong>, making the transition easier. In this post I take a look at how to fit a simple GAM in <strong>brms</strong> and compare it with the same model fitted using <strong>mgcv</strong>.
</p>
<p>
In this post we’ll use the following packages. If you don’t know <strong>schoenberg</strong>, it’s a package I’m writing to provide <code>ggplot</code> versions of plots that can be produced by <strong>mgcv</strong> from fitted GAM objects. <strong>schoenberg</strong> is in early development, but it currentl works well enough to plot the models we fit here. If you’ve never come across this package before, you can install it from Github using <code>devtools::install_github(‘gavinsimpson/schoenberg’)</code>
</p>
<figure>
<pre><code data-lang="r"><span>## packages
</span><span>library</span><span>(</span><span>'mgcv'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'brms'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'ggplot2'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'schoenberg'</span><span>)</span><span>
</span><span>theme_set</span><span>(</span><span>theme_bw</span><span>())</span></code></pre>
</figure>
<p>
To illustrate <strong>brms</strong> GAM-fittgin chops, we’ll use the <code>mcycle</code> data set that comes with the <strong>MASS</strong> package. It contains a set of measurements of the acceleration force on a rider’s head during a simulated motorcycle collision and the time, in milliseconds, post collision. The data are loaded using <code>data()</code> and we take a look at the first few rows
</p>
<figure>
<pre><code data-lang="r"><span>## load the example data mcycle
</span><span>data</span><span>(</span><span>mcycle</span><span>,</span><span> </span><span>package</span><span> </span><span>=</span><span> </span><span>'MASS'</span><span>)</span><span>
</span><span>## show data
</span><span>head</span><span>(</span><span>mcycle</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">  times accel
1   2.4   0.0
2   2.6  -1.3
3   3.2  -2.7
4   3.6   0.0
5   4.0  -2.7
6   6.2  -2.7</code></pre>
</figure>
<p>
The aim is to model the acceleration force (<code>accel</code>) as a function of time post collision (<code>times</code>). The plot below shows the data.
</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>mcycle</span><span>,</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>times</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>accel</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>()</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Miliseconds post impact&quot;</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Acceleration (g)&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Simulated Motorcycle Accident&quot;</span><span>,</span><span>
</span><span>subtitle</span><span> </span><span>=</span><span> </span><span>&quot;Measurements of head acceleration&quot;</span><span>)</span></code></pre>
</figure>
<figure>
<img src="https://i0.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/fitting-gams-with-brms-plot-data-1.png?w=456&ssl=1" data-recalc-dims="1"><br>
</figure>
<p>
We’ll model acceleration as a <em>smooth</em> function of time using a GAM and the default thin plate regression spline basis. This can be done using the <code>gam()</code> function in <strong>mgcv</strong> ad for comparison with the fully bayesian model we’ll fit shortly, we use `method = “REML” to estimate the smoothness parameter for the spline in mixed model form and REML
</p>
<figure>
<pre><code data-lang="r"><span>m</span><span>1</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>accel</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>times</span><span>),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>mcycle</span><span>,</span><span> </span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;REML&quot;</span><span>)</span><span>
</span><span>summary</span><span>(</span><span>m</span><span>1</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Family: gaussian 
Link function: identity 

Formula:
accel ~ s(times)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -25.546      1.951  -13.09   &lt;2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
           edf Ref.df    F p-value    
s(times) 8.625  8.958 53.4  &lt;2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.783   Deviance explained = 79.7%
-REML = 616.14  Scale est. = 506.35    n = 133</code></pre>
</figure>
<p>
As we can see from the model summary, the estimated smooth uses about 8.5 effective degrees of freedom and in the test of zero effect, the null hypothesis is strongly rejected. The fitted spline explains about 80% of the variance or deviance in the data.
</p>
<p>
To plot the fitted smooth we could use the <code>plot()</code> method provided by <strong>mgcv</strong>, but this uses base graphics. Instead we can use the <code>draw()</code> method from <strong>schoenberg</strong>, which can currently handle most of the univariate smooths in <strong>mgcv</strong> plus 2-d tensor product smooths
</p>
<figure>
<pre><code data-lang="r"><span>draw</span><span>(</span><span>m</span><span>1</span><span>)</span></code></pre>
</figure>
<figure>
<img src="https://i0.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/fitting-gams-with-brms-plot-mgcv-model-1.png?w=456&ssl=1" data-recalc-dims="1"><br>
</figure>
<p>
The equivalent model can be estimated using a fully-bayesian approach via the <code>brm()</code> function in the <strong>brms</strong> package. In fact, <code>brm()</code> will use the smooth specification functions from <strong>mgcv</strong>, making our lives much easier. The major difference though is that you can’t use <code>te()</code> or <code>ti()</code> smooths in <code>brm()</code> models; you need to use <code>t2()</code> tensor product smooths instead. This is because the smooths in the model are going to be treated as random effects and the model estimated as a GLMM, which exploits the duality of splines as random effects. In this representation, the wiggly parts of the spline basis are treated as a random effect and their associated variance parameter controls the degree of wiggliness of the fitted spline. The perfectly smooth parts of the basis are treated as a fixed effect. In this form, the GAM can be estimated using standard GLMM software; it’s what allows the <code>gamm4()</code> function for example to fit GAMMs using the <strong>lme4</strong> package for example. This is also the reason why we can’t use <code>te()</code> or <code>ti()</code> smooths; those smooths do not have nicely separable penalties which means they can’t be written in the form required to be fitted using typical mixed model software.
</p>
<p>
The <code>brm()</code> version of the GAM is fitted using the code below. Note that I have changed a few things from their default values as
</p>
<ol type="1">
<li>
the model required more than the default number of MCMC samples — <code>iter = 4000</code>,
</li>
<li>
the samples needed thinning to deal with some strong autocorrelation in the Markov chains — <code>thin = 10</code>,
</li>
<li>
the <code>adapt.delta</code> parameter, a tuning parameter in he NUTS sampler for Hamiltonian Monte Carlo, potentially needed raising — there was a warning about a potential divergent transition but I should have looked to see if it was one or not; instead I just increased the tuning parameter to <code>0.99</code>,
</li>
<li>
four chains fitted by default but I wanted these to be fitted using 4 CPU <code>cores</code>,
</li>
<li>
<code>seed</code> sets the internal random number generator seed, which allows reproducibility of models, and
</li>
<li>
for this post I didn’t want to print out the progress of the sampler — <code>refresh = 0</code> — typically you won’t want to do this so you can see how the sampling is progressing.
</li>
</ol>
<p>
The rest of the model is pretty similar to the <code>gam()</code> version we fitted earlier. The main difference is that I use the <code>bf()</code> function to create a special <strong>brms</strong> formula specifying the model. You don’t actually need to do this for such a simple model, but in a later post we’ll use this to fit distributional GAMs. Note that I’m leaving all the priors in the model at the default values. I’ll look at defining priors in a later post; for now I’m just going to use the default priors that <code>brm()</code> uses
</p>
<figure>
<pre><code data-lang="r"><span>m</span><span>2</span><span> </span><span>&lt;-</span><span> </span><span>brm</span><span>(</span><span>bf</span><span>(</span><span>accel</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>times</span><span>)),</span><span>
</span><span>data</span><span> </span><span>=</span><span> </span><span>mcycle</span><span>,</span><span> </span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>(),</span><span> </span><span>cores</span><span> </span><span>=</span><span> </span><span>4</span><span>,</span><span> </span><span>seed</span><span> </span><span>=</span><span> </span><span>17</span><span>,</span><span>
</span><span>iter</span><span> </span><span>=</span><span> </span><span>4000</span><span>,</span><span> </span><span>warmup</span><span> </span><span>=</span><span> </span><span>1000</span><span>,</span><span> </span><span>thin</span><span> </span><span>=</span><span> </span><span>10</span><span>,</span><span> </span><span>refresh</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span>
</span><span>control</span><span> </span><span>=</span><span> </span><span>list</span><span>(</span><span>adapt_delta</span><span> </span><span>=</span><span> </span><span>0.99</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Compiling the C++ model</code></pre>
</figure>
<figure>
<pre><code data-lang="text">Start sampling</code></pre>
</figure>
<p>
Once the model has finished compiling and sampling we can output the model summary
</p>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>m</span><span>2</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text"> Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: accel ~ s(times) 
   Data: mcycle (Number of observations: 133) 
Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 10; 
         total post-warmup samples = 1200
    ICs: LOO = NA; WAIC = NA; R2 = NA
 
Smooth Terms: 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sds(stimes_1)   722.44    198.12   450.17  1150.27       1180 1.00

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept   -25.54      2.02   -29.66   -21.50       1200 1.00
stimes_1     16.10     38.20   -61.46    90.91       1171 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma    22.78      1.47    19.94    25.68       1200 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).</code></pre>
</figure>
<p>
This outputs details of the model fitted plus parameter estimates (as posterior means), standard errors, (by default) 95% credible intervals and two other diagnostics:
</p>
<ol type="1">
<li>
<code>Eff.Sample</code> is the effective sample size of the posterior samples in the model, and
</li>
<li>
<code>Rhat</code> is the <em>potential scale reduction factor</em> or Gelman-Rubin diagnostic and is a measure of how well the chains have converged and ideally should be equal to <code>1</code>.
</li>
</ol>
<p>
The summary includes two entries for the smooth of <code>times</code>:
</p>
<ol type="1">
<li>
<code>sds(stimes_1)</code> is the variance parameter, which has the effect of controlling the wiggliness of the smooth — the larger this value the more wiggly the smooth. We can see that the credible interval doesn’t include 0 so there is evidence that a smooth is required over and above a linear parametric effect of <code>times</code>, details of which are given next,
</li>
<li>
<code>stimes_1</code> is the fixed effect part of the spline, which is the linear function that is perfectly smooth.
</li>
</ol>
<p>
The final parameter table includes information on the variance of the data about the conditional mean of the response.
</p>
<p>
How does this model compare with the one fitted using <code>gam()</code>? We can use the <code>gam.vcomp()</code> function to compute the variance component representation of the smooth estimated via <code>gam()</code>. To make it comparable with the value shown for the <strong>brms</strong> model, we don’t undo the rescaling of the penalty matrix that <code>gam()</code> performs to help with numeric stability during model fitting.
</p>
<figure>
<pre><code data-lang="r"><span>gam.vcomp</span><span>(</span><span>m</span><span>1</span><span>,</span><span> </span><span>rescale</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Standard deviations and 0.95 confidence intervals:

           std.dev     lower      upper
s(times) 807.88726 480.66162 1357.88215
scale     22.50229  19.85734   25.49954

Rank: 2/2</code></pre>
</figure>
<p>
This gives a posterior mean of 807.89 with 95% confidence interval of 480.66–1357.88, which compares well with posterior mean and credible interval of the <code>brm()</code> version of 722.44 (450.17 – 1150.27).
</p>
<p>
The <code>marginal_smooths()</code> function is used to extract the marginal effect of the spline.
</p>
<figure>
<pre><code data-lang="r"><span>msms</span><span> </span><span>&lt;-</span><span> </span><span>marginal_smooths</span><span>(</span><span>m</span><span>2</span><span>)</span></code></pre>
</figure>
<p>
This function extracts enough information about the estimated spline to plot it using the <code>plot()</code> method
</p>
<figure>
<pre><code data-lang="r"><span>plot</span><span>(</span><span>msms</span><span>)</span></code></pre>
</figure>
<figure>
<img src="https://i2.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/fitting-gams-with-brms-plot-marginal-smooths-1.png?w=456&ssl=1" data-recalc-dims="1"><br>
</figure>
<p>
Given the similarity in the variance components of the two models it is not surprising the two estimated smooth also look similar. The <code>marginal_smooths()</code> is effectively the equivalent of the <code>plot()</code> method for <strong>mgcv</strong>-based GAMs.
</p>
<p>
There’s a lot that we can and should do to check the model fit. For now, we’ll look at two posterior predictive check plots that <strong>brms</strong>, via the <strong>bayesplot</strong> package <span data-cites="bayesplot-2018">(Gabry and Mahr, 2018)</span>, makes very easy to produce using the <code>pp_check()</code> function.
</p>
<figure>
<pre><code data-lang="r"><span>pp_check</span><span>(</span><span>m</span><span>2</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Using 10 posterior samples for ppc type 'dens_overlay' by default.</code></pre>
</figure>
<figure>
<img src="https://i2.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/fitting-gams-with-brms-pp-check-density-1.png?w=456&ssl=1" data-recalc-dims="1"><br>
</figure>
<p>
The default produces a density plot overlay of the original response values (the thick black line) with 10 draws from the posterior distribution of the model. If the model is a good fit to the data, samples of data sampled from it at the observed values of the covariate(s) should be similar to one another.
</p>
<p>
Another type of posterior predictive check plot is the empirical cumulative distribution function of the observations and random draws from the model posterior, which we can produce with <code>type = “ecdf_overlay”</code>
</p>
<figure>
<pre><code data-lang="r"><span>pp_check</span><span>(</span><span>m</span><span>2</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>&quot;ecdf_overlay&quot;</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Using 10 posterior samples for ppc type 'ecdf_overlay' by default.</code></pre>
</figure>
<figure>
<img src="https://i2.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/fitting-gams-with-brms-pp-check-ecdf-1.png?w=456&ssl=1" data-recalc-dims="1"><br>
</figure>
<p>
Both plots show significant deviations between the the posterior simulations and the observed data. The poor posterior predictive check results are in large part due to the non-constant variance of the acceleration data conditional upon the covariate. Both models assumed that the observation are distributed Gaussian with means equal to the fitted values (estimated expectation of the response) with the same variance <span>(^2)</span>. The observations appear to have different variances, which we can model with a distributional model, which allow all parameters of the distribution of the response to be modelled with linear predictors. We’ll take a look at these models in a future post.
</p>
<h3>
References<br>
</h3>
<div>
<div>
<p>
Bürkner, P.-C. (2017). brms: An R package for bayesian multilevel models using Stan. <em>Journal of Statistical Software</em> 80, 1–28. doi:<a rel="nofollow" href="https://doi.org/10.18637/jss.v080.i01" target="_blank">10.18637/jss.v080.i01</a>.
</p>
</div>
<div>
<p>
Gabry, J., and Mahr, T. (2018). <em>Bayesplot: Plotting for bayesian models</em>. Available at: <a rel="nofollow" href="https://cran.r-project.org/package=bayesplot" target="_blank">https://CRAN.R-project.org/package=bayesplot</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="https://www.fromthebottomoftheheap.net/2018/04/21/fitting-gams-with-brms/"> 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 rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/_-h19fLVzJo" height="1"></a></div>