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, 13 variables have a Gelman-Rubin statistic > 1.1.
#> You may consider modifying the model run settings.
#> The variables with the poorest convergence are: PI[9,8], PI[8,6], PI[10,7], PI[6,2], PI[3,6], PI[8,7], PI[7,2], PI[10,8], PI[1,3], PI[9,3].
#> 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.182 minutes at time 2026-06-06 09:16:03.609788.
#> 
#>              mean     sd    2.5%     50%   97.5% overlap0 f  Rhat n.eff
#> eta[4,1]    0.654  0.085   0.482   0.654   0.810    FALSE 1 1.003  1039
#> eta[5,1]    0.588  0.088   0.411   0.590   0.748    FALSE 1 1.001   879
#> eta[3,2]    0.091  0.060   0.012   0.080   0.242    FALSE 1 1.000  1355
#> eta[6,2]    0.087  0.055   0.012   0.078   0.225    FALSE 1 1.001  1458
#> eta[7,2]    0.437  0.103   0.234   0.436   0.640    FALSE 1 1.005   367
#> eta[8,2]    0.221  0.084   0.080   0.213   0.404    FALSE 1 1.008   328
#> eta[9,2]    0.730  0.091   0.535   0.738   0.879    FALSE 1 1.004   486
#> eta[1,3]    0.386  0.086   0.232   0.386   0.556    FALSE 1 1.002   604
#> eta[4,3]    0.646  0.083   0.481   0.649   0.802    FALSE 1 1.007   283
#> eta[8,3]    0.809  0.068   0.666   0.817   0.925    FALSE 1 1.000  1500
#> eta[9,3]    0.814  0.068   0.661   0.820   0.924    FALSE 1 1.001  1500
#> eta[1,6]    0.122  0.059   0.033   0.115   0.259    FALSE 1 0.999  1500
#> eta[3,6]    0.782  0.072   0.630   0.789   0.902    FALSE 1 1.001  1099
#> eta[8,6]    0.155  0.062   0.054   0.151   0.291    FALSE 1 1.005   392
#> eta[9,6]    0.968  0.029   0.888   0.977   0.999    FALSE 1 1.003  1500
#> eta[10,6]   0.875  0.055   0.750   0.882   0.962    FALSE 1 1.004   566
#> eta[5,7]    0.487  0.088   0.324   0.487   0.657    FALSE 1 1.001   906
#> eta[8,7]    0.968  0.031   0.888   0.977   0.999    FALSE 1 1.003  1500
#> eta[10,7]   0.234  0.076   0.105   0.226   0.396    FALSE 1 1.006   413
#> eta[4,8]    0.609  0.105   0.392   0.613   0.798    FALSE 1 0.999  1500
#> eta[5,8]    0.513  0.105   0.307   0.515   0.716    FALSE 1 1.008   328
#> eta[9,8]    0.091  0.061   0.012   0.077   0.243    FALSE 1 1.006   662
#> eta[10,8]   0.199  0.085   0.057   0.192   0.386    FALSE 1 1.025    84
#> eta[4,9]    0.950  0.047   0.824   0.963   0.999    FALSE 1 1.002  1500
#> eta[5,10]   0.953  0.045   0.836   0.965   0.999    FALSE 1 1.001  1500
#> PI[4,1]     0.419  0.285   0.006   0.378   1.000    FALSE 1 1.014   185
#> PI[5,1]     0.581  0.285   0.000   0.622   0.994    FALSE 1 1.014   185
#> PI[3,2]     0.188  0.265   0.000   0.028   0.876    FALSE 1 1.040    89
#> PI[6,2]     0.090  0.197   0.000   0.003   0.765    FALSE 1 1.214    23
#> PI[7,2]     0.381  0.324   0.000   0.344   0.980    FALSE 1 1.188    15
#> PI[8,2]     0.217  0.270   0.000   0.096   0.975    FALSE 1 1.056    53
#> PI[9,2]     0.124  0.132   0.000   0.092   0.454    FALSE 1 1.105    25
#> PI[1,3]     0.145  0.218   0.000   0.002   0.700    FALSE 1 1.181    17
#> PI[4,3]     0.175  0.166   0.000   0.145   0.540    FALSE 1 1.113    22
#> PI[8,3]     0.226  0.189   0.000   0.182   0.661    FALSE 1 1.101    27
#> PI[9,3]     0.454  0.278   0.008   0.434   0.995    FALSE 1 1.120    23
#> PI[1,6]     0.009  0.030   0.000   0.000   0.083    FALSE 1 1.097   347
#> PI[3,6]     0.316  0.215   0.000   0.322   0.744    FALSE 1 1.201    14
#> PI[8,6]     0.080  0.173   0.000   0.000   0.617    FALSE 1 1.806     6
#> PI[9,6]     0.330  0.211   0.023   0.303   0.798    FALSE 1 1.010   325
#> PI[10,6]    0.266  0.205   0.000   0.239   0.710    FALSE 1 1.004  1379
#> PI[5,7]     0.242  0.180   0.000   0.249   0.590    FALSE 1 1.099    26
#> PI[8,7]     0.470  0.242   0.027   0.497   0.939    FALSE 1 1.196    14
#> PI[10,7]    0.288  0.320   0.000   0.140   0.924    FALSE 1 1.266    12
#> PI[4,8]     0.175  0.182   0.000   0.123   0.595    FALSE 1 1.064    87
#> PI[5,8]     0.221  0.211   0.000   0.182   0.762    FALSE 1 1.043   289
#> PI[9,8]     0.115  0.145   0.000   0.037   0.502    FALSE 1 1.832     6
#> PI[10,8]    0.489  0.265   0.007   0.493   0.997    FALSE 1 1.187    15
#> 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.683 11.334 845.458 865.331 890.262    FALSE 1 1.010   217
#> 
#> **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.7 and DIC = 929.406 
#> 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, 16 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,7], PI[8,7], PI[1,3], PI[4,1], PI[5,1], PI[5,8], PI[10,8], PI[9,3], PI[8,2], PI[6,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.182 minutes at time 2026-06-06 09:16:17.649832.
#> 
#>              mean     sd    2.5%     50%   97.5% overlap0 f  Rhat n.eff
#> eta[4,1]    0.667  0.085   0.489   0.668   0.826    FALSE 1 1.006   303
#> eta[5,1]    0.604  0.086   0.434   0.607   0.770    FALSE 1 1.001   942
#> eta[3,2]    0.359  0.084   0.206   0.355   0.527    FALSE 1 1.004   425
#> eta[6,2]    0.349  0.081   0.205   0.342   0.521    FALSE 1 1.006   330
#> eta[7,2]    0.615  0.082   0.449   0.619   0.767    FALSE 1 1.002  1106
#> eta[8,2]    0.447  0.085   0.286   0.445   0.613    FALSE 1 1.008   331
#> eta[9,2]    0.818  0.066   0.680   0.824   0.931    FALSE 1 1.000  1500
#> eta[1,3]    0.515  0.079   0.363   0.513   0.669    FALSE 1 1.024    95
#> eta[4,3]    0.721  0.072   0.569   0.724   0.847    FALSE 1 1.000  1500
#> eta[8,3]    0.847  0.054   0.723   0.852   0.938    FALSE 1 1.002   686
#> eta[9,3]    0.848  0.056   0.726   0.855   0.938    FALSE 1 1.007   322
#> eta[1,6]    0.157  0.065   0.055   0.150   0.303    FALSE 1 1.001  1365
#> eta[3,6]    0.788  0.072   0.636   0.794   0.910    FALSE 1 1.002   995
#> eta[8,6]    0.191  0.069   0.080   0.186   0.344    FALSE 1 1.000  1500
#> eta[9,6]    0.971  0.028   0.893   0.980   0.999    FALSE 1 1.006   482
#> eta[10,6]   0.879  0.055   0.759   0.886   0.964    FALSE 1 1.006   400
#> eta[5,7]    0.563  0.079   0.404   0.567   0.709    FALSE 1 0.999  1500
#> eta[8,7]    0.973  0.026   0.905   0.981   0.999    FALSE 1 1.000  1500
#> eta[10,7]   0.351  0.076   0.213   0.346   0.504    FALSE 1 1.026    82
#> eta[4,8]    0.673  0.090   0.488   0.677   0.837    FALSE 1 1.003   507
#> eta[5,8]    0.593  0.099   0.394   0.596   0.776    FALSE 1 1.000  1500
#> eta[9,8]    0.232  0.084   0.093   0.225   0.426    FALSE 1 1.000  1500
#> eta[10,8]   0.313  0.093   0.149   0.310   0.501    FALSE 1 1.005   342
#> eta[4,9]    0.954  0.044   0.836   0.967   0.999    FALSE 1 1.001  1236
#> eta[5,10]   0.958  0.039   0.856   0.970   0.999    FALSE 1 1.009   332
#> PI[4,1]     0.335  0.436   0.000   0.009   1.000    FALSE 1 1.822     6
#> PI[5,1]     0.665  0.436   0.000   0.991   1.000    FALSE 1 1.822     6
#> PI[3,2]     0.124  0.189   0.000   0.003   0.654    FALSE 1 1.206    18
#> PI[6,2]     0.080  0.160   0.000   0.002   0.559    FALSE 1 1.315    13
#> PI[7,2]     0.631  0.272   0.021   0.681   0.996    FALSE 1 1.100    25
#> PI[8,2]     0.046  0.104   0.000   0.001   0.366    FALSE 1 1.325    16
#> PI[9,2]     0.119  0.146   0.000   0.059   0.514    FALSE 1 1.142    23
#> PI[1,3]     0.168  0.262   0.000   0.000   0.920    FALSE 1 2.371     4
#> PI[4,3]     0.171  0.164   0.000   0.139   0.585    FALSE 1 1.077    36
#> PI[8,3]     0.096  0.150   0.000   0.018   0.522    FALSE 1 1.095    38
#> PI[9,3]     0.564  0.303   0.000   0.600   0.999    FALSE 1 1.536     7
#> PI[1,6]     0.015  0.058   0.000   0.000   0.162    FALSE 1 1.130    59
#> PI[3,6]     0.338  0.271   0.000   0.312   0.898    FALSE 1 1.026   134
#> PI[8,6]     0.283  0.298   0.000   0.154   0.942    FALSE 1 1.110    32
#> PI[9,6]     0.233  0.239   0.000   0.151   0.792    FALSE 1 1.058    44
#> PI[10,6]    0.131  0.207   0.000   0.015   0.747    FALSE 1 1.095    35
#> PI[5,7]     0.224  0.229   0.000   0.176   0.785    FALSE 1 1.279    12
#> PI[8,7]     0.440  0.375   0.000   0.527   1.000    FALSE 1 2.752     4
#> PI[10,7]    0.336  0.438   0.000   0.001   1.000    FALSE 1 4.486     3
#> PI[4,8]     0.140  0.196   0.000   0.030   0.677    FALSE 1 1.016   159
#> PI[5,8]     0.181  0.260   0.000   0.043   0.955    FALSE 1 1.745     6
#> PI[9,8]     0.359  0.283   0.000   0.318   0.963    FALSE 1 1.102    31
#> PI[10,8]    0.320  0.352   0.000   0.093   0.940    FALSE 1 1.598     7
#> 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  868.073 11.643 846.901 867.293 892.903    FALSE 1 1.019   108
#> 
#> **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 = 66.6 and DIC = 934.677 
#> 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)