ABSTRACT (original article): Dendrochronology, the study of annual rings formed by trees and woody plants, has important applications in research of climate and environmental phenomena of the past. Since its inception in the late 19th century, dendrochronology has not had a way to quantify uncertainty about the years assigned to each ring (dating). There are, however, many woody species and sites where it is difficult or impossible to delimit annual ring boundaries and verify them with crossdating, especially in the lowland tropics. Rather than ignoring dating uncertainty or discarding such samples as useless, we present for the first time a probabilistic approach to assign expected ages with a confidence interval. It is proven that the cumulative age in a tree-ring time series advances by an amount equal to the probability that a putative growth boundary is truly annual. Confidence curves for the tree stem radius as a function of uncertain ages are determined. A sensitivity analysis shows the effect of uncertainty of the probability that a recognizable boundary is annual, as well as of the number of expected missing boundaries. Furthermore, we derive a probabilistic version of the mean sensitivity of a dendrochronological time series, which quantifies a tree’s sensitivity to environmental variation over time, as well as probabilistic versions of the autocorrelation and process standard deviation. A computer code in Mathematica is provided, with sample input files, as supporting information. Further research is necessary to analyze frequency patterns of false and missing boundaries for different species and sites. CITATION (original article): Ricker M, Gutiérrez-García G, Juárez-Guerrero D, Evans MEK (2020) Statistical age determination of tree rings. PLoS ONE 15(9): e0239052. https://doi.org/10.1371/journal.pone.0239052

Introduction:

To provide the big picture of what this is about, it is useful to read first the article’s abstract: Dendrochronology, the study of annual rings formed by trees and woody plants, has important applications in research of climate and environmental phenomena of the past. Since its inception in the late 19th century, dendrochronology has not had a way to quantify uncertainty about the years assigned to each ring (dating). There are, however, many woody species and sites where it is difficult or impossible to delimit annual ring boundaries and verify them with crossdating, especially in the lowland tropics. Rather than ignoring dating uncertainty or discarding such samples as useless, we present for the first time a probabilistic approach to assign expected ages with a confidence interval. It is proven that the cumulative age in a tree-ring time series advances by an amount equal to the probability that a putative growth boundary is truly annual. Confidence curves for the tree stem radius as a function of uncertain ages are determined. A sensitivity analysis shows the effect of uncertainty of the probability that a recognizable boundary is annual, as well as of the number of expected missing boundaries. Furthermore, we derive a probabilistic version of the mean sensitivity of a dendrochronological time series, which quantifies a tree’s sensitivity to environmental variation over time, as well as probabilistic versions of the autocorrelation and process standard deviation. A computer code in Mathematica is provided, with sample input files, as supporting information. Further research is necessary to analyze frequency patterns of false and missing boundaries for different species and sites.

The input:

The notebook imports an Excel file (attached below) of type “csv” (“comma delimited”); this type is automatically a single-sheet file, and the Mathematica import of this format works better for large files than the “xlsx” type. The name of our input file is “S2 Alchornea Input.csv”, where “S2” refers to the second supplementary file, and “Alchornea” is the genus name of the tree, whose data was analyzed. The input file is imported from the directory of the Mathematica notebook:
In[]:=
SetDirectory[NotebookDirectory[]];​​inputFile="S2 Alchornea Input.csv";​​originalData=Import[inputFile];​​TableForm[originalData]
Out[]//TableForm=
Boundary number
Increment (cm)
BoundaRecognizabley type
Probability
0
0.124
Annual
1
1
0.099
Recognizable
0.83
2
0.256
Recognizable
0.83
3
0.129
Recognizable
0.83
4
0.167
Recognizable
0.83
5
0.148
Recognizable
0.83
6
0.22
Recognizable
0.83
7
0.3
Recognizable
0.83
8
0.172
Recognizable
0.83
9
0.143
Recognizable
0.83
10
0.119
Recognizable
0.83
11
0.282
Recognizable
0.83
12
0.092
Recognizable
0.83
13
0.196
Recognizable
0.83
14
0.131
Recognizable
0.83
15
0.179
Recognizable
0.83
16
0.136
Recognizable
0.83
17
0.444
Recognizable
0.83
18
0.034
Recognizable
0.83
19
0.482
Recognizable
0.83
20
0.146
Recognizable
0.83
21
0.458
Recognizable
0.83
22
0.067
Recognizable
0.83
23
0.271
Recognizable
0.83
24
0.217
Recognizable
0.83
25
0.16
Recognizable
0.83
26
0.219
Recognizable
0.83
27
0.173
Recognizable
0.83
28
0.133
Recognizable
0.83
29
0.149
Recognizable
0.83
30
0.213
Recognizable
0.83
31
0.318
Recognizable
0.83
32
0.154
Recognizable
0.83
33
0.194
Recognizable
0.83
34
0.103
Recognizable
0.83
35
0.241
Recognizable
0.83
36
0.33
Recognizable
0.83
37
0.094
Recognizable
0.83
38
0.384
Recognizable
0.83
39
0.111
Recognizable
0.83
40
0.239
Recognizable
0.83
41
0.196
Recognizable
0.83
42
0.267
Recognizable
0.83
43
0.092
Recognizable
0.83
44
0.181
Recognizable
0.83
45
0.115
Recognizable
0.83
46
0.136
Recognizable
0.83
47
0.256
Recognizable
0.83
48
0
Recognizable
0.83
As shown, the file consists of four columns and any number of rows. The first line contains the column titles; the exact column titles don’t matter, because they are only for orientation, and are not used further. Of the four columns, only the second (radius increments) and fourth (probabilities) are used by the algorithm, while the other two columns contain auxiliary information for labeling the output. The four columns contain the following data:
​
1) The boundary number is typically an integer, though for missing boundaries, we use for example 10, 10.5, 11, where 10.5 refers to a missing boundary after boundary 10.
​
2) The radius increment in centimeters (cm), though the unit is only used for labeling. The increments should be positive real number.
​
3) The boundary type indicates annual (probability = 1), recognizable (probability < 1), or missing boundaries. The column is only for orientation, and is not used further.
​
4) The probabilities should be positive real numbers larger than 0, and not larger than 1. For a recognizable boundary, the number represents the probability that a recognizable boundary is an annual, regular (and not a false) boundary. For a missing boundary, the number represents the probability that a missing annual boundary is expected at a specific location in a core sample. All probabilities after the base boundary are 0.83 here, but that must not always be so.
​
The first row is the base boundary, numbered as boundary 0, and having a probability of 1 (recognizable and regular). It contains the increment of the first tree ring. The probabilities of the subsequent boundaries may either vary or all be the same. The last boundary of the sample is distinct in that it has no increment (0 cm).

Calculate the probabilities for the different possible ages each tree-ring boundary could have:

As the first step, the column titles are deleted, and the data is sorted according to the chronology of the tree-ring boundaries (in our input file, the data is already ordered). Then the input probabilities are extracted in a list and rationalized. Rationalization causes the calculations to be carried out with exact numbers in Mathematica (0.83 becomes 83/100). This is recommended for two reasons: First, probabilities for some ages can become very small, and the sum of all probabilities for all possible ages of a given boundary may not be exactly 1 anymore. Second, interpolation of the discrete points in the cumulative probability function, for calculating confidence limits, may not work anymore with approximate numbers, for numerical reasons.
In[]:=
data=SortBy[Drop[originalData,1],First];​​probabilities=Rationalize[Transpose[data][[4]]];
The following code represents the key part to calculate iteratively the expected probabilities of different possible ages that each tree-ring boundary may have, starting at the base boundary (number 0). The base boundary is defined to have an age of 0 years and a probability of being an annual (regular) boundary of 1. In each iteration, an additional boundary is considered, and with it the number of possible ages augments, unless the input probability is 1. The boundary that is farthest away from the base boundary thus has the most possible ages, each with a certain probability. The sum of the probabilities of all possible ages for a given tree-ring boundary is always 1. In each iteration, the variables “calculatedProbsPerAge” is calculated for the corresponding tree-ring boundary. We will go through the first three iterations to see how the code works.
In[]:=
calculatedProbsPerAge={{-1},{1}};(*Age,probability*)​​ageResults=​​Table[​​If[probabilities[[i]]≠1,​​newAges=Append[calculatedProbsPerAge[[1]],calculatedProbsPerAge[[1,-1]]+1];newcalculatedProbs=Append[calculatedProbsPerAge[[2]]*(1-probabilities[[i]]),0]+Prepend[calculatedProbsPerAge[[2]]*probabilities[[i]],0],​​(*probabilities[[i]]==1*)​​newAges=calculatedProbsPerAge[[1]]+1;​​newcalculatedProbs=calculatedProbsPerAge[[2]]];(*endofIf*)​​calculatedProbsPerAge={newAges,newcalculatedProbs},​​{i,1,Length[data]}];
The input probability of the base boundary is 1 (for i = 1):
In[]:=
i=1;​​data[[i]]​​probabilities[[i]]
Out[]=
{0,0.124,Annual,1}
Out[]=
1
Given the If-clause, the first iteration reduces to the following:
In[]:=
calculatedProbsPerAge={{-1},{1}};(*Age,probability*)​​newAges=calculatedProbsPerAge[[1]]+1​​newcalculatedProbs=calculatedProbsPerAge[[2]]​​calculatedProbsPerAge={newAges,newcalculatedProbs}
Out[]=
{0}
Out[]=
{1}
Out[]=
{{0},{1}}
The starting values before the first iteration of “calculatedProbsPerAge” is set to a (meaningless) age of -1 year and a probability of 1. The addition of 1 in the first iteration results in an age of 0 years and 1 for the probability to represent an annual boundary, for the base boundary. The input data for the second iteration (i = 2) is:
In[]:=
i=2;​​data[[i]]​​probabilities[[i]]
Out[]=
{1,0.099,Recognizable,0.83}
Out[]=
83
100
With the input probability not being 1, the first option of the If-clause applies:
In[]:=
calculatedProbsPerAge={{0},{1}};(*fromthelastiteration*)​​i=2;​​​​newAges=Append[calculatedProbsPerAge[[1]],calculatedProbsPerAge[[1,-1]]+1]​​newcalculatedProbs=Append[calculatedProbsPerAge[[2]]*(1-probabilities[[i]]),0]+Prepend[calculatedProbsPerAge[[2]]*probabilities[[i]],0]​​calculatedProbsPerAge={newAges,newcalculatedProbs}
Out[]=
{0,1}
Out[]=

