Package 'analogue'

Description

analogue is a package for quantitative palaeoecology with a focus on analogue methods, transfer functions, and data handling and display.

Analogue methods

analogue provides functions for analogue matching and the modern analogue technique (MAT) via analog and mat. A wide range of dissimilarity coefficients are available via the distance function.

Additional analysis of modern and no-analogue problems is facilitated via a range of functions implementing many methods from the literature. In particular, the receiver operating characteric (ROC) curves method of Gavin et al (2003) is available in roc and a related method employing direct logistic regression modelling (Simpson & Birks, 2012) is available in logitreg.

Transfer function methods

Several approaches to fitting transfer function models are provided by analogue:

wa:

Simple and tolerance-downweighted weighted averaging with classical, inverse, and monotonic spline deshrinking.

mat:

The modern analogue technique (MAT).

pcr:

Principal components regression with ecologically meaningful transformations

A range of functions for working with and exploring training sets and palaeoenvironmental reconstructions is also included in analogue. These include

crossval

leave-one-out, repeated k-fold, and bootstrap cross-validation methods.

compare

compare properties of taxa or other proxies across modern and fossil data sets.

evenSample

are training set samples evenly distributed along the gradient of interest?

splitSample

splits a gradient into a set of bins or chunks and samples evenly from within each chunk to create a representative test set for cross-validation.

timetrack

overlays a fossil or secondary data set on to an (constrained) ordination of a modern or reference data set.

weightedCor

implements the weighted correlation test of a Weighted Averaging reconstruction as proposed by Telford & Birks (2011).

Utilities

analogue provides a range of utilities for working with palaeo data.

tran

a range of transformations applicable to or commonly used with palaeo data.

Stratiplot

draws stratigraphic diagrams using the Lattice package.

join

merging of modern/training set and fossil data sets.

chooseTaxa

selects taxa that meet certain abundance and occurrence criteria.

Documentation

A full tutorial and worked example for the main features of analogue matching and MAT is avilable in the vignette

analogue_methods

Analogue Methods in Palaeoecology

Author(s)

Gavin L. Simpson, Jari Oksanen

Maintainer: Gavin L. Simpson <[email protected]>

References

Gavin, D.G., Oswald, W.W., Wahl, E.R. and Williams, J.W. (2003) A statistical approach to evaluating distance metrics and analog assignments for pollen records. Quaternary Research 60, 356–367.

Simpson, G.L. & Birks H.J.B. (2012) Statistical Learning in Palaeolimnology. In Birks, H.J.B, Lotter, A.F. Juggins S., and Smol, J.P. (Eds) Tracking Environmental Change Using Lake Sediments, Volume 5: Data Handling and Numerical Techniques. Springer, Dordrecht.

Telford R.J. and Birks, H.J.B. (2011) A novel method for assessing the statistical significance of quantitative reconstructions inferred from biotic assemblages. Quanternary Science Reviews 30:1272-1278.

Description

The classic pollen data set from Abernethy Forest in the Scottish highlands, UK. The data originate from the work of Hilary Birks and Rolf Mathewes (1978) and have been analysed in several texts on quantitative numerical palaeoecology.

The data set consists of 36 pollen taxa from 49 levels, with two additional variables; Age, the age of each sample, and Depth the depth (in cm) below the surface of the peat sequence from which the core was taken.

Usage

data(abernethy)

Format

abernethy is a data frame with 49 samples on 36 species plus sample Age and Depth (in cm).

Source

These data were provided in electronic format by Prof. H. John B. Birks. The original source is Birks and Mathewes (1978).

References

Birks, H.H. and Mathewes, R.W. (1978) Studies in the vegetational history of Scotland. New Phytologist 80, 455-484.

Examples

data(abernethy)
head(abernethy)

(plt <- Stratiplot(Age ~ . - Depth,
                   data = chooseTaxa(abernethy, n.occ = 5, max.abun = 10),
                   type = "poly"))

Description

Analogue matching is a more general implementation of the modern analogue methodology than MAT, where we are only interested in identifying sufficiently similar samples from a modern training as being suitable modern analogues for one or more fossil samples.

Usage

analog(x, ...)

## Default S3 method:
analog(x, y,
       method = c("euclidean", "SQeuclidean", "chord", "SQchord",
                  "bray", "chi.square", "SQchi.square",
                  "information", "chi.distance", "manhattan",
                  "kendall", "gower", "alt.gower", "mixed"),
       keep.train = TRUE, ...)

## S3 method for class 'distance'
analog(x, train = NULL, keep.train = TRUE, ...)

Arguments

Details

analog implements analogue matching sensu Flower et al (1997) and Simpson et al (2005), where the aim is to identify suitable close analogues of fossil samples from a modern training set. These results are generally used within ecological restoration, but the identification of close modern analogues for fossil samples is also used as a technique for assessing transfer function reconstructions.

analog is a simple and very general function that generates a pairwise dissimilarity matrix for the modern training set, and a second matrix containing the pairwise dissimilarities between each fossil sample and each sample in the training set. These results can then be assessed using other functions and to extract the close modern analogues using function cma. See the See Also section below.

Analysis of the pairwise dissimilarity matrix for the modern training set can be used to decide on a suitable dissimilarity threshold for defining close modern analogues. By default this matrix is returned as part of the output from the analog function.

Value

A list of class "analog" with the following components:

Author(s)

Gavin L. Simpson

References

Flower, R.J., Juggins, S. and Battarbee, R.W. (1997) Matching diatom assemblages in lake sediment cores and modern surface sediment samples: the implications for lake conservation and restoration with special reference to acidified systems. Hydrobiologia 344; 27–40.

Simpson, G.L., Shilland, E.M., Winterbottom, J. M. and Keay, J. (2005) Defining reference conditions for acidified waters using a modern analogue approach. Environmental Pollution 137; 119–133.

See Also

distance for the function that calculates the dissimilarity matrices. cma for extraction of close modern analogues. dissimilarities and plot.dissimilarities for analysis of distribution of pairwise dissimilarity matrix for modern training set.

Examples

## Imbrie and Kipp example
## load the example data
data(ImbrieKipp)
data(SumSST)
data(V12.122)

## merge training and test set on columns
dat <- join(ImbrieKipp, V12.122, verbose = TRUE)

## extract the merged data sets and convert to proportions
ImbrieKipp <- dat[[1]] / 100
V12.122 <- dat[[2]] / 100

## Imbrie and Kipp foraminfera sea-surface temperature

## analog matching between SWAP and RLGH core
ik.analog <- analog(ImbrieKipp, V12.122, method = "chord")
ik.analog
summary(ik.analog)

## Can take pre-computed dissimilarity objects
d1 <- distance(ImbrieKipp, V12.122)
d2 <- distance(ImbrieKipp)
ik <- analog(d1, d2, keep.train = TRUE)
ik

Description

Calculates Bayes factors or likelihood ratios of analogue and no-analogue results.

Usage

bayesF(x, prior = rep(0.5, 2))

## S3 method for class 'bayesF'
plot(x, group = "all", xlab = NULL, ylab = "Pr (A+ | d)",
        col = "red", abline.col = "lightgrey", abline.lty = "dashed", ...)

Arguments

Details

LR(+), is the likelihood ratio of a positive test result, that the value of d assigns the sample to the group it belongs to. LR(-) is the likelihood ratio of a negative test result, that the value of d assigns the sample to the wrong group.

LR(+) is defined as $LR(+) = TPF / FPF$ (or sensitivity / (1 - specificity)), and LR(-) is defined as $LR(-) = FPF / TNF$ (or (1 - sensitivity) / specificity), in Henderson (1993).

The posterior probability of analogue given a dissimilarity is the LR(+) likelihood ratio values multiplied by the prior odds of analogue, for given values of the dissimilarity, and is then converted to a probability.

The plotting function currently only draws the posterior probability of analogue based on the Bayes factor or likelihood ratio of a positive event (analogue).

Value

For plot.bayesF a plot on the currently active device.

For bayesF, a list containing the results of computing Bayes factors for each group in x. Each component of this list is itself a list with the following components:

Author(s)

Gavin L. Simpson

References

Brown, C.D., and Davis, H.T. (2006) Receiver operating characteristics curves and related decision measures: A tutorial. Chemometrics and Intelligent Laboratory Systems 80, 24–38.

Gavin, D.G., Oswald, W.W., Wahl, E.R. and Williams, J.W. (2003) A statistical approach to evaluating distance metrics and analog assignments for pollen records. Quaternary Research 60, 356–367.

Henderson, A.R. (1993) Assessing test accuracy and its clinical consequences: a primer for receiver operating characteristic curve analysis. Annals of Clinical Biochemistry 30, 834–846.

See Also

roc and plot.bayesF.

Examples

