#############################################################################
# This script computes summary statistics reported in the manuscript. 
# These statistics are derived from the same underlying data presented in the
#  manuscript.
#############################################################################
# Baseline set ----
rm(list = ls())
library(tidyverse)

Locations <- c("Dubai", "Rio_de_Janeiro", "Rome", "Toronto")
delta_climate = -0.04
rho_k <- 0.1
sigma_beta_val = 0
rho_mean = 0.5
School.Terms = F
years = 10 
import = 1e-5

# Locations full names
location_labels <- c(
  "Rome" = "Rome, Italy",
  "Rio_de_Janeiro" = "Rio de Janeiro, Brazil",
  "Toronto" = "Toronto, Canada",
  "Dubai" = "Dubai, UAE"
)

# Variables full names
variable_labels <- c(
    "Te_est" = "Temperature",
    "RH_est" = "Relative Humidity"
)

# Models full names
model_labels <- c(CC_stand_smooth_AC = "Test",
                  CC_true = "Control")

df_main_fig <- NULL
for(roiLoc in Locations){
    data_Path <- paste("Results_Aggregated/R_sims_Te",
                      delta_climate, "_RH", delta_climate, 
                      "_",
                      "sigma", sigma_beta_val, "_",
                      paste("rho_mean", rho_mean, "_", sep=""),
                      roiLoc, "_",
                      sprintf("_rho_k%s", gsub("\\.", "-", as.character(rho_k))),
                      sprintf("_import_%s", as.character(import)),
                      ifelse(School.Terms == T, "_School-Forced", ""),
                      sep = "")
    final_pars <- readRDS(paste(data_Path, "/final_pars.rds", sep = ""))
    sims_reg_tot <- readRDS(paste(data_Path, "/sims_reg_tot.rds", sep = "")) |> 
      dplyr::filter(type %in% c("CC_true", "CC_stand_smooth_AC")) |> 
      dplyr::mutate(location = roiLoc) |>
      dplyr::select(-c(e_RH_low_CI, e_RH_high_CI,
                      log_SI_lag_VAR, mean_log_SI_lag,
                      R2, e_Te_low_CI, e_Te_high_CI,
                      TP_Te, TP_RH)) |> 
      dplyr::group_by(type, .id, location) |> 
      dplyr::summarise(
        MAB_Te = mean(abs(e_Te - unique(final_pars[,"True.e_Te"]))), # Computing the mean absolute bias for temperature
        MAB_RH = mean(abs(e_RH - unique(final_pars[,"True.e_RH"]))),
        sd_AB_Te = sd(abs(e_Te - unique(final_pars[,"True.e_Te"]))), # Computing the standard deviation in absolute biases for temperature
        sd_AB_RH = sd(abs(e_RH - unique(final_pars[,"True.e_RH"]))),
        MSE_Te = mean(e_Te_se),
        MSE_RH = mean(e_RH_se),
        SDSE_Te = sd(e_Te_se),
        SDSE_RH = sd(e_RH_se),
        MPOW_Te = mean(e_Te_pow),
        MPOW_RH = mean(e_RH_pow)
      ) |> dplyr::ungroup() |> 
      dplyr::mutate(true_val = if(unique(final_pars[,"True.e_Te"]) == unique(final_pars[,"True.e_RH"])){
                                unique(final_pars[,"True.e_Te"])      
                              })

    df_main_fig <- df_main_fig |>
      bind_rows(sims_reg_tot); 
    rm(sims_reg_tot)
}

df_main_fig |> dplyr::filter(type == "CC_stand_smooth_AC") |> dplyr::group_by(location) |> 
    dplyr::summarise(mean_MPOW_Te = mean(MPOW_Te) * 100,
                    mean_MPOW_RH = mean(MPOW_RH) * 100)

# Differences in rho_k ----
rm(list=ls())

Locations <- c("Dubai", "Rio_de_Janeiro", "Rome", "Toronto")
delta_climate = -0.04
Rho_ks <- c(0.1, 0.16)
sigma_beta_val = 0
Conds <- expand.grid(loc = Locations, rho_k = Rho_ks)
rho_mean = 0.5
School.Terms = F
years = 10
import = 1e-5

# Initialise the dataframe to store all
df_main = NULL
# For loop to generate big cumulative dataframe
for(i in 1:nrow(Conds)){
    data_Path = paste("Results_Aggregated/R_sims_Te",
                      delta_climate, "_RH", delta_climate,
                      "_",
                      "sigma", sigma_beta_val, "_",
                      paste("rho_mean", rho_mean, "_", sep=""),
                      Conds[i, "loc"], "_",
                      ifelse(School.Terms == T, "_School-Forced", ""),
                      sprintf("_rho_k%s", gsub("\\.", "-", as.character(Conds[i, "rho_k"]))),
                      sprintf("_import_%s", as.character(import)),
                      sep = "")
    R0s = readRDS(paste(data_Path, "/final_pars.rds", sep = "")) |> dplyr::select(.id, R0) |> distinct()
    sims_reg_tot = readRDS(paste(data_Path, "/sims_reg_tot.rds", sep = "")) |> 
      dplyr::filter(type == "CC_stand_smooth_AC") |> 
      dplyr::mutate(location = Conds[i, "loc"],
                    rho_k = Conds[i, "rho_k"]) |>
      dplyr::select(-c(e_RH_low_CI, e_RH_high_CI,
                      log_SI_lag_VAR, mean_log_SI_lag,
                      R2, e_Te_low_CI, e_Te_high_CI,
                      TP_Te, e_Te_pow, TP_RH, e_RH_pow)) |> 
      dplyr::group_by(.id, location, rho_k, type) |> 
      dplyr::full_join(
                      y = R0s, 
                      by = c(".id")
                      )

    df_main = df_main |>
      bind_rows(sims_reg_tot); 
    rm(sims_reg_tot)
}

df_diff <- df_main |> mutate(bias_Te = abs(e_Te - delta_climate), 
                             bias_RH = abs(e_RH - delta_climate)) |> 
                  dplyr::group_by(location, .id, rho_k, type, .groups = "keep") |> 
                  dplyr::summarise(MAB_id_Te = mean(bias_Te), 
                            MAB_id_RH = mean(bias_RH),
                            MSE_id_Te = mean(e_Te_se),
                            MSE_id_RH = mean(e_RH_se)) |> 
                            ungroup() |> 
                  mutate(RMAB_id_Te = MAB_id_Te / abs(delta_climate) * 100,
                         RMAB_id_RH = MAB_id_RH / abs(delta_climate) * 100,
                         RMSE_id_Te = MSE_id_Te / abs(delta_climate) * 100,
                         RMSE_id_RH = MSE_id_RH / abs(delta_climate) * 100) |> 
                  select(-c(MAB_id_Te, MAB_id_RH, MSE_id_Te, MSE_id_RH)) |> 
        pivot_wider(
          names_from = c(rho_k),
          values_from = c(RMAB_id_Te, RMAB_id_RH, RMSE_id_Te, RMSE_id_RH)
  ) |> 
  mutate( # New_Sensitivity - reference -> diff > 0 -> ref is better
    diff_RMAB_id_RH = `RMAB_id_RH_0.16` - `RMAB_id_RH_0.1`,
    diff_RMAB_id_Te = `RMAB_id_Te_0.16` - `RMAB_id_Te_0.1`,
    diff_RMSE_id_RH = `RMSE_id_RH_0.16` - `RMSE_id_RH_0.1`,
    diff_RMSE_id_Te = `RMSE_id_Te_0.16` - `RMSE_id_Te_0.1`
  ) |> 
  select(location, type, .id, diff_RMAB_id_RH, diff_RMAB_id_Te, diff_RMSE_id_RH, diff_RMSE_id_Te) |> 
  pivot_longer(cols = c(diff_RMAB_id_RH, diff_RMAB_id_Te, diff_RMSE_id_RH, diff_RMSE_id_Te),
               names_to= "variable", values_to= 'vals')