17
100
,
83
100

Out[]=
{0,1},
17
100
,
83
100

The second boundary thus is either 0 years old (it is not annual), with a probability of 0.17, or it is 1 year old with a probability of 0.83. The third iteration (i = 3) follows:
In[]:=
i=3;​​data[[i]]​​probabilities[[i]]
Out[]=
{2,0.256,Recognizable,0.83}
Out[]=
83
100
In[]:=
calculatedProbsPerAge={{0,1},{17/100,83/100}};(*fromthelastiteration*)​​i=3;​​​​newAges=Append[calculatedProbsPerAge[[1]],calculatedProbsPerAge[[1,-1]]+1]​​newcalculatedProbs=Append[calculatedProbsPerAge[[2]]*(1-probabilities[[i]]),0]+Prepend[calculatedProbsPerAge[[2]]*probabilities[[i]],0]​​calculatedProbsPerAge={newAges,newcalculatedProbs}
If the input probability is smaller than 1, then there is a certain probability that the additional tree-ring boundary is false, and thus no age would be added. However, if the additional tree-ring boundary is annual, then an additional age is attached in the list of “newAges”. As for the corresponding probabilities “newcalculatedProbs”, a 0 is appended at the end of the list of the previous probabilities, and all previous probabilities are multiplied with 1 minus the new input probability. The result is a list of the probabilities in the case that the boundary is false. A second list is generated by prepending a 0 to the list of previous probabilities, and all previous probabilities are multiplied with the new input probability. The result is a list of probabilities in the case that boundary is annual. The two lists are added to provide the final probabilities for each possible age, as shown below:
The list from "Transpose[ageResults[[-1]]]" represents a discrete probability density function that is both left and right-censured, shown in the following graph:

Calculate the additional statistical information, including expected ages and confidence intervals:

The result from the last section is a list of 49 sublists (corresponding to 49 boundaries for 48 tree rings), each consisting of 2 sub-sublists of variable length. One sub-sublist contains the possible ages and the other the corresponding probabilities.
​
This is easier to read in table form, and with probabilities as real numbers:
​
The fourth group corresponds to the third boundary after the base boundary. With a probability of 0.004913 all three boundaries after the base boundary are false, and the age of the wood at that part is still somewhere between 0 and 1 years. With a probability of 0.071961, the boundary is the first annual boundary after the base boundary (its age is 1 year), with 0.351339 the second (its age is 2 years), and with 0.571787 the third (its age is 3 years). The first parameters that are taken from that data are the minimum and maximum possible ages a boundary may have:
Of more interest is the expected age of each boundary. It is calculated for each boundary by multiplying each possible age with the corresponding probability, and taking the sum of the products (exact numbers are again converted to real numbers with “N”):
The formula for calculating the variances is slightly more complicated. The square of the age is taken before multiplication with the corresponding probability, the sum is taken again, and the square of the expected age is rested. The variances are presented together with the number of data points, so that subsequently one could also calculate the corresponding standard errors.
The standard deviation of the possible ages around the expected age is the square root of the variance:
To calculate confidence intervals, the inverse of the continuous Cumulative Distribution Function is needed. The points of ages at different probabilities are interpolated with straight lines, resulting in a piecewise linear (“PL”) Function:
The function values (converted to real numbers) for graphing the inverse Cumulative Distribution Function are:
The graph shows that on the left, the age is extremely sensitive to a change in cumulative probability, then there is a wide range where age is only moderately affected by a change in cumulative probability, and then there is a final range where age is again highly sensitive to a change in cumulative probability.
​
The previous explanations were all for the last tree-ring boundary. The procedure to calculate the continuous Cumulative Distribution Function and determine 95%-confidence limits has to be carried out for each boundary after the base boundary:
Following, the confidence intervals and age ranges are graphed as a function of boundary number, to be able to get an overall idea of their development. The expected age is shown with black circles, its 95% confidence intervals with red points, and the possible age range as black lines:
Next, three parameters are calculated to characterize the discrete confidence curves that are shown with red points in the previous graph. The first parameter is the mean width of all confidence intervals (years), without taking into account the base boundary that has always a zero confidence interval:
The confidence intervals are asymmetric around the expected age. The second parameter calculates the difference between the expected age and lower confidence limit, divided by the difference between upper confidence limit and expected age. Then the mean of all calculated ratios is taken:
This calculated ratio (1.95385) as an exact number is remarkably long: it has 2,357 integer digits in both the numerator and the denominator (without the possibility of simplification):
The confidence intervals also widen with more distance from the base boundary. The third parameter calculates how much (in years) the confidence intervals widen on average with each additional boundary:
A fundamental graph in forestry is a growth curve of trunk radius as a function of tree age. It can now be drawn, together with confidence intervals and overall possible-age range (function points are interpolated linearly). The radius increments are part of the original input. In contrast to the calculated probabilities (for different possible ages of tree-ring boundaries), where extreme differences of numbers resulted, and it was of great advantage to use Mathematica’s ability to calculate with exact numbers, there is no need to use exact numbers for the radius increments.
​
Note that the following graph of radius as a function of age implies inverted axes, because we originally calculated expected age and 95% confidence intervals (red lines) for given radius (as implied from boundary numbers). In contrast to the usual age-radius curves, the curves are delimited horizontally at the maximum measured radius, not vertically at the maximum measured age. Furthermore, there can be vertical lines, connecting several radii that are associated with the same age. The outer black lines show the possible age range. Points in each case are linearly interpolated.

Calculate the mean sensitivity of the tree-ring series:

The mean sensitivity of a tree-ring series was introduced by Andrew Douglass in 1936 (Climatic cycles and tree growth volume III: a study of cycles. Washington, D.C.: Carnegie Institution of Washington, publication No. 289). The formula is based on the tree-ring widths over the trunk radius that are part of the original input:
First, Douglass’ original formula is used (see (4) in Ricker et al. 2020):
Next, the probabilistic formulas from Ricker et al. (2020) is used, including the mean sensitivity’s variance, standard error, and approximate 95% confidence intervals:
Finally, the alternative indicators for mean sensitivity are calculated in their deterministic and probabilistic versions (see Table 6 in Ricker et al. 2020):

The overall output:

The overall output, presented below and in two Excel output files, starts with the four graphs that were already shown above. Next, two tables are shows. The first table presents five columns: the boundary number, the boundary type, and the assigned probability (all from the input file), the radius calculated cumulatively from the input increments, and the resulting range of possible ages for each boundary number. The second table repeats the boundary number, followed by the expected age, its standard deviation, the corresponding number of ages from which the two previous parameters were calculated, and the lower and upper 95% confidence limits. If the lower 95% confidence limit does not exist at early ages, then the lower age limit is given. All parameters are given from the base boundary to the last boundary. The tables are followed by three parameters that characterize the 95% confidence intervals, and were described above. Exact numbers are generally converted to real numbers.
​
Finally, the classical mean sensitivity is shown (with the number of widths for its calculation), followed by the probabilistic mean sensitivity (with its standard error, and confidence intervals), and the process standard deviation in classical and probabilistic form (with the corresponding first-order autocorrelations and time series’ variance).

In addition, all numerical values of the output (used also for drawing the graphs) are stored in two Excel files in the current directory. The two files are called “Probabilistic ages output.xlsx” (with the data from the two tables) and “Additional parameters.xlsx” (with all additional data). Note that files with these same name are overwritten without any warning (unless they are open). Also note that more decimals of the provided real numbers are given in these output files than in the above presentation.

CITE THIS NOTEBOOK

Statistical age determination of tree rings​
by Martin Ricker​
Wolfram Community, STAFF PICKS, January 17, 2021
​https://community.wolfram.com/t/17454