library(renv)
renv::load()
renv::restore() # Install missing librairies.
renv::status() # Check the project state.
library(targets)
library(here)
tar_config_set(
store = here::here("_targets")
)
library(foodwebbuilder)
library(dplyr)
library(tibble)
library(riversea)
library(ggplot2)
set_theme(theme_minimal())
library(INLA)Build food webs
Build the metaweb
To build the metaweb, we need the following.
The diet table for fish species considered.
diet_fish <- tar_read(diet_size)
head(diet_fish)# A tibble: 6 × 15
species length_min length_max stage biofilm crustacean detritus echinoderm
<chr> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
1 Abramis br… 0 2.00 larv… 0 0 0 0
2 Abramis br… 2.00 2.00 juv.… 0 1 1 0
3 Abramis br… 2.00 Inf adul… 0 1 1 0
4 Agonus cat… 0 2.00 larv… 0 0 0 0
5 Agonus cat… 2.00 Inf comb… 0 1 0 0
6 Alburnus a… 0 2.00 larv… 0 0 0 0
# ℹ 7 more variables: fish <dbl>, insect <dbl>, macrophyte <dbl>,
# mollusk <dbl>, phytoplankton <dbl>, worm <dbl>, zooplankton <dbl>The diet for the resources considered.
diet_resource <- tar_read(diet_resource)The size of each individual fish.
size <- tar_read(sea_data_imputed)$size
head(size) trait length year survey species_valid measured
1 61836630 6.3 2016 pomet Abramis brama TRUE
2 61836630 5.5 2016 pomet Abramis brama TRUE
3 61836630 7.0 2016 pomet Abramis brama TRUE
4 61836630 6.5 2016 pomet Abramis brama TRUE
5 61836630 19.3 2016 pomet Abramis brama TRUE
6 61836630 6.1 2016 pomet Abramis brama TRUEThe predation window.
predation_window <- tar_read(predation_window)
head(predation_window)# A tibble: 6 × 3
species beta_min beta_max
<chr> <dbl> <dbl>
1 Abramis brama 0.03 0.45
2 Alosa alosa 0.03 0.45
3 Alosa fallax 0.03 0.45
4 Ameiurus melas 0.03 0.45
5 Anguilla anguilla 0.03 0.45
6 Argyrosomus regius 0.03 0.45From this we can build the metaweb using helper functions of the foodwebbuilder package. We wrapped them in a riversea function so the metaweb can be build directly. We just need to specify the number of size class, and the selected resources, although default values are set if not specified.
For information, here is the structure that we consider between the different trophic groups on which fish can feed on.
tar_read(plot_net_groups)
We have already computed metaweb for different number of size classes, so we can check whether this parameter influences the metaweb structure.
metaweb_table <- tar_read(metaweb_table)
metaweb_table# A tibble: 7 × 2
num_classes metaweb
<int> <list>
1 3 <dbl [319 × 319]>
2 4 <dbl [422 × 422]>
3 5 <dbl [525 × 525]>
4 6 <dbl [628 × 628]>
5 7 <dbl [731 × 731]>
6 8 <dbl [834 × 834]>
7 9 <dbl [937 × 937]>We can plot the metaweb connectance against the number of size classes.
tar_read(plot_metaweb_connectance)
We observe little variation indicating that the number of size class does not change too much the structure of the food web reconstructed.
Get local food webs
web_list <- tar_read(web_list)
get_connectance(web_list$metaweb)[1] 0.1714286local_foodwebs <- enframe(web_list$local, name = "trait", value = "foodweb")
head(local_foodwebs)# A tibble: 6 × 2
trait foodweb
<chr> <list>
1 61836630 <dbl [21 × 21]>
2 61832572 <dbl [22 × 22]>
3 61835462 <dbl [21 × 21]>
4 61832074 <dbl [19 × 19]>
5 61836249 <dbl [20 × 20]>
6 61837171 <dbl [12 × 12]>local_foodwebs <- local_foodwebs |>
mutate(
log_richness = log(purrr::map_dbl(foodweb, get_richness)),
connectance = purrr::map_dbl(foodweb, get_connectance)
)
head(local_foodwebs)# A tibble: 6 × 4
trait foodweb log_richness connectance
<chr> <list> <dbl> <dbl>
1 61836630 <dbl [21 × 21]> 3.04 0.170
2 61832572 <dbl [22 × 22]> 3.09 0.209
3 61835462 <dbl [21 × 21]> 3.04 0.143
4 61832074 <dbl [19 × 19]> 2.94 0.158
5 61836249 <dbl [20 × 20]> 3.00 0.232
6 61837171 <dbl [12 × 12]> 2.48 0.25 trait_data <- tar_read(sea_data_tidy)$trait
local_foodwebs <- local_foodwebs |>
left_join(trait_data)
head(local_foodwebs)# A tibble: 6 × 9
trait foodweb log_richness connectance year month longitude latitude survey
<chr> <list> <dbl> <dbl> <int> <int> <dbl> <dbl> <chr>
1 61836… <dbl[…]> 3.04 0.170 2016 10 -0.164 49.2 pomet
2 61832… <dbl[…]> 3.09 0.209 2015 6 1.04 49.3 pomet
3 61835… <dbl[…]> 3.04 0.143 2006 7 -1.27 47.3 pomet
4 61832… <dbl[…]> 2.94 0.158 2011 9 -0.259 49.2 pomet
5 61836… <dbl[…]> 3.00 0.232 2007 5 -1.06 46.3 pomet
6 61837… <dbl[…]> 2.48 0.25 2006 11 -1.37 47.3 pomet We have computed the connectance and richness of local metaweb so we can plot the relationship between these two variables.
m1 <- inla(
connectance ~ 1 + log_richness + f(survey, model = "iid") + f(year, model = "iid"),
data = local_foodwebs |> select(-foodweb)
)
summary(m1)Time used:
Pre = 0.475, Running = 0.581, Post = 0.464, Total = 1.52
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.441 0.006 0.429 0.441 0.453 0.441 0
log_richness -0.096 0.001 -0.098 -0.096 -0.094 -0.096 0
Random effects:
Name Model
survey IID model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for the Gaussian observations 2238.61 31.84 2176.49 2238.41
Precision for survey 26060.92 19666.06 2539.33 20948.68
Precision for year 26741.54 6890.06 15526.11 25959.26
0.975quant mode
Precision for the Gaussian observations 2301.86 2238.08
Precision for survey 74797.11 7697.46
Precision for year 42458.72 24525.24
Marginal log-Likelihood: 24114.40
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')intercept <- m1$summary.fixed["(Intercept)", "mean"]
slope <- m1$summary.fixed["log_richness", "mean"]
ggplot(local_foodwebs, aes(x = log_richness, y = connectance, color = survey)) +
geom_point(alpha = 0.3) +
geom_abline(intercept = intercept, slope = slope) +
labs(x = "Log richness", y = "Connectance", color = "Survey")
As expected, we observe a negative trend between connectance and species richness. There is also a small random effect attributed to the survey.
Next, we investigate the effect of the year on species richness.
ggplot(local_foodwebs, aes(x = year, y = log_richness, color = survey)) +
geom_point(alpha = 0.1) +
stat_summary(
fun = mean,
geom = "line",
linewidth = 1
) +
labs(x = "Year", y = "Log richness", color = "Survey")
As we observe a positive trend between richness and year, and that we have observed a negative trend between richness and connectance, we expect a negative trend between connectance and year.
ggplot(local_foodwebs, aes(x = year, y = connectance, color = survey)) +
geom_point(alpha = 0.1) +
stat_summary(
fun = mean,
geom = "line",
linewidth = 1
) +
labs(x = "Year", y = "Connectance", color = "Survey")
Interestingly, we see that connectance values are in agreement between the two survey (when there is an overlap), suggesting that it makes sense to merge these two different data sources.
TODO
- Use life stage lengths from Anik dataset to refine diet table;
- Split zoobenthos in three categories: molluscs, crustaceans and polychaetes;
- Relate food web structure to environmental variables.