Exp-1 Met Results

1 Objective

Run ED using different climate to create an ensemble of runs

• How does that impact ED results?
• What impact does it make on dimensions of stability?
• Is there a relationship between the inputs & metrics?
• How do we want to calculate recovery?
  • When it reaches the control? When it reaches a steady state?

2 Compare Inputs

# Import the meterolical inputs. 
RESULTS_DIR <- file.path(BASE_DIR, "ED-outputs", "met_results")
# Create a data frame of all the meterolgical data used. 
list.files(RESULTS_DIR, pattern = 'yr', full.names = TRUE) %>%  
  lapply(function(x){
    
    name <- gsub(x = basename(x), pattern = '_yr.csv', replacement = '')
    data <- read.csv(x, stringsAsFactors = FALSE)
    data$scn <- name
    return(data)
    
  }) %>%  
  dplyr::bind_rows() %>%  
  as.data.table -> 
  data
# A mapping file of variable names, description, and units. 
# This will be used to add more information to the met data data frame. 
var_info <- data.table(variable = c('dlwrf', 'nbdsf', 'nddsf', 'vbdsf',
                                    'vddsf', 'prate', 'pres', 'hgt', 
                                    'ugrd', 'vgrd', 'sh', 'tmp'), 
                       description = c('downward long wave radiation', 
                                       'near infrared beam downward solar radiation', 
                                       'near IR diffuse downward solar radiation', 
                                       'visible beam downward solar radiation',
                                       'visible diffuse downward solar radiation', 
                                       'precipitation rate', 
                                       'atmospheric pressure', 
                                       'geopotential height', 
                                       'zonal wind',
                                       'meridional wind', 
                                       'specific humidity',
                                       'air temperature'),  
                       units = c('W m-2', 'W m-2', 'W m-2', 'W m-2',
                                 'W m-2', 'kgH2O m-2 s-1', 'Pa', 'm',
                                 'm s-1', 'm s-', 'kgH2O kgAir-1', 'K'), stringsAsFactors=FALSE)
data <- data[var_info, on = 'variable']
data <- na.omit(data)
data$label <- paste0(data$description, ' (', data$units, ')')

# Define function that generates the standard met input plot 
# Args 
#   data: data table of the met data to plot
#   yr1: the first year of met data to plot, default set to 1950
#   yr2: the final year of met data to plot, default set to 2030 
# Return: a faceted line plot of the met data, ggplot2 object
met_plot <- function(d, yr1 = 1950, yr2 = 2040){
  
  # Color the background 
  ALPHA <- 0.05
  BCK_GRND1 <- '#dcffcc'
  BCK_GRND2 <- '#9fdfcd'
  BCK_GRND3 <- '#baabda'
  
  d[year >= yr1 & year <= yr2] %>%  
    ggplot() + 
    geom_rect(aes(xmin=yr1, xmax = 1979,
                  ymin = -Inf, ymax = Inf), fill = BCK_GRND1, alpha = ALPHA) +
    geom_rect(aes(xmin=1979, xmax = 2010,
                  ymin = -Inf, ymax = Inf), fill = BCK_GRND2, alpha = ALPHA) +
    geom_rect(aes(xmin=2010, xmax = yr2,
                  ymin = -Inf, ymax = Inf), fill = BCK_GRND3, alpha = ALPHA) + 
    labs(caption = 'colored background: looped over obs., obs., 5 year obs. mean ') + 
    geom_line(aes(year, value, linetype = scn)) + 
    facet_wrap('label', scales = 'free', strip.position = "top",
             labeller = label_wrap_gen(width = 20, multi_line = TRUE)) + 
    THEME + 
    labs(y = NULL, x = 'Year', title = 'Annual Mean Meteorological Inputs') + 
    theme(legend.position = 'none') -> 
    plot 
  
  return(plot)
}
# Split the met data up into two different plots, this will make the plots easier to  look at. 
var_list <- unique(data$variable)
pannels_n1 <- ceiling(length(var_list)/2)

Up until 2010 the runs use the same meteorology inputs, starting in 2010, about a decade before the disturbance we start to use different meteorology inputs based on 5 year means sampled from the observational data set.

2.1 Plot Climate Input

# Plot the first half of the met input values 
met_plot(data[variable %in% var_list[1:pannels_n1]])

# Plot the second half of the met input values 
met_plot(data[variable %in% var_list[pannels_n1+1:length(var_list)]])

