knitr::opts_chunk$set(echo = TRUE, eval = FALSE, message = FALSE)
Transfer functions for palaeoenvironmental reconstructions can be fitted with the rioja
or analogue
packages in R. Other packages such as vegan
and palaeoSig
can help interpretation.
library("rioja") library("palaeoSig") library("tidyverse") library("ggpalaeo")
Packages rioja
and analogue
have functions with the same names - don't load both at the same time.
set.seed(42)#ensures results are same every time
For convenience, we are going to use some of the datasets available within the rioja
package.
#NB any text after a '#' is a comment and ignored by R #get help on any function with '?function_name' data("SWAP")#SWAP calibration set data("RLGH")#fossil data from the Round Loch of Glenhead
This has imported the data in lists. We can break these up into their component parts to simplify the code later.
# '<-' assigns the result on the right to the name on the left # '$' extracts a named component of a list spp <- SWAP$spec pH <- SWAP$pH spNames <- SWAP$names fos <- RLGH$spec ageDepth <- RLGH$depths
Plot calibration set species abundances against pH. Add a smoother to highlight the main trends.
ggplot(spp/100, aes(x = pH, y = TA003A)) + geom_point() + geom_smooth(method = "gam", method.args = list(family = "quasibinomial"))
The relationship between all species and the environment can be seen with an ordination.
ord <- cca(sqrt(spp) ~ pH) ord autoplot(ord) #learn more about ordinations in Bio303
mod <- WA(y = sqrt(spp), x = pH) #crossvalidate mod <- crossval(mod) performance(mod) autoplot(mod)
mod <- WA(y = sqrt(spp), x = pH) #crossvalidate mod <- crossval(mod)
performance(mod)
Question: do any of the other transfer function methods perform better - try WAPLS, MAT, MLRC? Look only at the cross-validated performance.
By default crossval
uses leave-one-out cross-validation. Bootstrap cross-validation is also useful as it gives sample specific errors. Other forms of cross-validation are appropriate if there is pseudoreplication or autocorrelation.
strat.plot(d = fos, yvar = ageDepth$Age, scale.percent = TRUE, y.rev = TRUE)
recon <- predict(mod, sqrt(fos)) tibble(age = ageDepth$Age, pH = recon$fit[, "WA.inv"]) %>% ggplot(aes(x = age, y = pH)) + geom_line()
Several diagnostics available - goodness of fit, distance to closest analogue - in analogue
package.
Timetrack of fossil samples in calibration set ordination
#use 'analogue::' to access functions without loading library - prevents name conflicts with rioja tt <- analogue::timetrack(X = sqrt(spp), passive = sqrt(fos), env = pH) plot(tt)
Residual length to constraining ordination axis
rlen <- analogue::residLen(X = sqrt(spp), env = pH, passive = sqrt(fos)) plot(rlen) autoplot(rlen, df = ageDepth, x_axis = "Age") + labs(x = "Age", y = "Squared residual distance", fill = "Goodness of fit")
Colours show how good the fit is.
Analogue distances
AD <- analogue_distances(spp, fos) autoplot(AD, df = ageDepth, x_axis = "Age") + labs(y = "Squared-chord distance", x = "Age")
Significance tests in palaeoSig
.
tst1 <- randomTF(spp = sqrt(spp), env = pH, fos = sqrt(fos), n = 99, col = 1, fun = WA) tst1$sig autoplot(tst1)
Juggins and Birks (2012) Quantitative Environmental Reconstructions from Biological Data. In Tracking Environmental Change Using Lake Sediments, Data Handling and Numerical Techniques https://www.springer.com/gp/book/9789400727441
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.