BI Change Point Detection
Hung Nguyen
Introduction
This document reproduces the paper “Correcting the heartwood-sapwood transition in blue intensity measurements with change point detection methods” by Nguyen et al. (2025, Dendrochronologia).
If you don’t have the following packages, please install them first.
install.packages(
c('data.table',
'changepoint',
'dplR',
'metR',
'ggplot2',
'cowplot',
'ggtext',
'ggprism',
'patchwork',
'magick',
'maps',
'remotes'))
remotes::install_github('eliocamp/tagger')We’re now ready to run the background setup for this project.
knitr::opts_chunk$set(echo = TRUE)
source('R/init.R') # Useful functions
source('R/read_tucson2.R') # A faster read_tucson()
theme_set(theme_prism(base_size = 12, base_line_size = 0.2, base_fontface = 'plain'))
illiniBlue <- '#13294B' # University of Illinois colors
illiniOrange <- '#FF5F05' # University of Illinois colorsData
# Rainfall data Chiang Mai 1951
chm <- fread('data/thai-precip-full.csv', key = 'ID')['CHM'][year(date) == 1951]
chm[, cumuP := cumsum(pre)] # Cumulative precipitation
chm[, doy := yday(date)] # Day of the year
# BI data
mckLW <- read.tucson2('data/MCKlwBI.rwl', verbose = FALSE) |> rwl_to_dt(value.name = 'lwbi')
mckEW <- read.tucson2('data/MCKewBI.rwl', verbose = FALSE) |> rwl_to_dt(value.name = 'ewbi')
# Fix some name errors produced by CooRecorder
mckLW[nchar(core) > 6, core := substr(core, 1, 5)]
mckLW[core %in% c('MCK20L', 'MCK24L', 'MCK36L'), core := substr(core, 1, 5)]
mckEW[nchar(core) > 6, core := substr(core, 1, 5)]
mckEW[core %in% c('MCK20E', 'MCK24E', 'MCK36E'), core := substr(core, 1, 5)]
mck <- merge(mckLW, mckEW, by = c('core', 'year'))
# Some cores have "MCK" prefix while the other cores have "MC" prefix.
# Standardize so that they appear in the right alphabetical order in plots
mck[substr(core, 3, 3) == 'K',
core := paste0('MC', substr(core, 4, nchar(core)))]
mck[core == 'MC200', core := 'MC20A']
mck[core == 'MC201', core := 'MC20B']
mck[core == 'MC230', core := 'MC23A']
mck[core == 'MC231', core := 'MC23B']
mck[core == 'MC245', core := 'MC24A']
mck[core == 'MC290', core := 'MC29A']
mck[core == 'MC291', core := 'MC29B']
mck[core == 'MC345', core := 'MC34A']
mck[core == 'MC390', core := 'MC39A']
mck[core == 'MC392', core := 'MC39B']
setkey(mck, core)
# cumulative BI
mck[, ':='(clw = cumsum(lwbi),
cew = cumsum(ewbi)),
by = core] Example HW-SW having different effects on EW and LW (Figure 1).
mck23B <- mck['MC23B']
mck23B[, dbi := lwbi - ewbi]
mck23B[, ':='(lwlp = pass.filt(lwbi, 100),
ewlp = pass.filt(ewbi, 100),
dblp = pass.filt(dbi, 100))]
p1 <- ggplot(mck23B) +
geom_line(aes(year, ewbi, color = 'EWBI'), linewidth = 0.4) +
# geom_line(aes(year, ewlp, color = 'EWBI'), linetype = 2) +
geom_line(aes(year, lwbi, color = 'LWBI'), linewidth = 0.4) +
# geom_line(aes(year, lwlp, color = 'LWBI'), linetype = 2) +
geom_line(aes(year, dbi, color = 'DBI'), linewidth = 0.4) +
# geom_line(aes(year, dblp, color = 'DBI'), linetype = 2) +
scale_color_manual(
values = c(EWBI = 'steelblue',
LWBI = illiniBlue,
DBI = 'darkorange'),
breaks = c('LWBI', 'EWBI', 'DBI')) +
scale_x_continuous(
guide = guide_prism_minor(),
breaks = seq(1750, 2000, 50),
minor_breaks = seq(1710, 2010, 10)) +
labs(x = NULL, y = 'BI [-]') +
theme(legend.position = 'top')
p2 <- image_ggplot(image_read('figures/MCK231.jpg'))
p1 / p2 +
plot_layout(heights = c(5, 1))
Change point detection
Code for the 3 methods
# Two-phase linear regression
# Calculate sums of squares on the left and right of a change point
ssq_left_right <- function(x, cp, required.direction = 'down') {
idx <- seq_along(x)
idx1 <- which(idx <= cp)
idx2 <- which(idx > cp)
t1 <- matrix(c(rep(1, length(idx1)), idx1), ncol = 2)
t2 <- matrix(c(rep(1, length(idx2)), idx2), ncol = 2)
left <- .lm.fit(t1, x[idx1])
right <- .lm.fit(t2, x[idx2])
if ((required.direction == 'down' & left$coefficients[2] > right$coefficients[2]) |
(required.direction == 'up' & left$coefficients[2] < right$coefficients[2]) |
(required.direction == 'none')) sum(left$residuals^2) + sum(right$residuals^2) else NA
}
# scan all time steps and calculate sums of squares
# x : data vector
# scan.range: indices of a subset of time.vector to be scanned for the change point
# by default 3 : N-3
scan_2phase <- function(x, scan.range = NULL, required.direction = 'down') {
if (is.null(scan.range)) scan.range <- 3:(length(x) - 3)
ssq <- sapply(scan.range, ssq_left_right, x = x, required.direction = required.direction)
which.min(ssq) + scan.range[1] - 1
}
# Radial growth averaging
rga_cp <- function(x, N, lambda) {
DT <- data.table(x = x)
DT[, meanBef := frollmean(x, N, align = 'right')]
DT[, meanAft := (frollsum(x, N + 1, align = 'left') - x) / N]
DT[, ratio := meanAft / meanBef]
idx <- which.min(DT$ratio)
if (DT[idx, ratio] <= lambda) idx else NULL
}
# Likelihood ratio test
likelihood_cp <- function(x, alpha = 0.01, required.direction = 'down') {
fit <- cpt.meanvar(x, penalty = 'SIC', method = 'AMOC', class = FALSE, minseglen = 20)
if (fit[2] > 1 - alpha) { # significant
cp <- fit[1]
left <- mean(x[1:cp])
right <- mean(x[(cp + 1) : length(x)])
if (required.direction == 'down') {
if (left > right) fit[1] else NULL
} else {
if (left < right) fit[1] else NULL
}
} else NULL
}Illustration of two-phase regression (Figure 3)
First, apply two-phase regression to Chiang Mai rainfall data
## [1] 144Day chmCpt was detected to be the change point.
Now we plot the precipitation data, BI data, and their cumulatives. On the precipitation data we will plot the trends before and after day 144.
p1 <- ggplot(chm) +
geom_col(aes(date, pre), fill = 'steelblue') +
labs(x = NULL, y = 'Precipitation [mm]') +
scale_x_date(date_labels = '%d %b') +
scale_y_continuous(expand = c(0, 0)) +
theme(plot.tag.position = c(0.2, 0.97),
plot.margin = margin(r = 10))
p2 <- ggplot(chm) +
geom_line(aes(date, cumuP), linewidth = 0.6, color = 'steelblue') +
geom_smooth(aes(date, cumuP), data = chm[1:144],
method = 'lm', formula = 'y ~ x', se = FALSE,
color = 'darkorange', linetype = 2) +
geom_smooth(aes(date, cumuP), data = chm[144:250],
method = 'lm', formula = 'y ~ x', se = FALSE,
color = 'darkorange', linetype = 2) +
geom_point(aes(date, cumuP + 30), chm[144]) +
labs(x = NULL, y = 'Cumulative precipitation [mm]') +
scale_x_date(date_labels = '%d %b') +
scale_y_continuous(expand = c(0, 0)) +
theme(plot.tag.position = c(0.24, 0.97))
p3 <- ggplot(mck['MC24N']) +
geom_line(aes(year, lwbi), color = '#8c510a') +
labs(x = NULL, y = 'LWBI [-]') +
theme(plot.tag.position = c(0.2, 0.97),
plot.margin = margin(r = 10))
p4 <- ggplot(mck['MC24N'][year >= 1900]) +
geom_point(aes(year, clw), color = '#8c510a', size = 0.1) +
labs(x = NULL, y = 'Cumulative LWBI [-]') +
theme(plot.tag.position = c(0.24, 0.97))
p1 + p2 + p3 + p4 +
plot_annotation(tag_levels = 'A') 
Apply the 3 methods
# Remove the five subfossil cores that do not have a transition
noTransition <- c('MC20', 'MC20A', 'MC20B', 'MC28S', 'MC39B')
mckAll <- copy(mck)
mck <- mck[core %ni% noTransition]
# Likelihood
cptLrLW <- mck[, tail(.SD, 100)[likelihood_cp(lwbi), .(year.lr = year)], by = core]
cptLrEW <- mck[, tail(.SD, 100)[likelihood_cp(ewbi), .(year.lr = year)], by = core]
# Two phase - apply on cumulative
cpt2pLW <- mck[, tail(.SD, 100)[scan_2phase(clw), .(year.2p = year)], by = core]
cpt2pEW <- mck[, tail(.SD, 100)[scan_2phase(cew), .(year.2p = year)], by = core]
# Radial growth averaging
cptRgaLW <- mck[, tail(.SD, 100)[rga_cp(lwbi, 20, 0.95), .(year.ga = year)], by = core]
cptRgaEW <- mck[, tail(.SD, 100)[rga_cp(ewbi, 20, 0.95), .(year.ga = year)], by = core]
# Merge
cptLW <- cpt2pLW |>
merge(cptLrLW, by = 'core', all = TRUE) |>
merge(cptRgaLW, by = 'core', all = TRUE)
setkey(cptLW, core)
cptEW <- cpt2pEW |>
merge(cptLrEW, by = 'core', all = TRUE) |>
merge(cptRgaEW, by = 'core', all = TRUE)
setkey(cptEW, core)Change points for latewood
## Key: <core>
## core year.2p year.lr year.ga
## <char> <int> <int> <int>
## 1: MC22B 1963 1958 1957
## 2: MC22N 1958 1959 1961
## 3: MC23A 1938 NA 1940
## 4: MC23B 1912 NA NA
## 5: MC23W 1958 1944 1951
## 6: MC24 1946 1948 1948
## 7: MC24A 1974 1980 1981
## 8: MC24N 1945 1948 1948
## 9: MC24S 1968 NA 1979
## 10: MC25E 1949 1944 1944
## 11: MC29A 1948 1953 1953
## 12: MC29B 1968 1964 1963
## 13: MC30N 1982 1986 1986
## 14: MC31N 1976 1982 1983
## 15: MC32S 1953 1951 NA
## 16: MC34A 1967 1973 1974
## 17: MC35N 1952 1951 1951
## 18: MC35W 1955 1960 1960
## 19: MC36 1946 1945 1945
## 20: MC36N 1945 1949 1949
## 21: MC37S 1971 1975 1975
## 22: MC39A 1945 1946 1948
## 23: MC41S 1968 1968 1975
## 24: MC42A 1940 1933 1937
## 25: MC42B 1933 1930 1930
## 26: MC900 1956 1953 1953
## 27: MC901 1964 1964 1963
## 28: MC902 1941 1943 1945
## 29: MC903 1959 1949 1987
## 30: MC904 1961 1961 1959
## 31: MC908 1953 1952 1952
## 32: MC910 1957 1957 1957
## 33: MC912 1946 1947 1947
## 34: MC916 1957 1952 1953
## 35: MC917 1956 1958 1987
## 36: MC918 1958 1957 1987
## 37: MC919 1964 1952 1977
## 38: MC920 1954 1952 1952
## 39: MC922 1979 1979 1986
## 40: MC923 1968 1960 1987
## core year.2p year.lr year.gaChange points for earlywood
## Key: <core>
## core year.2p year.lr year.ga
## <char> <int> <int> <int>
## 1: MC22B 1965 1959 1957
## 2: MC22N 1962 1962 1962
## 3: MC23A 1973 1982 1980
## 4: MC23B 1976 1980 1980
## 5: MC23W 1952 1953 1957
## 6: MC24 1946 1943 1949
## 7: MC24A 1975 1982 1982
## 8: MC24N 1945 1949 1947
## 9: MC24S 1984 1980 1979
## 10: MC25E 1940 1944 1943
## 11: MC29A 1960 1953 NA
## 12: MC29B 1965 1966 1958
## 13: MC30N 1986 1987 1987
## 14: MC31N 1967 1972 1972
## 15: MC32S 1983 1981 1981
## 16: MC34A 1978 1974 1974
## 17: MC34W 1979 1977 1977
## 18: MC35N 1951 1951 1951
## 19: MC35W 1967 1961 1961
## 20: MC36 1950 1944 1944
## 21: MC36N 1942 1942 1943
## 22: MC37S 1977 1975 1975
## 23: MC39A 1946 1948 1949
## 24: MC41S 1971 1970 1974
## 25: MC42A 1942 1932 1934
## 26: MC42B 1925 1928 1928
## 27: MC900 1955 1954 1955
## 28: MC901 1961 1959 1959
## 29: MC902 1935 1940 1939
## 30: MC903 1955 1945 1979
## 31: MC904 1962 1962 1962
## 32: MC908 1955 1953 1953
## 33: MC910 1958 1959 1960
## 34: MC912 1946 1947 1947
## 35: MC916 1959 1954 1955
## 36: MC917 1962 1958 1957
## 37: MC918 1961 1958 1957
## 38: MC919 1966 1965 1966
## 39: MC920 1954 1952 1953
## 40: MC922 1967 1973 1973
## 41: MC923 1958 1958 1958
## core year.2p year.lr year.gaMerge results into a long data set.
cptLong <- rbind(
melt(cptLW[, .(type = 'LWBI', core, year.2p, year.lr, year.ga)],
id.vars = c('core', 'type'),
variable.name = 'method', value.name = 'year') |>
merge(mck[, .(core, year, bi = lwbi)], by = c('core', 'year')),
melt(cptEW[, .(type = 'EWBI', core, year.2p, year.lr, year.ga)],
id.vars = c('core', 'type'), variable.name = 'method', value.name = 'year') |>
merge(mck[, .(core, year, bi = ewbi)], by = c('core', 'year')))
# Assign factors to the `method` column for plotting
cptLong[, method := factor(substr(method, 6, 7), levels = c('ga', '2p', 'lr'))]
setkey(cptLong, core)Comparing the methods
MC904 to illustrate (Figure 4).
mc904 <- mck['MC904'][year >= 1900]
p1 <- ggplot(mc904) +
geom_line(aes(year, lwbi, color = 'LWBI'), linewidth = 0.5) +
geom_line(aes(year, ewbi, color = 'EWBI'), linewidth = 0.5) +
scale_color_manual(values = c('steelblue', illiniBlue),
guide = guide_legend(override.aes = list(linewidth = 0.8), order = 1)) +
ggnewscale::new_scale_color() +
geom_point(aes(year, bi, color = method, shape = method, size = method), cptLong['MC904']) +
scale_color_manual(
values = c('2p' = 'darkorange', lr = 'steelblue', ga = 'darkviolet'),
breaks = c('lr', 'ga', '2p'),
labels = c('2p' = 'Two-phase linear regression',
lr = 'Likelihood ratio',
ga = 'Radial growth averaging')) +
scale_shape_manual(
values = c('2p' = 16, lr = 17, ga = 18),
breaks = c('lr', 'ga', '2p'),
labels = c('2p' = 'Two-phase linear regression',
lr = 'Likelihood ratio', ga = 'Radial growth averaging'),
guide = guide_legend(override.aes = list(size = 2))) +
scale_size_manual(
values = c('2p' = 4, lr = 3.2, ga = 3),
breaks = c('lr', 'ga', '2p'),
labels = c('2p' = 'Two-phase linear regression',
lr = 'Likelihood ratio',
ga = 'Radial growth averaging')) +
geom_text(aes(year, bi, label = year, color = method), cptLong['MC904'],
hjust = 0, nudge_x = 1, show.legend = FALSE) +
scale_x_continuous(
guide = guide_prism_minor(),
breaks = seq(1900, 2000, 10),
labels = skip_label(2),
minor_breaks = seq(1900, 2010, 1),
expand = c(0, 0)) +
labs(x = NULL, y = 'BI [-]') +
theme(
legend.position = 'top',
legend.box = 'vertical',
legend.spacing.y = unit(2, 'mm'),
legend.box.margin = margin(),
legend.margin = margin(t = 0),
plot.tag.position = c(0.01, 0.95))
p2 <- image_ggplot(image_read('figures/MC904-illustration.jpg')) +
theme(plot.tag.position = c(0.01, 0.9))
p1 / p2 + plot_layout(heights = c(3, 1)) +
plot_annotation(tag_levels = 'A')
Plot for all cores (Figure 5)
# Read the visually determined data and correct the core IDs to match with the BI daat
cpVis <- fread('data/mck-manual-cp-year.csv')
cpVis[substr(core, 3, 3) == 'K',
core := paste0('MC', substr(core, 4, nchar(core)))]
cpVis[core == 'MC200', core := 'MC20A']
cpVis[core == 'MC201', core := 'MC20B']
cpVis[core == 'MC230', core := 'MC23A']
cpVis[core == 'MC231', core := 'MC23B']
cpVis[core == 'MC245', core := 'MC24A']
cpVis[core == 'MC290', core := 'MC29A']
cpVis[core == 'MC291', core := 'MC29B']
cpVis[core == 'MC345', core := 'MC34A']
cpVis[core == 'MC390', core := 'MC39A']
cpVis[core == 'MC392', core := 'MC39B']
setkey(cpVis, core)
cptLong <- merge(cptLong, cpVis, by = 'core')
cptLong[, diff := year - year.vi]
# Plot
ggplot(mck[year >= 1900]) +
geom_vline(aes(xintercept = year.vi), cpVis, linewidth = 0.6, color = 'gray80') +
geom_line(aes(year, ewbi, color = 'EWBI')) +
geom_line(aes(year, lwbi, color = 'LWBI')) +
scale_color_manual(
values = c('steelblue', illiniBlue),
guide = guide_legend(override.aes = list(linewidth = 0.8), order = 1)) +
ggnewscale::new_scale_color() +
geom_point(aes(year, bi, color = method, shape = method), cptLong) +
scale_color_manual(
values = c('2p' = 'darkorange', lr = 'steelblue', ga = 'darkviolet'),
breaks = c('lr', 'ga', '2p'),
labels = c('2p' = 'Two-phase linear regression',
lr = 'Likelihood ratio',
ga = 'Radial growth averaging')) +
scale_shape_manual(
values = c('2p' = 16, lr = 17, ga = 18),
breaks = c('lr', 'ga', '2p'),
labels = c('2p' = 'Two-phase linear regression',
lr = 'Likelihood ratio',
ga = 'Radial growth averaging'),
guide = guide_legend(override.aes = list(size = 2))) +
facet_wrap(~core, ncol = 7) +
scale_x_continuous(expand = c(0, 0), breaks = c(1925, 1950, 1975, 2000), labels = skip_label(2)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = NULL, y = 'BI [-]') +
theme(
axis.text.x = element_text(size = 9),
legend.box = 'vertical',
legend.spacing.y = unit(-1, 'mm'),
legend.position = 'top') +
panel_border('black', 0.2)
Accuracy plot (Figure 6)
cptSpread <- cptLong[, .(yearMin = min(year, na.rm = TRUE),
yearMax = max(year, na.rm = TRUE)),
by = .(core, type)]
cptAccuracy <- cptLong[, .(mu = round(mean(diff), 1),
m = round(median(diff), 1),
sigma = round(sd(diff), 1)), by = .(method, type)]
ggplot(cptLong) +
geom_vline(xintercept = 0, color = 'gray90', linewidth = 0.2) +
geom_histogram(
aes(x = diff, fill = method),
binwidth = 1) +
geom_hline(yintercept = 0, linewidth = 0.4) +
geom_text(
aes(x = -62, y = 14, label = paste0('m = ', m)),
cptAccuracy, hjust = 0, size = 3.5) +
geom_text(
aes(x = -62, y = 12.4, label = paste0('\u03bc = ', mu)),
cptAccuracy, hjust = 0, size = 3.5) +
geom_text(
aes(x = -62, y = 10.8, label = paste0('\u03c3 = ', sigma)),
cptAccuracy, hjust = 0, size = 3.5) +
facet_grid(
type ~ method,
labeller = labeller(.cols = c('2p' = 'Two-phase linear regression',
lr = 'Likelihood ratio',
ga = 'Radial growth averaging'))) +
scale_x_continuous(
expand = c(0, 0),
guide = guide_prism_minor(),
minor_breaks = seq(-70, 25, 5)) +
scale_y_continuous(expand = c(0, 0)) +
scale_fill_manual(
values = c('2p' = 'darkorange', lr = 'steelblue', ga = 'darkviolet'),
breaks = c('2p', 'ga', 'lr')) +
labs(x = 'Difference [years]', y = 'Count [-]') +
theme(
axis.text = element_text(size = 11),
axis.title = element_text(size = 11),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
strip.text = element_text(size = 11),
panel.grid.major.y = element_line('gray90', 0.2),
panel.spacing.y = unit(5, 'mm'),
legend.position = 'none') +
panel_border('gray30', 0.2)
Assign final change point based on the spread among 3 methods
goodCases <- cptSpread[yearMax - yearMin <= 5]
cptEW[core %in% goodCases[type == 'EWBI', core],
year.fi := median(c(year.2p, year.lr, year.ga)), by = core]
cptLW[core %in% goodCases[type == 'LWBI', core],
year.fi := median(c(year.2p, year.lr, year.ga)), by = core]
cptEW['MC903', year.fi := 1979]
cptLW['MC903', year.fi := 1979]
cptEW['MC917', year.fi := 1987]
cptLW['MC917', year.fi := 1987]
cptLW['MC918', year.fi := 1958]
cptLW['MC919', year.fi := 1966]
cptLW['MC922', year.fi := 1973]
cptEW['MC922', year.fi := 1973]
cptLW['MC923', year.fi := 1960]
cptEW['MCK22B', year.fi := 1958]
cptLW['MCK22B', year.fi := 1959]
cptEW['MCK230', year.fi := 1980]
cptLW['MCK230', year.fi := NA]
cptLW['MCK231', year.fi := NA]
cptLW['MCK23W', year.fi := 1953]
cptEW['MCK24', year.fi := 1949]
cptLW['MCK245', year.fi := 1981]
cptEW['MCK245', year.fi := 1982]
cptLW['MCK24S', year.fi := 1979]
cptEW['MCK290', year.fi := 1953]
cptEW['MCK291', year.fi := 1964]
cptLW['MCK31N', year.fi := 1976]
cptLW['MCK32S', year.fi := 1981]
cptLW['MCK345', year.fi := 1974]
cptEW['MCK35W', year.fi := 1961]
cptEW['MCK36', year.fi := 1944]
cptLW['MCK41S', year.fi := 1968]
cptLW['MCK42A', year.fi := 1933]
cptEW['MCK42A', year.fi := 1932]
cptFinal <- merge(cptEW[, .(core, year.ew = year.fi)],
cptLW[, .(core, year.lw = year.fi)],
by = 'core', all = TRUE) |>
merge(cpVis, by = 'core', all = TRUE)Building the final chronologies
Color adjustment
Using 20 years
mck[, lwbiAdj := {
s <- .BY[[1]]
cpt <- cptFinal[core == s, year.lw]
if (is.na(cpt)) {
lwbi
} else {
nAfter <- length(which(year >= cpt))
nBefore <- 20
tailDT <- tail(.SD, nBefore + nAfter)
mean1 <- tailDT[year <= cpt, mean(lwbi)]
sd1 <- tailDT[year <= cpt, sd(lwbi)]
mean2 <- tailDT[year > cpt, mean(lwbi)]
sd2 <- tailDT[year > cpt, sd(lwbi)]
c(.SD[year <= cpt, lwbi], .SD[year > cpt, (lwbi - mean2) / sd2 * sd1 + mean1])
}
}, by = core]
mck[, ewbiAdj := {
s <- .BY[[1]]
cpt <- cptFinal[core == s, year.ew]
if (is.na(cpt)) {
ewbi
} else {
nAfter <- length(which(year >= cpt))
nBefore <- 20
tailDT <- tail(.SD, nBefore + nAfter)
mean1 <- tailDT[year <= cpt, mean(ewbi)]
sd1 <- tailDT[year <= cpt, sd(ewbi)]
mean2 <- tailDT[year > cpt, mean(ewbi)]
sd2 <- tailDT[year > cpt, sd(ewbi)]
c(.SD[year <= cpt, ewbi], .SD[year > cpt, (ewbi - mean2) / sd2 * sd1 + mean1])
}
}, by = core]
ggplot(mck[year >= 1900]) +
geom_line(aes(year, lwbi, linetype = 'Unadjusted', linewidth = 'Unadjusted', color = 'LWBI')) +
geom_line(aes(year, lwbiAdj, linetype = 'Adjusted', linewidth = 'Adjusted', color = 'LWBI')) +
geom_line(aes(year, ewbi, linetype = 'Unadjusted', linewidth = 'Unadjusted', color = 'EWBI')) +
geom_line(aes(year, ewbiAdj, linetype = 'Adjusted', linewidth = 'Adjusted', color = 'EWBI')) +
facet_wrap(~core, ncol = 7) +
scale_x_continuous(expand = c(0, 0), breaks = c(1925, 1950, 1975, 2000), labels = skip_label(2)) +
scale_y_continuous(expand = c(0, 0)) +
scale_color_manual(values = c('#a6611a', '#018571')) +
scale_linewidth_manual(name = 'type', values = c(0.3, 0.2)) +
scale_linetype_manual(name = 'type', values = c(1, 2)) +
labs(x = NULL, y = 'BI [-]') +
theme(
axis.text = element_text(size = 9),
legend.position = 'top') +
panel_border('black', 0.2)
Using long-term mean and variance
This analysis is done following a reviewer’s request.
mck[, lwbiAdj2 := {
s <- .BY[[1]]
cpt <- cptFinal[core == s, year.lw]
if (is.na(cpt)) {
lwbi
} else {
nAfter <- length(which(year >= cpt))
mean1 <- .SD[year <= cpt, mean(lwbi)]
sd1 <- .SD[year <= cpt, sd(lwbi)]
mean2 <- .SD[year > cpt, mean(lwbi)]
sd2 <- .SD[year > cpt, sd(lwbi)]
c(.SD[year <= cpt, lwbi], .SD[year > cpt, (lwbi - mean2) / sd2 * sd1 + mean1])
}
}, by = core]
mck[, ewbiAdj2 := {
s <- .BY[[1]]
cpt <- cptFinal[core == s, year.ew]
if (is.na(cpt)) {
ewbi
} else {
nAfter <- length(which(year >= cpt))
mean1 <- .SD[year <= cpt, mean(ewbi)]
sd1 <- .SD[year <= cpt, sd(ewbi)]
mean2 <- .SD[year > cpt, mean(ewbi)]
sd2 <- .SD[year > cpt, sd(ewbi)]
c(.SD[year <= cpt, ewbi], .SD[year > cpt, (ewbi - mean2) / sd2 * sd1 + mean1])
}
}, by = core]
ggplot(mck[year >= 1900]) +
geom_line(aes(year, lwbiAdj, linetype = '20-year', linewidth = '20-year', color = 'LWBI')) +
geom_line(aes(year, lwbiAdj2, linetype = 'Long-term', linewidth = 'Long-term', color = 'LWBI')) +
geom_line(aes(year, ewbiAdj, linetype = '20-year', linewidth = '20-year', color = 'EWBI')) +
geom_line(aes(year, ewbiAdj2, linetype = 'Long-term', linewidth = 'Long-term', color = 'EWBI')) +
facet_wrap(~core, ncol = 7) +
scale_x_continuous(expand = c(0, 0), breaks = c(1925, 1950, 1975, 2000), labels = skip_label(2)) +
scale_y_continuous(expand = c(0, 0)) +
scale_color_manual(values = c('#a6611a', '#018571')) +
scale_linewidth_manual(name = 'type', values = c(0.3, 0.2)) +
scale_linetype_manual(name = 'type', values = c(1, 2)) +
labs(x = NULL, y = 'BI [-]') +
theme(
axis.text = element_text(size = 9),
legend.position = 'top') +
panel_border('black', 0.2)
Use core 42B as an example.
ggplot(mck[core == 'MC42B']) +
geom_line(aes(year, lwbiAdj, linetype = '20-year', linewidth = '20-year', color = 'LWBI')) +
geom_line(aes(year, lwbiAdj2, linetype = 'Long-term', linewidth = 'Long-term', color = 'LWBI')) +
geom_line(aes(year, ewbiAdj, linetype = '20-year', linewidth = '20-year', color = 'EWBI')) +
geom_line(aes(year, ewbiAdj2, linetype = 'Long-term', linewidth = 'Long-term', color = 'EWBI')) +
facet_wrap(~core, ncol = 7) +
scale_color_manual(values = c('#a6611a', '#018571')) +
scale_linewidth_manual(name = 'type', values = c(0.3, 0.2)) +
scale_linetype_manual(name = 'type', values = c(1, 2)) +
labs(x = NULL, y = 'BI [-]') +
theme(
axis.text = element_text(size = 9),
legend.position = 'top') +
panel_border('black', 0.2)
Detrending
Figure S1.
mckAdj <- rbind(mck[, .(core, year, lwbi = lwbiAdj, ewbi = ewbiAdj)],
mckAll[noTransition][, .(core, year, lwbi, ewbi)])
ewbiCrnAdj <- dt_to_rwl(mckAdj[, .(core, year, ewbi)], 'ewbi') |>
ssf(method = 'AgeDepSpline', nyrs = 50, verbose = FALSE) |>
crn_to_dt()
ewbiCrnRaw <- dt_to_rwl(mckAll[, .(core, year, ewbi)], 'ewbi') |>
ssf(method = 'AgeDepSpline', nyrs = 50, verbose = FALSE) |>
crn_to_dt()
lwbiCrnAdj <- dt_to_rwl(mckAdj[, .(core, year, lwbi)], 'lwbi') |>
ssf(method = 'AgeDepSpline', nyrs = 50, verbose = FALSE) |>
crn_to_dt()
lwbiCrnRaw <- dt_to_rwl(mckAll[, .(core, year, lwbi)], 'lwbi') |>
ssf(method = 'AgeDepSpline', nyrs = 50, verbose = FALSE) |>
crn_to_dt()
dbiCrn <- dt_to_rwl(mckAll[, .(core, year, dbi = lwbi - ewbi)], 'dbi') |>
ssf(method = 'AgeDepSpline', nyrs = 50, verbose = FALSE) |>
crn_to_dt()
crnRawMerged <- rbind(ewbiCrnRaw[samp.depth >= 5, .(year, sfc, variable = 'EWBI')],
lwbiCrnRaw[samp.depth >= 5, .(year, sfc, variable = 'LWBI')],
dbiCrn[samp.depth >= 5, .(year, sfc, variable = 'DBI')])
crnMerged <- rbind(ewbiCrnAdj[samp.depth >= 5, .(year, sfc, variable = 'EWBI')],
lwbiCrnAdj[samp.depth >= 5, .(year, sfc, variable = 'LWBI')],
dbiCrn[samp.depth >= 5, .(year, sfc, variable = 'DBI')])
crnMerged[, variable := factor(variable, levels = c('EWBI', 'LWBI', 'DBI'))]
crnRawMerged[, variable := factor(variable, levels = c('EWBI', 'LWBI', 'DBI'))]
ggplot() +
geom_line(aes(year, sfc, color = variable, linetype = 'Adjusted'), crnMerged[variable != 'DBI']) +
geom_line(aes(year, sfc, color = variable, linetype = 'Unadjusted'), crnRawMerged) +
facet_wrap(~variable, ncol = 1, scales = 'free_y', strip.position = 'right') +
scale_color_manual(values = c(EWBI = 'steelblue',
LWBI = illiniBlue,
DBI = 'darkorange')) +
labs(x = NULL, y = 'BI [-]') +
theme(strip.text.y.right = element_text(angle = 0)) +
theme(legend.position = 'none')
Correlation analysis
With TAS averaged
First, we prepare the climate data.
# Mickey Creek coordinates
mckCoo <- data.table(name = 'Mickey Creek',
latitude = -41.48,
longitude = 146.42)
# Read CRU data
cruTmp <- fread('data/cruTmp.csv')
cruxy <- cruTmp[, .SD[1], by = .(lon, lat)]
cruTmp[, ':='(year = year(time), month = month(time))]
# Calculate Jan (+1) temperature
cruTmpJ1 <- cruTmp[month == 1, .(tmp = mean(tmp)), by = .(lon, lat, year)]
# Schulman-shift the year back to match the tree ring data
cruTmpJ1[, year := year - 1]
# TAS average
cruTmpJ1Tas <- cruTmpJ1[lon %between% c(144, 150) & lat %between% c(-44, -40.5), .(tmp = mean(tmp)), by = year]The following grid cells are used to calculate the TAS average.
ggplot() +
geom_point(aes(longitude, latitude), mckCoo, color = 'violet', shape = 'triangle', size = 2.5) +
geom_tile(aes(lon, lat), cruxy, fill = NA, color = 'gray', width = 0.5, height = 0.5, linewidth = 0.1) +
geom_polygon(aes(long, lat, group = group), map_data('world', 'Australia'), color = 'steelblue') +
scale_x_continuous(labels = degree_lon, name = NULL) +
scale_y_continuous(labels = degree_lat, name = NULL) +
coord_quickmap(xlim = c(144, 150), ylim = c(-44, -40.5)) +
panel_border('black', 0.2)
Now we visualize the chronologies and TAS-averaged Jan (+1) temperatures since 1901. Both temperatures and BI data are converted to z-scores to have a common scale for comparison. Because of the negative correlations, the BI data are inverted.
tasMerged <- merge(cruTmpJ1Tas, crnMerged, by = 'year')
ggplot(tasMerged) +
geom_line(aes(year, standardize(tmp)), linetype = 2) +
geom_line(aes(year, -standardize(sfc), color = variable)) +
facet_wrap(~variable, ncol = 1, strip.position = 'right') +
scale_color_manual(values = c(EWBI = 'steelblue',
LWBI = illiniBlue,
DBI = 'darkorange')) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = NULL, y = 'z-score') +
theme(strip.text.y.right = element_text(angle = 0),
legend.position = 'none')
Correlations without first differencing:
## variable V1
## <fctr> <num>
## 1: EWBI -0.5691474
## 2: LWBI -0.4228773
## 3: DBI -0.1751944Correlation with first differencing:
## variable V1
## <fctr> <num>
## 1: EWBI -0.59763960
## 2: LWBI -0.54459215
## 3: DBI -0.08063125Gridded correlations
With first-differencing
rhoDT <- merge(cruTmpJ1, crnMerged[, .(year, variable, sfc)], by = 'year', allow.cartesian = TRUE
)[, {
ct <- cor.test(diff(tmp), diff(sfc), 'two.sided', 'pearson')
list(rho = ct$estimate, p.value = ct$p.value)
}, by = .(variable, lon, lat)]
rhoDT[, variable := factor(variable, c('EWBI', 'LWBI', 'DBI'))]
crnMergedWide <- dcast(crnMerged, year ~ variable, value.var = 'sfc')
biCorr <- cor(crnMergedWide[, -'year']) |>
round(2) |>
as.data.table(keep.rownames = 'var1') |>
melt(id.vars = 'var1', variable.name = 'var2')
biCorr[, var1 := factor(var1, levels = levels(var2), ordered = TRUE)]
biCorr[, var2 := factor(var2, levels = levels(var2), ordered = TRUE)]
p1 <- ggplot(rhoDT[p.value < 0.05]) +
geom_tile(aes(lon, lat, fill = rho)) +
geom_polygon(aes(long, lat, group = group),
map_data('world', 'Australia'),
fill = NA, color = 'gray30') +
geom_point(aes(longitude, latitude), mckCoo, color = 'violet', shape = 'triangle', size = 2.5) +
facet_wrap(~variable, drop = FALSE, ncol = 2) +
scale_fill_distiller(name = 'Correlation',
palette = 'RdBu',
breaks = seq(-0.6, 0.6, 0.2) |> round(1),
limits = abs_range(rhoDT$rho)) +
scale_x_continuous(expand = c(0, 0),
labels = degree_lon, breaks = seq(120, 150, 10)) +
scale_y_continuous(expand = c(0, 0),
labels = degree_lat) +
labs(x = NULL, y = NULL) +
coord_quickmap(xlim = c(110, 157), ylim = c(-47, -8)) +
panel_border('black', 0.2) +
tag_facets(tag_levels = 'A', tag_suffix = '') +
theme(
strip.text = element_text(size = 12),
legend.title = element_text(),
legend.key.width = unit(5, 'mm'),
legend.key.height = unit(2.25, 'cm'))
p2 <- ggplot(biCorr[var1 < var2]) +
geom_text(aes(var1, var2, label = value)) +
scale_y_discrete(limits = c('DBI', 'LWBI')) +
geom_tile(aes(var1, var2), color = 'gray10', fill = NA) +
labs(x = NULL, y = NULL, tag = 'D',
subtitle = 'Cross-correlation') +
theme(axis.line = element_blank(),
axis.ticks = element_blank(),
plot.tag.position = c(-0.16, 1.23),
plot.tag = element_text(size = 11))
layout <- c(
area(1, 1, 16, 16),
area(12, 11, 15, 14)
)
p1 + p2 + plot_layout(design = layout)
Without first differencing
rhoDT2<- merge(cruTmpJ1, crnMerged[, .(year, variable, sfc)], by = 'year', allow.cartesian = TRUE
)[, {
ct <- cor.test(tmp, sfc, 'two.sided', 'pearson')
list(rho = ct$estimate, p.value = ct$p.value)
}, by = .(variable, lon, lat)]
rhoDT2[, variable := factor(variable, c('EWBI', 'LWBI', 'DBI'))]
rhoDT[, type := 'First-differenced']
rhoDT2[, type := 'No differencing']
rhoDTmerged <- rbind(rhoDT, rhoDT2)
ggplot(rhoDTmerged[p.value < 0.05]) +
geom_tile(aes(lon, lat, fill = rho)) +
geom_polygon(aes(long, lat, group = group),
map_data('world', 'Australia'),
fill = NA, color = 'gray30') +
geom_point(aes(longitude, latitude), mckCoo, color = 'violet', shape = 'triangle', size = 2.5) +
facet_grid(variable ~ type, drop = FALSE) +
scale_fill_distiller(name = 'Correlation',
palette = 'RdBu',
breaks = seq(-0.6, 0.6, 0.2) |> round(1),
limits = abs_range(rhoDTmerged$rho)) +
scale_x_continuous(expand = c(0, 0),
labels = degree_lon, breaks = seq(120, 150, 10)) +
scale_y_continuous(expand = c(0, 0),
labels = degree_lat) +
labs(x = NULL, y = NULL) +
coord_quickmap(xlim = c(110, 157), ylim = c(-47, -8)) +
panel_border('black', 0.2) +
tag_facets(tag_levels = 'A', tag_suffix = '') +
theme(
strip.text = element_text(size = 12),
legend.title = element_text(),
legend.key.width = unit(5, 'mm'),
legend.key.height = unit(2.25, 'cm'))
With ADS on DBI at different window lengths
This is in response to Reviewer 2.
dbiCrnRawADS <- lapplyrbind(c(10, 25, 50, 75),
\(k) {
out <- dt_to_rwl(mckAll[, .(core, year, dbi = lwbi - ewbi)], 'dbi') |>
detrend(method = 'AgeDepSpline', nyrs = k, verbose = FALSE) |>
chron() |>
crn_to_dt()
out[, nyrs := paste(k, 'years')]
})
rhoDT3<- merge(cruTmpJ1, dbiCrnRawADS[, .(nyrs, year, std)], by = 'year', allow.cartesian = TRUE
)[, {
ct <- cor.test(diff(tmp), diff(std), 'two.sided', 'pearson')
list(rho = ct$estimate, p.value = ct$p.value)
}, by = .(nyrs, lon, lat)]
ggplot(rhoDT3[p.value < 0.05]) +
geom_tile(aes(lon, lat, fill = rho)) +
geom_polygon(aes(long, lat, group = group),
map_data('world', 'Australia'),
fill = NA, color = 'gray30') +
geom_point(aes(longitude, latitude), mckCoo, color = 'violet', shape = 'triangle', size = 2.5) +
facet_wrap(~nyrs) +
scale_fill_distiller(name = 'Correlation',
palette = 'RdBu',
breaks = seq(-0.6, 0.6, 0.2) |> round(1),
limits = abs_range(rhoDT$rho)) +
scale_x_continuous(expand = c(0, 0),
labels = degree_lon, breaks = seq(120, 150, 10)) +
scale_y_continuous(expand = c(0, 0),
labels = degree_lat) +
labs(x = NULL, y = NULL) +
coord_quickmap(xlim = c(110, 157), ylim = c(-47, -8)) +
panel_border('black', 0.2) +
tag_facets(tag_levels = 'A', tag_suffix = '') +
theme(
strip.text = element_text(size = 12),
legend.title = element_text(),
legend.key.width = unit(5, 'mm'),
legend.key.height = unit(2.25, 'cm'))
dbiCrnRawSSF <- lapplyrbind(c(10, 25, 50, 75),
\(k) {
out <- dt_to_rwl(mckAll[, .(core, year, dbi = lwbi - ewbi)], 'dbi') |>
ssf(method = 'AgeDepSpline', nyrs = k, verbose = FALSE) |>
crn_to_dt()
out[, nyrs := paste(k, 'years')]
})
rhoDT4<- merge(cruTmpJ1, dbiCrnRawSSF[, .(nyrs, year, sfc)], by = 'year', allow.cartesian = TRUE
)[, {
ct <- cor.test(diff(tmp), diff(sfc), 'two.sided', 'pearson')
list(rho = ct$estimate, p.value = ct$p.value)
}, by = .(nyrs, lon, lat)]
rhoDT4[, nyrs := factor(nyrs)]
ggplot(rhoDT4[p.value < 0.05]) +
geom_tile(aes(lon, lat, fill = rho)) +
geom_polygon(aes(long, lat, group = group),
map_data('world', 'Australia'),
fill = NA, color = 'gray30') +
geom_point(aes(longitude, latitude), mckCoo, color = 'violet', shape = 'triangle', size = 2.5) +
facet_wrap(~nyrs, drop = FALSE) +
scale_fill_distiller(name = 'Correlation',
palette = 'RdBu',
breaks = seq(-0.6, 0.6, 0.2) |> round(1),
limits = abs_range(rhoDT$rho)) +
scale_x_continuous(expand = c(0, 0),
labels = degree_lon, breaks = seq(120, 150, 10)) +
scale_y_continuous(expand = c(0, 0),
labels = degree_lat) +
labs(x = NULL, y = NULL) +
coord_quickmap(xlim = c(110, 157), ylim = c(-47, -8)) +
panel_border('black', 0.2) +
tag_facets(tag_levels = 'A', tag_suffix = '') +
theme(
strip.text = element_text(size = 12),
legend.title = element_text(),
legend.key.width = unit(5, 'mm'),
legend.key.height = unit(2.25, 'cm'))
Graphical abstract
Core MC904 but with adjustments too.
p1 <- ggplot(mck['MC904'][year >= 1900]) +
geom_line(aes(year, lwbi, color = 'LWBI', linetype = 'Unadjusted'), linewidth = 0.5) +
geom_line(aes(year, lwbiAdj, color = 'LWBI', linetype = 'Adjusted'), linewidth = 0.5) +
geom_line(aes(year, ewbi, color = 'EWBI', linetype = 'Unadjusted'), linewidth = 0.5) +
geom_line(aes(year, ewbiAdj, color = 'EWBI', linetype = 'Adjusted'), linewidth = 0.5) +
scale_color_manual(values = c('steelblue', illiniBlue),
guide = guide_legend(override.aes = list(linewidth = 0.8), order = 1)) +
scale_linetype_manual(values = c(Adjusted = 1, Unadjusted = 2),
guide = guide_legend(order = 2)) +
ggnewscale::new_scale_color() +
geom_point(aes(year, bi, color = method, shape = method, size = method), cptLong['MC904']) +
scale_color_manual(
values = c('2p' = 'darkorange', lr = 'steelblue', ga = 'darkviolet'),
breaks = c('lr', 'ga', '2p'),
labels = c('2p' = 'Two-phase linear regression',
lr = 'Likelihood ratio',
ga = 'Radial growth averaging')) +
scale_shape_manual(
values = c('2p' = 16, lr = 17, ga = 18),
breaks = c('lr', 'ga', '2p'),
labels = c('2p' = 'Two-phase linear regression',
lr = 'Likelihood ratio',
ga = 'Radial growth averaging'),
guide = guide_legend(override.aes = list(size = 2))) +
scale_size_manual(
values = c('2p' = 4, lr = 3.2, ga = 3),
breaks = c('lr', 'ga', '2p'),
labels = c('2p' = 'Two-phase linear regression',
lr = 'Likelihood ratio',
ga = 'Radial growth averaging')) +
geom_text(aes(year, bi, label = year, color = method), cptLong['MC904'],
hjust = 0, nudge_x = 1, show.legend = FALSE) +
scale_x_continuous(
guide = guide_prism_minor(),
breaks = seq(1900, 2000, 10),
labels = skip_label(2),
minor_breaks = seq(1900, 2010, 1),
expand = c(0, 0)) +
labs(x = NULL, y = 'BI [-]') +
theme(
legend.position = 'top',
legend.box = 'vertical',
legend.key.width = unit(10, 'mm'),
legend.box.margin = margin(),
legend.margin = margin(t = 0),
plot.tag.position = c(0.01, 0.95))
p2 <- image_ggplot(image_read('figures/MC904-illustration.jpg')) +
theme(plot.tag.position = c(0.01, 0.9))
p1 / p2 + plot_layout(heights = c(3, 1)) +
plot_annotation(tag_levels = 'A')