---
title: "flower_plots"
output: html_document
date: "2024-02-26"
---
#===============================================================================================#
# GENERAL INFORMATION
This R script contains all statistical analyses presented in the manuscript
'Plant identity and traits determine pollinator, natural enemy, herbivore
and decomposer abundances in flower plantings'
by Neus Rodriguez-Gasol, Fabian A. Boetzl, Elodie Chapurlat, Johan A. Stenberg, Mattias Jonsson,
Ola Lundin and Maria Viketoft. For a detailed description of sampling protocols and inclusion
criteria, please see the main text.
All analyses were performed in R version 4.4.1 (2024-06-14 ucrt) on x86_64, mingw32 system.
#===============================================================================================#
# Set working directory
```{r setup, include=FALSE}
# remove junk first
rm(list=ls())
# library("knitr")
## setting working directory (here you need to change to the folder where the data are)
knitr::opts_knit$set(root.dir = ".")
setwd('.') # where '.' is your working directory
getwd()
```
#===============================================================================================#
# Loading required packages and defining plotting aesthetics for figures
```{r, echo = FALSE, message = FALSE, warning = FALSE, results = FALSE}
# packages & versions used
library(dplyr) # 1.1.4
library(tidyr) # 1.3.1
library(ggplot2) # 3.5.1
library(ggpubr) # 0.6.0
library(viridis) # 0.6.5
library(cowplot) # 1.1.3
library(glmmTMB) # 1.1.10
library(DHARMa) # 0.4.7
library(performance) # 0.12.4
library(car) # 3.1-3
library(emmeans) # 1.10.5
library(ggeffects) # 1.7.2
library(MASS) # 7.3-60.2
library(parameters) # 0.23.0
# plotting aesthetics
plot_theme <- theme(axis.title = element_text(color="black", face='bold', size=14),
axis.text.x = element_text(color="black", size=12),
axis.text.y = element_text(color="black", size=12),
panel.border = element_rect(colour = "black", fill=NA, size=1),
panel.background = element_rect(colour = NA, fill = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
legend.title = element_blank(),
legend.key = element_rect(fill = NA),
legend.text = element_text(size=18),
legend.background=element_blank(),
plot.title = element_text(size=20, face="bold",hjust = 0.5))
# car Anova computes wrong type III tests according to Sarraj & Forkman (2023) - change option for correct one
options(contrasts = c('contr.sum', 'contr.poly'))
```
#===============================================================================================#
# Load datasets & data handling
```{r, echo = FALSE, message = FALSE, warning = FALSE, results = FALSE}
# plant list for merging
plant_list <- read.csv("plant_list.csv", header=T, sep=",", na.strings="NA", fileEncoding = "ISO-8859-1")
plant_list$plot <- as.factor(plant_list$plot)
plant_list$plant_name <- as.factor(plant_list$plant_name)
# soil fauna
soilfauna <- read.csv("soil_fauna.csv", header=T, sep=",", na.strings="NA", fileEncoding = "ISO-8859-1")
soilfauna$year <- as.factor(soilfauna$year)
soilfauna$site <- as.factor(soilfauna$site)
soilfauna$plot <- as.factor(soilfauna$plot)
soilfauna$block <- as.factor(soilfauna$block)
soilfauna$exclude <- as.factor(soilfauna$exclude)
# leaf samples
leafsuction <- read.csv("leaf_suction.csv", header=T, sep=",", na.strings="NA", fileEncoding = "ISO-8859-1")
leafsuction$year <- as.factor(leafsuction$year)
leafsuction$site <- as.factor(leafsuction$site)
leafsuction$plot <- as.factor(leafsuction$plot)
leafsuction$block <- as.factor(leafsuction$block)
leafsuction$exclude <- as.factor(leafsuction$exclude)
# pollinators
pollinators <- read.csv("pollinators.csv", header=T, sep=",", na.strings="NA", fileEncoding = "ISO-8859-1")
pollinators$year <- as.factor(pollinators$year)
pollinators$site <- as.factor(pollinators$site)
pollinators$plot <- as.factor(pollinators$plot)
pollinators$block <- as.factor(pollinators$block)
pollinators$exclude <- as.factor(pollinators$exclude)
# peak bloom (and flower, species and weed cover at peak bloom)
peak_bloom <- read.csv("peakbloom.csv", header=T, sep=",", na.strings="NA", fileEncoding = "ISO-8859-1")
peak_bloom$year <- as.factor(peak_bloom$year)
peak_bloom$site <- as.factor(peak_bloom$site)
peak_bloom$plot <- as.factor(peak_bloom$plot)
peak_bloom$block <- as.factor(peak_bloom$block)
peak_bloom$exclude <- as.factor(peak_bloom$exclude)
# NOTE: 'Exclude' only refers to flower cover at peak bloom here!
peak_bloom$site_year <- interaction(peak_bloom$site, peak_bloom$year, sep='_')
# flower cover all season
flower_cover <- read.csv("flower_cover_season.csv", header=T, sep=",", na.strings="NA", fileEncoding = "ISO-8859-1")
flower_cover$year <- as.factor(flower_cover$year)
flower_cover$site <- as.factor(flower_cover$site)
flower_cover$plot <- as.factor(flower_cover$plot)
flower_cover$block <- as.factor(flower_cover$block)
flower_cover$exclude <- as.factor(flower_cover$exclude)
flower_cover$site_year <- interaction(flower_cover$site, flower_cover$year, sep='_')
# groundcover
groundcover <- read.csv("groundcover.csv", header=T, sep=",", na.strings="NA", fileEncoding = "ISO-8859-1")
groundcover$year <- as.factor(groundcover$year)
groundcover$plot <- as.factor(groundcover$plot)
groundcover$block <- as.factor(groundcover$block)
groundcover$exclude <- as.factor(groundcover$exclude)
groundcover$site_year <- interaction(groundcover$site, groundcover$year, sep='_')
```
#===============================================================================================#
# exclude samples (based on inclusion criteria, see main text)
```{r}
soilfauna <- droplevels(subset(soilfauna, exclude == FALSE))
soilfauna$site_year <- interaction(soilfauna$site, soilfauna$year, sep='_')
soilfauna$lg_soil_weight_g <- log(soilfauna$soil_weight_g)
leafsuction <- droplevels(subset(leafsuction, exclude == FALSE))
leafsuction$site_year <- interaction(leafsuction$site, leafsuction$year, sep='_')
pollinators <- droplevels(subset(pollinators, exclude == FALSE))
pollinators$site_year <- interaction(pollinators$site, pollinators$year, sep='_')
```
#===============================================================================================#
# Modelling part
## (a) hoverflies
```{r}
# restrict dataset
hoverflies <- pollinators %>% dplyr::select(year, site, site_year, plot, block, hoverflies, observations_is, observations_should)
hoverflies %>% group_by(plot) %>% tally()
# merge with plant list for names
hoverflies <- merge(hoverflies, plant_list, by = 'plot', all.x = T)
# calculate sample completeness
hoverflies$completeness <- hoverflies$observations_is / hoverflies$observations_should
# model
M_hoverflies <- glmmTMB(hoverflies ~ plant_name +
site +
year +
offset(log(completeness)) +
(1|site_year/block),
family = nbinom2(),
dispformula = ~ year,
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")),
data = hoverflies)
hist(residuals(M_hoverflies), breaks=100)
plot(simulateResiduals(M_hoverflies, n=1000))
testDispersion(M_hoverflies)
check_collinearity(M_hoverflies)
check_zeroinflation(M_hoverflies)
Anova(M_hoverflies, type=3)
# Chisq Df Pr(>Chisq)
# (Intercept) 37.7311 1 8.12e-10 ***
# plant_name 329.2606 26 < 2.2e-16 ***
# site 0.4702 1 0.4929
# year 0.3142 1 0.5751
# calculate residual df (based on https://www.rdocumentation.org/packages/glmmTMB/versions/1.1.9/topics/glmmTMB)
df.residual(M_hoverflies) * (sqrt(nobs(M_hoverflies) / (1+df.residual(M_hoverflies))))
# 288.0297
# summary(M_hoverflies)
performance(M_hoverflies)
# AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | RMSE | Sigma | Score_log | Score_spherical
# -------------------------------------------------------------------------------------------------------
# 1996.918 | 2004.683 | 2115.968 | | 0.720 | 16.357 | 1.683 | -3.179 | 0.031
# extracting marginal means - response scale
means_hoverflies2 <- as.data.frame(emmeans(M_hoverflies, specs = 'plant_name', type = "response"))
means_hoverflies2$group <- 'hoverflies'
means_hoverflies2 <- means_hoverflies2 %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
# calculate values relative to maximum (for radar plots)
means_hoverflies2$prop <- means_hoverflies2$response / max(means_hoverflies2$response)
```
## (b) wild bees
```{r}
# restrict dataset
wild_bees <- pollinators %>% dplyr::select(year, site, site_year, plot, block, wild_bees, observations_is, observations_should)
wild_bees %>% group_by(plot) %>% tally()
# merge with plant list for names
wild_bees <- merge(wild_bees, plant_list, by = 'plot', all.x = T)
# calculate sample completeness
wild_bees$completeness <- wild_bees$observations_is / wild_bees$observations_should
# model
M_wild_bees <- glmmTMB(wild_bees ~ plant_name +
site +
year +
offset(log(completeness)) +
(1|site_year/block),
family = nbinom2(),
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")),
data = wild_bees)
hist(residuals(M_wild_bees), breaks=100)
plot(simulateResiduals(M_wild_bees, n=1000))
testDispersion(M_wild_bees)
check_collinearity(M_wild_bees)
Anova(M_wild_bees, type=3)
# Chisq Df Pr(>Chisq)
# (Intercept) 237.9199 1 < 2e-16 ***
# plant_name 535.9687 26 < 2e-16 ***
# site 0.1505 1 0.69804
# year 6.5672 1 0.01039 *
# calculate residual df (based on https://www.rdocumentation.org/packages/glmmTMB/versions/1.1.9/topics/glmmTMB)
df.residual(M_wild_bees) * (sqrt(nobs(M_wild_bees) / (1+df.residual(M_wild_bees))))
# 288.0297
# summary(M_wild_bees)
performance(M_wild_bees)
# AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | RMSE | Sigma | Score_log | Score_spherical
# -------------------------------------------------------------------------------------------------------
# 1960.323 | 1968.088 | 2079.373 | | 0.833 | 32.852 | 1.316 | -3.135 | 0.031
# extracting marginal means - response scale
means_wild_bees2 <- as.data.frame(emmeans(M_wild_bees, specs = 'plant_name', type = "response"))
means_wild_bees2$group <- 'wild bees'
means_wild_bees2 <- means_wild_bees2 %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
# calculate values relative to maximum (for radar plots)
means_wild_bees2$prop <- means_wild_bees2$response / max(means_wild_bees2$response)
```
## (c) herbivorous arthropods
```{r}
# restrict dataset
leaf_herbivores <- leafsuction %>% dplyr::select(year, site, site_year, plot, block, herbivores, observations)
leaf_herbivores %>% group_by(plot) %>% tally()
# merge with plant list for names
leaf_herbivores <- merge(leaf_herbivores, plant_list, by = 'plot', all.x = T)
# model
M_leaf_herbivores <- glmmTMB(herbivores ~ plant_name +
site +
year +
offset(log(observations)) +
(1|site_year/block),
family = nbinom2(),
dispformula = ~site_year,
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")),
data = leaf_herbivores)
hist(residuals(M_leaf_herbivores), breaks=100)
plot(simulateResiduals(M_leaf_herbivores, n=1000))
testDispersion(M_leaf_herbivores)
check_collinearity(M_leaf_herbivores)
Anova(M_leaf_herbivores, type=3)
# Chisq Df Pr(>Chisq)
# (Intercept) 484.2957 1 < 2e-16 ***
# plant_name 274.7628 26 < 2e-16 ***
# site 0.0058 1 0.93926
# year 5.7669 1 0.01633 *
# calculate residual df (based on https://www.rdocumentation.org/packages/glmmTMB/versions/1.1.9/topics/glmmTMB)
df.residual(M_leaf_herbivores) * (sqrt(nobs(M_leaf_herbivores) / (1+df.residual(M_leaf_herbivores))))
# 282.0197
# summary(M_leaf_herbivores)
performance(M_leaf_herbivores)
# AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | Score_log | Score_spherical
# ----------------------------------------------------------------------------------------------------------------
# 3165.504 | 3173.444 | 3283.918 | 0.690 | 0.547 | 0.316 | 137.521 | 2.971 | -5.126 | 0.041
# extracting marginal means - response scale
means_leaf_herbivores2 <- as.data.frame(emmeans(M_leaf_herbivores, specs = 'plant_name', type = "response"))
means_leaf_herbivores2$group <- 'leaf herbivores'
means_leaf_herbivores2 <- means_leaf_herbivores2 %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
# calculate values relative to maximum (for radar plots)
means_leaf_herbivores2$prop <- means_leaf_herbivores2$response / max(means_leaf_herbivores2$response)
```
## (d) parasitic wasps
```{r}
# restrict dataset
leaf_parasitoids <- leafsuction %>% dplyr::select(year, site, site_year, plot, block, parasitoids, observations)
leaf_parasitoids %>% group_by(plot) %>% tally()
# merge with plant list for names
leaf_parasitoids <- merge(leaf_parasitoids, plant_list, by = 'plot', all.x = T)
# model
M_leaf_parasitoids <- glmmTMB(parasitoids ~ plant_name +
site +
year +
offset(log(observations)) +
(1|site_year/block),
family = nbinom2(),
data = leaf_parasitoids)
hist(residuals(M_leaf_parasitoids), breaks=100)
plot(simulateResiduals(M_leaf_parasitoids, n=1000))
testDispersion(M_leaf_parasitoids)
check_collinearity(M_leaf_parasitoids)
Anova(M_leaf_parasitoids, type=3)
# Chisq Df Pr(>Chisq)
# (Intercept) 154.6198 1 < 2e-16 ***
# plant_name 251.0564 26 < 2e-16 ***
# site 0.0092 1 0.92358
# year 3.7230 1 0.05367 .
# calculate residual df (based on https://www.rdocumentation.org/packages/glmmTMB/versions/1.1.9/topics/glmmTMB)
df.residual(M_leaf_parasitoids) * (sqrt(nobs(M_leaf_parasitoids) / (1+df.residual(M_leaf_parasitoids))))
# 282.0197
# summary(M_leaf_parasitoids)
performance(M_leaf_parasitoids)
# AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | Score_log | Score_spherical
# ---------------------------------------------------------------------------------------------------------------
# 2562.015 | 2569.955 | 2680.429 | 0.660 | 0.471 | 0.356 | 44.073 | 2.232 | -4.136 | 0.040
# extracting marginal means - response scale
means_leaf_parasitoids2 <- as.data.frame(emmeans(M_leaf_parasitoids, specs = 'plant_name', type = "response"))
means_leaf_parasitoids2$group <- 'leaf parasitoids'
means_leaf_parasitoids2 <- means_leaf_parasitoids2 %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
# calculate values relative to maximum (for radar plots)
means_leaf_parasitoids2$prop <- means_leaf_parasitoids2$response / max(means_leaf_parasitoids2$response)
```
## (e) predatory arthropods
```{r}
# restrict dataset
leaf_predators <- leafsuction %>% dplyr::select(year, site, site_year, plot, block, predators, observations) %>% droplevels()
leaf_predators %>% group_by(plot) %>% tally()
# merge with plant list for names
leaf_predators <- merge(leaf_predators, plant_list, by = 'plot', all.x = T)
# model
M_leaf_predators <- glmmTMB(predators ~ plant_name +
site +
year +
offset(log(observations)) +
(1|site_year/block),
family = nbinom2(),
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")),
data = leaf_predators)
hist(residuals(M_leaf_predators), breaks=100)
plot(simulateResiduals(M_leaf_predators, n=1000))
testDispersion(M_leaf_predators)
check_collinearity(M_leaf_predators)
Anova(M_leaf_predators, type=3)
# Chisq Df Pr(>Chisq)
# (Intercept) 990.6234 1 < 2.2e-16 ***
# plant_name 112.9236 26 8.336e-13 ***
# site 0.1013 1 0.7503
# year 28.3817 1 9.960e-08 ***
# calculate residual df (based on https://www.rdocumentation.org/packages/glmmTMB/versions/1.1.9/topics/glmmTMB)
df.residual(M_leaf_predators) * (sqrt(nobs(M_leaf_predators) / (1+df.residual(M_leaf_predators))))
# 282.0197
# summary(M_leaf_predators)
performance(M_leaf_predators)
# AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | RMSE | Sigma | Score_log | Score_spherical
# -------------------------------------------------------------------------------------------------------
# 2502.067 | 2510.007 | 2620.481 | | 0.459 | 20.596 | 2.230 | -4.093 | 0.045
# extracting marginal means - response scale
means_leaf_predators2 <- as.data.frame(emmeans(M_leaf_predators, specs = 'plant_name', type = "response"))
means_leaf_predators2$group <- 'leaf predators'
means_leaf_predators2 <- means_leaf_predators2 %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
# calculate values relative to maximum (for radar plots)
means_leaf_predators2$prop <- means_leaf_predators2$response / max(means_leaf_predators2$response)
```
## (f) herbivorous nematodes
```{r}
# restrict dataset
soil_plant_feeders <- soilfauna %>% dplyr::select(year, site, site_year, plot, block, plant_feeders, soil_weight_g, lg_soil_weight_g)
soil_plant_feeders %>% group_by(plot) %>% tally()
# merge with plant list for names
soil_plant_feeders <- merge(soil_plant_feeders, plant_list, by = 'plot', all.x = T)
# model
M_soil_plant_feeders <- glmmTMB(plant_feeders ~ plant_name +
site +
year +
offset(log(soil_weight_g)) +
(1|site_year/block),
family = nbinom2(),
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")),
data = soil_plant_feeders)
hist(residuals(M_soil_plant_feeders), breaks=100)
plot(simulateResiduals(M_soil_plant_feeders, n=1000))
testDispersion(M_soil_plant_feeders)
check_collinearity(M_soil_plant_feeders)
Anova(M_soil_plant_feeders, type=3)
# Chisq Df Pr(>Chisq)
# (Intercept) 0.2791 1 0.5973
# plant_name 122.3135 26 1.949e-14 ***
# site 21.6866 1 3.210e-06 ***
# year 16.5021 1 4.860e-05 ***
# calculate residual df (based on https://www.rdocumentation.org/packages/glmmTMB/versions/1.1.9/topics/glmmTMB)
df.residual(M_soil_plant_feeders) * (sqrt(nobs(M_soil_plant_feeders) / (1+df.residual(M_soil_plant_feeders))))
# 292.0361
summary(M_soil_plant_feeders)
performance(M_soil_plant_feeders)
# AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | RMSE | Sigma | Score_log | Score_spherical
# -------------------------------------------------------------------------------------------------------
# 2346.542 | 2354.194 | 2466.009 | | 0.390 | 42.427 | 1.034 | -Inf | 0.022
# extracting marginal means - response scale
means_soil_plant_feeders2 <- as.data.frame(emmeans(M_soil_plant_feeders, specs = 'plant_name', type = "response"))
means_soil_plant_feeders2$group <- 'soil plant feeders'
means_soil_plant_feeders2 <- means_soil_plant_feeders2 %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
# calculate values relative to maximum (for radar plots)
means_soil_plant_feeders2$prop <- means_soil_plant_feeders2$response / max(means_soil_plant_feeders2$response)
```
## (g) predatory nematodes
```{r}
# restrict dataset
soil_predators <- soilfauna %>% dplyr::select(year, site, site_year, plot, block, predators, soil_weight_g, lg_soil_weight_g)
soil_predators %>% group_by(plot) %>% tally()
# merge with plant list for names
soil_predators <- merge(soil_predators, plant_list, by = 'plot', all.x = T)
# model
M_soil_predators <- glmmTMB(predators ~ plant_name +
site + year +
offset(log(soil_weight_g)) +
(1|site_year/block),
zi =~1,
family = nbinom2(),
dispformula = ~ site_year,
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")),
data = soil_predators)
hist(residuals(M_soil_predators), breaks=100)
plot(simulateResiduals(M_soil_predators, n=1000))
testDispersion(M_soil_predators)
check_collinearity(M_soil_predators)
Anova(M_soil_predators, type=3)
# Chisq Df Pr(>Chisq)
# (Intercept) 22.8607 1 1.742e-06 ***
# plant_name 32.0298 26 0.1921375
# site 13.8618 1 0.0001968 ***
# year 3.1616 1 0.0753885 .
# calculate residual df (based on https://www.rdocumentation.org/packages/glmmTMB/versions/1.1.9/topics/glmmTMB)
df.residual(M_soil_predators) * (sqrt(nobs(M_soil_predators) / (1+df.residual(M_soil_predators))))
# 292.0361
# summary(M_soil_predators)
performance(M_soil_predators)
# AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | Score_log | Score_spherical
# --------------------------------------------------------------------------------------------------------------
# 1606.050 | 1613.703 | 1725.517 | 0.481 | 0.354 | 0.196 | 5.616 | 1.881 | -Inf | 0.030
# extracting marginal means - response scale
means_soil_predators2 <- as.data.frame(emmeans(M_soil_predators, specs = 'plant_name', type = "response"))
means_soil_predators2$group <- 'soil predators'
means_soil_predators2 <- means_soil_predators2 %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
# calculate values relative to maximum (for radar plots)
means_soil_predators2$prop <- means_soil_predators2$response / max(means_soil_predators2$response)
```
## (h) decomposer nematodes
```{r}
# restrict dataset
soil_decomposers <- soilfauna %>% dplyr::select(year, site, site_year, plot, block, decomposers, soil_weight_g, lg_soil_weight_g)
soil_decomposers %>% group_by(plot) %>% tally()
# merge with plant list for names
soil_decomposers <- merge(soil_decomposers, plant_list, by = 'plot', all.x = T)
# model
M_soil_decomposers <- glmmTMB(decomposers ~ plant_name +
site +
year +
offset(log(soil_weight_g)) +
(1|site_year/block),
family = nbinom2(link = 'log'),
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")),
data = soil_decomposers)
hist(residuals(M_soil_decomposers), breaks=100)
plot(simulateResiduals(M_soil_decomposers, n=1000))
testDispersion(M_soil_decomposers)
check_collinearity(M_soil_decomposers)
Anova(M_soil_decomposers, type=3)
# Chisq Df Pr(>Chisq)
# (Intercept) 1258.0443 1 < 2.2e-16 ***
# plant_name 70.2547 26 6.091e-06 ***
# site 0.6085 1 0.435355
# year 8.6732 1 0.003229 **
# calculate residual df (based on https://www.rdocumentation.org/packages/glmmTMB/versions/1.1.9/topics/glmmTMB)
df.residual(M_soil_decomposers) * (sqrt(nobs(M_soil_decomposers) / (1+df.residual(M_soil_decomposers))))
# 292.0361
# summary(M_soil_decomposers)
performance(M_soil_decomposers)
# AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | RMSE | Sigma | Score_log | Score_spherical
# --------------------------------------------------------------------------------------------------------
# 3772.025 | 3779.677 | 3891.492 | | 0.229 | 145.501 | 2.198 | -Inf | 0.025
# extracting marginal means - response scale
means_soil_decomposers2 <- as.data.frame(emmeans(M_soil_decomposers, specs = 'plant_name', type = "response"))
means_soil_decomposers2$group <- 'soil decomposers'
means_soil_decomposers2 <- means_soil_decomposers2 %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
# calculate values relative to maximum (for radar plots)
means_soil_decomposers2$prop <- means_soil_decomposers2$response / max(means_soil_decomposers2$response)
```
#===============================================================================================#
# Figure 1
```{r}
# bing model predictions together in one table
means <- means_soil_plant_feeders2
means <- rbind(means, means_soil_decomposers2)
means <- rbind(means, means_soil_predators2)
means <- rbind(means, means_leaf_predators2)
means <- rbind(means, means_leaf_parasitoids2)
means <- rbind(means, means_leaf_herbivores2)
means <- rbind(means, means_wild_bees2)
means <- rbind(means, means_hoverflies2)
# calculate percentages out of proportions
means$prop100 <- (means$prop * 100)
# reorder factor
means$group <- ordered(means$group, levels =c('hoverflies',
'wild bees',
'leaf herbivores',
'leaf parasitoids',
'leaf predators',
'soil plant feeders',
'soil predators',
'soil decomposers'))
# include dummy for legend
means_zzz <- data.frame(group = c('hoverflies', 'wild bees', 'leaf herbivores', 'leaf parasitoids',
'leaf predators', 'soil plant feeders', 'soil predators', 'soil decomposers'),
# insert additional 'dummy plant' with prop = 100 for legend
plant_name = c('zzz', 'zzz', 'zzz', 'zzz', 'zzz', 'zzz', 'zzz', 'zzz'),
response = c(NA, NA, NA, NA, NA, NA, NA, NA),
asymp.LCL = c(NA, NA, NA, NA, NA, NA, NA, NA),
asymp.UCL = c(NA, NA, NA, NA, NA, NA, NA, NA),
prop = c(NA, NA, NA, NA, NA, NA, NA, NA),
prop100 = 100)
means <- rbind(means, means_zzz)
# define plot aesthetics
plot_theme_radar <- theme(axis.title = element_blank(),
axis.text.x = element_text(color="black", size=12),
axis.text.y = element_text(color="black", size=12),
panel.border = element_blank(),
panel.background = element_rect(colour = NA, fill = NA),
panel.grid.major = element_blank(),
legend.position = "bottom",
legend.title = element_blank(),
legend.text = element_text(size=12),
legend.background=element_blank(),
plot.title = element_text(size=20, face="bold",hjust = 0.5))
# Define radarplot lines
segments <- c(0,25,50,75,100)
# actual plot
radarplot <- ggplot(data=means, aes(x = group, y = prop100, fill = group)) +
coord_polar() +
geom_hline(yintercept = segments, size=0.35, color = "grey80") +
geom_col(aes(x = group, y = prop100, fill = group), width = 0.85)+
scale_fill_viridis(discrete = TRUE,
labels=c('hoverflies',
'wild bees',
'herbivorous arthropods',
'parasitic wasps',
'predatory arthropods',
'herbivorous nematodes',
'predatory nematodes',
'decomposer nematodes')) +
scale_y_continuous('',breaks = seq(0, 100, by = 25), limits=c(-10,100), expand =c(0,0))+
facet_wrap(plant_name ~ ., ncol = 7) +
plot_theme_radar +
theme(axis.text.y=element_text(size=9),
axis.text.x = element_blank(),
strip.text = element_text(size = 10, face = 'italic'),
strip.text.y.left = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
strip.background = element_rect(colour=NA, fill="white"))
# export plot
tiff("Figure_1.tiff", res=800, width=340, height=240, units="mm", compression = "lzw")
radarplot
dev.off()
```
#===============================================================================================#
# Modelling traits
## (a) peak bloom
```{r}
# restrict dataset
peakbloom <- peak_bloom %>% dplyr::select(year, site, site_year, plot, block, peak_avg)
peakbloom %>% group_by(plot) %>% tally()
# merge with plant list for names
peakbloom <- merge(peakbloom, plant_list, by = 'plot', all.x = T)
# model
M_peakbloom <- glmmTMB(peak_avg ~ plant_name +
site +
year +
(1|site_year/block),
family = gaussian(),
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")),
data = peakbloom)
hist(residuals(M_peakbloom), breaks=100)
plot(simulateResiduals(M_peakbloom, n=1000))
testDispersion(M_peakbloom)
plotResiduals(simulateResiduals(M_peakbloom, n=1000), form=peakbloom$site)
check_collinearity(M_peakbloom)
plotResiduals(simulateResiduals(M_peakbloom, n=1000), form=peakbloom$site_year)
Anova(M_peakbloom, type=3)
# Chisq Df Pr(>Chisq)
# plant_name 106.6162 26 9.946e-12 ***
# site 0.2598 1 0.6103
# year 81.4294 1 < 2.2e-16 ***
performance(M_peakbloom)
# extracting marginal means - response scale
means_peakbloom <- as.data.frame(emmeans(M_peakbloom, specs = 'plant_name', type = "response"))
means_peakbloom$group <- 'peak_bloom'
means_peakbloom <- means_peakbloom %>% dplyr::select(group, plant_name, emmean, lower.CL, upper.CL)
```
## (b) species ground cover (at peak bloom)
```{r}
# restrict dataset
groundcover_pb <- peak_bloom %>% dplyr::select(year, site, site_year, plot, block, groundcover_species_avg_peak)
groundcover_pb %>% group_by(plot) %>% tally()
# merge with plant list for names
groundcover_pb <- merge(groundcover_pb, plant_list, by = 'plot', all.x = T)
hist(groundcover_pb$groundcover_species_avg_peak, breaks=100)
groundcover_pb$groundcover_species_avg_peak2 <- groundcover_pb$groundcover_species_avg_peak/100
groundcover_pb$groundcover_species_avg_peak2[which(groundcover_pb$groundcover_species_avg_peak2 == 1)] <- 0.99999
M_groundcover_pb <- glmmTMB(groundcover_species_avg_peak2 ~ plant_name +
site +
year +
(1|site_year/block),
family = beta_family(),
data = groundcover_pb)
hist(residuals(M_groundcover_pb), breaks=100)
plot(simulateResiduals(M_groundcover_pb, n=1000))
testDispersion(M_groundcover_pb)
check_collinearity(M_groundcover_pb)
Anova(M_groundcover_pb, type=3)
# Chisq Df Pr(>Chisq)
# plant_name 106.6162 26 9.946e-12 ***
# site 0.2598 1 0.6103
# year 81.4294 1 < 2.2e-16 ***
performance(M_groundcover_pb)
# extracting marginal means - response scale
means_groundcover_pb <- as.data.frame(emmeans(M_groundcover_pb, specs = 'plant_name', type = "response"))
means_groundcover_pb$group <- 'groundcover_pb'
means_groundcover_pb <- means_groundcover_pb %>% rename('emmean' = 'response',
'lower.CL' = 'asymp.LCL',
'upper.CL' = 'asymp.UCL') %>%
dplyr::select(group, plant_name, emmean, lower.CL, upper.CL)
```
## (c) flower cover (at peak bloom)
```{r}
# restrict dataset
flowercover_pb <- droplevels(subset(peak_bloom, exclude == FALSE))
flowercover_pb <- flowercover_pb %>% dplyr::select(year, site, site_year, plot, block, floral_area_peak)
flowercover_pb %>% group_by(plot) %>% tally()
# merge with plant list for names
flowercover_pb <- merge(flowercover_pb, plant_list, by = 'plot', all.x = T)
hist(flowercover_pb$floral_area_peak, breaks=100)
min(flowercover_pb$floral_area_peak)
# model
M_flowercover_pb <- glmmTMB(floral_area_peak ~ plant_name +
site +
year +
(1|site_year/block),
family = nbinom2(link='log'),
data = flowercover_pb)
hist(residuals(M_flowercover_pb), breaks=100)
plot(simulateResiduals(M_flowercover_pb, n=1000))
testDispersion(M_flowercover_pb)
check_collinearity(M_flowercover_pb)
Anova(M_flowercover_pb, type=2)
# Chisq Df Pr(>Chisq)
# plant_name 106.6162 26 9.946e-12 ***
# site 0.2598 1 0.6103
# year 81.4294 1 < 2.2e-16 ***
performance(M_flowercover_pb)
# extracting marginal means - response scale
means_flowercover_pb <- as.data.frame(emmeans(M_flowercover_pb, specs = 'plant_name', type = "response"))
means_flowercover_pb$group <- 'flowercover_pb'
means_flowercover_pb <- means_flowercover_pb %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
```
## (d) flower cover (whole season)
```{r}
# restrict dataset
flowercover <- flower_cover %>% dplyr::select(year, site, site_year, plot, block, observations_is, observation_should, floral_area)
flowercover %>% group_by(plot) %>% tally()
# merge with plant list for names
flowercover <- merge(flowercover, plant_list, by = 'plot', all.x = T)
# model
M_flowercover <- glmmTMB(floral_area ~ plant_name +
site +
year +
(1|site_year/block),
family = nbinom2(),
data = flowercover)
hist(residuals(M_flowercover), breaks=100)
plot(simulateResiduals(M_flowercover, n=1000))
testDispersion(M_flowercover)
check_collinearity(M_flowercover)
Anova(M_flowercover, type=2)
# Chisq Df Pr(>Chisq)
# plant_name 106.6162 26 9.946e-12 ***
# site 0.2598 1 0.6103
# year 81.4294 1 < 2.2e-16 ***
performance(M_flowercover)
# extracting marginal means - response scale
means_flowercover <- as.data.frame(emmeans(M_flowercover, specs = 'plant_name', type = "response"))
means_flowercover$group <- 'flowercover'
means_flowercover <- means_flowercover %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
means_flowercover <- means_flowercover %>% dplyr::select(group, plant_name, response) %>%
pivot_wider(names_from = group, values_from = response)
```
## (e) species ground cover (whole season)
```{r}
# restrict dataset
ground_cover <- groundcover %>% dplyr::select(year, site, site_year, plot, block, groundcover_species_avg_season)
ground_cover %>% group_by(plot) %>% tally()
# merge with plant list for names
ground_cover <- merge(ground_cover, plant_list, by = 'plot', all.x = T)
ground_cover$groundcover_species_avg_season_prop <- ground_cover$groundcover_species_avg_season /100
ground_cover$groundcover_species_avg_season_prop[which(ground_cover$groundcover_species_avg_season_prop == 1)] <- 0.9999999
# model
M_groundcover<- glmmTMB(groundcover_species_avg_season_prop ~ plant_name +
site +
year +
(1|site_year/block),
family = beta_family(),
data = ground_cover)
hist(residuals(M_groundcover), breaks=100)
plot(simulateResiduals(M_groundcover, n=1000))
testDispersion(M_groundcover)
check_collinearity(M_groundcover)
Anova(M_groundcover, type=2)
# Chisq Df Pr(>Chisq)
# plant_name 106.6162 26 9.946e-12 ***
# site 0.2598 1 0.6103
# year 81.4294 1 < 2.2e-16 ***
summary(M_groundcover)
performance(M_groundcover)
# extracting marginal means - response scale
means_groundcover <- as.data.frame(emmeans(M_groundcover, specs = 'plant_name', type = "response"))
means_groundcover$group <- 'groundcover'
means_groundcover <- means_groundcover %>% dplyr::select(group, plant_name, response, asymp.LCL, asymp.UCL)
means_groundcover <- means_groundcover %>% dplyr::select(group, plant_name, response) %>%
pivot_wider(names_from = group, values_from = response)
```
## (f) combined table
```{r}
# peak bloom
means_peakbloom <- means_peakbloom %>% dplyr::select(plant_name, emmean) %>%
rename('peak_bloom_timing' = 'emmean')
# groundcover
means_groundcover <- means_groundcover %>% dplyr::select(plant_name, groundcover) %>%
rename('groundcover_whole_season' = 'groundcover')
means_groundcover_pb <- means_groundcover_pb %>% dplyr::select(plant_name, emmean) %>%
rename('groundcover_at_peak_bloom' = 'emmean')
# flower cover
means_flowercover <- means_flowercover %>% dplyr::select(plant_name, flowercover) %>%
rename('flowercover_whole_season' = 'flowercover')
means_flowercover_pb <- means_flowercover_pb %>% dplyr::select(plant_name, response) %>%
rename('flowercover_at_peak_bloom' = 'response')
# merge for models
means_traits <- means_peakbloom
means_traits <- merge(means_traits, means_groundcover_pb, by = 'plant_name')
means_traits <- merge(means_traits, means_flowercover_pb, by = 'plant_name')
means_traits <- merge(means_traits, means_flowercover, by = 'plant_name')
means_traits <- merge(means_traits, means_groundcover, by = 'plant_name')
```
#===============================================================================================#
# Modelling groups vs. traits
```{r}
# create table with all preditions
# merge with predictions for peak bloom, ground cover and flower cover
predictions <- means %>% dplyr::select(plant_name, group, response) %>% pivot_wider(names_from = group, values_from = response)
predictions <- merge(predictions, means_traits, by = 'plant_name')
# merge with plant life history
plant_lh <- plant_list %>% dplyr::select(plant_name, lifecycle)
predictions <- merge(predictions, plant_lh, by = 'plant_name')
predictions$lifecycle <- as.factor(predictions$lifecycle)
# initiate table to write values in
coefficients <- data.frame(response=c(NA,NA,NA,NA), parameter=c('life_cycle','peak_bloom', 'ground_cover','flower_cover'))
```
## (a) hoverflies
```{r}
M_hoverflies2 <- glm.nb(predictions[,9] ~
lifecycle +
scale(peak_bloom_timing) +
scale(groundcover_at_peak_bloom) +
scale(log(flowercover_at_peak_bloom)),
data = predictions)
hist(residuals(M_hoverflies2), breaks=20)
plot(simulateResiduals(M_hoverflies2),quantreg = T)
check_collinearity(M_hoverflies2)
Anova(M_hoverflies2, type=3)
# LR Chisq Df Pr(>Chisq)
# lifecycle 8.9943 1 0.002708 **
# scale(peakbloom) 4.0178 1 0.045022 *
# scale(groundcover) 4.7287 1 0.029664 *
# scale(log(flowercover)) 0.1831 1 0.668734
hovfly <- as.data.frame(model_parameters(M_hoverflies2, type='response'))
coefficients_hovfly <- coefficients
coefficients_hovfly$response <- 'hoverflies'
coefficients_hovfly$coef <- hovfly$Coefficient[2:5]
coefficients_hovfly$CI_low <- hovfly$CI_low[2:5]
coefficients_hovfly$CI_high <- hovfly$CI_high[2:5]
```
## (b) wild bees
```{r}
M_wildbees2<- glm.nb(predictions[,8] ~
lifecycle +
scale(peak_bloom_timing) +
scale(groundcover_at_peak_bloom) +
scale(log(flowercover_at_peak_bloom)),
data = predictions)
hist(residuals(M_wildbees2), breaks=20)
plot(simulateResiduals(M_wildbees2),quantreg = T)
check_collinearity(M_wildbees2)
Anova(M_wildbees2, type=3)
# LR Chisq Df Pr(>Chisq)
# lifecycle 2.05923 1 0.1513
# scale(peakbloom) 0.04094 1 0.8397
# scale(groundcover) 1.11375 1 0.2913
# scale(log(flowercover)) 1.30720 1 0.2529
wildbee <- as.data.frame(model_parameters(M_wildbees2, type='response'))
coefficients_wildbee <- coefficients
coefficients_wildbee$response <- 'wild_bees'
coefficients_wildbee$coef <- wildbee$Coefficient[2:5]
coefficients_wildbee$CI_low <- wildbee$CI_low[2:5]
coefficients_wildbee$CI_high <- wildbee$CI_high[2:5]
```
## (c) herbivorous arthropods
```{r}
M_herbivorous_arthropods2 <- glm.nb(predictions[,7] ~
lifecycle +
scale(peak_bloom_timing) +
scale(groundcover_at_peak_bloom) +
scale(log(flowercover_at_peak_bloom)),
data = predictions)
hist(residuals(M_herbivorous_arthropods2), breaks=20)
plot(simulateResiduals(M_herbivorous_arthropods2),quantreg = T)
check_collinearity(M_herbivorous_arthropods2)
Anova(M_herbivorous_arthropods2, type=3)
# LR Chisq Df Pr(>Chisq)
# lifecycle 1.92617 1 0.1652
# scale(peakbloom) 0.41869 1 0.5176
# scale(groundcover) 1.54629 1 0.2137
# scale(log(flowercover)) 0.23973 1 0.6244
leafherb <- as.data.frame(model_parameters(M_herbivorous_arthropods2, type='response'))
coefficients_leafherb <- coefficients
coefficients_leafherb$response <- 'leaf_herbivores'
coefficients_leafherb$coef <- leafherb$Coefficient[2:5]
coefficients_leafherb$CI_low <- leafherb$CI_low[2:5]
coefficients_leafherb$CI_high <- leafherb$CI_high[2:5]
```
## (d) parasitic wasps
```{r}
M_parasitic_wasps2 <- glm.nb(predictions[,6] ~
lifecycle +
scale(peak_bloom_timing) +
scale(groundcover_at_peak_bloom) +
scale(log(flowercover_at_peak_bloom)),
data = predictions)
hist(residuals(M_parasitic_wasps2), breaks=20)
plot(simulateResiduals(M_parasitic_wasps2),quantreg = T)
check_collinearity(M_parasitic_wasps2)
Anova(M_parasitic_wasps2, type=3)
# LR Chisq Df Pr(>Chisq)
# lifecycle 0.0131 1 0.90896
# scale(peakbloom) 3.3635 1 0.06665 .
# scale(groundcover) 0.3522 1 0.55287
# scale(log(flowercover)) 2.1862 1 0.13925
leafpara <- as.data.frame(model_parameters(M_parasitic_wasps2, type='response'))
coefficients_leafpara <- coefficients
coefficients_leafpara$response <- 'leaf_parasitoids'
coefficients_leafpara$coef <- leafpara$Coefficient[2:5]
coefficients_leafpara$CI_low <- leafpara$CI_low[2:5]
coefficients_leafpara$CI_high <- leafpara$CI_high[2:5]
```
## (e) predatory arthropods
```{r}
M_predatory_arthropods2 <- glm.nb(predictions[,5] ~
lifecycle +
scale(peak_bloom_timing) +
scale(groundcover_at_peak_bloom) +
scale(log(flowercover_at_peak_bloom)),
data = predictions)
hist(residuals(M_predatory_arthropods2), breaks=20)
plot(simulateResiduals(M_predatory_arthropods2),quantreg = T)
check_collinearity(M_predatory_arthropods2)
Anova(M_predatory_arthropods2, type=3)
# LR Chisq Df Pr(>Chisq)
# lifecycle 3.1645 1 0.07525 .
# scale(peakbloom) 1.8209 1 0.17721
# scale(groundcover) 0.6179 1 0.43184
# scale(log(flowercover)) 0.1378 1 0.71047
leafpred <- as.data.frame(model_parameters(M_predatory_arthropods2, type='response'))
coefficients_leafpred <- coefficients
coefficients_leafpred$response <- 'leaf_predators'
coefficients_leafpred$coef <- leafpred$Coefficient[2:5]
coefficients_leafpred$CI_low <- leafpred$CI_low[2:5]
coefficients_leafpred$CI_high <- leafpred$CI_high[2:5]
```
## (f) herbivorous nematodes
```{r}
M_herbivorous_nematodes2 <- glm.nb(predictions[,2] ~
lifecycle +
scale(peak_bloom_timing) +
scale(groundcover_at_peak_bloom) +
scale(log(flowercover_at_peak_bloom)),
data = predictions)
hist(residuals(M_herbivorous_nematodes2), breaks=20)
plot(simulateResiduals(M_herbivorous_nematodes2),quantreg = T)
check_collinearity(M_herbivorous_nematodes2)
Anova(M_herbivorous_nematodes2, type=3)
# LR Chisq Df Pr(>Chisq)
# lifecycle 0.1217 1 0.72717
# scale(peakbloom) 1.7179 1 0.18997
# scale(groundcover) 3.3647 1 0.06661 .
# scale(log(flowercover)) 2.8513 1 0.09130 .
summary(M_herbivorous_nematodes2)
# extracting coefficients - model scale
soilpf <- as.data.frame(model_parameters(M_herbivorous_nematodes2, type='response'))
coefficients_soilpf <- coefficients
coefficients_soilpf$response <- 'soil_plant_feeders'
coefficients_soilpf$coef <- soilpf$Coefficient[2:5]
coefficients_soilpf$CI_low <- soilpf$CI_low[2:5]
coefficients_soilpf$CI_high <- soilpf$CI_high[2:5]
```
## (g) predatory nematodes
```{r}
M_predatory_nematodes2 <- glm(predictions[,4] ~
lifecycle +
scale(peak_bloom_timing) +
scale(groundcover_at_peak_bloom) +
scale(log(flowercover_at_peak_bloom)),
family = gaussian(link = 'log'),
data = predictions)
hist(residuals(M_predatory_nematodes2), breaks=20)
plot(simulateResiduals(M_predatory_nematodes2),quantreg = T)
check_collinearity(M_predatory_nematodes2)
Anova(M_predatory_nematodes2, type=3)
# LR Chisq Df Pr(>Chisq)
# lifecycle 2.20521 1 0.13754
# scale(peakbloom) 2.83597 1 0.09218 .
# scale(groundcover) 0.51757 1 0.47188
# scale(log(flowercover)) 1.17282 1 0.27882
soilpred <- as.data.frame(model_parameters(M_predatory_nematodes2, type='response'))
coefficients_soilpred <- coefficients
coefficients_soilpred$response <- 'soil_predators'
coefficients_soilpred$coef <- soilpred$Coefficient[2:5]
coefficients_soilpred$CI_low <- soilpred$CI_low[2:5]
coefficients_soilpred$CI_high <- soilpred$CI_high[2:5]
```
## (h) decomposer nematodes
```{r}
M_decomposer_nematodes2 <- glm.nb(predictions[,3] ~
lifecycle +
scale(peak_bloom_timing) +
scale(groundcover_at_peak_bloom) +
scale(log(flowercover_at_peak_bloom)),
data = predictions)
hist(residuals(M_decomposer_nematodes2), breaks=20)
plot(simulateResiduals(M_decomposer_nematodes2),quantreg = T)
check_collinearity(M_decomposer_nematodes2)
Anova(M_decomposer_nematodes2, type=3)
# lifecycle 15.5064 1 8.223e-05 ***
# scale(peakbloom) 0.2138 1 0.6438
# scale(groundcover) 2.3504 1 0.1253
# scale(log(flowercover)) 0.4869 1 0.4853
# extracting coefficients - model scale
soildec <- as.data.frame(model_parameters(M_decomposer_nematodes2, type='response'))
coefficients_soildec <- coefficients
coefficients_soildec$response <- 'soil_decomposers'
coefficients_soildec$coef <- soildec$Coefficient[2:5]
coefficients_soildec$CI_low <- soildec$CI_low[2:5]
coefficients_soildec$CI_high <- soildec$CI_high[2:5]
```
## (j) combined table
```{r}
coefficients_comb <- rbind(coefficients_soilpf,
coefficients_soildec,
coefficients_soilpred,
coefficients_leafpred,
coefficients_leafpara,
coefficients_leafherb,
coefficients_wildbee,
coefficients_hovfly)
```
# Figure 2
```{r}
# Define plotting aesthetics
plot_theme2 <- theme(axis.title = element_text(color="black", face='bold', size=14),
axis.text.x = element_text(color="black", size=12),
axis.text.y = element_text(color="black", size=12),
panel.border = element_rect(colour = "black", fill=NA, size=1),
panel.background = element_rect(colour = NA, fill = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "bottom",
legend.title = element_blank(),
legend.text = element_text(size=12),
legend.background=element_blank(),
plot.title = element_text(size=20, face="bold",hjust = 0.5))
# reorder factors to ensure correct order
coefficients_comb$response <- ordered(coefficients_comb$response, levels =c('hoverflies',
'wild_bees',
'leaf_herbivores',
'leaf_parasitoids',
'leaf_predators',
'soil_plant_feeders',
'soil_predators',
'soil_decomposers'))
# rename factor levels for plotting
coefficients_comb <- coefficients_comb %>%
mutate(parameter = as.factor(parameter),
parameter = dplyr::recode(parameter, 'flower_cover' = 'floral area',
'peak_bloom' = 'peak bloom',
'life_cycle' = 'life cycle',
'ground_cover' = 'plant cover'),
response = dplyr::recode(response, 'hoverflies' = 'hoverflies',
'wild_bees' = 'wild bees',
'leaf_herbivores' = 'herbivorous arthropods',
'leaf_parasitoids' = 'parasitic wasps',
'leaf_predators' = 'predatory arthropods',
'soil_plant_feeders' = 'herbivorous nematodes',
'soil_predators' = 'predatory nematodes',
'soil_decomposers' = 'decomposer nematodes'))
coefficients_comb$parameter <- ordered(coefficients_comb$parameter, levels =c('floral area', 'peak bloom', 'life cycle', 'plant cover'))
# actual plot
all_coef <- ggplot(coefficients_comb, aes(x = response, y = coef, color = response, fill = response)) +
geom_hline(yintercept = 0, linetype="dashed", color = "black", linewidth=0.5)+
geom_pointrange(aes(ymin=CI_low, ymax=CI_high), size=0.75, linewidth= 1.15, alpha = 1,
shape=21, key_glyph = 'rect') +
labs(y = " ", x = " ", limits=c(-0.75,1.0), fill = "parameter") +
scale_colour_viridis(discrete=TRUE, guide = "none") +
facet_grid( ~ parameter, scales = 'free_y', switch = 'y') +
scale_y_continuous(position="right") +
theme_bw() +
plot_theme2 +
theme(strip.text = element_text(size = 13),
axis.text.x=element_blank(),
axis.ticks.x = element_blank(),
axis.text.y=element_text(size=9),
strip.text.y.left = element_text(angle = 90),
strip.background = element_rect(colour=NA, fill="white"))
# export plot
tiff("Figure_2.tiff", res=600, width=220, height=100, units="mm", compression = "lzw")
all_coef
dev.off()
```
#===============================================================================================#
# Corrleations
```{r}
# load required package
library(metan) # 1.18.0
```
## Figure 3
```{r}
# select data for correlations
corrs <- predictions[,c(2:9)] %>%
rename('decomposer nematodes' = 'soil decomposers',
'herbivorous nematodes' = 'soil plant feeders',
'predatory nematodes' = 'soil predators',
'predatory arthropods' = 'leaf predators',
'parasitic wasps' = 'leaf parasitoids',
'herbivorous arthropods' = 'leaf herbivores',
'wild bees' = 'wild bees',
'hoverflies' = 'hoverflies')
# generate correlation matrix
CM <- corr_coef(corrs)
# generate plot
Figure_3 <- plot(CM, reorder = FALSE, digits.cor = 2) +
theme(axis.text.x = element_text(color = "black",size = 10, face = "bold"),
axis.text.y = element_text(color = "black",size = 10, face = "bold"),
plot.margin = margin(0.6,0.6,0.6,0.6, "cm"))+
scale_fill_gradient2(low = "#6D9EC1", mid = "white", high ="#E46726", midpoint = 0,
limits = c(-1,1), space = "Lab",
name="Pearson correlation\ncoefficient r")
# export plot
tiff("Figure_3.tiff", res=600, width=230, height=170, units="mm", compression = "lzw")
Figure_3
dev.off()
```
## Figure S1
```{r}
# select data for correlations
corrs2 <- means_traits[,-1] %>%
rename('floral area\nat peak bloom' = 'flowercover_at_peak_bloom',
'floral area' = 'flowercover_whole_season',
'plant cover\nat peak bloom' = 'groundcover_at_peak_bloom',
'plant cover' = 'groundcover_whole_season',
'peak bloom' = 'peak_bloom_timing')
# generate correlation matrix
CM2 <- corr_coef(corrs2)
# generate plot
Figure_S1 <- plot(CM2, reorder = FALSE, digits.cor = 2) +
theme(axis.text.x = element_text(color = "black",size = 10, face = "bold"),
axis.text.y = element_text(color = "black",size = 10, face = "bold"),
plot.margin = margin(0.6,0.6,0.6,0.6, "cm"))+
scale_fill_gradient2(low = "#6D9EC1", mid = "white", high ="#E46726", midpoint = 0,
limits = c(-1,1), space = "Lab",
name="Pearson correlation\ncoefficient r")
# export plot
tiff("Figure_S1.tiff", res=600, width=230, height=170, units="mm", compression = "lzw")
Figure_S1
dev.off()
```
####
```{r}
```