Created
October 7, 2019 21:18
-
-
Save philipp-baumann/75c11fb8084e096b88f474ddb863b5a9 to your computer and use it in GitHub Desktop.
Example script related to https://github.com/rstudio/renv/issues/208
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ## Wrangle results of nested group k-fold cross-validation ===================== | |
| # Add additional grouping IDs | |
| cubist_nested_results_ids <- cubist_nested_add_ids( | |
| object = cubist_nested_results_repgroupcv, x_vars = "spc_pre", | |
| id_cols = c("site_id", "depth_cm"), | |
| y_vars = param_names, | |
| custom_spc_averaging = FALSE) | |
| # Get predicted and observed values for each outer resample assessment set | |
| cubist_nested_results_predobs_cm <- merge_predobs_y_vars( | |
| object = cubist_nested_results_ids, y_vars = param_names) | |
| # Join predicted and observed parameter names for each modeled hydraulic | |
| # parameter: | |
| cubist_outer_predobs_cm <- join_predobs_y_vars( | |
| object = cubist_nested_results_predobs_cm, | |
| vg_theta_s, vg_theta_r, vg_alpha, vg_n, | |
| kosugi_theta_s, kosugi_theta_r, kosugi_sigma, kosugi_h_mi) | |
| ## Evaluate internal model performance for all model hydraulic parameters | |
| ## from outer repeated group v-fold cross-validation =========================== | |
| # Convert the inflection point psi_0 from cm to kPa | |
| cubist_nested_results_predobs <- cubist_nested_results_predobs_cm %>% | |
| mutate( | |
| kosugi_h_mi = map(kosugi_h_mi, | |
| ~ mutate(.x, | |
| obs = 0.0980665 * obs, | |
| pred = 0.0980665 * pred)) | |
| ) | |
| # Create the model evaluation summary and include rules and majority consensus | |
| # for committees and neighbors | |
| # Create the model evaluation summary | |
| param_eval_nocomm <- metrics_cols(cubist_nested_results_predobs) | |
| param_eval <- cubist_nested_results_repgroupcv %>% | |
| select(-c(splits, id, inner_resamples, id_fold)) %>% | |
| summarize_outer_tuning_vars() %>% | |
| inner_join(x = ., y = param_eval_nocomm) %>% | |
| group_by(variable) %>% | |
| slice(1L) %>% | |
| ungroup() | |
| # Modify columns for annotation for model evaluation summary graph | |
| param_eval_graph <- mutate(param_eval, | |
| rmse = as.character(as.expression(paste0("RMSE == ", | |
| sprintf("%.2f", rmse_mean), | |
| "%+-%", | |
| round(rmse_sd, abs(round(log10(rmse_sd))))) | |
| )), | |
| r2 = as.character(as.expression(paste0("italic(R)^2 == ", | |
| sprintf("%.2f", r2_mean), | |
| "%+-%", | |
| sprintf("%.3f", r2_sd)) | |
| )), | |
| bias = as.character(as.expression(paste0("bias == ", | |
| round(bias_mean, abs(round(log10(bias_sd)))), | |
| "%+-%", | |
| round(bias_sd, abs(round(log10(bias_sd) + 1)))) | |
| )), | |
| n = as.character(as.expression(paste0("italic(n) == ", n_mean))), | |
| one_one = "1:1" | |
| ) %>% | |
| select(variable, rmse, r2, bias, n, one_one) | |
| ## Model assessment plot for hydraulic parameters ============================== | |
| # Extract predicted and observed values for all resampling results for the | |
| # fitted hydraulic parameters stored as list-columns | |
| params_predobs <- extract_predobs_cols(.data = cubist_nested_results_predobs) | |
| # Calculate mean and standard deviation of hold-out predictions for | |
| # all site-depth combinations across resampling repeats | |
| params_predobs_aggr_nometa <- params_predobs %>% | |
| group_by(variable, site_id, depth_cm) %>% | |
| mutate_at( | |
| .vars = vars("pred"), | |
| .funs = list(pred_mean = mean, pred_sd = sd)) %>% | |
| slice(1L) %>% | |
| select(-resample_id, -id_repeat) | |
| ## Append additional data for observed theta that can be used for | |
| ## grouping/annotating model evaluation plots; add another column filtering | |
| ## for theta higher or lower 0.6 | |
| df_category <- df %>% | |
| mutate( | |
| theta_category = dplyr::case_when( | |
| suction_cm == 10 & moisture_vol > 0.6 ~ "theta > 0.6", | |
| TRUE ~ "theta < 0.6" | |
| ) | |
| ) %>% | |
| select(site_id, depth_cm, soil_class, suction, suction_cm, moisture_vol, | |
| theta_category, clay) | |
| # Join new theta category data set to aggregated hydraulic model parameter | |
| # data set | |
| params_predobs_aggr <- params_predobs_aggr_nometa %>% | |
| inner_join(x = ., y = df_category) %>% | |
| group_by(variable, site_id, depth_cm) %>% | |
| slice(1L) | |
| # Create labeller for facets: name hydraulic parameters | |
| lbl_params <- as_labeller( | |
| x = c( | |
| `kosugi_h_mi` = "Kosugi~-italic(psi)['0']~'[kPa]'", | |
| `kosugi_sigma` = "Kosugi~sigma", | |
| `kosugi_theta_s` = "Kosugi~theta[s]~'['*cm^{-1}*cm^{-1}*']'", | |
| `kosugi_theta_r` = "Kosugi~theta[r]~'['*cm^{-1}*cm^{-1}*']'", | |
| `vg_alpha` = "Van~Genuchten~alpha~'['*cm^{-1}*']'", | |
| `vg_n` = "Van~Genuchten~italic(n)", | |
| `vg_theta_s` = "Van~Genuchten~theta[s]~'[%]'", | |
| `vg_theta_r` = "Van~Genuchten~theta[r]~'[%]'" | |
| ), | |
| default = label_parsed | |
| ) | |
| # Reorder factor levels to order panels | |
| param_levels <- params_predobs_aggr$variable %>% unique() | |
| params_predobs_aggr$variable <- factor(params_predobs_aggr$variable, | |
| levels = param_levels, ordered = TRUE) | |
| ## Create separate subfigures for the parameters to have identical x and y | |
| ## ranges which differ for all parameters | |
| ## Kosugi model | |
| param_eval_graph_kosugi <- mutate(param_eval, | |
| rmse = as.character(as.expression(paste0("RMSE == ", | |
| "'", | |
| ifelse(rmse_mean > 100, sprintf("%.0f", rmse_mean), | |
| ifelse(rmse_mean > 0.1, sprintf("%.2f", rmse_mean), | |
| sprintf("%.3f", rmse_mean))), | |
| "'") | |
| )), | |
| r2 = as.character(as.expression(paste0("italic(R)^2 == ", | |
| "'", | |
| sprintf("%.2f", r2_mean), | |
| "'") | |
| )), | |
| n = as.character(as.expression(paste0("italic(n) == ", n_mean))), | |
| one_one = "1:1" | |
| ) %>% | |
| select(variable, rmse, r2, n, one_one) %>% | |
| filter(variable %in% | |
| c("kosugi_theta_s", "kosugi_theta_r", "kosugi_sigma", "kosugi_h_mi")) | |
| # Model assessment of Kosugi theta_s ------------------------------------------- | |
| # xy range | |
| xyrange_kosugi_theta_s <- params_predobs_aggr %>% | |
| filter(variable == "kosugi_theta_s") %>% | |
| xy_range(x = obs, y = pred_mean) | |
| p_param_eval_kosugi_theta_s <- params_predobs_aggr %>% | |
| filter(variable == "kosugi_theta_s") %>% | |
| ggplot(data = ., | |
| aes(x = obs, y = pred_mean)) + | |
| geom_abline(slope = 1, intercept = 0, color = "black") + | |
| geom_point(aes(colour = soil_class, fill = soil_class), | |
| shape = 21, size = 2.5, alpha = 0.65) + | |
| geom_errorbar(aes(ymin = pred_mean - pred_sd, ymax = pred_mean + pred_sd, | |
| colour = soil_class)) + | |
| facet_wrap(~ variable, labeller = lbl_params, | |
| ncol = 2) + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_theta_s"), | |
| aes(x = Inf, y = Inf, label = one_one), size = 4, | |
| hjust = 2.5, vjust = 2, colour = "black") + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_theta_s"), | |
| aes(x = Inf, y = -Inf, label = r2), size = 4, | |
| hjust = 1.09, vjust = -2.5, parse = TRUE) + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_theta_s"), | |
| aes(x = Inf, y = -Inf, label = rmse), size = 4, | |
| hjust = 1.08, vjust = -1.25, parse = TRUE) + | |
| guides( | |
| colour = guide_legend(title = "ASC\norder"), | |
| fill = guide_legend(title = "ASC\norder") | |
| ) + | |
| scale_colour_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", | |
| "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a")) + | |
| scale_fill_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", | |
| "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a")) + | |
| xlim(xyrange_kosugi_theta_s[1] - 0.02 * diff(range(xyrange_kosugi_theta_s)), | |
| xyrange_kosugi_theta_s[2] + 0.02 * diff(range(xyrange_kosugi_theta_s))) + | |
| ylim(xyrange_kosugi_theta_s[1] - 0.02 * diff(range(xyrange_kosugi_theta_s)), | |
| xyrange_kosugi_theta_s[2] + 0.02 * diff(range(xyrange_kosugi_theta_s))) + | |
| coord_fixed() + | |
| xlab("") + # use cowplot to create axis annotation for later | |
| ylab("") + | |
| theme_bw() + | |
| theme( | |
| strip.background = element_rect(colour = "black", fill = NA), | |
| panel.grid.major = element_blank(), | |
| panel.grid.minor = element_blank(), | |
| axis.title.x = element_text(size = 12), | |
| axis.text.x = element_text(size = 12), | |
| axis.title.y = element_text(size = 12), | |
| axis.text.y = element_text(size = 12), | |
| legend.title = element_text(size = 12), | |
| legend.text = element_text(size = 12), | |
| strip.text.x = element_text(size = 12)) | |
| # Model assessment of Kosugi theta_r ------------------------------------------- | |
| # xy range | |
| xyrange_kosugi_theta_r <- params_predobs_aggr %>% | |
| filter(variable == "kosugi_theta_r") %>% | |
| xy_range(x = obs, y = pred_mean) | |
| p_param_eval_kosugi_theta_r <- params_predobs_aggr %>% | |
| filter(variable == "kosugi_theta_r") %>% | |
| ggplot(data = ., | |
| aes(x = obs, y = pred_mean)) + | |
| geom_abline(slope = 1, intercept = 0, color = "black") + | |
| geom_point(aes(colour = soil_class, fill = soil_class), | |
| shape = 21, size = 2.5, alpha = 0.65) + | |
| geom_errorbar(aes(ymin = pred_mean - pred_sd, ymax = pred_mean + pred_sd, | |
| colour = soil_class)) + | |
| facet_wrap(~ variable, labeller = lbl_params, | |
| ncol = 2) + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_theta_r"), | |
| aes(x = Inf, y = Inf, label = one_one), size = 4, | |
| hjust = 2.5, vjust = 2, colour = "black") + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_theta_r"), | |
| aes(x = Inf, y = -Inf, label = r2), size = 4, | |
| hjust = 1.09, vjust = -2.5, parse = TRUE) + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_theta_r"), | |
| aes(x = Inf, y = -Inf, label = rmse), size = 4, | |
| hjust = 1.08, vjust = -1.25, parse = TRUE) + | |
| guides( | |
| colour = guide_legend(title = "ASC\norder"), | |
| fill = guide_legend(title = "ASC\norder") | |
| ) + | |
| scale_colour_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", | |
| "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a")) + | |
| scale_fill_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", | |
| "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a")) + | |
| xlim(xyrange_kosugi_theta_r[1] - 0.02 * diff(range(xyrange_kosugi_theta_r)), | |
| xyrange_kosugi_theta_r[2] + 0.02 * diff(range(xyrange_kosugi_theta_r))) + | |
| ylim(xyrange_kosugi_theta_r[1] - 0.02 * diff(range(xyrange_kosugi_theta_r)), | |
| xyrange_kosugi_theta_r[2] + 0.02 * diff(range(xyrange_kosugi_theta_r))) + | |
| coord_fixed() + | |
| xlab("") + # use cowplot to create axis annotation for later | |
| ylab("") + | |
| theme_bw() + | |
| theme( | |
| strip.background = element_rect(colour = "black", fill = NA), | |
| panel.grid.major = element_blank(), | |
| panel.grid.minor = element_blank(), | |
| axis.title.x = element_text(size = 12), | |
| axis.text.x = element_text(size = 12), | |
| axis.title.y = element_text(size = 12), | |
| axis.text.y = element_text(size = 12), | |
| legend.title = element_text(size = 12), | |
| legend.text = element_text(size = 12), | |
| strip.text.x = element_text(size = 12)) | |
| # Model assessment of Kosugi sigma --------------------------------------------- | |
| # xy range | |
| xyrange_kosugi_sigma <- params_predobs_aggr %>% | |
| filter(variable == "kosugi_sigma") %>% | |
| xy_range(x = obs, y = pred_mean) | |
| p_param_eval_kosugi_sigma <- params_predobs_aggr %>% | |
| filter(variable == "kosugi_sigma") %>% | |
| ggplot(data = ., | |
| aes(x = obs, y = pred_mean)) + | |
| geom_abline(slope = 1, intercept = 0, color = "black") + | |
| geom_point(aes(colour = soil_class, fill = soil_class), | |
| shape = 21, size = 2.5, alpha = 0.65) + | |
| geom_errorbar(aes(ymin = pred_mean - pred_sd, ymax = pred_mean + pred_sd, | |
| colour = soil_class)) + | |
| facet_wrap(~ variable, labeller = lbl_params, | |
| ncol = 2) + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_sigma"), | |
| aes(x = Inf, y = Inf, label = one_one), size = 4, | |
| hjust = 2.5, vjust = 2, colour = "black") + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_sigma"), | |
| aes(x = Inf, y = -Inf, label = r2), size = 4, | |
| hjust = 1.09, vjust = -2.5, parse = TRUE) + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_sigma"), | |
| aes(x = Inf, y = -Inf, label = rmse), size = 4, | |
| hjust = 1.08, vjust = -1.25, parse = TRUE) + | |
| guides( | |
| colour = guide_legend(title = "ASC\norder"), | |
| fill = guide_legend(title = "ASC\norder") | |
| ) + | |
| scale_colour_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", | |
| "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a")) + | |
| scale_fill_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", | |
| "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a")) + | |
| xlim(xyrange_kosugi_sigma[1] - 0.02 * diff(range(xyrange_kosugi_sigma)), | |
| xyrange_kosugi_sigma[2] + 0.02 * diff(range(xyrange_kosugi_sigma))) + | |
| ylim(xyrange_kosugi_sigma[1] - 0.02 * diff(range(xyrange_kosugi_sigma)), | |
| xyrange_kosugi_sigma[2] + 0.02 * diff(range(xyrange_kosugi_sigma))) + | |
| coord_fixed() + | |
| xlab("") + # use cowplot to create axis annotation for later | |
| ylab("") + | |
| theme_bw() + | |
| theme( | |
| strip.background = element_rect(colour = "black", fill = NA), | |
| panel.grid.major = element_blank(), | |
| panel.grid.minor = element_blank(), | |
| axis.title.x = element_text(size = 12), | |
| axis.text.x = element_text(size = 12), | |
| axis.title.y = element_text(size = 12), | |
| axis.text.y = element_text(size = 12), | |
| legend.title = element_text(size = 12), | |
| legend.text = element_text(size = 12), | |
| strip.text.x = element_text(size = 12)) | |
| # Model assessment of Kosugi h_mi ---------------------------------------------- | |
| # xy range | |
| xyrange_kosugi_h_mi <- params_predobs_aggr %>% | |
| filter(variable == "kosugi_h_mi") %>% | |
| xy_range(x = obs, y = pred_mean) | |
| p_param_eval_kosugi_h_mi <- params_predobs_aggr %>% | |
| filter(variable == "kosugi_h_mi") %>% | |
| ggplot(data = ., | |
| aes(x = obs, y = pred_mean)) + | |
| geom_abline(slope = 1, intercept = 0, color = "black") + | |
| geom_point(aes(colour = soil_class, fill = soil_class), | |
| shape = 21, size = 2.5, alpha = 0.65) + | |
| geom_errorbar(aes(ymin = pred_mean - pred_sd, ymax = pred_mean + pred_sd, | |
| colour = soil_class)) + | |
| facet_wrap(~ variable, labeller = lbl_params, | |
| ncol = 2) + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_h_mi"), | |
| aes(x = Inf, y = Inf, label = one_one), size = 4, | |
| hjust = 2.5, vjust = 2, colour = "black") + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_h_mi"), | |
| aes(x = Inf, y = -Inf, label = r2), size = 4, | |
| hjust = 1.09, vjust = -2.5, parse = TRUE) + | |
| geom_text(data = param_eval_graph_kosugi %>% | |
| filter(variable == "kosugi_h_mi"), | |
| aes(x = Inf, y = -Inf, label = rmse), size = 4, | |
| hjust = 1.08, vjust = -1.25, parse = TRUE) + | |
| guides( | |
| colour = guide_legend(title = "ASC\norder"), | |
| fill = guide_legend(title = "ASC\norder") | |
| ) + | |
| scale_colour_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", | |
| "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a")) + | |
| scale_fill_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c", | |
| "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a")) + | |
| xlim(xyrange_kosugi_h_mi[1] - 0.02 * diff(range(xyrange_kosugi_h_mi)), | |
| xyrange_kosugi_h_mi[2] + 0.02 * diff(range(xyrange_kosugi_h_mi))) + | |
| ylim(xyrange_kosugi_h_mi[1] - 0.02 * diff(range(xyrange_kosugi_h_mi)), | |
| xyrange_kosugi_h_mi[2] + 0.02 * diff(range(xyrange_kosugi_h_mi))) + | |
| coord_fixed() + | |
| xlab("") + # use cowplot to create axis annotation for later | |
| ylab("") + | |
| theme_bw() + | |
| theme( | |
| strip.background = element_rect(colour = "black", fill = NA), | |
| panel.grid.major = element_blank(), | |
| panel.grid.minor = element_blank(), | |
| axis.title.x = element_text(size = 12), | |
| axis.text.x = element_text(size = 12), | |
| axis.title.y = element_text(size = 12), | |
| axis.text.y = element_text(size = 12), | |
| legend.title = element_text(size = 12), | |
| legend.text = element_text(size = 12), | |
| strip.text.x = element_text(size = 12)) | |
| # Rearrange Kosugi single panel graphs into a combined graph ------------------- | |
| kosugi_margin <- unit(c(t = 0.1, r = 0, b = 0.1, l = 0.02), "cm") | |
| p_kosugi_params <- cowplot::plot_grid( | |
| p_param_eval_kosugi_theta_s + | |
| theme(legend.position = "none", | |
| plot.margin = kosugi_margin), | |
| p_param_eval_kosugi_theta_r + | |
| theme(legend.position = "none", | |
| plot.margin = kosugi_margin), | |
| p_param_eval_kosugi_sigma + | |
| theme(legend.position = "none", | |
| plot.margin = kosugi_margin), | |
| p_param_eval_kosugi_h_mi + | |
| theme(legend.position = "none", | |
| plot.margin = kosugi_margin), | |
| ncol = 2, label_x = "a") | |
| p_kosugi_params_legend <- get_legend(p_param_eval_kosugi_theta_s) | |
| p_kosugi_comb <- cowplot::plot_grid(p_kosugi_params, p_kosugi_params_legend, | |
| rel_widths = c(3, 0.5)) + | |
| draw_label(expression(paste(italic(y))), x = 0.465, y = 0, vjust = -0.3, | |
| angle = 0) + | |
| draw_label(expression(paste(hat(italic(y)))), x = -0.0175, y = 0.5, vjust = 1.5, | |
| angle = 90) | |
| p_kosugi_comb_pdf <- ggsave(filename = "params-spc-pred-kosugi-soilclass.pdf", | |
| plot = p_kosugi_comb, | |
| path = here("out", "figs", "spc-models-params"), width = 6.5, height = 6) | |
| ## Tabular model evaluation ==================================================== | |
| # Replace and format paramter values for xtable and LaTeX | |
| pats <- c("theta_r", "theta_s", "alpha", "vg_n", "sigma", "h_mi") | |
| repl <- c( | |
| "$\\\\theta_\\\\textnormal{r}\\\\,[\\\\textnormal{cm}^{3}\\\\textnormal{cm}^{-3}]$", | |
| "$\\\\theta_\\\\textnormal{s}\\\\,[\\\\textnormal{cm}^{3}\\\\textnormal{cm}^{-3}]$", | |
| "$\\\\alpha$", "$n$", "$\\\\sigma$", "$-\\\\psi_0\\\\,[\\\\textnormal{kPa}]$") | |
| pats_col <- rep(list(pats), nrow(param_eval)) | |
| repl_col <- rep(list(repl), nrow(param_eval)) | |
| param_eval_table_repl <- param_eval %>% | |
| add_column(pats_col, repl_col) %>% | |
| mutate( | |
| variable_new = pmap(list(variable, pats_col, repl_col), | |
| ~ str_replace(string = ..1, pattern = ..2, replacement = ..3)) | |
| ) %>% | |
| mutate( | |
| variable_extr = map(.x = variable_new, | |
| ~ str_extract(string = .x, pattern = "\\${1}[:graph:]+\\${1}")) | |
| ) %>% | |
| mutate( | |
| variable = map_chr(variable_extr, | |
| ~ .x[!is.na(.x)]) | |
| ) | |
| param_eval_table_nonm <- param_eval_table_repl %>% | |
| mutate( | |
| model = c( | |
| "Van Genuchten", rep("", nrow(.) / 2 - 1) , | |
| "Kosugi", rep("", nrow(.) / 2 - 1)) # before: 3 | |
| ) %>% | |
| rename( | |
| parameter = variable, | |
| n = n_mean, | |
| min = min_mean, | |
| max = max_mean, | |
| mean = mean_mean, | |
| median = median_mean, | |
| sdev = sdev_mean, | |
| skewness = skewness_b1_mean | |
| ) %>% | |
| paste_plusminus(rmse_mean, rmse_sd, .digit_lag = 2) %>% | |
| paste_plusminus(me_mean, me_sd) %>% | |
| paste_plusminus(sde_mean, sde_sd) %>% | |
| paste_plusminus(r2_mean, r2_sd, digits = 2) | |
| vars_sel_params <- c("model", "parameter", "mean", | |
| "committees", "neighbors", "rules", "rmse", "me", "sde", "r2") | |
| param_eval_table <- param_eval_table_nonm %>% | |
| select( | |
| one_of(vars_sel_params)) | |
| colnames(param_eval_table) <- c("Model", "$\\mathcal P$", "Mean", | |
| "\\#{}\\textit{C}", "\\#{}\\textit{N}", "\\#{}Rules", | |
| "\\textsc{rmse}", "\\textsc{me}", "\\textsc{sde}", "$R^2$") | |
| # Only show results for Kosugi (drop Van Genuchten) ---------------------------- | |
| param_eval_table_kosugi_raw <- param_eval[5:8, ] %>% | |
| rename( | |
| parameter = variable, | |
| n = n_mean, | |
| min = min_mean, | |
| max = max_mean, | |
| mean = mean_mean, | |
| median = median_mean, | |
| sdev = sdev_mean, | |
| skewness = skewness_b1_mean | |
| ) | |
| param_eval_table_kosugi_rows_1_3 <- param_eval_table_kosugi_raw[1:3, ] %>% | |
| paste_plusminus(rmse_mean, rmse_sd, digits = 3) %>% | |
| paste_plusminus(me_mean, me_sd, digits = 3) %>% | |
| paste_plusminus(sde_mean, sde_sd, digits = 3) %>% | |
| paste_plusminus(r2_mean, r2_sd, digits = 2) | |
| param_eval_table_kosugi_row_4 <- param_eval_table_kosugi_raw[4, ] %>% | |
| paste_plusminus(rmse_mean, rmse_sd, digits = 1) %>% | |
| paste_plusminus(me_mean, me_sd, digits = 1) %>% | |
| paste_plusminus(sde_mean, sde_sd, digits = 1) %>% | |
| paste_plusminus(r2_mean, r2_sd, digits = 2) | |
| param_eval_table_kosugi_rows <- bind_rows( | |
| param_eval_table_kosugi_rows_1_3, | |
| param_eval_table_kosugi_row_4 | |
| ) | |
| param_eval_table_kosugi <- param_eval_table_kosugi_rows %>% | |
| select( | |
| one_of(vars_sel_params[!vars_sel_params %in% "model"])) | |
| colnames(param_eval_table_kosugi) <- c("$\\mathcal P$", "Mean", | |
| "\\#{}\\textit{C}", "\\#{}\\textit{N}", "\\#{}Rules", | |
| "RMSE", "ME", "SDE", "$R^2$") | |
| # Set number of digits per column and row | |
| # Specifiy digits for each row and column | |
| # https://stackoverflow.com/questions/14404241/set-digits-row-wise-with-print-xtable-in-r | |
| mdat_kosugi <- matrix( | |
| c(0, 0, 2, 0, 0, 0, 2, 2, 2, 2), # use recycling | |
| ncol = ncol(param_eval_table_kosugi) + 1, | |
| nrow = nrow(param_eval_table_kosugi), | |
| byrow = TRUE # change values | |
| ) | |
| # Change digits for min and max columns of psi_0 | |
| mdat_kosugi[4, c(7:9)] <- 1 | |
| mdat_kosugi[4, 3] <- 1 | |
| param_eval_xtable_kosugi <- xtable(param_eval_table_kosugi, | |
| caption = "Evaluation of the vis--NIR estimates of the Kosugi parameters. | |
| The spectroscopic modelling was based on the averaged spectra of the samples | |
| with water contents at the matric potentials between -1 and -1500~kPa, | |
| for each of the 162 soils from 54 sites and 3 depths. | |
| Mean values $\\pm$ standard deviations were calculated from 3-repeated | |
| 10-fold cross-validation. \\#\\textit{C} and \\#\\textit{N} denote the best | |
| number of {\\sc cubist} committees and neighbours, respectively, used in the | |
| models. \\#Rules are the number of rules. Values for \\#\\textit{C} and | |
| \\#\\textit{N} as well as \\#Rules report majority consensus values over all | |
| hold-out data sets of the nested cross-validation.", | |
| label = "tab:spc_params_kosugi", | |
| digits = mdat_kosugi) | |
| param_eval_xtable_kosugi_prt <- print(param_eval_xtable_kosugi, | |
| digits = mdat_kosugi, | |
| type = "latex", # this uses tabular environment | |
| # tabular.environment = "tabularx", | |
| # floating.environment = "sidewaystable", | |
| size = "\\footnotesize", | |
| hline.after = c(0, nrow(param_eval_xtable_kosugi)), | |
| include.rownames = FALSE, | |
| caption.placement = "top", | |
| sanitize.colnames.function = identity, | |
| sanitize.text.function = identity, | |
| # hline.after = NULL, # see link above | |
| booktabs = TRUE, | |
| # Add a table header | |
| add.to.row = list( | |
| pos = list(-1), | |
| command = c("\\toprule{}") | |
| ), | |
| file = here("out", "tables", "kosugi-eval.txt") | |
| ) | |
| ## Below is used for the appendix | |
| # Van Genuchten model only, drop Kosugi model results -------------------------- | |
| param_eval_table_vg_raw <- param_eval[1:4, ] %>% | |
| rename( | |
| parameter = variable, | |
| n = n_mean, | |
| min = min_mean, | |
| max = max_mean, | |
| mean = mean_mean, | |
| median = median_mean, | |
| sdev = sdev_mean, | |
| skewness = skewness_b1_mean | |
| ) | |
| param_eval_table_vg_rows_1_3 <- param_eval_table_vg_raw[1:3, ] %>% | |
| paste_plusminus(rmse_mean, rmse_sd, digits = 3) %>% | |
| paste_plusminus(me_mean, me_sd, digits = 3) %>% | |
| paste_plusminus(sde_mean, sde_sd, digits = 3) %>% | |
| paste_plusminus(r2_mean, r2_sd, digits = 2) | |
| param_eval_table_vg_row_4 <- param_eval_table_vg_raw[4, ] %>% | |
| paste_plusminus(rmse_mean, rmse_sd, digits = 2) %>% | |
| paste_plusminus(me_mean, me_sd, digits = 2) %>% | |
| paste_plusminus(sde_mean, sde_sd, digits = 2) %>% | |
| paste_plusminus(r2_mean, r2_sd, digits = 2) | |
| param_eval_table_vg_rows <- bind_rows( | |
| param_eval_table_vg_rows_1_3, | |
| param_eval_table_vg_row_4 | |
| ) | |
| param_eval_table_vg <- param_eval_table_vg_rows %>% | |
| select( | |
| one_of(vars_sel_params[!vars_sel_params %in% "model"])) | |
| colnames(param_eval_table_vg) <- c("$\\mathcal P$", "Mean", | |
| "\\#{}\\textit{C}", "\\#{}\\textit{N}", "\\#{}Rules", | |
| "RMSE", "ME", "SDE", "$R^2$") | |
| # Modify parameter symbols and units | |
| param_eval_table_vg$`$\\mathcal P$` <- c( | |
| "$\\theta_\\textrm{s}$ [$\\textrm{cm}^3 \\textrm{cm}^{-3}$]", | |
| "$\\theta_\\textrm{r}$ [$\\textrm{cm}^3 \\textrm{cm}^{-3}$]", | |
| "$\\alpha$", | |
| "$n$" | |
| ) | |
| mdat_vg <- matrix( | |
| c(0, 0, 2, 0, 0, 0, 2, 2, 2, 2), # use recycling | |
| ncol = ncol(param_eval_table_vg) + 1, | |
| nrow = nrow(param_eval_table_vg), | |
| byrow = TRUE # change values | |
| ) | |
| # Change some digits | |
| mdat_vg[4, c(7:9)] <- 1 | |
| mdat_vg[4, 3] <- 1 | |
| param_eval_xtable_vg <- xtable(param_eval_table_vg, | |
| caption = "Evaluation of the vis--NIR estimates of the van Genuchten parameters. | |
| The spectroscopic modelling was based on the averaged spectra of the samples | |
| with water contents at the matric potentials between -1 and -1500~kPa, | |
| for each of the 162 soils from 54 sites and 3 depths. | |
| Mean values $\\pm$ standard deviations were calculated from 3-repeated | |
| 10-fold cross-validation. \\#\\textit{C} and \\#\\textit{N} denote the best | |
| number of {\\sc cubist} committees and neighbours, respectively, used in the | |
| models. \\#Rules are the number of rules. Values for \\#\\textit{C} and | |
| \\#\\textit{N} as well as \\#Rules report majority consensus values over all | |
| hold-out data sets of the nested cross-validation.", | |
| label = "tab:spc_params_vg", | |
| digits = mdat_vg) | |
| param_eval_xtable_vg_prt <- print(param_eval_xtable_vg, | |
| digits = mdat_kosugi, | |
| type = "latex", # this uses tabular environment | |
| # tabular.environment = "tabularx", | |
| # floating.environment = "sidewaystable", | |
| size = "\\footnotesize", | |
| hline.after = c(0, nrow(param_eval_xtable_vg)), | |
| include.rownames = FALSE, | |
| caption.placement = "top", | |
| sanitize.colnames.function = identity, | |
| sanitize.text.function = identity, | |
| # hline.after = NULL, # see link above | |
| booktabs = TRUE, | |
| # Add a table header | |
| add.to.row = list( | |
| pos = list(-1), | |
| command = c("\\toprule{}") | |
| ), | |
| file = here("out", "tables", "vg-eval.txt") | |
| ) | |
| ## Predict new hydraulic curves based on Monte Carlo cross-validated | |
| ## outer assessment predictions of hydraulic model parameters ================== | |
| # `r` and `suction_cm_fit` already defined in file 40_model-swr-curve.R | |
| # r <- range(df$suction_cm, na.rm = TRUE) | |
| # suction_kpa_fit <- seq(from = r[1], to = r[2], length.out = 100) | |
| # suction_cm_fit <- seq(from = r[1], to = r[2], length.out = 100) | |
| # Combine all Kosugi parameters and pressure data to predict water contents at | |
| cubist_outer_predobs2 <- cubist_outer_predobs_cm %>% | |
| group_by(site_id, depth_cm, resample_id) %>% | |
| nest() %>% | |
| # Append pressure values to predict the curve | |
| mutate(suction_cm_fit = rep(list(suction_cm_fit))) %>% | |
| mutate(suction_cm_fit = map(suction_cm_fit, | |
| ~ tibble::tibble(suction_cm_fit = .x))) %>% | |
| mutate(data = map2(.x = data, .y = suction_cm_fit, | |
| ~ replicate_df(x = .x, y = .y))) %>% | |
| mutate(data = map2(.x = suction_cm_fit, .y = data, | |
| ~ dplyr::bind_cols(.x, .y))) | |
| # Compute theta predictions ---------------------------------------------------- | |
| ## Kosugi simplified form with `theta_c` = 0 | |
| # Predict outer assessment sets with Kosugi | |
| cubist_outer_kosugi_pred <- cubist_outer_predobs2 %>% | |
| mutate(kosugi_theta_pred = map(.x = data, | |
| ~ predict_kosugi(df = .x, | |
| x = suction_cm_fit, | |
| theta_r = kosugi_theta_r_pred, | |
| theta_s = kosugi_theta_s_pred, | |
| h_mi = kosugi_h_mi_pred, | |
| sigma = kosugi_sigma_pred))) | |
| # Unnest theta predicted by spectroscopy modeled Kosugi parameters | |
| cubist_outer_kosugi_unnested <- cubist_outer_kosugi_pred %>% | |
| select(-c(data, suction_cm_fit)) %>% | |
| unnest(kosugi_theta_pred) %>% | |
| rename(theta_pred = kosugi_theta_pred) %>% | |
| mutate(suction_fit_MPa = suction_cm_fit / 10000) | |
| # Extract and mutate fold id into a new column | |
| cubist_outer_kosugi_nocomb <- add_cv_fold(data = cubist_outer_kosugi_unnested, | |
| repeat_fold_col = resample_id) | |
| cubist_outer_kosugi_nocomb_rep <- add_cv_repeat( | |
| data = cubist_outer_kosugi_nocomb, | |
| repeat_fold_col = resample_id) | |
| # Add combined `repeat` and `depth` id | |
| cubist_outer_kosugi <- cubist_outer_kosugi_nocomb_rep %>% | |
| mutate(depth_repeat = paste(depth_cm, id_repeat, sep = "_")) | |
| ## Van Genuchten | |
| # Predict outer assessment sets with Van Genuchten | |
| cubist_outer_vg_pred <- cubist_outer_predobs2 %>% | |
| mutate(vg_theta_pred = map(.x = data, | |
| ~ predict_van_genuchten(df = .x, | |
| x = suction_cm_fit, | |
| theta_s = vg_theta_s_pred, | |
| theta_r = vg_theta_r_pred, | |
| alpha = vg_alpha_pred, | |
| n = vg_n_pred))) | |
| # Unnest theta predicted by spectroscopy modeled Van Genuchten parameters | |
| cubist_outer_vg_unnested <- cubist_outer_vg_pred %>% | |
| select(-c(data, suction_cm_fit)) %>% | |
| unnest(vg_theta_pred) %>% | |
| rename(theta_pred = vg_theta_pred) %>% | |
| mutate(suction_fit_MPa = suction_cm_fit / 10000) | |
| # Extract and mutate fold id into a new column | |
| cubist_outer_vg_nocomb <- add_cv_fold(data = cubist_outer_vg_unnested, | |
| repeat_fold_col = resample_id) | |
| cubist_outer_vg_nocomb_rep <- add_cv_repeat(data = cubist_outer_vg_nocomb, | |
| repeat_fold_col = resample_id) | |
| # Add combined `repeat` and `depth` id | |
| cubist_outer_vg <- cubist_outer_vg_nocomb_rep %>% | |
| mutate(depth_repeat = paste(depth_cm, id_repeat, sep = "_")) | |
| ## Plot reconstructed water rentention curves using spectroscopically predicted | |
| ## hydraulic model parameters ================================================== | |
| # Overlap spectroscopy-based Kosugi predicted theta lines using outer assessment | |
| # samples with NLS fitted values based on observed data ------------------------ | |
| p_site_kosugi_spc <- df %>% | |
| filter(suction != "air_dry") %>% | |
| ggplot(aes(x = suction_kpa, y = moisture_vol)) + | |
| geom_point(aes(colour = depth_cm, shape = depth_cm)) + | |
| # Plot prediction curve based on spectroscopy predicted Kosugi paramters | |
| # for outer assessment sets | |
| geom_line(aes(x = suction_fit_kPa, y = theta_pred, | |
| colour = depth_cm, linetype = depth_repeat), | |
| data = cubist_outer_kosugi %>% | |
| mutate(suction_fit_kPa = 1000 * suction_fit_MPa), | |
| inherit.aes = FALSE) + | |
| facet_wrap(~ site_id) + | |
| xlab(expression(paste(- italic(psi), " [kPa]"))) + | |
| ylab(expression(paste(italic(theta), " [", cm^{-3}, cm^{-3}, "]"))) + | |
| scale_colour_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) + | |
| guides(colour = guide_legend(title = "Depth [cm]"), | |
| linetype = FALSE, | |
| shape = guide_legend(title = "Depth [cm]")) + | |
| scale_x_log10( | |
| breaks = scales::trans_breaks("log10", function(x) 10^x), | |
| labels = scales::trans_format("log10", scales::math_format(10^.x)) | |
| ) + | |
| theme_bw() + | |
| theme( | |
| strip.background = element_rect(colour = "black", fill = NA), | |
| panel.grid.major = element_blank(), | |
| panel.grid.minor = element_blank() | |
| ) | |
| p_site_kosugi_spc_pdf <- ggsave( | |
| filename = "h2o-ret-curve-per-site-kosugi-spc-pred-repgroupcv.pdf", | |
| plot = p_site_kosugi_spc, | |
| path = here("out", "figs", "swr-curves"), | |
| width = 10, height = 10) | |
| ## Show Kosugi re-prediction for exemplay sites | |
| p_site_kosugi_spc_example <- df_kosugi_example %>% | |
| ggplot(aes(x = suction_kpa, y = moisture_vol)) + | |
| geom_point(aes(colour = depth_cm, shape = depth_cm)) + | |
| # Plot prediction curve based on spectroscopy predicted Kosugi paramters | |
| # for outer assessment sets | |
| geom_line(aes(x = suction_fit_kPa, y = theta_pred, group = depth_repeat, | |
| colour = depth_cm, linetype = depth_cm), | |
| data = cubist_outer_kosugi %>% | |
| ungroup() %>% | |
| filter(site_id %in% c("apsru2", "apsru28", "apsru17")) %>% # 5 | |
| mutate(site_id = replace(site_id, site_id == "apsru2", "Dermosol")) %>% | |
| mutate(site_id = replace(site_id, site_id == "apsru28", "Vertosol")) %>% | |
| mutate(site_id = replace(site_id, site_id == "apsru17", "Chromosol")) %>% | |
| mutate(suction_fit_kPa = 1000 * suction_fit_MPa), | |
| inherit.aes = FALSE) + | |
| # geom_text_repel( | |
| # data = df_kosugi_example %>% filter(suction_cm == 50), | |
| # aes(label = df_kosugi_example %>% | |
| # filter(suction_cm == 50) %>% | |
| # pull(clay) %>% paste0("%")), | |
| # size = 2.5, | |
| # # nudge_x = 40, | |
| # # segment.color = NA | |
| # ) + | |
| facet_wrap(~ site_id) + | |
| xlab(expression(paste(- italic(psi), " [kPa]"))) + | |
| ylab(expression(paste(italic(theta), " [", cm^{-3}, cm^{-3}, "]"))) + | |
| scale_colour_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) + | |
| guides(colour = guide_legend(title = "Depth [cm]"), | |
| linetype = guide_legend(title = "Depth [cm]"), | |
| shape = guide_legend(title = "Depth [cm]")) + | |
| scale_x_log10( | |
| breaks = scales::trans_breaks("log10", function(x) 10^x), | |
| labels = scales::trans_format("log10", scales::math_format(10^.x)) | |
| ) + | |
| theme_bw() + | |
| theme( | |
| strip.background = element_rect(colour = "black", fill = NA), | |
| panel.grid.major = element_blank(), | |
| panel.grid.minor = element_blank(), | |
| axis.title.x = element_text(size = 9), | |
| axis.text.x = element_text(size = 9), | |
| axis.title.y = element_text(size = 9), | |
| axis.text.y = element_text(size = 9), | |
| legend.title = element_text(size = 9), | |
| legend.text = element_text(size = 9), | |
| strip.text.x = element_text(size = 9), | |
| aspect.ratio = 1) | |
| p_site_kosugi_spc_example_pdf <- ggsave( | |
| filename = "h2o-ret-curve-per-site-kosugi-spc-pred-repgroupcv-example.pdf", | |
| plot = p_site_kosugi_spc_example, | |
| path = here("out", "figs", "swr-curves"), | |
| width = 5, height = 2) | |
| # Recombine hydraulic model fits and spectroscopy-based curve reconstruction | |
| # This is the Example fit | |
| p_kosugi_nls_spc <- plot_grid( | |
| p_curve_site_kosugi_example + ylim(c(0, 0.9)) + | |
| theme(legend.position = "none", plot.margin = unit(c(0, 0, 0, 0.75), "cm")), | |
| get_legend(p_curve_site_kosugi_example + | |
| theme(legend.position = "right")), | |
| p_site_kosugi_spc_example + ylim(c(0, 0.9)) + | |
| theme(legend.position = "none", plot.margin = unit(c(0, 0, 0, 0.75), "cm")), | |
| labels = c("a", "", "b", ""), | |
| ncol = 2, nrow = 2, rel_widths = c(4, 0.8)) | |
| p_kosugi_nls_spc_comb <- plot_grid( | |
| p_curve_site_kosugi_example + ylim(c(0, 0.9)) + | |
| theme(legend.position = "none", plot.margin = unit(c(0, 0, 0, 0.75), "cm")), | |
| p_site_kosugi_spc_example + ylim(c(0, 0.9)) + | |
| theme(legend.position = "none", plot.margin = unit(c(0, 0, 0, 0.75), "cm")), | |
| labels = c("a", "b"), label_size = 12, | |
| nrow = 2) | |
| p_kosugi_nls_spc_grid <- plot_grid( | |
| p_kosugi_nls_spc_comb, | |
| get_legend(p_curve_site_kosugi_example + | |
| theme(legend.position = "right")), | |
| ncol = 2, rel_widths = c(4, 0.8)) | |
| p_kosugi_nls_spc_grid_pdf <- ggsave( | |
| filename = "swrc-kosugi-nls-spc-example.pdf", | |
| path = here("out", "figs", "swr-curves"), | |
| plot = p_kosugi_nls_spc_grid, | |
| width = 5, height = 3.75) | |
| # Overlap spectroscopy-based Van Genuchten predicted theta lines using | |
| # outer assessment samples with NLS fitted values based on observed data ------- | |
| p_site_vg_spc <- df %>% | |
| filter(suction != "air_dry") %>% | |
| ggplot(aes(x = suction_kpa, y = moisture_vol)) + | |
| geom_point(aes(colour = depth_cm, shape = depth_cm)) + | |
| # Plot prediction curve based on spectroscopy predicted Van Genuchten | |
| # paramters for outer assessment sets | |
| geom_line(aes(x = suction_fit_kPa, y = theta_pred, | |
| colour = depth_cm, linetype = depth_repeat), | |
| data = cubist_outer_vg %>% | |
| mutate(suction_fit_kPa = 1000 * suction_fit_MPa), | |
| inherit.aes = FALSE) + | |
| facet_wrap(~ site_id) + | |
| xlab(expression(paste(- italic(psi), " [kPa]"))) + | |
| ylab(expression(paste(italic(theta), " [", cm^{-3}, cm^{-3}, "]"))) + | |
| scale_colour_manual(values = c("#1b9e77", "#d95f02", "#7570b3")) + | |
| guides( | |
| colour = guide_legend(title = "Depth [cm]"), | |
| linetype = FALSE, | |
| shape = guide_legend(title = "Depth [cm]")) + | |
| scale_x_log10( | |
| breaks = scales::trans_breaks("log10", function(x) 10^x), | |
| labels = scales::trans_format("log10", scales::math_format(10^.x)) | |
| ) + | |
| theme_bw() + | |
| theme( | |
| strip.background = element_rect(colour = "black", fill = NA), | |
| panel.grid.major = element_blank(), | |
| panel.grid.minor = element_blank() | |
| ) | |
| p_site_vg_spc_pdf <- ggsave( | |
| filename = "h2o-ret-curve-per-site-vg-spc-pred-repgroupcv.pdf", | |
| plot = p_site_vg_spc, | |
| path = here("out", "figs", "swr-curves"), | |
| width = 10, height = 10) | |
| ## Assess available water (theta_a) from reconstructed water curve predictions | |
| ## between field capacity (theta_fc) and | |
| ## the permanent wilting point (theta_pwp) ===================================== | |
| # Prepare new `data` list-column that contains suction pressure values at | |
| # field capacity (- 0.33bar) and the permanent wilting point (-15bar) | |
| # predict based on hydraulic model equations and spectroscopically predicted | |
| # hydraulic model parameters --------------------------------------------------- | |
| suction_cm_awc <- c(300, 15295) | |
| # Combine all Kosugi parameters and pressure data to predict field capacity | |
| # and the permanent wilting point ---------------------------------------------- | |
| cubist_outer_predobs_awc <- cubist_outer_predobs_cm %>% | |
| group_by(site_id, depth_cm, resample_id) %>% | |
| nest() %>% | |
| # Append pressure values to predict the curve | |
| mutate(suction_cm_awc = rep(list(suction_cm_awc))) %>% | |
| mutate(suction_cm_awc = map(suction_cm_awc, | |
| ~ tibble::tibble(suction_cm_fit = .x))) %>% | |
| mutate(data = map2(.x = data, .y = suction_cm_awc, | |
| ~ replicate_df(x = .x, y = .y))) %>% | |
| mutate(data = map2(.x = suction_cm_awc, .y = data, | |
| ~ dplyr::bind_cols(.x, .y))) | |
| # Compute theta predictions at field capacity and | |
| # the permanent wilting point and ---------------------------------------------- | |
| ## Kosugi | |
| # Predict outer assessment sets with Kosugi | |
| cubist_outer_kosugi_awc_nested <- cubist_outer_predobs_awc %>% | |
| mutate(kosugi_theta_pred = map(.x = data, | |
| ~ predict_kosugi(df = .x, | |
| x = suction_cm_fit, | |
| theta_r = kosugi_theta_r_pred, | |
| theta_s = kosugi_theta_s_pred, | |
| h_mi = kosugi_h_mi_pred, | |
| sigma = kosugi_sigma_pred))) | |
| # Unnest theta predicted by spectroscopy modeled Kosugi parameters | |
| cubist_outer_kosugi_awc_unnested <- | |
| cubist_outer_kosugi_awc_nested %>% | |
| select(-c(data, suction_cm_awc)) %>% | |
| unnest(kosugi_theta_pred) %>% | |
| mutate(suction_fit_MPa = suction_cm_fit / 10000) | |
| ## Van Genuchten | |
| # Predict outer assessment sets with Van Genuchten | |
| cubist_outer_vg_awc_nested <- | |
| cubist_outer_predobs_awc %>% | |
| mutate(vg_theta_pred = map(.x = data, | |
| ~ predict_van_genuchten(df = .x, | |
| x = suction_cm_fit, | |
| theta_s = vg_theta_s_pred, | |
| theta_r = vg_theta_r_pred, | |
| alpha = vg_alpha_pred, | |
| n = vg_n_pred))) | |
| # Unnest theta predicted by spectroscopy modeled Van Genuchten parameters | |
| cubist_outer_vg_awc_unnested <- | |
| cubist_outer_vg_awc_nested %>% | |
| select(-c(data, suction_cm_awc)) %>% | |
| unnest(vg_theta_pred) %>% | |
| mutate(suction_fit_MPa = suction_cm_fit / 10000) | |
| # Rearrange predicted theta_fc and theta_pwc and calculate theta_a ------------- | |
| ## Van Genuchten | |
| cubist_outer_vg_awc <- cubist_outer_vg_awc_unnested %>% | |
| group_by(site_id, depth_cm, resample_id) %>% | |
| # creates list-column <data> nested by <site_depth>, | |
| # which contains every other column | |
| nest() %>% | |
| mutate( | |
| awc = map(data, | |
| ~ calculate_awc_from_cm(., suction = suction_cm_fit, | |
| theta = vg_theta_pred))) %>% | |
| unnest(awc) %>% | |
| # Rename `theta_pwp`, `theta_fc` and `theta_a` | |
| rename( | |
| vg_theta_fc = theta_fc, | |
| vg_theta_pwp = theta_pwp, | |
| vg_theta_a = theta_a | |
| ) %>% | |
| # remove data list-column | |
| select(-data) | |
| cubist_outer_kosugi_awc <- cubist_outer_kosugi_awc_unnested %>% | |
| group_by(site_id, depth_cm, resample_id) %>% | |
| # creates list-column <data> nested by <site_depth>, | |
| # which contains every other column | |
| nest() %>% | |
| mutate( | |
| awc = map(data, | |
| ~ calculate_awc_from_cm(., suction = suction_cm_fit, | |
| theta = kosugi_theta_pred))) %>% | |
| unnest(awc) %>% | |
| # Rename `theta_pwp`, `theta_fc` and `theta_a` | |
| rename( | |
| kosugi_theta_fc = theta_fc, | |
| kosugi_theta_pwp = theta_pwp, | |
| kosugi_theta_a = theta_a | |
| ) %>% | |
| # remove data list-column | |
| select(-data) | |
| # Merge theta_fc, theta_pwc, and theta_awc for Kosugi and Van Genuchten, and | |
| # the measured (observed) values ----------------------------------------------- | |
| awc_obs <- spc_awc %>% select(site_id, depth_cm, soil_class, | |
| theta_fc, theta_pwp, theta_a) | |
| awc_cubist_hydro_params <- Reduce(inner_join, | |
| list(awc_obs, cubist_outer_kosugi_awc, cubist_outer_vg_awc)) | |
| # Perform model evaluation for Cubist resampling-based theta_a spectroscopic | |
| # prediction ------------------------------------------------------------------- | |
| ## Plot average (by `site_id` and `depth_cm` outer resample predictions derived | |
| ## with 10 times repeated | |
| ## (0.67/0.33) Monte Carlo cross-validation | |
| awc_cubist_hydro_params_aggr <- awc_cubist_hydro_params %>% | |
| group_by(site_id, depth_cm) %>% | |
| # Average `theta_a` and predicted `vg_theta_a` and `kosugi_theta_a` | |
| # by `site_id` and `depth_cm` across resampling sets | |
| mutate( | |
| kosugi_theta_a_mean = mean(kosugi_theta_a), | |
| kosugi_theta_a_sd = sd(kosugi_theta_a), | |
| vg_theta_a_mean = mean(vg_theta_a), | |
| vg_theta_a_sd = sd(vg_theta_a) | |
| ) %>% | |
| slice(1L) %>% | |
| mutate( | |
| obs = theta_a, | |
| pred = kosugi_theta_a_mean, | |
| pred_sd = kosugi_theta_a_sd | |
| ) %>% | |
| mutate(method = c("awc_kosugi")) | |
| # Create model evaluation summaries on sample average hold-out predictions | |
| # across resamples ------------------------------------------------------------- | |
| awc_kosugi_eval <- awc_cubist_hydro_params %>% | |
| # Add the repeat id called `id_repeat` | |
| add_cv_repeat(repeat_fold_col = resample_id, repeat_col = id_repeat) %>% | |
| group_by(id_repeat) %>% | |
| nest() %>% | |
| mutate(metrics = map(data, | |
| ~ simplerspec::summary_df( | |
| df = as.data.frame(.x), x = "theta_a", y = "kosugi_theta_a")) | |
| ) %>% | |
| unnest(metrics) %>% | |
| ungroup() %>% | |
| mutate( | |
| n = mean(n), | |
| r2_sd = sd(r2), | |
| rmse_sd = sd(rmse), | |
| me_sd = sd(me), | |
| sde_sd = sd(sde), | |
| rmse = mean(rmse), | |
| r2 = mean(r2), | |
| me = mean(me), | |
| sde = mean(sde) | |
| ) %>% | |
| add_column(variable = "theta_a", .before = 1) %>% | |
| mutate(method = "awc_kosugi") %>% | |
| slice(1L) | |
| ## Kosugi model summary for available water prediction comparision across | |
| ## methods (concept of the paper) | |
| # Calculate model performance measures by Resample | |
| model_summary_awc_kosugi <- awc_cubist_hydro_params %>% | |
| # Add the repeat id called `id_repeat` | |
| add_cv_repeat(repeat_fold_col = resample_id, repeat_col = id_repeat) %>% | |
| group_by(id_repeat) %>% | |
| nest() %>% | |
| mutate(metrics = map(data, | |
| ~ simplerspec::summary_df( | |
| df = as.data.frame(.x), x = "theta_a", y = "kosugi_theta_a")) | |
| ) %>% | |
| unnest(metrics) %>% | |
| ungroup() %>% | |
| summarize( | |
| n = mean(n), | |
| r2_sd = sd(r2), | |
| rmse_sd = sd(rmse), | |
| rmse = mean(rmse), | |
| r2 = mean(r2), | |
| bias = mean(bias) | |
| ) %>% | |
| # Modify columns for plot annotation | |
| mutate( | |
| rmse = as.character(as.expression(paste0("RMSE == ", round(rmse, 3), | |
| "%+-%", round(rmse_sd, 3)))), | |
| r2 = as.character(as.expression(paste0("italic(R)^2 == ", round(r2, 2), | |
| "%+-%", round(r2_sd, 2)))), | |
| n = as.character(as.expression(paste0("italic(n) == ", n))), | |
| one_one = "1:1" | |
| ) %>% | |
| mutate(method = "awc_kosugi") | |
| ## Van Genuchten | |
| vg_hydro_spc_summary <- simplerspec::summary_df( | |
| df = as.data.frame(awc_cubist_hydro_params_aggr), | |
| x = "vg_theta_a", y = "vg_theta_a_mean") %>% | |
| mutate( | |
| rmse = as.character(as.expression(paste0("RMSE == ", | |
| round(rmse, 2)))), | |
| r2 = as.character(as.expression(paste0("italic(R)^2 == ", | |
| round(r2, 2)))), | |
| rpd = as.character(as.expression(paste("RPD == ", | |
| round(rpd, 2)))), | |
| n = as.character(as.expression(paste0("italic(n) == ", n))) | |
| ) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment