paper-code.Rmd

---
title: "Code for *Tree-Ring Evidence of Increasing Drought Risks amidst Projected Flood Intensification in the Kabul River Basin (Afghanistan and Pakistan)* by Khan et al. (2022)"
author: "Hung Nguyen"
date: "2022-05-14"
output: 
  rmdformats::readthedown: default 
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Introduction and Preparations

This document details the process of producing the results presented in **Tree-Ring Evidence of Increasing Drought Risks amidst Projected Flood Intensification in the Kabul River Basin (Afghanistan and Pakistan)** by Khan et al. (2022).

To reproduce the results, please do the following:

-   This code requires R 4.1.0 and above.
-   Download the code repository from the GitHub repo and extract the downloaded `.zip` file to your working folder.
-   Open `chitral-precip.Rproj` in RStudio (**It's important to open this first so that the file path is loaded properly**).
-   Install and load the following packages if you don't already have them. For package `ldsr`, please use the development version which can be installed from GitHub with

```{r, eval=FALSE}
remotes::install_github('ntthung/ldsr')
```

```{r, message=FALSE, warning=FALSE}
library(dplR)         # Tree ring data processing
library(ldsr)         # Tree ring data processing
library(data.table)   # Data handling 
library(missMDA)      # Imputation
library(qmap)         # Bias correction
library(quantreg)     # Quantile regression
library(ggplot2)      # Plotting
library(cowplot)      # Plotting
library(patchwork)    # Plotting
library(ggprism)      # Plotting
library(ggnewscale)   # Plotting
library(RColorBrewer) # Color palette
theme_set(theme_prism(base_size = 10, base_fontface = 'plain', base_line_size = 0.2))
```

-   Open `paper-code.Rmd`, which is the source code for this document.
-   Follow the written details below and run the code chunks one by one.

This R Markdown is set to render both HTML and PDF outputs. To do so, please run

```{r, eval=FALSE}
rmarkdown::render(
  'paper-code.Rmd', 
  output_format = 'all', 
  output_options = list(hightlight = 'tango', 
                        latex_engine = 'xelatex'))
```

For quick access to the final results please see the `.csv` file in `results/`. 

The code utilities to support the main code are stored in the folder `R/`. We need to load them first before running the main code.

```{r}
source('R/init.R')
source('R/correlation_functions.R')
source('R/drought_analysis_functions.R')
```

# Data

## Climate data

### Precipitation and Temperature

```{r climate data, warning=FALSE, message=FALSE}
# Monthly temperature
TmWide <- fread('data/chitral-monthly-T.csv')
Tm <- melt(TmWide, 
           id.var = 'year', 
           variable.name = 'month2', 
           variable.factor = TRUE, 
           value.name = 'Tm'
         )[, month := as.integer(month2)
         ][order(year, month)]
# Monthly precipitation
PmWide <- fread('data/chitral-monthly-P.csv')
Pm <- melt(PmWide, 
           id.var = 'year', 
           variable.name = 'month2', 
           variable.factor = TRUE, 
           value.name = 'Pm'
         )[, month := as.integer(month2)
         ][order(year, month)]
```

**Figure 1c**

```{r, fig.width=7, fig.height=3}
p1 <- ggplot(Pm) +
  geom_boxplot(aes(month2, Pm), fill = 'steelblue', alpha = 0.5) +
  labs(x = NULL, y = 'Precipitation [mm]', tag = 'c)') +
  scale_x_discrete(labels = \(x) substr(x, 1, 1)) +
  scale_y_continuous(guide = guide_prism_minor()) +
  panel_border('black', 0.2) 

p2 <- ggplot(Tm) +
  geom_boxplot(aes(month2, Tm), fill = 'darkorange', alpha = 0.5) +
  labs(x = NULL, y = 'Temperature [\u00b0C]') +
  scale_x_discrete(labels = \(x) substr(x, 1, 1)) +
  scale_y_continuous(guide = guide_prism_minor()) +
  panel_border('black', 0.2)
p1 + p2
```

### Wind vectors

Note: reanalysis and GPCP data are huge so they won't be included in this repository. You can download them from NOAA in order to run the code below, which produces Figure 1d.

```{r, eval=FALSE}
# Set map limit
xrange <- c(0, 80)
yrange <- c(0, 50)

world <- map('world', fill = TRUE, plot = FALSE) |> 
  st_as_sf(coords = c('x', 'y')) |> 
  st_buffer(0.01) |>
  st_union()
  
pak <- map('world', 'Pakistan', fill = TRUE, plot = FALSE) |> 
  st_as_sf(coords = c('xy', 'y')) 

## NCEP wind data
# Wind vector: read u-wind and v-wind, then merge
# Change the file path to your own
ncep <- merge(
  # u-wind
  ReadNetCDF('data/uwnd.mon.mean.nc', 'uwnd', 
             subset = list(level = 850, 
                           lon = xrange, lat = yrange,
                           time = c('1979-01-01', '2021-12-01'))),
  ReadNetCDF('data/vwnd.mon.mean.nc', 'vwnd', 
             subset = list(level = 850, 
                           lon = xrange, lat = yrange,
                           time = c('1979-01-01', '2021-12-01'))),
  by = c('time', 'level', 'lat', 'lon'))
# Create seasonal mean. 
ncepSeas <- ncep[month(time) %in% c(2:4, 6:8), 
                 .(uwnd = mean(uwnd), vwnd = mean(vwnd)), 
                 by = .(lon, lat, seas = fifelse(month(time) %in% 2:4, 'Feb-Apr', 'Jun-Aug'))]

## GPCP precipitation data
pre <- ReadNetCDF('data/precip.mon.mean.nc', 'precip',
                  subset = list(lon = xrange, lat = yrange,
                           time = c('1979-01-01', '2015-12-01')))
# Create seasonal mean. 
preSeas <- pre[month(time) %in% c(2:4, 6:8), 
               .(precip = mean(precip)), 
               by = .(lon, lat, seas = fifelse(month(time) %in% 2:4, 'Feb-Apr', 'Jun-Aug'))]
```

```{r wind plot, fig.width=7, fig.height=3.5, eval=FALSE}
maxSpeed <- ncepSeas[, sqrt(uwnd^2 + vwnd^2) |> max() |> ceiling()]
maxRain <- preSeas[, ceiling(max(precip))]
p3 <- ggplot() +
  # Precip
  geom_contour_fill(
    aes(lon, lat, z = precip, fill = stat(level_d)),
    data = preSeas,
    breaks = seq(0, maxRain, 2)) +
  # Land mass
  geom_sf(data = world, size = 0.2, fill = NA) +
  geom_sf(data = pak, size = 0.4, fill = NA, color = 'firebrick') +
  # Arrow vector
  geom_arrow(aes(lon, lat, dx = uwnd, dy = vwnd), 
             data = ncepSeas,
             # min.mag = 1,
             # skip = 1,
             color = blues9[7],
             size = 0.2,
             arrow.length = 0.3) +
  # Site locations
  geom_point(aes(Longitude, Latitude), crnMeta, color = 'darkviolet') +
  # Facet by seasons
  facet_wrap(vars(seas)) +
  # Set scales
  as.discretised_scale(scale_fill_distiller)(
    name = 'Precipitation rate [mm/day]',
    palette = 'Greens', direction = 1, breaks = seq(0, maxRain, 2),
    guide = guide_colorsteps(
      # override.aes = list(alpha = 0.5),
      title.position = 'top', 
      show.limits = TRUE)) + 
  scale_x_continuous(expand = c(0, 0), labels = paste_long) +
  scale_y_continuous(expand = c(0, 0), labels = paste_lat) +
  scale_mag(
    name = 'Wind speed [m/s]', 
    guide = guide_vector(title.position = 'top'),
    max = 20) +
  # guides(alpha = guide_none()) +
  # Axis labels
  labs(x = NULL, y = NULL, tag = 'd)') +
  # Map projection
  coord_sf(xlim = xrange + c(1, -1.25), ylim = yrange + c(1.25, 0)) +
  panel_border('black', 0.2) +
  theme(
    strip.text = element_text(face = 'bold'),
    legend.position = 'bottom',
    legend.key.width = unit(1.25, 'cm'),
    legend.title = element_text()) 
p3
# ggsave('precip-wind-pakistan.png', width = 7, height = 3.5)
```

## Tree ring data

```{r}
# Read ARSTAN outputs
crnRaw <- lapplyrbind(
  list.files('data/', '.tabs', full.names = TRUE),
  function(fn) {
    dt <- fread(fn)
    dt[, site := substr(fn, 6, nchar(fn) - 5)]
  })
setkey(crnRaw, site)
sssOut <- fread('data/sss.csv', key = 'site')
crn <- merge(crnRaw, sssOut, by = c('site', 'year'))

firstYear <- crn[sss > 0.6][, .N, by = year][N >= 4, first(year)]
lastYear <- 2018
# Use the residual chronology
crnWide <- crn[year %in% firstYear:lastYear, dcast(.SD, year ~ site, value.var = c('res'))]
```

**Figure S2 - Chronology and SSS plot**

```{r, fig.width=7, fig.height=6}
sf <- 2.2
ggplot(sssOut) +
  geom_hline(aes(yintercept = 1, linetype = 'Mean', color = 'Mean')) +
  # geom_hline(aes(yintercept = 0.85 * sf, 
                 # color = 'SSS = 0.85', linetype = 'SSS = 0.85')) +
  geom_hline(aes(yintercept = 0.6 * sf,
                 color = 'SSS = 0.6', linetype = 'SSS = 0.6')) +
  geom_line(aes(year, sss * sf, 
                col = 'Subsample signal strength', 
                linetype = 'Subsample signal strength'), size = 0.5) +
  geom_line(aes(year, res, 
                col = 'Ring width index', 
                linetype = 'Ring width index'), crn) +
  facet_wrap(~site, ncol = 2) +
  scale_x_continuous(
    name = NULL,
    breaks = seq(1300, 2000, 100),
    minor_breaks = seq(1300, 2000, 50),
    guide = guide_prism_minor(),
    expand = c(0, 0)) +
  scale_y_continuous(

    sec.axis = sec_axis(~ . / sf, name = 'Subsample signal strength',
                        breaks = seq(0, 1, 0.25)),
    expand = c(0, 0)) +
  scale_color_manual(
    name = NULL,
    breaks = c('Ring width index', 'Mean',
               'Subsample signal strength', 'SSS = 0.6'),
    values = c('steelblue', 'black', 'darkorange', 'black')) +
  scale_linetype_manual(
    name = NULL,
    breaks = c('Ring width index', 'Mean',
               'Subsample signal strength', 'SSS = 0.6'),
    values = c(1, 1, 1, 2)) +
  guides(color = guide_legend(override.aes = list(size = 0.50))) +
  theme(
    strip.background = element_rect('gray95', NA),
    legend.position = c(0.75, 0.1),
    legend.key.width = unit(1, 'cm'),
    panel.border = element_rect(NA, 'black', 0.2),
    axis.ticks.y.right = element_line(color = 'darkorange'),
    axis.text.y.right = element_text(color = 'darkorange', size = 9),
    axis.title.y.right = element_text(color = 'darkorange', size = 9),
    axis.line.y.right = element_line(color = 'darkorange'),
    axis.ticks.y.left = element_line(color = 'steelblue'),
    axis.text.y.left = element_text(color = 'steelblue', size = 9),
    axis.title.y.left = element_text(color = 'steelblue', size = 9),
    axis.line.y.left = element_line(color = 'steelblue'))
```

## Climate-growth relationship

First we infill the tree ring data.

```{r infilling}
X <- crnWide[, -'year'] |> 
  as.matrix() |> 
  imputePCA(ncp = 6) |> 
  {\(x) x$completeObs}() 
PC <- prcomp(X, scale. = TRUE)$x
Xraw <- as.matrix(crnWide[, -'year'])
impModel <- imputePCA(Xraw, ncp = 6)
Xfilled <- impModel$completeObs
Xfitted <- impModel$fittedX
colnames(Xfitted) <- colnames(Xfilled)

crnFitted <- as.data.table(as.data.frame(Xfitted))
crnFitted[, year := 1517:2018]
crnFittedLong <- melt(
  crnFitted, id.vars = 'year', variable.name = 'site', value.name = 'rwi')
crnMerge <- rbindlist(
  list(Infilled = crnFittedLong, Original = crn[, .(site, year, rwi = res)]),
  use.names = TRUE,
  idcol = 'type')
```

**Figure S5**

```{r, fig.width=7, fig.height=7}
ggplot(crnMerge[year >= 1517]) +
  geom_line(aes(year, rwi, color = type)) +
  facet_wrap(vars(site), ncol = 2, scales = 'free_y') +
  scale_color_manual(values = c('darkorange', 'steelblue')) +
  scale_x_continuous(breaks = seq(1500, 2000, 100)) +
  labs(x = NULL, y = 'RWI [-]') +
  panel_border('black', 0.2) +
  theme(legend.position = c(0.6, 0.1))
```

Calculate bootstrapped correlations.

```{r}
treeYears <- crnWide$year
# Current year
# Merge TR in 1965-2018 with precipitation in 1965-2018
instIndc <- which(treeYears %in% 1965:2018) 
XYc <- cbind(Xfilled[instIndc, ], as.matrix(PmWide[, -'year'])) 

# Tree rings and previous year precipitation
# Merge TR in 1966-2018 with precipitation in 1965-2017
instIndp <- which((treeYears - 1) %in% 1965:2017) # previous year streamflow
XYp <- cbind(Xfilled[instIndp, ], as.matrix(PmWide[-.N, -'year']))

# Tree rings and next year precipitation
# Merge TR in 1963-2016 with precipitation in 1965-2018
# instIndn <- which((treeYears + 1) %in% 1965:2018) # next year streamflow
# XYn <- cbind(Xfilled[instIndn, ], as.matrix(PmWide[, -'year']))

set.seed(2022)
corDTc <- cor_boot(XYc, 1:7, 8:19, groupNames = c('site', 'month')) # current year
corDTp <- cor_boot(XYp, 1:7, 8:19, groupNames = c('site', 'month')) # current year
# corDTn <- cor_boot(XYn, 1:7, 8:19, groupNames = c('site', 'month')) # current year

corDTc[, month := paste0(month, 'c')]
corDTp[, month := paste0(month, 'p')]
# corDTn[, month := paste0(month, 'n')]

corDTPm <- rbind(corDTp, corDTc)
corDTPm[, month := factor(month, c(paste0(month.abb, 'p'), paste0(month.abb, 'c')))]
```

**Figure 2**

```{r, fig.width=7, fig.height=6}
ggplot(corDTPm) +
  geom_linerange(
    aes(x = month, ymin = low, ymax = high, alpha = signif), color = 'steelblue') +
  geom_point(
    aes(x = month, y = rho0, alpha = signif), color = 'steelblue') +
  geom_hline(aes(yintercept = 0)) +
  geom_vline(xintercept = 12.5) +
  facet_wrap(vars(site), ncol = 2) +
  scale_x_discrete(name = NULL, labels = monthLab) +
  scale_y_continuous(name = 'Correlation [-]', 
                     breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6)) +
  scale_alpha_manual(
    values = c(0.25, 1), labels = c('Not significant', 'Significant')) +
  theme(
    legend.position = c(0.75, 0.1),
    axis.text.x = ggtext::element_markdown()) +
  panel_border('black', 0.2)
```

# September to August precipitation reconstruction

## Principal compnent analysis

```{r}
PCfit <- prcomp(X, scale. = TRUE)
PC <- PCfit$x
pcDT <- as.data.table(as.data.frame(summary(PCfit)$importance),
                      keep.rownames = 'var') |> 
  melt(id.vars = 'var', variable.name = 'PC')

# Extract PC1 and PC2 for plotting
pc12plot <- as.data.table(PC[, 1:2, drop = FALSE])[, year := crnWide$year]
```

**Figure S6**

```{r, fig.width=7, fig.height=6}
p1 <- ggplot(pcDT[var == 'Standard deviation']) +
  geom_col(aes(PC, value), color = 'black', fill = 'gray95') +
  geom_hline(yintercept = 1, size = 0.8) +
  scale_y_continuous(expand = expansion(add = c(0, 0.01)),
                     limits = c(0, 2)) +
  labs(x = NULL, y = 'Eigenvalue [-]', tag = 'a)') +
  theme(
    panel.grid.major.y = element_line('gray'),
    panel.ontop = TRUE,
    axis.ticks = element_blank(),
    axis.line = element_blank())
p2 <- ggplot() +
  geom_col(
    aes(PC, value * 100), 
    pcDT[var == 'Proportion of Variance'],
    fill = 'steelblue', color = 'black', alpha = 0.8) +
  geom_line(
    aes(PC, value * 100, group = 1),
    pcDT[var == 'Cumulative Proportion'],
    color = 'firebrick', size = 0.8) +
  geom_point(
    aes(PC, value * 100),
    pcDT[var == 'Cumulative Proportion'],
    color = 'firebrick', size = 1.6) +
  labs(x = NULL, y = 'Proportion of variance [%]', tag = 'b)') +
  scale_y_continuous(expand = expansion(add = c(0, 1))) +
  theme(
    panel.grid.major.y = element_line('gray'),
    panel.ontop = TRUE,
    axis.ticks = element_blank(),
    axis.line = element_blank())

p3 <- ggplot(pc12plot) +
  geom_line(aes(year, PC1), color = blues9[9]) +
  scale_x_continuous(
    expand = c(0, 0),
    minor_breaks = seq(1500, 2025, 25),
    limits = c(1500, 2025),
    guide = guide_prism_minor()) +
  labs(x = NULL, tag = 'c)') 

p4 <- ggplot(pc12plot) +
  geom_line(aes(year, PC2), color = blues9[9]) +
  scale_x_continuous(
    expand = c(0, 0),
    minor_breaks = seq(1500, 2025, 25),
    limits = c(1500, 2025),
    guide = guide_prism_minor()) +
  labs(x = NULL, tag = 'd)') 


# layout <- '
# AB
# CC
# DD
# '

layout <- c(
  patchwork::area(1, 1, 3, 1),
  patchwork::area(1, 2, 3, 2),
  patchwork::area(4, 1, 5, 2),
  patchwork:: area(6, 1, 7, 2)
)
p1 + p2 + p3 + p4 +
  plot_layout(design = layout)
```

## Reconstruction

```{r}
Pm[, year2 := fifelse(month %in% 9:12, year + 1, year)]
PSepAug <- Pm[, .(Qa = sum(Pm)), by = .(year = year2)][year %in% 1966:lastYear]
instInd <- which(firstYear:lastYear %in% PSepAug$year)

# Stepwise linear regression
DT <- cbind(PC[450:502, 1:3], Qa = PSepAug$Qa) |> as.data.frame()
step(lm(Qa ~ ., DT), direction = 'backward')
```

**Figure S1 - Autocorrelation**

```{r}
acf(PSepAug$Qa, main = NA, xlab = 'Lag [years]')
```

**Build final model**

```{r}
sPC <- as.data.table(PC[, 1:2])
lmFit <- PCR_reconstruction(PSepAug, sPC, firstYear, transform = 'none')
set.seed(42)
cvFolds <- make_Z(PSepAug$Qa, frac = 0.25, contiguous = TRUE)
lmCV <- cvPCR(PSepAug, sPC, firstYear, transform = 'none',
              Z = cvFolds, use.robust.mean = FALSE)
```

**Performance scores**

```{r}
round(lmCV$metrics[, 1:3], 2)
```
**Figure S7 - score distribution**

```{r, fig.width=7, fig.height=5}
scores <- melt(lmCV$metrics.dist[, .(R2, RE, CE)], 
               measure.vars = 1:3,
               variable.name = 'metric')
ggplot(scores) +
  geom_boxplot(aes(metric, value), fill = 'grey95') +
  geom_jitter(aes(metric, value, color = metric), width = 0.25) +
  scale_color_brewer(palette = 'Set2') +
  scale_x_discrete(
    name = NULL,
    guide = guide_prism_bracket(),
    labels = c('R\u00b2', 'RE', 'CE')) +
  scale_y_continuous(
    name = 'Metric value [-]',
    limits = c(0.25, 0.85),
    breaks = seq(0.2, 0.9, 0.1),
    minor_breaks = seq(0.25, 0.85, 0.05),
    guide = guide_prism_offset_minor()) +
  theme(legend.position = 'none')
```

Final reconstruction with bias correction.

```{r}
recFinal <- lmFit$rec
recFinal[, lp20 := dplR::pass.filt(Q, 20, 'low')]
recFinal[, lp50 := dplR::pass.filt(Q, 50, 'low')]
bcFit <- fitQmap(
  PSepAug$Qa,
  recFinal[year %in% PSepAug$year, Q],
  'RQUANT',
  wet.day = FALSE)
recFinal[, Qbc := doQmap(Q, bcFit)]
```

**Figure 3**

```{r, fig.width=7, fig.height=4}
p1 <- ggplot(recFinal[year %in% PSepAug$year]) +
  geom_line(aes(year, Q, colour = 'Reconstruction', linetype = 'Reconstruction'), size = 0.4) +
  geom_line(aes(year, Qbc, colour = 'Bias-corrected', linetype = 'Bias-corrected'), size = 0.4) +
  geom_line(aes(year, Qa, colour = 'Observation', linetype = 'Observation'), 
            PSepAug, size = 0.4) +
  scale_colour_manual(
    name = NULL, 
    values = c('Observation' = 'gray30',
               'Reconstruction' = 'steelblue',
               'Bias-corrected' = 'darkorange')) +
  scale_linetype_manual(
    name = NULL, 
    values = c('Observation' = 1, 
               'Reconstruction' = 1,
               'Bias-corrected' = 2)) +
  scale_x_continuous(
    minor_breaks = seq(1965, 2020, 5),
    limits = c(1966, 2018), 
    guide = guide_prism_offset_minor()) +
  scale_y_continuous(
    minor_breaks = seq(100, 850, 50),
    breaks = seq(100, 850, 100),
    limits = c(100, 850),
    guide = guide_prism_offset_minor()) +
  labs(x = 'Year', y = 'Precipitation [mm]', tag = 'a)') +
  theme(
    plot.margin = margin(r = 10),
    legend.box.margin = margin(),
    legend.margin = margin(),
    legend.key.width = unit(2, 'cm'),
    legend.position = 'top') 
p2 <- ggplot() +
  stat_density(aes(y = Q, colour = 'Reconstruction', linetype = 'Reconstruction'), 
               recFinal[year %in% PSepAug$year], geom = 'line', bw = 60, size = 0.4) +
  stat_density(aes(y = Qbc, colour = 'Bias-corrected', linetype = 'Bias-corrected'), 
               recFinal[year %in% PSepAug$year], geom = 'line', bw = 70, size = 0.4) +
  stat_density(aes(y = Qa, colour = 'Observation', linetype = 'Observation'), 
               PSepAug, geom = 'line', bw = 70, size = 0.4) +
  scale_colour_manual(
    name = NULL, 
    values = c('Observation' = 'gray30',
               'Reconstruction' = 'steelblue',
               'Bias-corrected' = 'darkorange')) +
  scale_linetype_manual(
    name = NULL, 
    values = c('Observation' = 1, 
               'Reconstruction' = 1,
               'Bias-corrected' = 2)) +
  scale_x_continuous(guide = guide_prism_offset()) +
  scale_y_continuous(
    minor_breaks = seq(100, 850, 50),
    breaks = seq(100, 850, 100),
    limits = c(100, 850),
    guide = guide_prism_offset_minor()) +
  labs(y = NULL, x = 'Density', colour = NULL, tag = 'b)') +
  theme(
    axis.text.y = element_blank(),
    legend.box.margin = margin(),
    legend.margin = margin(),
    legend.key.width = unit(2, 'cm'),
    legend.position = 'top') +
  theme(legend.position = 'top')

layout <- '
CCC
AAB
AAB
AAB
AAB
AAB
AAB
'

p1 + p2 + guide_area() +
  plot_layout(design = layout, guides = 'collect')
```

# Trend and drought analyses

## Quantile regression

```{r}
taus <- c(0.01, 0.05, 0.50, 0.95, 0.99)
# names(taus) <- sprintf(taus, fmt = '%.2f')
names(taus) <- paste0(taus * 100, '%')
qrDT <- taus |> 
  lapply(\(tau) {
    fit <- rq(Qbc ~ year, tau, recFinal)
    summ <- summary(fit, se = 'boot')
    data.table(
      b = fit$coefficients[1], 
      a = fit$coefficients[2], 
      signif = summ$coefficients[2, 4] < 0.1)
  }) |> 
  rbindlist(idcol = 'tau')
```

```{r}
recFinal[, period := fcase(
  year %in% 1517:1767, 1,
  default = 2)]
densCals <- recFinal[, {
  d <- density(Qbc, cut = 0, bw = 45)
  list(x = d$x, y = d$y / max(d$y) * 1.5)
}, by = period]

recFinal[, dP := Qbc - mean(Qbc),
       ][, type := classify_events(dP)
       ][, dp10 := pass.filt(dP, 10)]
droughts <- get_timing(recFinal$dP, recFinal$type)[type == 'drought']
droughts[, ':='(yearStart = recFinal[start, year],
                period = recFinal[start, period],
                yearFinal = recFinal[final, year])]
worstDroughts <- droughts[order(peak)][1:3]
```

**Figure 4**

```{r, fig.width=8, fig.height=6, message=FALSE}
reds <- brewer.pal(9, 'Reds')
p1 <- ggplot() +
  geom_line(aes(year, Qbc), recFinal, color = 'gray85') +
  geom_abline(aes(slope = a, intercept = b, color = tau, linetype = signif), 
              qrDT, show.legend = FALSE) +
  geom_line(aes(x = 1800, y = 400, color = tau, group = 1), qrDT) +
  geom_point(aes(year, Qbc), recFinal, color = 'gray50', shape = 1, size = 0.3) +
  scale_color_manual(
    values = c(reds[6], reds[3], 'darkviolet', blues9[4], blues9[7])) +
  scale_linetype_manual(values = c(2, 1)) +
  scale_x_continuous(
    expand = c(0, 0),
    minor_breaks = seq(1500, 2025, 25),
    limits = c(1500, 2025),
    guide = guide_prism_offset_minor()) +
  scale_y_continuous(
    breaks = seq(200, 800, 100),
    minor_breaks = seq(200, 850, 50),
    limits = c(200, 850),
    guide = guide_prism_offset_minor()) +
  guides(color = guide_legend(
    title.position = 'top',
    title.hjust = 0.5,
    override.aes = list(size = 0.6))) +
  labs(x = NULL, y = 'Precipitation [mm]', tag = 'a)') +
  theme(
    plot.tag.position = c(0, 0.94),
    legend.position = 'top')

p2 <- ggplot(densCals) +
  geom_ribbon(aes(x, ymin = period, ymax = y + period, group = factor(period),
                  fill = factor(period)), 
              alpha = 0.5) +
  geom_line(aes(x, y + period, group = factor(period))) +
  scale_fill_manual(
    labels = c('1517-1717', '1718-2018'),
    values = wesanderson::wes_palette('Royal1', 4)[c(1, 4)]) +
  scale_x_continuous(
    guide = guide_prism_offset_minor(),
    breaks = seq(200, 900, 100), 
    minor_breaks = seq(200, 900, 20)) +
  scale_y_continuous(
    expand = c(0, 0),
    breaks = c(1.5, 2.5),
    labels = c('1517-1717', '            1718-2018')) +
  labs(x = 'Precipitation [mm]', y = NULL, tag = 'b)') +
  theme(
    legend.position = 'none',
    axis.line.y = element_blank(),
    axis.text.y.left = element_text(angle = 90, hjust = 0.5),
    axis.ticks.y = element_blank())

p3 <- ggplot(recFinal) +
  geom_rect(
    aes(xmin = yearStart, xmax = yearFinal, ymin = -Inf, ymax = Inf),
    worstDroughts,
    fill = 'magenta4', alpha = 0.1) +
  geom_col(aes(year, dP, fill = dP), recFinal[dP > 0]) +
  scale_fill_distiller(
    palette = 'BrBG',
    limits = abs_range(recFinal$dP),
    # limits = abs_range(recFinal$dP),
    guide = guide_none(),
    direction = 1) +
  ggnewscale::new_scale_fill() +
  geom_col(aes(year, dP, fill = dP), recFinal[dP < 0]) +
  scale_fill_distiller(
    palette = 'BrBG', 
    # limits = abs_range(recFinal$dP),
    limits = c(min(recFinal$dP), -min(recFinal$dP)),
    guide = guide_none(),
    direction = 1) +
  scale_color_manual(values = 'black') +
  scale_x_continuous(
    expand = c(0, 0),
    minor_breaks = seq(1500, 2025, 25),
    limits = c(1500, 2025),
    guide = guide_prism_offset_minor()) +
  scale_y_continuous(
    breaks = seq(-300, 400, 100),
    minor_breaks = seq(-300, 400, 50),
    guide = guide_prism_offset_minor()) +
  labs(x = NULL, y = 'Precipitation anomaly [mm]', tag = 'c)') +
  theme(
    legend.position = 'top')

p4 <- ggplot(droughts) +
  geom_point(aes(dur, peak, color = as.character(period)), alpha = 0.6) +
  geom_text(
    aes(dur, peak, label = paste(yearStart, yearFinal, sep = '-')), 
    worstDroughts,
    size = 2, hjust = 0.8, nudge_y = -15) +
  scale_x_continuous(
    breaks = seq(2, 20, 2),
    guide = guide_prism_offset_minor()) +
  scale_y_continuous(
    breaks = seq(-300, 0, 50),
    limits = c(-275, -25),
    guide = guide_prism_offset_minor()) +
  scale_color_manual(
    labels = c('1517-1717', '1718-2018'),
    values = wesanderson::wes_palette('Royal1', 4)[c(1, 4)]) +
  labs(x = 'Drought duration [years]', y = 'Drought severity [mm]', tag = 'd)') +
  theme(legend.position = 'top')

p1 + p2 + p3 + p4 +
  plot_layout(widths = c(1.5, 1))
```

# Sensitivity analyses

We now rerun the reconstruction procedure with different SSS cut-off values.

## SSS = 0.75

```{r}
firstYear2 <- crn[sss > 0.75][, .N, by = year][N >= 4, first(year)]
lastYear <- 2018
crnWide2 <- crn[year %in% firstYear2:lastYear,
                dcast(.SD, year ~ site, value.var = c('res'))]
X2 <- crnWide2[, -'year'] |> 
  as.matrix() |> 
  imputePCA(ncp = 6) |> 
  {\(x) x$completeObs}() 
PC2 <- prcomp(X2, scale. = TRUE)$x
sv <- c(1, 2)
sPC2 <- as.data.table(PC2[, sv, drop = FALSE])
lmFit2 <- PCR_reconstruction(PSepAug, sPC2, firstYear2, transform = 'none')
set.seed(42)
lmCV2 <- cvPCR(PSepAug, sPC2, firstYear2, transform = 'none',
              Z = cvFolds, use.robust.mean = FALSE)
recFinal2 <- lmFit2$rec
bcFit2 <- fitQmap(
  PSepAug$Qa,
  recFinal2[year %in% PSepAug$year, Q],
  'RQUANT',
  wet.day = FALSE)
recFinal2[, Qbc := doQmap(Q, bcFit)]
```

## SSS = 0.85

```{r}
firstYear3 <- crn[sss > 0.85][, .N, by = year][N >= 4, first(year)]
lastYear <- 2018
crnWide3 <- crn[year %in% firstYear3:lastYear,
                dcast(.SD, year ~ site, value.var = c('res'))]
X3 <- crnWide3[, -'year'] |> 
  as.matrix() |> 
  imputePCA(ncp = 6) |> 
  {\(x) x$completeObs}() 
PC3 <- prcomp(X3, scale. = TRUE)$x
sv <- c(1, 2)
sPC3 <- as.data.table(PC3[, sv, drop = FALSE])
lmFit3 <- PCR_reconstruction(PSepAug, sPC3, firstYear3, transform = 'none')
lmCV3 <- cvPCR(PSepAug, sPC3, firstYear3, transform = 'none',
              Z = cvFolds, use.robust.mean = FALSE)
recFinal3 <- lmFit3$rec
bcFit3 <- fitQmap(
  PSepAug$Qa,
  recFinal3[year %in% PSepAug$year, Q],
  'RQUANT',
  wet.day = FALSE)
recFinal3[, Qbc := doQmap(Q, bcFit)]
```

## Compare reconstructins from different SSS

**Figure S3**

```{r, fig.width=7, fig.height=6}
p1 <- ggplot() +
  geom_line(aes(year, Qbc, color = 'SSS = 0.6', linetype = 'SSS = 0.6'), 
            recFinal, alpha = 1) +
  geom_line(aes(year, Qbc, color = 'SSS = 0.75', linetype = 'SSS = 0.75'), 
            recFinal2, alpha = 1) +
  # scale_color_manual(values = MetBrewer::met.brewer('Kandinsky', 2))
  scale_color_manual(name = NULL, 
                     values = wesanderson::wes_palettes$Cavalcanti1[2:1]) +
  scale_linetype_manual(name = NULL, values = c(1, 2)) +
  scale_x_continuous(
    expand = c(0, 0),
    minor_breaks = seq(1500, 2025, 25),
    limits = c(1500, 2025),
    guide = guide_prism_offset_minor()) +
  scale_y_continuous(
    breaks = seq(200, 800, 100),
    minor_breaks = seq(200, 850, 50),
    limits = c(200, 850),
    guide = guide_prism_offset_minor()) +
  labs(x = NULL, y = 'Precipitation [mm]') +
  theme(
    legend.key.height = unit(1, 'cm'),
    legend.position = 'top')
p2 <- ggplot() +
  geom_line(aes(year, Qbc, color = 'SSS = 0.6', linetype = 'SSS = 0.6'), 
            recFinal, alpha = 1) +
  geom_line(aes(year, Qbc, color = 'SSS = 0.85', linetype = 'SSS = 0.85'), 
            recFinal3, alpha = 1) +
  scale_color_manual(name = NULL, 
                     values = wesanderson::wes_palettes$Cavalcanti1[c(2,5)]) +
  scale_linetype_manual(name = NULL, values = c(1, 2)) +
  scale_x_continuous(
    expand = c(0, 0),
    minor_breaks = seq(1500, 2025, 25),
    limits = c(1500, 2025),
    guide = guide_prism_offset_minor()) +
  scale_y_continuous(
    breaks = seq(200, 800, 100),
    minor_breaks = seq(200, 850, 50),
    limits = c(200, 850),
    guide = guide_prism_offset_minor()) +
  labs(x = NULL, y = 'Precipitation [mm]') +
  theme(
    legend.key.height = unit(1, 'cm'),
    legend.position = 'top')
p1 / p2
```

## Compare skills

**Table S5**

```{r}
rbindlist(
  list(
    s1 = lmCV$metrics,
    s2 = lmCV2$metrics,
    s3 = lmCV3$metrics),
  idcol = 'cutoff')[, 1:4][, cutoff := c(0.6, 0.75, 0.85)][] |> 
  roundDT(3)
```

## Compare quantile regression results

```{r}
taus <- c(0.01, 0.05, 0.50, 0.95, 0.99)
names(taus) <- paste0(taus * 100, '%')
qrDT2 <- taus |> 
  lapply(\(tau) {
    fit <- rq(Qbc ~ year, tau, recFinal2)
    summ <- summary(fit, se = 'boot', R = 1000, bsmethod = 'wxy')
    data.table(
      b = fit$coefficients[1], 
      a = fit$coefficients[2], 
      signif = summ$coefficients[2, 4] < 0.1)
  }) |> 
  rbindlist(idcol = 'tau')

qrDT3 <- taus |> 
  lapply(\(tau) {
    fit <- rq(Qbc ~ year, tau, recFinal3)
    summ <- summary(fit, se = 'boot', R = 1000, bsmethod = 'wxy')
    data.table(
      b = fit$coefficients[1], 
      a = fit$coefficients[2], 
      signif = summ$coefficients[2, 4] < 0.1)
  }) |> 
  rbindlist(idcol = 'tau')
```

**Figure S4**

```{r, fig.width=7, fig.height=6}
p1 <- ggplot() +
  geom_line(aes(year, Qbc), recFinal2, color = 'gray85') +
  geom_abline(aes(slope = a, intercept = b, color = tau, linetype = signif), 
              qrDT2, show.legend = FALSE) +
  geom_line(aes(x = 1800, y = 400, color = tau, group = 1), qrDT2) +
  geom_point(aes(year, Qbc), recFinal2, color = 'gray50', shape = 1, size = 0.3) +
  scale_color_manual(
    values = c(reds[6], reds[3], 'darkviolet', blues9[4], blues9[7])) +
  scale_linetype_manual(values = c(2, 1)) +
  scale_x_continuous(
    expand = c(0, 0),
    minor_breaks = seq(1625, 2025, 25),
    limits = c(1625, 2025),
    guide = guide_prism_offset_minor()) +
  scale_y_continuous(
    breaks = seq(200, 800, 100),
    minor_breaks = seq(200, 850, 50),
    limits = c(200, 850),
    guide = guide_prism_offset_minor()) +
  guides(color = guide_legend(
    title.position = 'top',
    title.hjust = 0.5,
    override.aes = list(size = 0.6))) +
  labs(x = NULL, y = 'Precipitation [mm]') +
  theme(
    legend.position = 'top')
p2 <- ggplot() +
  geom_line(aes(year, Qbc), recFinal3, color = 'gray85') +
  geom_abline(aes(slope = a, intercept = b, color = tau, linetype = signif), 
              qrDT3, show.legend = FALSE) +
  geom_line(aes(x = 1800, y = 400, color = tau, group = 1), qrDT3) +
  geom_point(aes(year, Qbc), recFinal3, color = 'gray50', shape = 1, size = 0.3) +
  scale_color_manual(
    values = c(reds[6], reds[3], 'darkviolet', blues9[4], blues9[7])) +
  scale_linetype_manual(values = c(2, 1)) +
  scale_x_continuous(
    expand = c(0, 0),
    minor_breaks = seq(1625, 2025, 25),
    limits = c(1625, 2025),
    guide = guide_prism_offset_minor()) +
  scale_y_continuous(
    breaks = seq(200, 800, 100),
    minor_breaks = seq(200, 850, 50),
    limits = c(200, 850),
    guide = guide_prism_offset_minor()) +
  guides(color = guide_legend(
    title.position = 'top',
    title.hjust = 0.5,
    override.aes = list(size = 0.6))) +
  labs(x = NULL, y = 'Precipitation [mm]') +
  theme(
    legend.position = 'top')
p1 / p2
```

## Using ARSTAN chronologies

```{r}
crnWide4 <- crn[year %in% firstYear:lastYear,
                dcast(.SD, year ~ site, value.var = c('ars'))]
X4 <- crnWide4[, -'year'] |> 
  as.matrix() |> 
  imputePCA(ncp = 6) |> 
  {\(x) x$completeObs}() 
PC4 <- prcomp(X4, scale. = TRUE)$x
sv <- c(1, 2)
sPC4 <- as.data.table(PC4[, sv, drop = FALSE])
lmFit4 <- PCR_reconstruction(PSepAug, sPC4, firstYear, transform = 'none')
lmCV4 <- cvPCR(PSepAug, sPC4, firstYear, transform = 'none',
              Z = cvFolds, use.robust.mean = FALSE)
recFinal4 <- lmFit4$rec
bcFit4 <- fitQmap(
  PSepAug$Qa,
  recFinal4[year %in% PSepAug$year, Q],
  'RQUANT',
  wet.day = FALSE)
recFinal4[, Qbc := doQmap(Q, bcFit)]
```

**Table S3**

```{r}
rbindlist(list(
  residual = lmCV$metrics,
  ARSTAN = lmCV4$metrics),
  idcol = 'Chronology')[, 1:4] |> 
  roundDT(3)
```

Export the reconstruction.

```{r export, eval=FALSE}
recFinal[, period := NULL]
fwrite(
  recFinal[, .(year, P = Q, Plower = Ql, Pupper = Qu, Pbc = Qbc)], 
  'results/chitral-sep-aug-precip-reconst.csv')
```