The inputs for the two different met records starting in 2010, so that that ED has a steady climate for the ~10 yrs before the disturbance, the met ensemble captures a wide range of different climates.

3 ED Results

annual_data <- as.data.table(read.csv(file.path(BASE_DIR, "ED-outputs", "results",
                                                "exp-1-yr.csv"), stringsAsFactors = FALSE))

3.1 Plot Baseline/Control

annual_data$label <- paste0(annual_data$description, " (", annual_data$unit, ")")
annual_data[scn == 'harvest_0_1day_above', ] %>%  
  .[ , severity := '0'] %>% 
  ggplot() + 
      geom_rect(aes(xmin=1900, xmax = 1979,
                  ymin = -Inf, ymax = Inf), fill = BCK_GRND1, alpha = ALPHA) +
    geom_rect(aes(xmin=1979, xmax = 2010,
                  ymin = -Inf, ymax = Inf), fill = BCK_GRND2, alpha = ALPHA) +
    geom_rect(aes(xmin=2010, xmax = 2048,
                  ymin = -Inf, ymax = Inf), fill = BCK_GRND3, alpha = ALPHA) + 
    geom_vline(xintercept = 2019, color = 'black') + 
  geom_line(aes(year, value, linetype = met)) + 
 facet_wrap('label', scales = 'free',  
            labeller = label_wrap_gen(width = 50, multi_line = TRUE)) + 
  THEME + 
  theme(legend.position = 'none')+
  labs(title = 'Disturbance Severity 0%', 
       y = NULL, 
       x = "Year", 
       caption = 'vertical line at 2019, when the FoRTE treatment is applied\ncolored background: looped over obs., obs., 5 year obs. mean ') 

Yes changing the meteorology input files does make a difference in the baseline scenarios, as we can tell based on the vertical line by the time the FoRTE treatment is applied the ED forest stands diverged from one another.

4 Disturbance

What is the difference between the different treatments and the control (0%) runs?

baseline <- annual_data[scn == 'harvest_0_1day_above', ][ , .(year, variable, description, unit, met, baseline = value)]
disturbance <- annual_data[scn != 'harvest_0_1day_above', ]
setnames(x = disturbance, old = "value", new = "disturbance", skip_absent = TRUE)
dt <- baseline[disturbance, on = c( "year", "variable", "description", "unit", "met")]
dt[ , value := disturbance - baseline]
dt[ , severity := gsub(pattern = "harvest_|_1day_above", replacement = "", x = scn)]

4.1 Plot Disturbance

dt$`Disturbance Severity` <- paste0(dt$severity, ' %')
dt[year >= 2015, ] %>% 
  ggplot(aes(year, value, color = `Disturbance Severity`, linetype = met)) + 
  geom_vline(xintercept = 2019, color = 'black') + 
  geom_line() + 
   facet_wrap('label', scales = 'free',  
            labeller = label_wrap_gen(width = 30, multi_line = TRUE)) + 
  THEME + 
  labs(y = 'Disturbance - Baseline', 
       x = 'Year', 
       caption = 'vertical line 2019 disturbance event')  + 
  scale_color_manual(values = FORTE_SEVERITY_COLORS) + 
  guides(linetype = FALSE)

The over shoot in the met2 track for the carbon fluxes is really interesting! Do we think that would happen if the met1 run continued longer?

5 Dimensions of Stability

