Computes the correlation between each tree-ring series in a rwl object.
corr.rwl.seg(rwl, seg.length = 50, bin.floor = 100, n = NULL,
nyrs = NULL, prewhiten = TRUE, ar.order.max = NULL,
pcrit = 0.05, biweight = TRUE,
method = c("spearman", "pearson","kendall"),
make.plot = TRUE, label.cex = 1, floor.plus1 = FALSE,
master = NULL, lag.max = 0,
master.yrs = as.numeric(if (is.null(dim(master))) {
names(master)
} else {
rownames(master)
}),
...)A list containing matrices spearman.rho,
p.val, overall, bins,
rwi, vector avg.seg.rho,
numeric seg.lag, seg.length, pcrit,
label.cex, matrices best.lag,
best.rho and numeric lag.max. An additional character
flags is also returned if any segments fall below the
critical value. Matrix spearman.rho contains the
correlations for each series by bin. Matrix p.val
contains the p-values on the correlation for each series by
bin. Matrix overall contains the average correlation and
p-value for each series. Matrix bins contains the years
encapsulated by each bin. The vector avg.seg.rho
contains the average correlation for each bin. Matrix rwi
contains the detrended rwl data, the numerics seg.lag,
seg.length, pcrit, label.cex
are from the oroginal call and used to pass into plot.crs.
Matrices best.lag and best.rho have the same
shape and names as spearman.rho. best.lag
gives, for each series and bin, the lag (in years) at which the
segment correlates best with the master, and best.rho
gives that correlation. Both are NA where
spearman.rho is. With lag.max = 0,
best.lag is 0 and best.rho equals
spearman.rho. lag.max is from the original
call. Flags in flags depend on the p-value only and are
the same whatever lag.max is.
a data.frame with series as columns and years as
rows such as that produced by read.rwl.
an even integral value giving length of segments in years (e.g., 20, 50, 100 years).
a non-negative integral value giving the base for locating the first segment (e.g., 1600, 1700, 1800 AD). Typically 0, 10, 50, 100, etc.
NULL or an integral value giving the filter length
for the hanning filter used for removal of low
frequency variation.
NULL or a number greater than zero. If not
NULL, each series is divided by a smoothing spline with this
rigidity (see caps) to remove low-frequency variation,
as detrend does with method = "Spline". A value
of 1 or less is taken as a proportion of each series' length.
Unlike the hanning filter, the spline removes no years
from the ends of a series. Cannot be combined with n.
logical flag. If TRUE each series is
whitened using ar.
NULL or a positive integer giving the
maximum order of the ar model used to prewhiten.
Prewhitening removes as many years from the start of each series as
the order of the model. If NULL, the order is chosen by AIC
up to the ar default, which on long series can be 20
or more. Requires prewhiten = TRUE.
a number between 0 and 1 giving the critical value for the correlation test.
logical flag. If TRUE then a robust
mean is calculated using tbrm.
Can be either "pearson", "kendall", or
"spearman" which indicates the correlation coefficient to be
used. Defaults to "spearman". See cor.test.
logical flag indicating whether to make a
plot.
numeric scalar for the series labels on the
plot. Passed to axis.cex in axis.
logical flag. If TRUE, one year is
added to the base location of the first segment (e.g., 1601, 1701,
1801 AD).
NULL, a numeric vector or a
matrix-like object of numeric values, including a
data.frame. If NULL, a number of master chronologies,
one for each series in rwl, is built from
rwl using the leave-one-out principle. If a
vector, the function uses this as the master chronology. If
a matrix or data.frame, this object is used for
building the master chronology (no leave-one-out).
a non-negative whole number less than
seg.length. Each segment is also correlated against the
master shifted by every lag from -lag.max to
lag.max years, and the best position is returned in
best.lag and best.rho. The default, 0,
tests the dated position only, as before.
a numeric vector giving the years of
series. Defaults to names or rownames of
master coerced to numeric type.
other arguments passed to plot.
Andy Bunn. Patched and improved by Mikko Korpela.
This function calculates correlation serially between each tree-ring
series and a master chronology built from all the other series in the
rwl object (leave-one-out principle). Optionally, the
user may give a master chronology (a vector) as an
argument. In the latter case, the same master chronology is used for
all the series in the rwl object. The user can also
choose to give a master data.frame (series as
columns, years as rows), from which a single master chronology is
built.
Correlations are done for each segment of the series where segments
are lagged by half the segment length (e.g., 100-year segments would
be overlapped by 50-years). The first segment is placed according to
bin.floor. The minimum bin year is calculated as
ceiling(min.yr/bin.floor)*bin.floor where
min.yr is the first year in either the rwl
object or the user-specified master chronology, whichever
is smaller. For example if the first year is 626 and
bin.floor is 100 then the first bin would start in 700.
If bin.floor is 10 then the first bin would start in 630.
Correlations are calculated for the first segment, then the second
segment and so on. Correlations are only calculated for segments with
complete overlap with the master chronology. For now, correlations are
Spearman’s rho calculated via cor.test using
method = "spearman".
Each series (including those in the rwl object) is optionally
detrended as the residuals from a hanning filter with
weight n. The filter is not applied if n is
NULL. Detrending can also be done via prewhitening where the
residuals of an ar model are added to each series
mean. This is the default. The master chronology is computed as the
mean of the rwl object using tbrm if
biweight is TRUE and rowMeans if not. Note
that detrending can change the length of the series. E.g., a
hanning filter will shorten the series on either end by
floor(n/2). The prewhitening default will change the
series length based on the ar model fit. The effects of
detrending can be seen with series.rwl.plot.
As an alternative to the hanning filter, nyrs
divides each series by a smoothing spline, which removes no years from
the ends of the series. COFECHA uses a 32-year spline. The spline
alone does not usually save years at the start of a series, because
the ar model is then fitted to a high-pass filtered series,
and AIC tends to choose a high order for it. Setting
ar.order.max to a small value (e.g., 3) limits the years
lost to prewhitening, and after a spline a low-order model loses little
of the fit.
With neither filter, each series is divided by its mean before the
master is built, so every series counts equally. Data that can be
negative, such as indices from detrend with
difference = TRUE, log widths or isotope values, have
the mean subtracted instead, since dividing by a negative mean would
flip the series. The n and nyrs filters
divide by a smooth curve, so they refuse such data. The same applies
to the other crossdating functions.
The function is typically invoked to produce a plot where each segment for each series is colored by its correlation to the master chronology. Green segments are those that do not overlap completely with the width of the bin. Blue segments are those that correlate above the user-specified critical value. Red segments are those that correlate below the user-specified critical value and might indicate a dating problem.
A segment that correlates well where it is dated can still correlate
better somewhere else. With lag.max greater than 0, each
segment is also correlated against the master at every shift from
-lag.max to lag.max years. Each year of the
segment is paired with the master value k years away, so the
correlation at lag k is that of series[t] with
master[t + k]. The sign follows ccf.series.rwl
with series.x = FALSE: a negative lag means missing rings
in the series, a positive lag means false (extra) rings. The dated
position wins ties.
This gives the two kinds of flag COFECHA reports. An ‘A’
segment correlates below the critical value, but its dated position
is still the best one tested: it is weak, not misdated. A
‘B’ segment correlates better at some other position,
whether or not it is significant as dated. The letters are not
stored; they follow from the returned matrices as
best.lag != 0 for B and
best.lag == 0 & p.val >= pcrit for A (see Examples). When
lag.max is greater than 0, the plot draws B segments in
purple and leaves red for the rest of the segments below the critical
value.
A B flag is a hypothesis about the dating, not a verdict. A one-year
shift that gains 0.05 in correlation is noise as often as it is a
missing ring, which is why the size of the gain is returned
(best.rho - spearman.rho) alongside the lag, so that it can be
thresholded. Check any B segment against the wood, for example with
ccf.series.rwl and xskel.ccf.plot.
Read the lags along a series rather than one bin at a time. Dating
runs from the bark inward, so a missing ring moves every ring before
it one year late: every complete bin before the missing ring has
best.lag of -1 and every bin after it has 0. A false ring
does the same with +1. The error therefore lies where the lag
changes, not in the bin with the largest gain. The bin that
straddles it holds some shifted and some unshifted years, and may
show either lag with a small gain. A run of -1 with correctly dated
bins on both sides is a different fault: the ring count is right, but
a ring is missing at the later end of the run and there is an extra
ring at the earlier end. A locally absent ring entered as a zero in
too early a year looks exactly like this. A test on the whole series,
such as RWL_DATING_LAG in rwl.check, usually cannot
see it, because most of the series is dated correctly. The Examples
show both patterns.
The shifted correlations follow the same complete-overlap rule as the
dated one: a lag is tested only if the master has a value for every
year of the shifted window. Near the start and end of the record
(and of the master), the lags that run off the record cannot be
tested, and a segment that lies within lag.max years of
an end is only searched in one direction. A segment that is not
complete at its dated position is not tested at all. Floating series
and series with dating errors near their ends sit exactly where the
test is blind, so the absence of a flag there says little. Worse, a
segment misdated in the direction that runs off the record can still
be flagged B, at whichever of the remaining lags happens to correlate
best, so the lag reported for a B near either end can be noise even
when the segment really is misdated.
corr.series.seg, skel.plot,
series.rwl.plot, ccf.series.rwl,
plot.crs
library(utils)
data(co021)
crs <- corr.rwl.seg(co021, seg.length = 100, label.cex = 1.25)
names(crs)
## Average correlation and p-value for the first few series
head(crs$overall)
## Average correlation for each bin
crs$avg.seg.rho
## Plant a missing ring and follow it through the bins: delete the
## 1500 ring of series 641143, the same fault as in the examples for
## ccf.series.rwl(), corr.series.seg() and xskel.ccf.plot(). Dated
## from the bark, every ring before 1500 now sits one year late.
dat <- co021
x <- dat$"641143"
names(x) <- rownames(dat)
dat$"641143" <- delete.ring(x, year = 1500)
crs <- corr.rwl.seg(dat, lag.max = 5)
## -1 in every complete bin before 1500 (drawn in purple), 0 from 1500
## on, and a large gain in each shifted bin
ok <- !is.na(crs$best.lag["641143", ])
rbind(lag = crs$best.lag["641143", ok],
rho = round(crs$spearman.rho["641143", ok], 2),
best.rho = round(crs$best.rho["641143", ok], 2))
## The same missing ring with a false ring inserted at 1400. The ring
## count is right again, and only the rings between the two errors are
## one year late.
dat$"641143" <- insert.ring(delete.ring(x, year = 1500), year = 1400)
crs <- corr.rwl.seg(dat, lag.max = 5, make.plot = FALSE)
## -1 between the two errors only, 0 on both sides
crs$best.lag["641143", ok]
## COFECHA-style A and B flags for every series
flag <- ifelse(crs$best.lag != 0, "B",
ifelse(crs$p.val >= crs$pcrit, "A", ""))
flag[is.na(flag)] <- ""
## B segments with their lag and the gain in correlation
idx <- which(flag == "B", arr.ind = TRUE)
data.frame(series = rownames(flag)[idx[, 1]],
bin = colnames(flag)[idx[, 2]],
lag = crs$best.lag[idx],
rho = round(crs$spearman.rho[idx], 3),
best.rho = round(crs$best.rho[idx], 3))
Run the code above in your browser using DataLab