powt() now defaults to method = "cook", which fits a power to each
series as Cook and Peters (1997) did. The default was
method = "universal", which fits one power to every series in a
mixed model with no term for each series. Trees differ in how strongly
their spread depends on their level, and one power can leave that
dependence in some series while over-correcting others. In simulated
collections where half the trees have multiplicative noise and half
additive, the universal power under-corrects the first group and
over-corrects the second, and the per-series powers do well on both.
The universal method's year term does help on some data: on co021,
where climate moves the level and spread of neighboring rings
together, it leaves less of the dependence than the per-series fit
does. Try both. To get the old results, pass method = "universal". A single series now works with the default.
dplPy's powt() already defaulted to "cook".
powt(method = "cook") log transforms a series whose spread grows at
least as fast as its level (fitted power 1 - b <= 0), with a warning
that names it, as method = "universal" already did. It used the
absolute value, abs(1 - b), so a slope of 1.4 became a power of 0.4,
a weaker transform than the series needed, with no message. The power
returned with return.power = TRUE is now the fitted 1 - b, so a
value at or below zero marks a logged series. Zero rings in a logged
series become -Inf, and a second warning names those series.
powt() guessed the measurement precision, which it uses in place of
zero spreads and levels before taking logs, from the decimals of one
value padded by format(). A single value with many decimals, such as
a ring filled by fill.internal.NA(), moved the precision for the
whole collection. On wa082 with series 712011's missing 1900 ring
filled, the precision dropped from 0.01 to 1e-7, each of the 102 pairs
of equal rings got a log spread near -16 where it should have been
-4.6, and the median method = "cook" power fell from 0.28 to 0.11
(the universal power from 0.22 to 0). The precision now comes from the
most common number of decimals, so one odd value doesn't change it.
When most values have more than eight decimals, as transformed data
do, powt() warns that it can't tell the precision.
bai.out() warns, and names the series, when a series' ring widths
add up to more than the radius (diam / 2). Working in from the bark,
the radius then reaches zero before the innermost rings are used up,
and those rings get areas that are too small or negative. This used to
happen silently: on gp.rwl with the diameters inside the bark from
gp.dbh, 65 rings in five series came out negative, and two more
series had inner areas that were too small. The help page says so and
points to bai.in() and bakker().
bakker() returns data frames with the years and series of rwl, as
its help page said it did, and the ring areas as class bai, with the
dplR.bai record, like bai.in() and bai.out(). Both parts used to
be matrices.
caps(), and so detrend(method = "Spline"), takes a negative nyrs
as a percent of the series length, as ARSTAN does: nyrs = -67 is a
67% spline, the same as nyrs = 0.67. This is the rule dplPy uses.
The crossdating functions that take nyrs (corr.rwl.seg(),
corr.series.seg(), ccf.series.rwl(), interseries.cor(),
rwi.stats.running(), series.rwl.plot(), xskel.plot(),
xskel.ccf.plot(), xdate.floater() and xdate.report()) accept
it too, resolved against each series' own length. A series too short
for the spline is named in the error.
caps() now stops when nyrs, after converting a proportion or
percent, truncates to less than a 2-year spline. A 0-year spline used
to return NaN and a 1-year spline a nearly flat line, with no
message. The warning for an nyrs longer than the series now checks
the converted value.
xdate.report() stops on an invalid nyrs (0, for example). Its
per-series spline screen used to turn the error into an excluded
series, so every series was excluded and the report came back with
no segment correlations and only a note saying no series could be
crossdated.
morlet() has a lag1 argument, the lag-one autocorrelation of the
noise its significance levels test against. It was set to 0 inside the
function, so the contours wavelet.plot() draws have always been a
test against white noise, and the help page didn't say so. A
chronology with persistence has more power at long periods than white
noise, so long periods passed easily: on the co021 standard
chronology, 47% of the plane inside the cone of influence at periods
over 32 years passes at the 95% level against white noise and 23%
against red noise with the chronology's own lag-one autocorrelation.
The default is still 0, so existing plots don't change. lag1 = NULL
estimates it from the series. The help pages for morlet() and
wavelet.plot() say what the test is.
redfit() printed its warning for an AR(1) estimate above 1 as a
vector, c("redfitGetrho returned rho = %f, forced to zero", "1.0055..."), because the message was built with gettext() in
place of gettextf(). It now reads "redfitGetrho returned rho =
1.005554, forced to zero". A heavily smoothed series, such as a
32-year pass.filt() low-pass of a chronology, triggers it.
rwl.stats(), and so summary() of an rwl object, no longer stops on
a series with no values. Such series are what x[rows, ] leaves when a
collection is cut to years some series don't reach, which [ keeps on
purpose, and the whole table used to fail with "'ts' object must have
one or more observations", which named neither the series nor the
cause. The series now keeps its row, with a length of 0 and NA for
every statistic, and a warning names it and points to window() and
subset(). Of 14 common functions tried on co021 cut this way,
these two were the only ones that failed.
pass.filt() stops with a plain message when W can't be used.
A period longer than half the series used to fail with "only 0's may
be mixed with negative subscripts" (the series is padded at each end by
twice the longest period). A two-element W that mixed a frequency and
a period, such as c(0.05, 5), was read as two periods and failed
inside signal::butter(). Frequencies at or above 0.5 and periods of 2
or less now stop too, with a message naming the limit. The help page
says so.
summary() of ring-width indices lists any ratio index far above the
rest of its own series (a robust z-score of the log index over
zcrit, default 5), with its year and whether it is at the start, end
or middle of the series. A growth curve that comes close to zero at an
end of a series inflates the ratio there (Cook and Peters 1997), and
nothing in the summary showed it: detrend(wa082, method = "AgeDepSpline") ends 712082 at an index of 11.95 over a curve of
0.08 mm and the collection statistics hardly move. Every index the
check lists in the on-board data sets is in the first or last few
rings of a series. Differences are not tested. The per-series table
gains a column extreme.
xdate.report(prewhiten = FALSE) tested no segments at all unless
ar.order.max = NULL was also given, because the default cap of 3 had
nothing to cap. It now steps aside when prewhitening is off, as it does
in xdate.floater().
xdate.floater() now filters and correlates the way xdate.report()
does by default: a 32-year spline (nyrs = 32), an AR model of order
at most 3 (ar.order.max = 3) and Pearson's correlation. The old
defaults were no low-frequency filter, an AR order chosen by AIC with no
cap, and Spearman's rho, so a floater was dated one way and then checked
another. Pass nyrs = NULL, ar.order.max = NULL, method = "spearman"
to get the old behaviour. Giving n (the Hanning filter) without
nyrs still works: the spline default steps aside.
xdate.floater() placed a floater too early by half the Hanning filter
length whenever n was given: with n = 32, 16 years early. The filter
takes n/2 rings off each end of the series, and the end year was worked
out as if it took none off the end. Prewhitening and the spline were not
affected.
rwl.report() is deprecated and will be removed in a future release.
rwl.check() now opens its printed result with the same description
of the collection: the file and site, the number of series and
measurements, the span, the mean series length, the mean (SD) AR1 and
the mean (SD) correlation of each series with the master. Its findings
already covered the zero rings, internal NAs and small or big rings
the report listed, and the two did not always agree: the report counted
a run of filler zeros as missing rings, which rwl.check() reports as
RWL_ZERO_RUN. The correlation is the one the crossdating checks
compute, so it costs nothing extra and differs a little from the
interseries.cor() figure the report gave (0.718 against 0.675 for
gp.rwl). The per-series values are in the new series element of the
result. rwl.check() now warns when given indices or basal area
increment, as rwl.report() did.
rwl.check() on a file it could not read printed "0 error(s)": the
RWL_CHECK_ERROR finding had no severity and was not counted. It is
now an error. Printing that result, and summary() of it, also
stopped with an error.
rwl.check() no longer runs the file checks on a name it made up from
the call. rwl.check(r) checked the R/ folder (macOS ignores case),
and an object named like a file in the working directory was checked
against that file. The file checks run when x is a path or file is
given.
New function rwl.meta() returns the metadata an rwl, rwi or bai
object carries: the file and the reader; from a Tucson header, the
site, species, state or country, elevation, coordinates (as written)
and investigators; the span the header declares beside the span of the
data, the precision, interior gaps, series the reader renamed,
problems met while reading and, for indices, how they were made. It
used to be reachable only through attr(x, "dplR.provenance") and an
internal header parser. With no record it returns NULL with a
message saying why there might be none.
rwl.report() no longer parses the comment lines of a file read by
read.sheet() as an ITRDB header. They gave a site line made from
fragments of the comment, such as "Site: # Site urricane Ridge".
New data set me044: eastern hemlock at Howland Forest, Maine (ITRDB
ME044, Maxwell, Belmecheri and Taylor), 42 series from 21 trees,
1780-2012, measured to 0.001 mm. It has a full ITRDB header and was
read with the default read.tucson() settings, so it shows what the
provenance record holds. Many of the trees grow much faster from the
1910s on, and the file has 29 zero rings.
ssf(return.info = TRUE) returns a plain list. It had class "crn"
and "data.frame", so print() and head() showed 0 rows and
nrow() and summary() gave nonsense about a chronology that was
not there. plot.crn() on the list now plots the final chronology,
ssfCrn; it failed with "'x' and 'y' lengths differ". Code that
takes elements out of the list with $ or [[ is unaffected.
When the curve fitted in ssf() reaches zero or below, the error
now names the series, the iteration and the year, gives the lowest
value of the chronology and its year, and explains the cause: in sites
with many absent rings, dividing by a low chronology inflates the
signal-free measurements, and putting them back on the series mean
pushes the narrowest rings below zero. It says that recode.zeros
does not prevent this, and that removing the series named usually
moves the failure to others. The error used to say only that the
curve had values <= 0. Reported by a dplPy user.
ssf() stops before any fitting when a series has missing values
inside it, with error [7] naming each series and the years. It used to
stop partway through with "'NA's are not allowed in the middle of the
series", which named neither.
ssf() has a new argument, return.unconverged (default FALSE). By
default ssf() still stops with error [2] when the stopping criterion
is not met within max.iterations, and returns nothing. With
return.unconverged = TRUE it returns the last iteration anyway, marked
attr(, "converged") = FALSE (and converged in the return.info
list), so a runaway can be shown in a reproducible example. It warns
when the argument is set and again when the result is unconverged; an
unconverged chronology is not for analysis. CRUST returns unconverged
chronologies without a warning. The help page now lists the convergence
vector as MAD_Vec, the name the function returns (it said MAD_out).
ssf() has a new argument, curve.floor (default NULL, off). If set,
every fitted curve, the first fit included, is raised to at least that
value, as CRUST and dplPy's preset = "crust" do with 0.02 mm. It lets
ssf() run where a spline bends to zero at the end of a series
(nm046, series 644081). A floored curve can give very large indices
(19 for 644081 in 1969), so ssf() warns with the series, years and
largest index, and returns the floored years in attr(, "floored") and
flooredCurves. Error [1] now mentions the option.
ssf() has two more arguments for CRUST's steps, both off by default.
rescale = "ratio" puts each signal-free series back on its mean by
multiplying instead of adding a constant. crn.guard keeps the raw
measurements in years where the chronology is below it (CRUST uses
0.01), warns with the years, and lets a year with every ring absent
through instead of stopping with error [3] or [4] (ca533). With
rescale = "ratio", crn.guard = 0.01, curve.floor = 0.02 and
max.iterations = 40, ssf() matches dplPy's ssf(preset = "crust").
rescale = "ratio.masked" is CRUST's version of the ratio: the means
use only years with more than one series, and only those years are
scaled.
ssf() gains biweight (default TRUE); FALSE uses the arithmetic
mean for every chronology, as CRUST does by default. max.iterations = 1
is allowed with return.unconverged = TRUE. With these and the CRUST
arguments above, ssf() reproduces CRUST at a fixed number of
iterations to about 1e-5 or better (see the help page).
strip.rwl() accepts ring-width indices (class "rwi") and strips
them as they are, returning indices. A series is then judged on the
variation the chronology will keep. Indices used to be treated as
widths and detrended again, with a warning.
strip.rwl() gains nyrs, the stiffness of the spline that raw
widths are detrended with before stripping, and detrends once. A
message says which spline was used and that stripping indices is the
better route. It
used to apply a 20-year and then a 200-year spline to the same series.
That was a misreading of Fowler and Boswijk (2003), who made two
separate sets of indices and stripped each, and it came to the 20-year
spline alone, so the default of 20 keeps the old behaviour: the
indices differ in the third decimal place or beyond, and the series
left out are the same in the data sets tried. nyrs = 200 gives the
paper's other version. The warning about spline stiffness that the
second pass raised for every series under 200 years is gone.
strip.rwl() reports what stripping cost. EPS over the whole span can
rise while the years that few trees reach get worse, because a series
that was left out may be one of the few there. When a series is left
out the function now prints the years that lost a tree and the years
in which EPS for the trees present is lower than before, and the
result carries the record as the attribute "dplR.strip": the series
left out, rbar and EPS before and after, and the trees and EPS for
each year. The series returned are the same as before. The help page
gains a caution, after Buras (2017), that a higher EPS is not evidence
of a better chronology for reconstructing climate.
chron.ars() builds the ARSTAN chronology from the residual
chronology, as its help page said it did and as Cook (1985) describes.
It used to add the pooled AR model to each individually whitened
series and average those, which left out the second whitening that the
residual chronology gets. Whitening each series with its own
coefficients can remove more persistence than the series share, and
the mean is then negatively autocorrelated. Where persistence is
mostly within trees the ARSTAN chronology inherited that: for wa082
detrended with a spline its lag-one autocorrelation was -0.16 against
0.19 for the standard chronology, and is now 0.10. Collections with
shared persistence change less (co021 0.15 to 0.22, ca533 0.47 to
0.52). The ars column changes for every data set, and it is now
missing in the same years as res, a few more years at the start
than before. std and res are unchanged.
The help page for chron() says that the residual chronology from
prewhiten = TRUE need not be white. Each series is whitened
separately, and where persistence is mostly within trees their mean
is negatively autocorrelated (about -0.2 at lag one for wa082). It
points to chron.ars(), which whitens the mean a second time. No
change to the function.
treeMean() returns class "rwi" when it is given indices, with the
record of how they were detrended. It used to return class "rwl"
whatever it was given, so chron(), rwi.stats() and the other
functions that want indices warned that tree means of indices were
ring widths. The help page had na.rm backwards: na.rm = FALSE
is the one that keeps only the years where a tree's samples overlap.
plot.crn() gains xlab and ylab, and passes ... on to plot().
Axis labels and graphical parameters such as col were silently
ignored before. One ylab labels the chronology panels and leaves
sample depth and running rbar alone; a vector labels each panel. A
chronology column with a name dplR does not know is now plotted as a
chronology instead of stopping with an error.
rwl.check()'s RWL_UNCHECKABLE says why a series could not be
checked. It always blamed gaps inside the series, but a series with a
long run of zeros is refused because its spline reaches zero, and
mi021's PLK18A, with no gaps and 1,691 filler zeros, was reported as
having gaps. The message now says which it was.
rwl.check() gains RWL_SERIES_SCALE, a warning for a series whose
median ring width is five times or more, or a fifth or less, of the
median of the other series (scale.ratio in rwl.check.control()).
A 0.001 mm Tucson series that has lost its -9999 stop marker is read
at 0.01 mm, ten times too wide, and nothing else catches it: the reader
has nothing to go on, the collection mean can stay plausible, and
crossdating is blind to a change of scale. Across the ITRDB ring-width
files a factor of 5 flags series in about 2% of files.
rwl.check() gains RWL_ZERO_RUN, a warning for a run of ten or more
zeros in a series (zero.run in rwl.check.control()). A run that
long is filler for wood that was not measured, not absent rings. It
used to be reported as a note, RWL_CONSECUTIVE_ZERO, which now covers
shorter runs only.
New kalman.ar() and kalman.spec() for evolutive spectra. They are
experimental: arguments, defaults and the objects returned may change
in a later release.
kalman.ar() fits an autoregressive model whose coefficients change
through time, by the Kalman filter and a smoother (Kitagawa and Gersch
1985), and kalman.spec() gives the AR spectrum of every year, with a
plot method that shows it as a surface in time and period. The order has no
default. The freedom of the coefficients, lambda, is estimated by
maximum likelihood unless set, and the fit is reported beside the
ordinary AR model with fixed coefficients. profile() on a fit gives
the log-likelihood over a range of lambda, with a plot method, to
show how well the data determine it and what a hand-set value costs.
With variance = "varying" the innovation variance changes through
time as well, estimated from the prediction errors, so that a change
in year-to-year scatter is not read as a change in the coefficients.
Suggested by Ed Cook; written from the published method, not ported
from his programs.
i.detrend() and i.detrend.series() now point to the iDetrend app
(https://github.com/OpenDendro/iDetrend), with a message once in a
session and a note on their help pages. The app does interactive
detrending better: curves can be adjusted and compared on each series,
and it writes R code that reproduces the indices. The functions are
not deprecated and work as before.
read.tucson() warns when a series ends without a stop marker (999 or
-9999), naming the series. The marker is the only statement of a
series' precision, and without it the reader assumes 0.01 mm. A
0.001 mm series that lost its -9999 therefore came back ten times too
large with nothing said. The warning says what was assumed, says when
other series in the file are at 0.001 mm, and is recorded as the event
NO_STOP_MARKER, which rwl.check() reports as RWL_NO_STOP_MARKER.
strict = TRUE refuses such a file. The values read are unchanged.
It fires on one of the 6,807 ITRDB ring-width files.
read.tucson() honours a stop marker written in an eleventh field,
past column 72, after a line of ten measurements. In 1.8.0 the marker
was dropped, so a 0.001 mm series ending that way was read as 0.01 mm
and every value came back ten times too large, with no warning.
read.tucson.legacy() was not affected. No file in the ITRDB has this
layout.
write.tucson() warns when it rounds values to prec. A Tucson file
holds whole numbers of 0.01 or 0.001 mm, and widths finer than that
were rounded without a word: 0.123 written at the default
prec = 0.01 came back from read.tucson() as 0.12. The warning
gives the number of values and series, the largest change and the
first series affected, and suggests prec = 0.001. The file written
is unchanged.
insert.ring() and delete.ring() work on a whole rwl object as
well as a single series: delete.ring(rwl, series = "50A", year = 1950)
returns the edited rwl. No measurement is dropped; the object gains or
loses years at its ends as needed, and gaps keep their place.
insert.ring() on a vector warns when fix.length = TRUE drops a
measured ring to keep the length.
fill.internal.NA() gains series, to fill only some series, and now
returns the class it was given: an rwl stays an rwl, with its read
record. It used to return a plain data.frame.
fill.internal.NA() gains fill = "Chron", which estimates a missing
ring from the other series in the collection rather than from the
series' own neighbouring rings. It is in the spirit of the gap filling
in ARSTAN, not a copy of it: the numbers will differ. A year is filled
only if min.series other series (default 3) were measured in it;
otherwise the call stops and names the years. An estimate below zero is
returned as zero with a warning. The filled values are estimates built
from the other trees, so they inflate statistics of agreement among
series and should not be used for crossdating. See the help page for
the assumptions and limits.
xdate.floater() dates a series that is longer than the master. It
used to stop with "'x' and 'y' must have the same length". Results for
a series shorter than the master are unchanged. It also stops with a
clear message when min.overlap is more than the years in the master.
xdate.floater() stops when series has NA inside it and says
where. It used to drop every NA, which closed the gap and put the
rings on either side out of step, so the series was dated wrongly, or
only in part, with no warning. Fill the gap first or date the parts
separately. NA before and after the measurements is still dropped.
plot() and print() work on the object xdate.floater() returns.
The methods were documented but not registered, so plot(fo) stopped
with an error and print(fo) listed the whole object.
seg.plot(), and so plot(rwl, plot.type = "seg"), draws an rwl
that holds one series. It used to stop with "wrong sign in 'by'
argument".
corr.rwl.seg() is about three times faster on large collections
(45 s to 13 s on the 597-series chin067), with the same results.
interseries.cor() is four to six times faster on large collections,
and rwl.report(), which computed it twice, now computes it once. On
the 597-series chin067, rwl.report() drops from 72 s to 8 s. The
results are unchanged.
A large release. This file lists what users need to know; the ChangeLog has the full record, with the reasons for each change.
read.tucson() is new, and no longer fills gaps inside a series with
zeros. A gap is now NA, so detrend() and other functions that
reject NA inside a series stop on such data until the gap is filled
on purpose. fill.internal.NA = 0 reproduces the old values. The old
reader is kept as read.tucson.legacy().detrend() and detrend.series() now default to method = "Spline".
The old default ran every method and returned a list, not indices.
Code that relied on it must now name the methods.detrend() with one method, rcs(), cms() and i.detrend() now
return class "rwi" (ring-width indices). This used to be a plain
data.frame, or class "rwl" for rcs(), cms() and i.detrend().
The values are unchanged.bai.in() and bai.out() now return class "bai" (basal area
increment). They used to return class "rwl". The areas are
unchanged.detrend(), rcs(), bai.in(), rwl.report(),
...) warn on class "rwi". Those that want indices (chron(),
rwi.stats(), sss(), ...) warn on class "rwl". Indices read
from a file come back as "rwl"; label them with as.rwi().detrend(), rcs(), cms(), i.detrend(), interseries.cor(),
corr.rwl.seg() and subset() drop series with no values, with a
message naming them. Before, they stopped or kept a column of NA.detrend(difference = TRUE), are no longer sign-flipped when
normalized. That flipping gave wrong rbar, EPS and interseries
correlations (GitHub issue 22). Results for ring widths and ratio
indices are unchanged.detrend.series(difference = TRUE) now uses fitted curves that go
below zero as they are, instead of swapping them for a simpler fit.caps() stops on a series containing NA. Before, it returned all
NA.chron(x) and chron.ars(x) are
now rwi, and sgc(x) is now rwl. Calls that named the argument
x = must change.nyrs and ar.order.max were added to the crossdating and
statistics functions, next to n and prewhiten. Arguments after
them passed by position now shift. An old call that shifts stops with
an error; none silently changes meaning.csv2rwl() is deprecated in favour of read.sheet(), and may refuse
files it used to accept wrongly.rwl.check() looks for defects in a collection (duplicate or empty
series, gaps, zeros, units errors, dating errors) and returns them as
rows with stable IDs, for one file or a whole archive.xdate.report() and write.xdate.report() produce a COFECHA-style
crossdating report, with A and B flags and the lag and correlation
gain for each flagged segment.corr.rwl.seg() gains lag.max. It tests each segment at shifted
dates and returns the best lag. Negative lags mean missing rings, as
in COFECHA.read.sheet() and write.sheet() read and write spreadsheet-shaped
files (years down, series across, or long format). read.rwl() uses
them for csv.read.tucson() records what it saw in a provenance record that
travels with the data."rwi" class has summary(), which gives collection statistics
and flags series that don't fit, and plot(x, plot.type = "image"),
which shows trend left in by detrending."bai" class has as.bai() and methods for [, subset(),
time(), window(), summary() and plot(). chron(),
detrend() and the crossdating functions take it without a warning;
functions that want ring widths warn.window() methods for rwl, rwi, bai and crn objects:
window(x, 1800, 1900).[ and subset() methods for rwl and rwi objects keep the class
and the records, and trim years no series covers.