This function runs the EcoDiet model using a Markov chain Monte Carlo approximation through the 'jagsUI' package to provide an approximated distribution for the variables of interest.
Depending on the nb_iter entered, this function may take hours, or even days
to run. We advise you to first test whether your model is compiling properly with the by-default parameters,
as this should take 1-2 min to run depending on your data size.
To save time, this function can solicit several cores (if available) to parallelize chains. Note that progress bars won't be displayed if chains are parallelized.
A warning message is printed if the model has not converged in the end (if the Gelman-Rubin diagnostic of at least one variable is > 1.1). For each run, the default 'jagsUI' package messages summarize the '.txt' file used for the definition of the BUGS model, the configuration of the model (iteration, adaptation, burnin, thin rate), the time required to run the model, and main statistics for the variables.
You need to have run the preprocess_data and the write_model functions
before using this function, as their outputs are used as the inputs for run_model.
run_model(
model_file,
data,
inits = NULL,
run_param = "test",
variables_to_save = c("eta", "PI"),
parallelize = FALSE,
DIC.out = TRUE
)The file containing the BUGS definition of the EcoDiet model
output by the write_model function
The preprocessed data list output by the preprocess_data() function
A list containing the initial values of the variables. By default the initialisation values are NULL, which means that the chain initial values are drawn from the prior distributions.
A object that can be a list of the parameters to configure the JAGS model or a string acting as a shortcut characterizing the overall length of the run requested (e.g. "short" or "long"). If run_param is provided as a list, the user should provide at least nb_iter, i.e. the number of iterations to run (the more iterations, the better are the chances that the model will converge; very small by default to test if the model compiles properly), and nb_burnin, i.e. the number of burn-in steps to run (so that the variable approximations are not too influenced by the first initial random values). nb_thin, the thinning rate, is by default defined by the function. The number of adaptation steps nb_adapt can be specified but is not required (see jasgUI documentation for more details). If set manually, it should be at least set at 1000.
A vector of variable names defining the variables to output.
The number has a big number of variables but by default we only save the variables of interest
that are the trophic link probabilities eta and the diet proportions PI.
Only these saved variables are used to compute the Gelman-Rubin statistics that indicate whether
the model has converged or not.
Indicates whether chains should be parallelized using several cores. Recommended in case of complex models.
Indicates whether the DIC (Deviance Information Criterion) should be reported.
A MCMC output formatted as a jagsUI object.
preprocess_data to preprocess the data, and
write_model to define the model.
# \donttest{
realistic_biotracer_data <- read.csv(system.file("extdata", "realistic_biotracer_data.csv",
package = "EcoDiet"))
realistic_stomach_data <- read.csv(system.file("extdata", "realistic_stomach_data.csv",
package = "EcoDiet"))
data <- preprocess_data(biotracer_data = realistic_biotracer_data,
trophic_discrimination_factor = c(0.8, 3.4),
literature_configuration = FALSE,
stomach_data = realistic_stomach_data)
#> The model will investigate the following trophic links:
#> Bivalves Cod Crabs Detritus Phytoplankton Pout Sardine Shrimps
#> Bivalves 0 0 1 0 0 1 0 0
#> Cod 0 0 0 0 0 0 0 0
#> Crabs 0 1 0 0 0 1 0 0
#> Detritus 1 0 1 0 0 0 0 1
#> Phytoplankton 1 0 0 0 0 0 1 1
#> Pout 0 1 0 0 0 0 0 0
#> Sardine 0 1 0 0 0 0 0 0
#> Shrimps 0 1 1 0 0 1 1 0
#> Worms 0 1 1 0 0 1 0 1
#> Zooplankton 0 0 0 0 0 1 1 1
#> Worms Zooplankton
#> Bivalves 0 0
#> Cod 0 0
#> Crabs 0 0
#> Detritus 1 0
#> Phytoplankton 0 1
#> Pout 0 0
#> Sardine 0 0
#> Shrimps 0 0
#> Worms 0 0
#> Zooplankton 0 0
write_model(literature_configuration = FALSE)
mcmc_output <- run_model("EcoDiet_model.txt", data, run_param="test")
#>
#> Processing function input.......
#>
#> Done.
#>
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 316
#> Unobserved stochastic nodes: 125
#> Total graph size: 1104
#>
#> Initializing model
#>
#> Adaptive phase, 500 iterations x 3 chains
#> If no progress bar appears JAGS has decided not to adapt
#>
#>
#> Burn-in phase, 500 iterations x 3 chains
#>
#>
#> Sampling from joint posterior, 500 iterations x 3 chains
#>
#>
#> Calculating statistics.......
#>
#> Done.
#>
#> /!\ Convergence warning:
#> Out of the 51 variables, 14 variables have a Gelman-Rubin statistic > 1.1.
#> You may consider modifying the model run settings.
#> The variables with the poorest convergence are: PI[4,3], PI[10,8], PI[10,7], PI[8,6], PI[8,2], PI[6,2], PI[7,2], PI[9,8], PI[9,3], PI[8,7].
#> JAGS output for model 'EcoDiet_model.txt', generated by jagsUI.
#> Estimates based on 3 chains of 1000 iterations,
#> adaptation = 500 iterations (sufficient),
#> burn-in = 500 iterations and thin rate = 1,
#> yielding 1500 total samples from the joint posterior.
#> MCMC ran for 0.24 minutes at time 2024-07-10 10:07:40.642848.
#>
#> mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
#> eta[4,1] 0.648 0.084 0.481 0.648 0.800 FALSE 1 1.002 649
#> eta[5,1] 0.592 0.086 0.419 0.593 0.756 FALSE 1 1.001 1500
#> eta[3,2] 0.091 0.059 0.012 0.080 0.233 FALSE 1 1.000 1500
#> eta[6,2] 0.088 0.059 0.012 0.076 0.240 FALSE 1 1.011 200
#> eta[7,2] 0.425 0.101 0.233 0.425 0.624 FALSE 1 1.002 877
#> eta[8,2] 0.220 0.085 0.082 0.213 0.403 FALSE 1 1.013 153
#> eta[9,2] 0.734 0.092 0.536 0.742 0.888 FALSE 1 1.000 1500
#> eta[1,3] 0.389 0.087 0.225 0.384 0.565 FALSE 1 1.005 449
#> eta[4,3] 0.639 0.085 0.467 0.641 0.801 FALSE 1 1.019 107
#> eta[8,3] 0.808 0.069 0.664 0.813 0.924 FALSE 1 1.000 1500
#> eta[9,3] 0.811 0.069 0.657 0.819 0.927 FALSE 1 0.999 1500
#> eta[1,6] 0.126 0.058 0.036 0.119 0.257 FALSE 1 1.011 241
#> eta[3,6] 0.788 0.070 0.638 0.793 0.908 FALSE 1 1.000 1500
#> eta[8,6] 0.156 0.062 0.054 0.151 0.293 FALSE 1 1.005 404
#> eta[9,6] 0.969 0.029 0.893 0.978 0.999 FALSE 1 1.002 1385
#> eta[10,6] 0.876 0.056 0.751 0.884 0.963 FALSE 1 1.001 1500
#> eta[5,7] 0.491 0.087 0.335 0.487 0.667 FALSE 1 1.001 1161
#> eta[8,7] 0.970 0.029 0.893 0.979 0.999 FALSE 1 1.001 1500
#> eta[10,7] 0.225 0.076 0.096 0.220 0.379 FALSE 1 1.007 289
#> eta[4,8] 0.608 0.108 0.393 0.611 0.810 FALSE 1 1.000 1500
#> eta[5,8] 0.525 0.107 0.314 0.525 0.726 FALSE 1 1.003 540
#> eta[9,8] 0.100 0.065 0.014 0.089 0.251 FALSE 1 1.002 666
#> eta[10,8] 0.194 0.083 0.062 0.183 0.370 FALSE 1 1.005 384
#> eta[4,9] 0.948 0.050 0.815 0.964 0.998 FALSE 1 1.001 1034
#> eta[5,10] 0.953 0.045 0.836 0.967 0.999 FALSE 1 1.005 516
#> PI[4,1] 0.394 0.286 0.000 0.361 0.966 FALSE 1 1.014 138
#> PI[5,1] 0.606 0.286 0.034 0.639 1.000 FALSE 1 1.014 138
#> PI[3,2] 0.217 0.296 0.000 0.053 0.951 FALSE 1 1.059 60
#> PI[6,2] 0.167 0.231 0.000 0.024 0.766 FALSE 1 1.460 8
#> PI[7,2] 0.122 0.210 0.000 0.000 0.718 FALSE 1 1.364 11
#> PI[8,2] 0.284 0.335 0.000 0.124 0.998 FALSE 1 1.489 8
#> PI[9,2] 0.209 0.204 0.000 0.152 0.770 FALSE 1 1.032 281
#> PI[1,3] 0.163 0.232 0.000 0.040 0.795 FALSE 1 1.208 17
#> PI[4,3] 0.163 0.202 0.000 0.000 0.583 FALSE 1 2.878 4
#> PI[8,3] 0.277 0.210 0.000 0.253 0.719 FALSE 1 1.050 66
#> PI[9,3] 0.397 0.305 0.002 0.330 0.997 FALSE 1 1.291 11
#> PI[1,6] 0.069 0.133 0.000 0.001 0.493 FALSE 1 1.239 15
#> PI[3,6] 0.320 0.193 0.011 0.304 0.741 FALSE 1 1.069 36
#> PI[8,6] 0.051 0.120 0.000 0.001 0.467 FALSE 1 1.520 10
#> PI[9,6] 0.291 0.200 0.011 0.258 0.732 FALSE 1 1.003 1429
#> PI[10,6] 0.269 0.192 0.005 0.239 0.708 FALSE 1 1.033 75
#> PI[5,7] 0.328 0.180 0.000 0.340 0.653 FALSE 1 1.147 24
#> PI[8,7] 0.578 0.196 0.101 0.602 0.905 FALSE 1 1.279 13
#> PI[10,7] 0.094 0.234 0.000 0.000 0.882 FALSE 1 1.705 7
#> PI[4,8] 0.213 0.206 0.000 0.159 0.728 FALSE 1 1.021 143
#> PI[5,8] 0.301 0.216 0.001 0.283 0.777 FALSE 1 1.185 15
#> PI[9,8] 0.219 0.191 0.000 0.162 0.642 FALSE 1 1.331 10
#> PI[10,8] 0.266 0.271 0.000 0.169 0.857 FALSE 1 1.770 6
#> PI[4,9] 1.000 0.000 1.000 1.000 1.000 FALSE 1 NA 1
#> PI[5,10] 1.000 0.000 1.000 1.000 1.000 FALSE 1 NA 1
#> deviance 866.400 11.003 846.223 865.943 887.645 FALSE 1 1.019 106
#>
#> **WARNING** Some Rhat values could not be calculated.
#> **WARNING** Rhat values indicate convergence failure.
#> Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#> For each parameter, n.eff is a crude measure of effective sample size.
#>
#> overlap0 checks if 0 falls in the parameter's 95% credible interval.
#> f is the proportion of the posterior with the same sign as the mean;
#> i.e., our confidence that the parameter is positive or negative.
#>
#> DIC info: (pD = var(deviance)/2)
#> pD = 59.5 and DIC = 925.854
#> DIC is an estimate of expected predictive error (lower is better).
# }