#############################################################################
# This script produces the tables present in the Supplementary material.
#############################################################################
# Tables with the RMAB and RMSE for test model ----
rm(list=ls())
library(tidyverse)
library(kableExtra)
library(webshot)
## First import all necessary dataframes ----
Locations <- c("Dubai", "Rio_de_Janeiro", "Rome", "Toronto")
Delta_climates = c(-0.04, -.02)
Sigma_beta_vals = c(0, .05)
rho_mean = 0.5
School.Terms = c(T,F)
years = 10
Rho_ks <- c(.1, .16)
import = 1e-5
roi_model <- "CC_stand_smooth_AC"
Conds <- expand.grid(location = Locations, delta_climate = Delta_climates,
rho_k = Rho_ks, Term_time = School.Terms,
sigma_beta = Sigma_beta_vals)
model_labels <- c(CC_AC = "Autocorrelation",
CC_stand_smooth = "Time-smooth",
CC_stand_smooth_AC = "Time-smooth + Autocorrelation",
CC_true = "Control model")
cumulative_df = NULL
# For loop to generate big cumulative dataframe
for(i in 1:nrow(Conds)){
if(Conds[i,"delta_climate"] == -.02 & (Conds[i,"Term_time"] == T | Conds[i,"sigma_beta"] != 0)){next}
data_Path = paste("Results_Aggregated/R_sims_Te",
Conds[i, "delta_climate"], "_RH", Conds[i, "delta_climate"],
"_",
"sigma", Conds[i, "sigma_beta"], "_",
ifelse(rho_mean != 0, paste("rho_mean", rho_mean, "_", sep="") , ""),
Conds[i, "location"],"_",
sprintf("_rho_k%s", gsub("\\.", "-", as.character(Conds[i, "rho_k"]))),
sprintf("_import_%s", as.character(import)),
ifelse(Conds[i, "Term_time"] == T, "_School-Forced", ""),
sep = "")
final_pars = readRDS(paste(data_Path, "/final_pars.rds", sep = ""))
tmp_df = readRDS(paste(data_Path, "/sims_reg_tot.rds", sep = "")) |> filter(type == roi_model) |>
dplyr::ungroup() |>
dplyr::select(-type) |>
dplyr::mutate(.id = factor(.id),
rep = factor(rep)) |>
dplyr::rename(AB_Te = bias_Te, AB_RH = bias_RH) |>
dplyr::group_by(.id) |> dplyr::mutate(MAB_Te_id = mean(AB_Te),
MAB_RH_id = mean(AB_RH),
MSE_Te_id = mean(e_Te_se),
MSE_RH_id = mean(e_RH_se)
) |>
dplyr::ungroup() |>
dplyr::mutate(MAB_Te = mean(AB_Te),
MAB_RH = mean(AB_RH),
MSE_Te = mean(e_Te_se),
MSE_RH = mean(e_RH_se)) |>
dplyr::mutate(bias_v_id_Te = mean((MAB_Te_id - MAB_Te) ** 2),
bias_v_id_RH = mean((MAB_RH_id - MAB_RH) ** 2),
SE_v_id_Te = mean((MSE_Te_id - MSE_Te) ** 2),
SE_v_id_RH = mean((MSE_RH_id - MSE_RH) ** 2),
) |>
dplyr::mutate(bias_v_rep_Te = mean((AB_Te - MAB_Te_id) ** 2),
bias_v_rep_RH = mean((AB_RH - MAB_RH_id) ** 2),
SE_v_rep_Te = mean((e_Te_se - MSE_Te_id) ** 2),
SE_v_rep_RH = mean((e_RH_se - MSE_RH_id) ** 2)
) |>
dplyr::ungroup() |> dplyr::select(-c("rep",".id", "e_Te", "e_Te_se", "e_Te_low_CI",
"e_Te_high_CI", "e_RH", "e_RH_se", "e_RH_low_CI", "e_RH_high_CI",
"log_SI_lag_VAR", "mean_log_SI_lag", "R2", "AB_Te", "AB_RH",
"TP_Te", "e_Te_pow", "TP_RH", "e_RH_pow", "MAB_Te_id", "MAB_RH_id",
"MAB_Te_id", "MAB_RH_id", "MSE_Te_id", "MSE_RH_id")) |>
dplyr::distinct() |>
dplyr::mutate(
across(
.cols = c(bias_v_id_Te, bias_v_id_RH, bias_v_rep_Te, bias_v_rep_RH),
.fns = ~ sqrt(.) / abs(Conds[i,"delta_climate"]) * 100),
across(
.cols = c(MAB_Te, MAB_RH, MSE_Te, MSE_RH),
.fns = ~ . / abs(Conds[i,"delta_climate"]) * 100),
across(
.cols = c(SE_v_id_Te, SE_v_id_RH, SE_v_rep_Te, SE_v_rep_RH),
.fns = ~ sqrt(.) / abs(Conds[i,"delta_climate"]) * 100)
) |>
dplyr::rename(RMAB_Te = MAB_Te, RMAB_RH = MAB_RH, RMSE_Te = MSE_Te, RMSE_RH = MSE_RH) |>
dplyr::mutate(
across(matches("^(RMAB|bias)"), ~ round(.x, 0)),
across(matches("^(RMSE|SE)"), ~ round(.x, 0))
) |> cbind(Conds[i,])
cumulative_df <- rbind(cumulative_df, tmp_df)
}
clean_df <- cumulative_df |> dplyr::mutate(
location = recode(location,
Rome = "Rome, Italy",
Rio_de_Janeiro = "Rio de Janeiro, Brazil",
Toronto = "Toronto, Canada",
Dubai = "Dubai, UAE"
)
) |>
dplyr::mutate(
formatted_RMAB_Te = sprintf("%.f (%.f) [%.f]", RMAB_Te, bias_v_id_Te, bias_v_rep_Te),
formatted_RMAB_RH = sprintf("%.f (%.f) [%.f]", RMAB_RH, bias_v_id_RH, bias_v_rep_RH),
formatted_RMSE_Te = sprintf("%.f (%.f) [%.f]", RMSE_Te, SE_v_id_Te, SE_v_rep_Te),
formatted_RMSE_RH = sprintf("%.f (%.f) [%.f]", RMSE_RH, SE_v_id_RH, SE_v_rep_RH),
) |>
dplyr::select(-c("RMAB_Te", "RMAB_RH",
"RMSE_Te", "RMSE_RH",
"bias_v_id_Te", "bias_v_id_RH",
"SE_v_id_Te", "SE_v_id_RH",
"bias_v_rep_Te", "bias_v_rep_RH",
"SE_v_rep_Te", "SE_v_rep_RH")) |>
mutate(
Temperature = paste0(formatted_RMAB_Te, "<br/>", formatted_RMSE_Te),
`Relative Humidity` = paste0(formatted_RMAB_RH, "<br/>", formatted_RMSE_RH)
)
## Now create each single table -----
baseline <- clean_df |> dplyr::rename(Location = location) |> dplyr::filter(rho_k == 0.1,
delta_climate == -0.04,
sigma_beta == 0,
Term_time == F) |> dplyr::select(-c("delta_climate", "rho_k", "Term_time", "sigma_beta",
"formatted_RMAB_Te", "formatted_RMAB_RH",
"formatted_RMSE_Te", "formatted_RMSE_RH")) |>
kable(format = "html", escape = FALSE, caption = roi_model) |>
# caption = "Baseline parameter values") |>
column_spec(1:3, latex_valign = "m") |>
kable_styling(bootstrap_options = c("bordered"), html_font = "Cambria") |>
row_spec(0:4, align = "c") |>
save_kable(file="Figures/E_table.png", density = 1800, zoom = 1.5)
overdisp <- clean_df |> dplyr::rename(Location = location) |> dplyr::filter(rho_k == 0.16,
delta_climate == -0.04,
sigma_beta == 0,
Term_time == F) |> dplyr::select(-c("delta_climate", "rho_k", "Term_time", "sigma_beta",
"formatted_RMAB_Te", "formatted_RMAB_RH",
"formatted_RMSE_Te", "formatted_RMSE_RH")) |>
kable(format = "html", escape = FALSE,
caption = "With overdispersion set at 16%") |>
column_spec(1:3, latex_valign = "m") |>
kable_styling(bootstrap_options = c("bordered"), html_font = "Cambria") |>
row_spec(0:4, align = "c") |>
save_kable(file="Figures/F_table.png", density = 1800, zoom = 1.5)
moody <- clean_df |> dplyr::rename(Location = location) |> dplyr::filter(rho_k == 0.1,
delta_climate == -0.02,
sigma_beta == 0,
Term_time == F) |> dplyr::select(-c("delta_climate", "rho_k", "Term_time", "sigma_beta",
"formatted_RMAB_Te", "formatted_RMAB_RH",
"formatted_RMSE_Te", "formatted_RMSE_RH")) |>
kable(format = "html", escape = FALSE,
caption = "With stronger effect of climate") |>
column_spec(1:3, latex_valign = "m") |>
kable_styling(bootstrap_options = c("bordered"), html_font = "Cambria") |>
row_spec(0:4, align = "c") |>
save_kable(file="Figures/G_table.png", density = 1800, zoom = 1.5)
stochastic <- clean_df |> dplyr::rename(Location = location) |> dplyr::filter(rho_k == 0.1,
delta_climate == -0.04,
sigma_beta == 0.05,
Term_time == F) |> dplyr::select(-c("delta_climate", "rho_k", "Term_time", "sigma_beta",
"formatted_RMAB_Te", "formatted_RMAB_RH",
"formatted_RMSE_Te", "formatted_RMSE_RH")) |>
kable(format = "html", escape = FALSE,
caption = "With process noise set to 5%") |>
column_spec(1:3, latex_valign = "m") |>
kable_styling(bootstrap_options = c("bordered"), html_font = "Cambria") |>
row_spec(0:4, align = "c") |>
save_kable(file="Figures/H_table.png", density = 1800, zoom = 1.5)
forced <- clean_df |> dplyr::rename(Location = location) |> dplyr::filter(rho_k == 0.1,
delta_climate == -0.04,
sigma_beta == 0,
Term_time == T) |> dplyr::select(-c("delta_climate", "rho_k", "Term_time", "sigma_beta",
"formatted_RMAB_Te", "formatted_RMAB_RH",
"formatted_RMSE_Te", "formatted_RMSE_RH")) |>
kable(format = "html", escape = FALSE,
caption = "With Term-time forcing") |>
column_spec(1:3, latex_valign = "m") |>
kable_styling(bootstrap_options = c("bordered"), html_font = "Cambria") |>
row_spec(0:4, align = "c") |>
save_kable(file="Figures/I_table.png", density = 1800, zoom = 1.5)
# Table for baseline with the control model ----
rm(list=ls())
Locations <- c("Dubai", "Rio_de_Janeiro", "Rome", "Toronto")
Delta_climates = -0.04
Sigma_beta_vals = 0
rho_mean = 0.5
School.Terms = F
years = 10
Rho_ks <- .1
import = 1e-5
roi_model <- "CC_true"
Conds <- expand.grid(location = Locations, delta_climate = Delta_climates,
rho_k = Rho_ks, Term_time = School.Terms,
sigma_beta = Sigma_beta_vals)
model_labels <- c(CC_AC = "Autocorrelation",
CC_stand_smooth = "Time-smooth",
CC_stand_smooth_AC = "Time-smooth + Autocorrelation",
CC_true = "Control model")
cumulative_df = NULL
# For loop to generate big cumulative dataframe
for(i in 1:nrow(Conds)){
if(Conds[i,"delta_climate"] == -.02 & (Conds[i,"Term_time"] == T | Conds[i,"sigma_beta"] != 0)){next}
data_Path = paste("Results_Aggregated/R_sims_Te",
Conds[i, "delta_climate"], "_RH", Conds[i, "delta_climate"],
"_",
"sigma", Conds[i, "sigma_beta"], "_",
ifelse(rho_mean != 0, paste("rho_mean", rho_mean, "_", sep="") , ""),
Conds[i, "location"],"_",
sprintf("_rho_k%s", gsub("\\.", "-", as.character(Conds[i, "rho_k"]))),
sprintf("_import_%s", as.character(import)),
ifelse(Conds[i, "Term_time"] == T, "_School-Forced", ""),
sep = "")
final_pars = readRDS(paste(data_Path, "/final_pars.rds", sep = ""))
tmp_df = readRDS(paste(data_Path, "/sims_reg_tot.rds", sep = "")) |> filter(type == roi_model) |>
dplyr::ungroup() |>
dplyr::select(-type) |>
dplyr::mutate(.id = factor(.id),
rep = factor(rep)) |>
dplyr::rename(AB_Te = bias_Te, AB_RH = bias_RH) |>
dplyr::group_by(.id) |> dplyr::mutate(MAB_Te_id = mean(AB_Te),
MAB_RH_id = mean(AB_RH),
MSE_Te_id = mean(e_Te_se),
MSE_RH_id = mean(e_RH_se)
) |>
dplyr::ungroup() |>
dplyr::mutate(MAB_Te = mean(AB_Te),
MAB_RH = mean(AB_RH),
MSE_Te = mean(e_Te_se),
MSE_RH = mean(e_RH_se)) |>
dplyr::mutate(bias_v_id_Te = mean((MAB_Te_id - MAB_Te) ** 2),
bias_v_id_RH = mean((MAB_RH_id - MAB_RH) ** 2),
SE_v_id_Te = mean((MSE_Te_id - MSE_Te) ** 2),
SE_v_id_RH = mean((MSE_RH_id - MSE_RH) ** 2),
) |>
dplyr::mutate(bias_v_rep_Te = mean((AB_Te - MAB_Te_id) ** 2),
bias_v_rep_RH = mean((AB_RH - MAB_RH_id) ** 2),
SE_v_rep_Te = mean((e_Te_se - MSE_Te_id) ** 2),
SE_v_rep_RH = mean((e_RH_se - MSE_RH_id) ** 2)
) |>
dplyr::ungroup() |> dplyr::select(-c("rep",".id", "e_Te", "e_Te_se", "e_Te_low_CI",
"e_Te_high_CI", "e_RH", "e_RH_se", "e_RH_low_CI", "e_RH_high_CI",
"log_SI_lag_VAR", "mean_log_SI_lag", "R2", "AB_Te", "AB_RH",
"TP_Te", "e_Te_pow", "TP_RH", "e_RH_pow", "MAB_Te_id", "MAB_RH_id",
"MAB_Te_id", "MAB_RH_id", "MSE_Te_id", "MSE_RH_id")) |>
dplyr::distinct() |>
dplyr::mutate(
across(
.cols = c(bias_v_id_Te, bias_v_id_RH, bias_v_rep_Te, bias_v_rep_RH),
.fns = ~ sqrt(.) / abs(Conds[i,"delta_climate"]) * 100),
across(
.cols = c(MAB_Te, MAB_RH, MSE_Te, MSE_RH),
.fns = ~ . / abs(Conds[i,"delta_climate"]) * 100),
across(
.cols = c(SE_v_id_Te, SE_v_id_RH, SE_v_rep_Te, SE_v_rep_RH),
.fns = ~ sqrt(.) / abs(Conds[i,"delta_climate"]) * 100)
) |>
dplyr::rename(RMAB_Te = MAB_Te, RMAB_RH = MAB_RH, RMSE_Te = MSE_Te, RMSE_RH = MSE_RH) |>
dplyr::mutate(
across(matches("^(RMAB|bias)"), ~ round(.x, 0)),
across(matches("^(RMSE|SE)"), ~ round(.x, 3))
) |> cbind(Conds[i,])
cumulative_df <- rbind(cumulative_df, tmp_df)
}
clean_df <- cumulative_df |> dplyr::mutate(
location = recode(location,
Rome = "Rome, Italy",
Rio_de_Janeiro = "Rio de Janeiro, Brazil",
Toronto = "Toronto, Canada",
Dubai = "Dubai, UAE"
)
) |>
dplyr::mutate(
formatted_RMAB_Te = sprintf("%.f (%.f) [%.f]", RMAB_Te, bias_v_id_Te, bias_v_rep_Te),
formatted_RMAB_RH = sprintf("%.f (%.f) [%.f]", RMAB_RH, bias_v_id_RH, bias_v_rep_RH),
formatted_RMSE_Te = sprintf("%.f (%.f) [%.f]", RMSE_Te, SE_v_id_Te, SE_v_rep_Te),
formatted_RMSE_RH = sprintf("%.f (%.f) [%.f]", RMSE_RH, SE_v_id_RH, SE_v_rep_RH),
) |>
dplyr::select(-c("RMAB_Te", "RMAB_RH",
"RMSE_Te", "RMSE_RH",
"bias_v_id_Te", "bias_v_id_RH",
"SE_v_id_Te", "SE_v_id_RH",
"bias_v_rep_Te", "bias_v_rep_RH",
"SE_v_rep_Te", "SE_v_rep_RH")) |>
mutate(
Temperature = paste0(formatted_RMAB_Te, "<br/>", formatted_RMSE_Te),
`Relative Humidity` = paste0(formatted_RMAB_RH, "<br/>", formatted_RMSE_RH)
)
# Now create control table at baseline -----
baseline <- clean_df |> dplyr::rename(Location = location) |> dplyr::filter(rho_k == 0.1,
delta_climate == -0.04,
sigma_beta == 0,
Term_time == F) |> dplyr::select(-c("delta_climate", "rho_k", "Term_time", "sigma_beta",
"formatted_RMAB_Te", "formatted_RMAB_RH",
"formatted_RMSE_Te", "formatted_RMSE_RH")) |>
kable(format = "html", escape = FALSE, caption = roi_model) |>
# caption = "Baseline parameter values") |>
column_spec(1:3, latex_valign = "m") |>
kable_styling(bootstrap_options = c("bordered"), html_font = "Cambria") |>
row_spec(0:4, align = "c") |>
save_kable(file="Figures/D_table.png", density = 1800, zoom = 1.5)