# Find the recovery at specific time places. 
# Args 
#   d: a data frame of the annual ED outputs
#   x: a vectory of the years to calculate the recovery at. 
# Return: a data frame of the recovery values
recovery_at_X <- function(d, x){
  
# Check inputs 
req_names <- c('scn', 'year', 'variable', 'description', 'unit', 'value', 'met')
assert_that(has_name(d, req_names))
baseline_scn <- "harvest_0_1day_above" 
assert_that(baseline_scn %in% d$scn)
assert_that(all(x %in% d$year))
# Separate the control or baseline scneario values from the output from the other 
# treatment groups. 
baseline_values <- d[scn == baseline_scn][ , c("year", "variable", "value", "met")]
names(baseline_values) <- c("year", "variable", "baseline", "met")
treatment_values <- d[scn != baseline_scn][ , c("scn", "year", "variable", "description", 
                                                "unit", "treatment" = "value", "met", "label")]
names(treatment_values) <- c("scn", "year", "variable", "description", 
                                                "unit", "treatment", "met", "label")
# Combine the control and the treatment values as a wide df, then take the # ln(treatment/control) per year / variable / met realization. 
disturbance_d <- treatment_values[baseline_values, on = c("year", "variable", "met")][year %in% 2019:2050]
disturbance_d[ , value := log(treatment/baseline)]
# Now that the ln(treatment/control) has been calculated we can use these values to calculate
# the various dimensions of ecosystem stability. 
recovery_d <- disturbance_d[year %in% x]
recovery_d <- recovery_d[ , c("scn", "year", "variable", "description", 
                              "unit", "met",  "value")]
names(recovery_d) <- c("scn", "year", "variable", "description", 
                       "unit", "met",  "recovery")
dt <- recovery_d
# Add the disturbance severity column. 
dt$`Disturbance Severity` <- paste0(gsub(x = dt$scn,
                                         pattern = "harvest_|_1day_above", 
                                         replacement = ""), " %")
# Format and apply factor information to the variable name so the dimensions of stability plot 
# step through the carbon cycle in order GPP -> NPP -> NEP. 
dt$variable <- gsub(pattern = "_patch", replacement = "", x = dt$variable)
dt$variable <- factor(x = dt$variable, levels = c("GPP", "NPP", "NEP"), ordered = TRUE)
dt <- na.omit(dt)
return(dt)
}
recovery_at_X_data <- recovery_at_X(annual_data, c(2025, 2030, 2035, 2040, 2045))

# Calculate the various dimensiosn of stability 
# Args
#   d: data table of the annual ED results 
# Return: data table of the dimensions of stability. 
calculate_dimensions_stability <- function(d){
  
# Check inputs 
req_names <- c('scn', 'year', 'variable', 'description', 'unit', 'value', 'met')
assert_that(has_name(d, req_names))
baseline_scn <- "harvest_0_1day_above" 
assert_that(baseline_scn %in% d$scn)
# Separate the control or baseline scneario values from the output from the other 
# treatment groups. 
baseline_values <- d[scn == baseline_scn][ , c("year", "variable", "value", "met")]
names(baseline_values) <- c("year", "variable", "baseline", "met")
treatment_values <- d[scn != baseline_scn][ , c("scn", "year", "variable", "description", 
                                                "unit", "treatment" = "value", "met", "label")]
names(treatment_values) <- c("scn", "year", "variable", "description", 
                                                "unit", "treatment", "met", "label")
# Combine the control and the treatment values as a wide df, then take the # ln(treatment/control) per year / variable / met realization. 
disturbance_d <- treatment_values[baseline_values, on = c("year", "variable", "met")][year %in% 2019:2050]
disturbance_d[ , value := log(treatment/baseline)]
# Now that the ln(treatment/control) has been calculated we can use these values to calculate
# the various dimensions of ecosystem stability. 
treatment_yr <- 2019 # The year of the girdel 
# Resistance: ${r}_{t} = ln(\frac{{F}_{disturbance}}{{F}_{control}})$ at the year year 
# of disturabnce (2019). 
resitance_d <- disturbance_d[year == treatment_yr][ , c("scn", "variable",
                                                        "description", "unit", "met",  "value")]
names(resitance_d) <- c("scn", "variable", "description", "unit", "met",  "resitance")
# Resilience : ${r}_{i} = ln(\frac{{F}_{disturbance}}{{F}_{control}}) = i + {r}_{i} * t$
# Calculate the resilience (the slope of the recovery starting at the treatment_yr)
d_to_use <- disturbance_d[year >= 2019 & year <= 2030] # Do not include data from pre treatment years. 
# Get the result from the linear regression per variable / scn / met ensemble member. 
d_list <- split(d_to_use, interaction(d_to_use$variable, d_to_use$scn, d_to_use$met), 
                drop = TRUE) 
fits <- lapply(d_list, function(input){lm(value ~ year, data = input)})
# Get the slope from the fits, this is the resilience value. 
resilience_values <- sapply(fits, function(x){x$coefficients[2]})
resilience_d <- data.table(info = names(resilience_values), resilience = resilience_values)
resilience_d <- resilience_d[, c("variable", "scn", "met", "x") := tstrsplit(info, ".",
                                                                             fixed=TRUE)]
resilience_d <- resilience_d[, c("resilience", "variable", "scn", "met")]
# Now using the fits calculate the temporal stability. 
# ${s}_{t} = \frac{1}{ SE of resid}$
se_resid <- sapply(fits, function(model){sqrt(deviance(model)/df.residual(model))})
temp_stability_d <- data.table(info = names(se_resid), temp_stability = 1/se_resid)
temp_stability_d <- temp_stability_d[, c("variable", "scn", "met") := tstrsplit(info, ".",
                                                                             fixed=TRUE)]
temp_stability_d <- temp_stability_d[, c("temp_stability", "variable", "scn", "met")]
# Finally calculate the recovery, for now we will define recovery at year 2039
# two decades after the year of the girdel. 
recovery_d <- disturbance_d[year == 2030]
recovery_d <- recovery_d[ , c("scn", "variable", "description", "unit", "met",  "value")]
names(recovery_d) <- c("scn", "variable", "description", "unit", "met",  "recovery")
recovery_d$recovery <- recovery_d$recovery
# Combine all dimensions of stability into a single data table. 
resitance_d[resilience_d, on = c("variable", "scn", "met")] %>% 
  .[temp_stability_d, on = c("variable", "scn", "met")] %>% 
  .[recovery_d, on = c("variable", "scn", "met")] %>% 
  .[, c( "scn", "variable", "description", "unit", "met", "resitance", 
         "resilience", "temp_stability", "recovery")] -> 
  dt 
# Add the disturbance severity column. 
dt$`Disturbance Severity` <- paste0(gsub(x = dt$scn,
                                         pattern = "harvest_|_1day_above", 
                                         replacement = ""), " %")
# Format and apply factor information to the variable name so the dimensions of stability plot 
# step through the carbon cycle in order GPP -> NPP -> NEP. 
dt$variable <- gsub(pattern = "_patch", replacement = "", x = dt$variable)
dt$variable <- factor(x = dt$variable, levels = c("GPP", "NPP", "NEP"), ordered = TRUE)
dt <- na.omit(dt)
return(dt)
}
stability_metrics <- calculate_dimensions_stability(annual_data)

