## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.align = "center"
)
library(rsofun)
library(dplyr)
library(ggplot2)
library(tidyr)
library(sensitivity)
library(BayesianTools)
## -----------------------------------------------------------------------------
# Define log-likelihood function
ll_pmodel <- function(
par_v # a vector of all calibratable parameters including errors
){
rsofun::cost_likelihood_pmodel( # reuse likelihood cost function
par_v,
obs = rsofun::p_model_validation,
drivers = rsofun::p_model_drivers,
targets = "gpp"
)
}
# Compute log-likelihood for a given set of parameters
ll_pmodel( par_v = c(
kphio = 0.09423773, # setup ORG in Stocker et al. 2020 GMD
kphio_par_a = 0.0, # set to zero to disable temperature-dependence of kphio
kphio_par_b = 1.0,
soilm_thetastar = 0.6 * 240, # to recover old setup with soil moisture stress
soilm_betao = 0.0,
beta_unitcostratio = 146.0,
rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous
tau_acclim = 30.0,
kc_jmax = 0.41,
error_gpp = 0.9 # value from previous simulations
))
## -----------------------------------------------------------------------------
# best parameter values (from previous literature)
par_cal_best <- c(
kphio = 0.09423773,
kphio_par_a = -0.0025,
kphio_par_b = 20,
soilm_thetastar = 0.6*240,
soilm_betao = 0.2,
beta_unitcostratio = 146.0,
rd_to_vcmax = 0.014,
tau_acclim = 30.0,
kc_jmax = 0.41,
error_gpp = 1
)
# lower bound
par_cal_min <- c(
kphio = 0.03,
kphio_par_a = -0.004,
kphio_par_b = 10,
soilm_thetastar = 0,
soilm_betao = 0,
beta_unitcostratio = 50.0,
rd_to_vcmax = 0.01,
tau_acclim = 7.0,
kc_jmax = 0.2,
error_gpp = 0.01
)
# upper bound
par_cal_max <- c(
kphio = 0.15,
kphio_par_a = -0.001,
kphio_par_b = 30,
soilm_thetastar = 240,
soilm_betao = 1,
beta_unitcostratio = 200.0,
rd_to_vcmax = 0.1,
tau_acclim = 60.0,
kc_jmax = 0.8,
error_gpp = 4
)
## ----eval = FALSE-------------------------------------------------------------
# morris_setup <- BayesianTools::createBayesianSetup(
# likelihood = ll_pmodel,
# prior = BayesianTools::createUniformPrior(par_cal_min, par_cal_max, par_cal_best),
# names = names(par_cal_best)
# )
## ----eval = FALSE, echo = TRUE------------------------------------------------
# set.seed(432)
# morrisOut <- sensitivity::morris(
# model = morris_setup$posterior$density,
# factors = names(par_cal_best),
# r = 1000,
# design = list(type = "oat", levels = 20, grid.jump = 3),
# binf = par_cal_min,
# bsup = par_cal_max,
# scale = TRUE)
## ----simulate_morrisOut, include = FALSE--------------------------------------
# TODO: get rid of this and always fully run the vignettes
# fake output since Morris sensitivity isn't run
# saveRDS(morrisOut, "files/sensitivity_analysis.Rmd__morrisOut.RDS")
morrisOut <- readRDS("files/sensitivity_analysis.Rmd__morrisOut.RDS")
## ----eval = TRUE, fig.width=7, fig.height=4-----------------------------------
# Summarise the morris output
morrisOut.df <- data.frame(
parameter = names(par_cal_best),
mu.star = apply(abs(morrisOut$ee), 2, mean, na.rm = T),
sigma = apply(morrisOut$ee, 2, sd, na.rm = T)
) %>%
arrange( mu.star )
morrisOut.df |>
tidyr::pivot_longer( -parameter, names_to = "variable", values_to = "value") |>
ggplot(aes(
reorder(parameter, value),
value,
fill = variable),
color = NA) +
geom_bar(position = position_dodge(), stat = 'identity') +
scale_fill_brewer("", labels = c('mu.star' = expression(mu * "*"),
'sigma' = expression(sigma)),
palette = "Greys") +
theme_classic() +
theme(
axis.text = element_text(size = 6),
axis.title = element_blank(),
legend.position = c(0.05, 0.95), legend.justification = c(0.05, 0.95)
)
## ----eval = FALSE, echo = TRUE------------------------------------------------
# # Calibrates kphio, kphio_par_b, kc_jmax - top 3 model params
# set.seed(2023)
#
# # Define calibration settings
# settings_calib <- list(
# method = "BayesianTools",
# metric = rsofun::cost_likelihood_pmodel,
# control = list(
# sampler = "DEzs",
# settings = list(
# burnin = 3000,
# iterations = 9000,
# nrChains = 3, # number of independent chains
# startValue = 3 # number of internal chains to be sampled
# )),
# par = list(
# kphio = list(lower = 0.03, upper = 0.15, init = 0.05),
# #kphio_par_a = list(lower = -0.004, upper = -0.001, init = -0.0025),
# kphio_par_b = list(lower = 10, upper = 30, init =25),
# kc_jmax = list(lower = 0.2, upper = 0.8, init = 0.41),
# err_gpp = list(lower = 0.1, upper = 3, init = 0.8)
# )
# )
# par_fixed <- list(
# soilm_thetastar = 0.6*240,
# soilm_betao = 0.2,
# beta_unitcostratio = 146.0,
# rd_to_vcmax = 0.014,
# kphio_par_a = -0.0025,
# tau_acclim = 30.0)
#
# # Calibrate kphio-related parameters and err_gpp
# par_calib <- calib_sofun(
# drivers = p_model_drivers,
# obs = p_model_validation,
# settings = settings_calib,
# par_fixed = par_fixed,
# targets = "gpp"
# )
#
# # This code takes 15 minutes to run
## ----simulate_par_calib, include = FALSE--------------------------------------
# TODO: get rid of this and always fully run the vignettes
# fake output since calibration isn't run
# saveRDS(par_calib, "files/sensitivity_analysis.Rmd__par_calib.RDS")
par_calib <- readRDS("files/sensitivity_analysis.Rmd__par_calib.RDS")
## ----fig.height = 10, fig.width = 7-------------------------------------------
plot(par_calib$mod)
## ----fig.height = 10, fig.width = 7, eval = FALSE, echo = FALSE---------------
# # Define function for plotting chains separately
# plot_acf_mcmc <- function(chains, par_names){
# # chains: from the BayesianTools output
# n_chains <- length(chains)
# n_internal_chains <- length(chains[[1]]$chain)
# par(mfrow = c(length(par_names), n_chains))
# for(par_name in par_names){
# for(i in 1:n_chains){
# stopifnot(n_internal_chains<=3); color = c("blue", "red", "darkgreen")
# spacing = 0.5/n_internal_chains
# for(j in 1:n_internal_chains){
# autocorr_internal_chain <- pacf(getSample(chains[[i]]$chain[[j]])[, par_name], plot = FALSE)
# if(j==1){
# plot(autocorr_internal_chain, col = color[j],
# main = sprintf("Series of %s , chain (%i)", par_name, i))
# } else {
# lines(autocorr_internal_chain$lag + spacing*(j-1),
# autocorr_internal_chain$acf,
# col = color[j], type = "h")
# }
# }
# }
# }
# }
# plot_acf_mcmc(par_calib$mod, c("kphio", "kc_jmax", "kphio_par_b", "err_gpp"))
## ----fig.width=5, fig.height=5------------------------------------------------
correlationPlot(par_calib$mod, thin = 1) # use all samples, no thinning
## -----------------------------------------------------------------------------
gelmanDiagnostics(par_calib$mod)
## -----------------------------------------------------------------------------
summary(par_calib$mod)
## ----echo = TRUE, eval = FALSE------------------------------------------------
# # Evaluation of the uncertainty coming from the model parameters' uncertainty
#
# # Sample parameter values from the posterior distribution
# samples_par <- getSample(par_calib$mod, numSamples = 600) |>
# as.data.frame() |>
# dplyr::mutate(mcmc_id = 1:n()) |>
# tidyr::nest(.by = mcmc_id, .key = "pars")
#
# run_pmodel <- function(sample_par){
# # Function that runs the P-model for a sample of parameters
# # and also adds the new observation error
# out <- runread_pmodel_f(
# drivers = p_model_drivers,
# par = c(par_fixed, sample_par) # err_gpp is also here, but simply ignored
# )
#
# # return modelled GPP and prediction for a new GPP observation
# gpp <- out$data[[1]][, "gpp"]
# data.frame(gpp = gpp,
# gpp_pred = gpp + rnorm(n = length(gpp), mean = 0,
# sd = sample_par$err_gpp),
# date = out$data[[1]][, "date"])
# }
#
# set.seed(2023)
# # Run the P-model for each set of parameters
# pmodel_runs <- samples_par |>
# dplyr::mutate(sim = purrr::map(pars, ~run_pmodel(.x))) |>
# # format to obtain 90% credible intervals
# dplyr::select(mcmc_id, sim) |>
# tidyr::unnest(sim) |>
# dplyr::group_by(date) |>
# # compute quantiles for each day
# dplyr::summarise(
# gpp_q05 = quantile(gpp, 0.05, na.rm = TRUE),
# gpp_q50 = quantile(gpp, 0.5, na.rm = TRUE), # get median
# gpp_q95 = quantile(gpp, 0.95, na.rm = TRUE),
# gpp_pred_q05 = quantile(gpp_pred, 0.05, na.rm = TRUE),
# gpp_pred_q95 = quantile(gpp_pred, 0.95, na.rm = TRUE)
# )
## ----simulate_pmodel_runs, include = FALSE------------------------------------
# TODO: get rid of this and always fully run the vignettes
# fake output since calibration isn't run
# saveRDS(pmodel_runs, "files/sensitivity_analysis.Rmd__pmodel_runs.RDS")
pmodel_runs <- readRDS("files/sensitivity_analysis.Rmd__pmodel_runs.RDS")
## ----fig.width=7, fig.height=5------------------------------------------------
# Plot the credible intervals computed above
# for the first year only
data_to_plot <- pmodel_runs |>
# Plot only first year
dplyr::slice(1:365) |>
dplyr::left_join(
# Merge GPP validation data (first year)
p_model_validation$data[[1]][1:365, ] |>
dplyr::rename(gpp_obs = gpp),
by = "date")
# TODO: clean below
# plot_gpp_error <- ggplot(data = data_to_plot) +
# # geom_ribbon(
# # aes(ymin = gpp_pred_q05,
# # ymax = gpp_pred_q95,
# # x = date),
# # fill = 'grey40', alpha = 0.2) +
# geom_ribbon(
# aes(ymin = gpp_q05,
# ymax = gpp_q95,
# x = date),
# fill = 'blue', alpha = 0.5) +
# geom_ribbon(
# aes(ymin = gpp_pred_q05,
# ymax = gpp_pred_q95,
# x = date),
# fill = 'grey40', alpha = 0.2) +
# geom_line(
# aes(x = date,
# y = gpp_q50),
# colour = "grey40",
# alpha = 0.8) +
# # # Include observations in the plot
# geom_point(shape = 3,
# aes(x = date,
# y = gpp_obs),# color = "Observations"),
# alpha = 0.8) +
# theme_classic() +
# theme(panel.grid.major.y = element_line(),
# legend.position = "bottom") +
# labs(
# x = 'Date',
# y = expression(paste("GPP (g C m"^-2, "s"^-1, ")"))
# )
# plot_gpp_error
plot_gpp_error <- ggplot(data = data_to_plot) +
geom_ribbon(
aes(ymin = gpp_pred_q05,
ymax = gpp_pred_q95,
x = date,
fill = "Model uncertainty")) +
geom_ribbon(
aes(ymin = gpp_q05,
ymax = gpp_q95,
x = date,
fill = "Parameter uncertainty")) +
geom_line(
aes(x = date,
y = gpp_q50,
color = "Predictions")) +
# Include observations in the plot
geom_point(shape = 3,
aes(x = date,
y = gpp_obs,
color = "Observations")) +
theme_classic() +
theme(panel.grid.major.y = element_line(),
legend.position = "bottom") +
labs(
x = 'Date',
y = expression(paste("GPP (g C m"^-2, "s"^-1, ")"))
)
t_col <- function(color, percent = 50, name = NULL) {
# color = color name
# percent = % transparency
# name = an optional name for the color
## Get RGB values for named color
rgb.val <- col2rgb(color)
## Make new color using input color as base and alpha set by transparency
t.col <- rgb(rgb.val[1], rgb.val[2], rgb.val[3],
max = 255,
alpha = (100 - percent) * 255 / 100,
names = name)
## Save the color
invisible(t.col)
}
plot_gpp_error <- plot_gpp_error +
scale_color_manual(name = "",
breaks = c("Observations",
"Predictions",
"Non-calibrated predictions"),
values = c(t_col("black", 10),
t_col("#E69F00", 10),
t_col("#56B4E9", 10))) +
scale_fill_manual(name = "",
breaks = c("Model uncertainty",
"Parameter uncertainty"),
values = c(t_col("#E69F00", 60),
t_col("#009E73", 40)))
plot_gpp_error
#dir.create("./analysis/paper_results_files2")
## ----fig.width=7, fig.height=5, echo = FALSE, eval = FALSE--------------------
# #
# # Plot observed and predicted GPP, with a 95% confidence interval using err_gpp
# data_to_plot_MAP <- runread_pmodel_f(drivers = p_model_drivers,
# par = c(par_fixed, par_calib$par)) |>
# dplyr::select(sitename, data) |>
# tidyr::unnest(data) |>
# dplyr::slice(1:365)
#
# plot_gpp_error <- ggplot(data = data_to_plot_MAP) + # Plot only first year
# geom_ribbon(
# aes(ymin = gpp - 2*par_calib$par[4],
# ymax = gpp + 2*par_calib$par[4],
# x = date),
# fill = 'grey40', alpha = 0.2) +
# geom_line(
# aes(
# date,
# gpp
# ),
# colour = "grey40",
# alpha = 0.8
# ) +
# theme_classic() +
# theme(panel.grid.major.y = element_line()) +
# labs(
# x = 'Date',
# y = expression(paste("GPP (g C m"^-2, "s"^-1, ")"))
# )
#
# # Define GPP validation data (first year)
# validation_data <- p_model_validation$data[[1]][1:365, ]
#
# # Include observations in the plot
# plot_gpp_error +
# geom_point(
# data = validation_data,
# aes(
# date,
# gpp
# ),
# alpha = 0.8
# )