<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/introducing-gratia/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://www.fromthebottomoftheheap.net/2018/10/23/introducing-gratia/"> 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>
I use generalized additive models (GAMs) in my research work. I use them a lot! Simon Wood’s <strong>mgcv</strong> package is an excellent set of software for specifying, fitting, and visualizing GAMs for very large data sets. Despite recently dabbling with <strong>brms</strong>, <strong>mgcv</strong> is still my go-to GAM package. The only down-side to <strong>mgcv</strong> is that it is not very tidy-aware and the <strong>ggplot</strong>-verse may as well not exist as far as it is concerned. This in itself is no bad thing, though as someone who uses <strong>mgcv</strong> a lot but also prefers to do my plotting with <strong>ggplot2</strong>, this lack of awareness was starting to hurt. So, I started working on something to help bridge the gap between these two separate worlds that I inhabit. The fruit of that labour is <strong>gratia</strong>, and development has progressed to the stage where I am ready to talk a bit more about it.
</p>
<p>
<strong>gratia</strong> is an R package for working with GAMs fitted with <code>gam()</code>, <code>bam()</code> or <code>gamm()</code> from <strong>mgcv</strong> or <code>gamm4()</code> from the <strong>gamm4</strong> package, although functionality for handling the latter is not yet implement. <strong>gratia</strong> provides functions to replace the base-graphics-based <code>plot.gam()</code> and <code>gam.check()</code> that <strong>mgcv</strong> provides with <strong>ggplot2</strong>-based versions. Recent changes have also resulted in <strong>gratia</strong> being much more <strong>tidyverse</strong> aware and it now (mostly) returns outputs as tibbles.
</p>
<p>
In this post I wanted to give a flavour of what is currently possible with <strong>gratia</strong> and outline what still needs to be implemented. <span></span>
</p>
<p>
<strong>gratia</strong> currently lives on GitHub, so we need to install it from there using <code>devtools::install_github</code>:
</p>
<figure>
<pre><code data-lang="r"><span>devtools</span><span>::</span><span>install_github</span><span>(</span><span>'gavinsimpson/gratia'</span><span>)</span></code></pre>
</figure>
<p>
To do anything useful with <strong>gratia</strong> we need a GAM and for that we need <strong>mgcv</strong>
</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>'mgcv'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'gratia'</span><span>)</span></code></pre>
</figure>
<p>
and an old favourite example data set
</p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>20</span><span>)</span><span>
</span><span>dat</span><span> </span><span>&lt;-</span><span> </span><span>gamSim</span><span>(</span><span>1</span><span>,</span><span> </span><span>n</span><span> </span><span>=</span><span> </span><span>400</span><span>,</span><span> </span><span>dist</span><span> </span><span>=</span><span> </span><span>&quot;normal&quot;</span><span>,</span><span> </span><span>scale</span><span> </span><span>=</span><span> </span><span>2</span><span>,</span><span> </span><span>verbose</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>)</span></code></pre>
</figure>
<p>
The simulated data in <code>dat</code> are well-studied in GAM-related research and contain a number of covariates — labelled <code>x0</code> through <code>x3</code> — which have, to varying degrees, non-linear relationships with the response. We want to try to recover these relationships by approximating the true relationships between covariate and response using splines. To fit a purely additive model, we use
</p>
<figure>
<pre><code data-lang="r"><span>mod</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>y</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>x</span><span>0</span><span>)</span><span> </span><span>+</span><span> </span><span>s</span><span>(</span><span>x</span><span>1</span><span>)</span><span> </span><span>+</span><span> </span><span>s</span><span>(</span><span>x</span><span>2</span><span>)</span><span> </span><span>+</span><span> </span><span>s</span><span>(</span><span>x</span><span>3</span><span>),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>dat</span><span>,</span><span> </span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;REML&quot;</span><span>)</span></code></pre>
</figure>
<p>
<strong>mgcv</strong> provides a <code>summary()</code> method that is used to extract information about the fitted GAM
</p>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>mod</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Family: gaussian 
Link function: identity 

Formula:
y ~ s(x0) + s(x1) + s(x2) + s(x3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.7625     0.0959   80.95   &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(x0) 3.528  4.370  11.30 4.2e-09 ***
s(x1) 2.662  3.310 129.02 &lt; 2e-16 ***
s(x2) 8.146  8.799  84.72 &lt; 2e-16 ***
s(x3) 1.001  1.002   0.00   0.987    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.763   Deviance explained = 77.2%
-REML = 850.87  Scale est. = 3.6785    n = 400</code></pre>
</figure>
<p>
and the <code>k.check()</code> function for checking whether sufficient numbers of basis functions were used in each smooth in the model. (You may not have used <code>k.check()</code> directly — it is called by <code>gam.check()</code> which prints out other diagnostics and also produces four model diagnostic plots, which is one thing that <strong>gratia</strong> provides a replacement for.)
</p>
<h3>
Plotting smooths<br>
</h3>
<p>
To visualize estimated GAMs, <strong>mgcv</strong> provides the <code>plot.gam()</code> method and the <code>vis.gam()</code> function. <strong>gratia</strong> currently provides a <strong>ggplot2</strong>-based replacement for <code>plot.gam()</code>. Work is on-going to provide <code>vis.gam()</code>-like functionality within <strong>gratia</strong> — see <code>?gratia::data_slice</code> for early work in that direction. In <strong>gratia</strong>, we use the <code>draw()</code> generic to produce <strong>ggplot2</strong>-like plots from objects. To visualize the four estimated smooth functions in the GAM <code>mod</code>, we would use
</p>
<figure>
<pre><code data-lang="r"><span>draw</span><span>(</span><span>mod</span><span>)</span></code></pre>
</figure>
<figure>
<img alt="The result of draw(mod) is a plot of each of the four smooth functions in the mod GAM." src="https://i1.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/introducing-gratia-draw-mod-1.png?w=456&ssl=1" data-recalc-dims="1"><figcaption>
The result of <code>draw(mod)</code> is a plot of each of the four smooth functions in the <code>mod</code> GAM.<br>
</figcaption></figure>
<p>
Internally <code>draw()</code> uses the <code>plot_grid()</code> function from <strong>cowplot</strong> to draw multiple panels on the plot device, and to line up the individual plots.
</p>
<p>
There’s not an awful lot more you can do with this now, but the at least the plot is reasonably pretty. <strong>gratia</strong> includes tools for working with the underlying smooths represented in <code>mod</code>, and if you wanted to extract most of the data used to build the plot you’d use the <code>evaluate_smooth()</code> function.
</p>
<figure>
<pre><code data-lang="r"><span>evaluate_smooth</span><span>(</span><span>mod</span><span>,</span><span> </span><span>&quot;x1&quot;</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text"># A tibble: 100 x 5
   smooth fs_variable       x1   est    se
   </code></pre></figure>