This function operates a diagnostic of the fit EcoDiet model.
A message is printed to provide the number of variables for which the Gelman-Rubin diagnostic exceeds specific thresholds (> 1.01, > 1.05, >1.1). The list of the 10 worst variables in terms of convergence is also printed.
You need to have run the run_model function before using this function.
The design of this function is substantially inspired from a function with a similar objective
in the MixSIAR package [(Stock et al. 2018)](https://doi.org/10.7717/peerj.5096), for which code is available online on
the [MixSIAR GitHub repository](https://github.com/brianstock/MixSIAR).
The diagnostic plots are generated using the ggmcmc package [(Fernández-i-Marín, 2016)](https://CRAN.R-project.org/package=ggmcmc).
diagnose_model(jags_output, var.to.diag = "all", save = FALSE, save_path = ".")The MCMC output summarized in the class jagsUI object
output by run_model function
The list of variables for which diagnostic plots should be produced and save. By default, this argument is "all" hence is run for all the variables.
Indicates whether diagnostic plots should be produced and saved.
The path indicating where to save the diagnostic plots.
A matrix containing the Gelman diagnostic for all the variables monitored by the run_model function (variables_to_save argument).
run_model to run 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,1], PI[5,1], PI[8,6], PI[9,8], PI[10,7], PI[8,7], PI[3,2], PI[1,6], PI[7,2], PI[6,2].
#> 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.246 minutes at time 2024-07-10 10:07:01.742541.
#>
#> mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
#> eta[4,1] 0.648 0.086 0.474 0.650 0.808 FALSE 1 1.024 88
#> eta[5,1] 0.591 0.091 0.402 0.598 0.759 FALSE 1 1.033 68
#> eta[3,2] 0.085 0.058 0.012 0.075 0.229 FALSE 1 1.007 631
#> eta[6,2] 0.085 0.054 0.011 0.078 0.219 FALSE 1 1.002 781
#> eta[7,2] 0.442 0.102 0.246 0.439 0.645 FALSE 1 1.005 343
#> eta[8,2] 0.219 0.086 0.081 0.208 0.410 FALSE 1 1.004 425
#> eta[9,2] 0.729 0.092 0.530 0.733 0.892 FALSE 1 1.009 217
#> eta[1,3] 0.393 0.086 0.236 0.389 0.566 FALSE 1 1.002 820
#> eta[4,3] 0.646 0.083 0.480 0.649 0.805 FALSE 1 1.002 755
#> eta[8,3] 0.810 0.069 0.660 0.818 0.920 FALSE 1 1.000 1500
#> eta[9,3] 0.800 0.071 0.639 0.807 0.917 FALSE 1 1.002 743
#> eta[1,6] 0.122 0.055 0.038 0.115 0.253 FALSE 1 1.001 1500
#> eta[3,6] 0.785 0.071 0.624 0.792 0.902 FALSE 1 1.000 1500
#> eta[8,6] 0.156 0.062 0.054 0.151 0.303 FALSE 1 1.004 428
#> eta[9,6] 0.970 0.030 0.890 0.979 0.999 FALSE 1 1.002 1500
#> eta[10,6] 0.879 0.055 0.759 0.885 0.965 FALSE 1 1.001 1500
#> eta[5,7] 0.484 0.090 0.311 0.488 0.643 FALSE 1 1.006 510
#> eta[8,7] 0.968 0.031 0.886 0.977 0.999 FALSE 1 1.001 1500
#> eta[10,7] 0.221 0.073 0.096 0.214 0.384 FALSE 1 1.001 1065
#> eta[4,8] 0.619 0.103 0.415 0.622 0.801 FALSE 1 1.006 298
#> eta[5,8] 0.520 0.103 0.319 0.524 0.716 FALSE 1 1.005 349
#> eta[9,8] 0.091 0.060 0.013 0.078 0.230 FALSE 1 1.008 361
#> eta[10,8] 0.210 0.085 0.071 0.201 0.393 FALSE 1 1.003 566
#> eta[4,9] 0.950 0.045 0.830 0.962 0.999 FALSE 1 1.002 1500
#> eta[5,10] 0.955 0.042 0.845 0.968 0.999 FALSE 1 0.999 1500
#> PI[4,1] 0.368 0.424 0.000 0.113 1.000 FALSE 1 2.366 4
#> PI[5,1] 0.632 0.424 0.000 0.887 1.000 FALSE 1 2.366 4
#> PI[3,2] 0.097 0.202 0.000 0.001 0.763 FALSE 1 1.468 10
#> PI[6,2] 0.060 0.156 0.000 0.000 0.611 FALSE 1 1.305 18
#> PI[7,2] 0.508 0.368 0.000 0.572 1.000 FALSE 1 1.373 9
#> PI[8,2] 0.208 0.313 0.000 0.017 0.989 FALSE 1 1.179 17
#> PI[9,2] 0.127 0.166 0.000 0.038 0.562 FALSE 1 1.095 27
#> PI[1,3] 0.204 0.235 0.000 0.114 0.783 FALSE 1 1.022 201
#> PI[4,3] 0.183 0.171 0.000 0.154 0.573 FALSE 1 1.013 247
#> PI[8,3] 0.313 0.207 0.006 0.294 0.753 FALSE 1 1.027 94
#> PI[9,3] 0.300 0.260 0.000 0.239 0.902 FALSE 1 1.021 144
#> PI[1,6] 0.026 0.104 0.000 0.000 0.439 FALSE 1 1.397 19
#> PI[3,6] 0.340 0.199 0.001 0.332 0.744 FALSE 1 1.024 86
#> PI[8,6] 0.033 0.086 0.000 0.000 0.298 FALSE 1 1.647 8
#> PI[9,6] 0.307 0.206 0.016 0.271 0.749 FALSE 1 1.011 185
#> PI[10,6] 0.294 0.201 0.004 0.269 0.727 FALSE 1 1.003 722
#> PI[5,7] 0.256 0.202 0.000 0.275 0.643 FALSE 1 1.257 12
#> PI[8,7] 0.678 0.236 0.171 0.667 1.000 FALSE 1 1.557 7
#> PI[10,7] 0.067 0.180 0.000 0.000 0.709 FALSE 1 1.578 9
#> PI[4,8] 0.239 0.217 0.000 0.210 0.710 FALSE 1 1.077 35
#> PI[5,8] 0.135 0.187 0.000 0.031 0.669 FALSE 1 1.162 19
#> PI[9,8] 0.026 0.069 0.000 0.000 0.285 FALSE 1 1.582 9
#> PI[10,8] 0.600 0.271 0.022 0.610 1.000 FALSE 1 1.230 13
#> 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 867.613 11.535 846.523 866.875 891.767 FALSE 1 1.015 176
#>
#> **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 = 65.9 and DIC = 933.468
#> DIC is an estimate of expected predictive error (lower is better).
Gelman_test <- diagnose_model(mcmc_output)
#>
#> ################################################################################
#> # Gelman-Rubin Diagnostic
#> ################################################################################
#> In EcoDiet, we recommend a that Gelman diagnostic is < 1.1
#> Out of 51 variables:
#> 25 > 1.01
#> 16 > 1.05
#> 14 > 1.1
#> The worst variables are:
#> Point est. Upper C.I.
#> PI[4,1] 2.366201 4.775357
#> PI[5,1] 2.366201 4.775357
#> PI[8,6] 1.647353 6.946107
#> PI[9,8] 1.581949 6.255799
#> PI[10,7] 1.577810 6.098783
#> PI[8,7] 1.557001 2.426211
#> PI[3,2] 1.467843 3.218520
#> PI[1,6] 1.397272 4.057678
#> PI[7,2] 1.373138 1.970375
#> PI[6,2] 1.305271 2.345151
Gelman_test
#> $gelman
#> Point est. Upper C.I.
#> eta[4,1] 1.0238218 1.0840029
#> eta[5,1] 1.0326479 1.1110296
#> eta[3,2] 1.0073030 1.0160698
#> eta[6,2] 1.0024676 1.0093722
#> eta[7,2] 1.0050101 1.0205926
#> eta[8,2] 1.0042858 1.0169150
#> eta[9,2] 1.0086269 1.0332372
#> eta[1,3] 1.0022756 1.0088487
#> eta[4,3] 1.0018874 1.0090014
#> eta[8,3] 0.9996050 1.0007175
#> eta[9,3] 1.0017490 1.0089648
#> eta[1,6] 1.0009908 1.0026964
#> eta[3,6] 0.9998243 1.0017870
#> eta[8,6] 1.0043536 1.0168919
#> eta[9,6] 1.0017632 1.0032294
#> eta[10,6] 1.0005100 1.0037226
#> eta[5,7] 1.0064825 1.0172002
#> eta[8,7] 1.0009986 1.0015033
#> eta[10,7] 1.0014542 1.0065071
#> eta[4,8] 1.0064703 1.0244491
#> eta[5,8] 1.0051852 1.0205364
#> eta[9,8] 1.0076261 1.0227224
#> eta[10,8] 1.0026900 1.0121669
#> eta[4,9] 1.0016584 1.0036672
#> eta[5,10] 0.9992420 0.9998608
#> PI[4,1] 2.3662010 4.7753575
#> PI[5,1] 2.3662010 4.7753575
#> PI[3,2] 1.4678431 3.2185198
#> PI[6,2] 1.3052707 2.3451506
#> PI[7,2] 1.3731382 1.9703747
#> PI[8,2] 1.1787074 1.5153410
#> PI[9,2] 1.0952846 1.2997314
#> PI[1,3] 1.0223346 1.0510041
#> PI[4,3] 1.0129505 1.0353616
#> PI[8,3] 1.0268470 1.0850507
#> PI[9,3] 1.0210912 1.0595104
#> PI[1,6] 1.3972717 4.0576783
#> PI[3,6] 1.0239406 1.0852896
#> PI[8,6] 1.6473532 6.9461073
#> PI[9,6] 1.0105640 1.0394081
#> PI[10,6] 1.0026421 1.0101049
#> PI[5,7] 1.2566339 1.7154494
#> PI[8,7] 1.5570009 2.4262107
#> PI[10,7] 1.5778100 6.0987835
#> PI[4,8] 1.0772212 1.2342001
#> PI[5,8] 1.1618971 1.5347317
#> PI[9,8] 1.5819491 6.2557988
#> PI[10,8] 1.2298548 1.6498196
#> deviance 1.0146832 1.0456138
#>
# }