Fit model to data

The objective of this example is demonstrate how to users can use memc_modfit to fit a MEMC model configuration to observational data. In this example we will be fitting the MIMCS model to the data from CITATION EXPERIMENT.

Set Up

# This assumes that the MEMC package has already been installed. 
library(MEMC)

# Used to visualize results
library(ggplot2)  

# set a theme to use in all the plots
theme_set(theme_bw() + theme(legend.title = element_blank()))

Start by loading the comparison data, included in the package data.

comp_data <- read.csv(system.file("example/Ultisol_control.csv", package = "MEMC"))

Set up the the MIMCS model with different initial pool sizes to be consistent with the experiment step up.

# Save a copy of the mimics model configuration. 
my_mod_config <- MEMC::MIMCS_model

# Update the MEND model initial state conditions to reflect the experimental set up.
my_mod_config$state <- c(POM=4.71, MOM=17.67, QOM=0, MB=0.52, DOM=0.148, 
                         EP=0.052, EM=0.052, IC=0, Tot=23.484)

Run the model with the original parameter values.

original_rslts <- solve_model(my_mod_config, time = 0:max(comp_data$time))

Use memc_modfit to fit the model to the observed data.

Description of inputs to MEMC function memc_modfit are as follows:

• config: a memc model configuration object, either one of the pre-built configurations listed in model_configs or created using configure_model
• comp_data: a data frame containing the comparison data that the model will be fit to

fit <- memc_modfit(config = my_mod_config, 
                        comp_data = comp_data, 
                        x = c(V_d=2, V_p=5, V_m=0.01))

fit$par
#>        V_d        V_p        V_m 
#>  2.9079702 26.3228643  0.4300596

Re run the model with the parameter values.

new_rslts <- solve_model(my_mod_config, params = fit$par, time = 0:max(comp_data$time)) 

# Format into a single data table to plot
original_rslts$par <- "original values"
new_rslts$par <- "fitted values"
out <- rbind(original_rslts, new_rslts)

Plot the model results vs the comparison data.

ggplot(data = out[out$variable == "IC", ]) +
  geom_line(aes(time, value, color = par)) +
  geom_point(data = comp_data, aes(time, IC)) +
  facet_wrap("variable", scales = "free") + 
  labs(x = "Time", y = "mg C/g soil", title = "MIMCS")