5.1 Plot Resistance

\({r}_{t} = ln(\frac{{F}_{disturbance}}{{F}_{control}})\)

When \[{t}\] is the year of the disturbance, in this case 2019.

stability_metrics %>% 
  ggplot(aes(`Disturbance Severity`, resitance, color = `Disturbance Severity`)) + 
  geom_point(size = 2, alpha = 0.4) + 
  facet_wrap('variable') + 
  THEME  + 
  labs(title = 'Resistance', y = '2019 ln(F disturabnce / F baseline)') + 
  scale_color_manual(values = FORTE_SEVERITY_COLORS) 

5.2 Plot Resilience

\({r}_{i} = ln(\frac{{F}_{disturbance}}{{F}_{control}}) = i + {r}_{i} * t\)

stability_metrics %>%  
ggplot() + 
  geom_point(aes(`Disturbance Severity`, resilience, color = `Disturbance Severity`),
             size = 2, alpha = 0.4) + 
  facet_wrap('variable') + 
  THEME  + 
  labs(title = 'Resilience', y = 'Slope ln(dist/control)') + 
  scale_color_manual(values = FORTE_SEVERITY_COLORS)

5.3 Plot Temporal Stability

\({s}_{t} = \frac{1}{ SE of resid}\)

stability_metrics %>%  
ggplot() + 
  geom_point(aes(`Disturbance Severity`, temp_stability, color = `Disturbance Severity`),
             size = 2, alpha = 0.4) + 
  facet_wrap('variable') + 
  THEME  + 
  labs(title = 'Temporal Stability', 
       y = ' 1 / (SE of resid)') + 
  scale_color_manual(values = FORTE_SEVERITY_COLORS)

5.4 Plot Recovery

Take a look at the recovery at separate time steps.

ggplot(data = recovery_at_X_data[ , year := as.character(year)], 
       aes(x = year, y = recovery, fill = `Disturbance Severity`)) + 
  geom_boxplot() + 
  THEME + 
  labs(x = 'Year', y = "ln(disturbance/baseline)", title = "Recovery at Year") + 
  facet_wrap("variable", ncol = 1, scales = "free")

Recovery 2 decades after the treatment.

stability_metrics %>%  
ggplot() + 
  geom_point(aes(`Disturbance Severity`, recovery, color = `Disturbance Severity`),
             size = 2, alpha = 0.4) + 
  facet_wrap('variable') + 
  THEME  + 
  labs(title = 'Recovery', 
       y = 'Recovery at 2030') + 
  scale_color_manual(values = FORTE_SEVERITY_COLORS)

