EcoDiet explored on a realistic example

The introduction employed a simplistic expemple of food web to familiarize the user with the basic commands and options of the EcoDiet package. Here we will use a more realistic example (although still artificial!) to run the different EcoDiet configurations, compare their results and hence higlight the complementarity in the different data used.

The data corresponds to 10 trophic groups with stomach content data, and very distinct isotopic measures.

realistic_stomach_data_path <- system.file("extdata", "realistic_stomach_data.csv",
                                           package = "EcoDiet")
realistic_stomach_data <- read.csv(realistic_stomach_data_path)
knitr::kable(realistic_stomach_data)
X Cod Pout Sardine Shrimps Crabs Bivalves Worms Zooplankton Phytoplankton Detritus
Cod 0 0 0 0 0 0 0 0 0 0
Pout 1 0 0 0 0 0 0 0 0 0
Sardine 9 0 0 0 0 0 0 0 0 0
Shrimps 4 4 29 0 24 0 0 0 0 0
Crabs 1 24 0 0 0 0 0 0 0 0
Bivalves 0 3 0 0 11 0 0 0 0 0
Worms 16 30 0 1 24 0 0 0 0 0
Zooplankton 0 27 6 3 0 0 0 0 0 0
Phytoplankton 0 0 14 10 0 16 0 20 0 0
Detritus 0 0 0 12 19 18 18 0 0 0
full 21 30 29 19 29 27 18 20 0 0 
realistic_biotracer_data_path <- system.file("extdata", "realistic_biotracer_data.csv",
                                           package = "EcoDiet")
realistic_biotracer_data <- read.csv(realistic_biotracer_data_path)
knitr::kable(realistic_biotracer_data[c(1:3, 31:33, 61:63), ])
group d13C d15N
1 Cod -12.94144 19.18913
2 Cod -14.96070 20.23939
3 Cod -13.77822 19.48809
31 Pout -13.47127 18.57353
32 Pout -13.16888 17.58714
33 Pout -14.23085 17.38938
61 Sardine -14.56111 16.95231
62 Sardine -15.04729 17.15197
63 Sardine -14.63688 16.90906 
library(EcoDiet)

plot_data(biotracer_data = realistic_biotracer_data,
          stomach_data = realistic_stomach_data)

#> Warning: Use of `biotracer_data$group` is discouraged.
#> ℹ Use `group` instead.

Yes, we are aware that isotopic data is usually messier, but isn’t it a beautiful plot?

The configuration without literature data

We define the configuration we are in, and preprocess the data:

literature_configuration <- FALSE

