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, 15 variables have a Gelman-Rubin statistic > 1.1.
#> You may consider modifying the model run settings.
#> The variables with the poorest convergence are: PI[1,6], PI[9,8], PI[1,3], PI[10,7], PI[10,8], PI[8,6], PI[8,2], PI[5,8], PI[8,7], PI[3,2].
#> 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.179 minutes at time 2026-09-20 07:45:11.547371.
#> 
#>              mean     sd    2.5%     50%   97.5% overlap0 f  Rhat n.eff
#> eta[4,1]    0.650  0.088   0.473   0.654   0.808    FALSE 1 1.003   695
#> eta[5,1]    0.592  0.090   0.408   0.592   0.760    FALSE 1 1.000  1500
#> eta[3,2]    0.089  0.057   0.014   0.078   0.234    FALSE 1 1.003   716
#> eta[6,2]    0.091  0.059   0.012   0.080   0.231    FALSE 1 1.002  1424
#> eta[7,2]    0.443  0.100   0.251   0.440   0.635    FALSE 1 1.000  1500
#> eta[8,2]    0.217  0.084   0.078   0.211   0.396    FALSE 1 1.015   132
#> eta[9,2]    0.736  0.090   0.536   0.740   0.887    FALSE 1 1.002  1067
#> eta[1,3]    0.389  0.087   0.230   0.386   0.566    FALSE 1 1.016   153
#> eta[4,3]    0.649  0.083   0.482   0.654   0.805    FALSE 1 0.999  1500
#> eta[8,3]    0.810  0.068   0.665   0.817   0.921    FALSE 1 1.002   724
#> eta[9,3]    0.806  0.071   0.651   0.813   0.926    FALSE 1 1.002  1263
#> eta[1,6]    0.124  0.058   0.034   0.116   0.256    FALSE 1 1.003   503
#> eta[3,6]    0.788  0.073   0.624   0.794   0.908    FALSE 1 1.001  1340
#> eta[8,6]    0.159  0.064   0.054   0.152   0.307    FALSE 1 1.004   459
#> eta[9,6]    0.970  0.030   0.891   0.979   0.999    FALSE 1 0.999  1500
#> eta[10,6]   0.876  0.056   0.743   0.884   0.966    FALSE 1 1.004   706
#> eta[5,7]    0.485  0.088   0.318   0.485   0.652    FALSE 1 1.002   779
#> eta[8,7]    0.967  0.032   0.882   0.976   0.999    FALSE 1 1.005  1500
#> eta[10,7]   0.225  0.074   0.093   0.222   0.377    FALSE 1 1.004   661
#> eta[4,8]    0.606  0.103   0.407   0.609   0.790    FALSE 1 1.008   241
#> eta[5,8]    0.511  0.108   0.302   0.513   0.717    FALSE 1 1.004   554
#> eta[9,8]    0.093  0.063   0.011   0.077   0.257    FALSE 1 0.999  1500
#> eta[10,8]   0.204  0.087   0.065   0.194   0.407    FALSE 1 1.004   383
#> eta[4,9]    0.948  0.049   0.819   0.963   0.999    FALSE 1 0.999  1500
#> eta[5,10]   0.953  0.048   0.813   0.968   0.999    FALSE 1 1.000  1500
#> PI[4,1]     0.273  0.280   0.000   0.182   0.975    FALSE 1 1.142    19
#> PI[5,1]     0.727  0.280   0.025   0.818   1.000    FALSE 1 1.142    19
#> PI[3,2]     0.106  0.218   0.000   0.004   0.854    FALSE 1 1.231    21
#> PI[6,2]     0.161  0.217   0.000   0.046   0.726    FALSE 1 1.081    58
#> PI[7,2]     0.309  0.277   0.000   0.243   0.897    FALSE 1 1.019   297
#> PI[8,2]     0.209  0.268   0.000   0.058   0.885    FALSE 1 1.324    10
#> PI[9,2]     0.215  0.192   0.001   0.170   0.722    FALSE 1 1.143    27
#> PI[1,3]     0.183  0.233   0.000   0.040   0.770    FALSE 1 1.397     9
#> PI[4,3]     0.190  0.149   0.000   0.172   0.521    FALSE 1 1.023   509
#> PI[8,3]     0.354  0.222   0.004   0.347   0.785    FALSE 1 1.197    14
#> PI[9,3]     0.272  0.251   0.000   0.219   0.834    FALSE 1 1.033    90
#> PI[1,6]     0.051  0.114   0.000   0.000   0.388    FALSE 1 1.755     7
#> PI[3,6]     0.318  0.198   0.007   0.306   0.723    FALSE 1 1.010   297
#> PI[8,6]     0.069  0.139   0.000   0.000   0.515    FALSE 1 1.366    12
#> PI[9,6]     0.299  0.204   0.016   0.266   0.762    FALSE 1 1.013   252
#> PI[10,6]    0.263  0.204   0.000   0.231   0.719    FALSE 1 1.005   367
#> PI[5,7]     0.269  0.176   0.000   0.281   0.609    FALSE 1 1.027   131
#> PI[8,7]     0.589  0.232   0.112   0.595   1.000    FALSE 1 1.259    12
#> PI[10,7]    0.141  0.231   0.000   0.004   0.765    FALSE 1 1.376    10
#> PI[4,8]     0.209  0.213   0.000   0.169   0.717    FALSE 1 1.172    19
#> PI[5,8]     0.192  0.231   0.000   0.092   0.786    FALSE 1 1.266    12
#> PI[9,8]     0.124  0.177   0.000   0.012   0.600    FALSE 1 1.480     8
#> PI[10,8]    0.475  0.270   0.014   0.459   0.995    FALSE 1 1.370     9
#> 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.495 11.809 845.786 865.616 891.512    FALSE 1 1.002  1446
#> 
#> **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 = 69.7 and DIC = 936.223 
#> 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, 17 variables have a Gelman-Rubin statistic > 1.1.
#> You may consider modifying the model run settings.
#> The variables with the poorest convergence are: PI[10,8], PI[9,8], PI[1,3], PI[7,2], PI[6,2], PI[5,8], PI[8,3], PI[9,3], PI[1,6], PI[8,7].
#> 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-09-20 07:45:25.335932.
#> 
#>              mean     sd    2.5%     50%   97.5% overlap0 f  Rhat n.eff
#> eta[4,1]    0.668  0.083   0.504   0.668   0.817    FALSE 1 1.005   334
#> eta[5,1]    0.605  0.087   0.438   0.607   0.768    FALSE 1 1.001  1500
#> eta[3,2]    0.357  0.087   0.198   0.350   0.551    FALSE 1 1.000  1500
#> eta[6,2]    0.359  0.083   0.201   0.357   0.528    FALSE 1 1.013   154
#> eta[7,2]    0.599  0.084   0.435   0.599   0.761    FALSE 1 1.008   228
#> eta[8,2]    0.442  0.086   0.281   0.443   0.611    FALSE 1 1.006   315
#> eta[9,2]    0.813  0.068   0.667   0.818   0.925    FALSE 1 1.002   862
#> eta[1,3]    0.518  0.080   0.355   0.519   0.677    FALSE 1 1.019   121
#> eta[4,3]    0.717  0.071   0.571   0.722   0.844    FALSE 1 1.000  1500
#> eta[8,3]    0.851  0.057   0.718   0.856   0.945    FALSE 1 1.000  1500
#> eta[9,3]    0.851  0.056   0.731   0.856   0.943    FALSE 1 1.003  1500
#> eta[1,6]    0.157  0.062   0.055   0.152   0.296    FALSE 1 1.001   946
#> eta[3,6]    0.787  0.069   0.645   0.792   0.907    FALSE 1 1.001   866
#> eta[8,6]    0.183  0.064   0.076   0.178   0.312    FALSE 1 1.000  1500
#> eta[9,6]    0.968  0.031   0.881   0.978   0.999    FALSE 1 1.007   602
#> eta[10,6]   0.880  0.055   0.747   0.888   0.964    FALSE 1 1.000  1500
#> eta[5,7]    0.559  0.083   0.394   0.558   0.715    FALSE 1 1.006   322
#> eta[8,7]    0.973  0.026   0.899   0.980   0.999    FALSE 1 1.004  1500
#> eta[10,7]   0.365  0.078   0.222   0.361   0.525    FALSE 1 1.000  1500
#> eta[4,8]    0.675  0.094   0.481   0.678   0.852    FALSE 1 1.002   905
#> eta[5,8]    0.596  0.094   0.406   0.597   0.773    FALSE 1 1.002   749
#> eta[9,8]    0.226  0.082   0.087   0.220   0.401    FALSE 1 1.001  1500
#> eta[10,8]   0.316  0.091   0.151   0.310   0.509    FALSE 1 1.003   760
#> eta[4,9]    0.953  0.043   0.838   0.964   0.999    FALSE 1 0.999  1500
#> eta[5,10]   0.958  0.039   0.854   0.970   0.999    FALSE 1 1.002  1500
#> PI[4,1]     0.258  0.348   0.000   0.053   1.000    FALSE 1 1.157    24
#> PI[5,1]     0.742  0.348   0.000   0.947   1.000    FALSE 1 1.157    24
#> PI[3,2]     0.127  0.206   0.000   0.006   0.772    FALSE 1 1.029   105
#> PI[6,2]     0.285  0.317   0.000   0.117   0.915    FALSE 1 1.558     7
#> PI[7,2]     0.467  0.392   0.000   0.507   0.998    FALSE 1 1.607     7
#> PI[8,2]     0.043  0.115   0.000   0.000   0.402    FALSE 1 1.266    19
#> PI[9,2]     0.078  0.111   0.000   0.022   0.371    FALSE 1 1.026    80
#> PI[1,3]     0.275  0.265   0.000   0.251   0.792    FALSE 1 1.667     6
#> PI[4,3]     0.122  0.124   0.000   0.089   0.443    FALSE 1 1.038   145
#> PI[8,3]     0.022  0.076   0.000   0.000   0.246    FALSE 1 1.421    16
#> PI[9,3]     0.582  0.264   0.102   0.578   0.997    FALSE 1 1.353     9
#> PI[1,6]     0.135  0.211   0.000   0.022   0.778    FALSE 1 1.296    13
#> PI[3,6]     0.332  0.252   0.000   0.313   0.809    FALSE 1 1.004   420
#> PI[8,6]     0.209  0.293   0.000   0.076   0.989    FALSE 1 1.149    19
#> PI[9,6]     0.182  0.218   0.000   0.080   0.712    FALSE 1 1.005   370
#> PI[10,6]    0.143  0.209   0.000   0.024   0.700    FALSE 1 1.009   541
#> PI[5,7]     0.076  0.107   0.000   0.025   0.375    FALSE 1 1.133    28
#> PI[8,7]     0.070  0.177   0.000   0.000   0.719    FALSE 1 1.274    20
#> PI[10,7]    0.854  0.198   0.179   0.929   1.000    FALSE 1 1.136    29
#> PI[4,8]     0.148  0.211   0.000   0.040   0.779    FALSE 1 1.271    12
#> PI[5,8]     0.174  0.215   0.000   0.074   0.699    FALSE 1 1.478     8
#> PI[9,8]     0.157  0.237   0.000   0.024   0.811    FALSE 1 1.727     6
#> PI[10,8]    0.521  0.374   0.000   0.627   0.998    FALSE 1 2.221     5
#> 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.907 11.323 845.944 865.213 890.194    FALSE 1 1.014   136
#> 
#> **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 = 63.2 and DIC = 929.145 
#> 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)