objective: test/demonstrate how the MEMC package can be used to calibrate a the SOM model with multiple output stream variables.

knitr::opts_chunk$set(echo = TRUE)
devtools::load_all("/Users/dorh012/Documents/2023/memc2/MEMC")
#library(MEMC)
library(dplyr)
library(ggplot2)

theme_set(theme_bw())

Generate Data

In this experiment we will be calibrating MEND to MEND output. To do this we want to run the model with non-default parameter values.

# Call the MEND model configuration file 
data(MEND_model)
# Define some new parameters, in theory the calibration process will return these parameter values! 
new_params <- c("V_p" = 28, "f_d" = 0.6)
full_comp_data <- solve_model(mod = MEND_model, time = 0:365, params = new_params)

Using IC Only

Fit MEND only using IC as comparison data.

full_comp_data %>% 
  dplyr::filter(variable == "IC") %>% 
  select(time, IC = value) ->
  IC_comp

fit1 <- memc_modfit(config = MEND_model, x = c("V_p" = 1, "f_d" = 0.4), lower = c(0, 0), upper = c(Inf, 1) , comp_data = IC_comp)
fit1$par
##        V_p        f_d 
## 20.4051574  0.8171424
out1 <- solve_model(mod = MEND_model, time = 0:365, params = fit1$par)
out1$name <- "IC only"

Let’s take a look at how well the IC values fit

ggplot() + 
  geom_point(data = IC_comp, aes(time, IC)) + 
  geom_line(data = out1 %>% filter(variable == "IC"), aes(time, value), color = "green") + 
  labs(title = "MEND fit using only IC")

Now consider the other variables

ggplot() + 
  geom_line(data = full_comp_data, aes(time, value, color = "truth")) + 
  geom_line(data = out1, aes(time, value, color = "fit with IC"), linetype = 2) + 
  facet_wrap("variable", scales = "free") + 
  scale_color_manual(values = c("truth" = "black","fit with IC" = "green"))

From the persepctive of IC, this model fit does well, however when we consider the other variables this model fit does not do well.

Using IC & DOM

Now try to fit MEND using multiple variables

full_comp_data %>% 
  dplyr::filter(variable %in% c("IC", "DOM")) %>% 
  select(time, variable, value) %>%  
  tidyr::spread(variable, value) -> 
  IC_DOM_comp

fit2 <- memc_modfit(config = MEND_model, x = c("V_p" = 1, "f_d" = 0.4), lower = c(0, 0), upper = c(Inf, 1) , comp_data = IC_DOM_comp)
fit2$par
##  V_p  f_d 
## 28.0  0.6
out2 <- solve_model(mod = MEND_model, time = 0:365, params = fit2$par)
out1$name <- "IC & DOM"

Let’s take a look at how well the IC and DOM values fit

ggplot() + 
  geom_point(data = IC_DOM_comp, aes(time, IC)) + 
  geom_line(data = out2 %>% filter(variable == "IC"), aes(time, value), color = "blue") + 
  labs(title = "MEND fit using IC & DOM")

ggplot() + 
  geom_point(data = IC_DOM_comp, aes(time, DOM)) + 
  geom_line(data = out2 %>% filter(variable == "DOM"), aes(time, value), color = "blue") + 
  labs(title = "MEND fit using IC & DOM")

ggplot() + 
  geom_line(data = full_comp_data, aes(time, value, color = "truth")) + 
  geom_line(data = out2, aes(time, value, color = "fit with IC & DOM"), linetype = 2) + 
  facet_wrap("variable", scales = "free") + 
  scale_color_manual(values = c("truth" = "black","fit with IC & DOM" = "blue"))

The optimization routine returns the correct values! And as expected now all the variables match!