dplR Basics
Using R for Tree-Ring Analysis
Andy Bunn, Mikko Korpela
25-June-2018
1 Abstract
This document describes basic features of dplR by following the initial steps that an analyst might follow when working with a new tree-ring data set. The vignette starts with reading in ring widths and plotting them. We describe a few of the available methods for detrending and then show how to extract basic descriptive statistics. We show how to build and plot a simple mean-value chronology. We also show how to build a chronology using the expressed population signal from the detrended ring widths as an example of how more complicated analysis can be done using dplR.
2 What is Covered
The Dendrochronology Program Library in R (dplR) is a package for dendrochronologists to handle data processing and analysis. This document gives just a brief introduction of some of the most commonly used functions in dplR. There is more detailed information available in the help files and in the literature (Bunn, 2008).
In this vignette, we will walk through the most basic activities of working with tree-ring data in roughly the order that a user might follow. E.g., reading data, detrending, chronology building, and doing preliminary exploratory data analysis via descriptive statistics.
2.1 Load dplR
We will be using dplR in here. Load it:
library(dplR)3 Working with Ring-Width Data
3.1 Reading Data
There are many ways that tree-ring data are digitally stored. These range in sophistication from the simple (and commonly used) Tucson/decadal format file of ring widths to the more complex (but richer) TRiDaS format. The type of data used most often by dendrochronologists is a series of ring widths from tree cores. We generally refer to these as rwl objects for “ring width length” but there is no reason these cannot be other types of tree-ring data (e.g., density).
The workhorse function for getting tree-ring data into R is dplR’s read.rwl function. This function reads files in "tucson", "compact", "tridas", and "heidelberg" formats. The on-board rwl data sets in dplR (i.e., co021, ca533, gp.rwl) were all imported into R using this function.
Throughout this vignette we will use the on-board data set ca533 which gives the raw ring widths for bristlecone pine Pinus longaeva at Campito Mountain in California, USA. There are 34 series spanning 1358 years.
These objects are structured very simply as a data.frame with the series in columns and the years as rows. The series IDs are the column names and the years are the row names (both stored as characters). For instance, using the Campito Mountain ring widths we can load the data and learn some basic things about it:
data(ca533) # the result of ca533 <- read.rwl('ca533.rwl')
nrow(ca533) # 1358 years## [1] 1358ncol(ca533) # 34 series## [1] 34We can look a little deeper at this object (ca533) and get the series names as well as look at the time values associates with the data:
colnames(ca533) # the series IDs## [1] "CAM011" "CAM021" "CAM031" "CAM032" "CAM041" "CAM042" "CAM051"
## [8] "CAM061" "CAM062" "CAM071" "CAM072" "CAM081" "CAM082" "CAM091"
## [15] "CAM092" "CAM101" "CAM102" "CAM111" "CAM112" "CAM121" "CAM122"
## [22] "CAM131" "CAM132" "CAM141" "CAM151" "CAM152" "CAM161" "CAM162"
## [29] "CAM171" "CAM172" "CAM181" "CAM191" "CAM201" "CAM211"head(time(ca533),n = 10) # the first 10 years## [1] 626 627 628 629 630 631 632 633 634 6353.2 Describing and Plotting Ring-Width Data
Once a rwl data set has been read into R, there are a variety of ways to describe and visualize those data. Take note that this object is stored both as a generic data.frame in R but it also is part of a special class called rwl which will let R know how to do some special things with it like summarize and plot the data:
class(ca533) ## [1] "rwl" "data.frame"Thus, we can plot a rwl object by showing either the segments arranged over time as straight lines or as a “spaghetti plot”. The rwl objects have generic S3 methods for plot and summary meaning that R knows how to do something special when plot or summary are invoked on an object with class rwl. E.g.,:
plot(ca533, plot.type="spag")
4 Descriptive Statistics
The simplest report on a rwl object can be print to the screen via:
rwl.report(ca533)## Number of dated series: 34
## Number of measurements: 23276
## Avg series length: 684.6
## Range: 1358
## Span: 626 - 1983
## Mean (Std dev) series intercorrelation: 0.6294 (0.08593)
## Mean (Std dev) AR1: 0.7093 (0.09812)
## -------------
## Years with absent rings listed by series
## Series CAM011 -- 1753 1782
## Series CAM031 -- 1497 1500 1523 1533 1540 1542 1545 1578 1579 1580 1655 1668 1670 1681
## Series CAM032 -- 1497 1523 1579 1654 1670 1681 1782
## Series CAM051 -- 1475
## Series CAM061 -- 1497 1523 1542 1545 1547 1579 1654 1655 1668 1670 1672 1782 1858 1960
## Series CAM062 -- 1542 1545 1547 1548 1579 1654 1655 1670 1672 1782 1836 1857 1858 1929
## Series CAM071 -- 1269 1497 1498 1523 1542 1547 1578 1579 1612 1655 1656 1668 1670 1672 1674 1690 1707 1708 1756 1782 1795 1820 1836 1845 1857 1858 1924 1948 1960
## Series CAM072 -- 1218 1497 1498 1523 1533 1538 1542 1545 1546 1547 1571 1579 1580 1590 1654 1655 1668 1670 1672 1675 1690
## Series CAM081 -- 1218 1336
## Series CAM082 -- 1362 1858 1865
## Series CAM091 -- 1655 1669 1670 1782 1858
## Series CAM092 -- 1624 1654 1655 1670 1672 1675 1677 1690 1703 1705 1707 1708 1710 1733 1753 1756 1757 1774 1777 1781 1782 1783 1784 1795 1807 1824 1829 1836 1845 1857 1858 1899 1904 1929 1936 1961
## Series CAM101 -- 1782 1783 1899 1929
## Series CAM102 -- 1669 1690 1782 1858 1899 1929
## Series CAM111 -- 1542
## Series CAM112 -- 1542
## Series CAM121 -- 1093 1218 1254 1361 1365 1460 1462 1468 1473 1475 1492 1497 1542 1544 1545 1547 1600 1899 1960
## Series CAM122 -- 1117 1133 1147 1177 1218 1254 1361 1475 1497 1670
## Series CAM131 -- 1361
## Series CAM151 -- 1670 1703
## Series CAM161 -- 1523
## Series CAM162 -- 1618 1624 1641
## Series CAM181 -- 1450 1523
## Series CAM191 -- 1475 1497 1523 1533 1542 1558 1571 1578 1618 1655 1668 1670 1675 1677 1690 1705 1777 1929
## Series CAM201 -- 1523
## Series CAM211 -- 645 762 809 847 924 957 1014 1118 1123 1133 1147 1189 1350 1384 1468 1571 1641
## 234 absent rings (1.005%)
## -------------
## Years with internal NA values listed by series
## NoneThat’s pretty basic information. We can look at some common (and not-so common) descriptive statistics of a rwl object:
ca533.stats <- summary(ca533) # same as calling rwl.stats(ca533)
head(ca533.stats,n=5) # look at the first five series## series first last year mean median stdev skew gini ar1
## 1 CAM011 1530 1983 454 0.440 0.40 0.222 1.029 0.273 0.696
## 2 CAM021 1433 1983 551 0.424 0.40 0.185 0.946 0.237 0.701
## 3 CAM031 1356 1983 628 0.349 0.29 0.214 0.690 0.341 0.808
## 4 CAM032 1435 1983 549 0.293 0.26 0.163 0.717 0.309 0.661
## 5 CAM041 1683 1983 301 0.526 0.53 0.223 0.488 0.238 0.690These are common summary statistics like mean, median, etc. but also statistics that are more specific to dendrochronology like the first-order autocorrelation (ar1), gini (gini), and mean sensitivity (sens1 and sens2). We would be remiss if we did not here mention that mean sensitivity is actually a terrible statistic that should rarely, if ever, be used (Bunn et al., 2013). Note that output object ca533.stats is itself a data.frame and its data can be used to plot, etc.
For instance, we can look at the spread of the first-order autocorrelation via summary(ca533.stats$ar1) or make a plot to show the data. Here we will demonstrate a somewhat involved plot to get you an idea of how to layer plotting commands:
boxplot(ca533.stats$ar1,ylab=expression(phi[1]),col = "lightblue")
stripchart(ca533.stats$ar1, vertical = TRUE,
method = "jitter", jitter = 0.02,add = TRUE, pch = 20, col = 'darkblue',cex=1.25)
ar1Quant <- quantile(ca533.stats$ar1,probs = c(0.25,0.5,0.75))
abline(h=ar1Quant,lty="dashed",col="grey")
mtext(text = names(ar1Quant),side = 4,at = ar1Quant,las=2)
5 Detrending
Analysts often (but not always) detrend a rwl data set to create an object containing ring-width index (rwi) values. The dplR package contains most standard detrending methods including detrending via splines, fitting negative exponential curves, and so on. There are also dplR functions for less commonly used detrending methods like regional curve standardization. There has been much hue and cry for having the “signal-free” method implemented in dplR and that should be done soon (he said optimistically).
A rwi object has the same basic properties as the rwl object from which it is made. I.e., it has the same number of rows and columns, the same names, and so on. The difference is that each series has been standardized by dividing the ring widths against a growth model (e.g., a stiff spline, a negative exponential, etc.). This gives each series a mean of one (thus referred to as “indexed”) and allows a chronology to be built (next section). As read.rwl is the primary function for getting data into R, detrend is the primary function for standardizing rwl objects (but see cms, rcs, bai.in, and bai.out as well).
5.1 Common Detrending Methods
As any dendrochronologist will tell you, detrending is a dark art. In dplR we have implemented some of the standard tools for detrending but all have drawbacks. In all of the methods, the detrending is the estimation and removal of the low frequency variability that is due to biological or stand effects. The standardization is done by dividing each series by the growth trend to produce units in the dimensionless ring-width index (RWI). Much of the text that follows is modified from the help page of detrend.
Probably the most common method for detrending is what is often called the “conservative” approach of attempting to fit a negative exponential curve to a series. In the dplR implementation the "ModNegExp" method of detrend attempts to fit a classic nonlinear model of biological growth of the form \((f(t) = a \times \mathrm{e}^{bt} + k)\), where the argument of the function is time, using nls. See Fritts (2001) for details about the parameters. If a suitable nonlinear model cannot be fit (function is non-decreasing or some values are not positive) then a linear model is fit using lm. That linear model can have a positive slope unless pos.slope is FALSE in which case the series is standardized by its mean (method "Mean" in detrend).
For instance, every series in the ca533 object can be detrended at once via:
ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")This saves the results in ca533.rwi which is a data.frame with the same dimensions as the rwl object ca533 and each series standardized.
nrow(ca533.rwi) # 1358 years## [1] 1358ncol(ca533.rwi) # 34 series## [1] 34colMeans(ca533.rwi, na.rm=TRUE)## CAM011 CAM021 CAM031 CAM032 CAM041 CAM042 CAM051 CAM061 CAM062 CAM071
## 0.9996 1.0000 1.0000 1.0000 1.0000 1.0012 1.0002 0.9999 1.0000 1.0000
## CAM072 CAM081 CAM082 CAM091 CAM092 CAM101 CAM102 CAM111 CAM112 CAM121
## 1.0000 1.0000 1.0000 1.0000 0.9996 1.0000 1.0000 1.0000 1.0000 1.0000
## CAM122 CAM131 CAM132 CAM141 CAM151 CAM152 CAM161 CAM162 CAM171 CAM172
## 1.0000 0.9998 0.9985 0.9999 0.9995 0.9999 1.0004 0.9994 0.9997 0.9998
## CAM181 CAM191 CAM201 CAM211
## 1.0000 0.9953 1.0000 0.9998When detrend is run on a rwl object the function loops through each series. It does this by calling a different function (detrend.series) for each column in the rwl object. But, a user can also call detrend.series and it is useful to do so here for educational purposes.
Let us detrend a single series and apply more than one detrending method when we call the detrend function. We will call detrend.series using the verbose mode so that we can see the parameters applied for each method. The detrend.series function produces a plot by default.
CAM011.rwi <- detrend.series(y = ca533[, "CAM011"],verbose=TRUE)## Verbose output:
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Options
## make.plot TRUE
## method(s)1 c("Spline", "ModNegExp", "Mean", "Ar", "Friedman", "ModHugershoff"
## method(s)2 )
## nyrs NULL
## f 0.5
## pos.slope FALSE
## constrain.nls never
## verbose TRUE
## return.info FALSE
## wt default
## span cv
## bass 0
## difference FALSE
##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Zero indices in input series:
## 1128 1157
##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Detrend by ModNegExp.
## Trying to fit nls model...
## nls coefs
## a: 0.66110
## b: -0.01184
## k: 0.31793
##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Detrend by ModHugershoff.
## Trying to fit nls model...
## nls coefs
## a: 0.45550
## b: 0.15420
## g: 0.01532
## d: 0.32392
##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Detrend by spline.
## Spline parameters
## nyrs = 304, f = 0.5
##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Detrend by mean.
## Mean = 0.4396
##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Detrend by prewhitening.
## Call:
## ar(x = y[idx.goody])
##
## Coefficients:
## 1 2 3 4 5 6 7 8 9
## 0.388 0.139 0.000 0.084 0.132 0.061 0.038 -0.126 0.037
## 10 11 12 13 14 15 16 17 18
## -0.100 -0.010 0.015 0.088 0.010 0.064 -0.013 0.015 -0.004
## 19 20 21 22 23
## -0.054 0.124 -0.030 -0.054 0.137
##
## Order selected 23 sigma^2 estimated as 0.0209
##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Detrend by FriedMan's super smoother.
## Smoother parameters
## span = cv, bass = 0
## wt = default
##
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Zero indices in Ar series:
## 993
Note that advanced users can use return.info=TRUE to have all the curve-fitting parameters returned in a list. Having access to these curves is occassionally desireable. Also. sometimes, a user will want to interactively detrend each series and fit a negative exponential curve to one series, a spline to another, and the mean to a third. This can be done via the i.detrend and i.detrend.series functions. See their help pages for details. Advanced users can also power transform tree-ring data (powt) and detrend via substraction instead of division using detrend and detrend.series.
5.2 Other Detrending Methods
There are other detrending methods that are less commonly used but might have theoretical advantages. These include regional curve standardization (function rcs), C-Method Standardization (function cms), and converting measurements of ring widths to basal area increment (functions bai.in and bai.out). See their help pages for further information.
6 Descriptive Statistics for Detrended Data
It is also easy in dplR to compute commonly used descriptive statistics that describe the correlation between series (both within and between tree correlations) as well as the expressed population signal and signal-to-noise ratio for a data set. These are done in dplR using the rwi.stats function so-named because these statistics are typically (but not always) carried out on detrended and standardized ring widths (rwi). If a data set has more than one core taken per tree this information can be used in the calculations to calculate within vs. between tree correlation. The function read.ids is used to identify which trees have multiple cores.
ca533.ids <- read.ids(ca533, stc = c(3, 2, 1))
rwi.stats(ca533.rwi, ca533.ids, prewhiten=TRUE)## n.cores n.trees n n.tot n.wt n.bt rbar.tot rbar.wt rbar.bt c.eff
## 1 34 21 11.7 523 13 510 0.444 0.603 0.439 1.448
## rbar.eff eps snr
## 1 0.501 0.922 11.75There is (at least) one other way of looking at the average interseries correlation of a data set. The interseries.cor function in dplR gives a measure of average interseries correlation that is different from the rbar statistics from rwi.stats. In this function, correlations are calculated serially between each tree-ring series and a master chronology built from all the other series in the rwl object (leave-one-out principle). The average of those correlations is sometimes called the “overall interseries correlation” or even the “COFECHA correlation” in reference to commonly used crossdating software COFECHA. This number is typically higher than the various rbar values given by rwi.stats. We are showing just the first five series and the mean for all series here:
ca533.rho <- interseries.cor(ca533.rwi, prewhiten=TRUE,
method="spearman")
ca533.rho[1:5, ]## res.cor p.val
## CAM011 0.5358 0
## CAM021 0.6760 0
## CAM031 0.5258 0
## CAM032 0.6265 0
## CAM041 0.4907 0mean(ca533.rho[, 1])## [1] 0.6368Again, if these concepts are unknown to you statistically look at some of the canonical works in dendrochronology like Cook et al. (1990), Fritts (2001), and Hughes et al. (2011).
7 Building a Mean Value Chronology
After detrending, a user will typically build a chronology by averaging across the years of the rwi object. In dplR the function for doing this is chron which by default uses Tukey’s biweight robust mean (an average that is unaffected by outliers).
ca533.crn <- chron(ca533.rwi, prefix = "CAM")This object has the same number of rows as the rwi object that was used as the input and two columns. The first gives the chronology and the second the sample depth (the number of series available in that year).
dim(ca533.rwi)## [1] 1358 34dim(ca533.crn)## [1] 1358 2An object produced by chron has a generic S3 method for plotting which calls the crn.plot function (which has many arguments for customization). Here we will just make a simple plot of the chronology with a smoothing spline (function ffcsaps) added.
plot(ca533.crn, add.spline=TRUE, nyrs=20)
8 Conclusion
In general this page aims to give a very cursory overview of basic tasks that most dendrochronologists will want to be aware of. Know that we are just scratching the surface of what dplR is capable of. As a small example, here is a way that a user might decide to truncate a chronology based on the subsample signal strength.
ca533.sss <- sss(ca533.rwi,ca533.ids)
yr <- time(ca533)
par(mar = c(2, 2, 2, 2), mgp = c(1.1, 0.1, 0), tcl = 0.25, xaxs='i')
plot(yr, ca533.crn[, 1], type = "n", xlab = "Year",
ylab = "RWI", axes=FALSE)
cutoff <- max(yr[ca533.sss < 0.85])
xx <- c(500, 500, cutoff, cutoff)
yy <- c(-1, 3, 3, -1)
polygon(xx, yy, col = "grey80")
abline(h = 1, lwd = 1.5)
lines(yr, ca533.crn[, 1], col = "grey50")
lines(yr, ffcsaps(ca533.crn[, 1], nyrs = 32), col = "red", lwd = 2)
axis(1); axis(2); axis(3);
par(new = TRUE)
## Add SSS
plot(yr, ca533.sss, type = "l", xlab = "", ylab = "",
axes = FALSE, col = "blue")
abline(h=0.85,col="blue",lty="dashed")
axis(4, at = pretty(ca533.sss))
mtext("SSS", side = 4, line = 1.1, lwd=1.5)
box()
We hope that this helps users cover introductory data handling and processing using dplR and R. As we noted above we are just providing a short introduction as to what is possible in dplR. There are many other functions in dplR that will help users analyze tree rings. These include a host of functions for statistical cross dating as well as spectral and wavelet analysis. You can see other vignettes for dplR via:
vignette(package="dplR")And if you were interested in ARMA models, wavelets, and spectral analysis for tree rings you might want to display a PDF of the vignette on time-series analysis via:
vignette(topic="timeseries-dplR",package="dplR")Or just extract the code via:
edit(vignette(topic="timeseries-dplR",package="dplR"))9 Modules
In the following pages we will cover the basics of dplR, crossdating, and some limited time-series analysis. Users should be familiar with the basics of dendrochronology and concepts like detrending, autocorrelation, spectral analysis and so on.
If this is all new to you – you should proceed immediately to a good primer on dendrochronology like Fritts (2001) or the Cook Book (1990). These pages are not intended to teach you about how to do tree-ring analysis. They are intended to teach you how to use R for dendro.
9.1 Pages
References
Bunn AG (2008) A dendrochronology program library in R (dplR). Dendrochronologia, 26, 115–124.
Bunn AG, Jansma E, Korpela M, Westfall RD, Baldwin J (2013) Using simulations and data to evaluate mean sensitivity (\(\zeta\)) as a useful statistic in dendrochronology. Dendrochronologia, 31, 250–254.
Cook E, Kairiukstis L (1990) Methods of dendrochronology: Applications in the environmental sciences. Kluwer, pp.
Cook E, Briffa K, Shiyatov S, Mazepa V (1990) Tree-Ring Standardization and Growth-Trend Estimation. In: Methods of dendrochronology: Applications in the environmental sciences (eds Cook E, Kairiukstis L), pp. 104–123. Kluwer, Dordrecht.
Fritts HC (2001) Tree Rings and Climate. The Blackburn Press, pp.
Hughes MK, Swetnam TW, Diaz HF (eds.) (2011) Dendroclimatology, Vol. 11. Springer Netherlands, Dordrecht, pp.