df_scatter <- df_diff |> dplyr::filter(type == "CC_stand_smooth_AC") |> 
  separate(variable, into = c("metric", "var_type"), sep = "_id_") |>
  pivot_wider(
    names_from = metric,
    values_from = vals
  ) |>   mutate(var_type = factor(var_type, levels = c("Te", "RH")))

df_scatter |> dplyr::group_by(var_type) |> 
    dplyr::summarise(mean_diff_RMAB = mean(diff_RMAB),
                    sd_diff_RMAB = sd(diff_RMAB),
                    mean_diff_RMSE = mean(diff_RMSE),
                    sd_diff_RMSE = sd(diff_RMSE))

# Differneces in delta_climate ----
rm(list=ls())

Locations <- c("Dubai", "Rio_de_Janeiro", "Rome", "Toronto")
Delta_climates = c(-0.04, -.02)
Conds <- expand.grid(loc = Locations, delta_climate = Delta_climates)
sigma_beta_val = 0
rho_mean = 0.5
School.Terms = F
years = 10 
Pois <- F
rho_k <- 0.1
import = 1e-5
roi_model <- "CC_stand_smooth_AC"

df_effects = NULL
# For loop to generate big cumulative dataframe 
for(i in 1:nrow(Conds)){
    data_Path = paste("Results_Aggregated/R_sims_Te",
                      Conds[i, "delta_climate"], "_RH", Conds[i, "delta_climate"], 
                      "_",
                      "sigma", sigma_beta_val, "_",
                      paste("rho_mean", rho_mean, "_", sep=""),
                      Conds[i, "loc"], "_",
                      ifelse(School.Terms == T, "_School-Forced", ""),
                      sprintf("_rho_k%s", gsub("\\.", "-", as.character(rho_k))),
                      sprintf("_import_%s", as.character(import)),
                      sep = "")
    climate_effects <- readRDS(paste(data_Path, "/final_pars.rds", sep = "")) |> dplyr::select(True.e_Te, True.e_RH) |> distinct()
    sims_reg_tot <- readRDS(paste(data_Path, "/sims_reg_tot.rds", sep = "")) |>  filter(type == roi_model) |> 
      dplyr::mutate(location = Conds[i, "loc"],
             True_e_Te = climate_effects$True.e_Te,
             True_e_RH = climate_effects$True.e_RH) |>
      dplyr::ungroup() |> 
      dplyr::select(c(location, rep, .id, True_e_RH, True_e_Te, e_Te, e_RH, e_Te_se, e_RH_se))

    df_effects <- df_effects |>
      bind_rows(sims_reg_tot); 
    rm(sims_reg_tot)
}

df_diff <- df_effects |> mutate(bias_Te = abs(e_Te - True_e_Te), 
                                bias_RH = abs(e_RH - True_e_RH)) |> 
                  group_by(location, .id, True_e_Te, True_e_RH) |> 
                  dplyr::summarise(MAB_id_Te = mean(bias_Te), 
                            MAB_id_RH = mean(bias_RH),
                            MSE_id_Te = mean(e_Te_se),
                            MSE_id_RH = mean(e_RH_se)) |> 
                            ungroup() |> 
                  mutate(RMAB_id_Te = MAB_id_Te / abs(True_e_Te) * 100,
                         RMAB_id_RH = MAB_id_RH / abs(True_e_RH) * 100,
                         RMSE_id_Te = MSE_id_Te / abs(True_e_Te) * 100,
                         RMSE_id_RH = MSE_id_RH / abs(True_e_RH) * 100) |> 
                  select(-c(MAB_id_Te, MAB_id_RH, MSE_id_Te, MSE_id_RH)) |> 
        pivot_wider(
          names_from = c(True_e_Te, True_e_RH),
          values_from = c(RMAB_id_Te, RMAB_id_RH, RMSE_id_Te, RMSE_id_RH)
  ) |> 
  mutate( # New_Sensitivity - reference -> diff > 0 -> ref is better
    diff_RMAB_id_RH = `RMAB_id_RH_-0.02_-0.02` - `RMAB_id_RH_-0.04_-0.04`,
    diff_RMAB_id_Te = `RMAB_id_Te_-0.02_-0.02` - `RMAB_id_Te_-0.04_-0.04`,
    diff_RMSE_id_RH = `RMSE_id_RH_-0.02_-0.02` - `RMSE_id_RH_-0.04_-0.04`,
    diff_RMSE_id_Te = `RMSE_id_Te_-0.02_-0.02` - `RMSE_id_Te_-0.04_-0.04`
  ) |> select(location, .id, diff_RMAB_id_RH, diff_RMAB_id_Te, diff_RMSE_id_RH, diff_RMSE_id_Te) |> 
  pivot_longer(cols = c(diff_RMAB_id_RH, diff_RMAB_id_Te, diff_RMSE_id_RH, diff_RMSE_id_Te),
               names_to = "variable", values_to= 'vals')

  df_scatter <- df_diff |>
  separate(variable, into = c("metric", "var_type"), sep = "_id_") |>
  pivot_wider(
    names_from = metric,
    values_from = vals
  ) |>   mutate(var_type = factor(var_type, levels = c("Te", "RH")))

