shiggo/jezioro
Tools for Paleolimnological Analyses at the Paleoecological Environmental Assessment and Research Laboratory

Package index
Search the shiggo/jezioro package
Vignettes
Functions
30
Source code
20
Man pages
20
Home
/
GitHub
/
shiggo/jezioro
/

In shiggo/jezioro: Tools for Paleolimnological Analyses at the Paleoecological Environmental Assessment and Research Laboratory 

library(jezioro)

Introduction

The jezioro package is a collection of functions and datasets intended for use by members of the Paleoecological Environmental Assessment and Research Laboratory (PEARL). The package is maintained by Adam Jeziorski, and contains contributions from Adam Jeziorski, Joshua Thienpont, and Andrew Labaj. If you have any suggestions for its improvement please let us know.

The package includes the following functions (along with some example data sets to illustrate their use):

The package also contains some commonly used calibration sets to aid the construction of the relevant transfer functions for application to fossil assemblages:

The chirInfo data set tabulates several aspects of chironomid taxonomy/ecology to facilitate manipulation of chironomid data (e.g. sorting by littoral/profundal habitat preferences)

Functions

cladCount

cladCount determines the maximum number of individuals represented in raw counts of cladoceran subfossils.

In paleolimnological cladoceran analyses, all remains (carapaces, headshields, ephippia, postabdominal claws, etc.) should be tabulated seperately, with only the most frequently encountered remain for each taxon used to estimate its abundance (Korhola and Rautio 2001). cladCount is a function to quickly calculate the maximum number of individuals represented in appropriately formatted raw counts of subfossils.

Input data must contain the taxon name and subfossil name in the first two columns (respectively), with each subsequent column a sample/interval. If a taxon can be represented by more than one subfossil type, the 'Taxon' cell should be left blank from the second row onwards (these blank cells are how the function identifies the number of subfossils present for each taxon).

The required format of cladCount input data is illustrated by cladCountInput:

data(cladCountInput)

data(cladCountInput)
knitr::kable(head(cladCountInput[,1:6], 10), row.names=F)

The largest number of individuals possibly represented by all the subfossils for a particular taxon is identified by attributing the following numbers of each subfossil to an individual (rounding up in the case of 'half-individuals'; Korhola and Rautio, 2001):

The output of cladCount is controlled by three arguments (percCutoff, sampleCutoff, and outputType):

For example, to return the number of individuals for all taxa present in the dataset:

data(cladCountInput)
cladCount(cladCountInput)

Similarly, to return the relative abundances of only those taxa with greater than 4% abundance in at least 2 samples:

data(cladCountInput)
cladCount(cladCountInput, percCutoff=4, sampleCutoff=2, outputType='gt')


References

Korhola A, Rautio M (2001) 2. Cladocera and other branchiopod crustaceans. In: Smol JP, Birks HJB, Last WM (eds.) Tracking Environmental Change Using Lake Sediments. Volume 4: Zoological Indicators. Kluwer Academic Publishers, Dordrecht, The Netherlands, pp 4-41

dipteranVWHO

dipteranVWHO harmonizes chironomid/chaoborid count data with the taxonomy used in the Quinlan and Smol (2001, 2010) calibration set, returns the chironomid inferred volume weighted hypolimnetic oxygen concentrations (VWHO) from the sample data using the inference model described in Quinlan and Smol (2001, 2010), and also performs some analog matching between the sample data and the calibration set.

The required format of the input data is provided in an example data set dipteranVWHOInput

The output of dipteranVWHO is modified by three arguments (evaluate, percentileCut, and lowCount):

The function returns a list of four lists:

In addition, three plots are produced:

For example, to perform the analyses on the example dataset using a modern analog cutoff of 10, and only flagging counts with less than 40 individuals as low:

data(dipteranVWHOInput)
dipteranVWHO(dipteranVWHOInput, evaluate=TRUE, percentileCut=10, lowCount=40)


References

