<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/comparing-smooths-in-factor-smooth-interactions-i/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.fromthebottomoftheheap.net/2017/10/10/difference-splines-i/"> 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>
One of the really appealing features of the <strong>mgcv</strong> package for fitting GAMs is the functionality it exposes for fitting quite complex models, models that lie well beyond what many of us may have learned about what GAMs can do. One of those features that I use a lot is the ability to model the smooth effects of some covariate <span>(x)</span> in the different levels of a factor. Having estimated a separate smoother for each level of the factor, the obvious question is, which smooths are different? In this post I’ll take a look at one way to do this using <code>by</code>-variable smooths.
</p>
<p>
With <strong>mgcv</strong>, smooths are included in model formulae using the <code>s()</code> function. If you want to have the smooth equivalent of a continuous-factor interaction, one way to achieve this is via the <code>by</code> argument to <code>s()</code>. If you pass a factor to <code>by</code>, <strong>mgcv</strong> sets up the model matrix in such a way that you get a separate smoother for each level of the <code>by</code> factor. Each of these smoothers gets its own smoothness parameter — so you can fit a wiggly function in level <em>foo</em> and a smooth function in level <em>bar</em>, with each level’s function being learned from the data associated with that level.
</p>
<p>
I used this technique in a <a rel="nofollow" href="http://doi.org/10.1016/j.gca.2010.12.026" target="_blank">paper</a> I wrote with my colleagues at UCL, Neil Rose, Handong Yang, and Simon Turner <span data-cites="Rose2012-pl">(Rose et al., 2012)</span>. Neil, Handong, and Simon had collected sediment cores from several Scottish lochs and measured metal concentrations, especially of lead (Pb) and mercury (Hg), in sediment slices covering the last 200 years. The aim of the study was to investigate sediment profiles of these metals in three regions of Scotland; north east, north west, and south west. A pair of lochs in each region was selected, one in a catchment with visibly eroding peat/soil, and the other in a catchment without erosion. The different regions represented variations in historical deposition levels, whilst the hypothesis was that cores from eroded and non-eroded catchments would show differential responses to reductions in emissions of Pb and Hg to the atmosphere. The difference, it was hypothesised, was that the eroding soil acts as a secondary source of pollutants to the lake. You can read more about it in the <a rel="nofollow" href="http://doi.org/10.1016/j.gca.2010.12.026" target="_blank">paper</a> — if you’re interested but don’t have access to the journal, send me an email and I’ll pass on a pdf.
</p>
<p>
It was relatively simple to fit splines to each sediment profile, but once I’d done this, how were we going to estimate the difference between the fitted trends? Thankfully, I already had the answer as Simon Wood had supplied code to do it to an OP on the R-Help listserver some years previous. That answer involved <code>by</code>-variable smoothers, which I was already using, and the use of the <span>(Xp)</span> matrix of the fitted GAM.
</p>
<p>
Readers of this blog will have heard about the <span>(Xp)</span> matrix before; it’s used a lot when we want to simulate from the posterior of the estimated model. Importantly, for our purposes, it allows for the creation of derived quantities, from the fitted model, and the assignment of uncertainty to those quantities.
</p>
<p>
In this post I’ll illustrate how to do the required comparison using some of the data from that study on Scottish lochs.
</p>
<p>
In this post I’ll use the the following packages
</p>
<ul>
<li>
<strong>readr</strong>
</li>
<li>
<strong>dplyr</strong>
</li>
<li>
<strong>ggplot2</strong>, and
</li>
<li>
<strong>mgcv</strong>
</li>
</ul>
<p>
You’ll more than likely have these installed, but if you get errors about missing packages when you run the code chunk below, install any missing packages and run the chunk again
</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>'readr'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'dplyr'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'ggplot2'</span><span>)</span><span>
</span><span>theme_set</span><span>(</span><span>theme_bw</span><span>())</span><span>
</span><span>library</span><span>(</span><span>'mgcv'</span><span>)</span></code></pre>
</figure>
<p>
Next, load the data set and convert the <code>SiteCode</code> variable to a factor for use in fitting the GAM with <code>gam()</code>
</p>
<figure>
<pre><code data-lang="r"><span>uri</span><span> </span><span>&lt;-</span><span> </span><span>'https://gist.githubusercontent.com/gavinsimpson/eb4ff24fa9924a588e6ee60dfae8746f/raw/geochimica-metals.csv'</span><span>
</span><span>metals</span><span> </span><span>&lt;-</span><span> </span><span>read_csv</span><span>(</span><span>uri</span><span>,</span><span> </span><span>skip</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>col_types</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>'ciccd'</span><span>))</span><span>
</span><span>metals</span><span> </span><span>&lt;-</span><span> </span><span>mutate</span><span>(</span><span>metals</span><span>,</span><span> </span><span>SiteCode</span><span> </span><span>=</span><span> </span><span>factor</span><span>(</span><span>SiteCode</span><span>))</span></code></pre>
</figure>
<p>
This is a subset of the data used in <span data-cites="Rose2012-pl">Rose et al. (2012)</span> — the Hg concentrations in the sediments for just three of the lochs are included here in the interests of simplicity. The data set contains 5 variables
</p>
<figure>
<pre><code data-lang="r"><span>metals</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text"># A tibble: 44 x 5
   SiteCode  Date SoilType Region        Hg
     </code></pre></figure>