df_scatter |> dplyr::group_by(var_type) |> 
    dplyr::summarise(mean_diff_RMAB = mean(diff_RMAB),
                    sd_diff_RMAB = sd(diff_RMAB),
                    mean_diff_RMSE = mean(diff_RMSE),
                    sd_diff_RMSE = sd(diff_RMSE))

# Differences in sigma_beta ----
rm(list=ls())

Locations <- c("Dubai", "Rio_de_Janeiro", "Rome", "Toronto")
delta_climate = -0.04
rho_k <- 0.1
Sigma_beta_vals = c(0, 0.05)
Conds <- expand.grid(loc = Locations, sigma_beta = Sigma_beta_vals)
rho_mean = 0.5
School.Terms = F
years = 10 
import = 1e-5

# Initialise the dataframe to store all
df_main = NULL
# For loop to generate big cumulative dataframe 
for(i in 1:nrow(Conds)){
    data_Path = paste("Results_Aggregated/R_sims_Te",
                      delta_climate, "_RH", delta_climate, 
                      "_",
                      "sigma", Conds[i, "sigma_beta"], "_",
                      paste("rho_mean", rho_mean, "_", sep=""),
                      Conds[i, "loc"], "_",
                      ifelse(School.Terms == T, "_School-Forced", ""),
                      sprintf("_rho_k%s", gsub("\\.", "-", as.character(rho_k))),
                      sprintf("_import_%s", as.character(import)),
                      sep = "")
    R0s = readRDS(paste(data_Path, "/final_pars.rds", sep = "")) |> dplyr::select(.id, R0) |> distinct()
    sims_reg_tot = readRDS(paste(data_Path, "/sims_reg_tot.rds", sep = "")) |> 
    dplyr::filter(type == "CC_stand_smooth_AC") |> 
      dplyr::mutate(location = Conds[i, "loc"],
                    sigma_beta = Conds[i, "sigma_beta"]) |>
      dplyr::select(-c(e_RH_low_CI, e_RH_high_CI,
                      log_SI_lag_VAR, mean_log_SI_lag,
                      R2, e_Te_low_CI, e_Te_high_CI,
                      TP_Te, e_Te_pow, TP_RH, e_RH_pow)) |> 
      dplyr::group_by(.id, location, sigma_beta) |> 
      dplyr::full_join(
                      y = R0s, 
                      by = c(".id")
                      )

    df_main = df_main |>
      bind_rows(sims_reg_tot); 
    rm(sims_reg_tot)
}

df_diff <- df_main |> mutate(bias_Te = abs(e_Te - delta_climate), 
                             bias_RH = abs(e_RH - delta_climate)) |> 
                  dplyr::group_by(location, .id, sigma_beta, .groups = "keep") |> 
                  dplyr::summarise(MAB_id_Te = mean(bias_Te), 
                            MAB_id_RH = mean(bias_RH),
                            MSE_id_Te = mean(e_Te_se),
                            MSE_id_RH = mean(e_RH_se)) |> 
                            ungroup() |> 
                  mutate(RMAB_id_Te = MAB_id_Te / abs(delta_climate) * 100,
                         RMAB_id_RH = MAB_id_RH / abs(delta_climate) * 100,
                         RMSE_id_Te = MSE_id_Te / abs(delta_climate) * 100,
                         RMSE_id_RH = MSE_id_RH / abs(delta_climate) * 100) |> 
                  select(-c(MAB_id_Te, MAB_id_RH, MSE_id_Te, MSE_id_RH)) |> 
        pivot_wider(
          names_from = c(sigma_beta),
          values_from = c(RMAB_id_Te, RMAB_id_RH, RMSE_id_Te, RMSE_id_RH)
  ) |> 
  mutate( # New_Sensitivity - reference -> diff > 0 -> ref is better
    diff_RMAB_id_RH = `RMAB_id_RH_0.05` - `RMAB_id_RH_0`,
    diff_RMAB_id_Te = `RMAB_id_Te_0.05` - `RMAB_id_Te_0`,
    diff_RMSE_id_RH = `RMSE_id_RH_0.05` - `RMSE_id_RH_0`,
    diff_RMSE_id_Te = `RMSE_id_Te_0.05` - `RMSE_id_Te_0`
  ) |> select(location, .id, diff_RMAB_id_RH, diff_RMAB_id_Te, diff_RMSE_id_RH, diff_RMSE_id_Te) |> 
  pivot_longer(cols = c(diff_RMAB_id_RH, diff_RMAB_id_Te, diff_RMSE_id_RH, diff_RMSE_id_Te),
               names_to= "variable", values_to= 'vals')

df_scatter <- df_diff |>
  tidyr::separate(variable, into = c("metric", "var_type"), sep = "_id_") |>
  tidyr::pivot_wider(
    names_from = metric,
    values_from = vals
  ) |> dplyr::mutate(var_type = factor(var_type, levels = c("Te", "RH")))

df_scatter |> dplyr::group_by(var_type) |> 
    dplyr::summarise(mean_diff_RMAB = mean(diff_RMAB),
                    sd_diff_RMAB = sd(diff_RMAB),
                    mean_diff_RMSE = mean(diff_RMSE),
                    sd_diff_RMSE = sd(diff_RMSE))

# Differences in Term-time ----
rm(list=ls())

Locations <- c("Dubai", "Rio_de_Janeiro", "Rome", "Toronto")
delta_climate = -0.04
rho_k <- 0.1
sigma_beta_val = 0
School.Terms = c(T,F)
Conds <- expand.grid(loc = Locations, term_time = School.Terms)
rho_mean = 0.5
years = 10 
import = 1e-5

# Initialise the dataframe to store all
df_main = NULL
# For loop to generate big cumulative dataframe 
for(i in 1:nrow(Conds)){
    data_Path = paste("Results_Aggregated/R_sims_Te",
                      delta_climate, "_RH", delta_climate, 
                      "_",
                      "sigma", sigma_beta_val, "_",
                      paste("rho_mean", rho_mean, "_", sep=""),
                      Conds[i, "loc"], "_",
                      sprintf("_rho_k%s", gsub("\\.", "-", as.character(rho_k))),
                      sprintf("_import_%s", as.character(import)),
                      ifelse(Conds[i, "term_time"] == T, "_School-Forced", ""),
                      sep = "")
    sims_reg_tot = readRDS(paste(data_Path, "/sims_reg_tot.rds", sep = "")) |> 
    dplyr::filter(type == "CC_stand_smooth_AC") |> 
      dplyr::mutate(location = Conds[i, "loc"],
                    term_time = Conds[i, "term_time"]) |>
      dplyr::ungroup() |> 
      dplyr::select(-c(e_RH_low_CI, e_RH_high_CI, type,
                      log_SI_lag_VAR, mean_log_SI_lag,
                      R2, e_Te_low_CI, e_Te_high_CI,
                      TP_Te, e_Te_pow, TP_RH, e_RH_pow))

    df_main = df_main |>
      bind_rows(sims_reg_tot); 
    rm(sims_reg_tot)
}