## load the example data
data(swapdiat, swappH, rlgh)

## merge training and test set on columns
dat <- join(swapdiat, rlgh, verbose = TRUE)

## extract the merged data sets and convert to proportions
swapdiat <- dat[[1]] / 100
rlgh <- dat[[2]] / 100

## fit an analogue matching (AM) model using the squared chord distance
## measure - need to keep the training set dissimilarities
swap.ana <- analog(swapdiat, rlgh, method = "SQchord",
                   keep.train = TRUE)

## fit the ROC curve to the SWAP diatom data using the AM results
## Generate a grouping for the SWAP lakes
METHOD <- if (getRversion() < "3.1.0") {"ward"} else {"ward.D"}
clust <- hclust(as.dist(swap.ana$train), method = METHOD)
grps <- cutree(clust, 12)

## fit the ROC curve
swap.roc <- roc(swap.ana, groups = grps)
swap.roc

## calculate the Bayes factors of analogue and no-analogue
## (uses observed probabilities of analogue/no-analogue
swap.bayes <- bayesF(swap.roc)
swap.bayes

## plot the probability of analogue
plot(swap.bayes)

## Not run: 
## calculate the Bayes factors of analogue and no-analogue
## with prior probabilities c(0.5, 0.05)
swap.bayes2 <- bayesF(swap.roc, prior = c(0.5, 0.05))
swap.bayes

## plot the probability of analogue
plot(swap.bayes2)

## End(Not run)

Description

Function to calculate bootstrap statistics for transfer function models such as bootstrap estimates, model RMSEP, sample specific errors for predictions and summary statistics such as bias and $R^2$ between oberved and estimated environment.

residuals method for objects of class "bootstrap.mat".

Usage

bootstrap(object, ...)

## Default S3 method:
bootstrap(object, ...)

## S3 method for class 'mat'
bootstrap(object, newdata, newenv, k,
          weighted = FALSE, n.boot = 1000, ...)

## S3 method for class 'bootstrap.mat'
fitted(object, k, ...)

## S3 method for class 'bootstrap.mat'
residuals(object, which = c("model", "bootstrap"), ...)

Arguments

Details

bootstrap is a fairly flexible function, and can be called with or without arguments newdata and newenv.

If called with only object specified, then bootstrap estimates for the training set data are returned. In this case, the returned object will not include component predictions.

If called with both object and newdata, then in addition to the above, bootstrap estimates for the new samples are also calculated and returned. In this case, component predictions will contain the apparent and bootstrap derived predictions and sample-specific errors for the new samples.

If called with object, newdata and newenv, then the full bootstrap object is returned (as described in the Value section below). With environmental data now available for the new samples, residuals, RMSE(P) and $R^2$ and bias statistics can be calculated.

The individual components of predictions are the same as those described in the components relating to the training set data. For example, returned.object$predictions$bootstrap contains the components as returned.object$bootstrap.

It is not usual for environmental data to be available for the new samples for which predictions are required. In normal palaeolimnological studies, it is more likely that newenv will not be available as we are dealing with sediment core samples from the past for which environmental data are not available. However, if sufficient training set samples are available to justify producing a training and a test set, then newenv will be available, and bootstrap can accomodate this extra information and calculate apparent and bootstrap estimates for the test set, allowing an independent assessment of the RMSEP of the model to be performed.

Typical usage of residuals is

    resid(object, which = c("model", "bootstrap"), \dots)

Value

For bootstrap.mat an object of class "bootstrap.mat" is returned. This is a complex object with many components and is described in bootstrapObject.

For residuals, a list containg the requested residuals and metadata, with the following components:

Author(s)

Gavin L. Simpson

References

Birks, H.J.B., Line, J.M., Juggins, S., Stevenson, A.C. and ter Braak, C.J.F. (1990). Diatoms and pH reconstruction. Philosophical Transactions of the Royal Society of London; Series B, 327; 263–278.

See Also

mat, plot.mat, summary.bootstrap.mat, residuals

Examples

## Imbrie and Kipp example
## load the example data
data(ImbrieKipp)
data(SumSST)
data(V12.122)

## merge training and test set on columns
dat <- join(ImbrieKipp, V12.122, verbose = TRUE)

## extract the merged data sets and convert to proportions
ImbrieKipp <- dat[[1]] / 100
V12.122 <- dat[[2]] / 100

## Imbrie and Kipp foraminfera sea-surface temperature 
## fit the MAT model using the squared chord distance measure
ik.mat <- mat(ImbrieKipp, SumSST, method = "SQchord")

## bootstrap training set
## IGNORE_RDIFF_BEGIN
ik.boot <- bootstrap(ik.mat, n.boot = 100)
ik.boot
summary(ik.boot)
## IGNORE_RDIFF_END

## Bootstrap fitted values for training set
## IGNORE_RDIFF_BEGIN
fitted(ik.boot)
## IGNORE_RDIFF_END

## residuals
resid(ik.boot) # uses abbreviated form

Description

Function to calculate bootstrap statistics for transfer function models such as bootstrap estimates, model RMSEP, sample specific errors for predictions and summary statistics such as bias and $R^2$ between oberved and estimated environment.

Usage

## S3 method for class 'wa'
bootstrap(object, n.boot = 1000, verbose = TRUE, ...)

Arguments

Details

See bootstrap.mat for further details. This method is not as feature packed as bootstrap.mat but can be used to evaluate the model performance of WA transfer function models.

Value

An object with the same components as predict.wa.

Author(s)

Gavin L. Simpson

References

Birks, H.J.B., Line, J.M., Juggins, S., Stevenson, A.C. and ter Braak, C.J.F. (1990). Diatoms and pH reconstruction. Philosophical Transactions of the Royal Society of London; Series B, 327; 263–278.

See Also

wa, plot.wa.

Examples

## Imbrie and Kipp
data(ImbrieKipp)
data(SumSST)
ik.wa <- wa(SumSST ~ ., data = ImbrieKipp, tol.dw = TRUE,
            min.tol = 2, small.tol = "min")
ik.wa

## compare actual tolerances to working values
with(ik.wa, rbind(tolerances, model.tol))

## bootstrap the WA model
ik.boot <- bootstrap(ik.wa, n.boot = 100)

## performance statistics
performance(ik.boot)

Description

Objects of class bootstrap.mat are a complex containing many sub-components. This object is described here in more detail.

Details

A large object is returned with some or all of the following depending on whether newdata and newenv are supplied or not.

observed:

vector of observed environmental values.

model:

a list containing the apparent or non-bootstrapped estimates for the training set. With the following components:

estimated:

estimated values for the response

residuals:

model residuals

r.squared:

Apparent $R^2$ between observed and estimated values of response

avg.bias:

Average bias of the model residuals

max.bias:

Maximum bias of the model residuals

rmse:

Apparent error (RMSE) for the model.

k:

numeric; indicating the size of model used in estimates and predictions

bootstrap:

a list containing the bootstrap estimates for the training set. With the following components:

estimated:

Bootstrap estimates for the response

residuals:

Bootstrap residuals for the response

r.squared:

Bootstrap derived $R^2$ between observed and estimated values of the response

avg.bias:

Average bias of the bootstrap derived model residuals

max.bias:

Maximum bias of the bootstrap derived model residuals

rmsep:

Bootstrap derived RMSEP for the model

s1:

Bootstrap derived S1 error component for the model

s2:

Bootstrap derived S2 error component for the model

k:

numeric; indicating the size of model used in estimates and predictions

sample.errors:

a list containing the bootstrap-derived sample specific errors for the training set. With the following components:

rmsep:

Bootstrap derived RMSEP for the training set samples

s1:

Bootstrap derived S1 error component for training set samples

s2:

Bootstrap derived S2 error component for training set samples

weighted:

logical; whether the weighted mean was used instead of the mean of the environment for k-closest analogues

auto:

logical; whether "k" was choosen automatically or user-selected

n.boot:

numeric; the number of bootstrap samples taken

call:

the matched call

type:

model type

predictions:

a list containing the apparent and bootstrap-derived estimates for the new data, with the following components:

observed:

the observed values for the new samples — only if newenv is provided

model:

a list containing the apparent or non-bootstrapped estimates for the new samples. A list with the same components as model, above

bootstrap:

a list containing the bootstrap estimates for the new samples, with some or all of the same components as bootstrap, above

sample.errors:

a list containing the bootstrap-derived sample specific errors for the new samples, with some or all of the same components as sample.errors, above

Author(s)

Gavin L. Simpson

See Also

mat, plot.mat, summary.bootstrap.mat, residuals

Description

Draws a caterpillar plot of the weighted average optima and tolerance range for each of the species in a data set.

Usage

## Default S3 method:
caterpillarPlot(x, tol, mult = 1, decreasing = TRUE,
                labels, xlab = NULL, pch = 21, bg = "white",
                col = "black", lcol = col, lwd = 2, frame.plot = FALSE, ...)

