<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/pangaea-and-r-and-open-palaeo-data/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.fromthebottomoftheheap.net/2016/12/16/pangaea-r-open-palaeo-data/"> 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>
For a while now, I’ve been wanting to experiment with rOpenSci’s <strong>pangaear</strong> package <span data-cites="pangaear-024">(Chamberlain et al. 2016)</span>, which allows you to search, and download data from, the Pangaea, a major data repository for the earth and environmental sciences. Earlier in the year, as a member of the editorial board of <a target="_blank" rel="nofollow" href="http://www.nature.com/sdata/">Scientific Data</a>, Springer Nature’s open data journal I was handling a data descriptor submission that described a new 2,200-year foraminiferal δ<sup><sup>18</sup></sup>O record from the Gulf of Taranto in the Ionian Sea <span data-cites="Taricco2016-pv">(Carla Taricco et al. 2016)</span>. The data descriptor was recently <a target="_blank" rel="nofollow" href="http://doi.org/10.1038/sdata.2016.42">published</a> and as part of the submission Carla Taricco deposited the data set in Pangaea. So, what better opportunity to test out <strong>pangaear</strong>? (Oh and to fit a GAM to the data while I’m at it!)
</p>
<p>
The post makes use of the following packages: <strong>pangaear</strong> (obviously), <strong>mgcv</strong> and <strong>ggplot2</strong> for modelling and plotting, and <strong>tibble</strong> because <strong>pangaear</strong> returns search results and data sets in tibbles that I need to manipulate before I can fit a GAM to the δ<sup><sup>18</sup></sup>O record.
</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>&quot;pangaear&quot;</span><span>)</span><span>
</span><span>library</span><span>(</span><span>&quot;tibble&quot;</span><span>)</span><span>
</span><span>library</span><span>(</span><span>&quot;ggplot2&quot;</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>&quot;mgcv&quot;</span><span>)</span></code></pre>
</figure>
<p>
To download a data set from Pangaea you need to know the DOI of the deposit. If you don’t know the Pangaear DOI, you can search the data records held by Pangaea for specific terms. In <strong>pangaear</strong> searching is done using the <code>pg_search()</code> function. To find the data set I want, I’m going to search for records that have the string <code>“Taricco”</code> in the citation.
</p>
<figure>
<pre><code data-lang="r"><span>recs</span><span> </span><span>&lt;-</span><span> </span><span>pg_search</span><span>(</span><span>query</span><span> </span><span>=</span><span> </span><span>&quot;citation:Taricco&quot;</span><span>)</span><span>
</span><span>recs</span><span>$</span><span>citation</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">[1] &quot;Versteegh, GJM; de Leeuw, JW; Taricco, C et al. (2007): Alkenone-derived UK'37 data and sea surface temperatures (SST) of sediment cores from the Gulf of Taranto (Italy)&quot;
[2] &quot;Taricco, C; Alessio, S; Rubinetti, S et al. (2016): A foraminiferal d18O record of sediment core GT90-3 covering the last 2,200 years&quot;                                    
[3] &quot;Versteegh, GJM; de Leeuw, JW; Taricco, C et al. (2007): Alkenone-derived UK'37 data and sea surface temperatures (SST) of sediment core GT89-3&quot;                           
[4] &quot;Versteegh, GJM; de Leeuw, JW; Taricco, C et al. (2007): Alkenone-derived UK'37 data and sea surface temperatures (SST) of a combined sediment core&quot;                       
[5] &quot;Versteegh, GJM; de Leeuw, JW; Taricco, C et al. (2007): Alkenone-derived UK'37 data and sea surface temperatures (SST) of sediment core GT91-1&quot;                           
[6] &quot;Versteegh, GJM; de Leeuw, JW; Taricco, C et al. (2007): Alkenone-derived UK'37 data and sea surface temperatures (SST) of sediment core GT90-3&quot;                           </code></pre>
</figure>
<p>
Assuming that the query didn’t time out (Pangaea can be a little slow to respond on occasion, so you might find increasing the timeout on the query helps), <code>recs</code> shoud contain 6 records with <code>“Taricco”</code> in the citation. The one we want is the second entry.
</p>
<p>
To download the data object(s) associated with a record in Pangaea, we use the <code>pg_data()</code> function, supplying it with a single DOI.
</p>
<figure>
<pre><code data-lang="r"><span>res</span><span> </span><span>&lt;-</span><span> </span><span>pg_data</span><span>(</span><span>doi</span><span> </span><span>=</span><span> </span><span>recs</span><span>$</span><span>doi</span><span>[</span><span>2</span><span>])</span><span> </span><span>#</span><span> </span><span>doi</span><span> </span><span>=</span><span> </span><span>&quot;10.1594/PANGAEA.857573&quot;</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Downloading 1 datasets from 10.1594/PANGAEA.857573</code></pre>
</figure>
<figure>
<pre><code data-lang="text">Processing 1 files</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>res</span><span>
</span><span>str</span><span>(</span><span>res</span><span>[[</span><span>1</span><span>]],</span><span> </span><span>max</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">[[1]]
&lt;Pangaea data> 10.1594/PANGAEA.857573
# A tibble: 560 × 4
   `Depth [m]` `Age [a AD]` `Age [ka BP]` `G. ruber d18O [per mil PDB]`
         &lt;dbl>        &lt;dbl>         &lt;dbl>                         &lt;dbl>
1       1.4000      -188.20       2.13820                         0.742
2       1.3975      -184.33       2.13433                         0.290
3       1.3950      -180.46       2.13046                         0.706
4       1.3925      -176.59       2.12659                         0.356
5       1.3900      -172.72       2.12272                         0.558
6       1.3875      -168.85       2.11885                         0.746
7       1.3850      -164.98       2.11498                         0.346
8       1.3825      -161.11       2.11111                         0.554
9       1.3800      -157.24       2.10724                         0.510
10      1.3775      -153.37       2.10337                         0.543
# ... with 550 more rows

List of 4
 $ doi     : chr &quot;10.1594/PANGAEA.857573&quot;
 $ citation:List of 1
  ..- attr(*, &quot;class&quot;)= chr &quot;citation&quot;
 $ meta    :List of 1
  ..- attr(*, &quot;class&quot;)= chr &quot;meta&quot;
 $ data    :Classes 'tbl_df', 'tbl' and 'data.frame':   560 obs. of  4 variables:
 - attr(*, &quot;class&quot;)= chr &quot;pangaea&quot;</code></pre>
</figure>
<p>
In Pangaea, a DOI might refer to a collection of data objects, in which case the object returned by <code>pg_data()</code> would be a list with as many components as objects in the collection. In this instance there is but a single data object associated with the requested DOI, but for consistency it is returned in a list with a single component.
</p>
<p>
Rather than work with the <code>pangaea</code> object directly, for modelling or plotting it is, for the moment at least, going to be simpler if we extract out the data object, which is stored in the <code><span>\(data&lt;/code> component. We’ll also want to tidy up those variable/column names&lt;/p> &lt;figure class=&quot;highlight&quot;> &lt;pre>&lt;code class=&quot;language-r&quot; data-lang=&quot;r&quot;>&lt;span class=&quot;n&quot;>foram&lt;/span>&lt;span class=&quot;w&quot;> &lt;/span>&lt;span class=&quot;o&quot;>&amp;lt;-&lt;/span>&lt;span class=&quot;w&quot;> &lt;/span>&lt;span class=&quot;n&quot;>res&lt;/span>&lt;span class=&quot;p&quot;>[[&lt;/span>&lt;span class=&quot;m&quot;>1&lt;/span>&lt;span class=&quot;p&quot;>]]&lt;/span>&lt;span class=&quot;o&quot;>\)</span><span>data</span><span> </span><span>names</span><span>(</span><span>foram</span><span>)</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>“Depth”</span><span>,</span><span> </span><span>“Age_AD”</span><span>,</span><span> </span><span>“Age_kaBP”</span><span>,</span><span> </span><span>“d18O”</span><span>)</span><span> </span><span>foram</span></code></p>

<figure>
<pre><code data-lang="text"># A tibble: 560 × 4
    Depth  Age_AD Age_kaBP  d18O
    &lt;dbl>   &lt;dbl>    &lt;dbl> &lt;dbl>
1  1.4000 -188.20  2.13820 0.742
2  1.3975 -184.33  2.13433 0.290
3  1.3950 -180.46  2.13046 0.706
4  1.3925 -176.59  2.12659 0.356
5  1.3900 -172.72  2.12272 0.558
6  1.3875 -168.85  2.11885 0.746
7  1.3850 -164.98  2.11498 0.346
8  1.3825 -161.11  2.11111 0.554
9  1.3800 -157.24  2.10724 0.510
10 1.3775 -153.37  2.10337 0.543
# ... with 550 more rows</code></pre>
</figure>
<p>
Now that’s done, we can take a look at the data set
</p>
<figure>
<pre><code data-lang="r"><span>ylabel</span><span> </span><span>&lt;-</span><span> </span><span>expression</span><span>(</span><span>delta</span><span>^</span><span>{</span><span>18</span><span>}</span><span> </span><span>*</span><span> </span><span>O</span><span> </span><span>~</span><span> </span><span>&quot;[‰ VPDB]&quot;</span><span>)</span><span>
</span><span>xlabel</span><span> </span><span>&lt;-</span><span> </span><span>&quot;Age [ka BP]&quot;</span><span>
</span><span>ggplot</span><span>(</span><span>foram</span><span>,</span><span> </span><span>aes</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>d</span><span>18</span><span>O</span><span>,</span><span> </span><span>x</span><span> </span><span>=</span><span> </span><span>Age_kaBP</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_path</span><span>()</span><span> </span><span>+</span><span>
</span><span>scale_x_reverse</span><span>(</span><span>sec.axis</span><span> </span><span>=</span><span> </span><span>sec_axis</span><span>(</span><span> </span><span>~</span><span> </span><span>1950</span><span> </span><span>-</span><span> </span><span>(</span><span>.</span><span> </span><span>*</span><span> </span><span>1000</span><span>),</span><span> </span><span>name</span><span> </span><span>=</span><span> </span><span>&quot;Age [AD]&quot;</span><span>))</span><span> </span><span>+</span><span>
</span><span>scale_y_reverse</span><span>()</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>ylabel</span><span>,</span><span> </span><span>x</span><span> </span><span>=</span><span> </span><span>xlabel</span><span>)</span></code></pre>
</figure>
<figure>
<img alt="The δ18O record of Taricco et al (2016)" src="https://i1.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/simultaneous-intervals-revisited-plot-data-1.png?w=456" data-recalc-dims="1">
<figcaption>
The δ<sup><sup>18</sup></sup>O record of Taricco <em>et al</em> (2016)<br>
</figcaption>
</figure>
<p>
Notice that the x-axis has been reversed on this plot so that as we move from left to right the observations become younger, as is standard for a time series. In the code block above I’ve used <code>sec_axis()</code> to add an AD scale to the x-axis. This is a new feature in version 2.2.0 of <strong>ggplot2</strong> which allows secondary axis that is a one-to-one transformation of the main scale. This isn’t quite right as the two scales don’t map in a fully one-to-one fashion; as there is no year 0AD (or 0&lt;abbrv, title = &quot;Before Common Era&quot;>BCE), the scale will be a year out for the BCE period.
</p>
<p>
Note too that the y-axis has also been reversed, to match the published versions of the data. This is done in those publications because δ<sup><sup>18</sup></sup>O has an interpretation as temperature, with lower δ<sup><sup>18</sup></sup>O indicating higher temperatures. As is common for data from proxies that have a temperature interpretation, the values are plotted in a way that <em>up</em> on the plot means <em>warmer</em> and <em>down</em> means colder.
</p>
<p>
To model the data in the same time ordered way using the year BP variable we need to create a variable that is the negative of <code>Age-kaBP</code>.
</p>
<figure>
<pre><code data-lang="r"><span>foram</span><span> </span><span>&lt;-</span><span> </span><span>with</span><span>(</span><span>foram</span><span>,</span><span> </span><span>add_column</span><span>(</span><span>foram</span><span>,</span><span> </span><span>Age</span><span> </span><span>=</span><span> </span><span>-</span><span> </span><span>Age_kaBP</span><span>))</span></code></pre>
</figure>
<p>
Note that we don’t want to use the <code>Age_AD</code> scale for this as this has the problem of having a discontinuity at 0AD (which doesn’t exist).
</p>
<p>
Now we can fit a GAM to the δ<sup><sup>18</sup></sup>O record
</p>
<figure>
<pre><code data-lang="r"><span>m</span><span> </span><span>&lt;-</span><span> </span><span>gam</span><span>(</span><span>d</span><span>18</span><span>O</span><span> </span><span>~</span><span> </span><span>s</span><span>(</span><span>Age</span><span>,</span><span> </span><span>k</span><span> </span><span>=</span><span> </span><span>100</span><span>,</span><span> </span><span>bs</span><span> </span><span>=</span><span> </span><span>&quot;ad&quot;</span><span>),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>foram</span><span>,</span><span> </span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;REML&quot;</span><span>,</span><span> </span><span>select</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)</span></code></pre>
</figure>
<p>
In this instance I used an adaptive spline basis <code>bs = “ad”</code>, which allows the degree of wiggliness to vary along the fitted function. With a relatively large data set like this, which has over 500 observations, using an adaptive smoother can provide a better fit to the observations, and it is especially useful in situations where it is plausible that the response will vary more over some time periods than others. Adaptive smooths aren’t going to work well in short time series; there just isn’t the information available to estimate what in effect can be thought of as several separate splines over small chunks of the data. That said, I’ve had success with data sets with about 100–200 observations. Also note that fitting an adaptive smoother requires cranking the CPU over a lot more calculations; be aware of that if you throw a very large data set at it.
</p>
<p>
Also note that the model was fitted using REML — in most cases this is the default you want to be using as GCV can undersmooth in some circumstances. The double penalty approach of <span data-cites="Marra2011-sf">Marra and Wood (2011)</span> is used here too (<code>select = TRUE</code>), which in this instance is being used to apply a bit of shrinkage to the fitted trend; it’s good to be a little conservative at times.
</p>
<p>
The model diagnostics look OK for this model and the check of sufficient dimensionality in the basis doesn’t indicate anything to worry about (partly because we used a large basis in the first place: <code>99 = k - 1 = 100 - 1</code>)
</p>
<figure>
<pre><code data-lang="r"><span>gam.check</span><span>(</span><span>m</span><span>)</span><span>
</span><span>## RStudio users might need
## layout(matrix(1:4, ncol = 2, byrow = TRUE))
## gam.check(m)
## layout(1)
</span><span>##</span><span> </span><span>to</span><span> </span><span>see</span><span> </span><span>all</span><span> </span><span>the</span><span> </span><span>plots</span><span> </span><span>on</span><span> </span><span>one</span><span> </span><span>device</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Method: REML   Optimizer: outer newton
full convergence after 26 iterations.
Gradient range [-0.0003244869,0.000148452]
(score -128.435 &amp; scale 0.03324845).
Hessian positive definite, eigenvalue range [6.249446e-06,279.6116].
Model rank =  100 / 100 