final_pars <- readRDS(paste(data_Path, "/final_pars.rds", sep = ""))
True_e_Te <- unique(final_pars$True.e_Te)
True_e_RH <- unique(final_pars$True.e_RH)

df_diff <- df_main |> mutate(bias_Te = abs(e_Te - delta_climate), 
                  bias_RH = abs(e_RH - delta_climate)) |> 
                  select(-c(e_Te, e_RH)) |> 
                  group_by(location, .id, term_time) |> 
                  dplyr::summarise(MAB_id_Te = mean(bias_Te), 
                            MAB_id_RH = mean(bias_RH),
                            MSE_id_Te = mean(e_Te_se),
                            MSE_id_RH = mean(e_RH_se)) |> 
                            ungroup() |> 
                  mutate(RMAB_id_Te = MAB_id_Te / abs(True_e_Te) * 100,
                         RMAB_id_RH = MAB_id_RH / abs(True_e_RH) * 100,
                         RMSE_id_Te = MSE_id_Te / abs(True_e_Te) * 100,
                         RMSE_id_RH = MSE_id_RH / abs(True_e_RH) * 100) |> 
                  select(-c(MAB_id_Te, MAB_id_RH, MSE_id_Te, MSE_id_RH)) |> 
        pivot_wider(
          names_from = c(term_time),
          values_from = c(RMAB_id_Te, RMAB_id_RH, RMSE_id_Te, RMSE_id_RH)) |> 
  mutate(
    diff_RMAB_id_RH = `RMAB_id_RH_TRUE` - `RMAB_id_RH_FALSE`,
    diff_RMAB_id_Te = `RMAB_id_Te_TRUE` - `RMAB_id_Te_FALSE`,
    diff_RMSE_id_RH = `RMSE_id_RH_TRUE` - `RMSE_id_RH_FALSE`,
    diff_RMSE_id_Te = `RMSE_id_Te_TRUE` - `RMSE_id_Te_FALSE`
  ) |> select(location, .id, diff_RMAB_id_RH, diff_RMAB_id_Te, diff_RMSE_id_RH, diff_RMSE_id_Te) |> 
  pivot_longer(cols = c(diff_RMAB_id_RH, diff_RMAB_id_Te, diff_RMSE_id_RH, diff_RMSE_id_Te),
               names_to= "variable", values_to= 'vals')

  df_scatter <- df_diff |>
  separate(variable, into = c("metric", "var_type"), sep = "_id_") |>
  pivot_wider(
    names_from = metric,
    values_from = vals
  ) |>   mutate(var_type = factor(var_type, levels = c("Te", "RH"))) 
  
df_scatter |> dplyr::group_by(var_type) |> 
    dplyr::summarise(mean_diff_RMAB = mean(diff_RMAB),
                    sd_diff_RMAB = sd(diff_RMAB),
                    mean_diff_RMSE = mean(diff_RMSE),
                    sd_diff_RMSE = sd(diff_RMSE))

# Differences in performance with 52 and 75 knots -----
rm(list=ls())

Locations <- c("Dubai", "Rio_de_Janeiro", "Rome", "Toronto")
delta_climate = -0.04
Rho_ks <- 0.1
sigma_beta_val = 0
Knots = 52
Conds <- expand.grid(loc = Locations, rho_k = Rho_ks, k_max = Knots)
rho_mean = 0.5
School.Terms = F
years = 10
import = 1e-5