## S3 method for class 'data.frame'
caterpillarPlot(x, env, useN2 = TRUE, xlab, ...)

## S3 method for class 'wa'
caterpillarPlot(x, type = c("observed","model"), ...)

Arguments

Details

The function may also be called using the short form name caterpillar:

caterpillar(x, ...)
  

Value

The function results in a plot on the currently active device. A data frame with components Optima and Tolerance is returned invisibly.

Author(s)

Gavin L. Simpson

See Also

For the underlying computations optima and tolerance.

Examples

data(ImbrieKipp)
data(SumSST)

## default plot
caterpillar(ImbrieKipp, SumSST)

## customisation
opttol <-
    caterpillar(ImbrieKipp, SumSST, col = "red2",
                bg = "yellow", lcol = "blue",
                xlab = expression(Summer ~ Sea ~ Surface ~
                                 Temperature~(degree*C)))

## invisibly returns the optima and tolerances
head(opttol)

Description

Select taxa (variables) from an object on the basis of one or both of maximum abundance and number of occurrences greater than user-specified values. This is a simple utility function to encapsulate this common task in filtering palaeoecological data sets.

Usage

chooseTaxa(object, ...)

## Default S3 method:
chooseTaxa(object, n.occ = 1, max.abun = 0,
           type = c("AND","OR"), value = TRUE, na.rm = FALSE,
           ...)

Arguments

Value

If value = TRUE, returns the supplied data frame or matrix with a subset of columns (taxa) that meet the criteria chosen. If value = FALSE, a logical vector is returned.

Author(s)

Gavin L. Simpson

Examples

data(ImbrieKipp)
IK2 <- chooseTaxa(ImbrieKipp, n.occ = 5)
dim(ImbrieKipp)
dim(IK2)

## return a logical vector to select species/columns
chooseTaxa(ImbrieKipp, n.occ = 5, value = FALSE)

Description

Extracts and formats close modern analogue samples from a modern reference set that are closer than a defined cut off threshold.

Usage

cma(object, ...)

## Default S3 method:
cma(object, ...)

## S3 method for class 'analog'
cma(object, cutoff, prob = c(0.01, 0.025, 0.05), ...)

## S3 method for class 'mat'
cma(object, k, cutoff, prob = c(0.01, 0.025, 0.05), ...)

## S3 method for class 'predict.mat'
cma(object, k, cutoff, prob = c(0.01, 0.025,
0.05), ...)

## S3 method for class 'cma'
plot(x, method = c("overplot", "jitter", "stack"),
   jitter = 0.1, vertical = FALSE,
   draw.quant = TRUE, xlab = NULL, ylab = "",
   main = "", cex.axis = NULL, ...,
   col.quant = "red", lty.quant= "dashed")

Arguments

Details

The plot method is simply a wrapper to stripchart.

The methods for mat and predict.mat objects allow the user to select the k-closest analogues (argument k) or those samples as close or closer than a stated threshold of dissimilarity (argument cutoff). Only one of k and cutoff may be specified. If neither is specified, getK is used to extract the value for k stored within object. As such, the default is to return the automatically selected set of k closest samples, behaviour that is consistent with other functions in the package.

Value

For the plot method, a plot on the current device. Invisibly the plotted data are returned; see Note for further details.

A list of class "cma" with the following components:

Note

Only objects of classes analog, mat, and predict.mat are supported.

The plot method invisibly returns a list with the following components:

distances

a vector of stacked distances extracted from object.

groups

a factor listing the fossil sample for which the distances are the distances to the close modern analogues for the training set.

Author(s)

Gavin L. Simpson

References

Flower, R.J., Juggins, S. and Battarbee, R.W. (1997) Matching diatom assemblages in lake sediment cores and modern surface sediment samples: the implications for lake conservation and restoration with special reference to acidified systems. Hydrobiologia 344; 27–40.

Simpson, G.L., Shilland, E.M., Winterbottom, J. M. and Keay, J. (2005) Defining reference conditions for acidified waters using a modern analogue approach. Environmental Pollution 137; 119–133.

See Also

analog, stripchart, or boxplot for an alternative representation.

Examples

## Imbrie and Kipp example
## load the example data
data(ImbrieKipp)
data(SumSST)
data(V12.122)

## merge training and test set on columns
dat <- join(ImbrieKipp, V12.122, verbose = TRUE)

## extract the merged data sets and convert to proportions
ImbrieKipp <- dat[[1]] / 100
V12.122 <- dat[[2]] / 100

## analog matching between SWAP and RLGH reference samples
(ik.ana <- analog(ImbrieKipp, V12.122, method = "chord"))

## close modern analogues
(ik.cma <- cma(ik.ana, cutoff = 0.4))
summary(ik.cma)

## plot the results
plot(ik.cma)

Description

compare() compares a proxy dataset with a training set or other data set that is considered to be the master. A range of metrics is returned, currently for samples only.

Usage

compare(x, ...)

## Default S3 method:
compare(x, y, env,
        by = c("sites", "species"),
        ordination = "rda",
        method = "chord",
        transform = NULL,
        n2limit = 5L,
        ...)

Arguments

Details

ToDo

Value

If by = "species" a data frame of diagnostics for each species (proxy) in y relative to x. If by = "sites", the diagnostics are for each sample (row) in y. Depending on the value of by some of the following columns will be returned

Author(s)

Gavin L. Simpson

Examples

data(ImbrieKipp, V12.122, SumSST)
compare(ImbrieKipp, V12.122, env = SumSST, ordination = "rda",
        method = "chord")

Description

Performs leave-one-out, k-fold, n k-fold and bootstrap cross-validation of palaeoecological transfer function models.

Usage

crossval(obj, ...)

## S3 method for class 'wa'
crossval(obj, method = c("LOO","kfold","bootstrap"),
         nboot = 100, nfold = 10, folds = 5,
         verbose = getOption("verbose"), ...)

## S3 method for class 'pcr'
crossval(obj, method = c("LOO","kfold","bootstrap"),
         ncomp, nboot = 100, nfold = 10, folds = 5,
         verbose = getOption("verbose"), ...)

Arguments

Value

Returns an object of class "crossval", a list with the following components:

Author(s)

Gavin L. Simpson

See Also

wa

Examples

## Load the Imbrie & Kipp data and
## summer sea-surface temperatures
data(ImbrieKipp)
data(SumSST)
     
## fit the WA model
mod <- wa(SumSST ~., data = ImbrieKipp)
mod

## Leave one out CV
cv.loo <- crossval(mod)
cv.loo

## k-fold CV (k == 10)
cv.kfold <- crossval(mod, method = "kfold", kfold = 10, folds = 1)
cv.kfold

## n k-fold CV (k == 10, n = 10)
cv.nkfold <- crossval(mod, method = "kfold", kfold = 10, folds = 10)
cv.nkfold

## bootstrap with 100 bootstrap samples
cv.boot <- crossval(mod, method = "bootstrap", nboot = 100)
cv.boot

## extract fitted values and residuals
fitted(cv.boot)
resid(cv.boot)

## Principal Components Regression
mpcr <- pcr(SumSST ~., data = ImbrieKipp, ncomp = 10)
crossval(mpcr, method = "kfold", kfold = 10, folds = 2, ncomp = 10)

crossval(mpcr, method = "bootstrap", nboot = 100, ncomp = 10)

Description

Lattice densityplot method for residLen objects.

Usage

## S3 method for class 'residLen'
densityplot(x, ..., xlab = NULL, ylab = NULL)

Arguments

Value

Returns an object of class "trellis". See densityplot for details.

Author(s)

Gavin L. Simpson

See Also

residLen, plot.residLen, hist.residLen, histogram.residLen.

Examples

## load the Imbrie and Kipp example data
data(ImbrieKipp, SumSST, V12.122)

## squared residual lengths for Core V12.122
rlens <- residLen(ImbrieKipp, SumSST, V12.122)
rlens

## plot the density functions of the residual distances
densityplot(rlens)

Description

In Weighted Averaging models averages are taken twice and thus WA estimates shrink towards the training set mean and need to be deshrunk.deshrink performs this deshrinking using several techniques, whilst deshrinkPred will deshrink WA estimates for new samples given a set of deshrinking coefficients.

Usage

deshrink(env, wa.env,
         type = c("inverse", "classical", "expanded", "none",
                  "monotonic"))

deshrinkPred(x, coef,
         type = c("inverse", "classical", "expanded", "none",
                  "monotonic"))

Arguments

Value

For deshrinkPred a numeric vector of deshrunk estimates.

For an object of class "deshrink", inheriting from class "list", with two components. The type of deshrinking performed is stroed within attribute "type". The componets of the returned object are:

Warning

deshrinkPred, does not currently check that the correct coefficients have been supplied in the correct order.

Author(s)

Gavin L. Simpson & Jari Oksanen

References

Birks, H.J.B. (1995) Quantitative environmental reconstructions. In Statistical modelling of Quaternary science data (eds.~D.Maddy & J.S. Brew). Quaternary Research Association technical guide 5. Quaternary Research Association, Cambridge.

See Also

wa

Description

Extracts a vector of dissimilarity coefficients from an object for further analysis.

Usage

dissimilarities(object, ...)
dissim(object, ...)

## S3 method for class 'analog'
dissimilarities(object, which = c("train", "analogs"),
                 ...)

## S3 method for class 'mat'
dissimilarities(object, ...)

Arguments

Details

The function can be called using the much shorter name "dissim".

Value

A vector of dissimilarities.

Author(s)

Gavin L. Simpson

See Also

analog, plot.dissimilarities

Examples

## Imbrie and Kipp example
## load the example data
data(ImbrieKipp)
data(SumSST)
data(V12.122)

## merge training and test set on columns
dat <- join(ImbrieKipp, V12.122, verbose = TRUE)

## extract the merged data sets and convert to proportions
ImbrieKipp <- dat[[1]] / 100
V12.122 <- dat[[2]] / 100

## analog matching between SWAPImbrie & Kipp and V12.122 core
ik.analog <- analog(ImbrieKipp, V12.122, method = "chord")
ik.analog
summary(ik.analog)

## compare training set dissimilarities with normals
## and derive cut-offs
ik.dij <- dissim(ik.analog)
plot(ik.dij)

Description

Flexibly calculates distance or dissimilarity measures between a training set x and a fossil or test set y. If y is not supplied then the pairwise dissimilarities between samples in the training set, x, are calculated.

Usage

distance(x, ...)

## Default S3 method:
distance(x, y, method = "euclidean", weights = NULL,
         R = NULL, dist = FALSE, double.zero = FALSE, ...)

## S3 method for class 'join'
distance(x, ...)

oldDistance(x, ...)
## Default S3 method:
oldDistance(x, y, method = c("euclidean", "SQeuclidean",
            "chord", "SQchord", "bray", "chi.square",
            "SQchi.square", "information", "chi.distance",
            "manhattan", "kendall", "gower", "alt.gower",
            "mixed"),
            fast = TRUE,
            weights = NULL, R = NULL, ...)
## S3 method for class 'join'
oldDistance(x, ...)

Arguments

Details

A range of dissimilarity coefficients can be used to calculate dissimilarity between samples. The following are currently available:

Argument fast determines whether fast C versions of some of the dissimilarity coefficients are used. The fast versions make use of dist for methods "euclidean", "SQeuclidean", "chord", "SQchord", and vegdist for method == "bray". These fast versions are used only when x is supplied, not when y is also supplied. Future versions of distance will include fast C versions of all the dissimilary coefficients and for cases where y is supplied.

Value

A matrix of dissimilarities where columns are the samples in y and the rows the samples in x. If y is not provided then a square, symmetric matrix of pairwise sample dissimilarities for the training set x is returned, unless argument dist is TRUE, in which case an object of class "dist" is returned. See dist.

