Package 'dplR'

Description

This package contains functions for performing some standard tree-ring analyses.

Details

Main Functions

read.rwl

reads rwl files

detrend

detrends raw ring widths

chron

builds chronologies

corr.rwl.seg

crossdating function

Author(s)

Andy Bunn andy.bunn@wwu.edu with major additions from Mikko Korpela and other significant contributions from Franco Biondi, Filipe Campelo, Pierre Mérian, Fares Qeadan and Christian Zang. Function redfit is an improved translation of program REDFIT which is original work of Manfred Mudelsee and Michael Schulz. Jacob Cecile contributed a bug fix to detrend.series. Allan Buras came up with the revised formula of glk in dplR >= 1.6.1.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Fritts, H. C. (2001) Tree Rings and Climate. Blackburn. ISBN-13: 978-1-930665-39-2.

Description

Applies an age-dependent smoothing spline to y.

Usage

ads(y, nyrs0 = 50, pos.slope = TRUE)

Arguments

Details

This implements the age-dependent smoothing spline similar to that described by Melvin (2004). In this implementation a cubic smoothing spline (caps) is fit to y with an initial stiffness of nyrs0. For each succesive measurement, the siffness is incremented by that ring index. This results in a spline is nyrs0 flexible at the start of the series and grows progressively stiffer. This can help capture the initial fast growth of a juvinielle tree with a flexible spline that then progresses to a stiffer spline that can better model the constant growth commonly found in mature trees. In its details, the cubic smoothing spline follows the Cook and Peters (1981) spline with a 50% frequency cutoff. See Cook and Kairiukstis (1990) for more information.

The default setting for nyrs0 is 50 years which is approprite for trees with a classic growth model but a value of 10 or 20 might be more appropriate for a Hugershoff-like initial increase in growth. Cook (pers comm) suggests a value of 20 for RCS.

If pos.slope is FALSE, the function will attempt to prevent a positive slope at the end of the series. In some cases when ads is used for detrending, a positive slope can be considered biologically unlikely. In those cases, the user can constrain the positive slope in the spline. This works by calculating the spline and taking the first difference. Then the function finds the last (outside) index where the spline changes slope and fixes the spline values from the that point to the end. Finally the spline is rerun along this constrained curve. See examples for details. The wisdom of constraining the slope in this manner depends very much on expert knowledge of the system.

Value

A filtered vector.

Author(s)

Fortran code provided by Ed Cook. Ported and adapted for dplR by Andy Bunn.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Melvin, T. M. (2004) Historical Growth Rates and Changing Climatic Sensitivity of Boreal Conifers. PhD Thesis, Climatic Research Unit, School of Environmental Sciences, University of East Anglia.

Cook, E. R. and Peters, K. (1981) The Smoothing Spline: A New Approach to Standardizing Forest Interior Tree-Ring Width Series for Dendroclimatic Studies. Tree-Ring Bulletin, 41, 45-53.

See Also

caps, detrend

Examples

# fit a curve
data(co021)
aSeries <- na.omit(co021$`641114`)
plot(aSeries,type="l",col="grey50")
lines(ads(y = aSeries),col="blue",lwd=2)

# show an artificial series with a Hugershoff-like curve.
a <- 0.5
b <- 1
g <- 0.1
d <- 0.25

n <- 300
x <- 1:n
y <- I(a*x^b*exp(-g*x)+d)
# add some noise
y <- y + runif(n=length(y),min = 0,max = 0.5)
# Plot with two different splines.
plot(y,type="l",col="grey50")
lines(ads(y,50),col="darkgreen",lwd=2) # bad
lines(ads(y,10),col="darkblue",lwd=2) # good

# now repeat with a positive slope to constrain
y <- I(a*x^b*exp(-g*x)+d)
y[251:300] <- y[251:300] + seq(0,0.25,length.out=50)
y <- y + runif(n=length(y),min = 0,max = 0.5)
plot(y,type="l",col="grey50")
lines(ads(y,10),col="darkgreen",lwd=2) #bad?
lines(ads(y,10,pos.slope=FALSE),col="darkgreen",lwd=2,lty="dashed")

Description

This data set gives the raw ring widths for Norway spruce Picea abies at Rothenburg ob der Tauber, Bavaria, Germany. There are 20 series from 10 trees. Data set was created using read.rwl and saved to an .rda file using save.

Usage

data(anos1)

Format

A data.frame containing 20 tree-ring series from 10 trees in columns and 98 years in rows. The correct stc mask for use with read.ids is c(5, 2, 1).

References

Zang, C. (2010) Growth reaction of temperate forest tree species to summer drought – a multispecies tree-ring network approach. Ph.D. thesis, Technische Universität München.

Description

Class "bai" holds basal area increment: the area of each ring, with one series per column and one year per row, as made by bai.in and bai.out. as.bai labels a data.frame or matrix of ring areas as one.

Usage

as.bai(x)

## S3 method for class 'bai'
x[i, j, drop]

## S3 method for class 'bai'
subset(x, subset, select, drop = FALSE, ...)

## S3 method for class 'bai'
time(x, ...)

## S3 method for class 'bai'
summary(object, ...)

## S3 method for class 'bai'
plot(x, plot.type = c("seg", "spag"), ...)

Arguments

Details

A "bai" object has the same shape as an "rwl" object and makes the same promise: the row names are the years, consecutive and increasing. It inherits from neither "rwl" nor "rwi", so that areas, widths and indices can be told apart.

dplR's functions use the class to check what they are given. Those that want ring-width indices, such as chron and rwi.stats, take a "bai" object without a word: a mean of basal area increment by year is a BAI chronology. So do the crossdating functions, the plots, common.interval, and detrend, since fitting a curve to basal area increment is established practice. Functions that want ring widths, such as bai.in, rcs, cms and rwl.report, warn and treat the areas as widths. If the values really are widths, relabel them with as.rwl.

Two attributes travel with it. attr(x, "dplR.bai") says how the areas were made: the function (fun) and whether the distance to pith (d2pith, for bai.in) or the diameters (diam, for bai.out) were given rather than estimated from the widths. attr(x, "dplR.provenance") is the provenance record of the ring widths the areas were made from, if they had one (see read.tucson). Neither is set by as.bai, which only labels the object.

Subsetting works as for "rwl" objects ([.rwl), and window takes a span of years (window.rwl). summary gives the per-series statistics of rwl.stats, in the units of the areas.

Value

as.bai: an object of class c("bai", "data.frame"). x is returned unchanged if it already has that class. time: a numeric vector of years. summary: a data.frame, as from rwl.stats.

Author(s)

Andy Bunn

See Also

bai.in, bai.out, as.rwl, as.rwi, chron

Examples

library(utils)
data(gp.rwl)
data(gp.d2pith)
gp.bai <- bai.in(rwl = gp.rwl, d2pith = gp.d2pith)
class(gp.bai)
str(attr(gp.bai, "dplR.bai"))
head(summary(gp.bai))
## A BAI chronology
gp.bai.crn <- chron(gp.bai)
plot(gp.bai.crn, add.spline = TRUE, nyrs = 32)
plot(gp.bai[, 1:5])

Description

Class "rwi" holds ring-width indices: detrended series with one series per column and one year per row, as made by detrend, rcs and cms. as.rwi labels a data.frame or matrix of indices as one.

Usage

as.rwi(x)

## S3 method for class 'rwi'
x[i, j, drop]

## S3 method for class 'rwi'
subset(x, subset, select, drop = FALSE, ...)

## S3 method for class 'rwi'
time(x, ...)

## S3 method for class 'rwi'
summary(object, ids = NULL, pcrit = 0.05, ...)

## S3 method for class 'summary.rwi'
print(x, max.print = 10, ...)

## S3 method for class 'summary.rwi'
as.data.frame(x, ...)

## S3 method for class 'rwi'
plot(x, plot.type = c("spag", "seg", "image"), ...)

Arguments

Details

An "rwi" object has the same shape as an "rwl" object and makes the same promise: the row names are the years, consecutive and increasing. It does not inherit from "rwl", so that indices and ring widths can be told apart.

dplR's functions use the class to check what they are given. Those that want ring widths – detrend, rcs, cms, i.detrend, bai.in, bai.out, pointer, strip.rwl, ssf and rwl.report – warn when given an "rwi" object, and say what goes wrong: detrend on indices, for instance, divides out a growth curve that has already been removed. Those that want indices – chron, chron.ars, chron.stabilized, rwi.stats, rwi.stats.running and sss – warn when given an "rwl" object: rwi.stats(ca533) gives an rbar.eff of 0.350 from the widths against 0.423 from the Spline indices, and nothing about the number says it is wrong. They still take a plain data.frame or matrix without a word, and a "bai" object of basal area increment (see as.bai). The crossdating functions, the plots, common.interval, sgc and rwl.stats take either class quietly. Each warning is only a warning, and the function goes on: if the class is what is wrong – indices read from a file with read.rwl come back as class "rwl" – relabel them with as.rwi, or relabel widths with as.rwl.

Two attributes travel with it. attr(x, "dplR.detrend") is a list that says how the indices were made: the function (fun), the method and its settings, and always difference, which is TRUE when the indices are differences from the fitted curve (centred on 0) rather than ratios (centred on 1). attr(x, "dplR.provenance") is the provenance record of the ring widths the indices were made from, if they had one (see read.tucson). Neither is set by as.rwi, which only labels the object.

Subsetting works as for "rwl" objects ([.rwl): the class and both records are kept, dropping series trims the years that none of the remaining series cover, and a row subset that leaves the years out of order or with holes returns a plain data.frame, with a warning. subset and window.rwi drop series left with no values and name them in a message; [ keeps them, on purpose, so that its columns stay lined up with anything indexed alongside them.

summary describes the indices as a collection. It gives how they were made (from attr(x, "dplR.detrend")), their span and common interval, the collection statistics from rwi.stats (including rbar.eff, EPS and SNR), and for each series its first and last years, length, mean, standard deviation, first-order autocorrelation, and its correlation with the mean of the other series, with the p-value, from interseries.cor with its defaults. Printing it shows the collection and lists the series that do not correlate with the others at pcrit; as.data.frame gives the per-series table. The collection statistics count each series as its own tree unless ids is given, and the print says which. The print ends by pointing to summary with ids (when it was not given), rwi.stats.running (rbar and EPS through time, which usually fall off where sample depth does) and corr.rwl.seg (where in a series the fit breaks down). A series with no values at all is listed and left out of the collection statistics, the correlations and the common interval. With fewer than two series with values there is no collection, and the statistics and correlations are NULL and NA.

plot draws, by default, spag.plot, which draws indices around 1 (or 0 for differences) and not around each series' mean, so the grey line under each series is the value the indices should sit at. plot.type = "image" draws each index as a coloured cell, with years across and series down, earliest-starting at the bottom. Colours are brown below 1 (or 0) and green above it. Each side is scaled separately and clipped at the clip quantile of its departures, since ratio indices cannot fall below 0 but can run well above 2; the key gives the clips. A run of one colour at the start or end of a series is a growth trend the detrending did not remove, and a vertical stripe is a year the series agree on.

Arithmetic on an "rwi" object (e.g. x + 1) returns a plain data.frame, as it does for any data.frame subclass. Use as.rwi on the result if it is still a set of indices.

Value

as.rwi: an object of class c("rwi", "data.frame"). x is returned unchanged if it already has that class. time: a numeric vector of years. summary: an object of class "summary.rwi", a list with how (the dplR.detrend record), n.series, first, last, common (the first and last years in which every series has a value, or NA), stats (from rwi.stats), ids.given, empty (the names of series with no values), series (the per-series data.frame) and pcrit.

Author(s)

Andy Bunn

See Also

detrend, rcs, cms, as.rwl, chron, rwi.stats

Examples

library(utils)
data(ca533)
ca533.rwi <- detrend(ca533, method = "Spline")
class(ca533.rwi)
str(attr(ca533.rwi, "dplR.detrend"))
summary(ca533.rwi)
## Group cores by tree, from the series IDs
summary(ca533.rwi, ids = autoread.ids(ca533))
head(as.data.frame(summary(ca533.rwi)))
plot(ca533.rwi[, 1:5])
data(co021)
## The growth trend "Mean" leaves in: a run of green at the start of
## nearly every series.
plot(detrend(co021, method = "Mean"), plot.type = "image")

Description

Attempts to turn its argument into a rwl object.

Usage

as.rwl(x)

Arguments

Details

This tries to coerce x into class c("rwl","data.frame"). Failable.

Value

An object of class c("rwl", "data.frame") with the series in columns and the years as rows. The series IDs are the column names and the years are the row names.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

Examples

library(graphics)
library(stats)
library(utils)
## Toy
n <- 100
## Make a data.frame that is tree-ring like
base.series <- 0.75 + exp(-0.2 * 1:n)
foo <- data.frame(x1 = base.series + abs(rnorm(n, 0, 0.25)),
                  x2 = base.series + abs(rnorm(n, 0, 0.25)),
                  x3 = base.series + abs(rnorm(n, 0, 0.25)),
                  x4 = base.series + abs(rnorm(n, 0, 0.25)),
                  x5 = base.series + abs(rnorm(n, 0, 0.25)),
                  x6 = base.series + abs(rnorm(n, 0, 0.25)))
# coerce to rwl and use plot and summary methods
foo <- as.rwl(foo)
class(foo)
plot(foo, plot.type="spag")
summary(foo)

Description

Convert multiple ring-width series to basal area increment (i.e., ring area) going from the pith to the bark.

Usage

bai.in(rwl, d2pith = NULL)

Arguments

Details

This converts ring-width series (mm) to ring-area series (mm squared) (aka basal area increments) based on the distance between the innermost measured ring and the pith of the tree. It is related to bai.out, which calculates each ring area starting from the outside of the tree and working inward. Both methods assume a circular cross section (Biondi 1999). See the references below for further details.

Value

An object of class c("bai", "data.frame") containing the ring areas for each series with the column names, row names and dimensions of rwl. It carries the provenance record of rwl, if it has one, and records how it was made in attr(x, "dplR.bai"): the function (fun) and whether d2pith was given (d2pith). See as.bai for which functions take it.

Note

DendroLab website: https://dendrolaborg.wordpress.com/

Author(s)

Code by Andy Bunn based on work from DendroLab, University of Nevada Reno, USA. Patched and improved by Mikko Korpela.

References

Biondi, F. (1999) Comparing tree-ring chronologies and repeated timber inventories as forest monitoring tools. Ecological Applications, 9(1), 216–227.

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

See Also

bai.out, as.bai

Examples

library(graphics)
library(stats)
library(utils)
## Toy
n <- 100
## Make three fake tree-ring series to show that these funcs work on rwl objects
base.series <- 0.75 + exp(-0.2 * 1:n)
rwl <- data.frame(x1 = base.series + abs(rnorm(n, 0, 0.05)),
                  x2 = base.series + abs(rnorm(n, 0, 0.05)),
                  x3 = base.series + abs(rnorm(n, 0, 0.05)))

## The inside out method
foo <- bai.in(rwl = rwl)
## The outside in method
bar <- bai.out(rwl = rwl)

## Identical
head(bar)
head(foo)

## Use gp data
data(gp.rwl)
data(gp.d2pith)
foo <- bai.in(rwl = gp.rwl, d2pith = gp.d2pith)
foo.crn <- chron(foo)
yrs <- time(foo.crn)
plot(yrs, foo.crn[, 1], type = "n",
     xlab = "Year", ylab = expression(mm^2))
lines(yrs, foo.crn[, 1], col = "grey", lty = "dashed")
lines(yrs, caps(foo.crn[, 1], nyrs = 32), col = "red", lwd = 2)

Description

Convert multiple ring-width series to basal area increment (i.e., ring area) going from the bark to the pith.

Usage

bai.out(rwl, diam = NULL)

Arguments

Details

This converts ring-width series (mm) to ring-area series (mm squared) (aka basal area increments) based on the diameter of the tree and the width of each ring moving towards the pith of the tree. It is related to bai.in, which calculates each ring area starting from the inside of the tree and working outward. Both methods assume a circular cross section (Biondi 1999). See the references below for further details.

Value

An object of class c("bai", "data.frame") containing the ring areas for each series with the column names, row names and dimensions of rwl. It carries the provenance record of rwl, if it has one, and records how it was made in attr(x, "dplR.bai"): the function (fun) and whether diam was given (diam). See as.bai for which functions take it.

Note

DendroLab website: https://dendrolaborg.wordpress.com/

Author(s)

Code by Andy Bunn based on work from DendroLab, University of Nevada Reno, USA. Patched and improved by Mikko Korpela.

References

Biondi, F. (1999) Comparing tree-ring chronologies and repeated timber inventories as forest monitoring tools. Ecological Applications, 9(1), 216–227.

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

See Also

bai.in, as.bai

Examples

library(graphics)
library(utils)
## Not run: 
library(stats)
## Toy
n <- 100
## Make three fake tree-ring series to show that these funcs work on rwl objects
base.series <- 0.75 + exp(-0.2 * 1:n)
rwl <- data.frame(x1 = base.series + abs(rnorm(n, 0, 0.05)),
                  x2 = base.series + abs(rnorm(n, 0, 0.05)),
                  x3 = base.series + abs(rnorm(n, 0, 0.05)))

## The inside out method
foo <- bai.in(rwl = rwl)
## The outside in method
bar <- bai.out(rwl = rwl)

## Identical
head(bar)
head(foo)

## End(Not run)
## Use gp data
data(gp.rwl)
data(gp.dbh)
## dbh (minus the bark) from cm to mm 
gp.dbh2 <- gp.dbh[, 1:2]
gp.dbh2[, 2] <- (gp.dbh[, 2] - gp.dbh[, 3]) * 10
bar <- bai.out(rwl = gp.rwl, diam = gp.dbh2)
bar.crn <- chron(bar)
yrs <- time(bar.crn)
plot(yrs, bar.crn[, 1], type = "n",
     xlab = "Year", ylab = expression(mm^2))
lines(yrs, bar.crn[, 1], col = "grey", lty = "dashed")
lines(yrs, caps(bar.crn[, 1], nyrs = 32), col = "red", lwd = 2)

Description

Convert multiple ring-width series to basal area increment (i.e., ring area) following the proportional method of Bakker (2005).

Usage

bakker(rwl, ancillary)

Arguments

Details

This converts ring-width series (mm) to ring-area series (mm squared) (aka basal area increments) based on the diameter of the tree, the missing distance to the pith and the missing number of rings to the pith, following the proportional method for reconstructing historical tree diameters by Bakker (2005). It prevents bai.out transformations from producing negative increments when the sum of all ring widths in a series is larger than DBH/2. It prevents bai.in transformations from producing too small values when the sum of all ring widths in a series is smaller than DBH/2.

Value

A list containing the following objects:

Author(s)

Code by Stefan Klesse. Adapted for dplR by Andy Bunn.

References

Bakker, J.D., 2005. A new, proportional method for reconstructing historical tree diameters. Canadian Journal of Forest Research 35, 2515–2520. https://doi.org/10.1139/x05-136

Examples

data(zof.rwl)
data(zof.anc)
zof.bakker <- bakker(rwl = zof.rwl,ancillary = zof.anc)
zof.bai <- zof.bakker$baiBakker_raw

# first series bai
yrs <- time(zof.rwl)
plot(yrs,zof.bai[,1],type="l",
     xlab="Year",
     ylab=expression(BAI~(mm^2)),
     main = colnames(zof.bai)[1])

Description

This data set gives the raw ring widths for bristlecone pine Pinus longaeva at Campito Mountain in California, USA. There are 34 series. Data set was created using read.rwl and saved to an .rda file using save.

Usage

data(ca533)

Format

A data.frame containing 34 tree-ring series in columns and 1358 years in rows.

Source

International tree-ring data bank, Accessed on 20-April-2021 at https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/usa/ca533.rwl

References

Graybill, D. A. and LaMarche, Jr., V. C. (1983) Campito Mountain Data Set. IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series 1983-CA533.RWL. NOAA/NCDC Paleoclimatology Program, Boulder, Colorado, USA.

Examples

library(utils)
data(ca533)

## Where the data came from. read.tucson() records what it saw and attaches
## it to what it returns; this data set carries that record.
prov <- attr(ca533, "dplR.provenance")
prov$file          # the archived file it was read from
prov$precision     # the precision each series was measured at
prov$gaps          # interior gaps, if any -- this file has none
prov$events        # problems found while reading -- none here either

## rwl.report() puts the file and precision at the head of its report
rwl.report(ca533)

Description

This data set gives the standard chronology for white spruce Picea glauca at Twisted Tree Heartrot Hill in Yukon, Canada. Data set was created using read.crn and saved to an .rda file using save.

Usage

data(cana157)

Format

A data.frame containing the standard chronology in column one and the sample depth in column two. There are 463 years (1530–1992) in the rows.

Source

International tree-ring data bank, Accessed on 20-April-2021 at https://www.ncei.noaa.gov/pub/data/paleo/treering/chronologies/northamerica/canada/cana157.crn

References

Jacoby, G., D’Arrigo, R. and Buckley, B. (1992) Twisted Tree Heartrot Hill Data Set. IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series 1992-CANA157.CRN. NOAA/NCDC Paleoclimatology Program, Boulder, Colorado, USA.

Description

Applies a smoothing spline to y with rigidity determined by two parameters: frequency response f at a wavelength of nyrs years.

Usage

caps(y, nyrs = 32, f = 0.5)

Arguments

Details

This applies the classic smoothing spline from Cook and Peters (1981). The rigidity of the spline has a frequency response of 50% at a wavelength of nyrs. The references, of course, have more information.

This function was introduced to dplR in version 1.7.3 and replaces the now defunct ffcsaps. Where ffcsaps was written entirely in R, caps is a wrapper for a Fortran subroutine from Ed Cook's ARSTAN program that is hundreds of times faster.

The default value of nyrs was changed to 32 from 2/3 length of y in version 1.7.8 based on a suggestion from Klesse (2021).

Note: nyrs is passed to the Fortran subroutine as an integer, so a fractional nyrs is truncated. This matters when nyrs is given as a proportion of series length, which is rarely a whole number. The effect on the fitted curve is small but not zero; see the “Mathematical Details” vignette for a worked comparison.

Note: caps stops if there are any NA values in y. A smoothing spline through a series with missing values is not defined, and the caller has to decide what to do about the gap rather than have one chosen silently. Before version 1.8.0 such a call returned a vector of NA, which looked like an answer and was not. See examples, and fill.internal.NA.

Value

A filtered vector.

Note

The mathematical relationship between nyrs, f and the spline's smoothing parameter, together with an empirical check of the frequency response and a demonstration that caps reproduces the deprecated ffcsaps, is given in the vignette: vignette("math-dplR", package = "dplR").

Author(s)

Fortran code provided by Ed Cook and adapted for dplR by Andy Bunn.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Cook, E. R. and Peters, K. (1981) The Smoothing Spline: A New Approach to Standardizing Forest Interior Tree-Ring Width Series for Dendroclimatic Studies. Tree-Ring Bulletin, 41, 45-53.

Klesse, S. (2021) Critical Note on the Application of the "Two-Third"" Spline. Dendrochronologia Volume 65: 125786

See Also

ads

Examples

library(graphics)
library(utils)

## Use first series from the Mesa Verde data set
data(co021)
series <- co021[, 1]
series <- series[!is.na(series)]
plot(series, type = "l", ylab = "Ring Width (mm)", col = "grey")
lines(caps(series, nyrs = 10), col = "red", lwd = 2)
lines(caps(series, nyrs = 100), col = "green", lwd = 2)
# The default, a 32-year spline since version 1.7.8
lines(caps(series), col = "blue", lwd = 2)
legend("topright",
       c("Series", "nyrs=10", "nyrs=100", "Default nyrs (32)"),
       fill=c("grey", "red", "green", "blue"))


## caps() needs a series with no missing values. Note that the example
## above drops the NA padding before fitting. A gap inside the measured
## span has to be dealt with deliberately: drop the series, or close the
## gap with fill.internal.NA().
y <- c(NA, NA, rnorm(100))
try(caps(y))
## fit over the measured values instead
head(caps(y[!is.na(y)]))

Description

Computes cross-correlations between a tree-ring series and a master chronology built from a rwl object at user-specified lags and segments.

Usage

ccf.series.rwl(rwl, series, series.yrs = as.numeric(names(series)),
               seg.length = 50, bin.floor = 100, n = NULL,
               nyrs = NULL, prewhiten = TRUE, ar.order.max = NULL,
               biweight = TRUE, pcrit = 0.05,
               lag.max = 5, make.plot = TRUE,
               floor.plus1 = FALSE, series.x = FALSE, ...)

Arguments

Details

This function calculates the cross-correlation function between a tree-ring series and a master chronology built from rwl looking at correlations lagged positively and negatively using ccf at overlapping segments set by seg.length. For instance, with lag.max set to 5, cross-correlations would be calculated at for each segment with the master lagged at k = -5:5 years.

The cross correlations are calculated calling ccf as
ccf(x=master, y=series, lag.max=lag.max, plot=FALSE) if series.x is FALSE and as ccf(x=series, y=master, lag.max=lag.max, plot=FALSE) if series.x is TRUE. This argument was introduced in dplR version 1.7.0. Different users have different expectations about how missing or extra rings are notated. If series.x = FALSE the behavior will be like COFECHA where a missing ring in a series produces a negative lag in the plot rather than a positive lag.

Correlations are calculated for the first segment, then the second segment and so on. Correlations are only calculated for segments with complete overlap with the master chronology.

Each series (including those in the rwl object) is optionally detrended as the residuals from a hanning filter with weight n. The filter is not applied if n is NULL. Detrending can also be done via prewhitening where the residuals of an ar model are added to each series mean. This is the default. The master chronology is computed as the mean of the rwl object using tbrm if biweight is TRUE and rowMeans if not. Note that detrending typically changes the length of the series. E.g., a hanning filter will shorten the series on either end by floor(n/2). The prewhitening default will change the series length based on the ar model fit. The effects of detrending can be seen with series.rwl.plot.

Value

A list containing matrices ccf and bins. Matrix ccf contains the correlations between the series and the master chronology at the lags window given by lag.max. Matrix bins contains the years encapsulated by each bin.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

References

Bunn, A. G. (2010) Statistical and visual crossdating in R using the dplR library. Dendrochronologia, 28(4), 251–258.

See Also

corr.rwl.seg, corr.series.seg, skel.plot, series.rwl.plot

Examples

library(utils)
data(co021)
dat <- co021
## Create a missing ring: delete the 1500 ring of series 641143 and
## drop the original from the master.  Dated from the bark, every ring
## before 1500 now sits one year late.
flagged <- dat$"641143"
names(flagged) <- rownames(dat)
flagged <- delete.ring(flagged, year = 1500)
dat$"641143" <- NULL
ccf.100 <- ccf.series.rwl(rwl = dat, series = flagged, seg.length = 100)
## Not run: 
flagged2 <- co021$"641143"
names(flagged2) <- rownames(dat)
ccf.100.1 <- ccf.series.rwl(rwl = dat, seg.length = 100,
                            series = flagged2)
## Select series by name or column position
ccf.100.2 <- ccf.series.rwl(rwl = co021, seg.length = 100,
                            series = "641143")
ccf.100.3 <- ccf.series.rwl(rwl = co021, seg.length = 100,
                            series = which(colnames(co021) == "641143"))
identical(ccf.100.1, ccf.100.2) # TRUE
identical(ccf.100.2, ccf.100.3) # TRUE

## End(Not run)

Description

Checks that an object is a valid rwl and, if not, attempts coercion via as.rwl. Warns about internal NA values found within any series.

Usage

check.rwl(rwl, why = NULL, bai.ok = FALSE)

Arguments

Details

check.rwl is for functions that want ring widths. If rwl is an "rwi" object of ring-width indices (see as.rwi), a warning names the calling function, says that it wants widths, adds why, and suggests as.rwl if the values really are widths. The object is then relabelled as class "rwl" and checked as below. A "bai" object of basal area increment (see as.bai) is treated the same way, unless bai.ok is TRUE, when it is relabelled without a warning.

If rwl is not already of class "rwl", coercion is attempted using as.rwl. If coercion succeeds, a warning is issued. If coercion fails (e.g., because row names are not consecutive integers, columns are not numeric, or the input is not a data.frame or matrix), the function stops with a single informative message that describes both the class problem and the reason coercion failed.

After class validation, each series is checked for internal NA values, defined as NA values that are sandwiched between non-NA values within a series (as opposed to leading or trailing NAs, which are standard in rwl objects). If any are found, a warning names the affected series and suggests fill.internal.NA.

This function is called at the top of the dplR functions that want ring widths, providing uniform validation across the package. Functions that take widths or indices alike, and those that want indices, use internal checks of their own (see as.rwi).

Value

An object of class c("rwl", "data.frame") with series in columns and years as row names. The series IDs are the column names. The returned object is identical to the input when the input is already a valid rwl.

Author(s)

Andy Bunn

See Also

as.rwl, as.rwi, read.rwl, fill.internal.NA

Examples

library(utils)
data(ca533)

## A proper rwl passes silently and is returned unchanged
ca533.checked <- check.rwl(ca533)
identical(ca533, ca533.checked)

## A plain data.frame is coerced with a warning.
## suppressWarnings() is used here to keep R CMD check clean since the
## warning is expected. Remove it to see the coercion warning.
ca533.df <- ca533
class(ca533.df) <- "data.frame"
ca533.rwl <- suppressWarnings(check.rwl(ca533.df))
#ca533.rwl <- check.rwl(ca533.df)
class(ca533.rwl)

Description

This function builds a mean value chronology, typically from a data.frame of detrended ring widths as produced by detrend.

Usage

chron(rwi, biweight = TRUE, prewhiten = FALSE, ...)

Arguments

Details

This either averages the rows of the data.frame using a mean or a robust mean (the so-called standard chronology) or can do so from the residuals of an AR process (the residual chronology).

The residual chronology here is the mean of series that were each whitened with their own AR model, the order chosen for each series by ar. Each whitened series is free of autocorrelation but their mean need not be. Where much of the persistence belongs to the individual trees and is not shared among them, whitening each series removes more than the series have in common, and the residual chronology comes out negatively autocorrelated. For wa082 detrended with a spline its lag-one autocorrelation is about -0.2. It is worth checking, for example with acf.

The residual chronology from chron.ars is built differently and will not have the same values. It uses one AR order for all the series, taken from a pooled model, and whitens the mean a second time, which removes the autocorrelation that averaging leaves behind. Use it when the residual chronology from this function is not close to white.

Value

An object of of class crn and data.frame with the standard chronology, residual chronology (if prewhitening was performed), and the sample depth. The years are stored as row numbers.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Fritts, H. C. (2001) Tree Rings and Climate. Blackburn. ISBN-13: 978-1-930665-39-2.

See Also

read.rwl, detrend, chron.ars, ar, crn.plot

Examples

library(graphics)
library(utils)
data(ca533)
ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")
ca533.crn <- chron(ca533.rwi)
plot(ca533.crn,xlab="Year",ylab="RWI")
## With residual chron
ca533.crn2 <- chron(ca533.rwi, prewhiten = TRUE)
plot(ca533.crn2,xlab="Year",ylab="RWI")

Description

This function builds three varieties of the mean-value chronology, including the ARSTAN chronology, typically from a data.frame of detrended ring widths as produced by detrend.

Usage

chron.ars(rwi, biweight=TRUE, maxLag=10, firstAICmin=TRUE,
  verbose=TRUE, prewhitenMethod=c("ar.yw","arima.CSS-ML"))

Arguments

Details

This produces three mean-value chronologies: standard, residual, and ARSTAN. Users unfamiliar with the concept behind the ARSTAN method should look to Cook (1985) for background and inspiration.

The standard chronology is the (biweight) mean value across rows and identical to chron.

The residual chronology is the prewhitened chronology as described by Cook (1985) and uses multivariate autoregressive modeling to determine the order of the AR process. It is important to note that the residual chronology produced here is different from the simple residual chronology produced by chron, which returns the residuals of an AR process using a naive call to ar. In practice the results will be similar. For more on the residual chronology in this function, see pp. 153-154 in Cook's 1985 dissertation.

The ARSTAN chronology is the residual chronology with the pooled AR coefficients from the multivariate autoregressive modeling put back in. The residual chronology is whitened a second time after the whitened series are averaged, because whitening each series with its own coefficients can remove more persistence than the series share and leave their mean negatively autocorrelated. The pooled model is added to that twice-whitened chronology, so the ARSTAN chronology has missing values in the same years as the residual chronology. Through dplR 1.8.0 the pooled model was instead added to each whitened series before averaging, which skipped the second whitening and could give an ARSTAN chronology with negative autocorrelation where the persistence is mostly within trees. See references for details.

The order of the AR model is selected from the pooled AR coefficients by AIC using either the first local AIC minimum (firstAICmin = TRUE, the default) or the overall minimum across all lags up to maxLag (firstAICmin = FALSE). When firstAICmin = TRUE and the AIC decreases monotonically across all tested lags without reaching a local minimum, the function stops with an error. This can occur when maxLag is too small for the data or when the series have persistence structures that are not well-suited to this approach (e.g., non-stationary or very long-memory series). In such cases, consider increasing maxLag or reviewing the data preparation steps described in Cook (1985). When firstAICmin = FALSE the function always selects an order and does not stop in this way.

If the selected AR order is zero (no detectable common autocorrelation among the series), no prewhitening or post-reddening is applied and all three returned chronologies are identical to the standard chronology.

Once the AR order is determined an AR(p) model is fit to each series using either ar via the Yule-Walker method or by arima via conditional-sum-of-squares to find starting values, then maximum likelihood. It is possible that the model will not converge in which case a warning is produced. The AR fitting is determined via prewhitenMethod and defaults to using ar.

Value

A data.frame of class "crn" with the following columns:

Row names are the years taken from rownames(rwi).

Author(s)

Andy Bunn with contributions from Kevin Achukaitis and Ed Cook. Much of the function is a port of Cook's FORTRAN code.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Cook, E. R. (1985). A Time Series Analysis Approach to Tree Ring Standardization. PhD thesis, The University of Arizona.

See Also

chron, crn.plot, ar, arima

Examples

library(graphics)
library(utils)
data(co021)
co021.rwi <- detrend(rwl = co021, method = "AgeDepSpline")
co021.crn <- chron.ars(co021.rwi)
plot(co021.crn,xlab="Year",ylab="RWI",add.spline=TRUE,nyrs=20)
cor(co021.crn)

Description

This function builds a mean value chronology with bootstrapped confidence intervals, typically from a data.frame of detrended ring widths as produced by detrend.

Usage

chron.ci(x, biweight=TRUE, conf=0.95, R=100)

Arguments

Details

This either averages the rows of the data.frame using a mean or a robust mean (the so-called standard chronology) and calculates boostrapped confidence intervals using the normal approximation. The function will fail if there are any rows in x that contain only one sample and in practice there should be several samples in a row. One of the guiding principles of bootstrapping is that the population is to the sample as the sample is to the bootstrap samples.

Value

An object of of class data.frame with the standard chronology, the upper and lower confidence interval, and the sample depth. The years are stored as row numbers.

Author(s)

Andy Bunn.

See Also

read.rwl, detrend, boot, boot.ci

Examples

library(graphics)
library(utils)
data(wa082)
# truncate to a sample depth of five
wa082Trunc <- wa082[rowSums(!is.na(wa082))>4,]
# detrend
wa082RWI <- detrend(wa082Trunc, method="AgeDepSpline")
# bootstrap the chronology and
wa082Crn <- chron.ci(wa082RWI, biweight = TRUE, R = 100, conf = 0.99)
head(wa082Crn)
# plot (this is so much easier in ggplot!)
wa082Crn$yrs <- time(wa082Crn)
xx <- c(wa082Crn$yrs,rev(wa082Crn$yrs))
yy <- c(wa082Crn$lowerCI,rev(wa082Crn$upperCI))
plot(wa082Crn$yrs,wa082Crn$std,type="n",ylim=range(yy),
     ylab="RWI",xlab="Year",main="Chronology with CI")
polygon(x=xx,y=yy,col = "grey",border = NA)
lines(wa082Crn$yrs,wa082Crn$std)

Description

This function builds a variance stabilized mean-value chronology, typically from a data.frame of detrended ring widths as produced by detrend.

Usage

chron.stabilized(x, winLength, biweight = TRUE, running.rbar = FALSE)

Arguments

Details

The variance of a mean chronology depends on the variance of the individual samples, the number of series averaged together, and their interseries correlation (Wigley et al. 1984). As the number of series commonly decreases towards the beginning of a chronology averaging introduces changes in variance that are a solely an effect of changes in sample depth.

Additionally, time-dependent changes in interseries correlation can cause artificial variance changes of the final mean chronology. The function chron.stabilized accounts for both temporal changes in the interseries correlation and sample depth to produce a mean value chronology with stabilized variance.

The basic correction centers around the use of the effective independent sample size, Neff, which considers sample replication and mean interseries correlation between the samples at every time. This is defined as: Neff = n(t) / 1+(n(t)-1)rbar(t)

where n(t) is the number of series at time t, and rbar is the interseries correlation (see interseries.cor). Multiplication of the mean time series with the square root of Neff at every time t theoretically results in variance that is independent of sample size. In the limiting cases, when the rbar is zero or unity, Neff obtains values of the true sample size and unity, respectively.

Value

An object of of class crn and data.frame with the variance stabilized chronology, running interseries correlation ('if running.rbar=TRUE), and the sample depth.

Author(s)

Original code by David Frank and adapted for dplR by Stefan Klesse. Patched and improved by Andy Bunn.

References

Frank, D, Esper, J, Cook, E, (2006) On variance adjustments in tree-ring chronology development. Tree rings in archaeology, climatology and ecology, TRACE 4, 56–66

Frank, D, Esper, J, Cook, E, (2007) Adjustment for proxy number and coherence in a large-scale temperature reconstruction. Geophysical Research Letters 34

Wigley, T, Briffa K, Jones P (1984) On the Average Value of Correlated Time Series, with Applications in Dendroclimatology and Hydrometeorology. J. Climate Appl. Meteor., 23, 201–213

See Also

chron

Examples

library(graphics)
library(utils)
data(co021)
co021.rwi <- detrend(co021,method = "Spline")
co021.crn <- chron(co021.rwi)
co021.crn2 <- chron.stabilized(co021.rwi,
                                winLength=101,
                                biweight = TRUE,
                                running.rbar = FALSE)
yrs <- time(co021)
plot(yrs,co021.crn$std,type="l",col="grey")
lines(yrs,co021.crn2$adj.crn,col="red")

Description

Detrend multiple ring-width series simultaneously using the C-method.

Usage

cms(rwl, po, c.hat.t = FALSE, c.hat.i = FALSE)

Arguments

Details

This method detrends and standardizes tree-ring series by calculating a growth curve based on constant annual basal area increment. The method is based on the “assumption that constant growth is expressed by a constant basal area increment distributed over a growing surface” (Biondi and Qeadan 2008). The detrending is the estimation and removal of the tree’s natural biological growth trend. The standardization is done by dividing each series by the growth trend to produce units in the dimensionless ring-width index (RWI).

This attempts to remove the low frequency variability that is due to biological or stand effects.

A series with no values at all is dropped, with a message naming it, and the rest are standardized as if it had never been there. po must still have a row for every series in rwl, empty or not.

See the reference below for further details.

Value

An object of class c("rwi", "data.frame") (see as.rwi) containing the dimensionless and detrended ring-width indices with column names, row names and dimensions of rwl if c.hat.t is FALSE and c.hat.i is FALSE.

Otherwise a list of length 2 or 3 containing the RWI data.frame, a data.frame containing the C-curves for each tree (c.hat.t), and/or a vector containing the C-values for each tree (c.hat.i) depending on the output flags. See Eq. 12 in Biondi and Qeadan (2008) for more detail on c.hat.t, and c.hat.i.

Note

DendroLab website: https://dendrolaborg.wordpress.com/

Author(s)

Code provided by DendroLab based on programming by F. Qeadan and F. Biondi, University of Nevada Reno, USA and adapted for dplR by Andy Bunn. Patched and improved by Mikko Korpela.

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

See Also

detrend, chron, rcs

Examples

library(graphics)
library(utils)
data(gp.rwl)
data(gp.po)
gp.rwi <- cms(rwl = gp.rwl, po = gp.po)
gp.crn <- chron(gp.rwi)
crn.plot(gp.crn, add.spline = TRUE)
## c.hat
gp.rwi <- cms(rwl = gp.rwl, po = gp.po, c.hat.t = TRUE, c.hat.i = TRUE)
dotchart(gp.rwi$c.hat.i, ylab = "Series", xlab = expression(hat(c)[i]))
tmp <- gp.rwi$c.hat.t
plot(tmp[, 1], type = "n", ylim = range(tmp, na.rm = TRUE),
     xlab = "Cambial Age", ylab = expression(hat(c)[t]))
apply(tmp, 2, lines)

Description

This data set gives the raw ring widths for Douglas fir Pseudotsuga menziesii at Mesa Verde in Colorado, USA. There are 35 series. Data set was created using read.rwl and saved to an .rda file using save.

Usage

data(co021)

Format

A data.frame containing 35 tree-ring series in columns and 788 years in rows.

Source

International tree-ring data bank, Accessed on 20-April-2021 at https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/usa/co021.rwl

References

Schulman, E. (1963) Schulman Old Tree No. 1 Data Set. IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series 1983-CO021.RWL. NOAA/NCDC Paleoclimatology Program, Boulder, Colorado, USA.

Examples

library(utils)
data(co021)

## Where the data came from. read.tucson() records what it saw and attaches
## it to what it returns; this data set carries that record.
prov <- attr(co021, "dplR.provenance")
prov$file          # the archived file it was read from
prov$precision     # the precision each series was measured at
prov$gaps          # interior gaps, if any -- this file has none
prov$events        # problems found while reading -- none here either

## rwl.report() puts the file and precision at the head of its report
rwl.report(co021)

Description

This function combines any number of data.frames of tree-ring data into one data.frame.

Usage

combine.rwl(x, y = NULL)

Arguments

Details

The sequence of years in each data.frame must be increasing and continuous. The output produced by the function also fulfills this condition. If the input is differently formatted, the result will be wrong.

Value

An object of class c("rwl", "data.frame") with the series in columns and the years as rows. The keycodes are the column names and the years are the row names.

Author(s)

Christian Zang. Patched by Mikko Korpela.

Examples

library(utils)
data(ca533)
data(co021)
combi1 <- combine.rwl(list(ca533, co021))
## or alternatively for data.frames to combine
combi2 <- combine.rwl(ca533, co021)
identical(combi1, combi2) # TRUE

Description

This function finds the common interval on a set of tree-ring widths such as that produced by read.rwl.

Usage

common.interval(rwl, type=c("series", "years", "both"),
                make.plot=TRUE)

Arguments

Details

This trims an rwl object to a common interval that maximizes the number of series (type="series"), the number of years (type="years"), or a compromise between the two (type="both"). A modified seg.plot can be drawn as well.

Series with no values are left out. A single series is its own common interval, and its measured years are returned. If there are two or more series and no year in which any two of them overlap, or no series has any values, there is no common interval and that is an error.

Value

A data.frame with colnames(x) and rownames(x), with no missing values. It has the class of rwl: ring-width indices (class "rwi", see as.rwi) come back as indices. If the years kept are not consecutive, which type = "years" can give when a series has an interior gap, the result is a plain data.frame, with a warning (see [.rwl).

Author(s)

Filipe Campelo, Andy Bunn and Mikko Korpela

See Also

seg.plot

Examples

library(utils)
data(co021)
co021.s <- common.interval(co021, type="series", make.plot=TRUE)
co021.y <- common.interval(co021, type="years", make.plot=TRUE)
co021.b <- common.interval(co021, type="both", make.plot=TRUE)

dim(co021)
dim.s <- dim(co021.s)
dim.s       # the highest number of series
prod(dim.s) #   (33 series x 288 years = 9504)
dim.y <- dim(co021.y)
dim.y       # the highest number of years
prod(dim.y) #   (27 series x 458 years = 12366)
dim.b <- dim(co021.b)
dim.b       # compromise solution
prod(dim.b) #   (28 series x 435 years = 12180)

Description

Computes the correlation between each tree-ring series in a rwl object.

Usage

corr.rwl.seg(rwl, seg.length = 50, bin.floor = 100, n = NULL,
             nyrs = NULL, prewhiten = TRUE, ar.order.max = NULL,
             pcrit = 0.05, biweight = TRUE,
             method = c("spearman", "pearson","kendall"),
             make.plot = TRUE, label.cex = 1, floor.plus1 = FALSE,
             master = NULL, lag.max = 0,
             master.yrs = as.numeric(if (is.null(dim(master))) {
                              names(master)
                          } else {
                              rownames(master)
                          }),
             ...)

Arguments

Details

This function calculates correlation serially between each tree-ring series and a master chronology built from all the other series in the rwl object (leave-one-out principle). Optionally, the user may give a master chronology (a vector) as an argument. In the latter case, the same master chronology is used for all the series in the rwl object. The user can also choose to give a master data.frame (series as columns, years as rows), from which a single master chronology is built.

Correlations are done for each segment of the series where segments are lagged by half the segment length (e.g., 100-year segments would be overlapped by 50-years). The first segment is placed according to bin.floor. The minimum bin year is calculated as ceiling(min.yr/bin.floor)*bin.floor where min.yr is the first year in either the rwl object or the user-specified master chronology, whichever is smaller. For example if the first year is 626 and bin.floor is 100 then the first bin would start in 700. If bin.floor is 10 then the first bin would start in 630.

Correlations are calculated for the first segment, then the second segment and so on. Correlations are only calculated for segments with complete overlap with the master chronology. For now, correlations are Spearman’s rho calculated via cor.test using method = "spearman".

Each series (including those in the rwl object) is optionally detrended as the residuals from a hanning filter with weight n. The filter is not applied if n is NULL. Detrending can also be done via prewhitening where the residuals of an ar model are added to each series mean. This is the default. The master chronology is computed as the mean of the rwl object using tbrm if biweight is TRUE and rowMeans if not. Note that detrending can change the length of the series. E.g., a hanning filter will shorten the series on either end by floor(n/2). The prewhitening default will change the series length based on the ar model fit. The effects of detrending can be seen with series.rwl.plot.

As an alternative to the hanning filter, nyrs divides each series by a smoothing spline, which removes no years from the ends of the series. COFECHA uses a 32-year spline. The spline alone does not usually save years at the start of a series, because the ar model is then fitted to a high-pass filtered series, and AIC tends to choose a high order for it. Setting ar.order.max to a small value (e.g., 3) limits the years lost to prewhitening, and after a spline a low-order model loses little of the fit.

With neither filter, each series is divided by its mean before the master is built, so every series counts equally. Data that can be negative, such as indices from detrend with difference = TRUE, log widths or isotope values, have the mean subtracted instead, since dividing by a negative mean would flip the series. The n and nyrs filters divide by a smooth curve, so they refuse such data. The same applies to the other crossdating functions.

The function is typically invoked to produce a plot where each segment for each series is colored by its correlation to the master chronology. Green segments are those that do not overlap completely with the width of the bin. Blue segments are those that correlate above the user-specified critical value. Red segments are those that correlate below the user-specified critical value and might indicate a dating problem.

A segment that correlates well where it is dated can still correlate better somewhere else. With lag.max greater than 0, each segment is also correlated against the master at every shift from -lag.max to lag.max years. Each year of the segment is paired with the master value k years away, so the correlation at lag k is that of series[t] with master[t + k]. The sign follows ccf.series.rwl with series.x = FALSE: a negative lag means missing rings in the series, a positive lag means false (extra) rings. The dated position wins ties.

This gives the two kinds of flag COFECHA reports. An ‘A’ segment correlates below the critical value, but its dated position is still the best one tested: it is weak, not misdated. A ‘B’ segment correlates better at some other position, whether or not it is significant as dated. The letters are not stored; they follow from the returned matrices as best.lag != 0 for B and best.lag == 0 & p.val >= pcrit for A (see Examples). When lag.max is greater than 0, the plot draws B segments in purple and leaves red for the rest of the segments below the critical value.

A B flag is a hypothesis about the dating, not a verdict. A one-year shift that gains 0.05 in correlation is noise as often as it is a missing ring, which is why the size of the gain is returned (best.rho - spearman.rho) alongside the lag, so that it can be thresholded. Check any B segment against the wood, for example with ccf.series.rwl and xskel.ccf.plot.

Read the lags along a series rather than one bin at a time. Dating runs from the bark inward, so a missing ring moves every ring before it one year late: every complete bin before the missing ring has best.lag of -1 and every bin after it has 0. A false ring does the same with +1. The error therefore lies where the lag changes, not in the bin with the largest gain. The bin that straddles it holds some shifted and some unshifted years, and may show either lag with a small gain. A run of -1 with correctly dated bins on both sides is a different fault: the ring count is right, but a ring is missing at the later end of the run and there is an extra ring at the earlier end. A locally absent ring entered as a zero in too early a year looks exactly like this. A test on the whole series, such as RWL_DATING_LAG in rwl.check, usually cannot see it, because most of the series is dated correctly. The Examples show both patterns.

The shifted correlations follow the same complete-overlap rule as the dated one: a lag is tested only if the master has a value for every year of the shifted window. Near the start and end of the record (and of the master), the lags that run off the record cannot be tested, and a segment that lies within lag.max years of an end is only searched in one direction. A segment that is not complete at its dated position is not tested at all. Floating series and series with dating errors near their ends sit exactly where the test is blind, so the absence of a flag there says little. Worse, a segment misdated in the direction that runs off the record can still be flagged B, at whichever of the remaining lags happens to correlate best, so the lag reported for a B near either end can be noise even when the segment really is misdated.

Value

A list containing matrices spearman.rho, p.val, overall, bins, rwi, vector avg.seg.rho, numeric seg.lag, seg.length, pcrit, label.cex, matrices best.lag, best.rho and numeric lag.max. An additional character flags is also returned if any segments fall below the critical value. Matrix spearman.rho contains the correlations for each series by bin. Matrix p.val contains the p-values on the correlation for each series by bin. Matrix overall contains the average correlation and p-value for each series. Matrix bins contains the years encapsulated by each bin. The vector avg.seg.rho contains the average correlation for each bin. Matrix rwi contains the detrended rwl data, the numerics seg.lag, seg.length, pcrit, label.cex are from the oroginal call and used to pass into plot.crs.

Matrices best.lag and best.rho have the same shape and names as spearman.rho. best.lag gives, for each series and bin, the lag (in years) at which the segment correlates best with the master, and best.rho gives that correlation. Both are NA where spearman.rho is. With lag.max = 0, best.lag is 0 and best.rho equals spearman.rho. lag.max is from the original call. Flags in flags depend on the p-value only and are the same whatever lag.max is.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

corr.series.seg, skel.plot, series.rwl.plot, ccf.series.rwl, plot.crs

Examples

library(utils)
data(co021)
crs <- corr.rwl.seg(co021, seg.length = 100, label.cex = 1.25)
names(crs)
## Average correlation and p-value for the first few series
head(crs$overall)
## Average correlation for each bin
crs$avg.seg.rho

## Plant a missing ring and follow it through the bins: delete the
## 1500 ring of series 641143, the same fault as in the examples for
## ccf.series.rwl(), corr.series.seg() and xskel.ccf.plot().  Dated
## from the bark, every ring before 1500 now sits one year late.
dat <- co021
x <- dat$"641143"
names(x) <- rownames(dat)
dat$"641143" <- delete.ring(x, year = 1500)
crs <- corr.rwl.seg(dat, lag.max = 5)
## -1 in every complete bin before 1500 (drawn in purple), 0 from 1500
## on, and a large gain in each shifted bin
ok <- !is.na(crs$best.lag["641143", ])
rbind(lag = crs$best.lag["641143", ok],
      rho = round(crs$spearman.rho["641143", ok], 2),
      best.rho = round(crs$best.rho["641143", ok], 2))

## The same missing ring with a false ring inserted at 1400.  The ring
## count is right again, and only the rings between the two errors are
## one year late.
dat$"641143" <- insert.ring(delete.ring(x, year = 1500), year = 1400)
crs <- corr.rwl.seg(dat, lag.max = 5, make.plot = FALSE)
## -1 between the two errors only, 0 on both sides
crs$best.lag["641143", ok]

## COFECHA-style A and B flags for every series
flag <- ifelse(crs$best.lag != 0, "B",
               ifelse(crs$p.val >= crs$pcrit, "A", ""))
flag[is.na(flag)] <- ""
## B segments with their lag and the gain in correlation
idx <- which(flag == "B", arr.ind = TRUE)
data.frame(series = rownames(flag)[idx[, 1]],
           bin = colnames(flag)[idx[, 2]],
           lag = crs$best.lag[idx],
           rho = round(crs$spearman.rho[idx], 3),
           best.rho = round(crs$best.rho[idx], 3))

Description

Compute correlation between a tree-ring series and a master chronology by segment.

Usage

corr.series.seg(rwl, series, series.yrs = as.numeric(names(series)),
                seg.length = 50, bin.floor = 100, n = NULL,
                nyrs = NULL, prewhiten = TRUE, ar.order.max = NULL,
                biweight = TRUE,
                method = c("spearman", "pearson","kendall"),
                pcrit = 0.05,
                make.plot = TRUE, floor.plus1 = FALSE, ...)

Arguments

Details

This function calculates the correlation between a tree-ring series and a master chronology built from a rwl object. Correlations are done by segment (see below) and with a moving correlation with length equal to the seg.length. The function is typically invoked to produce a plot.

Value

A list containing matrices bins, moving.rho, and vectors spearman.rho, p.val, and overall.

Matrix bins contains the years encapsulated by each bin (segments). Matrix moving.rho contains the moving correlation and p-value for a moving average equal to seg.length. Vector spearman.rho contains the correlations by bin and p.val contains the p-values. Vector overall contains the average correlation and p-value.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

corr.series.seg, skel.plot, series.rwl.plot, ccf.series.rwl

Examples

library(utils)
data(co021)
dat <- co021
## Create a missing ring: delete the 1500 ring of series 641143 and
## drop the original from the master.  Dated from the bark, every ring
## before 1500 now sits one year late.
flagged <- dat$"641143"
names(flagged) <- rownames(dat)
flagged <- delete.ring(flagged, year = 1500)
dat$"641143" <- NULL
seg.100 <- corr.series.seg(rwl = dat, series = flagged,
                           seg.length = 100, biweight = FALSE)
## Not run: 
flagged2 <- co021$"641143"
names(flagged2) <- rownames(dat)
seg.100.1 <- corr.series.seg(rwl=dat, seg.length=100, biweight=FALSE,
                             series = flagged2)
## Select series by name or column position
seg.100.2 <- corr.series.seg(rwl=co021, seg.length=100, biweight=FALSE,
                             series = "641143")
seg.100.3 <- corr.series.seg(rwl=co021, seg.length=100, biweight=FALSE,
                             series = which(colnames(co021) == "641143"))
identical(seg.100.1, seg.100.2) # TRUE
identical(seg.100.2, seg.100.3) # TRUE

## End(Not run)

Description

This is a wrapper for detrend.series to detrend many ring-width series at once.

Usage

detrend(rwl, y.name = names(rwl), make.plot = FALSE,
        method = "Spline",
        nyrs = NULL, f = 0.5, pos.slope = FALSE,
        constrain.nls = c("never", "when.fail", "always"),
        verbose = FALSE, return.info = FALSE,
        wt, span = "cv", bass = 0, difference = FALSE)

Arguments

Details

See detrend.series for details on detrending methods. Setting make.plot = TRUE will cause plots of each series to be produced. These could be saved using Devices if desired.

A series with no values at all is dropped, with a message naming it, and the rest are detrended as if it had never been there; its name is dropped from y.name too. x[rows, ] and head can leave such a series (see [.rwl). If no series has any values, that is an error. A series with NA inside its measured span stops the call, with a message naming the series; see fill.internal.NA.

Value

If one detrending method is used, an object of class c("rwi", "data.frame") containing the dimensionless detrended ring widths with the column names, row names and dimensions of rwl, less any series with no values. It records how it was made in attr(x, "dplR.detrend") and carries the provenance record of rwl, if it has one; see as.rwi. If more methods are used, a list with ncol(rwl) elements each containing a data.frame with the detrended ring widths in each column.

If return.info is TRUE, the return value is a list with four parts:

Note

This function uses the foreach looping construct with the %dopar% operator. For parallel computing and a potential speedup, a parallel backend must be registered before running the function. If verbose is TRUE, parallel computation is disabled.

Author(s)

Andy Bunn. Improved by Mikko Korpela.

See Also

detrend.series, as.rwi

Examples

library(utils)
data(ca533)
## Detrend using modified exponential decay. Returns an rwi object
ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")
## Detrend using splines and compute
## residuals via subtraction
ca533.rwi <- detrend(rwl = ca533, method = "Spline",
                     difference = TRUE)

## Detrend using modified Hugershoff curve and return info on the model
## fits. Returns a list with: series, curves, modelinfo and data.info
data(co021)
co021.rwi <- detrend(rwl = co021, method = "ModHugershoff",
                     return.info=TRUE)

## Not run: 
library(grDevices)
## Detrend using all methods. Returns a list
ca533.rwi <- detrend(rwl = ca533,
                     method = c("Spline", "ModNegExp", "Mean", "Ar",
                                "Friedman", "ModHugershoff", "AgeDepSpline"))
## Save a pdf of all series
fname <- tempfile(fileext=".pdf")
print(fname) # tempfile used for output
pdf(fname)
ca533.rwi <- detrend(rwl = ca533, method = c("Spline", "ModNegExp"),
                     make.plot = TRUE)
dev.off()

unlink(fname) # remove the file

## End(Not run)

Description

Detrend a tree-ring series by one of two methods, a smoothing spline or a statistical model. The series and fits are plotted by default.

Usage

detrend.series(y, y.name = "", make.plot = TRUE,
               method = "Spline",
               nyrs = NULL, f = 0.5, pos.slope = FALSE,
               constrain.nls = c("never", "when.fail", "always"),
               verbose = FALSE, return.info = FALSE,
               wt, span = "cv", bass = 0, difference = FALSE)

Arguments

Details

This detrends and standardizes a tree-ring series. The detrending is the estimation and removal of the tree’s natural biological growth trend. The default standardization is done by dividing each series by the growth trend to produce units in the dimensionless ring-width index (RWI). If difference is TRUE, the index is calculated by subtraction. Values of zero (typically missing rings) in y are set to 0.001 when dividing.

A curve that is divided by must be positive, so when difference is FALSE a fit that is not all positive is replaced by a straight line or the mean (or, for "Ar", its negative values are set to zero) with a warning. Subtraction has no such need, so with difference = TRUE the fits are used as they are, zeros are left alone, and y may be negative, as log widths or isotope values are.

There are currently seven methods available for detrending although more are certainly possible. The methods implemented are an age-dependent spline via ads (method = "AgeDepSpline"), the residuals of an AR model (method = "Ar"), Friedman's Super Smoother (method = "Friedman"), a simple horizontal line (method = "Mean"), or a modified Hugershoff curve (method = "ModHugershoff"), a modified negative exponential curve (method = "ModNegExp"), and a smoothing spline via caps (method = "Spline").

The "AgeDepSpline" approach uses an age-dependent spline with an initial stiffness of 50 (nyrs=50). See ads. If some of the fitted values are not positive then method "Mean" is used. However, this is extremely unlikely.

The "Ar" approach is also known as prewhitening where the detrended series is the residuals of an ar model divided by the mean of those residuals to yield a series with white noise and a mean of one. This method removes all but the high frequency variation in the series and should only be used as such.

The "Friedman" approach uses Friedman’s ‘super smoother’ as implemented in supsmu. The parameters wt, span and bass can be adjusted, but periodic is always set to FALSE. If some of the fitted values are not positive then method "Mean" is used.

The "Mean" approach fits a horizontal line using the mean of the series. This method is the fallback solution in cases where the "Spline" or the linear fit (also a fallback solution itself) contains zeros or negative values, which would lead to invalid ring-width indices.

The "ModHugershoff" approach attempts to fit a Hugershoff model of biological growth of the form $f(t) = a t^b e^{-g t} + d$, where the argument of the function is time, using nls. See Fritts (2001) for details about the parameters. Option constrain.nls gives a possibility to constrain the parameters of the modified negative exponential function. If the constraints are enabled, the nonlinear optimization algorithm is instructed to keep the parameters in the following ranges: $a \ge 0$, $b \ge 0$ and $d \ge 0$. The default is to not constrain the parameters (constrain.nls = "never") for nls but warn the user when the parameters go out of range.

If a suitable nonlinear model cannot be fit (function is non-decreasing or some values are not positive) then a linear model is fit. That linear model can have a positive slope unless pos.slope is FALSE in which case method "Mean" is used.

The "ModNegExp" approach attempts to fit a classic nonlinear model of biological growth of the form $f(t) = a e^{b t} + k$, where the argument of the function is time, using nls. See Fritts (2001) for details about the parameters. Option constrain.nls gives a possibility to constrain the parameters of the modified negative exponential function. If the constraints are enabled, the nonlinear optimization algorithm is instructed to keep the parameters in the following ranges: $a \ge 0$, $b \le 0$ and $k \ge 0$. The default is to not constrain the parameters (constrain.nls = "never") for nls but warn the user when the parameters go out of range.

If a suitable nonlinear model cannot be fit (function is non-decreasing or some values are not positive) then a linear model is fit. That linear model can have a positive slope unless pos.slope is FALSE in which case method "Mean" is used.

The "Spline" approach uses a spline where the frequency response is 0.50 at a wavelength of 0.67 * “series length in years”, unless specified differently using nyrs and f in the function caps. If some of the fitted values are not positive then method "Mean" is used.

These methods are chosen because they are commonly used in dendrochronology. There is a rich literature on detrending and many researchers are particularly skeptical of the use of the classic nonlinear model of biological growth ($f(t) = a e^{b t} + k$) for detrending. It is, of course, up to the user to determine the best detrending method for their data.

Note that the user receives a warning if a series has negative values in the fitted curve. This happens fairly commonly with with the ‘Ar’ method on high-order data. It happens less often with method ‘Spline’ but isn't unheard of (see ‘Examples’). If this happens, users should look carefully at their data before continuing. Automating detrending and not evaluating each series individually is folly. Remember, frustration over detrending is the number one cause of dendros going to live as hermits in the tallgrass prairie, where there are no trees to worry about.

See the references below for further details on detrending. It's a dark art.

Value

If several methods are used, returns a data.frame containing the detrended series (y) according to the methods used. The columns are named and ordered to match method. If only one method is selected, returns a vector.

If return.info is TRUE, the return value is a list with four parts:

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela. A bug fix related to negative output values is based on work by Alice Cecile.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Fritts, H. C. (2001) Tree Rings and Climate. Blackburn. ISBN-13: 978-1-930665-39-2.

See Also

detrend

Examples

library(stats)
library(utils)
## Use series CAM011 from the Campito data set
data(ca533)
series <- ca533[, "CAM011"]
names(series) <- rownames(ca533)
# all seven methods
series.rwi <- detrend.series(y = series, y.name = "CAM011", verbose=TRUE,
                             method = c("Spline", "ModNegExp", "Mean", "Ar",
                                        "Friedman", "ModHugershoff",
                                        "AgeDepSpline"))
# see plot with three methods
series.rwi <- detrend.series(y = series, y.name = "CAM011",
                             method=c("Spline", "ModNegExp","Friedman"),
                             difference=TRUE)
# see plot with two methods
# interesting to note difference from ~200 to 250 years 
# in terms of what happens to low frequency growth
series.rwi <- detrend.series(y = series, y.name = "CAM011",
                             method=c("Spline", "ModNegExp"))
# see plot with just one method and change the spline
# stiffness to 50 years (which is not neccesarily a good choice!)
series.rwi <- detrend.series(y = series, y.name = "CAM011",
                             method="Spline",nyrs=50)
                             
# note that method "Ar" doesn't get plotted in first panel
# since this approach doesn't approximate a growth curve.
series.rwi <- detrend.series(y = series, y.name = "CAM011",
                             method="Ar")
                             
# note the difference between ModNegExp and ModHugershoff at the 
# start of the series. Also notice how curves, etc. are returned
# via return.info
data(co021)
series <- co021[, 4]
names(series) <- rownames(co021)
series.rwi <- detrend.series(y = series, y.name = names(co021)[4],
                             method=c("ModNegExp", "ModHugershoff"),
                             verbose = TRUE, return.info = TRUE, 
                             make.plot = TRUE)
                             
# A dirty dog.
# In the case of method=="Spline" the function carries-on
# and applies method=="Mean" as an alternative. 
data(nm046)
series <- nm046[,8]
names(series) <- rownames(nm046)
series.rwi <- detrend.series(y = series, y.name = names(nm046)[8],
                             method="Spline",
                             make.plot = FALSE)

Description

These are functions no longer included in dplR

Details

Description

Subset a rwl object by series (columns) and years (rows), keeping the class and the provenance record of the data.

Usage

## S3 method for class 'rwl'
x[i, j, drop]

## S3 method for class 'rwl'
subset(x, subset, select, drop = FALSE, ...)

Arguments

Details

Subsetting an rwl object works as it does for a data.frame, with i indexing years and j indexing series. Four things are different.

Dropping series returns the years the remaining series cover. A collection is ragged, and the years at the ends of an rwl object usually belong to a few long series, so dropping some of the series ordinarily leaves years in which nothing that remains was measured. The result runs from the first year to the last year in which any of the selected series has a measurement. Empty years between those two are kept: an rwl object holds one row per year, and an interior year cannot be dropped without breaking the sequence.

This applies to calls that drop series without indexing years – x[, j], x[j], and subset given only select. Two kinds of call leave the years alone. A call that indexes years returns exactly the years it asks for, measured or not: x[i, ], x[i, j], head, tail, common.interval, and subset with its subset argument. A year window is how an rwl object is aligned with something else, such as a climate series or a second collection, and returning fewer years than were asked for would put that alignment out. And a call that keeps every series, such as the reordering x[, order(names(x))] or x[], has dropped nothing to trim for, whatever empty years the object already held.

Two results are not trimmed. A subset holding no measurement in any year is returned whole rather than emptied, and rwl.check reports it. A single series taken with drop = TRUE is a numeric vector, which carries no years at all, so it is returned at full length and stays aligned with the object it came from; use drop = FALSE for a one-series rwl object trimmed to its own years.

subset drops empty series; [ keeps them. Taking years can leave a series with no values at all: in ca533[time(ca533) %in% 1800:1899, ], the series CAM132, CAM152, CAM161 and CAM201 end before 1800 and are left as columns of nothing but NA. An empty column is not a series. chron, rwi.stats and the plots skip it, and summary.rwi, interseries.cor, corr.rwl.seg and detrend either drop it or fail on it.

subset and window.rwl therefore drop any series left with no values, and a message names the series dropped. If no series is left with values, that is an error. Use one of these to take a span of years for analysis.

[ keeps every series it is asked for, empty or not, on purpose. x[i, j] always returns the columns that j names, in that order, so the result stays lined up with anything indexed alongside it: a vector of series names, an ids table from read.ids, or a second object with the same series. If [ dropped columns, those would fall out of step without a word, and code that indexes with [ – dplR's own included – could no longer count on names(x[, j]) being the series it asked for. So row subsets made with [, and head and tail, can return empty series. Drop them with window or subset, or find them with colSums(!is.na(x)) == 0.

Years must stay consecutive. An rwl object holds one row per year in increasing order, and dplR reads the year of a measurement from its row name, so time.rwl, plot.rwl, detrend, chron, rwl.stats and the crossdating functions all take the next row to be the next year. A subset such as x[c(1, 5, 9), ] leaves that untrue. Where row subsetting leaves years that are not consecutive and increasing, because rows were taken from the middle, repeated, or reordered, the result is returned as a plain data.frame, with a warning, and without the provenance record. The measurements are unchanged; what is withdrawn is the claim that they are a set of dated series.

The provenance record is carried and cut to fit. The "dplR.provenance" attribute attached by read.tucson follows the subset. What describes the file – its name, the reader, the header lines – is kept. What describes series – the precision of each, the renames, the interior gaps, the problems found while parsing – is kept for the series that remain, and gaps are dropped where they fall outside the years that remain. mixed.precision is recomputed, since a file measured at two precisions can be subset down to series measured at one. A subset element is added, holding the series in the result, its first and last year, and all.series, which is TRUE when every series the file held is still present. A record that cannot be matched to a series in the result is dropped rather than kept, so an object whose series have been renamed since it was read carries a record that says nothing about them in preference to one that says something wrong.

common.interval cuts an object down to the years its series share, which is a different question from the years they cover between them.

Value

An object of class c("rwl", "data.frame"), or a numeric vector when one column is selected and drop is TRUE, or a data.frame when row subsetting has left the years no longer consecutive. When series are dropped and years are not indexed, the result runs only over the years the remaining series cover. [ returns every series asked for, including any with no values; subset drops those, with a message.

Author(s)

Andy Bunn

See Also

read.rwl, read.tucson, rwl.check, common.interval, time.rwl, subset.data.frame

Examples

library(utils)
data(ca533)

## Series are columns and years are rows. Selecting series drops the years
## that none of them covers: ca533 begins in 626, but these five series do
## not.
ca533.sub <- ca533[, 1:5]
class(ca533.sub)
range(time(ca533))
range(time(ca533.sub))

## Naming years returns those years, whether or not they hold data.
dim(ca533[time(ca533) %in% 1000:1099, 1:5])

## subset() works the same way.
range(time(subset(ca533, select = 1:5)))

## Years are subset by name or by time().
ca533.1800s <- ca533[time(ca533) %in% 1800:1899, ]
range(time(ca533.1800s))

## Four series end before 1800. `[` keeps them as empty columns, so the
## result still has the 34 series asked for; subset() and window() drop
## them and say which.
ncol(ca533.1800s)
names(ca533.1800s)[colSums(!is.na(ca533.1800s)) == 0]
ncol(subset(ca533, time(ca533) %in% 1800:1899))
ncol(window(ca533, 1800, 1899))

## A single series comes back as a vector unless drop is FALSE.
class(ca533[, 1])
class(ca533[, 1, drop = FALSE])

## The provenance record of the data follows the subset.
data(wa082)
gap.series <- attr(wa082, "dplR.provenance")$gaps$series
wa082.sub <- wa082[, gap.series, drop = FALSE]
attr(wa082.sub, "dplR.provenance")$gaps
attr(wa082.sub, "dplR.provenance")$subset$all.series

## Rows taken out of the middle break the one-row-per-year promise that
## dplR relies on. This warns, and gives back a data.frame.
not.rwl <- ca533[c(1, 5, 9), ]
class(not.rwl)

Description

This function fills internal NA values (i.e., those with numeric data above and below a small data gap) in each column of a data.frame such as a data set of ring widths as produced by read.rwl.

Usage

fill.internal.NA(x, fill = c("Mean", "Spline", "Linear", "Chron"),
                 series = names(x), nyrs = 20, min.series = 3)

Arguments

Details

There are occasionally data gaps within a tree-ring series. Some of the functions in dplR will fail when an internal NA is encountered (e.g. caps). This function fills internal NA values with either a given numeric value (e.g., 0) or through crude imputation. The latter can be calculated as the mean of the series (fill="Mean") or calculated by fitting a cubic spline (fill="Spline") using the spline function or calculated by linear approximation (fill="Linear") using the function approx. Those three look only at the series being filled. fill="Chron" takes the value from the other series in x and is described in its own section below.

Editorial: Having internal NA in a tree-ring series is often bad practice and filling those values should be done with caution. For instance, some users code missing rings as NA instead of 0. And missing values (i.e., NA) are sometimes present in maximum latewood density data when the rings are small. A common, but not recommended, practice is to leave stretches of NA values in places where it has been impossible to accurately measure rings (perhaps because of a break in the core). It is often better to treat that core as two separate series (e.g., "01A" and "01B" rather than have internal NA values. As with all processing, the analyst should make a decision based on their experience with the wood and not rely on software to make a choice for them!

Value

x, with the internal NAs of the series named in series filled as above. It keeps the class and attributes of x: an rwl object comes back as an rwl object, with the record read.tucson attached to it. Before dplR 1.8.1 the result was a plain data.frame.

Filling from the other series (fill = "Chron")

A mean, a line or a spline through one series can say what size of ring to expect in a gap, but nothing about whether the missing years were wide or narrow ones. The other trees at the site can. fill = "Chron" estimates each missing ring from them.

The method is in the spirit of the one in the program ARSTAN (Cook and Krusic), whose subroutine fillin fills every gap in a data set before detrending. It follows the same outline:

1. A spline of nyrs years (caps) is fitted to every series in x and each series is divided by its spline, giving indices. The default of 20 years is ARSTAN's.

2. For the series being filled, the indices of all the other series are averaged year by year with a robust mean (tbrm).

3. In each gap year, the departure of that mean from its own average is expressed in standard deviations, and then given the mean and standard deviation of the indices of the series being filled.

4. The result is multiplied by the spline of the series being filled, to return to ring width. An estimate below zero is set to zero.

It is not a port of fillin and will not reproduce ARSTAN's numbers. The differences are deliberate:

• ARSTAN fits a weighted spline that ignores the gap years. Here the gap is first bridged with a straight line and the spline is fitted to the bridged series.

• ARSTAN's mean includes the series being filled. Here it is left out.

• ARSTAN detrends the mean a second time and stabilises its variance. Neither is done here.

• ARSTAN fills from a single series if that is all there is and leaves the value missing if there is none. Here a year with fewer than min.series other series stops the call.

• ARSTAN treats any negative value as a gap, and zero as well when the data are transformed. Here only NA is a gap.

Assumptions:

• The data are ring widths or another quantity for which a ratio to a growth curve makes sense: no negative values, and a curve that stays above zero.

• The series share a common year-to-year signal. The method assumes the series being filled would have followed the others in the years it is missing. Where the series do not crossdate, the filled values are noise with a convincing shape.

• The series are correctly dated. A dating error in the series being filled, or in several of the others, puts the wrong year's signal in the gap.

• The gap is short relative to nyrs. The growth curve across a gap comes from a straight line between the rings either side of it, so over a long gap the level of the filled values is a guess, whatever their year-to-year pattern.

Limitations:

• The filled values are estimates, not measurements, and nothing in the returned object marks which they are. Keep the unfilled data.

• The fill is circular. The values come from the other series and then sit beside them, so statistics that measure agreement among series (rwi.stats, interseries.cor, corr.rwl.seg) are biased upwards in the filled years, and sample depth counts series that were not measured there. Do not crossdate with filled values.

• The filled values are given the full year-to-year variance of the series, as if it tracked the others perfectly. That keeps the variance of the series from dropping in the gap, at the cost of filled values that vary more than a regression estimate would.

• An estimate that falls below zero is returned as zero, with a warning naming the series and years. A zero ring width reads as a locally absent ring everywhere else in dplR.

• A year can be filled only if at least min.series other series were measured in it. If any gap year falls short, the call stops, names the series and years, and fills nothing. With min.series of 1 or 2 the “mean” is one tree or the average of two.

• The series being filled and the mean of the others need at least three years in common to set the scale, and a scale set from a few years is poor. Short series are filled badly.

• A series with fewer than three measurements, or whose spline reaches zero, has no usable growth curve. If it is one of the series to be filled the call stops. If it is one of the others it is left out of the mean without comment.

• Gaps are not filled from each other: a year that is missing in several series is estimated for each from the series that have it.

• Filled values do not fall on the measurement grid of the data (0.01 or 0.001 mm). They are not rounded.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

References

Cook, E. R. and Krusic, P. J. Program ARSTAN: a tree-ring standardization program based on detrending and autoregressive time series modeling, with interactive graphics. Tree-Ring Laboratory, Lamont-Doherty Earth Observatory of Columbia University. The "Chron" method follows the outline of subroutine fillin in version 49v1 of the Fortran source.

See Also

spline, approx, caps, tbrm

Examples

library(graphics)
library(stats)
foo <- data.frame(x1=c(rnorm(5), NA, NA, rnorm(3)),
                  x2=c(rnorm(10)),
                  x3=c(NA, NA, rnorm(3), NA, rnorm(4)),
                  x4=c(NA, NA, rnorm(3), NA, rnorm(3), NA),
                  x5=c(NA, NA, rnorm(8)),
                  x6=c(NA, rnorm(9)),
                  x7=c(NA, rnorm(5), NA, rnorm(3)),
                  x8=c(rnorm(8), NA, NA),
                  x9=c(rnorm(5), NA, rnorm(3), NA))
row.names(foo) <- 1901:1910
class(foo) <- c("rwl","data.frame")
fill.internal.NA(foo, fill=0)

bar <- fill.internal.NA(foo, fill="Spline")
baz <- fill.internal.NA(foo, fill="Linear")

## note differences in method "Spline" vs. "Linear"
yrs <- time(foo)
plot(yrs, foo$x7, type="b", lwd=3)
lines(yrs, bar$x7, col="red", lwd=2)
lines(yrs, baz$x7, col="green", lwd=1)

## fill = "Chron" uses the other series. Remove ten rings from one
## series of co021 and compare the estimates with what was measured.
data(co021)
gap <- 1701:1710
gap2 <- as.character(gap)
dat <- co021
dat[gap2, "641114"] <- NA
chr <- fill.internal.NA(dat, fill="Chron", series="641114")
lin <- fill.internal.NA(dat, fill="Linear", series="641114")
plot(gap, co021[gap2, "641114"], type="b", lwd=3,
     xlab="Year", ylab="Ring width (mm)", ylim = c(0.1,0.75))
lines(gap, chr[gap2, "641114"], col="darkred", lwd=2)
lines(gap, lin[gap2, "641114"], col="darkgreen", lwd=2)

Description

This function calculates the Gini coefficient on raw or detrended ring-width series.

Usage

gini.coef(x)

Arguments

Details

This calculates the Gini coefficient of inequality which is used as an all-lag measure of diversity in tree-ring records – typically detrended series. Lower values indicate lower diversity. The use of the Gini coefficient in dendrochronology is described by Biondi and Qeadan (2008). See Handcock and Morris (1999) for more information.

Value

the Gini coefficient.

Note

The equivalence between the Lorenz-curve formula used in the C code and the relative-mean-difference definition of the Gini coefficient is derived in the vignette: vignette("math-dplR", package = "dplR").

Author(s)

Mikko Korpela, based on original by Andy Bunn

References

Biondi, F. and Qeadan, F. (2008) Inequality in Paleorecords. Ecology, 89(4), 1056–1067.

Handcock, M. S. and Morris, M. (1999) Relative Distribution Methods in the Social Sciences. Springer. ISBN: 0-387-98778-9.

See Also

rwl.stats

Examples

library(utils)
data(ca533)
ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")
ca533.crn <- chron(ca533.rwi)
gini.coef(ca533.crn)

Description

This function calculates the Gleichläufigkeit and related measures for a given set of tree-ring records.

Usage

glk(x, overlap = 50, prob = TRUE)
  glk.legacy(x)

Arguments

Details

Gleichläufigkeit is a classical agreement test based on sign tests (Eckstein and Bauch, 1969). This function implements Gleichläufigkeit as the pairwise comparison of all records in data set. This vectorized implementation is faster than the previous version and follows the original definition (Huber 1942), instead of the incorrect interpretation that has been used in the past (Schweingruber 1988, see Buras/Wilmking 2015 for the correction).

The probability of exceedence (p) for the Gleichläufigkeit expresses the chance that the Gleichläufigkeit is incorrect. The observed value of the Gleichläufigkeit is converted to a z-score and based on the standard normal curve the probability of exceedence is calculated. The result is a matrix of all p-values (Jansma 1995, 60-61, see also Visser 2020).

Note that prior to dplR version 1.7.2, glk did not have the overlap or prob and returned a matrix with just the Gleichläufigkeit for all possible combinations of records. That function can still be accessed via glk.legacy.

Value

The funtions returns a named list of two or three matrices (p_mat is optional if prob = TRUE):

1. glk_mat: matrix with Gleichläufigkeit

2. overlap: matrix with number of overlapping growth changes.This is the number of overlapping years minus one.

3. p_mat: matrix of all probabilities of exceedence for all observed Gleichläufigkeit values.

The matrices can be extracted from the list by selecting the name or index number. If two curves have less than 3 years of overlap, Gleichläufigkeit cannot be computed, and NA is returned. To calculate the global glk of the dataset mean(x$glk_mat, na.rm = TRUE). See examples for calculating the the global glk without self-similarities (i.e., the ones on the diagonal).

Author(s)

Christian Zang. Patched and improved by Mikko Korpela. Improved by Allan Buras. Further improved and expanded by Ronald Visser and Andy Bunn

References

Buras, A. and Wilmking, M. (2015) Correcting the calculation of Gleichläufigkeit, Dendrochronologia 34, 29-30. DOI: https://doi.org/10.1016/j.dendro.2015.03.003

Eckstein, D. and Bauch, J. (1969) Beitrag zur Rationalisierung eines dendrochronologischen Verfahrens und zur Analyse seiner Aussagesicherheit. Forstwissenschaftliches Centralblatt, 88(1), 230-250.

Huber, B. (1943) Über die Sicherheit jahrringchronologischer Datierung. Holz als Roh- und Werkstoff 6, 263-268. DOI: https://doi.org/10.1007/BF02603303

Jansma, E., 1995. RemembeRINGs; The development and application of local and regional tree-ring chronologies of oak for the purposes of archaeological and historical research in the Netherlands, Nederlandse Archeologische Rapporten 19, Rijksdienst voor het Oudheidkundig Bodemonderzoek, Amersfoort

Schweingruber, F. H. (1988) Tree rings: basics and applications of dendrochronology, Kluwer Academic Publishers, Dordrecht, Netherlands, 276 p.

Visser, R.M. (2020) On the similarity of tree-ring patterns: Assessing the influence of semi-synchronous growth changes on the Gleichläufigkeit for big tree-ring data sets,Archaeometry, 63, 204-215 DOI: https://doi.org/10.1111/arcm.12600

See Also

sgc sgc is an alternative for glk)

Examples

library(utils)
  data(ca533)
  ca533.glklist <- glk(ca533)
  ca533.glk_mat <- ca533.glklist$glk_mat
  # calculating the mean GLK including self-similarities
  mean(ca533.glk_mat, na.rm = TRUE) 
  # calculating the mean GLK excluding self-similarities
  mean(ca533.glk_mat[upper.tri(ca533.glk_mat)], na.rm = TRUE)

Ponderosa Pine Distance to Pith Corresponding to gp.rwl

Description

This data set gives the distance to pith for each series (in mm) that matches the ring widths for gp.rwl – a data set of ponderosa pine (Pinus ponderosa) from the Gus Pearson Natural Area (GPNA) in northern Arizona, USA. Data are further described by Biondi and Qeadan (2008) and references therein.

Usage

data(gp.d2pith)

Format

A data.frame with 58 rows, one for each series of gp.rwl and in the order of its columns, containing series IDs in column 1 (series) and the distance (in mm) from the innermost ring to the pith of the tree (d2pith). This can be used together with the ring widths to calculate the area of each ring.

Source

DendroLab, University of Nevada Reno, USA. https://dendrolaborg.wordpress.com/

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

See Also

gp.rwl, gp.po, gp.dbh, bai.in

Ponderosa Pine Stem Diameters and Bark Thickness (gp.rwl)

Description

This data set gives the diameter at breast height for each series that matches the series in gp.rwl – a data set of ponderosa pine (Pinus ponderosa) from the Gus Pearson Natural Area (GPNA) in northern Arizona, USA. Data are further described by Biondi and Qeadan (2008) and references therein.

Usage

data(gp.dbh)

Format

A data.frame with 58 rows, one for each series of gp.rwl and in the order of its columns, containing series IDs in column 1 (ID), tree diameter (in cm) at breast height (DBH), and the bark thickness (in cm, BARK). The two cores of a tree have the same diameter. This can be used together with the ring widths to calculate the area of each ring: see the example in bai.out.

Source

DendroLab, University of Nevada Reno, USA. https://dendrolaborg.wordpress.com/

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

See Also

gp.rwl, gp.po, gp.d2pith, bai.out

Ponderosa Pine Pith Offsets Corresponding to gp.rwl

Description

This data set gives the pith offsets that match the ring widths for gp.rwl – a data set of ponderosa pine (Pinus ponderosa) from the Gus Pearson Natural Area (GPNA) in northern Arizona, USA. Data are further described by Biondi and Qeadan (2008) and references therein.

Usage

data(gp.po)

Format

A data.frame with 58 rows, one for each series of gp.rwl and in the order of its columns, and two variables:

series

the series ID, a factor. It matches names(gp.rwl).

pith.offset

an integer: the number of years between the beginning of that series in gp.rwl and the pith of the tree. The values run from 1 to 249.

Details

A pith offset says how old the tree was, at the height of the core, when the first measured ring was laid down. An offset of 1 means the first ring of the series is the first ring of the tree. Added to the position of a ring in its series, it gives the cambial age of that ring: the first ring of a series with an offset of 79 has a cambial age of 79.

This is the form of po that rcs and cms take, and the data can be used as they are with gp.rwl. See the examples.

The two cores of a tree have offsets of their own.

Source

DendroLab, University of Nevada Reno, USA. https://dendrolaborg.wordpress.com/

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

See Also

gp.rwl, gp.d2pith, gp.dbh, rcs, cms

Examples

data(gp.rwl)
data(gp.po)
head(gp.po)
summary(gp.po$pith.offset)

## the offsets are in the order of the series
identical(as.character(gp.po$series), names(gp.rwl))

## cambial age of each ring of the first series
age <- gp.po$pith.offset[1] - 1 + seq_len(sum(!is.na(gp.rwl[[1]])))
range(age)

## as the po argument of cms() and rcs()
gp.rwi <- cms(rwl = gp.rwl, po = gp.po)

Description

This data set includes ring-width measurements for ponderosa pine (Pinus ponderosa) increment cores collected at the Gus Pearson Natural Area (GPNA) in northern Arizona, USA. There are 58 series from 29 trees (2 cores per tree). Data are further described by Biondi and Qeadan (2008) and references therein.

Usage

data(gp.rwl)

Format

An object of class c("rwl", "data.frame") containing 58 ring-width series in columns and 421 years in rows, 1570 to 1990.

Each series is named by a two-digit tree number and a letter for the core, A or B: "01A" and "01B" are the two cores of tree 1. The tree numbers run from 1 to 59 with gaps.

Ring widths are in mm, to 0.01 mm. The series hold 16,408 ring widths in all, of which 194 are zero. Series are between 92 and 418 rings long. All end in 1990 but one, which ends in 1987. No series has a year without a measurement inside its span.

Details

Three other data sets hold one row for each series of gp.rwl, in the order of its columns:

• gp.po, the pith offsets, used by rcs and cms;

• gp.d2pith, the distance from the innermost ring to the pith, used by bai.in;

• gp.dbh, stem diameter and bark thickness, used by bai.out.

With two cores for every tree, these data are useful for anything that separates agreement within trees from agreement between them. read.ids(gp.rwl, stc = c(0, 2, 1)) gives the tree and core of each series for rwi.stats, sss and summary.rwi.

The Gus Pearson Natural Area is a permanent plot of 800 by 400 m, established in 1908 in ponderosa pine forest near Flagstaff, Arizona. Its history, and the growth of its trees through the twentieth century, are described by Biondi (1996, 1999).

Source

DendroLab, University of Nevada Reno, USA. https://dendrolaborg.wordpress.com/

References

Biondi, F. (1996) Decadal-scale dynamics at the Gus Pearson Natural Area: evidence for inverse (a)symmetric competition? Canadian Journal of Forest Research, 26(8), 1397–1406.

Biondi, F. (1999) Comparing tree-ring chronologies and repeated timber inventories as forest monitoring tools. Ecological Applications, 9(1), 216–227.

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

See Also

gp.po, gp.d2pith, gp.dbh, read.ids

Examples

data(gp.rwl)
summary(gp.rwl)[1:6, ]

## which series are cores of the same tree
gp.ids <- read.ids(gp.rwl, stc = c(0, 2, 1))
head(gp.ids)
## two cores for each of 29 trees
table(table(gp.ids$tree))

## the pith offsets, distances to pith and diameters are in the
## order of the series
data(gp.po)
identical(as.character(gp.po$series), names(gp.rwl))

Description

Applies a Hanning filter of length n to x.

Usage

hanning(x, n = 7)

Arguments

Details

This applies a low frequency Hanning (a.k.a. Hann) filter to x with weight set to n.

Value

A filtered vector.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

References

Oppenheim, A. V., Schafer, R. W., and Buck, J. R. (1999) Discrete-Time Signal Processing. Prentice-Hall, 2nd edition. ISBN-13: 978-0-13-754920-7.

See Also

filter

Examples

library(graphics)
library(utils)
data(ca533)
yrs <- time(ca533)
y <- ca533[, 1]
not.na <- !is.na(y)
yrs <- yrs[not.na]
y <- y[not.na]
plot(yrs, y, xlab = "Years", ylab = "Series1 (mm)",
     type = "l", col = "grey")
lines(yrs, hanning(y, n = 9), col = "red", lwd = 2)
lines(yrs, hanning(y, n = 21), col = "blue", lwd = 2)
legend("topright", c("Series", "n=9", "n=21"),
       fill=c("grey", "red", "blue"))

Description

Interactively detrend multiple tree-ring series by one of two methods, a smoothing spline or a statistical model. This is a wrapper for detrend.series.

Usage

i.detrend(rwl, y.name = names(rwl), nyrs = NULL, f = 0.5,
          pos.slope = FALSE)

Arguments

Details

This function allows a user to choose detrending curves based on plots that are produced by detrend.series for which it is essentially a wrapper. The user enters their choice of detrended method via keyboard at a prompt for each ring width series in rwl. See detrend.series for examples and details on the detrending methods.

A series with no values at all is dropped, with a message naming it, before any series is shown, and its name is dropped from y.name too.

Value

An object of class c("rwi", "data.frame") (see as.rwi) containing each detrended series according to the method used as columns and rownames set to colnames(y). These are typically years. The method chosen for each series is recorded, named by series, in attr(x, "dplR.detrend")$method. Plots are also produced as the user chooses the detrending methods through keyboard input.

Note

This function is kept for old scripts. For interactive detrending use the iDetrend app (https://github.com/OpenDendro/iDetrend), part of openDendro. There the curve for each series can be adjusted and compared with others on screen, series that need attention are pointed out, and the app writes R code that reproduces the indices. Here the choice is made once per series from a static plot, the parameters of the curves cannot be changed series by series, and only the name of the method chosen is recorded. The function says so, once in a session, when it is called.

Author(s)

Andy Bunn

See Also

detrend.series

Description

Interactively detrend a tree-ring series by one of three methods, a smoothing spline, a linear model, or the mean. This is a wrapper for detrend.series.

Usage

i.detrend.series(y, y.name = NULL, nyrs = NULL, f = 0.5,
                 pos.slope = FALSE)

Arguments

Details

This function allows a user to choose a detrending method based on a plot that is produced by detrend.series for which it is essentially a wrapper. The user enters their choice of detrended method via keyboard at a prompt. See detrend.series for examples and details on the detrending methods.

Value

A vector containing the detrended series (y) according to the method used with names set to colnames(y). These are typically years. The name of the method chosen is attached as attr(x, "method"). A plot is also produced and the user chooses a method through keyboard input.

Note

This function is kept for old scripts. For interactive detrending use the iDetrend app (https://github.com/OpenDendro/iDetrend), part of openDendro. There the curve for each series can be adjusted and compared with others on screen, series that need attention are pointed out, and the app writes R code that reproduces the indices. Here the choice is made once per series from a static plot, the parameters of the curves cannot be changed series by series, and only the name of the method chosen is recorded. The function says so, once in a session, when it is called.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

detrend.series

Description

Insert or delete rings from a ring-width series, given either as a vector or as one series of an rwl object.

Usage

insert.ring(rw.vec, ...)
delete.ring(rw.vec, ...)

## Default S3 method:
insert.ring(rw.vec, rw.vec.yrs = as.numeric(names(rw.vec)),
            year, ring.value = mean(rw.vec, na.rm = TRUE),
            fix.last = TRUE, fix.length = TRUE, ...)
## Default S3 method:
delete.ring(rw.vec, rw.vec.yrs = as.numeric(names(rw.vec)),
            year, fix.last = TRUE, fix.length = TRUE, ...)

## S3 method for class 'rwl'
insert.ring(rw.vec, series, year,
            ring.value = mean(rw.vec[[series]], na.rm = TRUE),
            fix.last = TRUE, ...)
## S3 method for class 'rwl'
delete.ring(rw.vec, series, year, fix.last = TRUE, ...)

Arguments