Movatterモバイル変換

Dendrochronology Program Library inR

Description

This package contains functions for performing some standard tree-ringanalyses.

Details

Main Functions

read.rwl

reads rwl files

detrend

detrends raw ring widths

chron

builds chronologies

corr.rwl.seg

crossdating function

Author(s)

Andy Bunnandy.bunn@wwu.edu with major additions from MikkoKorpela and other significant contributions from Franco Biondi, FilipeCampelo, Pierre Mérian, Fares Qeadan and ChristianZang. Functionredfit is an improved translation ofprogram REDFIT which is original work of Manfred Mudelsee and MichaelSchulz. Jacob Cecile contributed a bug fix todetrend.series. Allan Buras came up with the revisedformula ofglk in dplR >= 1.6.1.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990)Methods ofDendrochronology: 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.

Age-Dependent Spline

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

This implements the age-dependent smoothing spline similar to that described by Melvin (2004). In this implementation a cubic smoothing spline (caps) is fit toy with an initial stiffness ofnyrs0. For each succesive measurement, the siffness is incremented by that ring index. This results in a spline isnyrs0 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 fornyrs0 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.

Ifpos.slope is TRUE, the function will attempt to prevent a postive slope at the end of the series. In some cases whenads is used for detrending, a postive slope can be considered considered biologically unlikely. In those cases, the user can contrain the positive slope in the slpine. 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 contrained curve. See examples for details. The wisdom of contraining 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 curvedata(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.5b <- 1g <- 0.1d <- 0.25n <- 300x <- 1:ny <- I(a*x^b*exp(-g*x)+d)# add some noisey <- 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) # badlines(ads(y,10),col="darkblue",lwd=2) # good# now repeat with a positive slope to constrainy <- 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")

Rothenburg Tree Ring Widths

Description

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

Usage

data(anos1)

Format

Adata.frame containing 20 tree-ring series from 10 treesin columns and 98 years in rows. The correct stc mask for use withread.ids isc(5, 2, 1).

References

Zang, C. (2010)Growth reaction of temperate forest tree speciesto summer drought – a multispecies tree-ring networkapproach. Ph.D. thesis, Technische UniversitätMünchen.

as.rwl

Description

Attempts to turn its argument into a rwl object.

Usage

as.rwl(x)

Arguments

Details

This tries to coercex into classc("rwl","data.frame"). Failable.

Value

An object of classc("rwl", "data.frame") with the series incolumns and the years as rows. The seriesIDs are thecolumn 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)## Toyn <- 100## Make a data.frame that is tree-ring likebase.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 methodsfoo <- as.rwl(foo)class(foo)plot(foo, plot.type="spag")summary(foo)

Basal Area Increment (Inside Out)

bai.in(rwl, d2pith = NULL)

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

Value

Adata.frame containing the ring areas for each series withcolumn names, row names and dimensions ofrwl.

Note

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

Author(s)

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

References

Biondi, F. (1999) Comparing tree-ring chronologies and repeated timberinventories as forest monitoring tools.EcologicalApplications,9(1), 216–227.

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

See Also

bai.out

Examples

library(graphics)library(stats)library(utils)## Toyn <- 100## Make three fake tree-ring series to show that these funcs work on rwl objectsbase.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 methodfoo <- bai.in(rwl = rwl)## The outside in methodbar <- bai.out(rwl = rwl)## Identicalhead(bar)head(foo)## Use gp datadata(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)

Basal Area Increment (Outside In)

Description

Convert multiple ring-width series to basal area increment (i.e., ringarea) 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 thewidth of each ring moving towards the pith of the tree. It is relatedtobai.in, which calculates each ring area starting fromthe inside of the tree and working outward. Both methods assume acircular cross section (Biondi 1999). See the references below forfurther details.

Value

Adata.frame containing the ring areas for each series withcolumn names, row names and dimensions ofrwl.

Note

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

Author(s)

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

References

Biondi, F. (1999) Comparing tree-ring chronologies and repeated timberinventories as forest monitoring tools.EcologicalApplications,9(1), 216–227.

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

See Also

bai.in

Examples

library(graphics)library(utils)## Not run: library(stats)## Toyn <- 100## Make three fake tree-ring series to show that these funcs work on rwl objectsbase.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 methodfoo <- bai.in(rwl = rwl)## The outside in methodbar <- bai.out(rwl = rwl)## Identicalhead(bar)head(foo)## End(Not run)## Use gp datadata(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]) * 10bar <- 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)

Basal Area Increment (Bakker)

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 preventsbai.out transformations from producing negative increments when the sum of all ring widths in a series is larger than DBH/2. It preventsbai.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 baiyrs <- time(zof.rwl)plot(yrs,zof.bai[,1],type="l",     xlab="Year",     ylab=expression(BAI~(mm^2)),     main = colnames(zof.bai)[1])

Campito Mountain Tree Ring Widths

Description

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

Usage

data(ca533)

Format

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

Source

International tree-ring data bank, Accessed on 20-April-2021 athttps://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 DataSet.IGBPPAGES/World Data Center forPaleoclimatology Data Contribution Series 1983-CA533.RWL.NOAA/NCDC Paleoclimatology Program, Boulder,Colorado,USA.

Twisted Tree Heartrot Hill Standard Chronology

This data set gives the standard chronology for white sprucePicea glauca at Twisted Tree Heartrot Hill in Yukon,Canada. Data set was created usingread.crn and saved toan .rda file usingsave.

Usage

data(cana157)

Format

Adata.frame containing the standard chronology in column oneand 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 athttps://www.ncei.noaa.gov/pub/data/paleo/treering/chronologies/northamerica/canada/cana157.crn

References

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

Cook and Peters Smoothing Spline with User-Specified Rigidity and Frequency Cutoff

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

Cook, E. R. and Kairiukstis, L. A., editors (1990)Methods ofDendrochronology: 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 setdata(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)# A 2/3 spline, the default, 462.6667 yrs herelines(caps(series), col = "blue", lwd = 2)legend("topright",       c("Series", "nyrs=10", "nyrs=100",         paste("Default nyrs (", floor(length(series) * 2/3), ")", sep="")),       fill=c("grey", "red", "green", "blue"))## Note behavior when NA is encountered and## take appropriate measures as demonstrated abovey <- c(NA,NA,rnorm(100))caps(y)

Cross-Correlation between a Series and a Master Chronology

Description

Computes cross-correlations between a tree-ring series and a masterchronology 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,               prewhiten = TRUE, 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 atree-ring series and a master chronology built fromrwllooking at correlations lagged positively and negatively usingccf at overlapping segments set byseg.length. For instance, withlag.max setto 5, cross-correlations would be calculated at for each segment withthe master lagged atk = -5:5 years.

The cross correlations are calculated callingccf as
ccf(x=master, y=series, lag.max=lag.max, plot=FALSE) ifseries.x isFALSE and asccf(x=series, y=master, lag.max=lag.max, plot=FALSE) ifseries.x isTRUE. This argument was introduced in dplR version 1.7.0.Different users have different expectations about how missing or extra rings are notated. Ifswitch.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 thesecond segment and so on. Correlations are only calculated for segments withcomplete overlap with the master chronology.

Each series (including those in therwl object) isoptionally detrended as the residuals from ahanningfilter with weightn. The filter is not applied ifn isNULL. Detrending can also be done viaprewhitening where the residuals of anar model areadded to each series mean. This is the default. The master chronologyis computed as the mean of therwl object usingtbrm ifbiweight isTRUE androwMeans if not. Note that detrending typically changes thelength of the series. E.g., ahanning filter willshorten the series on either end byfloor(n/2). Theprewhitening default will change the series length based on thear model fit. The effects of detrending can be seen withseries.rwl.plot.

Value

Alist containing matricesccf andbins. Matrixccf contains the correlationsbetween the series and the master chronology at the lags window givenbylag.max. Matrixbins contains the yearsencapsulated 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 by deleting a year of growth in a random seriesflagged <- dat$"641143"flagged <- c(NA, flagged[-325])names(flagged) <- rownames(dat)dat$"641143" <- NULLccf.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 positionccf.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) # TRUEidentical(ccf.100.2, ccf.100.3) # TRUE## End(Not run)

Build Mean Value Chronology

Description

This function builds a mean value chronology, typically from adata.frame of detrended ring widths as produced bydetrend.

Usage

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

Arguments

Details

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

Note that the residual chronology in this function will return differentvalues than the residual chronology fromchron.ars which usesa slightly different method for determining AR order.

Value

An object of of classcrn anddata.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 ofDendrochronology: 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,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 chronca533.crn2 <- chron(ca533.rwi, prewhiten = TRUE)plot(ca533.crn2,xlab="Year",ylab="RWI")

Build ARSTAN Chronology

Description

This function builds three varieties of the mean-value chronology, including the ARSTAN chronology, typically from adata.frame of detrended ring widths as produced bydetrend.

Usage

chron.ars(x, 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 tochron.

The residual chronology is the prewhitened chronology as described by Cook (1985) and uses uses multivariate autoregressive modeling to determine the order of the AR process. It's important to note that residual chronology produced here is different than the simple residual chronology produced bychron which returns the residuals of an AR process using a naive call toar. But 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 builds on the residual chronology but returns a re-whitened chronology where the pooled AR coefficients from the multivariate autoregressive modeling are reintroduced. 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 otherwise or the overall minimum considering the maximum lag (argumentmaxLag).

Once the AR order is determined an AR(p) model is fit using eitherar via the Yule-Walker method or byarima 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 viaprewhitenMethod and defaults to usingar.

Value

Adata.frame with the standard, residual, and ARSTAN chronologies. The sample depth is also included.

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 ofDendrochronology: 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)

Build Mean Value Chronology with Confidence Intervals

Description

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

Usage

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

Arguments

Details

This either averages the rows of thedata.frame using a mean ora 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 inx 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 classdata.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 fivewa082Trunc <- wa082[rowSums(!is.na(wa082))>4,]# detrendwa082RWI <- detrend(wa082Trunc, method="AgeDepSpline")# bootstrap the chronology andwa082Crn <- 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)

Build Mean Value Chronology with Stabilized Variance

Description

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

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 functionchron.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)

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

Value

An object of of classcrn anddata.frame with the variance stabilized chronology, running interseries correlation ('ifrunning.bar=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")

C-Method Standardization

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 calculatinga growth curve based on constant annual basal area increment. Themethod is based on the “assumption that constant growth isexpressed by a constant basal area increment distributed over agrowing surface” (Biondi and Qeadan 2008). The detrending is theestimation and removal of the tree’s natural biologicalgrowth trend. The standardization is done by dividing each series bythe growth trend to produce units in the dimensionless ring-widthindex (RWI).

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

See the reference below for further details.

Value

Adata.frame containing the dimensionless and detrendedring-width indices with column names, row names and dimensions ofrwl ifc.hat.t isFALSE andc.hat.i isFALSE.

Otherwise alist of length 2 or 3 containing theRWIdata.frame, adata.frame containing the C-curves foreach tree (c.hat.t), and/or a vector containing theC-values for each tree (c.hat.i) depending on the outputflags. See Eq. 12 in Biondi and Qeadan (2008) for more detail onc.hat.t, andc.hat.i.

Note

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

Author(s)

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

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ringstandardization: Defining the biological trend from expected basal areaincrement.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.hatgp.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.tplot(tmp[, 1], type = "n", ylim = range(tmp, na.rm = TRUE),     xlab = "Cambial Age", ylab = expression(hat(c)[t]))apply(tmp, 2, lines)

Schulman Old Tree No. 1, Mesa Verde

Description

This data set gives the raw ring widths for Douglas firPseudotsuga menziesii at Mesa Verde in Colorado,USA.There are 35 series. Data set was created usingread.rwland saved to an .rda file usingsave.

Usage

data(co021)

Format

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

Source

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

References

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

Combine Tree-Ring Data Sets

combine.rwl(x, y = NULL)

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

This function finds the common interval on a set of tree-ring widthssuch as that produced byread.rwl.

Usage

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

Arguments

Details

This trims anrwl object to a common interval that maximizesthe number of series (type="series"), the number of years(type="years"), or a compromise between the two(type="both"). A modifiedseg.plot can be drawnas well.

Value

Adata.frame withcolnames(x) andrownames(x).

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 seriesprod(dim.s) #   (33 series x 288 years = 9504)dim.y <- dim(co021.y)dim.y       # the highest number of yearsprod(dim.y) #   (27 series x 458 years = 12366)dim.b <- dim(co021.b)dim.b       # compromise solutionprod(dim.b) #   (28 series x 435 years = 12180)

Compute Correlations between Series

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,             prewhiten = TRUE, pcrit = 0.05, biweight = TRUE,             method = c("spearman", "pearson","kendall"),             make.plot = TRUE, label.cex = 1, floor.plus1 = FALSE,             master = NULL,             master.yrs = as.numeric(if (is.null(dim(master))) {                              names(master)                          } else {                              rownames(master)                          }),             ...)

Arguments

Details

This function calculates correlation serially between each tree-ringseries and a master chronology built from all the other series in therwl object (leave-one-out principle). Optionally, theuser may give amaster chronology (avector) as anargument. In the latter case, the same master chronology is used forall the series in therwl object. The user can alsochoose to give amasterdata.frame (series ascolumns, years as rows), from which a single master chronology isbuilt.

Correlations are done for each segment of the series where segmentsare lagged by half the segment length (e.g., 100-year segments wouldbe overlapped by 50-years). The first segment is placed according tobin.floor. The minimum bin year is calculated asceiling(min.yr/bin.floor)*bin.floor wheremin.yr is the first year in either therwlobject or the user-specifiedmaster chronology, whicheveris smaller. For example if the first year is 626 andbin.floor is 100 then the first bin would start in 700.Ifbin.floor is 10 then the first bin would start in 630.

Correlations are calculated for the first segment, then the secondsegment and so on. Correlations are only calculated for segments withcomplete overlap with the master chronology. For now, correlations areSpearman’s rho calculated viacor.test usingmethod = "spearman".

Each series (including those in the rwl object) is optionallydetrended as the residuals from ahanning filter withweightn. The filter is not applied ifn isNULL. Detrending can also be done via prewhitening where theresiduals of anar model are added to each seriesmean. This is the default. The master chronology is computed as themean of therwl object usingtbrm ifbiweight isTRUE androwMeans if not. Notethat detrending can change the length of the series. E.g., ahanning filter will shorten the series on either end byfloor(n/2). The prewhitening default will change theseries length based on thear model fit. The effects ofdetrending can be seen withseries.rwl.plot.

The function is typically invoked to produce a plot where each segmentfor each series is colored by its correlation to the masterchronology. Green segments are those that do not overlap completelywith the width of the bin. Blue segments are those that correlateabove the user-specified critical value. Red segments are those thatcorrelate below the user-specified critical value and might indicate adating problem.

Value

Alist containing matricesspearman.rho,p.val,overall,bins,rwi, vectoravg.seg.rho,numericseg.lag,seg.length,pcrit,label.cex. An additionalcharacterflags is also returned if any segments fall below thecritical value. Matrixspearman.rho contains thecorrelations for each series by bin. Matrixp.valcontains the p-values on the correlation for each series bybin. Matrixoverall contains the average correlation andp-value for each series. Matrixbins contains the yearsencapsulated by each bin. The vectoravg.seg.rhocontains the average correlation for each bin. Matrixrwicontains the detrended rwl data, the numericsseg.lag,seg.length,pcrit,label.cex are from the oroginal call and used to pass intoplot.crs.

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 serieshead(crs$overall)## Average correlation for each bincrs$avg.seg.rho

Compute Correlation between a Series and a Master Chronology

Description

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

Usage

corr.series.seg(rwl, series, series.yrs = as.numeric(names(series)),                seg.length = 50, bin.floor = 100, n = NULL,                prewhiten = TRUE, 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 amaster chronology built from a rwl object. Correlations are done bysegment (see below) and with a moving correlation with length equal totheseg.length. The function is typically invoked toproduce a plot.

Value

Alist containing matricesbins,moving.rho, and vectorsspearman.rho,p.val, andoverall.

Matrixbins contains the years encapsulated by each bin(segments). Matrixmoving.rho contains the movingcorrelation and p-value for a moving average equal toseg.length. Vectorspearman.rho containsthe correlations by bin andp.val containsthe p-values. Vectoroverall contains the averagecorrelation 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 by deleting a year of growth in a random seriesflagged <- dat$"641143"flagged <- c(NA, flagged[-325])names(flagged) <- rownames(dat)dat$"641143" <- NULLseg.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 positionseg.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) # TRUEidentical(seg.100.2, seg.100.3) # TRUE## End(Not run)

Read Ring Width File from CSV

Description

This function reads in a file of ring widths (.rwl) from a text file with comma separated values (csv).

Usage

csv2rwl(fname,...)

Arguments

Details

This is a simple wrapper toread.table that reads in a text file with ring-width data in "spreadsheet" format. I.e., with series in columns and the the years as rows. The first column should contain the years and each subsequent column should contain a tree-ring series. The series names should be in the first row of the file. The deafult forNA values are empty cells or as the character string"NA" but can also be set using thena.strings argument passed toread.table. E.g.,:

Note that this is a rudimentary convenience function that isn't doing anything sophisticated. It reads in a file, assigns the years to the row names and sets the class of the object toc("rwl","data.frame") which allowsdplR to recognize it.

Although arguments can be passed toread.table, this is not designed to be a flexible function. As is the philosophy withdplR, modifying the code is easy should the user wish to read in a different style of text file (e.g., tab delimited). Typingcsv2rwl at theR prompt will get the user started.

Value

Function returns an object of classc("rwl", "data.frame") with the series in columns and the yearsas rows. The seriesIDs are the column names and the yearsare the row names.

Author(s)

Andy Bunn

See Also

read.rwl,read.table

Examples

library(utils)data(ca533)# write out a rwl file in a format that csv2rwl will understandtm <- time(ca533)foo <- data.frame(tm,ca533)# this is the temp file where foo will be writtentmpName <- tempfile()write.csv(foo,file=tmpName,row.names=FALSE)# read it back in using csv2rwlbar <- csv2rwl(tmpName)# check to see if identicalidentical(ca533,bar)# delete temp fileunlink(tmpName)

Detrend Multiple Ring-Width Series Simultaneously

Description

This is a wrapper fordetrend.series to detrend manyring-width series at once.

Usage

detrend(rwl, y.name = names(rwl), make.plot = FALSE,        method = c("Spline", "ModNegExp", "Mean", "Ar", "Friedman",                   "ModHugershoff", "AgeDepSpline"),        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

Seedetrend.series for details on detrendingmethods. Settingmake.plot = TRUE will cause plots ofeach series to be produced. These could be saved usingDevices if desired.

Value

If one detrending method is used, adata.frame containing thedimensionless detrended ring widths with column names, row names anddimensions ofrwl. If more methods are used, a list withncol(rwl) elements each containing adata.framewith the detrended ring widths in each column.

Ifreturn.info isTRUE, the return value is alist with four parts:

Note

This function uses theforeach loopingconstruct with the%dopar% operator.For parallel computing and a potential speedup, a parallel backendmust be registered before running the function. Ifverbose isTRUE, parallel computation is disabled.

Author(s)

Andy Bunn. Improved by Mikko Korpela.

See Also

detrend.series

Examples

library(utils)data(ca533)## Detrend using modified exponential decay. Returns a data.frameca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")## Detrend using splines and compute## residuals via subtractionca533.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.infodata(co021)co021.rwi <- detrend(rwl = co021, method = "ModHugershoff",                     return.info=TRUE)## Not run: library(grDevices)## Detrend using all methods. Returns a listca533.rwi <- detrend(rwl = ca533)## Save a pdf of all seriesfname <- tempfile(fileext=".pdf")print(fname) # tempfile used for outputpdf(fname)ca533.rwi <- detrend(rwl = ca533, method = c("Spline", "ModNegExp"),                     make.plot = TRUE)dev.off()unlink(fname) # remove the file## End(Not run)

Detrend a Ring-Width Series

Description

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

Usage

detrend.series(y, y.name = "", make.plot = TRUE,               method = c("Spline", "ModNegExp", "Mean",                "Ar", "Friedman", "ModHugershoff","AgeDepSpline"),               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). Ifdifference is TRUE, the index is calculated by subtraction. Values of zero (typically missing rings) iny are set to 0.001.

There are currently seven methods available fordetrending although more are certainly possible. The methodsimplemented are an age-dependent spline viaads (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 Hugershoffcurve (method = "ModHugershoff"), a modified negative exponentialcurve (method = "ModNegExp"), and a smoothing spline viacaps (method = "Spline").

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

The"Ar" approach is also known as prewhitening where the detrended series is the residuals of anar model divided by themean 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 seriesand should only be used as such.

The"Friedman" approach uses Friedman’s ‘supersmoother’ as implemented insupsmu. The parameterswt,span andbass can beadjusted, butperiodic is always set toFALSE. If some of the fitted values are not positive then method"Mean" is used.

The"Mean" approach fits a horizontal line using the mean ofthe 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 invalidring-width indices.

The"ModHugershoff" approach attempts to fit a Hugershoffmodel of biological growth of the formf(t) = a t^b e^{-g t} + d, where the argument of the function is time, usingnls. See Fritts (2001) for details about theparameters. Optionconstrain.nls gives apossibility to constrain the parameters of the modified negativeexponential function. If the constraints are enabled, the nonlinearoptimization algorithm is instructed to keep the parameters in thefollowing ranges:a \ge 0,b \ge 0 andd \ge 0. The default is to not constrain the parameters(constrain.nls = "never") fornls butwarn 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 alinear model is fit. That linear model can have a positive slopeunlesspos.slope isFALSE in which case method"Mean" is used.

The"ModNegExp" approach attempts to fit a classic nonlinearmodel of biological growth of the formf(t) = a e^{b t} + k, where the argument of the function is time, usingnls. See Fritts (2001) for details about theparameters. Optionconstrain.nls gives apossibility to constrain the parameters of the modified negativeexponential function. If the constraints are enabled, the nonlinearoptimization algorithm is instructed to keep the parameters in thefollowing ranges:a \ge 0,b \le 0 andk \ge 0. The default is to not constrain the parameters(constrain.nls = "never") fornls butwarn 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 alinear model is fit. That linear model can have a positive slopeunlesspos.slope isFALSE in which case method"Mean" is used.

The"Spline" approach uses a spline where the frequencyresponse is 0.50 at a wavelength of 0.67 * “series length inyears”, unless specified differently usingnyrs andf in the functioncaps. If some of the fitted values are not positive then method"Mean" is used.

These methods are chosen because they are commonly used indendrochronology. There is a rich literature on detrendingand 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’ butisn't unheard of (see ‘Examples’). If this happens, users should lookcarefully 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 adata.frame containingthe detrended series (y) according to the methods used.The columns are named and ordered to matchmethod. Ifonly one method is selected, returns a vector.

Ifreturn.info isTRUE, the return value is alist with four parts:

Author(s)

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

References

Cook, E. R. and Kairiukstis, L. A., editors (1990)Methods ofDendrochronology: 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 setdata(ca533)series <- ca533[, "CAM011"]names(series) <- rownames(ca533)# defaults to all six methodsseries.rwi <- detrend.series(y = series, y.name = "CAM011", verbose=TRUE)# see plot with three methodsseries.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 growthseries.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.infodata(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)

Defunct functions in dplR

Description

These are functions no longer included in dplR

Details

Fill Internal NA

Description

This function fills internalNA values (i.e., those with numeric dataabove and below a small data gap) in each column of adata.frame such as a data set of ring widths as produced byread.rwl.

Usage

fill.internal.NA(x, fill = c("Mean", "Spline", "Linear"))

Arguments

Details

There are occasionally data gaps within a tree-ring series. Some ofthe functions in dplR will failwhen an internalNA is encountered (e.g.caps). Thisfunction fills internalNA 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 thespline function or calculatedby linear approximation (fill="Linear") using the functionapprox.

Editorial: Having internalNA in a tree-ring series isoften bad practice and filling those values should be done withcaution. For instance, some users code missing rings asNAinstead of0. And missing values (i.e.,NA) aresometimes present in maximum latewood density data when the rings aresmall. A common, but not recommended, practice is to leave stretchesofNA values in places where it has been impossible toaccurately measure rings (perhaps because of a break in the core). Itis often better to treat that core as two separate series (e.g., "01A"and "01B" rather than have internalNA values. As with allprocessing, the analyst should make a decision based on theirexperience with the wood and not rely on software to make a choice forthem!

Value

Adata.frame withcolnames(x) andrownames(x). InternalNAsfilled as above.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

spline,approx

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:1910class(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)

Calculate the Gini Coefficient

Description

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

Usage

gini.coef(x)

Arguments

Details

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

Value

the Gini coefficient.

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 DistributionMethods 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)

Calculate Gleichläufigkeit

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 theoverlap orprob and returned amatrix with just the Gleichläufigkeit for all possible combinations of records. That function can still be accessed viaglk.legacy.

Value

The funtions returns a namedlist of two or three matrices (p_mat is optional ifprob = 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, andNA is returned.To calculate the global glk of the datasetmean(x$glk_mat, na.rm = TRUE).

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,Dendrochronologia34, 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 Werkstoff6, 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

sgcsgc is an alternative forglk)

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 togp.rwl

Description

This data set gives the distance to pith for each series (in mm) thatmatches the ring widths forgp.rwl – a data set ofponderosa pine (Pinus ponderosa) from the Gus Pearson NaturalArea (GPNA) in northern Arizona,USA. Data arefurther described by Biondi and Qeadan (2008) and references therein.

Usage

data(gp.d2pith)

Format

Adata.frame containing seriesIDs in column 1(series) and the distance (in mm) from the innermost ringto the pith of the tree (d2pith). This can be usedtogether 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-ringstandardization: Defining the biological trend from expected basal areaincrement.Tree-Ring Research,64(2), 81–96.

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

Description

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

Usage

data(gp.dbh)

Format

Adata.frame containing seriesIDs in column 1(series), tree diameter (in cm) at breast height(dbh), and the bark thickness (in cm). This can be usedtogether 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-ringstandardization: Defining the biological trend from expected basal areaincrement.Tree-Ring Research,64(2), 81–96.

Ponderosa Pine Pith Offsets Corresponding togp.rwl

Description

This data set gives the pith offsets that match the ring widths forgp.rwl – a data set of ponderosa pine (Pinusponderosa) from the Gus Pearson Natural Area (GPNA) innorthern Arizona,USA. Data are further described by Biondiand Qeadan (2008) and references therein.

Usage

data(gp.po)

Format

Adata.frame containing seriesIDs in column 1(series) and the number of years between the beginning ofthat series ingp.rwl and the pith of the tree(pith.offset). This can be used together with the ringwidths to calculate the cambial age 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-ringstandardization: Defining the biological trend from expected basal areaincrement.Tree-Ring Research,64(2), 81–96.

Ponderosa Pine Ring Widths from Gus Pearson Natural Area

Description

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

Usage

data(gp.rwl)

Format

Adata.frame containing 58 ring-width series in columns and 421years in rows.

Source

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

References

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

Hanning Filter

Description

Applies a Hanning filter of lengthn tox.

Usage

hanning(x, n = 7)

Arguments

Details

This applies a low frequency Hanning (a.k.a. Hann) filter tox with weight set ton.

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"))

Interactively Detrend Multiple Ring-Width Series

Description

Interactively detrend multiple tree-ring series by one of two methods, asmoothing spline or a statistical model. This is a wrapper fordetrend.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 plotsthat are produced bydetrend.series for which it isessentially a wrapper. The user enters their choice of detrendedmethod via keyboard at a prompt for each ring width series inrwl. Seedetrend.series for examples anddetails on the detrending methods.

Value

Adata.frame containing each detrended series according to themethod used as columns and rownames set tocolnames(y). These are typically years. Plots are alsoproduced as the user chooses the detrending methods through keyboardinput.

Author(s)

Andy Bunn

See Also

detrend.series

Interactively Detrend a Ring-Width Series

Description

Interactively detrend a tree-ring series by one of three methods, asmoothing spline, a linear model, or the mean. This is a wrapper fordetrend.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 aplot that is produced bydetrend.series for which it isessentially a wrapper. The user enters their choice of detrendedmethod via keyboard at a prompt. Seedetrend.series forexamples and details on the detrending methods.

Value

A vector containing the detrended series (y) according tothe method used with names set tocolnames(y). These aretypically years. A plot is also produced and the user chooses a methodthrough keyboard input.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

detrend.series

Edit a Ring-Width Series

Description

Insert or delete rings from a ring-width series

Usage

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)delete.ring(rw.vec,rw.vec.yrs=as.numeric(names(rw.vec)),            year,fix.last=TRUE,fix.length=TRUE)

Arguments

Details

Simple editing of ring widths.

Value

A named vector.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

dplR

Examples

library(utils)data(gp.rwl)dat <- gp.rwl# insert a value of zero for the year 1950 in series 50A# fix the last year of growth and maintain the length of the seriestmp <- insert.ring(rw.vec=dat$"50A",rw.vec.yrs = time(dat),                   year=1950,ring.value=0,fix.length = TRUE)# with fix.length=TRUE this can be merged back into the rwl object:data.frame(dat$"50A",tmp)dat$"50A" <- tmp# note that if fix.last = FALSE and fix.length = FALSE inserting a ring causes the# ending year of the series to be pushed forward and the length of the output to# be longer than the original series.tmp <- insert.ring(rw.vec=dat$"50A",rw.vec.yrs = time(dat),                   year=1950,ring.value=0, fix.last = FALSE,                    fix.length = FALSE)# with fix.length=FALSE this can't be merged back into the rwl object the# same way as abovetail(tmp)length(tmp)nrow(dat)# the same logic applies to deleting rings.dat <- gp.rwl# delete the year 1950 in series 50A# fix the last year of growth and maintain the length of the seriestmp <- delete.ring(rw.vec=dat$"50A",rw.vec.yrs = time(dat),                   year=1950,fix.last = TRUE, fix.length = TRUE)# with fix.length=TRUE this can be merged back into the rwl object:data.frame(dat$"50A",tmp)dat$"50A" <- tmp# note that if fix.last = FALSE and fix.length = FALSE inserting a ring causes the# ending year of the series to be pushed forward and the length of the output to# be longer than the original series.tmp <- delete.ring(rw.vec=dat$"50A", rw.vec.yrs = time(dat),                   year=1950, fix.last = FALSE,                    fix.length = FALSE)# with fix.length=FALSE this can't be merged back into the rwl object the# same way as abovetail(tmp)length(tmp)nrow(dat)

Individual Series Correlation Against a Master Chronology

Description

This function calculates the correlation between a series and a masterchronology.

Usage

  interseries.cor(rwl,n=NULL,prewhiten=TRUE,biweight=TRUE,    method = c("spearman", "pearson", "kendall"))

Arguments

Details

This function calculates correlation serially between each tree-ringseries and a master chronology built from all the other series in therwl object (leave-one-out principle).

Each series in the rwl object is optionallydetrended as the residuals from ahanning filter withweightn. The filter is not applied ifn isNULL. Detrending can also be done via prewhitening where theresiduals of anar model are added to each seriesmean. This is the default. The master chronology is computed as themean of therwl object usingtbrm ifbiweight isTRUE androwMeans if not. Notethat detrending can change the length of the series. E.g., ahanning filter will shorten the series on either end byfloor(n/2). The prewhitening default will change theseries length based on thear model fit. The effects ofdetrending can be seen withseries.rwl.plot.

This function produces the same output of theoverall portion ofcorr.rwl.seg. The mean correlation value given is sometimes referred to as the “overall interseries correlation” or the “COFECHA interseries correlation”. This output differs from therbar statistics given byrwi.stats in thatrbar is the average pairwise correlation between series where this is the correlation between a series and a master chronology.

Value

adata.frame with correlation values and p-values given fromcor.test

Author(s)

Andy Bunn, patched and improved by Mikko Korpela

See Also

rwl.stats,rwi.stats

Examples

library(utils)data(gp.rwl)foo <- interseries.cor(gp.rwl)# compare to: # corr.rwl.seg(rwl=gp.rwl,make.plot=FALSE)$overall# using pearson's rfoo <- interseries.cor(gp.rwl,method="pearson")# two measures of interseries correlation# compare interseries.cor to rbar from rwi.statsgp.ids <- read.ids(gp.rwl, stc = c(0, 2, 1))bar <- rwi.stats(gp.rwl, gp.ids, prewhiten=TRUE)bar$rbar.effmean(foo[,1])

Date Conversion to Character in LaTeX Format

Description

This is a simple convenience function that returns a date in theformat used by ‘⁠\today⁠’ in LaTeX. A possible use case is fixingthe date shown in a vignette at weaving time.

Usage

latexDate(x = Sys.Date(), ...)

Arguments

Value

Acharacter vector

Author(s)

Mikko Korpela

Examples

latexDate()                              # todaylatexDate(Sys.Date() + 5)                # today + 5 dayslatexDate(c("2013-12-06", "2014-09-19")) # fixed dates## [1] "December 6, 2013"   "September 19, 2014"latexDate(5*60*60*24, origin=Sys.Date()) # today + 5 days

Convert Character Strings for Use with LaTeX

Description

Some characters cannot be entered directly into a LaTeX document.This function converts the inputcharacter vector to a formsuitable for inclusion in a LaTeX document in text mode. It can beused together with ‘⁠\Sexpr⁠’ in vignettes.

Usage

latexify(x, doublebackslash = TRUE, dashdash = TRUE,         quotes = c("straight", "curved"),         packages = c("fontenc", "textcomp"))

Arguments

Details

The function is intended for use with unformatted inline text.Newlines, tabs and other whitespace characters ("[:space:]" inregex) are converted to spaces. Control characters("[:cntrl:]") that are not whitespace are removed. Other moreor less special characters in theASCII set are ‘{’,‘}’, ‘\’, ‘#’, ‘$’, ‘%’,‘^’, ‘&’, ‘_’, ‘~’, double quote,‘/’, single quote, ‘<’, ‘>’, ‘|’, graveand ‘-’. They are converted to the corresponding LaTeXcommands. Some of the conversions are affected by user options,e.g.dashdash.

Before applying the substitutions described above, input elements withEncoding set to"bytes" are printed and theoutput is stored usingcaptureOutput. The result ofthis intermediate stage isASCII text where some charactersare shown as their byte codes using a hexadecimal pair prefixed with"\x". This set includes tabs, newlines and controlcharacters. The substitutions are then applied to the intermediateresult.

The quoting functionssQuote anddQuotemay use non-ASCII quote characters, depending on the locale.Also these quotes are converted to LaTeX commands. This means thatthe quoting functions are safe to use with any LaTeX input encoding.Similarly, some other non-ASCII characters, e.g. letters,currency symbols, punctuation marks and diacritics, are converted tocommands.

Adding"eurosym" topackages enables the use of theeuro sign as provided by the"eurosym" package (‘⁠\euro⁠’).

The result is converted to UTF-8 encoding, Normalization Form C (NFC).Note that this function will not add any non-ASCIIcharacters that were not already present in the input. On thecontrary, some non-ASCII characters, e.g. all characters inthe"latin1" (ISO-8859-1)Encoding(character set), are removed when converted to LaTeX commands. Anyremaining non-ASCII character has a good chance of workingwhen the document is processed with XeTeX or LuaTeX, but the Unicodesupport available with pdfTeX is limited.

Assuming that ‘⁠pdflatex⁠’ is used for compilation, suggestedpackage loading commands in the document preamble are:

\usepackage[T1]{fontenc}    % no '"' in OT1 font encoding\usepackage{textcomp}       % some symbols e.g. straight single quote\usepackage[utf8]{inputenx} % UTF-8 input encoding\input{ix-utf8enc.dfu}      % more supported characters

Value

Acharacter vector

Author(s)

Mikko Korpela

References

INRIA. Tralics: a LaTeX to XML translator, HTML documentation of allTeX commands.https://www-sop.inria.fr/marelle/tralics/.

Levitt, N., Persch, C., and Unicode, Inc. (2013) GNOME Character Map,application version 3.10.1.

Mittelbach, F., Goossens, M., Braams, J., Carlisle, D., andRowley, C. (2004)The LaTeX Companion. Addison-Wesley,second edition.ISBN-13: 978-0-201-36299-2.

Pakin, S. (2009) The Comprehensive LaTeX SymbolList.https://www.ctan.org/tex-archive/info/symbols/comprehensive.

The Unicode Consortium. The Unicode Standard.https://home.unicode.org/.

Examples

x1 <- "clich\xe9\nma\xf1ana"Encoding(x1) <- "latin1"x1x2 <- x1Encoding(x2) <- "bytes"x2x3 <- enc2utf8(x1)testStrings <-    c("different     kinds\nof\tspace",      "control\a characters \ftoo",      "{braces} and \\backslash",      '#various$ %other^ &characters_ ~escaped"/coded',      x1,      x2,      x3)latexStrings <- latexify(testStrings, doublebackslash = FALSE)## All should be "unknown"Encoding(latexStrings)cat(latexStrings, sep="\n")## Input encoding does not matteridentical(latexStrings[5], latexStrings[7])

Perform a Continuous Morlet Wavelet Transform

Description

This function performs a continuous wavelet transform on a time series.

Usage

morlet(y1, x1 = seq_along(y1), p2 = NULL, dj = 0.25, siglvl = 0.95)

Arguments

Details

This performs a continuous wavelet transform of a time series. Thisfunction is typically invoked withwavelet.plot.

Value

Alist containing:

Note

This is a port of Torrence’sIDL code,which can be accessed through theInternet Archive Wayback Machine.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

References

Torrence, C. and Compo, G. P. (1998) A practical guide to waveletanalysis.Bulletin of the American Meteorological Society,79(1), 61–78.

See Also

wavelet.plot

Examples

library(utils)data(ca533)ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")ca533.crn <- chron(ca533.rwi, prewhiten = FALSE)Years <- time(ca533.crn)CAMstd <- ca533.crn[, 1]out.wave <- morlet(y1 = CAMstd, x1 = Years, dj = 0.1, siglvl = 0.99)

Calculate NET

Description

Computes the\mathit{NET} parameter for a set of tree-ringrecords or other time-series data.

Usage

net(x, weights = c(v = 1, g = 1))

Arguments

Details

This function computes the\mathit{NET} parameter (Esper etal., 2001). The overall\mathit{NET} is an average of all(non-NA) yearly values\mathit{NET_j}, which arecomputed as follows:

\mathit{NET_j}=v_j+(1-G_j)

The yearly variationv_j is the standard deviation of themeasurements of a single year divided by their mean.Gegenläufigkeit1-G_j is basedon one definition of GleichläufigkeitG_j, similar to but not the same as whatglkcomputes. Particularly, in the formula used by this function (Esperet al., 2001), simultaneous zero differences in two series are notcounted as a synchronous change.

The weights ofv_j and1-G_j in the sum canbe adjusted with the argumentweights (see above). As arather extreme example, it is possible to isolate variation orGegenläufigkeit by setting one of the weightsto zero (see ‘Examples’).

Value

Alist with the following components, in the same order asdescribed here:

Author(s)

Mikko Korpela

References

Esper, J., Neuwirth, B., and Treydte, K. (2001) A new parameter toevaluate temporal signal strength of tree-ring chronologies.Dendrochronologia,19(1), 93–102.

Examples

library(utils)data(ca533)ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")ca533.net <- net(ca533.rwi)tail(ca533.net$all)ca533.net$average## Not run: ## Isolate the components of NETca533.v <- net(ca533.rwi, weights=c(v=1,0))ca533.g <- net(ca533.rwi, weights=c(g=1,0))## End(Not run)

Los Alamos Tree Ring Widths

Description

This data set gives the raw ring widths for Douglas firPseudotsuga menziesii near Los Alamos New Mexico,USA. There are 8 series. Data set was created usingread.rwl and saved to an .rda file usingsave.

Usage

data(nm046)

Format

Adata.frame containing 8 tree-ring series in columns and 289years in rows.

Source

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

References

O'Brien, D (2002) Los Alamos DataSet.IGBPPAGES/World Data Center forPaleoclimatology Data Contribution Series 2002-NM046.RWL.NOAA/NCDC Paleoclimatology Program, Boulder,Colorado,USA.

Low-pass, high-pass, band-pass, and stop-pass filtering

pass.filt(y, W, type = c("low", "high", "stop", "pass"),          method = c("Butterworth", "ChebyshevI"),          n = 4, Rp = 1)

This function applies either a Butterworth or a Chebyshev type I filter of ordern to a signal and is nothing more than a wrapper for functions in thesignal package. The filters are designed viabutter andcheby1. The filter is applied viafiltfilt.

The input data (y) has the mean value subtracted and is then padded via reflection at the start and the end to a distance of twice the maximum period. The padded data and the filter are passed tofiltfilt after which the data are unpadded and returned afer the mean is added back.

The argumementW can be given in either frequency between 0 and 0.5 or, for convenience, period (minimum value of 2). For low-pass and high-pass filters,W must have a length of one. For low-pass and high-pass filtersW must be a two-element vector (c(low, high)) specifying the lower and upper boundaries of the filter.

Because this is just a wrapper for casual use with tree-ring data the frequencies and periods assume a sampling frequency of one. Users are encouraged to build their own filters using thesignal package.

Value

A filtered vector.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

hanning,detrend

Examples

data("co021")x <- na.omit(co021[, 1])# 20-year low-pass filter -- note freq is passed inbSm <- pass.filt(x, W=0.05, type="low", method="Butterworth")cSm <- pass.filt(x, W=0.05, type="low", method="ChebyshevI")plot(x, type="l", col="grey")lines(bSm, col="red")lines(cSm, col="blue")# 20-year high-pass filter -- note period is passed inbSm <- pass.filt(x, W=20, type="high")plot(x, type="l", col="grey")lines(bSm, col="red")# 20 to 100-year band-pass filter -- note freqs are passed inbSm <- pass.filt(x, W=c(0.01, 0.05), type="pass")cSm <- pass.filt(x, W=c(0.01, 0.05), type="pass", method="ChebyshevI")plot(x, type="l", col="grey")lines(bSm, col="red")lines(cSm, col="blue")# 20 to 100-year stop-pass filter -- note periods are passed incSm <- pass.filt(x, W=c(20, 100), type="stop", method="ChebyshevI")plot(x, type="l", col="grey")lines(cSm, col="red")

Plot a Tree-Ring Chronology

Description

This function makes a default plot of a tree-ring chronology from adata.frame of the type produced bychron,chron.ars,chron.stabilized,ssf.

Usage

## S3 method for class 'crn'plot(x, add.spline = FALSE, nyrs = NULL, ...)

Arguments

Details

This makes a crude plot of one or more tree-ring chronologies.

Value

None. Invoked for side effect (plot).

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

chron

Examples

library(graphics)data(wa082)# Truncate the RW data to a sample depth at least 5wa082Trunc <- wa082[rowSums(!is.na(wa082))>4,]# Detrend with age-dependent splinewa082RWI <- detrend(wa082Trunc,method = "AgeDep")# make several chronologieswa082CRN1 <- chron(wa082RWI)wa082CRN2 <- chron.stabilized(wa082RWI,                              winLength=51,                              biweight = TRUE,                              running.rbar = TRUE)wa082CRN3 <- chron.ars(wa082RWI)wa082CRN4 <- ssf(wa082Trunc)# and plotplot.crn(wa082CRN1,add.spline = TRUE,nyrs=20)plot.crn(wa082CRN2,add.spline = TRUE,nyrs=20)plot(wa082CRN3,add.spline = TRUE,nyrs=20)plot(wa082CRN4,add.spline = TRUE,nyrs=20)# a custom crnfoo <- data.frame(wa082CRN1,sfc=wa082CRN4$sfc)foo <- foo[,c(1,3,2)]class(foo) <- c("crn","data.frame")plot.crn(foo,add.spline = TRUE,nyrs=20)

Plotting crs Objects

Description

Plots objects returned fromcorr.rwl.seg.

Usage

## S3 method for class 'crs'plot(x, ...)

Arguments

Value

None. A plot is produced.

Author(s)

Andy Bunn

See Also

corr.rwl.seg

Examples

library(graphics)data(co021)foo <- corr.rwl.seg(co021, make.plot = FALSE)plot(foo)

Plotting Rwl Objects

Description

Plots rwl objects.

Usage

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

Arguments

Value

None. A plot is produced.

Author(s)

Andy Bunn

See Also

read.rwl

Examples

library(graphics)library(utils)data(co021)plot(co021, plot.type="seg")plot(co021, plot.type="spag")plot(co021, plot.type="spag", zfac=2)

Convert Pith Offset to Wood Completeness

Description

This function creates a partial wood completeness data structure basedon pith offset data.

Usage

po.to.wc(po)

Arguments

Details

Usespith.offset - 1 as the number of missing heartwoodrings.

Value

Adata.frame containing one variable of wood completeness data:n.missing.heartwood (integer type). This can beused as input towrite.tridas.

Author(s)

Mikko Korpela

See Also

wc.to.po,rcs,write.tridas

Examples

## Not run: library(utils)data(gp.po)all(wc.to.po(po.to.wc(gp.po)) == gp.po)## End(Not run)

Calculates Pointer Years from a Group of Ring-Width Series

Description

This function calculates pointer years on adata.frame ofring-width series using the Becker algorithm. The pointer years arecomputed with adjustable thresholds of relative radial growthvariation and number of series displaying similar growth pattern(i.e. positive or negative variations).

Usage

pointer(rwl, rgv.thresh = 10, nseries.thresh = 75, round.decimals = 2)

Arguments

Details

This calculates pointer years from ring-width series for each yeart of the time period covered by the series using theBecker algorithm. This algorithm is based on, first, the calculationof the individual relative radial growth variation by comparison ofring-width of yeart to that of yeart-1 foreach series, and second, the inter-series comparison of both sign andmagnitude of these variations.

For example, ifrgv.thresh andnseries.thresh are set at 10 and 75 respectively, pointeryears will be defined as those years when at least 75% of the seriespresent an absolute relative radial growth variation higher than 10%.

Users unfamiliar with the Becker algorithm should refer to Becker etal. (1994) and Mérian and Lebourgeois (2011) for furtherdetails.

Value

Adata.frame containing the following columns (each rowcorresponds to one position of the window):

Author(s)

Pierre Mérian. Improved by Mikko Korpela and AndyBunn.

References

Becker, M., Nieminen, T. M., and Gérémia, F. (1994)Short-term variations and long-term changes in oak productivity innortheastern France – the role of climate and atmosphericCO2.Annals of Forest Science,51(5),477–492.

Mérian, P. and Lebourgeois, F. (2011) Size-mediatedclimate–growth relationships in temperate forests: A multi-speciesanalysis.Forest Ecology and Management,261(8), 1382–1391.

See Also

skel.plot

Examples

## Pointer years calculation on ring-width series. Returns a data.frame.library(utils)data(gp.rwl)py <- pointer(rwl=gp.rwl, rgv.thresh=10, nseries.thresh=75,              round.decimals=2)tail(py)

Power Transformation of Tree-Ring Data

Description

Power transformation of tree-ring width.

Usage

powt(rwl, method = "universal", rescale = FALSE,     return.power=FALSE)

Arguments

Details

In dendrochronology, ring width series are sometimes power transformed to address heteroscedasticity.

The classic procedure used bymethod="cook" is a variance stabilization technique implemented afterCook & Peters (1997): for each series a linear model is fitted on thelogs of level and spread, where level is defined as the local meanM_t = \left(R_t + R_{t-1}\right)/2 withring widths R, and spread S is the local standard deviation defined asS_t = \left|R_t - R_{t-1}\right|. Theregression coefficient b from a linear model\log S = k + b \log M is then used for the power transform\star{R}_t = R_t^{1-b}.

The procedure above is modified withmethod="universal" where all samples are used simultaneously in a linear mixed-effects model with time (year)as a random effect:lmer(log S ~ log M + (1|year). This "universal" or "signal free" approach accounts for the common year effect across all of the series inrwl and should address that not every year has the same change in environmental conditions to the previous year.

Therescale argument will return the series with a mean and standarddeviation that matches the input data. While this is a common convention,users should note that this can produce negative values which can be confusingif thought of as "ring widths."

Value

Either an object of classc("rwl", "data.frame") containing the power transformed ring width series with the series in columns and the years as rows or in the case of a single series, a possibly named vector of the same. With classrwl, the seriesIDs are the column names and the years are the row names.

Ifreturn.power=TRUE the returned object is alist containing the power transformed data and anumeric with the power estimate(s) used to transform the data.

Author(s)

Christian Zang implemented the Cook and Peters method. Stefan Klesse conceivedand wrote the universal method. Patched and improved by Mikko Korpela and Andy Bunn.

References

Cook, E. R. and Peters, K. (1997) Calculating unbiasedtree-ring indices for the study of climatic and environmentalchange.The Holocene,7(3), 361–370.

Examples

library(utils)data(zof.rwl)powtUniversal <- powt(zof.rwl, method = "universal")powtCook <- powt(zof.rwl, method = "cook")op <- par(no.readonly = TRUE)par(mfcol = c(1, 3))hist(summary(zof.rwl)$skew,     breaks = seq(-2.25,2.25,by=0.25),     main="Raw Data",xlab="Skew")hist(summary(powtUniversal)$skew,     breaks = seq(-2.25,2.25,by=0.25),     main="Universal POWT",xlab="Skew")hist(summary(powtCook)$skew,     breaks = seq(-2.25,2.25,by=0.25),     main="Cook POWT",xlab="Skew")par(op) # restore graphical parameters

Printing Redfit Results

## S3 method for class 'redfit'print(x, digits = NULL, csv.out = FALSE, do.table = FALSE,      prefix = "", row.names = FALSE, file = "", ...)

This function is based on the Fortran programREDFIT,which is in the public domain.

Schulz, M. and Mudelsee, M. (2002) REDFIT: estimating red-noisespectra directly from unevenly spaced paleoclimatic time series.Computers & Geosciences,28(3), 421–426.

See Also

redfit

Examples

library(utils)data(ca533)tm <- time(ca533)x <- ca533[[1]]idx <- which(!is.na(x))redf <- redfit(x[idx], tm[idx], "time",               nsim = 100, iwin = 0, ofac = 1, n50 = 1)print(redf)fname <- tempfile(fileext=".csv")print(fname) # tempfile used for outputprint(redf, csv.out = TRUE, file = fname)redftable <- read.csv(fname)unlink(fname) # remove the file

Do some reporting on a RWL object

Description

This function prints the results ofrwl.report

Usage

## S3 method for class 'rwl.report'print(x, ...)

Arguments

Details

This function formats thelist fromrwl.report for the user to have a summary report of the number of series, the mean lengthof all the series, the first year, last year, the mean first-order autocorrelation (viasummary.rwl), the mean interseries correlation (viainterseries.cor), the years where a series has a missing ring (zero), internal NA, or a very small ring (<0.005).

Value

Invisible

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

rwl.report,summary.rwl,interseries.cor

Examples

data("gp.rwl")rwl.report(gp.rwl)foo <- gp.rwlfoo[177,1] <- NA foo[177:180,3] <- NA foo[185,4] <- 0.001 rwl.report(foo)

Add Raster Elements to Plot

Description

This function takes plotting commands and uses a temporarybitmap graphics device to capture their output. Theresulting raster image is drawn in the plot or figure region of theactive high-level plot. A new plot is started if one does not exist.

Usage

rasterPlot(expr, res = 150, region = c("plot", "figure"), antialias,           bg = "transparent", interpolate = TRUE, draw = TRUE,           Cairo = FALSE, ...)

Arguments

Details

The appropriate graphical parameters of the current graphics deviceare copied to the temporary bitmap device. Therefore theappearance of the raster contents should be almost the same as whendirectly drawn.

The call or expressionexpr is evaluated in theenvironment of the caller.

It is possible that the raster contents will maintain a constant sizewhen the graphics device is resized. If resizing works, however, theimage may become distorted. For example, circle symbols will turninto ellipses if the width to height ratio is not maintained (see‘Examples’). This is in contrast to a standard plot in adisplay graphics device, e.g.x11, where text andsymbols maintain their size when the device is resized.

Value

Ifdraw isTRUE, there is no return value. Thefunction is used for the side effects.

Ifdraw isFALSE, an object of class"nativeRaster" is returned. The object can be used as inputforrasterImage orgrid.raster. SeereadPNG. If no bitmap device is available(see ‘Note’),NULL is returned.

Note

• The graphics device used for the output must have support forincluding raster images. See"rasterImage" indev.capabilities.

• TheR build must have a functionalpng device,which requires one of the followingcapabilities:"png","aqua" or"cairo". Alternatively, aCairo device from package Cairo must beavailable withCairo.capabilities"raster" or"png".

If either of these requirements is not met, at least onemessage is generated and the function reverts to regularplotting. Thebg argument is then handled by drawing afilled rectangle. Alsoregion is honored, but the othersettings do not apply.

Author(s)

Mikko Korpela

Examples

library(graphics)library(stats)## Picture with various graphical elementsx <- 1:100y0 <- quote(sin(pi * x / 20) + x / 100 + rnorm(100, 0, 0.2))y <- eval(y0)ylab <- deparse(y0)spl <- smooth.spline(y)plot(x, y, type = "n", axes = FALSE, ylab = ylab)usr <- par("usr")xrange <- usr[2] - usr[1]xsize <- xrange * 0.4nsteps <- 8xmar <- xsize / 20yrange <- usr[4] - usr[3]ysize <- yrange / 20ymar <- 0.5 * ysizeX <- seq(usr[1] + xmar, by = xsize / nsteps, length.out = nsteps + 1)xleft <- X[-(nsteps + 1)]xright <- X[-1]pin <- par("pin")maxrad <- xsize / 3 * min(1, pin[2] / pin[1])nrad <- 16minrad <- maxrad / nradRad <- seq(maxrad, by = (minrad - maxrad) / (nrad - 1), length.out=nrad)xmar2 <- xmar + maxradymar2 <- (xmar2 / xrange) * pin[1] / pin[2] * yrangeexpr <- quote({    rect(xleft, usr[4] - 1.5 * ysize, xright, usr[4] - ymar,         col = rainbow(8), border = NA)    symbols(rep(usr[2] - xmar2, nrad), rep(usr[3] + ymar2, nrad),            circles = Rad, inches = FALSE, add = TRUE, fg = NA,            bg = gray.colors(nrad + 1, 1, 0)[-1])    points(y)    lines(spl)})rasterPlot(expr, res = 50)box()axis(1)axis(2)## The same picture with higher resolution but no antialiasingplot(y, type = "n", axes = FALSE, ann = FALSE)## No content in margin, but region = "figure" and bg = "white"## paints margin whiterasterPlot(expr, antialias = "none", interpolate = FALSE,           region = "figure", bg = "white")## Draw box, axes, labelsparnew <- par(new = TRUE)plot(x, y, type = "n", ylab = ylab)par(parnew)## Draw plot(1:5) with adjusted margins and additional axes.  Some parts## are drawn with rasterPlot, others normally.  Resize to see stretching.op <- par(no.readonly = TRUE)par(mar = c(5.1, 4.1, 2.1, 2.1))plot(1:5, type = "n", axes = FALSE, ann = FALSE)expr2 <- quote({    points(c(2, 4), c(2, 4))    axis(2)    axis(3)})rasterPlot(expr2, region = "figure", bg = "white")points(c(1, 3, 5), c(1, 3, 5))box()axis(1)axis(4)title(xlab = "Index", ylab = "1:5")par(op)

Regional Curve Standardization

Description

Detrend multiple ring-width series simultaneously using a regionalcurve.

Usage

rcs(rwl, po, nyrs = NULL, f = 0.5, biweight = TRUE, ratios = TRUE,    rc.out = FALSE, make.plot = TRUE, ..., rc.in = NULL, check = TRUE)

Arguments

Details

This method detrends and standardizes tree-ring series by calculatingan age-related growth curve specific to therwl. Thedetrending is the estimation and removal of the tree’s naturalbiological growth trend. The standardization is done by eitherdividing each series by the growth trend or subtracting the growthtrend from each series to produce units in the dimensionlessring-width index (RWI). The option to produce indices bysubtraction is intended to be used on series that have been subject tovariance stabilization (e.g., usingpowt).

The spline approach uses an n-year spline where the frequency responseis 0.50 at a wavelength of 10 percent of the maximum cambial ageunless specified differently usingnyrs andf in the functioncaps.

This attempts to remove the low frequency variability that is due tobiological or stand effects. See the references below for furtherdetails on detrending in general, and Biondi and Qeadan (2008) for anexplanation ofRCS.

Value

Adata.frame containing the dimensionless and detrendedring-width indices with column names, row names and dimensions ofrwl. Ifrc.out isTRUE then alist will be returned with adata.frame containing thedetrended ring widths as above and avector containing theregional curve.

Note

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

Author(s)

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

References

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

Cook, E. R. and Kairiukstis, L. A., editors (1990)Methods ofDendrochronology: 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,chron,cms,caps

Examples

library(utils)data(gp.rwl)data(gp.po)gp.rwi <- rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,              rc.out = TRUE, make.plot = FALSE)str(gp.rwi)gp.rwi <- rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,              make.plot = TRUE, main = "Regional Curve")

ReadDPL Compact Format Ring Width File

This function reads in aDPL compact format file of ringwidths.

Usage

read.compact(fname)

Arguments

Details

This function should be able to read files written by theDendrochronology Program Library (DPL) in its compactformat.

Value

An object of classc("rwl", "data.frame") with the series incolumns and the years as rows. The seriesIDs are thecolumn names and the years are the row names.

Author(s)

Mikko Korpela

See Also

read.rwl,read.tucson,read.tridas,read.fh,write.compact

Examples

## Not run: data(co021)fname <- write.compact(rwl.df = co021,                       fname = tempfile(fileext=".rwl"),                       append = FALSE, prec = 0.001)foo <- read.compact(fname)str(foo)str(co021)all.equal(foo,co021)unlink(fname)## End(Not run)

Read Tucson Format Chronology File

Description

This function reads in a Tucson (decadal) format file of tree-ringchronologies (.crn).

Usage

read.crn(fname, header = NULL, encoding = getOption("encoding"),         long = TRUE)

Arguments

Details

This reads in a standard crn file as defined according to thestandards of theITRDB athttps://www1.ncdc.noaa.gov/pub/data/paleo/treering/treeinfo.txt. Despite thestandards at theITRDB, this occasionally fails due toformatting problems.

Value

Adata.frame with each chronology in columns and the years asrows. The chronologyIDs are the column names and the yearsare the row names. If the file includes sample depth that is includedas the last column (samp.depth). The output class isclass "crn" and "data.frame"

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

Read Heidelberg Format Ring Width File

read.fh(fname, BC_correction = FALSE)

Rinn, F. (2003)TSAP-Win User Reference Manual.Rinntech, Heidelberg.https://rinntech.info/products/tsap-win/.

See Also

read.rwl

Read Site-Tree-CoreIDs

These functions try to read site, tree, and coreIDs from arwldata.frame.

Usage

read.ids(rwl, stc = c(3, 2, 3), ignore.site.case = FALSE,         ignore.case = FALSE, fix.typos = FALSE, typo.ratio = 5,         use.cor = TRUE)autoread.ids(rwl, ignore.site.case = TRUE, ignore.case = "auto",             fix.typos = TRUE, typo.ratio = 5, use.cor = TRUE)

Arguments

The following parameters affect the handling of suspected typingerrors. Some have different default values inread.ids andautoread.ids.

Details

Because dendrochronologists often take more than one core per tree, itis occasionally useful to calculate within vs. between tree variance.The International Tree Ring Data Bank (ITRDB) allows thefirst eight characters in an rwl file for seriesIDs butthese are often shorter. Typically the creators of rwl files use alogical labeling method that can allow the user to determine the treeand coreID from the label.

Argumentstc tells how each series separate into site,tree, and coreIDs. For instance a series code might be"ABC011" indicating site"ABC", tree 1, core 1. If thisformat is consistent then thestc mask would bec(3, 2, 3) allowing up to three characters for the coreID (i.e., pad to the right). If it is not possible todefine the scheme (and often it is not possible to machine readIDs), then the outputdata.frame can be builtmanually. See Value for format.

The functionautoread.ids is a wrapper toread.ids withstc="auto", i.e. automatic detection of the site / tree / corescheme, and different default values of some parameters. In automaticmode, the names in the samerwl can even follow differentsite / tree / core schemes. As there are numerous possible encodingschemes for naming measurement series, the function cannot alwaysproduce the correct result.

Withstc="auto", the site part can be one of the following.

• In names mostly consisting of numbers, the longest commonprefix is the site part

• Alphanumeric site part ending with alphabet, when followed bynumbers and alphabets

• Alphabetic site part (quite complicated actualdefinition). Settingignore.case to"auto"allows the function to try to guess when a case change in the middleof a sequence of alphabets signifies a boundary between the sitepart and the tree part.

• The characters before the first sequence of space /punctuation characters in a name that contains at least two suchsequences

These descriptions are somewhat general, and the detailscan be found in regular expressions inside the function. If a namedoes not match any of the descriptions, it is matched against apreviously found site part, starting from the longest.

The followingID schemes are detected and supported in thetree / core part. The detection is done per site.

• Numbers in tree part, core part starts with something else

• Alphabets in tree part, core part starts with something else

• Alphabets, either tree part all lower case and core part allupper case or vice versa. For this to work,ignore.case must be set to"auto" orFALSE.

• All digits. In this case, the number of characters belongingto the tree and core parts is detected with one of the followingmethods.

  • If numeric tree parts were found before, it is assumed thatthe core part is missing (one core per tree).

  • It the series are numbered continuously, one core per treeis assumed.

  • Otherwise, try to find a core part as the suffix so that thecores are numbered continuously.

  If none of the above fits, the tree / core split of the all-digitnames will be decided with the methods described further down thelist, or finally with the fallback mechanism.

• The combined tree / core part is empty or one character.In this case, the core part is assumed to be missing.

• Tree and core parts separated by a punctuation or white spacecharacter

If the split of a tree / core part cannot be found with any of themethods described above, the prefix of the string is matched against apreviously found tree part, starting from the longest. The fallbackmechanism for the still undecided tree / core parts is one of thefollowing. The first one is used ifuse.cor isTRUE, number two if it isFALSE.

1. Pairwise correlation coefficients are computed between allremaining series. Pairs of series with above median correlationare flagged as similar, and the other pairs are flagged asdissimilar. Each possible number of characters (minimum 1) isconsidered for the share of the treeID. Thecorresponding unique would-be treeIDs determine a set ofclusterings where one cluster is formed by all the measurementseries of a single tree. For each clustering (allocation ofcharacters), an agreement score is computed. The agreement score isdefined as the sum of the number of similar pairs with matchingcluster number and the number of dissimilar pairs with non-matchingcluster number. The number of characters with the maximum agreementis chosen.

2. If the majority of the names in the site usekcharacters for the tree part, that number is chosen. Otherwise, onecore per tree is assumed. Parametertypo.ratio has adouble meaning as it also defines what is meant by majority here: atleasttypo.ratio / (typo.ratio + 1) *n.tot, wheren.tot is the number of names in the site.

In both fallback mechanisms, the number of characters allocated forthe tree part will be increased until all trees have a non-zeroID or there are no more characters.

Suspected typing errors will be fixed by the function iffix.typos isTRUE. The parametertypo.ratio affects the eagerness to fix typos, i.e. thenumber of counterexamples required to declare a typo. The followingmain typo fixing mechanisms are implemented:

SiteIDs.

If a rare site string resembles an atleasttypo.ratio times more frequent alternative, andif fixing it would not create any name collisions, make the fix.The alternative string must be unique, or if there is more thanone alternative, it is enough if only one of them is a look-alikestring. Any kind of substitution in one character place isallowed if the alternative string has the same length as theoriginal string. The alternative string can be one characterlonger or one character shorter than the original string, but onlyif it involves interpreting one digit as the look-alike alphabetor vice versa. There are requirements to how long a site stringmust be in order to be eligible for replacement / typo fixing,i.e. cannot be shortened to zero length, cannot change the onlycharacter of a site string. The parametersignore.case andignore.site.case havesome effect on this typo fixing mechanism.

Tree and coreIDs.

If all tree / core parts of asite have the same length, each character position is inspectedindividually. If the characters in thei:th position arepredominantly digits (alphabets), any alphabets (digits) arechanged to the corresponding look-alike digit (alphabet) if thereis one. The look-alike groups are {0, O, o}, {1, I, i}, {5,S, s} and {6, G}. The parametertypo.ratiodetermines the decision threshold of interpreting the type of eachcharacter position as alphabet (digit): the ratio of alphabets(digits) to the total number of characters must be at leasttypo.ratio / (typo.ratio + 1). If a namediffers from the majority type in more than one characterposition, it is not fixed. Also, no fixes are performed if any ofthem would cause a possible monotonic order of numeric prefixes tobreak.

The function attempts to convert the tree and core substrings tointegral values. When this succeeds, the converted values are copiedto the output without modification. When non-integral substrings areobserved, each unique tree is assigned a unique integral value. Thesame applies to cores within a tree, but there are some subtletieswith respect to the handling of duplicates. Substrings are sortedbefore assigning thenumericIDs.

The order of columns inrwl, in most cases, does notaffect the tree and coreIDs assigned to each series.

Value

Adata.frame with column one named"tree" giving anID for each tree and column two named"core" givinganID for each core. The original seriesIDs arecopied from rwl as rownames. The order of the rows in the outputmatches the order of the series inrwl. If more than onesite is detected, an additional third column named"site" willcontain a siteID. All columns have integral valuednumeric values.

Author(s)

Andy Bunn (original version) and Mikko Korpela (patches,stc="auto",fix.typos, etc.).

See Also

rwi.stats,read.rwl

Examples

library(utils)data(ca533)read.ids(ca533, stc = c(3, 2, 3))autoread.ids(ca533)

Read Ring Width File

Description

This function reads in a file of ring widths (.rwl) in one of theavailable formats.

Usage

read.rwl(fname,         format = c("auto", "tucson", "compact", "tridas", "heidelberg", "csv"),         ...)

Arguments

Details

This is a simple wrapper to the functions actually implementing theread operation.

Value

If a"tucson","compact","heidelberg","csv" file isread (even through"auto"), returns an object of classc("rwl", "data.frame") with the series in columns and the yearsas rows. The seriesIDs are the column names and the yearsare the row names.

If a"tridas" file is read (even through"auto"),returns a list of results. Seeread.tridas for moreinformation.

Author(s)

Mikko Korpela

See Also

read.tucson,read.compact,read.tridas,read.fh,csv2rwl,write.rwl

Read Tree Ring Data Standard (TRiDaS) File

Description

This function reads in a TRiDaS formatXML file.Measurements, derived series and various kinds of metadata aresupported.

Usage

read.tridas(fname, ids.from.titles = FALSE,            ids.from.identifiers = TRUE, combine.series = TRUE,            trim.whitespace = TRUE, warn.units = TRUE)

Arguments

Details

The Tree Ring Data Standard (TRiDaS) is described in Jansma et. al(2010).

The parameters used for rearranging (ids.from.titles,ids.from.identifiers) and combining(combine.series) measurement series only affect the fourlowest levels of document structure: element, sample, radius,measurementSeries. Series are not reorganized or combined at theupper structural levels (project, object).

Value

A list with a variable number of components according to the contentsof the input file. The possible list components are:

Note

This is an early version of the function. Bugs are likely toexist, and parameters and return values are subject to change. Notall metadata defined in the TRiDaS specification is supported –unsupported elements are quietly ignored.

Author(s)

Mikko Korpela

References

Jansma, E., Brewer, P. W., and Zandhuis, I. (2010) TRiDaS 1.1: Thetree-ring data standard.Dendrochronologia,28(2),99–130.

See Also

read.rwl,read.tucson,read.compact,read.fh,write.tridas

Read Tucson Format Ring Width File

Description

This function reads in a Tucson (decadal) format file of ring widths (.rwl).

Usage

read.tucson(fname, header = NULL, long = FALSE,            encoding = getOption("encoding"),             edge.zeros = TRUE, verbose = TRUE)

Arguments

Details

This reads in a standard rwl file as defined according to thestandards of theITRDB athttps://www1.ncdc.noaa.gov/pub/data/paleo/treering/treeinfo.txt. Despite thestandards at theITRDB, this occasionally fails due toformatting problems.

Value

An object of classc("rwl", "data.frame") with the series incolumns and the years as rows. The seriesIDs are thecolumn names and the years are the row names.

Author(s)

Andy Bunn. Patched and greatly improved by Mikko Korpela.

See Also

read.rwl,read.compact,read.tridas,read.fh,write.tucson

Red-Noise Spectra of Time-Series

Description

Estimate red-noise spectra from a possibly unevenly spaced time-series.

Usage

redfit(x, t, tType = c("time", "age"), nsim = 1000, mctest = TRUE,       ofac = 4, hifac = 1, n50 = 3, rhopre = NULL,       p = c(0.10, 0.05, 0.02), iwin = 2,       txOrdered = FALSE, verbose = FALSE, seed = NULL,       maxTime = 10, nLimit = 10000)runcrit(n, p = c(0.10, 0.05, 0.02), maxTime = 10, nLimit = 10000)

Arguments

Details

Functionredfit computes the spectrum of a possibly unevenlysampled time-series by using the Lomb-Scargle Fourier transform. Thespectrum is bias-corrected using spectra computed from simulated AR1series and the theoretical AR1 spectrum.

The function duplicates the functionality of program REDFIT by Schulzand Mudelsee. See the manual of that program for more information.The results of this function should be very close to REDFIT. However,some changes have been made:

• More precision is used in some constants and computations.

• All the data are used: the last segment always contains thelast pair of (t,x). There may be small differencesbetweenredfit and REDFIT with respect to the number ofpoints per segment and the overlap of consecutive segments.

• The critical values of the runs test (see the description ofruncrit below) differ betweenredfit and REDFIT. Theapproximate equations in REDFIT produce values quite far off fromthe exact values when the number of frequencies is large.

• The user can select the significance levels of the runs test.

• Most of the window functions have been adjusted.

• 6 dB bandwidths have been computed for discrete-time windows.

Functionruncrit computes the limits of the acceptance regionof a number of runs test: assuming a sequence ofn i.i.d.discrete random variables with two possible valuesa andbof equal probability (0.5), we are examining the distribution of thenumber of runs. A run is an uninterrupted sequence of onlya oronlyb. The minimum number of runs is 1 (a sequence with onlya or onlyb) while the maximum number isn(alternatinga andb). See Bradley, p. 253–254,259–263. The function is also called fromredfit;seercnt in ‘Value’ for the interpretation. Inthis case the argumentsp,maxTime andnLimit are passed fromredfit toruncrit,andn is the number of output frequencies.

The results ofruncrit have been essentially precomputed forsome values ofp andn. If a precomputedresult is not found andn is not too large(nLimit,maxTime), the exact results arecomputed on-demand. Otherwise, or if package"gmp" is notinstalled, the normal distribution is used for approximation.

Value

Functionruncrit returns alist containingrcritlo,rcrithi andrcritexact(see below). Functionredfit returns alist with thefollowing elements:

Author(s)

Mikko Korpela. Examples by Andy Bunn.

References

Functionredfit is based on the Fortran programREDFIT(version 3.8e), which is in the public domain.

Bradley, J. V. (1968)Distribution-Free StatisticalTests. Prentice-Hall.

Schulz, M. and Mudelsee, M. (2002) REDFIT: estimating red-noisespectra directly from unevenly spaced paleoclimatic time series.Computers & Geosciences,28(3), 421–426.

See Also

print.redfit

Examples

# Create a simulated tree-ring width series that has a red-noise# background ar1=phi and sd=sigma and an embedded signal with # a period of 10 and an amplitude of have the rednoise sd.library(graphics)library(stats)runif(1)rs <- .Random.seedset.seed(123)nyrs <- 500yrs <- 1:nyrs# Here is an ar1 time series with a mean of 2mm,# an ar1 of phi, and sd of sigmaphi <- 0.7sigma <- 0.3sigma0 <- sqrt((1 - phi^2) * sigma^2)x <- arima.sim(list(ar = phi), n = nyrs, sd = sigma0) + 2# Here is a sine wave at f=0.1 to add in with an amplitude# equal to half the sd of the red noise backgroundper <- 10amp <- sigma0 / 2wav <- amp * sin(2 * pi / per * yrs)# Add them together so we have signal and noisex <- x + wav# Here is the redfit specredf.x <- redfit(x, nsim = 500)# Acceptance region of number of runs test# (not useful with default arguments of redfit())runcrit(length(redf.x[["freq"]]))op <- par(no.readonly = TRUE) # Save to reset on exitpar(tcl = 0.5, mar = rep(2.2, 4), mgp = c(1.1, 0.1, 0))plot(redf.x[["freq"]], redf.x[["gxxc"]],     ylim = range(redf.x[["ci99"]], redf.x[["gxxc"]]),     type = "n", ylab = "Spectrum", xlab = "Frequency (1/yr)",     axes = FALSE)grid()lines(redf.x[["freq"]], redf.x[["gxxc"]], col = "#1B9E77")lines(redf.x[["freq"]], redf.x[["ci99"]], col = "#D95F02")lines(redf.x[["freq"]], redf.x[["ci95"]], col = "#7570B3")lines(redf.x[["freq"]], redf.x[["ci90"]], col = "#E7298A")freqs <- pretty(redf.x[["freq"]])pers <- round(1 / freqs, 2)axis(1, at = freqs, labels = TRUE)axis(3, at = freqs, labels = pers)mtext(text = "Period (yr)", side = 3, line = 1.1)axis(2); axis(4)legend("topright", c("x", "CI99", "CI95", "CI90"), lwd = 2,       col = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A"),       bg = "white")box()## Not run: # Second example with tree-ring data# Note the long-term low-freq signal in the data. E.g.,# crn.plot(cana157)library(utils)data(cana157)yrs <- time(cana157)x <- cana157[, 1]redf.x <- redfit(x, nsim = 1000)plot(redf.x[["freq"]], redf.x[["gxxc"]],     ylim = range(redf.x[["ci99"]], redf.x[["gxxc"]]),     type = "n", ylab = "Spectrum", xlab = "Frequency (1/yr)",     axes = FALSE)grid()lines(redf.x[["freq"]], redf.x[["gxxc"]], col = "#1B9E77")lines(redf.x[["freq"]], redf.x[["ci99"]], col = "#D95F02")lines(redf.x[["freq"]], redf.x[["ci95"]], col = "#7570B3")lines(redf.x[["freq"]], redf.x[["ci90"]], col = "#E7298A")freqs <- pretty(redf.x[["freq"]])pers <- round(1 / freqs, 2)axis(1, at = freqs, labels = TRUE)axis(3, at = freqs, labels = pers)mtext(text = "Period (yr)", side = 3, line = 1.1)axis(2); axis(4)legend("topright", c("x", "CI99", "CI95", "CI90"), lwd = 2,       col = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A"),       bg = "white")box()par(op)## End(Not run).Random.seed <- rs

(Running Window) Statistics on Detrended Ring-Width Series

Description

These functions calculate descriptive statistics on adata.frame of (usually) ring-width indices. The statistics areoptionally computed in a running window with adjustable length andoverlap. The data can be filtered so that the comparisons are made toon just high-frequency data.

Usage

rwi.stats.running(rwi, ids = NULL, period = c("max", "common"),                  method = c("spearman", "pearson","kendall"),                  prewhiten=FALSE,n=NULL,                  running.window = TRUE,                  window.length = min(50, nrow(rwi)),                  window.overlap = floor(window.length / 2),                  first.start = NULL,                  min.corr.overlap = min(30, window.length),                  round.decimals = 3,                  zero.is.missing = TRUE)rwi.stats(rwi, ids=NULL, period=c("max", "common"),           method = c("spearman", "pearson","kendall"), ...)rwi.stats.legacy(rwi, ids=NULL, period=c("max", "common"))

Arguments

Details

This calculates a variety of descriptive statistics commonly used indendrochronology.

The functionrwi.stats is a wrapper that callsrwi.stats.running withrunning.window = FALSE.The results may differ from those prior to dplR 1.5.3, where theformerrwi.stats (now renamed torwi.stats.legacy) wasreplaced with a call torwi.stats.running.

For correctly calculating the statistics on within and between seriesvariability, an appropriate mask (parameterids) must beprovided that identifies each series with a tree as it is common fordendrochronologists to take more than one core per tree. The functionread.ids is helpful for creating a mask based on theseriesID.

Ifids has duplicate tree/core combinations, thecorresponding series are averaged before any statistics are computed.The value of the parameterzero.is.missing is relevant inthe averaging:TRUE ensures that zeros don’tcontribute to the average. The default value ofzero.is.missing isTRUE. The default prior todplR 1.5.3 wasFALSE. If the parameter is set toFALSE,the user will be warned in case zeros are present. Duplicatetree/core combinations are not detected byrwi.stats.legacy.

Row names ofids may be used for matching theIDs with series inrwi. In this case, thenumber of rows inids is allowed to exceed the number ofseries. If some names ofrwi are missing from the rownames ofids, the rows ofids are assumed tobe in the same order as the columns ofrwi, and thedimensions must match. The latter is also the way thatrwi.stats.legacy handlesids, i.e. names areignored and dimensions must match.

Note thatperiod = "common" can produceNaN for many ofthe stats if there is no common overlap period among the cores. Thishappens especially in chronologies with floating subfossil samples(e.g.,ca533).

Some of the statistics are specific to dendrochronology (e.g., theeffective number of cores or the expressed population signal). Usersunfamiliar with these should see Cook and Kairiukstis (1990) andFritts (2001) for further details for computational details on theoutput. The signal-to-noise ratio is calculated following Cook and Pederson (2011).

Note that Buras (2017) cautions against using the expressed population signal as a statistic to determine the whether a chronology correctly represents the population signal of a data set. He reccomends the use of subsample signal strength (sss) over EPS.

If desired, therwi can be filtered in the same manneras the family of cross-dating functions usingprewhiten andn. See the help page forcorr.rwl.seg formore details.

Value

Adata.frame containing the following columns (each rowcorresponds to one position of the window):

Note

This function uses theforeach loopingconstruct with the%dopar% operator.For parallel computing and a potential speedup, a parallel backendmust be registered before running the function.

Author(s)

Mikko Korpela, based onrwi.stats.legacy by AndyBunn

References

Buras, A. (2017) A comment on the Expressed Population Signal. Dendrochronologia 44:130-132.

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

Cook, E. R. and Pederson, N. (2011) Uncertainty, Emergence, andStatistics in Dendrochronology. In Hughes, M. K., Swetnam, T. W., andDiaz, H. F., editors,Dendroclimatology: Progress andProspects, pages 77–112. Springer.ISBN-13:978-1-4020-4010-8.

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

See Also

detrend,cor,read.ids,rwi.stats,corr.rwl.seg

Examples

library(utils)data(gp.rwl)data(gp.po)gp.rwi <- cms(rwl = gp.rwl, po = gp.po)gp.ids <- read.ids(gp.rwl, stc = c(0, 2, 1))# On a running windowrwi.stats.running(gp.rwi, gp.ids)## With no running window (i.e. running.window = FALSE)rwi.stats(gp.rwi, gp.ids)## Restrict to common overlap (in this case 1899 to 1987)rwi.stats(gp.rwi, gp.ids, period="common")rwi.stats.legacy(gp.rwi, gp.ids) # rwi.stats prior to dplR 1.5.3

Do some reporting on a RWL object

Description

This function generates a small report on arwl (orrwi) objectthat gives the user some basic information on the data including the number ofseries, the span of the data, the mean interseries correlation, the number ofmissing rings (zeros), internalNA values, and rings that are very small,or very large.

Usage

rwl.report(rwl,small.thresh=NA,big.thresh=NA)

Arguments

Details

This generates information about arwl object including the number of series, the mean length of all the series, the first year, last year, the mean first-order autocorrelation (viasummary.rwl), the mean interseries correlation (viainterseries.cor), the years where a series has a missing ring (zero), internal NA, very small ring, very large rings, etc.

This output of this function is not typically meant for the user to access but has aprint method for the user.

Value

Alist with elements containing descriptive information on therwl object. Specifically:

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

read.rwl,summary.rwl,interseries.cor

Examples

data("gp.rwl")rwl.report(rwl = gp.rwl)# list very small (smallest 1pct) of rings as wellone.pct <- quantile(gp.rwl[gp.rwl != 0], na.rm=TRUE, probs=0.01)rwl.report(rwl = gp.rwl, small.thresh = one.pct)

Calculate Descriptive Summary Statistics on Ring-Width Series

Description

This function calculates descriptive statistics on arwl objectof raw or detrended ring-width series.

Usage

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

Arguments

Details

This calculates a variety of descriptive statistics commonly used indendrochronology (see below). Users unfamiliar with these should seeCook and Kairiukstis (1990) and Fritts (2001) for further details.

Thesummary method for class"rwl" is a wrapperforrwl.stats.

Value

Adata.frame containing descriptive stats on each"series". These are the first and last year of the series aswell as the length of the series ("first","last","year"). The mean, median, standard deviation are given("mean","median","stdev") as are the skewness,the excess kurtosis (calculated as Pearson’s kurtosis minus 3), the Gini coefficient, and first orderautocorrelation ("skew","kurtosis","gini.coef","ar1").

Note that prior to version 1.6.8, two measures of sensitivity were also included. However mean sensitivity is not a robust statistic that should rarely, if ever, be used (Bunn et al. 2013). Those sensitivity functions ("sens1" and"sens2") are still available for continuity. Users should consider the coef of variation in lieu of mean sensitivity.

Author(s)

Andy Bunn. Slightly improved by Mikko Korpela.

References

Bunn, A. G., Jansma, E., Korpela, M., Westfall, R. D., and Baldwin,J. (2013) Using simulations and data to evaluate mean sensitivity(\zeta) as a useful statistic in dendrochronology.Dendrochronologia,31(3), 250–254.

Cook, E. R. and Kairiukstis, L. A., editors (1990)Methods ofDendrochronology: 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

rwi.stats,read.rwl

Examples

library(utils)data(ca533)rwl.stats(ca533)summary(ca533)

Superposed Epoch Analysis

Description

This function calculates the significance of the departure from themean for a given set of key event years and lagged years.

Usage

sea(x, key, lag = 5, resample = 1000)

Arguments

Details

Superposed epoch analysis (SEA) is used to test the significance of a meantree growth response to certain events (such as droughts). Departuresfrom the meanRWI values for the specified years prior toeach event year, the event year, and the specified years immediatelyafter each event are averaged to a superposed epoch. To determine ifRWI for these years was significantly different fromrandomly selected sets oflag+1 other years, bootstrapresampling is used to randomly select sets oflag+1 yearsfrom the data set and to estimate significances for the departuresfrom the meanRWI.

SEA computation is based on scaledRWI values, and95%-confidence intervals are computed for the scaled values for eachyear in the superposed epoch.

Value

Adata.frame with

Author(s)

Christian Zang. Patched and improved by Mikko Korpela.

References

Lough, J. M. and Fritts, H. C. (1987) An assessment of the possibleeffects of volcanic eruptions on North American climate usingtree-ring data, 1602 to 1900AD.Climatic Change,10(3), 219–239.

Examples

library(graphics)library(utils)data(cana157)event.years <- c(1631, 1742, 1845)cana157.sea <- sea(cana157, event.years)foo <- cana157.sea$se.unscalednames(foo) <- cana157.sea$lagbarplot(foo, col = ifelse(cana157.sea$p < 0.05, "grey30", "grey75"),         ylab = "RWI", xlab = "Superposed Epoch")

Segment Plot

seg.plot(rwl, ...)

spag.plot

Examples

library(utils)data(co021)seg.plot(co021)

Calculate Mean Sensitivity

Description

This function calculates mean sensitivity of a detrended ring-widthseries.

Usage

  sens1(x)

Arguments

Details

This calculates mean sensitivity according to Eq. 1 in Biondi andQeadan (2008). This is the standard measure of sensitivity indendrochronology and is typically calculated on detrended series.However, note that mean sensitivity is not a robust statistic and should rarely, if ever, be used (Bunn et al. 2013).

Value

the mean sensitivity.

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.

Bunn, A. G., Jansma, E., Korpela, M., Westfall, R. D., and Baldwin,J. (2013) Using simulations and data to evaluate mean sensitivity(\zeta) as a useful statistic in dendrochronology.Dendrochronologia,31(3), 250–254.

See Also

sens2,rwl.stats

Examples

library(utils)data(ca533)ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")sens1(ca533.rwi[, 1])

Calculate Mean Sensitivity on Series with a Trend

Description

This function calculates mean sensitivity of a raw or detrendedring-width series.

Usage

  sens2(x)

Arguments

Details

This calculates mean sensitivity according to Eq. 2 in Biondi andQeadan (2008). This is a measure of sensitivity in dendrochronologythat is typically used in the presence of a trend. However, note that mean sensitivity is not a robust statistic and should rarely, if ever, be used (Bunn et al. 2013).

Value

the mean sensitivity.

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.

Bunn, A. G., Jansma, E., Korpela, M., Westfall, R. D., and Baldwin,J. (2013) Using simulations and data to evaluate mean sensitivity(\zeta) as a useful statistic in dendrochronology.Dendrochronologia,31(3), 250–254.

See Also

sens1,rwl.stats

Examples

library(utils)data(ca533)ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")sens2(ca533.rwi[, 1])

Plot Series and a Master

Description

Plots a tree-ring series with a master chronology and displaystheir fit, segments, and detrending options in support of thecross-dating functions.

Usage

series.rwl.plot(rwl, series, series.yrs = as.numeric(names(series)),                seg.length = 100, bin.floor = 100, n = NULL,                prewhiten = TRUE, biweight = TRUE, floor.plus1 = FALSE)

Arguments

Details

The function is typically invoked to produce four plots showing theeffect of the detrending optionsn andprewhiten and the binning optionsseg.lengthandbin.floor.

Plot 1

Time series plot of the filtered series and themaster

Plot 2

Scatterplot of series vs. master

Plot 3

Segments that would be used in the other cross-datingfunctions (e.g.,corr.series.seg)

Plot 4

Text giving the detrending options and the time spanof the raw and filtered series and master

The series and master are returned as well.

See help pages forcorr.rwl.seg,corr.series.seg, andccf.series.rwl formore information on these arguments.

Value

Alist containing the filtered vectorsseries andmaster.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

corr.rwl.seg,corr.series.seg,ccf.series.rwl

Examples

library(utils)data(co021)foo <- series.rwl.plot(rwl = co021, series = "646244", seg.length = 100,                       n = 5)## note effect of n on first year in the seriesfoo <- series.rwl.plot(rwl = co021, series = "646244", seg.length = 100,                       n = 13, prewhiten = FALSE)bar <- series.rwl.plot(rwl = co021, series = "646244", seg.length = 100,                       n = 7, prewhiten = FALSE)head(foo$series)head(bar$series)

Synchronous Growth Changes

Description

This function calculates the synchronous growth changes (sgc), semi synchronous growth changes (ssgc) and the length of the compared overlap for a given set of tree-ring records. Optionally the probability of exceedence is calculated.

Usage

  sgc(x,overlap = 50, prob = TRUE)

Arguments

Details

The sgc is a non parametric test based on sign tests.The synchronous growth changes (sgc) and semi synchronous growth changes (ssgc) are meant to replace the Gleichläufigkeit (glk()), since the Gleichläufigkeit can be (strongly) influenced by years when one of the compared series shows no growth change. The sgc gives a better description of the similarity (Visser, 2020). The ssgc gives the percentage years that one of the compared series shows no growth change. This function implements sgc and ssgc as the vectorized pairwise comparison of all records in data set.

The probability of exceedence (p) for the sgc expresses the chance that the sgc is incorrect. The observed value of the sgc is converted to a z-score and based on the standard normal curve the probability of exceedence is calculated (Visser 2020). The result is a matrix of all p-values.

Value

Alist with three or four matrices (p_mat is optional ifprob = TRUE):