The dissimilarity coefficient used (method) is returned as attribute "method". Attribute "type" indicates whether the object was computed on a single data matrix ("symmetric") or across two matrices (i.e. the dissimilarties between the rows of two matrices; "asymmetric".

Warning

For method = "mixed" it is essential that a factor in x and y have the same levels in the two data frames. Previous versions of analogue would work even if this was not the case, which will have generated incorrect dissimilarities for method = "mixed" for cases where factors for a given species had different levels in x to y.

distance now checks for matching levels for each species (column) recorded as a factor. If the factor for any individual species has different levels in x and y, an error will be issued.

Note

The dissimilarities are calculated in native R code. As such, other implementations (see See Also below) will be quicker. This is done for one main reason - it is hoped to allow a user defined function to be supplied as argument "method" to allow for user-extension of the available coefficients.

The other advantage of distance over other implementations, is the simplicity of calculating only the required pairwise sample dissimilarities between each fossil sample (y) and each training set sample (x). To do this in other implementations, you would need to merge the two sets of samples, calculate the full dissimilarity matrix and then subset it to achieve similar results.

Author(s)

Gavin L. Simpson and Jari Oksanen (improvements leading to method "metric.mixed" and proper handling of ordinal data via Podani's (1999) modification of Gower's general coefficient in method "mixed").

References

Faith, D.P., Minchin, P.R. and Belbin, L. (1987) Compositional dissimilarity as a robust measure of ecological distance. Vegetatio 69, 57–68.

Gavin, D.G., Oswald, W.W., Wahl, E.R. and Williams, J.W. (2003) A statistical approach to evaluating distance metrics and analog assignments for pollen records. Quaternary Research 60, 356–367.

Kendall, D.G. (1970) A mathematical approach to seriation. Philosophical Transactions of the Royal Society of London - Series B 269, 125–135.

Legendre, P. and Legendre, L. (1998) Numerical Ecology, 2nd English Edition. Elsevier Science BV, The Netherlands.

Overpeck, J.T., Webb III, T. and Prentice I.C. (1985) Quantitative interpretation of fossil pollen spectra: dissimilarity coefficients and the method of modern analogues. Quaternary Research 23, 87–108.

Podani, J. (1999) Extending Gower's General Coefficient of Similarity to Ordinal Characters. Taxon 48, 331–340).

Prentice, I.C. (1980) Multidimensional scaling as a research tool in Quaternary palynology: a review of theory and methods. Review of Palaeobiology and Palynology 31, 71–104.

See Also

vegdist in package vegan, daisy in package cluster, and dist provide comparable functionality for the case of missing y.

Examples

## simple example using dummy data
train <- data.frame(matrix(abs(runif(200)), ncol = 10))
rownames(train) <- LETTERS[1:20]
colnames(train) <- as.character(1:10)
fossil <- data.frame(matrix(abs(runif(100)), ncol = 10))
colnames(fossil) <- as.character(1:10)
rownames(fossil) <- letters[1:10]

## calculate distances/dissimilarities between train and fossil
## samples
test <- distance(train, fossil)

## using a different coefficient, chi-square distance
test <- distance(train, fossil, method = "chi.distance")

## calculate pairwise distances/dissimilarities for training
## set samples
test2 <- distance(train)

## Using distance on an object of class join
dists <- distance(join(train, fossil))
str(dists)

## calculate Gower's general coefficient for mixed data
## first, make a couple of variables factors

## fossil[,4] <- factor(sample(rep(1:4, length = 10), 10))
## train[,4] <- factor(sample(rep(1:4, length = 20), 20))
## ## now fit the mixed coefficient
## test3 <- distance(train, fossil, "mixed")

## ## Example from page 260 of Legendre & Legendre (1998)
x1 <- t(c(2,2,NA,2,2,4,2,6))
x2 <- t(c(1,3,3,1,2,2,2,5))
Rj <- c(1,4,2,4,1,3,2,5) # supplied ranges

## 1 - distance(x1, x2, method = "mixed", R = Rj)

## note this gives ~0.66 as Legendre & Legendre describe the
## coefficient as a similarity coefficient. Hence here we do
## 1 - Dij here to get the same answer.

## Tortula example from Podani (1999)
data(tortula)
Dij <- distance(tortula[, -1], method = "mixed") # col 1 includes Taxon ID

## Only one ordered factor
data(mite.env, package = "vegan")
Dij <- distance(mite.env, method = "mixed")

## Some variables are constant
data(BCI.env, package = "vegan")
Dij <- distance(BCI.env, method = "mixed")

Description

The number of samples in sections along the gradient is a useful diagnostic as to the quality of reconstructions at gradient values within those sections.

Usage

evenSample(grad, n = 10)

Arguments

Details

The sampling design of a training set, i.e. the number of samples taken at points along the gradient, can influence the uncertainty in the transfer function predictions at those values of the gradient. Poorly sampled sections of the gradient may have far larger RMSEP than the overall model RMSEP.

Value

Numeric vector of length n containing the numbers of samples per gradient segment.

Author(s)

Gavin L. Simpson

Examples

data(SumSST)
ev <- evenSample(SumSST) ## not an even sample...
plot(ev)

Description

Extracts fitted values for training set samples from logistic regression models fitted to each group of samples that describe the probability two samples are analogues (from the same group) as a function of dissimilarity between the paired samples.

Usage

## S3 method for class 'logitreg'
fitted(object, combined = FALSE, ...)

Arguments

Value

If combined == FALSE (the default) then a matrix of fitted probabilities, where the rows are the training set samples and the columns the groupings, is returned. If combined == TRUE, then a list with components "group" and "combined". "group" is a matrix of fitted probabilities as above. "combined" is a vector of fitted values for the entire set of pairwise comparisons considered.

Author(s)

Gavin L. Simpson

See Also

See logitreg for example usage.

Description

Combines dissimilarities from two or more dissimilarity objects into a single dissimilarity object so that both original dissimilarities contribute equally. Weighted combinations of the original objects can also be created.

Usage

fuse(..., weights = NULL)

## S3 method for class 'matrix'
fuse(..., weights = NULL)

## S3 method for class 'dist'
fuse(..., weights = NULL)

Arguments

Details

Fuses, or combines, dissimilarity objects in a very flexible way to create a single dissimilarity object that incorporates the separate dissimilarities. In analogue matching, we may wish to combine information from two or more proxies, such as diatoms and cladocera, or from biological and chemical or physical data in the case of matching modern samples.

The function can also be used to fuse dissimilarity objects created from a single data set but using different dissimilarity coefficients. In this way one could create a new dissimilarity object combining dissimilarity based on abundance data and presence absence data into a single measure.

fuse uses the method of Melssen et al. (2006) to combine dissimilarities. The dissimilarities in each dissimilarity object are scaled so that the maximum dissimilarity in each object is 1. The scaled dissimilarity objects are then weighted according to the supplied weights. If no weights are supplied (the default) the dissimilarity objects are weighted equally; weights = rep(1/N, N), where N is the number of dissimilarity objects fused.

$D_{fused}(j, k) = \sum_{i = 1}^N w_i D_{ijk}$

where $D_{fused}(j, k)$ is the fused dissimilarity between samples $j$ and $k$, $w_i$ is the weight assigned to the $i$th dissimilarity object and $D_{ijk}$ is the dissimilarity between $j$ and $k$ for the $i$th dissimilarity object.

Value

fuse returns an object of class "dist" with the attribute "method" set to "fuse".

This is the case even if the supplied objects are full dissimilarity matrices. If you want a full dissimilarity object, use as.matrix.dist on the returned object.

The returned object contains an extra attribute "weights", which records the weights used in the fusing.

Author(s)

Gavin L. Simpson

References

Melssen W., Wehrens R. and Buydens L. (2006) Supervised Kohonen networks for classification problems. Chemometrics and intelligent laboratory systems 83, 99–113.

See Also

dist, vegdist, distance.

Examples

train1 <- data.frame(matrix(abs(runif(100)), ncol = 10))
train2 <- data.frame(matrix(sample(c(0,1), 100, replace = TRUE),
                     ncol = 10))
rownames(train1) <- rownames(train2) <- LETTERS[1:10]
colnames(train1) <- colnames(train2) <- as.character(1:10)

d1 <- vegdist(train1, method = "bray")
d2 <- vegdist(train2, method = "jaccard")

dd <- fuse(d1, d2, weights = c(0.6, 0.4))
dd
str(dd)

Description

An extractor function to access the number of analogues used in particular models. The stored value of $k$ can be updated using setK.

Usage

getK(object, ...)

## S3 method for class 'mat'
getK(object, weighted = FALSE, ...)

## S3 method for class 'bootstrap.mat'
getK(object, which = c("bootstrap", "model"),
     prediction = FALSE, ...)

## S3 method for class 'predict.mat'
getK(object, which = c("model", "bootstrap"),
     ...)

setK(object, weighted = FALSE) <- value

## S3 replacement method for class 'mat'
setK(object, weighted = FALSE) <- value

Arguments

Details

getK is a generic accessor function, and setK<- is a generic replacement function.

Objects of class bootstrap.mat contain several different k's. If no predictions are performed, there will be two k's, one for the model and one from bootstrapping the model. Where predictions are performed with newenv supplied, in addition to the k's above, there will be two k' for the predictions, one for the model-based and one for the bootstrap-based predictions. To select k for the predictions, use prediction = TRUE. Argument which determines whether the model-based or the bootstrap-based k is returned.

Value

For getK, an integer value that is the number of analogues stored for use. The returned object has attributes “auto” and “weighted”. “auto” refers to whether the extracted value of $k$ was set automatically (TRUE) or by the user (FALSE). “weighted” states if the returned value is for a weighted analysis or an un-weighted analysis (FALSE).

For setK<-, the updated object.

Author(s)

Gavin L. Simpson

See Also

mat

Examples

## Imbrie and Kipp Sea Surface Temperature
data(ImbrieKipp)
data(SumSST)
data(V12.122)

## merge training set and core samples
dat <- join(ImbrieKipp, V12.122, verbose = TRUE)

## extract the merged data sets and convert to proportions
ImbrieKipp <- dat[[1]] / 100
ImbrieKippCore <- dat[[2]] / 100

## fit a MAT model
ik.mat <- mat(ImbrieKipp, SumSST, method = "chord")

## How many analogues gives lowest RMSE?
getK(ik.mat)
## note that this value was chosen automatically

## Now set k to be 10
setK(ik.mat) <- 10

## check
getK(ik.mat)

Description

Extracts information as to the locations of samples along an ordination gradient. gradientDist() standardises the entire gradient to the interval 0, ..., 1, to allow comparison between methods or data sets.

Usage

gradientDist(object, ...)

## Default S3 method:
gradientDist(object, na.rm = TRUE, ...)

## S3 method for class 'cca'
gradientDist(object, na.rm = TRUE, axis = 1L,
             scaling = 0, ...)

## S3 method for class 'prcurve'
gradientDist(object, na.rm = TRUE, ...)

Arguments

Value

A numeric vector of positions along the gradient, scaled to the range 0, ..., 1.

Author(s)

Gavin L. Simpson

See Also

See cca and prcurve for functions that produce objects that gradientDist() can work with.

Examples

data(abernethy)

## Remove the Depth and Age variables
abernethy2 <- abernethy[, -(37:38)]

## Fit PCA
aber.pca <- rda(abernethy2)

## Distance along the first PCA axis
gradientDist(aber.pca)

Description

Base graphics histogram plot method for residLen objects.

Usage

## S3 method for class 'residLen'
hist(x, breaks = "Sturges", freq = TRUE,
     probs = c(0.9, 0.95, 0.99), ncol = 1, lcol = "red",
     llty = "dashed", xlab = NULL, ylab = NULL,
     main = "Residual distances", rug = TRUE, ...)

Arguments

Value

A plot on the current device.

Returns a list with two components (train and passive), each of which is an object returned by hist.

Author(s)

Gavin L. Simpson

See Also

residLen, plot.residLen, histogram.residLen, densityplot.residLen.

Examples

## load the Imbrie and Kipp example data
data(ImbrieKipp, SumSST, V12.122)

## squared residual lengths for Core V12.122
rlens <- residLen(ImbrieKipp, SumSST, V12.122)
rlens

## plot a histogram of the residual distances
hist(rlens)

Description

Lattice histogram method for residLen objects.

Usage

## S3 method for class 'residLen'
histogram(x, ..., xlab = NULL, ylab = NULL,
          type = c("percent", "count", "density"))

Arguments

Value

Returns an object of class "trellis". See histogram for details.

Author(s)

Gavin L. Simpson

See Also

residLen, plot.residLen, hist.residLen, densityplot.residLen.

Examples

## load the Imbrie and Kipp example data
data(ImbrieKipp, SumSST, V12.122)

## squared residual lengths for Core V12.122
rlens <- residLen(ImbrieKipp, SumSST, V12.122)
rlens

## plot a histogram of the residual distances
histogram(rlens)

Description

The classic Imbrie and Kipp (1971) training set of counts on 27 species of foraminifera from 61 ocean sediment core surface samples and associated measures of summer and winter sea-surface temperatures and salinity at each location.

110 sediment cores samples from core V12-122 are also supplied in V12.122.

Usage

data(SumSST)
data(WinSST)
data(Salinity)
data(V12.122)

Format

ImbrieKipp is a data frame with 61 observations on the following 27 species:

O.univ

a numeric vector

G.cglob

a numeric vector

G.ruber

a numeric vector

G.tenel

a numeric vector

G.saccu

a numeric vector

G.rubes

a numeric vector

G.pacL

a numeric vector

G.pacR

a numeric vector

G.bullo

a numeric vector

G.falco

a numeric vector

G.calid

a numeric vector

G.aequi

a numeric vector

G.gluti

a numeric vector

G.duter

a numeric vector

G.infla

a numeric vector

G.trnL

a numeric vector

G.trnR

a numeric vector

G.crasf

a numeric vector

G.scitu

a numeric vector

G.mentu

a numeric vector

P.obliq

a numeric vector

C.nitid

a numeric vector

S.dehis

a numeric vector

G.digit

a numeric vector

Other

a numeric vector

G.quin

a numeric vector

G.hirsu

a numeric vector

Summer and Winter sea-surface temperatures, and salinity values for the 61 sites in the Imbrie and Kipp training set (ImbrieKipp):

SumSST

a numeric vector of summer sea-surface water temperatures. Values are in degrees C.

WinSST

a numeric vector of winter sea-surface water temperatures. Values are in degrees C.

Salinity

a numeric vector of sea water salinity values.

V12.122 is a data frame with 110 observations from core V12-122 on the following 28 species:

O.univ

a numeric vector

G.cglob

a numeric vector

G.ruber

a numeric vector

G.tenel

a numeric vector

G.saccu

a numeric vector

G.rubes

a numeric vector

G.pacL

a numeric vector

G.pacR

a numeric vector

G.bullo

a numeric vector

G.falco

a numeric vector

G.calid

a numeric vector

G.aequi

a numeric vector

G.gluti

a numeric vector

G.duter

a numeric vector

G.infla

a numeric vector

G.trnL

a numeric vector

G.trnR

a numeric vector

G.crasf

a numeric vector

G.scitu

a numeric vector

G.mentu

a numeric vector

P.obliq

a numeric vector

C.nitid

a numeric vector

S.dehis

a numeric vector

G.digit

a numeric vector

G.hexag

a numeric vector

G.cglom

a numeric vector

cfH.pel

a numeric vector

Other

a numeric vector

Source

Imbrie and Kipp (1971) TODO: Get the full reference.

These data were provided in electronic format by Prof. H. John B. Birks.

Examples

data(ImbrieKipp)
head(ImbrieKipp)

data(SumSST)
data(WinSST)
data(Salinity)

plot(cbind(SumSST, WinSST, Salinity))

data(V12.122)
head(V12.122)

Description

Merges any number of species matrices on their common columns to create a new data set with number of columns equal to the number of unqiue columns across all data frames. Needed for analysis of fossil data sets with respect to training set samples.

Usage

join(..., verbose = FALSE, na.replace = TRUE, split = TRUE, value = 0,
     type = c("outer", "left", "inner"))

## S3 method for class 'join'
head(x, ...)

## S3 method for class 'join'
tail(x, ...)

Arguments

Details

When merging multiple data frames the set of variables in the merged data can be determined via a number of routes. join provides for two (currently) join types; the outer join and the left outer (or simply the left) join. Which type of join is performed is determined by the argument type.

The outer join returns the union of the set of variables found in the data frames to be merged. This means that the resulting data frame(s) contain columns for all the variable observed across all the data frames supplied for merging.

With the left outer join the resulting data frame(s) contain only the set of variables found in the first data frame provided.

The inner join returns the intersection of the set of variables found in the supplied data frames. The resulting data frame(s) contains the variables common to all supplied data frames.

Value

If split = TRUE, an object of class "join", a list of data frames, with as many components as the number of data frames originally merged.

Otherwise, an object of class c("join", "data.frame"), a data frame containing the merged data sets.

head.join and tail.join return a list, each component of which is the result of a call to head or tail on each data set compont of the joined object.

Author(s)

Gavin L. Simpson

See Also

merge

Examples

## load the example data
data(swapdiat, swappH, rlgh)

## merge training and test set on columns
dat <- join(swapdiat, rlgh, verbose = TRUE)

## extract the merged data sets and convert to proportions
swapdiat <- dat[[1]] / 100
rlgh <- dat[[2]] / 100

## merge training and test set using left join
head(join(swapdiat, rlgh, verbose = TRUE, type = "left"))

## load the example data
data(ImbrieKipp, SumSST, V12.122)

## merge training and test set on columns
dat <- join(ImbrieKipp, V12.122, verbose = TRUE)

## extract the merged data sets and convert to proportions
ImbrieKipp <- dat[[1]] / 100
V12.122 <- dat[[2]] / 100

## show just the first few lines of each data set
head(dat, n = 4)

## show just the last few lines of each data set
tail(dat, n = 4)

## merge training and test set using inner join
head(join(ImbrieKipp, V12.122, verbose = TRUE, type = "inner"))

## merge training and test set using outer join and replace
## NA with -99.9
head(join(ImbrieKipp, V12.122, verbose = TRUE, value = -99.9))

Description

Fits logistic regression models to each level of group to model the probability of two samples being analogues conditional upon the dissimilarity between the two samples.

Usage

logitreg(object, groups, k = 1, ...)

## Default S3 method:
logitreg(object, groups, k = 1,
         biasReduced = FALSE, ...)

## S3 method for class 'analog'
logitreg(object, groups, k = 1, ...)

## S3 method for class 'logitreg'
summary(object, p = 0.9, ...)

Arguments

Details

Fits logistic regression models to each level of group to model the probability of two samples being analogues (i.e. in the same group) conditional upon the dissimilarity between the two samples.

This function can be seen as a way of directly modelling the probability that two sites are analogues, conditional upon dissimilarity, that can also be done less directly using roc and bayesF.

Often, the number of true analogues in the training set is small, both in absolute terms and as a proportion of comparisons. Logistic regression is known to suffer from a small-sample bias. Firth's method of bias reduction is a general solution to this problem and is implemented in logitreg through the brglm package of Ioannis Kosmidis.

Value

logitreg returns an object of class "logitreg"; a list whose components are objects returned by glm. See glm for further details on the returned objects.

The components of this list take their names from group.

For summary.logitreg an object of class "summary.logitreg", a data frame with summary statistics of the model fits. The components of this data frame are:

Note

The function may generate warnings from function glm.fit. These should be investigated and not simply ignored.

If the message is concerns fitted probabilities being numerically 0 or 1, then check the fitted values of each of the models. These may well be numerically 0 or 1. Heed the warning in glm and read the reference cited therein which may indicate problems with the fitted models, such as (quasi-)complete separation.

Author(s)

Gavin L. Simpson

References

Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika 80, 27-38.

See Also

roc, bayesF, glm, and brglm.

Examples

## load the example data
data(swapdiat, swappH, rlgh)

## merge training and test set on columns
dat <- join(swapdiat, rlgh, verbose = TRUE)

## extract the merged data sets and convert to proportions
swapdiat <- dat[[1]] / 100
rlgh <- dat[[2]] / 100

## fit an analogue matching (AM) model using the squared chord distance
## measure - need to keep the training set dissimilarities
swap.ana <- analog(swapdiat, rlgh, method = "SQchord",
                   keep.train = TRUE)

## fit the ROC curve to the SWAP diatom data using the AM results
## Generate a grouping for the SWAP lakes
METHOD <- if (getRversion() < "3.1.0") {"ward"} else {"ward.D"}
clust <- hclust(as.dist(swap.ana$train), method = METHOD)
grps <- cutree(clust, 6)

## fit the logit models to the analog object
swap.lrm <- logitreg(swap.ana, grps)
swap.lrm

## summary statistics
summary(swap.lrm)

## plot the fitted logit curves
plot(swap.lrm, conf.type = "polygon")

## extract fitted posterior probabilities for training samples
## for the individual groups
fit <- fitted(swap.lrm)
head(fit)

## compute posterior probabilities of analogue-ness for the rlgh
## samples. Here we take the dissimilarities between fossil and
## training samples from the `swap.ana` object rather than re-
## compute them
pred <- predict(swap.lrm, newdata = swap.ana$analogs)
head(pred)

## Bias reduction
## fit the logit models to the analog object
swap.brlrm <- logitreg(swap.ana, grps, biasReduced = TRUE)
summary(swap.brlrm)

Description

Modern Analogue Technique (MAT) transfer function models for palaeoecology. The fitted values are the, possibly weighted, averages of the environment for the k-closest modern analogues. MAT is a k-NN method.

Usage

mat(x, ...)

## Default S3 method:
mat(x, y,
    method = c("euclidean", "SQeuclidean", "chord", "SQchord",
               "bray", "chi.square", "SQchi.square",
               "information", "chi.distance", "manhattan",
               "kendall", "gower", "alt.gower", "mixed"),
    kmax, ...)

## S3 method for class 'formula'
mat(formula, data, subset, na.action,
    method = c("euclidean", "SQeuclidean", "chord", "SQchord",
               "bray", "chi.square", "SQchi.square",
               "information", "chi.distance", "manhattan",
               "kendall", "gower", "alt.gower", "mixed"),
    model = FALSE, ...)

## S3 method for class 'mat'
fitted(object, k, weighted = FALSE, ...)

## S3 method for class 'mat'
residuals(object, k, weighted = FALSE, ...)

Arguments

Details

The Modern Analogue Technique (MAT) is perhaps the simplest of the transfer function models used in palaeoecology. An estimate of the environment, $x$, for the response for a fossil sample, $y$, is the, possibly weighted, mean of that variable across the k-closest modern analogues selected from a modern training set of samples. If used, weights are the reciprocal of the dissimilarity between the fossil sample and each modern analogue.

A typical model has the form response ~ terms where response is the (numeric) response data frame and terms is a series of terms which specifies a linear predictor for response. A typical form for terms is ., which is shorthand for "all variables" in data. If . is used, data must also be provided. If specific species (variables) are required then terms should take the form spp1 + spp2 + spp3.

Pairwise sample dissimilarity is defined by dissimilarity or distance coefficients. A variety of coefficients are supported — see distance for details of the supported coefficients.

k is chosen by the user. The simplest choice for k is to evaluate the RMSE of the difference between the predicted and observed values of the environmental variable of interest for the training set samples for a sequence of models with increasing k. The number of analogues chosen is the value of k that has lowest RMSE. However, it should be noted that this value is biased as the data used to build the model are also used to test the predictive power.

An alternative approach is to employ an optimisation data set on which to evaluate the size of $k$ that provides the lowest RMSEP. This may be impractical with smaller sample sizes.

A third option is to bootstrap re-sample the training set many times. At each bootstrap sample, predictions for samples in the bootstrap test set can be made for $k = 1, ..., n$, where $n$ is the number of samples in the training set. $k$ can be chosen from the model with the lowest RMSEP. See function bootstrap.mat for further details on choosing $k$.

The output from summary.mat can be used to choose $k$ in the first case above. For predictions on an optimsation or test set see predict.mat. For bootstrap resampling of mat models, see bootstrap.mat.

The fitted values are for the training set and are taken as the, possibly weighted, mean of the environmental variable in question across the k-closest analogues. The fitted value for each sample does not include a contribution from itself — it is the closest analogue, having zero dissimilarity. This spurious distance is ignored and analogues are ordered in terms of the non-zero distances to other samples in the training set, with the k-closest contributing to the fitted value.

Typical usages for residuals.mat are:

    resid(object, k, weighted = FALSE, \dots)
  

Value

mat returns an object of class mat with the following components:

If model = TRUE then additional components "terms" and "model" are returned containing the terms object and model frame used.

fitted.mat returns a list with the following components:

Note

The object returned by mat contains lists "standard" and "weighted" both with the following elements:

est

a matrix of estimated values for the training set samples for models using $k$ analogues, where $k = 1, ..., n$. $n$ is the number of smaples in the training set. Rows contain the values for each model of size $k$, with colums containing the estimates for each training set sample.

resid

matrix; as for "est", but containing the model residuals.

rmsep

vector; containing the leave-one-out root mean square error or prediction.

avg.bias

vector; contains the average bias (mean of residuals) for models using k analogues, where $k = 1, ..., n$. $n$ is the number of smaples in the training set.

max.bias

vector; as for "avg.bias", but containing the maximum bias statistics.

r.squared

vector; as for "avg.bias", but containing the $R^2$ statistics.

Author(s)

Gavin L. Simpson

References

Gavin, D.G., Oswald, W.W., Wahl, E.R. and Williams, J.W. (2003) A statistical approach to evaluating distance metrics and analog assignments for pollen records. Quaternary Research 60, 356–367.

Overpeck, J.T., Webb III, T. and Prentice I.C. (1985) Quantitative interpretation of fossil pollen spectra: dissimilarity coefficients and the method of modern analogues. Quaternary Research 23, 87–108.

Prell, W.L. (1985) The stability of low-latitude sea-surface temperatures: an evaluation of the CLIMAP reconstruction with emphasis on the positive SST anomalies, Report TR 025. U.S. Department of Energy, Washington, D.C.

Sawada, M., Viau, A.E., Vettoretti, G., Peltier, W.R. and Gajewski, K. (2004) Comparison of North-American pollen-based temperature and global lake-status with CCCma AGCM2 output at 6 ka. Quaternary Science Reviews 23, 87–108.

See Also

summary.mat, bootstrap.mat for boostrap resampling of MAT models, predict.mat for making predictions from MAT models, fitted.mat and resid.mat for extraction of fitted values and residuals from MAT models respectively. plot.mat provides a plot.lm-like plotting tool for MAT models.

Examples

## Imbrie and Kipp Sea Surface Temperature
data(ImbrieKipp)
data(SumSST)
data(V12.122)

## merge training set and core samples
dat <- join(ImbrieKipp, V12.122, verbose = TRUE)

## extract the merged data sets and convert to proportions
ImbrieKipp <- dat[[1]] / 100
ImbrieKippCore <- dat[[2]] / 100

## fit the MAT model using the squared chord distance measure
ik.mat <- mat(ImbrieKipp, SumSST, method = "chord")
ik.mat

## model summary
summary(ik.mat)

## fitted values
fitted(ik.mat)

## model residuals
resid(ik.mat)

## draw summary plots of the model
par(mfrow = c(2,2))
plot(ik.mat)
par(mfrow = c(1,1))

## reconstruct for the V12.122 core data
coreV12.mat <- predict(ik.mat, V12.122, k = 3)
coreV12.mat
summary(coreV12.mat)

## draw the reconstruction
reconPlot(coreV12.mat, use.labels = TRUE, display.error = "bars",
          xlab = "Depth", ylab = "SumSST")

## fit the MAT model using the squared chord distance measure
## and restrict the number of analogues we fit models for to 1:20
ik.mat2 <- mat(ImbrieKipp, SumSST, method = "chord", kmax = 20)
ik.mat2

Description

Permutations and Monte Carlo simulations to define critical values for dissimilarity coefficients for use in MAT reconstructions.

Usage

mcarlo(object, ...)

## Default S3 method:
mcarlo(object, nsamp = 10000,
       type = c("paired", "complete", "bootstrap", "permuted"),
       replace = FALSE, 
       method = c("euclidean", "SQeuclidean", "chord", "SQchord",
                  "bray", "chi.square", "SQchi.square",
                  "information", "chi.distance", "manhattan",
                  "kendall", "gower", "alt.gower", "mixed"),
       is.dcmat = FALSE, diag = FALSE, ...)

## S3 method for class 'mat'
mcarlo(object, nsamp = 10000,
       type = c("paired", "complete", "bootstrap", "permuted"),
       replace = FALSE, diag = FALSE, ...)

## S3 method for class 'analog'
mcarlo(object, nsamp = 10000,
       type = c("paired", "complete", "bootstrap", "permuted"),
       replace = FALSE, diag = FALSE, ...)

Arguments

Details

Only "type" "paired" and "bootstrap" are currently implemented.

distance produces square, symmetric dissimilarity matrices for training sets. The upper triangle of these matrices is a duplicate of the lower triangle, and as such is redundant. mcarlo works on the lower triangle of these dissimilarity matrices, representing all pairwise dissimilarity values for training set samples. The default is not to include the diagonal (zero) values of the dissimilarity matrix. If you feel that these diagonal (zero) values are part of the population of dissimilarities then use "diag = TRUE" to include them in the permutations.

Value

A vector of simulated dissimilarities of length "nsamp". The "method" used is stored in attribute "method".

Note

The performance of these permutation and simulation techniques still needs to be studied. This function is provided for pedagogic reasons. Although recommended by Sawada et al (2004), sampling with replacement ("replace = TRUE") and including diagonal (zero) values ("diag = TRUE") simulates too many zero distances. This is because the same training set sample can, on occasion be drawn twice leading to a zero distance. It is impossible to find in nature two samples that will be perfectly similar, and as such sampling with replacement and "diag = TRUE" seems undesirable at best.

Author(s)

Gavin L. Simpson

References

Sawada, M., Viau, A.E., Vettoretti, G., Peltier, W.R. and Gajewski, K. (2004) Comparison of North-American pollen-based temperature and global lake-status with CCCma AGCM2 output at 6 ka. Quaternary Science Reviews 23, 87–108.

See Also

mat for fitting MAT models and analog for analogue matching. roc as an alternative method for determining critical values for dissimilarity measures when one has grouped data.

plot.mcarlo provides a plotting method to visualise the distribution of simulated dissimilarities.

Examples

## Imbrie and Kipp example
## load the example data
data(ImbrieKipp)
data(SumSST)
data(V12.122)

## merge training and test set on columns
dat <- join(ImbrieKipp, V12.122, verbose = TRUE)

## extract the merged data sets and convert to proportions
ImbrieKipp <- dat[[1]] / 100
V12.122 <- dat[[2]] / 100

## perform the modified method of Sawada (2004) - paired sampling,
## with replacement
ik.mcarlo <- mcarlo(ImbrieKipp, method = "chord", nsamp = 1000,
                    type = "paired", replace = FALSE)
ik.mcarlo

## plot the simulated distribution
layout(matrix(1:2, ncol = 1))
plot(ik.mcarlo)
layout(1)

Description

Minimum dissimilarity is a useful indicator of reliability of reconstructions performed via MAT and other methods, and for analogue matching. Minimum dissimilarity for a sample is the smallest dissimilarity between it and the training set samples.

Usage

minDC(x, ...)

## Default S3 method:
minDC(x, ...)

## S3 method for class 'predict.mat'
minDC(x, ...)

## S3 method for class 'analog'
minDC(x, probs = c(0.01, 0.02, 0.05, 0.1), ...)

## S3 method for class 'wa'
minDC(x, y,
      method = c("euclidean", "SQeuclidean", "chord", "SQchord",
                 "bray", "chi.square", "SQchi.square", "information",
                 "chi.distance", "manhattan", "kendall", "gower",
                 "alt.gower", "mixed"),
      percent = FALSE, probs = c(0.01, 0.025, 0.05, 0.1), ...)

Arguments

Value

minDC returns an object of class "minDC".

An object of class minDC is a list with some or all of the following components:

Note

The "default" method of minDC will attempt to extract the relevant component of the object in question. This may be useful until a specific minDC method is written for a given class.

Author(s)

Gavin L. Simpson

See Also

predict.mat, and plot.minDC for a plotting method.

Examples

## Imbrie and Kipp example
## load the example data
data(ImbrieKipp)
data(SumSST)
data(V12.122)

## merge training and test set on columns
dat <- join(ImbrieKipp, V12.122, verbose = TRUE)

## extract the merged data sets and convert to proportions
ImbrieKipp <- dat[[1]] / 100
V12.122 <- dat[[2]] / 100

## fit the MAT model using the squared chord distance measure
ik.mat <- mat(ImbrieKipp, SumSST, method = "SQchord")
ik.mat

## reconstruct for the V12-122 core data
v12.mat <- predict(ik.mat, V12.122)

## extract the minimum DC values
v12.mdc <- minDC(v12.mat)
v12.mdc

## draw a plot of minimum DC by time
plot(v12.mdc, use.labels = TRUE, xlab = "Depth (cm.)")

Description

Hills N2 is a measure of species diversity, commonly referred to as "effective" diversity. If computed on the rows (samples) then the "effective" number of species in each row is returned, whereas, if computed on the columns (species) then the "effective" number of occurences of each species in the data set is returned.

Usage

n2(x, ...)

## Default S3 method:
n2(x, which = c("species", "sites"), ...)

Arguments

Value

A numeric vector of N2 values.

Author(s)

Gavin L. Simpson

Examples

data(swapdiat)
sppN2 <- n2(swapdiat, "species")
head(sppN2)
sampN2 <- n2(swapdiat, "sites")
head(sampN2)

Description

Computes weighted average optima and tolerance ranges from species abundances and values of the environment.

Usage

optima(x, ...)

## Default S3 method:
optima(x, env, boot = FALSE, nboot = 1000,
       alpha = 0.05, ...)

## Default S3 method:
tolerance(x, env, useN2 = TRUE, ...)

Arguments

Value

Both functions return a named vector containing the WA optima or tolerances for the environmental gradient specified by env.

Note

Objects of class "optima" or "tolerance" can be coerced to data frames using methods for as.data.frame.

Author(s)

Gavin L. Simpson

See Also

wa

Examples

## Load the Imbrie & Kipp data and
## summer sea-surface temperatures
data(ImbrieKipp)
data(SumSST)

## WA optima
(opt <- optima(ImbrieKipp, SumSST))

## WA tolerances
(tol <- tolerance(ImbrieKipp, SumSST, useN2 = TRUE))

## caterpillar plot
caterpillarPlot(opt, tol)

## convert to data frame
as.data.frame(opt)
as.data.frame(tol)

## bootstrap WA optima - 100 resamples too low for SD & pCI
bopt <- optima(ImbrieKipp, SumSST, boot = TRUE, nboot = 100)
head(bopt)

Description

A modified version of panel.loess, for drawing Loess smooths on stratigraphic diagrams.

Usage

panel.Loess(x, y,
            span = 1/3, degree = 1,
            family = c("symmetric","gaussian"),
            evaluation = 50,
            lwd = plot.line$lwd,
            lty = plot.line$lty,
            col,
            col.line = plot.line$col,
            type,
            ...)

Arguments

Details

The standard panel function panel.loess treats the data as the x-axis acting as the time component. For stratigraphic plots where time flows along the y-axis, we want the smoother to be fitted with the x-axis data as the response and the time component (y-axis) as the predictor.

This modified version of panel.loess flips the two axes to produce the desired effect. Note also that it does not have argument horizontal as this is not required or supported by Stratiplot. In other respects, panel.Loess is equivalent to the lattice panel function panel.loess.

User should note that warnings can be generated by the fitting function if span is set too small for species with few observations. In such cases, the user is directed to the help page for loess.smooth, but increasing span slightly can often stop the warnings.

Author(s)

Gavin L. Simpson, slightly modified from the Lattice function panel.loess by Deepayan Sarkar.

See Also

loess.smooth, panel.loess.

Description

A Lattice panel function for drawing individuals panels on stratigraphic diagrams using a range of plot types commonly used within palaeoecology.

Usage

panel.Stratiplot(x, y,
                 type = "l",
                 col,
                 pch = plot.symbol$pch,
                 cex = plot.symbol$cex,
                 col.line = plot.line$col,
                 col.symbol = plot.symbol$col,
                 col.refline = ref.line$col,
                 col.smooth = "red",
                 col.poly = plot.line$col,
                 lty = plot.line$lty,
                 lwd = plot.line$lwd,
                 lty.smooth = plot.line$lty,
                 lwd.smooth = 2,
                 lwd.h = 3,
                 fill = plot.symbol$fill,
                 zones = NULL,
                 col.zones = plot.line$col,
                 lty.zones = plot.line$lty,
                 lwd.zones = plot.line$lwd,
                 gridh = -1, gridv = -1,
                 ...)

Arguments

Details

Creates stratigraphic scatter plots of x and y, with various modifications possible via the type argument.

Zones can be drawn on the panels by supplying the numeric vector of zone boundaries as argument zones. The panel function will then draw horizontal lines across the panels at the desired y-axis locations. Note that the panel function does not attempt to identify the zone boundaries automatically; these must be determined via a chronological (constrained) cluster analysis function or similar.

Note that all the arguments controlling the display can be supplied directly to a high-level call of the function Stratiplot.

Note

Histogram-like bars (type = "h") are drawn with lineend = "butt" to improve their appearance. This can not be changed by the user and you can't include that grid parameter in any call to panel.Stratiplot that uses type = "h".

Author(s)

Gavin L. Simpson

See Also

Stratiplot, panel.Loess, panel.xyplot.

Description

Fits a palaeoecological transfer function model using principal component regression, using an optional transformation of the matrix of predictor variables when these are species abundance data.

Usage

## Default S3 method:
pcr(x, y, ncomp, tranFun, ...)

## S3 method for class 'formula'
pcr(formula, data, subset, na.action, ..., model = FALSE)

Hellinger(x, ...)

ChiSquare(x, apply = FALSE, parms)

## S3 method for class 'pcr'
performance(object, ...)

## S3 method for class 'pcr'
residuals(object, comps = NULL, ...)

## S3 method for class 'pcr'
fitted(object, comps = NULL, ...)

## S3 method for class 'pcr'
coef(object, comps = NULL, ...)

## S3 method for class 'pcr'
screeplot(x, restrict = NULL,
          display = c("RMSE","avgBias","maxBias","R2"),
          xlab = NULL, ylab = NULL, main = NULL, sub = NULL, ...)

## S3 method for class 'pcr'
eigenvals(x, ...)

Arguments

Details

When applying cross-validation (CV) to transfer function models, any transformation of the predictors must be applied separately during each iteration of the CV procedure to the part of the data used in fitting the model. In the same way, any samples to be predicted from the model must use any meta-parameters derived from the training data only. For examle, centring is appled to the training data only and the variables means used to centre the training data are used to centre the test samples. The variable means should not be computed on a combination of the training and test samples.

When using PCR, we might wish to apply a transformation to the species data predictor variables such that the PCA of those data preserves a dissimilarity coefficient other than the Euclidean distance. This transformation is applied to allow PCA to better describe patterns in the species data (Legendre & Gallagher 2001).

How this is handled in pcr is to take a user-supplied function that takes a single argument, the matrix of predictor variables. The function should return a matrix of the same dimension as the input. If any meta-parameters are required for subsequent use in prediction, these should be returned as attribute "parms", attached to the matrix.

Two example transformation functions are provided implementing the Hellinger and Chi Square transformations of Legendre & Gallagher (2001). Users can base their transformation functions on these. ChiSquare() illustrates how meta-parameters should be returned as the attribute "parms".

Value

Returns an object of class "pcr", a list with the following components: