Getting Started with dplR

Reading ring widths

Ring widths are read with read.rwl(), which handles Tucson (decadal), compact, Heidelberg, CSV-style spreadsheets and TRiDaS files. It guesses the format, but naming it is safer:

dat <- read.rwl("mysite.rwl", format = "tucson")

The result is an rwl object: a data frame with one column per series and one row per year, with the years as row names and NA where a series has no ring. Here we use co021, 35 Douglas-fir series from Mesa Verde, Colorado, which was read from the International Tree-Ring Data Bank this way.

data(co021)
class(co021)

## [1] "rwl"        "data.frame"

dim(co021)

## [1] 788  35

co021[1:5, 1:4]

##      641114 641121 641132 641143
## 1176     NA     NA     NA     NA
## 1177     NA     NA     NA     NA
## 1178     NA     NA     NA     NA
## 1179     NA     NA     NA     NA
## 1180     NA     NA     NA     NA

Describing the data

rwl.check() gives an overview: where the file came from, number of series, span, mean length, mean correlation with the master, and then anything worth a look, from missing rings to suspicious values and series that may be misdated.

rwl.check(co021)

## rwl.check: co021
## File: co021.rwl
## Site: not given; the file carries no header
## Precision: 0.01 mm
## -------------
## 35 series, 19772 measurements, 1176-1963, 0.01 mm precision
## Mean series length: 564.9
## Mean (SD) AR1: 0.604 (0.157)
## Mean (SD) correlation with the master: 0.892 (0.02)
## -------------
## 0 error(s), 0 warning(s), 47 note(s)
## -------------
## NOTE
##   RWL_CONSECUTIVE_ZERO (45)
##     641114: 2 consecutive absent rings, 1505-1506
##     641114: 2 consecutive absent rings, 1584-1585
##     641121: 2 consecutive absent rings, 1584-1585
##     641132: 2 consecutive absent rings, 1584-1585
##     641143: 2 consecutive absent rings, 1584-1585
##     642114: 2 consecutive absent rings, 1505-1506
##     642114: 2 consecutive absent rings, 1584-1585
##     642121: 2 consecutive absent rings, 1584-1585
##     642143: 2 consecutive absent rings, 1584-1585
##     642211: 2 consecutive absent rings, 1584-1585
##     ... and 35 more; see as.data.frame()
##   RWL_LOW_DEPTH (1)
##     1176-1199 rests on fewer than 2 series
##   RWL_ZERO_RING (1)
##     716 ring(s) of zero width in 35 series

summary() returns per-series statistics, and plot() shows where each series sits in time.

head(summary(co021))

##   series first last year  mean median stdev  skew kurtosis  gini   ar1
## 1 641114  1270 1963  694 0.287   0.23 0.231 2.884   13.963 0.372 0.686
## 2 641121  1250 1963  714 0.328   0.26 0.315 3.306   15.199 0.410 0.744
## 3 641132  1256 1963  708 0.357   0.29 0.337 4.741   33.324 0.373 0.686
## 4 641143  1237 1963  727 0.344   0.27 0.287 2.341    7.637 0.397 0.708
## 5 642114  1243 1963  721 0.281   0.24 0.219 2.848   12.385 0.358 0.673
## 6 642121  1260 1963  704 0.313   0.21 0.416 4.399   22.962 0.474 0.865

plot(co021, plot.type = "spag")

Detrending

Raw ring widths carry an age-related growth trend and differences in mean growth between trees. Detrending fits a curve to each series and divides the ring widths by it, giving dimensionless indices with a mean of about one. detrend() does this for every series; detrend.series() does one and plots the fit, which is a good way to choose a method.

x <- co021[, "641114"]
names(x) <- rownames(co021)
x.rwi <- detrend.series(x, method = c("Spline", "ModNegExp"),
                        make.plot = TRUE)

Here we use a cubic smoothing spline with a 50% frequency cutoff at two-thirds of each series’ length (the default for "Spline").

co021.rwi <- detrend(co021, method = "Spline")

The result is an rwi object with the same shape as the rwl, which records how it was made. summary() describes the indices as a collection. rbar.eff is the mean interseries correlation and EPS the expressed population signal; an EPS above about 0.85 is the usual rule of thumb for a chronology that represents the population. These come from rwi.stats(), and here count each core as its own tree; pass ids to group cores by tree. The summary also correlates each series with the mean of the others and lists any that do not fit.

summary(co021.rwi)

## Ring-width indices: 35 series, 1176-1963
## Made by detrend(), method "Spline", as ratios (centred on 1).
## Common interval: none (no year has every series)
## rbar.eff 0.737, EPS 0.990, SNR 97.85 (rwi.stats(), one tree per series)
## Series vs the others (interseries.cor()): mean r 0.881, range 0.803 to 0.923
## Every series correlates with the others at p < 0.05.
## No index is far above the rest of its series (robust z > 5, log scale).
## More:
##   as.data.frame(summary(x))  every series: span, mean, sd, ar1, cor, p, extreme
##   summary(x, ids = )         rbar and EPS with cores grouped by tree
##   rwi.stats.running(x)       rbar and EPS through time
##   corr.rwl.seg() on widths   where in a series the fit breaks down

plot() with plot.type = "image" shows every index at once, years across and series down, brown below 1 and green above. Vertical stripes are years the trees agree on, which is the signal a chronology is built from.

plot(co021.rwi, plot.type = "image")

The same plot is a quick check on the detrending. Dividing each series by its mean leaves the age trend in, and it shows as a green run at the start of nearly every series.

plot(detrend(co021, method = "Mean"), plot.type = "image")

Building a chronology

chron() averages the indices by year, using Tukey’s biweight robust mean by default. The result is a crn object: the chronology and the number of series behind each year.

co021.crn <- chron(co021.rwi)
tail(co021.crn)

##            std samp.depth
## 1958 1.1284148         32
## 1959 0.1459389         32
## 1960 1.1804609         32
## 1961 0.7720185         32
## 1962 0.6232634         32
## 1963 0.6440584         29

plot(co021.crn, add.spline = TRUE, nyrs = 32)

The early part of this chronology rests on few series (fewer than five before 1234), so check samp.depth before trusting a given year.

Checking the crossdating

Crossdating assigns each ring its exact calendar year. dplR does not replace visual crossdating, but it can check it statistically, the way COFECHA does. To see what an error looks like, we plant one: the 1500 ring of series 641143 is deleted, so every ring before 1500 is now dated one year too late.

dat <- co021
x <- dat[, "641143"]
names(x) <- rownames(dat)
dat[, "641143"] <- delete.ring(x, year = 1500)

corr.rwl.seg() correlates overlapping segments of each series against a master built from all the other series. Segments that do not correlate significantly are flagged, and with lag.max it also finds the shift at which each segment correlates best.

crs <- corr.rwl.seg(dat, seg.length = 50, pcrit = 0.01,
                    lag.max = 10, label.cex = 0.7)

Every tested segment of 641143 that ends before 1500 correlates best at a lag of -1, and every later one at a lag of 0. In dplR, as in COFECHA, a negative lag means the series is missing a ring, so this points to a missing ring near 1500.

ok <- !is.na(crs$best.lag["641143", ])
crs$best.lag["641143", ok]

## 1275.1324 1300.1349 1325.1374 1350.1399 1375.1424 1400.1449 1425.1474 1450.1499 
##        -1        -1        -1        -1        -1        -1        -1        -1 
## 1475.1524 1500.1549 1525.1574 1550.1599 1575.1624 1600.1649 1625.1674 1650.1699 
##         0         0         0         0         0         0         0         0 
## 1675.1724 1700.1749 1725.1774 1750.1799 1775.1824 1800.1849 1825.1874 1850.1899 
##         0         0         0         0         0         0         0         0 
## 1875.1924 1900.1949 
##         0         0

ccf.series.rwl() looks at one series in more detail, plotting the full cross-correlation with the master for each segment.

ccf <- ccf.series.rwl(rwl = dat[, colnames(dat) != "641143"],
                      series = dat[, "641143"],
                      series.yrs = as.numeric(rownames(dat)),
                      seg.length = 50, bin.floor = 50)

## Negative lags suggest a missing ring in the series

xdate.report() puts all of this into a COFECHA-style report. It lists the flagged segments with their best lag and the gain in correlation, and records the settings and file checksum so the report can be reproduced. Printing it shows the report in the console.

rpt <- xdate.report(dat, title = "co021 with a planted fault")
rpt

The report opens with a summary of the collection:

Summary

The flagged segments section is where to start. All nine flags are on 641143, all are B flags at a lag of -1, and all end before 1500:

Flagged segments

The same table is in the report as a data frame, rpt$flagged, ready for further work:

rpt$flagged[, c("series", "from", "to", "flag", "best.lag", "gain")]

##   series from   to flag best.lag      gain
## 1 641143 1250 1299    B       -1 1.0043412
## 2 641143 1275 1324    B       -1 0.7509562
## 3 641143 1300 1349    B       -1 0.6912424
## 4 641143 1325 1374    B       -1 0.7945452
## 5 641143 1350 1399    B       -1 0.6579816
## 6 641143 1375 1424    B       -1 0.7378189
## 7 641143 1400 1449    B       -1 0.8823181
## 8 641143 1425 1474    B       -1 1.0267152
## 9 641143 1450 1499    B       -1 0.9913002

A B flag means the segment correlates better at another lag than at its dated position. An A flag, which does not appear here, means the dated position is the segment’s best but falls short of significance. The report also has a correlation table for every segment of every series, descriptive statistics, the output of rwl.check(), and notes on how each number was computed:

write.xdate.report() saves the report as text, Markdown or HTML:

write.xdate.report(rpt, "co021-report.html")

The statistics say where to look; the wood says what happened. Having found the likely missing ring, you would go back to the sample, and once confirmed, fix the series with insert.ring().

Where to go next

• Learning to Love dplR covers these steps and more in depth.
• vignette("math-dplR") gives the mathematics behind the smoothing splines, detrending and other functions.
• help(package = "dplR") lists every function. Good next stops are ?detrend, ?chron, ?ssf (signal-free chronologies), ?corr.rwl.seg and ?xdate.report.
• dplR is part of openDendro.