Quinlan R, Smol JP (2001) Setting minimum head capsule abundance and taxa deletion criteria in chironomid-based inference models. Journal of Paleolimnology 26: 327-342

Quinlan R, Smol JP (2001) Chironomid-based inference models for estimating end-of-summer hypolimnetic oxygen from south-central Ontario shield lakes. Freshwater Biology 46: 1529-1551

Quinlan R, Smol JP (2010) Use of Chaoborus subfossil mandibles in models for inferring past hypolimnetic oxygen. Journal of Paleolimnology 44: 43-50

interpDates

interpDates interpolates dates for any undated intervals in a sediment core from the midpoint and age of the dated intervals, as well as the sectioning resolution.

Input data must contain two columns, and the argument intervalWidth is used to specifiy the sectioning resolution:

The required format of interpDates input data:

interpDatesInput <- cbind(c(0.25, 4.25, 8.25, 18.25), c(2017, 2000, 1950, 1850))
colnames(interpDatesInput) <- c("Midpt", "Age")
interpDatesInput

The output of interpDates is modified by a single argument (intervalWidth):

Three different interpolation methods are used to determine the dates:

For example, to return the interpolated dates from all 3 methods (as well as a plot comparing the two fitted lines) for a core with four dated intervals and a constant sectioning resolution of 0.5cm:

interpDates.input <- cbind(c(0.25, 4.25, 8.25, 18.25), c(2017, 2000, 1950, 1850))
interpDates(interpDates.input)

interpDates.input <- cbind(c(0.25, 4.25, 8.25, 18.25), c(2017, 2000, 1950, 1850))
temp <- interpDates(interpDates.input)

Similarly, to return interpolated dates for a core with four dated intervals sectioned at a 0.5cm resolution for the first 10 intervals, and then a 1.0cm resolution for the next 10 intervals:

interpDates.input <- cbind(c(0.25, 4.25, 8.5, 13.5), c(2017, 2000, 1950, 1850))
intervalWidth <- c(rep(0.5, 10), rep(1.0, 10))
interpDates(interpDates.input, intervalWidth)

For more involved approaches to the estimation of age-depth relationships consult Blaauw and Heegaard (2012).


References

Blaauw M, Heegaard E (2012) 12. Estimation of age-depth relationships. In: Birks HJB, Lotter AF, Juggins S, Smol JP (eds.) Tracking Environmental Change Using Lake Sediments. Volume 5: Data Handling and Numerical Techniques. Springer, Netherlands, pp 379-413

quickGAM

quickGAM fits a GAM along with the related analyses described in Simpson 2018. This function applies the code provided in the Simpson 2018 Supplementary Information.

Input data must be composed of two vectors:

The output of 'quickGAM' is modified by three optional arguments, that can be used to ensure the axis labels and y-axis range remain consistent across all the output plots.

quickGAM returns a series of six plots: - A simple scatterplot of the input data - Estimated CAR1 process from the GAM fitted to the input data ( i.e. an autocorrelation check) - GAM-based trends fitted to the input data - Estimated trends (w 20 random draws of the GAM fit to the input data) - 95% simultaneous confidence intervals (light grey) and across-the-function confidence intervals (dark grey) on the estimated trends. - Estimated first derivatives (black lines) and 95% simultaneous confidence intervals of the GAM trends fitted to the data. Where the simultaneous interval does not include 0, the models detect significant temporal change in the response.

For example, to fit a GAM between sediment depth (x) and the inferred-chlorophyll a values (y) provided in the example data set vrsChlaInput:

Read in the example data set:

data(vrsChlaInput)

The optional arguments allow for consistent output plots.

indepVarLabel <- "Core Depth (cm)"
depenVarLabel <- expression("VRS-Inferred Chlorophyll "~italic("a")~" (mg"%.%"g"^"-1"*textstyle(")"))
depenVarRange <- c(0, 0.0825)

Note, that to use the default values for the plots, using 'drop=FALSE' when reading in data can preserve column names for use in the plots.

quickGAM(plotData[,1, drop=FALSE],plotData[,2, drop=FALSE], xLabel = indepVarLabel, yLabel = depenVarLabel, yRange = depenVarRange)


References

Simpson GL (2018) Modelling Palaeoecological Time Series Using Generalised Additive Models. Frontiers in Ecology and Evolution [https://dx.doi.org/10.3389/fevo.2018.00149]

vrsChla

vrsChla infers chlorophyll a concentrations of sediments from spectral measurements of absorbance at wavelengths between 650-700 nm, following the approach described in Wolfe et al. (2006) and Michelutti et al. (2010), and reviewed in Michelutti and Smol (2016).

The technique uses a simple linear predictive model to infer sedimentary chlorophyll a concentrations (along with its primary degradation products, pheophytin a and pheophorbide a) from the absorbance peak centered on 675nm following the equation:

$$[\mbox{chlorophyll }\textit{a } + \mbox{ derivatives}] = 0.0919 · \mbox{peak area}_{\textit{ 650-700 nm}} + 0.0011$$

Typical lake sediment spectrum showing the distinct peak centred near 675nm (modified from Michelutti et al. 2010)

Input data must contain 27 columns:

The required format of vrsChla input data is illustrated by vrsChlaInput:

data(vrsChlaInput)

data(vrsChlaInput)
head(vrsChlaInput[,1:6], 10)

NOTE: When using the Model 6500 series Rapid Content Analyzer at PEARL, the necessary values are contained in the 'spectra' tab of the excel file output (although they must be transposed). Ensure cells are formatted to 15 decimal places to avoid small rounding errors.

For example, to determine the chlorophyll a concentrations of vrsChlaInput and plot sediment depth vs chlorophyll a (red line denotes the estimated 0.01 mg·g^-1^ lower detection limit of the method):

data(vrsChlaInput)
vrsChla(vrsChlaInput)


References

Michelutti N, Blais JM, Cumming BF, Paterson AM, Ruhland K, Wolfe AP, Smol JP (2010) Do spectrally inferred determinations of chlorophyll a reflect trends in lake trophic status? Journal of Paleolimnology 43: 208-217

Michelutti N, Smol JP (2016) Visible spectroscopy reliably tracks trends in paleo production. Journal of Paleolimnology 56: 253-265

Wolfe AP, Vinebrooke RD, Michelutti N, Rivard B, Das B (2006) Experimental calibration of lake-sediment spectral reflectance to chlorophyll a concentrations: methodology and paleolimnological validation. Journal of Paleolimnology 36: 91-100

the binford series

The binford series are three functions to calculate ^210^Pb dates via gamma spectroscopy using the gamma counters at PEARL.

The functions are intended to be ran in sequence with the output of binfordRho and binfordActivity providing the input for binfordDates. In practice, binfordRho and binfordDates each be run once, while binfordActivity will be run repeatedly as individual samples pass through the gamma counter, until background activities are obtained.

Flow chart illustrating the intended use of the binford series.

NOTE: These functions were rendered obsolete by the lab's adoption of ScienTissiME. The principal use of these functions now is to rexamine cores originally dated prior to 2013.

binfordRho

binfordRho determines the water content and dry sediment density for the intervals of a sediment core freeze-dried in preparation for gamma dating, from the dry masses of the freeze-dried intervals and the bag weights of wet sediment. The values for intervening (i.e. non-prepared samples) are interpolated.

Input data must contain seven columns:

The required format of binfordRho input data is illustrated by binfordRhoInput:

data(binfordRhoInput)

data(binfordRhoInput)
head(binfordRhoInput, 10)

NOTE: Columns 1-3 must have values for the entire core length. Columns 4-7 will only have values for the intervals prepared for dating by freeze-drying, with the values for non-prepared intervals left blank. If the entire "Whirlpak" bag was freeze-dried (i.e. subsample=FALSE), then Columns 4-5 should be left blank and the mass (g) of the bag of sediment after freeze drying entered into Column 6.

The output of binfordRho is controlled by three arguments (bagwt, subsample, and coreD):

For example, to calculate water content and sediment densities for the example dataset, assuming a larger "Whirlpak" bag was used along with core tubes 6.8cm in diameter.

data(binfordRhoInput)
binfordRho(binfordRhoInput, bagwt=5, coreD=6.8)

data(binfordRhoInput)
head(binfordRho(binfordRhoInput, bagwt=5, coreD=6.8), 10)

The output table from binfordRho is intended to be used as input for binfordDates and has seven columns:

binfordActivity

binfordActivity calculates the activity of the ^210^Pb, ^137^Cs, and ^214^Bi isotopes, corrected for the efficiency of the gamma counter device, and the time between sediment core collection and analysis, using the approaches described in Schelske et al. (1994).

Input data must contain 14 columns:

The required format of binfordActivity input data is illustrated by binfordActivityInput:

data(binfordActivityInput)

data(binfordActivityInput)
head(binfordActivityInput, 10)

NOTE: Columns 1 and 2 must have values for all samples prepared for gamma analysis (the same values entered into Columns 1 and 2 of binfordRho), but no values for columns 3-9 in cases where samples have not yet been run or are not required (in cases where background has been reached).

The output of binfordActivity is controlled by the following arguments:

NOTE: It is very important to use the correct background and error measurements (e.g. PbBk, PbBkErr, etc.). That is, the values specific to both the gamma counter used and the time period when the samples were measured.

binfordActivity is intended to be run in an iterative fashion, following analysis of each sample. Allowing accurate and timely determination of when supported levels of ^210^Pb have been reached, in order to avoid the superfluous analysis of samples beyond background.

For example, to calculate activities for the example dataset, using the default error values (i.e. 0), label the plot as "Example Lake", and use Bq/g for y-axis units:

data(binfordActivityInput)
binfordActivity(binfordActivityInput, LakeName="Example Lake", graph.units="dpm")

data(binfordActivityInput)
head(binfordActivity(binfordActivityInput, LakeName="Example Lake", graph.units="bq"), 10)

The output table from binfordActivity is intended to be used as input for binfordDates and has eight columns:


References

Schelske CL, Peplow A, Brenner M, Spencer CN (1994) Low-background gamma counting: applications for ^210^Pb dating of sediments Journal of Paleolimnology 10: 115-128

binfordDates

binfordDates uses the unsupported fraction of ^210^Pb activity in sediment samples to calculate sediment ages via the Constant Rate of Supply (CRS) model following the procedures outlined in Binford (1990).

binfordDates uses the output of binfordRho and binfordActivity as its input data.

The ouput of binfordDates is controlled by six arguments:

For example, to calculate sediment ages for the example datasets using only default values.

data(binfordRhoInput)
data(binfordActivityInput)
binfordDates(binfordRho(binfordRhoInput), binfordActivity(binfordActivityInput))

data(binfordRhoInput)
data(binfordActivityInput)
binfordRhoOutput <- binfordRho(binfordRhoInput)
binfordActivityOutput <- binfordActivity(binfordActivityInput)

binfordDates(binfordRhoOutput, binfordActivityOutput)

The output table from binfordDates has 15 columns:

NOTE: The outpuy from binfordDates lends itself to use with interpDates:

data(binfordRhoInput)
data(binfordActivityInput)
dates <- binfordDates(binfordRho(binfordRhoInput), binfordActivity(binfordActivityInput))
dates <- dates[,c(3,13)]
interpDates(dates, 0.5)


References

Binford MW (1990) Calculation and uncertainty analysis of 210PB dates for PIRLA project lake sediment cores. Journal of Paleolimnology 3: 253-267

wiltseBrien/Decompose

The Wiltse functions wiltseBrien and wiltseDecompose test for homogeneity and coherence within a dataset. They are slightly modified (so that they no longer write to the global environment) versions of the 'brien' and 'decompose' functions described in Brendan Wiltse’s 2014 PhD thesis.

wiltseBrien

wiltseBrien performs the Brien test (Brien et al. 1984) on variables contained within a data frame, calculating the grand mean, main effects interactions and equal correlation for a correlation matrix of the z-scored data and returns the degrees of freedom, chi-squared value and p-value for each test.

The input data should be a data frame with the variables to be tested arranged in columns with appropriate column names. An example data set is provided by `wiltseInput' that contains dates in column 1, and the rest of the columns are PCA axis-1 scores for the dominant diatom taxa in 8 lakes in the ELA. A plot of the data shows the trends through time.

data(wiltseInput)
matplot(wiltseInput$Date, wiltseInput[,-match("Date", colnames(wiltseInput))], ylim=c(-4,4),bty="n", type="l", xlab="Year", ylab="PCA Axis 1 Score")

An initial run of the Brien test (on just the PCA Axis 1 scores) with wiltseBrien:

wiltseBrien(wiltseInput[,-match("Date", colnames(wiltseInput))])



wiltseDecompose

wiltseDecompose builds upon the Brien's test by identifying homogenous and coherent subsets within the correlation matrix.

An initial test for homogeneity and coherence within the correlation matrix is performed. If heterogeneity is detected (a homogenous matrix fails to reject the null-hypothesis for the tests of main effects, interactions and equal correlations), then the lake with lowest mean correlation within the matrix is deleted and the test is rerun on the reduced matrix. This process is repeated until the remaining variables are both homogenous and synchronous (i.e. fail to reject the null-hypothesis for the grand mean).

To run wiltseDecompose on the wiltseInput data:

wiltseDecompose(wiltseInput[,-match("Date", colnames(wiltseInput))], print.detail=T)

Through the decomposition process, first Lake 224 was removed, then Lake 468, then Lake 129 and finally Lake 199:

par(mfrow=c(2,2),mar=c(5,5,1,2), oma=c(0.1,0.1,0.1,0.1))
matplot(wiltseInput$Date, wiltseInput[,-match(c("Date", "ELA224"), colnames(wiltseInput))], ylim=c(-4,4),bty="n", type="l", xlab="Year", ylab="PCA Axis 1 Score", main="Lake 224 Removed")

matplot(wiltseInput$Date, wiltseInput[,-match(c("Date", "ELA224", "ELA468"), colnames(wiltseInput))], ylim=c(-4,4),bty="n", type="l", xlab="Year", ylab="PCA Axis 1 Score", main="Lake 468 Removed")

matplot(wiltseInput$Date, wiltseInput[,-match(c("Date", "ELA224", "ELA468", "ELA129"), colnames(wiltseInput))], ylim=c(-4,4),bty="n", type="l", xlab="Year", ylab="PCA Axis 1 Score", main="Lake 129 Removed")

matplot(wiltseInput$Date, wiltseInput[,-match(c("Date", "ELA224", "ELA468", "ELA129", "ELA99"), colnames(wiltseInput))], ylim=c(-4,4),bty="n", type="l", xlab="Year", ylab="PCA Axis 1 Score", main="Lake 99 Removed")

To test for synchrony among just the removed lakes (224, 468, 129 and 99):

wiltseBrien(wiltseInput[,match(c("ELA224", "ELA468", "ELA129", "ELA99"), colnames(wiltseInput))])
matplot(wiltseInput$Date, wiltseInput[,match(c("ELA224", "ELA468", "ELA129", "ELA99"), colnames(wiltseInput))], ylim=c(-4,4),bty="n", type="l", xlab="Year", ylab="PCA Axis 1 Score")

Even though the matrix was determined to be both synchronous and homogenous, Lake 99 looks like it might be doing something different. Therefore to run the test without Lake 99:

wiltseBrien(wiltseInput[,match(c("ELA224", "ELA468", "ELA129"), colnames(wiltseInput))])
matplot(wiltseInput$Date, wiltseInput[,match(c("ELA224", "ELA468", "ELA129"), colnames(wiltseInput))], ylim=c(-4,4),bty="n", type="l", xlab="Year", ylab="PCA Axis 1 Score")


References

Brien CJ, Venables WN, James AT, Mayo O (1984) An analysis of correlation matrices: equal correlations. Biometrika 71: 545-54

Rusak JA, Yan ND, Somers KM, McQueen DJ (1999) The temporal coherence of zooplankton population abundances in neighboring north-temperate lakes. The American Naturalist 153: 46–58

Wiltse B (2014) The response of Discostella species to climate change at the Experimental Lakes Area, Canada. PhD Thesis Queen's University

Datasets

tpHall1996

tpHall1996 contains the total phosphorus (TP) measurements and sedimentary diatom assemblage data for the 54 lake south-central Ontario calibration set from Hall and Smol (1996). The purpose of this dataset is to facilitate reproducing the described diatom-inferred TP model for application to down-core diatom sedimentary assemblage data.

The 425x57 data frame is composed of the following:

NOTE: The genus/species names included in the data frame are those used in Hall and Smol (1996). Due to changes in taxonomic nomenclature since its publication, applying the model to more recent data will likely require some taxonomic harmonization.


References

Hall RI, Smol JP (1996) Paleolimnological assessment of long-term water-quality changes in south-central Ontario lakes affected by cottage development and acidification. Canadian Journal of Fisheries and Aquatic Sciences 53: 1-17

vwhoQuinlan2010

vwhoQuinlan2010 contains the the volume-weighted hypolimnetic oxygen (VWHO) concentrations and sedimentary midge (chironomid and chaoborid) assemblage data for the 54 lake south-central Ontario calibration set from Quinlan and Smol (2001, 2010). The purpose of this dataset is to facilitate reproducing the described midge-inferred VWHO model for application to down-core midge sedimentary assemblage data.

The 54x47 data frame is composed of the following:

Information on each lake (i.e. latitude, longitude, etc.) and the 44 taxa included in the model are provided in Quinlan and Smol (2001, 2010).

NOTE: The taxonomic information included in the data frame are those used in Quinlan and Smol (2001,2010). Due to changes in taxonomic nomenclature since its publication, applying the model to more recent data will likely require some taxonomic harmonization.


References

Quinlan R, Smol JP (2001) Chironomid-based inference models for estimating end-of-summer hypolimnetic oxygen from south-central Ontario shield lakes. Freshwater Biology 46: 1529-1551

Quinlan R, Smol JP (2010) Use of Chaoborus subfossil mandibles in models for inferring past hypolimnetic oxygen. Journal of Paleolimnology 44: 43-50

chirInfo

chirInfo is a data frame tabulating several aspects of chironomid ecology into a single file to facilitate sorting/manipulating chironomid data files. The data frame is currently organized as follows:


References

Brooks SJ, Langdon PG, Heiri O (2007) The identification and use of Palaearctic Chironomidae larvae in paleoecology. Quat Res (London): Technical Guide no. 10.

Coffman WP (1978) Chironomidae. pp. 345-376. In: Merrit, RW, Cummins KW (Eds). An introduction of aquatic insects of North America. Kendall Hunt Publishing, Dubuque, 441p

Quinlan R, Smol JP (2001) Chironomid-based inference models for estimating end-of-summer hypolimnetic oxygen from south-central Ontario shield lakes. Freshwater Biology 46: 1529-1551

Quinlan R, Smol JP (2010) Use of Chaoborus subfossil mandibles in models for inferring past hypolimnetic oxygen. Journal of Paleolimnology 44: 43-50

Version History

v0.20 (Sept 2020)

General

Additions

v0.19 (Sept 2019)

General

Additions

Modifications

v0.18 (Dec 2018)

Reorganization

Additions

v0.17 (Dec 2017)

Reorganization

Additions

Modifications

v0.16 (May 2017)

Additions

Modifications

v0.14 (Nov 2015)

Modifications



shiggo/jezioro documentation built on Sept. 7, 2020, 7:34 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close