<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-ii/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://www.fromthebottomoftheheap.net/2017/12/14/difference-splines-ii/"> 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>
In a <a rel="nofollow" href="https://www.fromthebottomoftheheap.net/2017/10/10/difference-splines-i/" target="_blank">previous post</a> I looked at an approach for computing the differences between smooths estimated as part of a factor-smooth interaction using <code>s()</code>’s <code>by</code> argument. When a common-or-garden factor variable is passed to <code>by</code>, <code>gam()</code> estimates a separate smooth for each <em>level</em> of the <code>by</code> factor. Using the <span>(Xp)</span> matrix approach, we previously saw that we can post-process the model to generate estimates for pairwise differences of smooths. However, the <code>by</code> variable approach of estimating a separate smooth for each level of the factor my be quite inefficient in terms of degrees of freedom used by the model. This is especially so in situations where the estimated curves are quite similar but wiggly; why estimate many separate wiggly smooths when one, plus some simple difference smooths, will do the job just as well? In this post I look at an alternative to estimating separate smooths using an <em>ordered</em> factor for the <code>by</code> variable.
</p>
<p>
When an <em>ordered</em> factor is passed to <code>by</code>, <strong>mgcv</strong> does something quite different to the model I described previously, although the end results should be similar. What <strong>mgcv</strong> does in the <em>ordered</em> factor case is to fit <span>(L-1)</span> <em>difference smooths</em>, where <span>(l = 1, , L)</span> are the levels of the factor and <span>(L)</span> the number of levels. These smooths model the difference between the smooth estimated for the reference level and the <span>(l)</span>th level of the factor. Additionally, the <code>by</code> variable smooth doesn’t itself estimate the smoother for the reference level; so we are required to add a second smooth to the model that estimates that particular smooth.
</p>
<p>
In pseudo code our model would be something like, for ordered factor <code>of</code>,
</p>
<pre><code>model &lt;-<span> </span><span>gam</span>(y ~<span> </span>of +<span> </span><span>s</span>(x) +<span> </span><span>s</span>(x, <span>by =</span> of), <span>data =</span> df)</code></pre>
<p>
As with any <code>by</code> factor smooth we are required to include a parametric term for the factor because the individual smooths are centered for identifiability reasons. The first <code>s(x)</code> in the model is the smooth effect of <code>x</code> on the <em>reference</em> level of the ordered factor <code>of</code>. The second smoother, <code>s(x, by = of)</code> is the set of <span>(L-1)</span> <em>difference</em> smooths, which model the smooth differences between the reference level smoother and those of the individual levels (excluding the reference one).
</p>
<p>
Note that this model still estimates a separator smoother for each level of the ordered factor, it just does it in a different way. The smoother for the reference level is estimated via contribution from <code>s(x)</code> <em>only</em>, whilst the smoothers for the other levels are formed from the additive combination of <code>s(x)</code> and the relevant difference smoother from the set created by <code>s(x, by = of)</code>. This is analogous to the situation we have when estimating an ANOVA using the default contrasts and <code>lm()</code>; the intercept is then an estimate of the mean response for the reference level of the factor, and the remaining model coefficients estimate the <em>differences</em> between the mean response of the reference level and that of the other factor levels.
</p>
<p>
This <em>ordered-factor-smooth interaction</em> is most directly applicable to situations where you have a reference category and you are interested in difference between that category and the other levels. If you are interested in pair-wise comparison of smooths you could use the ordered factor approach — it may be more parsimonious than estimating separate smoothers for each level — but you will still need to post-process the results in a manner similar to that described in the previous post<a rel="nofollow" href="https://www.fromthebottomoftheheap.net/2017/12/14/difference-splines-ii/#fn1" target="_blank"><sup>1</sup></a>.
</p>
<p>
To illustrate the ordered factor difference smooths, I’ll reuse the example from the <em>Geochimica</em> <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>, and which formed the basis for the previous post.
</p>
<p>
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 hypothesized, 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>
Below I make use of 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>library</span><span>(</span><span>'dplyr'</span><span>)</span>
<span>library</span><span>(</span><span>'ggplot2'</span><span>)</span>
<span>theme_set</span><span>(</span><span>theme_bw</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
</p>
<figure>
<pre><code data-lang="r"><span>uri</span> <span>&lt;-</span> <span>'https://gist.githubusercontent.com/gavinsimpson/eb4ff24fa9924a588e6ee60dfae8746f/raw/geochimica-metals.csv'</span>
<span>metals</span> <span>&lt;-</span> <span>read_csv</span><span>(</span><span>uri</span><span>,</span> <span>skip</span> <span>=</span> <span>1</span><span>,</span> <span>col_types</span> <span>=</span> <span>c</span><span>(</span><span>'ciccd'</span><span>))</span>
<span>metals</span> <span>&lt;-</span> <span>mutate</span><span>(</span><span>metals</span><span>,</span> <span>SiteCode</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>