Basis dimension (k) checking results. Low p-value (k-index&lt;1) may
indicate that k is too low, especially if edf is close to k'.

           k'    edf k-index p-value
s(Age) 99.000 15.539   0.993    0.44</code></pre>
</figure>
<figure>
<img alt="Diagnostic plots for the fitted GAM" src="https://i0.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/simultaneous-intervals-revisited-gam-check-1.png?w=456" data-recalc-dims="1">
<figcaption>
Diagnostic plots for the fitted GAM<br>
</figcaption>
</figure>
<p>
and the fitted trend is <em>inconsistent</em> with a null-model of no trend
</p>
<figure>
<pre><code data-lang="r"><span>summary</span><span>(</span><span>m</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Family: gaussian 
Link function: identity 

Formula:
d18O ~ s(Age, k = 100, bs = &quot;ad&quot;)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.455329   0.007705   59.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(Age) 15.54     99 3.357  &lt;2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.373   Deviance explained =   39%
-REML = -128.44  Scale est. = 0.033248  n = 560</code></pre>
</figure>
<p>
There is a lot of variation about the fitted trend, but a model with about 15 degrees of freedom explains about 40% of the variance in the data set, which is pretty good.
</p>
<p>
While we could use the provided <code>plot()</code> method for <code>“gam”</code> objects to draw the fitted function, I now find myself prefering plotting with <strong>ggplot2</strong>. To recreate the sort of plot that <code>plot.gam()</code> would produce, we first need to predict for a fine grid of values, here 200 values, over the observed time interval. <code>predict.gam()</code> is used to generate predictions and standard errors; the standard errors requested here use a new addition to <strong>mgcv</strong> which includes the extra uncertainty in the model because we are also estimating the smoothness parameters (the parameters that control the degree of wiggliness in the spline). This is achieved through the use of <code>unconditional = TRUE</code> in the call to <code>predict()</code>. The standard errors you get with the default, <code>unconditional = FALSE</code>, assume that the smoothness parameters, and therefore the amount of wiggliness, is known before fitting, which is rarely the case. This doesn’t make much difference in this example, but I thought I’d mention it as it is a relatively new addition to <strong>mgcv</strong>.
</p>
<figure>
<pre><code data-lang="r"><span>pred</span><span> </span><span>&lt;-</span><span> </span><span>with</span><span>(</span><span>foram</span><span>,</span><span> </span><span>data.frame</span><span>(</span><span>Age</span><span> </span><span>=</span><span> </span><span>-</span><span>seq</span><span>(</span><span>min</span><span>(</span><span>Age_kaBP</span><span>),</span><span> </span><span>max</span><span>(</span><span>Age_kaBP</span><span>),</span><span> </span><span>length</span><span> </span><span>=</span><span> </span><span>200</span><span>)))</span><span>
</span><span>pred</span><span> </span><span>&lt;-</span><span> </span><span>cbind</span><span>(</span><span>pred</span><span>,</span><span> </span><span>as.data.frame</span><span>(</span><span>predict</span><span>(</span><span>m</span><span>,</span><span> </span><span>pred</span><span>,</span><span> </span><span>se.fit</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>unconditional</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)))</span><span>
</span><span>pred</span><span> </span><span>&lt;-</span><span> </span><span>transform</span><span>(</span><span>pred</span><span>,</span><span>
</span><span>Fitted</span><span> </span><span>=</span><span> </span><span>fit</span><span>,</span><span>
</span><span>Upper</span><span> </span><span>=</span><span> </span><span>fit</span><span> </span><span>+</span><span> </span><span>(</span><span>2</span><span> </span><span>*</span><span> </span><span>se.fit</span><span>),</span><span>
</span><span>Lower</span><span> </span><span>=</span><span> </span><span>fit</span><span> </span><span>-</span><span> </span><span>(</span><span>2</span><span> </span><span>*</span><span> </span><span>se.fit</span><span>),</span><span>
</span><span>Age_kaBP</span><span> </span><span>=</span><span> </span><span>-</span><span> </span><span>Age</span><span>)</span></code></pre>
</figure>
<p>
The code above uses these standard errors to create an approximate 95% point-wise confidence on the fitted function, and prepares this in tidy format for plotting with <code>ggplot()</code>.
</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>foram</span><span>,</span><span> </span><span>aes</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>d</span><span>18</span><span>O</span><span>,</span><span> </span><span>x</span><span> </span><span>=</span><span> </span><span>Age_kaBP</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>()</span><span> </span><span>+</span><span>
</span><span>geom_ribbon</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>pred</span><span>,</span><span> </span><span>mapping</span><span> </span><span>=</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>Age_kaBP</span><span>,</span><span> </span><span>ymin</span><span> </span><span>=</span><span> </span><span>Lower</span><span>,</span><span> </span><span>ymax</span><span> </span><span>=</span><span> </span><span>Upper</span><span>),</span><span>
</span><span>fill</span><span> </span><span>=</span><span> </span><span>&quot;grey&quot;</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>NA</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.7</span><span>,</span><span> </span><span>inherit.aes</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>)</span><span> </span><span>+</span><span>
</span><span>geom_path</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>pred</span><span>,</span><span> </span><span>mapping</span><span> </span><span>=</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>Age_kaBP</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>Fitted</span><span>),</span><span> </span><span>inherit.aes</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>,</span><span>
</span><span>size</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span> </span><span>+</span><span>
</span><span>scale_x_reverse</span><span>(</span><span>sec.axis</span><span> </span><span>=</span><span> </span><span>sec_axis</span><span>(</span><span> </span><span>~</span><span> </span><span>1950</span><span> </span><span>-</span><span> </span><span>(</span><span>.</span><span> </span><span>*</span><span> </span><span>1000</span><span>),</span><span> </span><span>name</span><span> </span><span>=</span><span> </span><span>&quot;Age [AD]&quot;</span><span>))</span><span> </span><span>+</span><span>
</span><span>scale_y_reverse</span><span>()</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>ylabel</span><span>,</span><span> </span><span>x</span><span> </span><span>=</span><span> </span><span>xlabel</span><span>)</span></code></pre>
</figure>
<figure>
<img alt="Observed δ18O values with the fitted trend and 95% point-wise confidence interval superimposed" src="https://i1.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/simultaneous-intervals-revisited-final-plot-1.png?w=456" data-recalc-dims="1">
<figcaption>
Observed δ<sup><sup>18</sup></sup>O values with the fitted trend and 95% point-wise confidence interval superimposed<br>
</figcaption>
</figure>
<p>
<span data-cites="Taricco2009-eh">C Taricco et al. (2009)</span> used singular spectrum analysis (SSA), among other spectral methods, to decompose the δ<sup><sup>18</sup></sup>O time series into components of variability with a range of periodicities. A visual comparison with the SSA components and the fitted GAM trend, suggests that the GAM trend maps on to the sum of the long-term trend component plus the ~600 year and (potentially) the 350 year frequency components of the SSA. This does make we wonder a little about how real the higher frequency components identified in the SSA are? No matter how hard I tried (even setting the basis dimension of the GAM to <code>k = 500</code>) I couldn’t get it to be more wiggly than shown in the plots above). Figure 4 of <span data-cites="Taricco2009-eh">C Taricco et al. (2009)</span> also showed the spectral power for the 4 non-trend components from the SSA. The power associated with the 200-year and the 125-year components is substantially less than that of the two longer-frequency components. The significance of the SSA components was determined using a Monte Carlo approach <span data-cites="Allen1996-gc">(Allen and Smith 1996)</span>, where surrogate time series are generate using AR(1) noise. It’s reasonable to ask whether this is a reasonable null model for these data? It’s also reasonable to ask whether the GAM approach I used above has sufficient statistical power to detect higher-freqency components if they actually exist? This warrants further study.
</p>
<p>
I started this post with some details on why I was prompted to look at this particular data set. Palaeo-scientists have a long record of sharing data — less so in some specific fields: yes, I’m looking at you (&amp; me), palaeolimnologists — but, and perhaps this is just me, I’m seeing more of an open-data culture within palaeoecology and palaeoclimatology. This is great to see, and avenues for publishing and hence generating traditional academic merit for the data we generate will only help foster this. With my “editorial board member” hat on, I would encourage people to consider writing a data paper and submitting it to Scientific Data or one of the other data journals that are springing up. But, if you can’t or don’t want to do that, depositing your data in an open repository like Pangaea brings with it many benefits and is something that we the palaeo community should be supportive of. I wouldn’t have been writing this post if Tarrico and co-authors hadn’t chosen to make their data openly available.
</p>
<p>
And that brings me on to my final point for this post; having access to an excellent data repository like Pangaea from within a data analysis platform like R makes it so much easier to engaged with the literature and ask new and interesting questions. I’ve highlighted Pangaea here, but other initiatives are doing a great job of making palaeo data available and also deserve our recognition and support, like the <a target="_blank" rel="nofollow" href="http://www.neotomadb.org/">Neotoma</a> database; we might take access to these reources for granted, but implementing and maintaining web servers and APIs requires a lot of time, effort, and resources. Also, this post wouldn’t have been possible without the work of the wonderful <a target="_blank" rel="nofollow" href="https://ropensci.org/">rOpenSci</a> community that make available R packages to query the APIs of online repositories like Pangaea and Neotoma. Thank you!
</p>
<h3>
References<br>
</h3>
<div>
<div>
<p>
Allen, Myles R, and Leonard A Smith. 1996. “Monte Carlo SSA: Detecting Irregular Oscillations in the Presence of Colored Noise.” <em>Journal of Climate</em> 9 (12): 3373–3404. doi:<a target="_blank" rel="nofollow" href="https://doi.org/10.1175/1520-0442(1996)009%3c3373:MCSDIO%3e2.0.CO;2">10.1175/1520-0442(1996)009&lt;3373:MCSDIO>2.0.CO;2</a>.
</p>
</div>
<div>
<p>
Chamberlain, Scott, Kara Woo, Andrew MacDonald, Naupaka Zimmerman, and Gavin Simpson. 2016. <em>Pangaear: Client for the ’Pangaea’ Database</em>. <a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=pangaear">https://CRAN.R-project.org/package=pangaear</a>.
</p>
</div>
<div>
<p>
Marra, Giampiero, and Simon N Wood. 2011. “Practical Variable Selection for Generalized Additive Models.” <em>Computational Statistics &amp; Data Analysis</em> 55 (7): 2372–87. doi:<a target="_blank" rel="nofollow" href="https://doi.org/10.1016/j.csda.2011.02.004">10.1016/j.csda.2011.02.004</a>.
</p>
</div>
<div>
<p>
Taricco, C, M Ghil, S Alessio, and G Vivaldo. 2009. “Two Millennia of Climate Variability in the Central Mediterranean.” <em>Climate of the Past</em> 5 (2): 171–81. doi:<a target="_blank" rel="nofollow" href="https://doi.org/10.5194/cp-5-171-2009">10.5194/cp-5-171-2009</a>.
</p>
</div>
<div>
<p>
Taricco, Carla, Silvia Alessio, Sara Rubinetti, Gianna Vivaldo, and Salvatore Mancuso. 2016. “A Foraminiferal <span>()</span>18O Record Covering the Last 2,200 Years.” <em>Scientific Data</em> 3 (June). Nature Publishing Group: 160042. doi:<a target="_blank" rel="nofollow" href="https://doi.org/10.1038/sdata.2016.42">10.1038/sdata.2016.42</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/2016/12/16/pangaea-r-open-palaeo-data/"> 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/GHhqtXJC1kU"></a></div>