6 Input & Metric Interactions

# Create a data frame of the disturbance period annual met inputs, note that it should be 
# a wide data frame in preperation for the pair wise plots. 
data[year == 2019][ , c("variable", "value", "scn")] %>%
  unique() %>% 
  setnames(old = c("scn", "value", "variable"), 
           new = c("met", "met_vavlue", "met_variable")) -> 
  dist_met_data
# Combine the met and the dimensions of stability data into a single data frame. 
stability_metrics[ , c("scn", "variable", "met", "resitance", 
                       "resilience", "temp_stability", "recovery", "Disturbance Severity")] %>%
  melt(id.vars = c("scn", "variable", "met", "Disturbance Severity"),
       measure.vars = c("resitance", "resilience", "temp_stability", "recovery"), 
       variable.name = "metric_name", value.name = "metric_value") %>% 
  unique() -> 
  metric_long
full_df <- dist_met_data[metric_long, on = "met", allow.cartesian=TRUE]

var_list <- unique(full_df$variable)
metric   <- unique(full_df$metric_name)
full_df  <- full_df[var_info, on = .(met_variable = variable)] 
lapply(metric, function(m){
  #  & met_variable %in% c("prate", "tmp", "sh")
  full_df[metric_name == m] %>% 
  ggplot(aes(met_vavlue, metric_value, color = `Disturbance Severity`)) + 
  geom_point() + 
  facet_grid(variable ~ description, scales = "free",  labeller = label_wrap_gen(width = 20, multi_line = TRUE)) + 
  THEME + 
  labs(y = NULL, x = NULL, title = m) + 
  theme(legend.position = "bottom") + 
  scale_color_manual(values = FORTE_SEVERITY_COLORS)
}) -> 
  plots
names(plots) <- metric

6.1 Plot Input x Resitance

plots$resitance + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=0.5))

6.2 Plot Input x Resilience

plots$resilience + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=0.5))

6.3 Plot Input x Temporal Stability

plots$temp_stability + 
  labs(title = "Temporal Stability") + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=0.5))

6.4 Plot Input x Recovery

plots$recovery + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=0.5))

7 Dynamic Recovery

# Calculate the ln(disturbance/control)
# Args 
#   d: a data frame of the anual data 
# Return: a data frame of the ln(disturbance/control) per year/variable/met/disturbance severity 
calc_ln_disturbance <- function(d){
  
  # Check inputs 
  req_names <- c('scn', 'year', 'variable', 'description', 'unit', 'value', 'met')
  assert_that(has_name(d, req_names))
  
  baseline_scn <- "harvest_0_1day_above" 
  assert_that(baseline_scn %in% d$scn)
  
  # Separate the control or baseline scneario values from the output from the other 
  # treatment groups. 
  baseline_values <- d[scn == baseline_scn][ , c("year", "variable", "value", "met")]
  names(baseline_values) <- c("year", "variable", "baseline", "met")
  treatment_values <- d[scn != baseline_scn][ , c("scn", "year", "variable", "description", 
                                                "unit", "treatment" = "value", "met")]
  names(treatment_values) <- c("scn", "year", "variable", "description", "unit", "treatment", "met")
  
  # Combine the control and the treatment values as a wide df, then take the # ln(treatment/control) per year / variable / met
  # realization. 
  disturbance_d <- treatment_values[baseline_values, on = c("year", "variable", "met")][year %in% 2019:2050]
  disturbance_d <- disturbance_d[ , value := log(treatment/baseline)]
  disturbance_d <- disturbance_d[ , disturbance_severity := paste0(gsub(x = scn, pattern = "harvest_|_1day_above", replacement = ""), " %")]
  disturbance_d <- disturbance_d[ , variable :=  gsub(pattern = "_patch", replacement = "", x = variable)]
  disturbance_d$variable <- factor(x = disturbance_d$variable, levels = c("GPP", "NPP", "NEP"), ordered = TRUE)
  disturbance_d <- disturbance_d[ , .(year, variable, met, value, disturbance_severity)]
  disturbance_d <- na.omit(disturbance_d)
  disturbance_d <- disturbance_d[year %in% 2019:2050]
  
  return(disturbance_d)
}
disturbance <- calc_ln_disturbance(annual_data)
# Define the functions to be used by frollapply to get the fits from the rolling slopes 
slope_func <- function(value){
  
  assert_that(is.vector(value))
  sq <- seq(1:length(value))
  fit <- lm(value~sq)
  
  slope <- coef(fit)[2]
  names(slope) <- "slope"
  return(slope)
  
}
intercept_func <- function(value){
  
  assert_that(is.vector(value))
  sq <- seq(1:length(value))
  fit <- lm(value~sq)
  
  inter <- coef(fit)[1]
  names(inter) <- "y-intercept"
  return(inter)
  
}
# Apply the slope and intercept functions to each variable / met / severity disturbance values. 
split(disturbance,
      interaction(disturbance$variable, disturbance$met, disturbance$disturbance_severity)) %>% 
  lapply(function(x){
    
    window <- 5
    slope <- frollapply(x = x[["value"]], n = window, FUN = slope_func, align = "left")
    slope2 <- frollapply(x = slope, n = window, FUN = slope_func, align = "left")
    
    out <- cbind(x, slope = slope, slope2 = slope2)
    return(out)
    
  }) %>% 
  rbindlist %>%  
  na.omit -> 
  dist_slope_y

7.1 Plot of Disturbance Rate of Change

This is the first derivative, the rate of change of disturabnce, we were intrested in when this is 0, if that would be the recovery point becacuse that would refelect when the ecosystem has reached a new steady state…

ggplot(data = dist_slope_y[year %in% 2019:2050]) +
 geom_hline(aes(yintercept = 0)) +
 geom_line(aes(year, slope, group = interaction(disturbance_severity, met)), color = "grey", alpha = 0.3) +
 geom_point(aes(year, slope, color = disturbance_severity, group = interaction(disturbance_severity, met))) +
 facet_wrap("variable", scales = "free", ncol = 1) +
 theme_bw() +
 labs(title = "5 year running slope")

# Selet the recovery years as the years where the mangnitude of the rate of change is less than 1e-3 or set a 
# hard limit of 2030, which is unideal but for now it is not not to avoid some sort of arbitriaty cut off. 
recovery_yr <- dist_slope_y[abs(slope) <= 1e-3][ , .(recovery_yr = min(year)), by = .(variable, met, disturbance_severity)] 
recovery_yr <- recovery_yr[recovery_yr <= 2030]
# Now join with the disurbance values. 
recovery_yr[disturbance, on = c("variable", "met", "disturbance_severity"), nomatch = NA] %>% 
  .[ , recovery_yr:= ifelse(test = is.na(recovery_yr), yes = 2030, no = recovery_yr)] %>% 
  .[year <= recovery_yr] -> 
  disturbance_recovery_subset
# Find recovery value aka the disturbance value at the time of recovery. 
recovery_values <- disturbance_recovery_subset[ year == recovery_yr]

Using when the rate of change or running slope is less than or equal to 0 or with a cut off there are only a handful of runs that recovery or reach a new stead state.

recovery_yr[disturbance, on = c("variable", "met", "disturbance_severity"), nomatch = NA] %>% 
 # .[ , recovery_yr:= ifelse(test = is.na(recovery_yr), yes = 2030, no = recovery_yr)] %>% 
  .[year <= recovery_yr] %>% 
  ggplot(aes(year, value, color = disturbance_severity, group = interaction(met, disturbance_severity))) + 
  geom_line() + 
  facet_grid(variable~disturbance_severity, scales = "free") + 
  labs(title = "Dynamic Recovery", y = "ln(treatment/control)") + 
  theme_bw() + 
  theme(legend.position = "none")

What if the cut off value is when the rate off change is less than or equal to 0 or at most 2030? If that is the case then our results look something like.

recovery_yr[disturbance, on = c("variable", "met", "disturbance_severity"), nomatch = NA] %>%
  .[ , recovery_yr:= ifelse(test = is.na(recovery_yr), yes = 2030, no = recovery_yr)] %>% 
  .[year <= recovery_yr] %>% 
  ggplot(aes(year, value, color = disturbance_severity, group = interaction(met, disturbance_severity))) + 
  geom_line() + 
  facet_grid(variable~disturbance_severity, scales = "free") + 
  labs(title = "Dynamic Recovery", y = "ln(treatment/control)") + 
  theme_bw() + 
  theme(legend.position = "none")

Is a linear fit the right fit here? How do we deal with these different recovery dates?