# Initialise the dataframe to store all
df_52 = NULL
# For loop to generate big cumulative dataframe
for(i in 1:nrow(Conds)){
        data_Path = paste("Results_Aggregated/R_sims_Te",
                      delta_climate, "_RH", delta_climate,
                      "_",
                      "sigma", sigma_beta_val, "_",
                      paste("rho_mean", rho_mean, "_", sep=""),
                      Conds[i, "loc"], "_",
                      sprintf("_rho_k%s", gsub("\\.", "-", as.character(Conds[i, "rho_k"]))),
                      sprintf("_import_%s", as.character(import)),
                      ifelse(School.Terms == T, "_School-Forced", ""),
                      sep = "")
    sims_reg_tot <- readRDS(paste(data_Path, "/sims_reg_tot.rds", sep = "")) |> 
      dplyr::filter(type == "CC_stand_smooth_AC") |> 
      dplyr::mutate(loc = Conds[i, "loc"],
                    rho_k = Conds[i, "rho_k"],
                    k_max = Conds[i, "k_max"]) |>
        dplyr::ungroup() |> 
      dplyr::select(loc, rep, .id, k_max, e_Te, e_RH, e_Te_se, e_RH_se)
    df_52 = df_52 |>
      bind_rows(sims_reg_tot); 
    rm(sims_reg_tot)
}

df_75 <- readRDS("Results_Aggregated/Regressions_75Knots.rds") |> dplyr::select(-c(k_index, edf))

df_main <- rbind(df_75,df_52)

df_diff <- df_main |> mutate(bias_Te = abs(e_Te - delta_climate), 
                             bias_RH = abs(e_RH - delta_climate)) |> 
                  dplyr::group_by(loc, .id, k_max, .groups = "keep") |> 
                  dplyr::summarise(MAB_id_Te = mean(bias_Te), 
                            MAB_id_RH = mean(bias_RH),
                            MSE_id_Te = mean(e_Te_se),
                            MSE_id_RH = mean(e_RH_se)) |> 
                            ungroup() |> 
                  mutate(RMAB_id_Te = MAB_id_Te / abs(delta_climate) * 100,
                         RMAB_id_RH = MAB_id_RH / abs(delta_climate) * 100,
                         RMSE_id_Te = MSE_id_Te / abs(delta_climate) * 100,
                         RMSE_id_RH = MSE_id_RH / abs(delta_climate) * 100) |> 
                  select(-c(MAB_id_Te, MAB_id_RH, MSE_id_Te, MSE_id_RH)) |> 
        pivot_wider(
          names_from = c(k_max),
          values_from = c(RMAB_id_Te, RMAB_id_RH, RMSE_id_Te, RMSE_id_RH)
  ) |> 
  mutate( # New_Sensitivity - reference -> diff > 0 -> ref is better
    diff_RMAB_id_RH = `RMAB_id_RH_75` - `RMAB_id_RH_52`,
    diff_RMAB_id_Te = `RMAB_id_Te_75` - `RMAB_id_Te_52`,
    diff_RMSE_id_RH = `RMSE_id_RH_75` - `RMSE_id_RH_52`,
    diff_RMSE_id_Te = `RMSE_id_Te_75` - `RMSE_id_Te_52`
  ) |> 
  select(loc, .id, diff_RMAB_id_RH, diff_RMAB_id_Te, diff_RMSE_id_RH, diff_RMSE_id_Te) |> 
  pivot_longer(cols = c(diff_RMAB_id_RH, diff_RMAB_id_Te, diff_RMSE_id_RH, diff_RMSE_id_Te),
               names_to= "variable", values_to= 'vals')

  df_diff |> 
  separate(variable, into = c("metric", "var_type"), sep = "_id_") |>
  pivot_wider(
    names_from = metric,
    values_from = vals
  ) |> mutate(var_type = factor(var_type, levels = c("Te", "RH"))) |> dplyr::group_by(var_type) |> 
    dplyr::summarise(mean_diff_RMAB = mean(diff_RMAB),
                    sd_diff_RMAB = sd(diff_RMAB),
                    mean_diff_RMSE = mean(diff_RMSE),
                    sd_diff_RMSE = sd(diff_RMSE))

# Additional lagged effect of weahter ----
file_path <- "Results/R_sims_Te-0.016_Te_Lag-0.027_RH-0.016_RH_Lag-0.027_sigma0_rho_mean0.5_Toronto__rho_k0-1_import_1e-05_LAGGED/sims_reg_tot.rds"
delta_Now <- -0.016
delta_Lag <- -0.027

sims_reg_tot_lag <- readRDS(file = file_path) |>
      dplyr::filter(type == "CC_stand_smooth_AC") |> 
      dplyr::mutate(type = "CC_stand_smooth_AC_lag") |> 
      dplyr::select(-c(e_RH_low_CI, e_RH_high_CI,
                      log_SI_lag_VAR, mean_log_SI_lag,
                      R2, e_Te_low_CI, e_Te_high_CI,
                      TP_Te, TP_RH)) |> 
      dplyr::group_by(type, .id, location) |> 
      dplyr::summarise(
        # First for the t-1 effect estimates
        MAB_Te = mean(abs(e_Te - delta_Now)),
        MAB_RH = mean(abs(e_RH - delta_Now)),
        sd_AB_Te = sd(abs(e_Te - delta_Now)),
        sd_AB_RH = sd(abs(e_RH - delta_Now)),
        MSE_Te = mean(e_Te_se),
        MSE_RH = mean(e_RH_se),
        MA_e_Te = mean(abs(e_Te)),
        MA_e_RH = mean(abs(e_RH)),
        SDSE_Te = sd(e_Te_se),
        SDSE_RH = sd(e_RH_se),
        MPOW_Te = mean(e_Te_pow),
        MPOW_RH = mean(e_RH_pow),
        # Now for the t-2 effect estimates
        MAB_Te_Lag = mean(abs(e_Te_Lag - delta_Lag)),
        MAB_RH_Lag = mean(abs(e_RH_Lag - delta_Lag)),
        sd_AB_Te_Lag = sd(abs(e_Te_Lag - delta_Lag)),
        sd_AB_RH_Lag = sd(abs(e_RH_Lag - delta_Lag)),
        MSE_Te_Lag = mean(e_Te_se_Lag),
        MSE_RH_Lag = mean(e_RH_se_Lag),
        MA_e_Te_Lag = mean(abs(e_Te_Lag)),
        MA_e_RH_Lag = mean(abs(e_RH_Lag)),
        SDSE_Te_Lag = sd(e_Te_se_Lag),
        SDSE_RH_Lag = sd(e_RH_se_Lag),
        MPOW_Te_Lag = mean(e_Te_pow_Lag),
        MPOW_RH_Lag = mean(e_RH_pow_Lag)
      ) |> dplyr::ungroup() |> 
      dplyr::mutate(true_val_Now = delta_Now,
                    true_val_Lag = delta_Lag)


summ_lag <- sims_reg_tot_lag |> dplyr::filter(type == "CC_stand_smooth_AC_lag") |> 
    dplyr::mutate(RMAB_Te_Now = MAB_Te / abs(true_val_Now) * 100,
        RMAB_RH_Now = MAB_RH / abs(true_val_Now) * 100,
        RMSE_Te_Now = MSE_Te / abs(true_val_Now) * 100,
        RMSE_RH_Now = MSE_RH / abs(true_val_Now) * 100,
        RMAB_Te_Lag = MAB_Te_Lag / abs(true_val_Lag) * 100,
        RMAB_RH_Lag = MAB_RH_Lag / abs(true_val_Lag) * 100,
        RMSE_Te_Lag = MSE_Te_Lag / abs(true_val_Lag) * 100,
        RMSE_RH_Lag = MSE_RH_Lag / abs(true_val_Lag) * 100,
        MPOW_Te = MPOW_Te * 100,
        MPOW_Te_Lag = MPOW_Te_Lag * 100,
        MPOW_RH = MPOW_RH * 100,
        MPOW_RH_Lag = MPOW_RH_Lag * 100
    ) |> dplyr::select(MPOW_Te,MPOW_Te_Lag,MPOW_RH,MPOW_RH_Lag,RMAB_Te_Now,
                        RMAB_RH_Now,RMSE_Te_Now,RMSE_RH_Now,RMAB_Te_Lag,RMAB_RH_Lag,RMSE_Te_Lag,RMSE_RH_Lag)

range_table <- data.frame(t(sapply(summ_lag, range, na.rm = TRUE)))
colnames(range_table) <- c("min","max")

range_table |>
  dplyr::mutate(across(everything(), ~ round(.x, 0)))