data <- preprocess_data(biotracer_data = realistic_biotracer_data,
                        trophic_discrimination_factor = c(0.8, 3.4),
                        literature_configuration = literature_configuration,
                        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

In this configuration, priors are set for each trophic link identified as plausible by the user but the priors are not informed by literature data, and are thus uninformative:

plot_prior(data, literature_configuration)

The marginal prior distributions have different shape depending on the variables:

  • it is flat or uniform for \(\eta\), the probabilities that a trophic link exists (all the probabilities of existence are thus equiprobable),

  • the marginal distributions for each diet proportion \(\Pi\) are peaking at zero, although the joint distribution for \(\Pi\)s is a flat Dirichlet prior, because all the diet proportions must sum to one.

plot_prior(data, literature_configuration, pred = "Pout")

We define the model, and test if it compiles well with a few iterations and adaptation steps:

filename <- "mymodel.txt"
write_model(file.name = filename, literature_configuration = literature_configuration, print.model = F)
mcmc_output <- run_model(filename, 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[7,2], PI[10,8], PI[8,2], PI[4,8], PI[1,6], PI[10,7], PI[3,2], PI[6,2], PI[8,6], PI[9,8].
#> JAGS output for model 'mymodel.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.18 minutes at time 2026-10-14 09:09:19.488828.
#> 
#>              mean     sd    2.5%     50%   97.5% overlap0 f  Rhat n.eff
#> eta[4,1]    0.654  0.088   0.474   0.654   0.818    FALSE 1 1.002  1086
#> eta[5,1]    0.589  0.090   0.411   0.591   0.763    FALSE 1 1.000  1500
#> eta[3,2]    0.086  0.059   0.012   0.074   0.235    FALSE 1 1.006   478
#> eta[6,2]    0.090  0.059   0.012   0.077   0.234    FALSE 1 1.007   400
#> eta[7,2]    0.443  0.102   0.252   0.441   0.654    FALSE 1 1.018   112
#> eta[8,2]    0.217  0.083   0.079   0.210   0.389    FALSE 1 1.010   193
#> eta[9,2]    0.736  0.089   0.540   0.742   0.896    FALSE 1 1.004   528
#> eta[1,3]    0.384  0.087   0.221   0.381   0.559    FALSE 1 1.002   769
#> eta[4,3]    0.647  0.081   0.484   0.649   0.798    FALSE 1 1.004   479
#> eta[8,3]    0.807  0.070   0.650   0.814   0.923    FALSE 1 1.001  1156
#> eta[9,3]    0.805  0.070   0.652   0.811   0.920    FALSE 1 0.999  1500
#> eta[1,6]    0.126  0.059   0.037   0.117   0.258    FALSE 1 1.001  1250
#> eta[3,6]    0.785  0.072   0.631   0.791   0.903    FALSE 1 1.000  1500
#> eta[8,6]    0.154  0.063   0.052   0.144   0.294    FALSE 1 1.002  1078
#> eta[9,6]    0.968  0.031   0.886   0.977   0.999    FALSE 1 1.005  1001
#> eta[10,6]   0.877  0.054   0.758   0.884   0.961    FALSE 1 1.000  1500
#> eta[5,7]    0.488  0.085   0.329   0.485   0.656    FALSE 1 1.002   686
#> eta[8,7]    0.968  0.032   0.876   0.977   0.999    FALSE 1 1.010   401
#> eta[10,7]   0.226  0.074   0.102   0.220   0.385    FALSE 1 1.020   103
#> eta[4,8]    0.617  0.104   0.407   0.623   0.810    FALSE 1 1.002   725
#> eta[5,8]    0.523  0.105   0.321   0.518   0.728    FALSE 1 1.006   355
#> eta[9,8]    0.095  0.063   0.014   0.081   0.256    FALSE 1 1.002   827
#> eta[10,8]   0.207  0.089   0.067   0.196   0.404    FALSE 1 1.016   135
#> eta[4,9]    0.952  0.048   0.829   0.967   0.999    FALSE 1 1.001  1319
#> eta[5,10]   0.956  0.042   0.848   0.967   0.999    FALSE 1 1.004  1009
#> PI[4,1]     0.417  0.281   0.002   0.386   1.000    FALSE 1 1.013   196
#> PI[5,1]     0.583  0.281   0.000   0.614   0.998    FALSE 1 1.013   196
#> PI[3,2]     0.083  0.172   0.000   0.003   0.675    FALSE 1 1.519     9
#> PI[6,2]     0.138  0.236   0.000   0.002   0.792    FALSE 1 1.411    10
#> PI[7,2]     0.425  0.375   0.000   0.413   0.997    FALSE 1 2.545     4
#> PI[8,2]     0.178  0.296   0.000   0.001   0.935    FALSE 1 1.691     7
#> PI[9,2]     0.176  0.187   0.000   0.123   0.619    FALSE 1 1.038    58
#> PI[1,3]     0.096  0.198   0.000   0.000   0.709    FALSE 1 1.205    18
#> PI[4,3]     0.201  0.182   0.000   0.171   0.615    FALSE 1 1.019   107
#> PI[8,3]     0.291  0.229   0.000   0.274   0.777    FALSE 1 1.011   193
#> PI[9,3]     0.412  0.304   0.000   0.374   0.996    FALSE 1 1.012   207
#> PI[1,6]     0.030  0.079   0.000   0.000   0.284    FALSE 1 1.593     9
#> PI[3,6]     0.334  0.189   0.007   0.335   0.713    FALSE 1 1.005   502
#> PI[8,6]     0.010  0.030   0.000   0.000   0.082    FALSE 1 1.386    16
#> PI[9,6]     0.337  0.210   0.014   0.315   0.777    FALSE 1 1.020   126
#> PI[10,6]    0.289  0.195   0.013   0.256   0.721    FALSE 1 1.013   154
#> PI[5,7]     0.244  0.203   0.000   0.249   0.646    FALSE 1 1.169    16
#> PI[8,7]     0.549  0.234   0.083   0.565   0.998    FALSE 1 1.217    13
#> PI[10,7]    0.207  0.296   0.000   0.006   0.900    FALSE 1 1.546     7
#> PI[4,8]     0.340  0.367   0.000   0.207   1.000    FALSE 1 1.594     7
#> PI[5,8]     0.218  0.226   0.000   0.150   0.722    FALSE 1 1.107    25
#> PI[9,8]     0.023  0.087   0.000   0.000   0.366    FALSE 1 1.363    23
#> PI[10,8]    0.419  0.311   0.000   0.429   0.985    FALSE 1 1.835     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  867.391 11.069 848.486 866.563 891.123    FALSE 1 1.017   127
#> 
#> **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 = 60.4 and DIC = 927.755 
#> DIC is an estimate of expected predictive error (lower is better).

You should now try to run the model until it converges (it should take around half an hour to run, so we won’t do it in this vignette):

mcmc_output <- run_model(filename, data, run_param="normal", parallelize = T)

Here are the figures corresponding to the results that have converged:

plot_results(mcmc_output, data)

plot_results(mcmc_output, data, pred = "Pout")

You can also plot the results for specific prey if you want a clearer figure:

plot_results(mcmc_output, data, pred = "Pout", 
             variable = "PI", prey = c("Bivalves", "Worms"))

The configuration with literature data

We now change the configuration to add literature data to the model:

literature_configuration <- TRUE
realistic_literature_diets_path <- system.file("extdata", "realistic_literature_diets.csv",
                                               package = "EcoDiet")
realistic_literature_diets <- read.csv(realistic_literature_diets_path)
knitr::kable(realistic_literature_diets)
X Cod Pout Sardine Shrimps Crabs Bivalves Worms Zooplankton Phytoplankton Detritus
Cod 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0 0.0 0 0
Pout 0.4275065 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0 0.0 0 0
Sardine 0.3603675 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0 0.0 0 0
Shrimps 0.0300737 0.5295859 0.0002143 0.0000000 0.0082107 0.0000000 0.0 0.0 0 0
Crabs 0.1410430 0.3332779 0.0000000 0.0000000 0.0000000 0.0000000 0.0 0.0 0 0
Bivalves 0.0000000 0.0006130 0.0000000 0.0000000 0.3441081 0.0000000 0.0 0.0 0 0
Worms 0.0410093 0.1023676 0.0000000 0.0171336 0.4435377 0.0000000 0.0 0.0 0 0
Zooplankton 0.0000000 0.0341557 0.7381375 0.9121505 0.0000000 0.0000000 0.0 0.0 0 0
Phytoplankton 0.0000000 0.0000000 0.2616482 0.0000610 0.0000000 0.9966847 0.0 1.0 0 0
Detritus 0.0000000 0.0000000 0.0000000 0.0706550 0.2041434 0.0033153 1.0 0.0 0 0
pedigree 0.8000000 0.1000000 0.5000000 0.3000000 0.7000000 0.1000000 0.2 0.2 1 1 
data <- preprocess_data(biotracer_data = realistic_biotracer_data,
                        trophic_discrimination_factor = c(0.8, 3.4),
                        literature_configuration = literature_configuration,
                        stomach_data = realistic_stomach_data,
                        literature_diets = realistic_literature_diets,
                        nb_literature = 12,
                        literature_slope = 0.5)
#> 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

Now we see that the prior distributions are informed by the literature data:

  • when the literature diet input is > 0, the trophic link probabilities \(\eta\) are shifted toward one. Here this is the case for all prey but we could imagine that the user identify a species as a plausible prey whereas it has not been observed being consumed by the predator in the literature. In that case, the literature diet of 0 would drive \(\eta\) toward 0.

  • the average prior for the diet proportions \(\Pi\) is directly the literature diet input.

plot_prior(data, literature_configuration)

plot_prior(data, literature_configuration, pred = "Pout")

Again, we verify that the model compiles well:

filename <- "mymodel_literature.txt"
write_model(file.name = filename, literature_configuration = literature_configuration, print.model = F)
mcmc_output <- run_model(filename, 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: 1594
#> 
#> 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, 11 variables have a Gelman-Rubin statistic > 1.1.
#> You may consider modifying the model run settings.
#> The variables with the poorest convergence are: PI[8,7], PI[8,6], PI[10,7], PI[8,3], PI[1,6], PI[7,2], PI[5,8], PI[10,8], PI[9,6], PI[3,2].
#> JAGS output for model 'mymodel_literature.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.181 minutes at time 2026-10-14 09:09:33.709737.
#> 
#>              mean     sd    2.5%     50%   97.5% overlap0 f  Rhat n.eff
#> eta[4,1]    0.665  0.085   0.483   0.666   0.813    FALSE 1 1.001  1105
#> eta[5,1]    0.603  0.086   0.420   0.608   0.759    FALSE 1 1.000  1500
#> eta[3,2]    0.367  0.083   0.213   0.366   0.543    FALSE 1 1.001  1500
#> eta[6,2]    0.355  0.082   0.204   0.353   0.521    FALSE 1 1.001  1132
#> eta[7,2]    0.603  0.084   0.426   0.607   0.763    FALSE 1 1.003   569
#> eta[8,2]    0.452  0.086   0.292   0.451   0.622    FALSE 1 1.009   201
#> eta[9,2]    0.813  0.069   0.663   0.821   0.922    FALSE 1 1.000  1500
#> eta[1,3]    0.523  0.077   0.371   0.524   0.677    FALSE 1 1.003   599
#> eta[4,3]    0.723  0.070   0.574   0.726   0.847    FALSE 1 0.999  1500
#> eta[8,3]    0.850  0.056   0.724   0.857   0.939    FALSE 1 1.002   764
#> eta[9,3]    0.854  0.054   0.729   0.860   0.944    FALSE 1 1.002  1012
#> eta[1,6]    0.154  0.061   0.055   0.148   0.293    FALSE 1 1.003   781
#> eta[3,6]    0.789  0.070   0.639   0.792   0.910    FALSE 1 1.004   598
#> eta[8,6]    0.182  0.063   0.076   0.175   0.322    FALSE 1 1.014   149
#> eta[9,6]    0.970  0.029   0.895   0.979   0.999    FALSE 1 1.001  1500
#> eta[10,6]   0.879  0.056   0.751   0.886   0.965    FALSE 1 1.001   953
#> eta[5,7]    0.564  0.080   0.401   0.564   0.725    FALSE 1 1.004   463
#> eta[8,7]    0.972  0.026   0.903   0.980   1.000    FALSE 1 1.004  1432
#> eta[10,7]   0.362  0.077   0.224   0.360   0.526    FALSE 1 1.004   401
#> eta[4,8]    0.671  0.094   0.474   0.676   0.836    FALSE 1 1.001  1013
#> eta[5,8]    0.592  0.095   0.402   0.595   0.772    FALSE 1 1.000  1500
#> eta[9,8]    0.227  0.085   0.085   0.219   0.412    FALSE 1 1.001  1123
#> eta[10,8]   0.318  0.094   0.150   0.312   0.519    FALSE 1 1.001  1500
#> eta[4,9]    0.955  0.043   0.839   0.967   0.999    FALSE 1 1.007   395
#> eta[5,10]   0.961  0.038   0.861   0.972   0.999    FALSE 1 1.006  1500
#> PI[4,1]     0.260  0.361   0.000   0.017   1.000    FALSE 1 1.057    54
#> PI[5,1]     0.740  0.361   0.000   0.983   1.000    FALSE 1 1.057    54
#> PI[3,2]     0.250  0.253   0.000   0.166   0.834    FALSE 1 1.163    18
#> PI[6,2]     0.204  0.262   0.000   0.054   0.863    FALSE 1 1.046    71
#> PI[7,2]     0.402  0.294   0.000   0.420   0.906    FALSE 1 1.248    13
#> PI[8,2]     0.052  0.075   0.000   0.023   0.275    FALSE 1 1.077   233
#> PI[9,2]     0.091  0.093   0.000   0.063   0.324    FALSE 1 1.013   180
#> PI[1,3]     0.330  0.224   0.001   0.322   0.798    FALSE 1 1.066    42
#> PI[4,3]     0.109  0.115   0.000   0.080   0.393    FALSE 1 1.028   142
#> PI[8,3]     0.074  0.144   0.000   0.013   0.494    FALSE 1 1.563     9
#> PI[9,3]     0.487  0.246   0.030   0.472   0.993    FALSE 1 1.064    51
#> PI[1,6]     0.098  0.187   0.000   0.004   0.681    FALSE 1 1.330    13
#> PI[3,6]     0.395  0.277   0.000   0.380   0.924    FALSE 1 1.045    50
#> PI[8,6]     0.187  0.297   0.000   0.012   0.975    FALSE 1 1.885     6
#> PI[9,6]     0.283  0.264   0.000   0.236   0.847    FALSE 1 1.187    15
#> PI[10,6]    0.037  0.091   0.000   0.001   0.310    FALSE 1 1.009  1473
#> PI[5,7]     0.131  0.192   0.000   0.046   0.625    FALSE 1 1.137    24
#> PI[8,7]     0.143  0.221   0.000   0.018   0.740    FALSE 1 2.156     5
#> PI[10,7]    0.726  0.326   0.000   0.860   1.000    FALSE 1 1.767     6
#> PI[4,8]     0.161  0.202   0.000   0.063   0.681    FALSE 1 1.045    84
#> PI[5,8]     0.083  0.158   0.000   0.003   0.569    FALSE 1 1.245    15
#> PI[9,8]     0.070  0.155   0.000   0.003   0.588    FALSE 1 1.097    49
#> PI[10,8]    0.686  0.270   0.067   0.738   1.000    FALSE 1 1.213    14
#> 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  865.753 10.945 845.952 865.146 888.261    FALSE 1 1.005   659
#> 
#> **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.8 and DIC = 925.55 
#> DIC is an estimate of expected predictive error (lower is better).

You should now try to run the model until it converges (it should take around half an hour to run, so we won’t do it in this vignette):

mcmc_output <- run_model(filename, data, run_param=list(nb_iter=100000, nb_burnin=50000, nb_thin=50, nb_adapt=50000), parallelize = T)

Here are the figures corresponding to the results that have converged:

plot_results(mcmc_output, data)

plot_results(mcmc_output, data, pred = "Pout")

You can save the figures as PNG using:

plot_results(mcmc_output, data, pred = "Pout", save = TRUE, save_path = ".")

Last, if you want to explore further in detail the a posteriori distribution of your parameters \(\Pi\) and \(\eta\), you can run the following code line, which will store the values for all iterations into a data frame.

reshape_mcmc(mcmc_output, data)