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, 12 variables have a Gelman-Rubin statistic > 1.1.
#> You may consider modifying the model run settings.
#> The variables with the poorest convergence are: PI[6,2], PI[9,8], PI[5,8], PI[1,6], PI[10,8], PI[3,2], PI[7,2], PI[10,7], PI[8,6], PI[8,7].
#> 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.197 minutes at time 2026-06-22 09:16:11.23618.
#> 
#>              mean     sd    2.5%     50%   97.5% overlap0 f  Rhat n.eff
#> eta[4,1]    0.653  0.085   0.491   0.655   0.809    FALSE 1 1.010   184
#> eta[5,1]    0.594  0.088   0.419   0.595   0.755    FALSE 1 1.000  1500
#> eta[3,2]    0.092  0.060   0.010   0.080   0.236    FALSE 1 1.019   136
#> eta[6,2]    0.093  0.059   0.013   0.082   0.238    FALSE 1 1.031    77
#> eta[7,2]    0.440  0.102   0.249   0.441   0.641    FALSE 1 1.013   158
#> eta[8,2]    0.216  0.083   0.079   0.208   0.392    FALSE 1 1.003   537
#> eta[9,2]    0.739  0.087   0.558   0.744   0.896    FALSE 1 1.000  1500
#> eta[1,3]    0.389  0.086   0.232   0.385   0.560    FALSE 1 1.000  1500
#> eta[4,3]    0.648  0.086   0.472   0.651   0.801    FALSE 1 1.000  1500
#> eta[8,3]    0.810  0.070   0.653   0.816   0.927    FALSE 1 1.001  1500
#> eta[9,3]    0.808  0.069   0.653   0.814   0.922    FALSE 1 1.000  1427
#> eta[1,6]    0.124  0.056   0.036   0.118   0.253    FALSE 1 1.008   291
#> eta[3,6]    0.782  0.070   0.637   0.786   0.901    FALSE 1 1.000  1500
#> eta[8,6]    0.156  0.062   0.054   0.150   0.295    FALSE 1 1.001  1500
#> eta[9,6]    0.969  0.029   0.887   0.978   0.999    FALSE 1 1.006   594
#> eta[10,6]   0.878  0.059   0.737   0.886   0.968    FALSE 1 1.002  1500
#> eta[5,7]    0.485  0.088   0.320   0.488   0.655    FALSE 1 1.000  1500
#> eta[8,7]    0.968  0.031   0.891   0.978   0.999    FALSE 1 1.001  1500
#> eta[10,7]   0.229  0.073   0.106   0.222   0.386    FALSE 1 1.007   272
#> eta[4,8]    0.609  0.103   0.403   0.611   0.799    FALSE 1 1.005   383
#> eta[5,8]    0.524  0.105   0.321   0.522   0.733    FALSE 1 1.006   427
#> eta[9,8]    0.093  0.060   0.012   0.081   0.234    FALSE 1 1.000  1500
#> eta[10,8]   0.189  0.082   0.058   0.181   0.369    FALSE 1 1.012   185
#> eta[4,9]    0.950  0.046   0.828   0.964   0.998    FALSE 1 1.003   739
#> eta[5,10]   0.955  0.043   0.842   0.967   0.999    FALSE 1 1.002  1500
#> PI[4,1]     0.308  0.286   0.000   0.252   0.967    FALSE 1 1.031    77
#> PI[5,1]     0.692  0.286   0.033   0.748   1.000    FALSE 1 1.031    77
#> PI[3,2]     0.174  0.226   0.000   0.047   0.764    FALSE 1 1.468     8
#> PI[6,2]     0.160  0.286   0.000   0.003   0.902    FALSE 1 2.446     4
#> PI[7,2]     0.335  0.300   0.000   0.317   0.930    FALSE 1 1.392     9
#> PI[8,2]     0.164  0.229   0.000   0.033   0.790    FALSE 1 1.094    29
#> PI[9,2]     0.167  0.175   0.000   0.115   0.590    FALSE 1 1.036   121
#> PI[1,3]     0.115  0.200   0.000   0.001   0.715    FALSE 1 1.097    68
#> PI[4,3]     0.209  0.167   0.000   0.194   0.568    FALSE 1 1.028   126
#> PI[8,3]     0.313  0.219   0.001   0.293   0.791    FALSE 1 1.096    28
#> PI[9,3]     0.362  0.259   0.000   0.338   0.900    FALSE 1 1.025   104
#> PI[1,6]     0.040  0.094   0.000   0.000   0.332    FALSE 1 1.672     8
#> PI[3,6]     0.304  0.196   0.000   0.300   0.685    FALSE 1 1.055    43
#> PI[8,6]     0.037  0.084   0.000   0.000   0.298    FALSE 1 1.180    21
#> PI[9,6]     0.325  0.205   0.022   0.298   0.760    FALSE 1 1.025    94
#> PI[10,6]    0.293  0.201   0.008   0.267   0.738    FALSE 1 1.016   138
#> PI[5,7]     0.240  0.202   0.000   0.240   0.651    FALSE 1 1.122    21
#> PI[8,7]     0.509  0.243   0.042   0.543   0.936    FALSE 1 1.152    18
#> PI[10,7]    0.250  0.317   0.000   0.078   0.941    FALSE 1 1.268    12
#> PI[4,8]     0.225  0.240   0.000   0.157   0.842    FALSE 1 1.108    25
#> PI[5,8]     0.374  0.349   0.000   0.270   1.000    FALSE 1 1.813     6
#> PI[9,8]     0.118  0.191   0.000   0.006   0.657    FALSE 1 1.992     5
#> PI[10,8]    0.284  0.281   0.000   0.203   0.879    FALSE 1 1.489     8
#> 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.644 11.421 847.294 867.314 891.195    FALSE 1 1.006   304
#> 
#> **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 = 64.9 and DIC = 932.52 
#> 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,8], PI[8,7], PI[10,7], PI[1,3], PI[9,8], PI[9,3], PI[5,8], PI[4,8], PI[7,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.199 minutes at time 2026-06-22 09:16:26.198494.
#> 
#>              mean     sd    2.5%     50%   97.5% overlap0 f  Rhat n.eff
#> eta[4,1]    0.675  0.082   0.503   0.681   0.822    FALSE 1 1.001  1500
#> eta[5,1]    0.609  0.088   0.435   0.612   0.773    FALSE 1 1.000  1500
#> eta[3,2]    0.357  0.083   0.201   0.354   0.528    FALSE 1 1.001  1138
#> eta[6,2]    0.357  0.081   0.207   0.355   0.522    FALSE 1 1.018   116
#> eta[7,2]    0.605  0.085   0.432   0.607   0.759    FALSE 1 1.011   211
#> eta[8,2]    0.450  0.087   0.294   0.448   0.621    FALSE 1 1.001  1307
#> eta[9,2]    0.816  0.064   0.678   0.823   0.924    FALSE 1 1.013   173
#> eta[1,3]    0.518  0.078   0.360   0.521   0.669    FALSE 1 1.006   324
#> eta[4,3]    0.722  0.071   0.570   0.727   0.845    FALSE 1 1.002   801
#> eta[8,3]    0.849  0.056   0.722   0.856   0.942    FALSE 1 1.003  1500
#> eta[9,3]    0.850  0.058   0.712   0.859   0.942    FALSE 1 1.003  1270
#> eta[1,6]    0.159  0.063   0.054   0.152   0.296    FALSE 1 1.001  1321
#> eta[3,6]    0.790  0.069   0.643   0.797   0.911    FALSE 1 1.004   721
#> eta[8,6]    0.189  0.066   0.073   0.182   0.329    FALSE 1 1.004   732
#> eta[9,6]    0.969  0.029   0.891   0.977   0.999    FALSE 1 1.002   880
#> eta[10,6]   0.880  0.055   0.761   0.887   0.967    FALSE 1 0.999  1500
#> eta[5,7]    0.566  0.080   0.413   0.564   0.718    FALSE 1 1.004   474
#> eta[8,7]    0.972  0.027   0.897   0.980   1.000    FALSE 1 1.003  1500
#> eta[10,7]   0.361  0.078   0.218   0.358   0.524    FALSE 1 1.007   259
#> eta[4,8]    0.673  0.092   0.478   0.681   0.825    FALSE 1 1.002   712
#> eta[5,8]    0.585  0.096   0.397   0.585   0.773    FALSE 1 0.999  1500
#> eta[9,8]    0.225  0.082   0.084   0.220   0.396    FALSE 1 1.009   218
#> eta[10,8]   0.315  0.094   0.147   0.309   0.510    FALSE 1 1.022   124
#> eta[4,9]    0.955  0.043   0.839   0.966   0.999    FALSE 1 1.004  1500
#> eta[5,10]   0.961  0.038   0.856   0.972   0.999    FALSE 1 1.003  1500
#> PI[4,1]     0.297  0.376   0.000   0.056   1.000    FALSE 1 1.213    16
#> PI[5,1]     0.703  0.376   0.000   0.944   1.000    FALSE 1 1.213    16
#> PI[3,2]     0.094  0.133   0.000   0.030   0.469    FALSE 1 1.075    36
#> PI[6,2]     0.210  0.262   0.000   0.023   0.794    FALSE 1 1.581     7
#> PI[7,2]     0.503  0.310   0.000   0.521   0.984    FALSE 1 1.590     7
#> PI[8,2]     0.125  0.179   0.000   0.053   0.702    FALSE 1 1.345    12
#> PI[9,2]     0.069  0.098   0.000   0.025   0.324    FALSE 1 1.008  1500
#> PI[1,3]     0.242  0.246   0.000   0.215   0.798    FALSE 1 1.938     5
#> PI[4,3]     0.114  0.128   0.000   0.080   0.392    FALSE 1 1.030   243
#> PI[8,3]     0.020  0.052   0.000   0.001   0.166    FALSE 1 1.057   302
#> PI[9,3]     0.624  0.273   0.000   0.636   0.999    FALSE 1 1.671     6
#> PI[1,6]     0.165  0.207   0.000   0.083   0.748    FALSE 1 1.037    59
#> PI[3,6]     0.376  0.248   0.001   0.367   0.891    FALSE 1 1.134    21
#> PI[8,6]     0.153  0.239   0.000   0.014   0.801    FALSE 1 1.188    22
#> PI[9,6]     0.204  0.226   0.000   0.108   0.721    FALSE 1 1.025   104
#> PI[10,6]    0.102  0.149   0.000   0.023   0.506    FALSE 1 1.011   349
#> PI[5,7]     0.069  0.094   0.000   0.026   0.309    FALSE 1 1.105    28
#> PI[8,7]     0.251  0.375   0.000   0.011   1.000    FALSE 1 2.669     4
#> PI[10,7]    0.680  0.367   0.000   0.862   1.000    FALSE 1 2.315     5
#> PI[4,8]     0.241  0.298   0.000   0.130   1.000    FALSE 1 1.612     7
#> PI[5,8]     0.085  0.157   0.000   0.000   0.523    FALSE 1 1.619     7
#> PI[9,8]     0.198  0.282   0.000   0.035   0.920    FALSE 1 1.729     6
#> PI[10,8]    0.475  0.380   0.000   0.551   1.000    FALSE 1 2.843     4
#> 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.222 11.798 848.277 867.425 893.518    FALSE 1 1.082    30
#> 
#> **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 = 64.9 and DIC = 933.148 
#> 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)