In today’s post, I’m going to critically examine another widely used tree ring chronology: the Icefields (Alberta) MXD RCS chronology of Luckman and Wilson (2005 pdf), used most recently in Wilson et al 2016.
I’ll show that the RCS technique used in the LW2005 MXD chronology eliminated high medieval values as a tautology of their method, not as a property of the data and that the Icefields data provides equal (or greater) justification for MXD RCS chronologies with elevated medieval values. The measure of potential difference is previewed in Figure 1 below, which compares the reported LW2005 chronology (top panel) to an MXD chronology with elevated medieval values, calculated using equally (or more plausible) RCS methodology, with several other similar variations shown in the post.
Figure 1. Top – MXD RCS chronology from Luckman and Wilson 2005 Figure 2; bottom – MXD RCS chronology calculated under alternative assumptions – see Figure 3 below.
I will use random effects statistical techniques both to analyse the data and to place prior analysis in a more formal statistical context. Because LW2005 coauthor Rob Wilson stands alone for civility in the paleoclimate world and because the present post is critical of past analysis, in some ways, I would have preferred to use another example. However, on the other hand, it is also possible that Rob will recognize that the techniques applied in the present post – techniques that are unfamiliar to dendros – can yield fresh and important insight into the data and will look at the Icefields data a little differently.
Although the article in consideration was published more than a decade ago, the analysis in today’s article was impossible until relatively recently, because coauthor Luckman withheld the relevant data for over a decade.
Background
In 1997, Brian Luckman (Luckman 1997) published a tree ring chronology from the Athabasca Glacier, Alberta area that was used in several contemporary multiproxy studies (Jones et al 1998; Crowley and Lowery 2000).
In the late 1990s and early 2000s, the data was updated and expanded. In 2005, Luckman and a former student well known to CA readers (Rob Wilson) combined data from the Athabasca Glacier location with data from near Peyto and Robson glaciers to form “regional chronologies” for both MXD (density) and RW (ring widths), settling on “RCS” and “standard” methods, respectively, though, as we will see below, the eventual “RCS” method was idiosyncratic. The MXD and RW chronologies were then used in a regression to calculate a temperature reconstruction, which, in turn, has been used in multiple subsequent multiproxy and/or multi-site studies, including the recent Wilson et al 2016 multi-site temperature reconstruction. reported in Luckman and Wilson (2005 ).
As is far too common in paleoclimate, Luckman did not archive his data when the article was published. Luckman also turned down/ignored requests for measurement data that I made many years ago. In May 2014, long after the original articles had faded from view, Luckman finally archived the underlying measurement data. (By coincidence, this was almost exactly the same time as the decades overdue archiving by Jacoby and D’Arrigo.) I didn’t notice this long overdue archiving at the time. However, during my re-examination of series used in Wilson et al 2016, I noticed that the long refused data for the Icefields sites had finally been archived, making the present assessment possible for the first time.
Unsurprisingly, there are some surprising results – if I may put it that way.
Today’s post will be limited to the MXD data; I hope to get to the RW data on another occasion. I’m looking first at the MXD data since CA readers will recall some recent assertions by specialists favoring MXD (density) data over RW (ring-width) data.
The MXD Data
Because of the long delay between original publication (2005) and archiving (2014), I took particular care to try to ensure that the data in my present collation matches (at least as closely as possible) the data that was used in the original article. Figure 1 of Luckman and Wilson shows that data for the regional chronology came from three sites: Athabasca Glacier, Peyto Glacier and Robson Glacier.
Figure 2 of Luckman and Wilson showed core counts by year for four groups: Icefields living and “snags” (subfossil), Peyto living and snags from Peyto and Robson snags.
I was able to closely replicate LW2005 Figure 2 (See bottom panel at right) by collating the following six MXD datasets at the NOAA archive:
- three datasets from the Athabasca Glacier area (living: cana096x, sampled in 1983, and cana430x, sampled in 1994 and archived in 2014; subfossil: cana431x, archived in 2014.
- two datasets from Peyto Glacier area (living: cana097x, sampled in 1983; subfossil: cana445x, archived in 2014.
- one dataset from Robson Glacier area (subfossil: cana448x, archived in 2014).
There are two Athabasca Glacier MXD datasets at NOAA that were not used (cana170x, cana171x – both from Schweingruber.) My collation of the MXD data from the six sites into a “tidy” file (tidy, as defined in R usage) is at http://www.climateaudit.info/data/proxy/tree/wilson/2005/collated_icefields_mxd.csv.
RCS Chronology Variations
“Standard” (as opposed to “RCS”) tree ring chronologies fit a growth curve to each individual tree and then take an average of ratios (rather than residuals, though the effect is similar) for each year to form a “chronology”. As I’ve discussed from time to time, this is formally equivalent to the assignment of random effect to each tree. Since the trees have different life spans, this technique necessarily transfers variance from annual effects (the chronology) to inter-tree variance and effectively eliminates centennial variability – a phenomenon familiar to dendros (the “segment-length curse”) though they don’t explain it in the statistical terminology preferred at CA.
“RCS” methods are homemade attempts by dendros to mitigate this difficult problem. One common technique is to use single regional curve to “standardize” cores for juvenile growth , as, for example, in the Yamal regional chronology of Briffa 2000 or the Gulf of Alaska chronology of Wiles et al 2014. Luckman and Wilson 2005 briefly discussed use of a “single regional curve”, but dismissed it summarily on the grounds that it might “possibly introduc[e] bias in the earliest part of the reconstruction”:
Using a single, regional curve would inflate the final indexed values for the Peyto-Robson data, possibly introducing bias in the earliest part of the reconstruction.
I cannot help but comment that one seldom if ever hears corresponding concerns about inflation of final indexed values when the bias pertains to the latest part of the reconstruction. The results of an MXD RCS calculation using a single regional curve are illustrated in the bottom panel in Figure 1 above: it has a very “warm” medieval period. Perhaps this is the “bias” that concerned the authors of LW2005.
To implement a single regional curve RCS using a well-known statistical model, I fitted the data using the nls function of the nlme package (and then taking annual averages of the residuals). I used a negative exponential, as a Hugershoff function is not required here – in a Hugershoff fit, the coefficient is close to zero. The single curve model is:
fmnls=nls(x~ A + B* exp(-C*age),data=tree,start=c(A=50,B=30,C=.02))
To examine potential random effects from the three different sites, I used the lmer function of the lme4 package to concurrently estimate crossed random effects for (1) site and (2) year (as discussed in a prior post here ^). As in my prior discussion, I simplified the calculation by keeping the negative exponential coefficient C fixed (using the coefficient from the nls fit).
nm= lmer(x~ exp(-0.0030179 * age)+(1|site)+(1|yearf),data=tree)
While I have not exhaustively studied the matter, the chronology of random effects do not appear to be materially impacted by slight changes in this coefficient.
I also did a model in which I included one more crossed factor distinguishing between subfossil/living – as well as for site.
nm1= lmer(x~ exp(-0.0030179 * age)+(1|site)+(1|sf)+ (1|yearf),data=tree)
I am more than a little troubled about introducing a random effect for a subfossil/living distinction without evidence of a physical difference, because the introduction of such a factor will necessarily transfer variance from the chronology to the factor and, in the process, compromise the chronology – an issue that I will discuss further below.
The next figure shows the RCS chronologies resulting from each of these three statistical models: top – nls; middle – lmer with a random effects for year and site; bottom – lmer with random effects for year, site and subfossil/living. In each of the above three cases, there is a pronounced MWP period, most pronounced in the nls case (see discussion in caption), but evident in all three. Each has a very different medieval-modern relationship as compared with the reported LW2005 chronology.
Figure 3. RCS chronology variations: top – nls fit; middle – lmer fit with random factors for year and site; bottom – lmer fit with random factors for year, site and living/subfossil. As we previously observed with the Polar Urals MXD chronologies, the lmer chronologies damp chronology values when measurements are sparse – observable here in the earliest values. When measurements are sparse, the lmer method keeps some of the variance with the residuals.
The LW2005 RCS Chronology
So how did Luckman and Wilson arrive at their RCS chronology without an MWP?
Rather than calculating a random effect for each site or crossed random effects for site and living/subfossil (as the two examples above), they divided the data into three somewhat idiosyncratic groups:
- living trees (combining Athabasca and Peyto);
- Icefields/Athabasca snags (subfossil); and,
- a combination of Peyto and Robson snags.
As an experiment, I calculated a chronology using lmer with crossed factors for (1) year, and (2) the three LW2005 groups. This did indeed yield a chronology that no longer had the high medieval values of the other RCS choices shown above.
Figure 3. Top – MXD RCS chronology from Luckman and Wilson 2005 Figure 2; bottom- lmer chronology of annual effects, using LW2005 grouping as a random effect. As in other examples, the lmer technique damps early values when measurements are sparse. The model used is:
nm3= lmer(x~ exp(-0.0030179 * age)+(1|alias)+ (1|yearf),data=tree).
However, there is a serious hidden cost to the grouping selected in LW2005.
In a grouped RCS calculation, each of the three groups is effectively centered at the same mean, even though their temporal coverage is entirely different. The average dates of the Peyto/Robson snags, Icefields snags and living trees are 1109, 1484 and 1855 respectively. The average date of the Peyto and Robson subfossil material is in the heart of the medieval period, while the average date of the living tree datasets is in the very cold 19th century. Nonetheless, the LW2005 grouping centers all three groups at the same level, as can be seen in the table below showing the average ring width for each group, together with the random effect of the above statistical model: after subtraction of the random effect for each group, the group values become virtually identical:
By choosing to separately standardize living and medieval data (Peyto/Robson subfossil), the possibility of drawing conclusions on medieval-modern data is rendered impossible, in effect defeating the purpose of the RCS calculation.
I think that the problem is neatly illustrated by the next diagram, which shows average annual MXD values from the two Peyto datasets (living and subfossil). In a single regional curve – or even when a single factor is assigned to the Peyto site – the high medieval MXD values yield high values in a chronology. However, if Peyto subfossil/snags are treated as a different site than Peyto living (as, effectively in LW2005) the high MXD values are allocated as a random effect of the site and cancelled out of the chronology.
Figure 4. Annual average MXD values from Peyto glacier data, combining data on living trees and snags.
LW2005 purported to justify a technique that cancelled the high medieval values on the basis that this technique permitted the “better assessment of potential bias from the RCS method”:
To address the potential occurrence of differing ‘populations’ within the MXD data, the MXD series were divided into three groups, to allow better assessment of potential bias from the RCS method (Table 4).
However, I fail to see how the LW2005 permutation allowed a “better assessment of potential bias”. Even if this permutation proved better in some sense than other potential permutations (including the ones presented in this post), LW2005 merely presented this result as fait accompli, rather than as a measured recommendation following analysis of alternatives. To actually carry out such an “assessment”, one would need to examine results from other plausible groupings, such as those presented above – all of which had pronounced MWPs – and explain precisely why the grouping that happened to eliminate the MWP was “better” on a criterion more meaningful than a Calvinball elimination of high MWP values.
Conclusion
The RCS technique used in the LW2005 MXD chronology eliminated high medieval values as a tautology of their method, not as a property of the data. The Icefields data provides equal (or greater) justification for MXD RCS chronologies with elevated medieval values.
As I work through these examples, I am increasingly convinced that a statistical model using crossed random effects for site and for year provides a relatively rational basis for accommodating site inhomogeneity without sacrificing “low frequency” variability. It certainly seems more effective in each of the recent examples. This does not mean that the chronology of annual effects is necessarily a temperature “proxy” – only that such a chronology is a better estimate of annual effects (which may be due to temperature, precipitation, insects, fire etc.)
However, the analysis of today’s post was impossible until recently because coauthor Luckman withheld relevant measurement data for over a decade. One of the retorts to the present analysis will undoubtedly be that I’m wasting my time dissecting a decade-old paper, because the field has “moved on” (all too often to still unarchived data (e.g. most of the component series of Wilson et al 2016) which, if the past is any guide, may not be available for another decade.
But while new datasets have arisen during this period and while “RCS” chronologies have become fashionable among dendroclimatologists, dendroclimatologists have limited understanding of the statistical issues involved in the construction of a chronology, and, in the absence of a rational statistical theory for “RCS” chronologies, the construction of such chronologies all too often seems to become prone to a peculiar form of Calvinball in which arbitrary decisions invariably seem to inflate modern values relative to medieval values (as appears to have taken place in the present and other cases.)
References:
Luckman et al 1997. Tree-ring based reconstruction of summer temperatures at the Columbia Icefield, Alberta, Canada, AD 1073-1983. The Holocene.
Luckman and Wilson 2005. Summer temperatures in the Canadian Rockies during the last millennium: a revised record. Climate Dynamics pdf.
Wilson et al 2016. Last millennium northern hemisphere summer temperatures from tree rings: Part I: The long term context. QSR. pdf