flower_plots
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
#===============================================================================================# # Loading required packages and defining plotting aesthetics for figures
#===============================================================================================# # Load datasets & data handling
#===============================================================================================# # exclude samples (based on inclusion criteria, see main text)
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
# restrict dataset
hoverflies <- pollinators %>% dplyr::select(year, site, site_year, plot, block, hoverflies, observations_is, observations_should)
hoverflies %>% group_by(plot) %>% tally()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 12
## 4 6 8
## 5 7 12
## 6 8 11
## 7 9 11
## 8 10 12
## 9 11 15
## 10 12 12
## # ℹ 17 more rows# 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) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.56454, p-value = 0.776
## alternative hypothesis: two.sidedcheck_collinearity(M_hoverflies)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.03 [1.00, 2.11] 1.01 0.97 [0.47, 1.00]
## site 1.00 [1.00, 7.16e+07] 1.00 1.00 [0.00, 1.00]
## year 1.02 [1.00, 2.73] 1.01 0.98 [0.37, 1.00]check_zeroinflation(M_hoverflies)## # Check for zero-inflation
##
## Observed zeros: 46
## Predicted zeros: 43
## Ratio: 0.93## Model seems ok, ratio of observed and predicted zeros is within the tolerance range (p = 0.792).Anova(M_hoverflies, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: hoverflies
## Chisq Df Pr(>Chisq)
## (Intercept) 39.0604 1 4.109e-10 ***
## plant_name 262.8988 26 < 2.2e-16 ***
## site 0.6197 1 0.4311
## year 0.0380 1 0.8454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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))))## [1] 287.4998# 288.0297
# summary(M_hoverflies)
performance(M_hoverflies)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | Score_log | Score_spherical
## ---------------------------------------------------------------------------------------------------------------
## 1992.139 | 2000.419 | 2114.909 | 1.000 | 0.640 | 1.000 | 17.355 | | -3.194 | 0.031# 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
# 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()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 12
## 4 6 8
## 5 7 12
## 6 8 11
## 7 9 11
## 8 10 12
## 9 11 15
## 10 12 12
## # ℹ 17 more rows# 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) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.60123, p-value = 0.6
## alternative hypothesis: two.sidedcheck_collinearity(M_wild_bees)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.21 [1.10, 1.41] 1.10 0.83 [0.71, 0.91]
## site 1.04 [1.00, 1.65] 1.02 0.96 [0.61, 1.00]
## year 1.16 [1.07, 1.37] 1.08 0.86 [0.73, 0.93]Anova(M_wild_bees, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: wild_bees
## 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 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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))))## [1] 288.0297# 288.0297
# summary(M_wild_bees)
performance(M_wild_bees)## # Indices of model performance
##
## AIC | AICc | BIC | RMSE | Sigma | Score_log | Score_spherical
## -----------------------------------------------------------------------------
## 1960.323 | 1968.088 | 2079.373 | 32.852 | 1.316 | -3.135 | 0.031# 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
# restrict dataset
leaf_herbivores <- leafsuction %>% dplyr::select(year, site, site_year, plot, block, herbivores, observations)
leaf_herbivores %>% group_by(plot) %>% tally()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 12
## 4 6 8
## 5 7 11
## 6 8 11
## 7 9 11
## 8 10 12
## 9 11 15
## 10 12 12
## # ℹ 17 more rows# 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) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.2088, p-value = 0.472
## alternative hypothesis: two.sidedcheck_collinearity(M_leaf_herbivores)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.05 [1.01, 1.46] 1.02 0.95 [0.69, 0.99]
## site 1.01 [1.00, 21.78] 1.01 0.99 [0.05, 1.00]
## year 1.03 [1.00, 1.74] 1.02 0.97 [0.57, 1.00]Anova(M_leaf_herbivores, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: herbivores
## Chisq Df Pr(>Chisq)
## (Intercept) 535.2252 1 < 2.2e-16 ***
## plant_name 265.8028 26 < 2.2e-16 ***
## site 0.0380 1 0.845526
## year 8.1479 1 0.004311 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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))))## [1] 280.4249# 282.0197
# summary(M_leaf_herbivores)
performance(M_leaf_herbivores)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | Score_log | Score_spherical
## ----------------------------------------------------------------------------------------------------------------
## 3156.605 | 3166.187 | 3286.121 | 1.000 | 0.809 | 1.000 | 132.920 | | -5.175 | 0.041# 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
# restrict dataset
leaf_parasitoids <- leafsuction %>% dplyr::select(year, site, site_year, plot, block, parasitoids, observations)
leaf_parasitoids %>% group_by(plot) %>% tally()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 12
## 4 6 8
## 5 7 11
## 6 8 11
## 7 9 11
## 8 10 12
## 9 11 15
## 10 12 12
## # ℹ 17 more rows# 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) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.2681, p-value = 0.432
## alternative hypothesis: two.sidedcheck_collinearity(M_leaf_parasitoids)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.03 [1.00, 2.06] 1.01 0.97 [0.48, 1.00]
## site 1.01 [1.00, 59.78] 1.01 0.99 [0.02, 1.00]
## year 1.02 [1.00, 7.57] 1.01 0.98 [0.13, 1.00]Anova(M_leaf_parasitoids, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: parasitoids
## Chisq Df Pr(>Chisq)
## (Intercept) 154.6180 1 < 2e-16 ***
## plant_name 251.0650 26 < 2e-16 ***
## site 0.0092 1 0.92370
## year 3.7247 1 0.05361 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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))))## [1] 282.0197# 282.0197
# summary(M_leaf_parasitoids)
performance(M_leaf_parasitoids)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | Score_log | Score_spherical
## ---------------------------------------------------------------------------------------------------------------
## 2562.015 | 2569.955 | 2680.429 | 0.663 | 0.474 | 0.360 | 44.075 | 2.232 | -4.136 | 0.040# 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
# restrict dataset
leaf_predators <- leafsuction %>% dplyr::select(year, site, site_year, plot, block, predators, observations) %>% droplevels()
leaf_predators %>% group_by(plot) %>% tally()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 12
## 4 6 8
## 5 7 11
## 6 8 11
## 7 9 11
## 8 10 12
## 9 11 15
## 10 12 12
## # ℹ 17 more rows# 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) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.73606, p-value = 0.304
## alternative hypothesis: two.sidedcheck_collinearity(M_leaf_predators)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.26 [1.14, 1.47] 1.12 0.79 [0.68, 0.87]
## site 1.09 [1.02, 1.34] 1.04 0.92 [0.75, 0.98]
## year 1.16 [1.07, 1.37] 1.08 0.86 [0.73, 0.94]Anova(M_leaf_predators, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: predators
## 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 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # 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))))## [1] 282.0197# 282.0197
# summary(M_leaf_predators)
performance(M_leaf_predators)## # Indices of model performance
##
## AIC | AICc | BIC | RMSE | Sigma | Score_log | Score_spherical
## -----------------------------------------------------------------------------
## 2502.067 | 2510.007 | 2620.481 | 20.596 | 2.230 | -4.093 | 0.045# 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
# 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()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 11
## 4 6 16
## 5 7 12
## 6 8 12
## 7 9 11
## 8 10 12
## 9 11 14
## 10 12 11
## # ℹ 17 more rows# 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)## Warning in glmmTMB(plant_feeders ~ plant_name + site + year + offset(log(soil_weight_g)) + : non-integer counts in a nbinom2 modelhist(residuals(M_soil_plant_feeders), breaks=100)
plot(simulateResiduals(M_soil_plant_feeders, n=1000)) 
testDispersion(M_soil_plant_feeders) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 2.1518, p-value = 0.072
## alternative hypothesis: two.sidedcheck_collinearity(M_soil_plant_feeders)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.35 [1.21, 1.56] 1.16 0.74 [0.64, 0.83]
## site 1.14 [1.06, 1.35] 1.07 0.88 [0.74, 0.95]
## year 1.21 [1.10, 1.41] 1.10 0.83 [0.71, 0.91]Anova(M_soil_plant_feeders, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: plant_feeders
## 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 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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))))## [1] 292.0361# 292.0361
summary(M_soil_plant_feeders)## Family: nbinom2 ( log )
## Formula: plant_feeders ~ plant_name + site + year + offset(log(soil_weight_g)) + (1 | site_year/block)
## Data: soil_plant_feeders
##
## AIC BIC logLik deviance df.resid
## 2346.5 2466.0 -1141.3 2282.5 277
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## block:site_year (Intercept) 4.791e-02 0.218875
## site_year (Intercept) 9.981e-06 0.003159
## Number of obs: 309, groups: block:site_year, 16; site_year, 4
##
## Dispersion parameter for nbinom2 family (): 1.03
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.04566 0.08643 -0.528 0.597300
## plant_name1 0.07318 0.29697 0.246 0.805355
## plant_name2 -1.06244 0.32324 -3.287 0.001013 **
## plant_name3 0.12442 0.29512 0.422 0.673324
## plant_name4 0.00012 0.49976 0.000 0.999808
## plant_name5 0.08456 0.29953 0.282 0.777718
## plant_name6 -0.73672 0.28215 -2.611 0.009025 **
## plant_name7 -0.43181 0.26312 -1.641 0.100767
## plant_name8 -0.48095 0.31845 -1.510 0.130979
## plant_name9 -0.40653 0.29861 -1.361 0.173376
## plant_name10 1.17211 0.35143 3.335 0.000852 ***
## plant_name11 1.03211 0.29702 3.475 0.000511 ***
## plant_name12 0.23327 0.29975 0.778 0.436437
## plant_name13 -0.11660 0.29556 -0.395 0.693196
## plant_name14 -0.17434 0.31665 -0.551 0.581910
## plant_name15 -0.89503 0.30781 -2.908 0.003640 **
## plant_name16 0.27657 0.29334 0.943 0.345756
## plant_name17 0.75002 0.25095 2.989 0.002801 **
## plant_name18 0.35402 0.25294 1.400 0.161622
## plant_name19 -1.34808 0.27994 -4.816 1.47e-06 ***
## plant_name20 1.03455 0.35178 2.941 0.003273 **
## plant_name21 -0.09667 0.37069 -0.261 0.794270
## plant_name22 1.56021 0.30261 5.156 2.52e-07 ***
## plant_name23 -0.30673 0.37263 -0.823 0.410423
## plant_name24 0.07189 0.26465 0.272 0.785902
## plant_name25 -0.62876 0.37512 -1.676 0.093713 .
## plant_name26 0.28381 0.38717 0.733 0.463541
## site1 -0.41060 0.08817 -4.657 3.21e-06 ***
## year1 -0.36599 0.09009 -4.062 4.86e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1performance(M_soil_plant_feeders)## # Indices of model performance
##
## AIC | AICc | BIC | RMSE | Sigma | Score_log | Score_spherical
## -----------------------------------------------------------------------------
## 2346.542 | 2354.194 | 2466.009 | 42.427 | 1.034 | -Inf | 0.022# 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
# 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()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 11
## 4 6 16
## 5 7 12
## 6 8 12
## 7 9 11
## 8 10 12
## 9 11 14
## 10 12 11
## # ℹ 17 more rows# 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)## Warning in glmmTMB(predators ~ plant_name + site + year + offset(log(soil_weight_g)) + : non-integer counts in a nbinom2 modelhist(residuals(M_soil_predators), breaks=100)
plot(simulateResiduals(M_soil_predators, n=1000)) 
testDispersion(M_soil_predators) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.65386, p-value = 0.88
## alternative hypothesis: two.sidedcheck_collinearity(M_soil_predators)## # Check for Multicollinearity
##
## * conditional component:
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.03 [1.00, 1.95] 1.02 0.97 [0.51, 1.00]
## site 1.01 [1.00, 19.11] 1.01 0.99 [0.05, 1.00]
## year 1.02 [1.00, 2.56] 1.01 0.98 [0.39, 1.00]Anova(M_soil_predators, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: predators
## Chisq Df Pr(>Chisq)
## (Intercept) 21.0437 1 4.489e-06 ***
## plant_name 51.3991 26 0.0021321 **
## site 12.4786 1 0.0004116 ***
## year 2.2502 1 0.1336003
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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))))## [1] 289.9123# 292.0361
# summary(M_soil_predators)
performance(M_soil_predators)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | Score_log | Score_spherical
## --------------------------------------------------------------------------------------------------------------
## 1563.526 | 1573.320 | 1697.926 | 1.000 | 0.737 | 1.000 | 5.656 | | -Inf | 0.031# 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
# 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()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 11
## 4 6 16
## 5 7 12
## 6 8 12
## 7 9 11
## 8 10 12
## 9 11 14
## 10 12 11
## # ℹ 17 more rows# 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)## Warning in glmmTMB(decomposers ~ plant_name + site + year + offset(log(soil_weight_g)) + : non-integer counts in a nbinom2 modelhist(residuals(M_soil_decomposers), breaks=100)
plot(simulateResiduals(M_soil_decomposers, n=1000)) 
testDispersion(M_soil_decomposers) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.1157, p-value = 0.488
## alternative hypothesis: two.sidedcheck_collinearity(M_soil_decomposers)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.22 [1.12, 1.43] 1.11 0.82 [0.70, 0.89]
## site 1.10 [1.03, 1.33] 1.05 0.91 [0.75, 0.97]
## year 1.13 [1.05, 1.34] 1.06 0.89 [0.75, 0.96]Anova(M_soil_decomposers, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: decomposers
## 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 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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))))## [1] 292.0361# 292.0361
# summary(M_soil_decomposers)
performance(M_soil_decomposers)## # Indices of model performance
##
## AIC | AICc | BIC | RMSE | Sigma | Score_log | Score_spherical
## ------------------------------------------------------------------------------
## 3772.025 | 3779.677 | 3891.492 | 145.501 | 2.198 | -Inf | 0.025# 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
# 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"))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.# export plot
tiff("Figure_1.tiff", res=800, width=340, height=240, units="mm", compression = "lzw")
radarplot
dev.off()## png
## 2#===============================================================================================# # Modelling traits
(a) peak bloom
# restrict dataset
peakbloom <- peak_bloom %>% dplyr::select(year, site, site_year, plot, block, peak_avg)
peakbloom %>% group_by(plot) %>% tally()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 12
## 4 6 8
## 5 7 8
## 6 8 12
## 7 9 12
## 8 10 12
## 9 11 16
## 10 12 12
## # ℹ 17 more rows# 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) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.0454, p-value = 0.648
## alternative hypothesis: two.sidedplotResiduals(simulateResiduals(M_peakbloom, n=1000), form=peakbloom$site)
check_collinearity(M_peakbloom)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.03 [1.00, 2.02] 1.01 0.97 [0.50, 1.00]
## site 1.01 [1.00, 1687.62] 1.00 0.99 [0.00, 1.00]
## year 1.02 [1.00, 3.53] 1.01 0.98 [0.28, 1.00]plotResiduals(simulateResiduals(M_peakbloom, n=1000), form=peakbloom$site_year)
Anova(M_peakbloom, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: peak_avg
## Chisq Df Pr(>Chisq)
## (Intercept) 7658.0326 1 < 2.2e-16 ***
## plant_name 338.9022 26 < 2.2e-16 ***
## site 0.5613 1 0.4537
## year 45.9359 1 1.222e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## --------------------------------------------------------------------------------
## 1260.229 | 1267.881 | 1379.696 | 0.792 | 0.758 | 0.143 | 1.639 | 1.649# 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)
# restrict dataset
groundcover_pb <- peak_bloom %>% dplyr::select(year, site, site_year, plot, block, groundcover_species_avg_peak)
groundcover_pb %>% group_by(plot) %>% tally()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 12
## 4 6 8
## 5 7 8
## 6 8 12
## 7 9 12
## 8 10 12
## 9 11 16
## 10 12 12
## # ℹ 17 more rows# 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) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.051, p-value = 0.592
## alternative hypothesis: two.sidedcheck_collinearity(M_groundcover_pb)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.49 [1.33, 1.74] 1.22 0.67 [0.58, 0.75]
## site 1.10 [1.03, 1.33] 1.05 0.91 [0.75, 0.97]
## year 1.37 [1.23, 1.59] 1.17 0.73 [0.63, 0.81]Anova(M_groundcover_pb, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: groundcover_species_avg_peak2
## Chisq Df Pr(>Chisq)
## (Intercept) 182.1048 1 < 2.2e-16 ***
## plant_name 123.2613 26 1.329e-14 ***
## site 1.6387 1 0.2005
## year 0.1052 1 0.7456
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## Random effect variances not available. Returned R2 does not account for random effects.## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | RMSE | Sigma
## ------------------------------------------------------------------------
## -327.509 | -319.857 | -208.042 | | 0.860 | 0.203 | 3.470# 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)
# 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()## # A tibble: 27 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 3 4
## 3 4 12
## 4 6 8
## 5 7 8
## 6 8 11
## 7 9 11
## 8 10 12
## 9 11 15
## 10 12 12
## # ℹ 17 more rows# 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)## [1] 0.38465# model
M_flowercover_pb <- glmmTMB(floral_area_peak ~ plant_name +
site +
year +
(1|site_year/block),
family = nbinom2(link='log'),
data = flowercover_pb)## Warning in glmmTMB(floral_area_peak ~ plant_name + site + year + (1 | site_year/block), : non-integer counts in a nbinom2 modelhist(residuals(M_flowercover_pb), breaks=100)
plot(simulateResiduals(M_flowercover_pb, n=1000)) 
testDispersion(M_flowercover_pb) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.39339, p-value = 0.84
## alternative hypothesis: two.sidedcheck_collinearity(M_flowercover_pb)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.03 [1.00, 1.75] 1.02 0.97 [0.57, 1.00]
## site 1.01 [1.00, 1.06e+05] 1.00 0.99 [0.00, 1.00]
## year 1.03 [1.00, 2.15] 1.01 0.97 [0.46, 1.00]Anova(M_flowercover_pb, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: floral_area_peak
## Chisq Df Pr(>Chisq)
## plant_name 431.3469 26 <2e-16 ***
## site 1.7417 1 0.1869
## year 0.3352 1 0.5626
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## # Indices of model performance
##
## AIC | AICc | BIC | RMSE | Sigma | Score_log
## --------------------------------------------------------------
## 5283.970 | 5291.735 | 5403.020 | 16022.439 | 0.840 | -Inf# 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)
# restrict dataset
flowercover <- flower_cover %>% dplyr::select(year, site, site_year, plot, block, observations_is, observation_should, floral_area)
flowercover %>% group_by(plot) %>% tally()## # A tibble: 30 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 2 5
## 3 3 9
## 4 4 12
## 5 5 8
## 6 6 8
## 7 7 12
## 8 8 11
## 9 9 11
## 10 10 12
## # ℹ 20 more rows# 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)## Warning in glmmTMB(floral_area ~ plant_name + site + year + (1 | site_year/block), : non-integer counts in a nbinom2 modelhist(residuals(M_flowercover), breaks=100)
plot(simulateResiduals(M_flowercover, n=1000)) 
testDispersion(M_flowercover) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.35249, p-value = 0.928
## alternative hypothesis: two.sidedcheck_collinearity(M_flowercover)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.02 [1.00, 3.13] 1.01 0.98 [0.32, 1.00]
## site 1.00 [1.00, 1.62e+09] 1.00 1.00 [0.00, 1.00]
## year 1.02 [1.00, 5.31] 1.01 0.98 [0.19, 1.00]Anova(M_flowercover, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: floral_area
## Chisq Df Pr(>Chisq)
## plant_name 529.7056 29 <2e-16 ***
## site 1.6651 1 0.1969
## year 0.3459 1 0.5564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## # Indices of model performance
##
## AIC | AICc | BIC | RMSE | Sigma | Score_log
## --------------------------------------------------------------
## 6723.212 | 6730.989 | 6859.225 | 26993.826 | 0.901 | -Inf# 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)
# restrict dataset
ground_cover <- groundcover %>% dplyr::select(year, site, site_year, plot, block, groundcover_species_avg_season)
ground_cover %>% group_by(plot) %>% tally()## # A tibble: 30 × 2
## plot n
## <fct> <int>
## 1 1 12
## 2 2 12
## 3 3 16
## 4 4 12
## 5 5 12
## 6 6 16
## 7 7 12
## 8 8 12
## 9 9 16
## 10 10 12
## # ℹ 20 more rows# 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) 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.1134, p-value = 0.152
## alternative hypothesis: two.sidedcheck_collinearity(M_groundcover)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## plant_name 1.04 [1.00, 1.44] 1.02 0.96 [0.70, 1.00]
## site 1.02 [1.00, 3.69] 1.01 0.98 [0.27, 1.00]
## year 1.03 [1.00, 2.00] 1.01 0.98 [0.50, 1.00]Anova(M_groundcover, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: groundcover_species_avg_season_prop
## Chisq Df Pr(>Chisq)
## plant_name 303.5391 29 <2e-16 ***
## site 0.4997 1 0.4796
## year 0.9658 1 0.3257
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## Family: beta ( logit )
## Formula: groundcover_species_avg_season_prop ~ plant_name + site + year + (1 | site_year/block)
## Data: ground_cover
##
## AIC BIC logLik deviance df.resid
## -249.2 -108.8 159.6 -319.2 373
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## block:site_year (Intercept) 0.02737 0.1654
## site_year (Intercept) 0.06500 0.2549
## Number of obs: 408, groups: block:site_year, 16; site_year, 4
##
## Dispersion parameter for beta family (): 3.39
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.25960 0.14534 1.786 0.074071 .
## plant_name1 0.22968 0.26819 0.856 0.391772
## plant_name2 0.23878 0.22639 1.055 0.291550
## plant_name3 -2.27310 0.29619 -7.675 1.66e-14 ***
## plant_name4 0.15721 0.26806 0.586 0.557557
## plant_name5 -1.49951 0.24761 -6.056 1.40e-09 ***
## plant_name6 1.69511 0.30627 5.535 3.12e-08 ***
## plant_name7 0.31899 0.22862 1.395 0.162923
## plant_name8 0.31326 0.23998 1.305 0.191773
## plant_name9 -0.40572 0.27789 -1.460 0.144293
## plant_name10 1.01073 0.28434 3.555 0.000379 ***
## plant_name11 -0.16391 0.24729 -0.663 0.507435
## plant_name12 2.25666 0.31358 7.196 6.18e-13 ***
## plant_name13 1.91986 0.31074 6.178 6.48e-10 ***
## plant_name14 -1.64857 0.28018 -5.884 4.01e-09 ***
## plant_name15 -0.03670 0.28356 -0.129 0.897033
## plant_name16 -0.01750 0.22320 -0.078 0.937521
## plant_name17 0.05114 0.22733 0.225 0.821996
## plant_name18 0.24236 0.23397 1.036 0.300254
## plant_name19 0.36838 0.24011 1.534 0.124972
## plant_name20 0.87536 0.25595 3.420 0.000626 ***
## plant_name21 0.82656 0.24124 3.426 0.000612 ***
## plant_name22 -0.80280 0.25792 -3.113 0.001855 **
## plant_name23 -0.61346 0.27557 -2.226 0.026002 *
## plant_name24 0.34927 0.23245 1.503 0.132950
## plant_name25 -0.56022 0.27249 -2.056 0.039789 *
## plant_name26 0.18576 0.23106 0.804 0.421424
## plant_name27 -0.26367 0.26878 -0.981 0.326608
## plant_name28 -1.20208 0.27026 -4.448 8.68e-06 ***
## plant_name29 -1.75479 0.29288 -5.991 2.08e-09 ***
## site1 0.10248 0.14497 0.707 0.479612
## year1 -0.14298 0.14549 -0.983 0.325728
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1performance(M_groundcover)## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## --------------------------------------------------------------------------------
## -249.239 | -242.465 | -108.845 | 0.865 | 0.791 | 0.351 | 0.220 | 3.392# 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
# 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
# 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
M_hoverflies2 <- glm.nb(predictions[,9] ~
lifecycle +
scale(peak_bloom_timing) +
scale(groundcover_at_peak_bloom) +
scale(log(flowercover_at_peak_bloom)),
data = predictions)## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.718217## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.702049## Warning in dpois(y, mu, log = TRUE): non-integer x = 6.755978## Warning in dpois(y, mu, log = TRUE): non-integer x = 3.928072## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.791931## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.581773## Warning in dpois(y, mu, log = TRUE): non-integer x = 9.706595## Warning in dpois(y, mu, log = TRUE): non-integer x = 38.536064## Warning in dpois(y, mu, log = TRUE): non-integer x = 4.838491## Warning in dpois(y, mu, log = TRUE): non-integer x = 31.384490## Warning in dpois(y, mu, log = TRUE): non-integer x = 50.961580## Warning in dpois(y, mu, log = TRUE): non-integer x = 7.399208## Warning in dpois(y, mu, log = TRUE): non-integer x = 0.592068## Warning in dpois(y, mu, log = TRUE): non-integer x = 7.327761## Warning in dpois(y, mu, log = TRUE): non-integer x = 4.199072## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.397437## Warning in dpois(y, mu, log = TRUE): non-integer x = 3.116650## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.353379## Warning in dpois(y, mu, log = TRUE): non-integer x = 14.394604## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.228516## Warning in dpois(y, mu, log = TRUE): non-integer x = 1.683211## Warning in dpois(y, mu, log = TRUE): non-integer x = 5.509222## Warning in dpois(y, mu, log = TRUE): non-integer x = 1.705361## Warning in dpois(y, mu, log = TRUE): non-integer x = 3.844724## Warning in dpois(y, mu, log = TRUE): non-integer x = 6.094923## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.813701## Warning in dpois(y, mu, log = TRUE): non-integer x = 3.281113hist(residuals(M_hoverflies2), breaks=20)
plot(simulateResiduals(M_hoverflies2),quantreg = T) 
check_collinearity(M_hoverflies2)## Warning: Using indexed data frames, such as `df[, 5]`, as model response can produce unexpected results. Specify your model using the literal name of
## the response variable instead.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## lifecycle 1.06 [1.00, 12.74] 1.03 0.94 [0.08, 1.00]
## scale(peak_bloom_timing) 1.23 [1.04, 2.34] 1.11 0.81 [0.43, 0.96]
## scale(groundcover_at_peak_bloom) 1.10 [1.00, 3.95] 1.05 0.91 [0.25, 1.00]
## scale(log(flowercover_at_peak_bloom)) 1.21 [1.03, 2.39] 1.10 0.83 [0.42, 0.97]Anova(M_hoverflies2, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: predictions[, 9]
## LR Chisq Df Pr(>Chisq)
## lifecycle 3.9108 1 0.04798 *
## scale(peak_bloom_timing) 3.1549 1 0.07570 .
## scale(groundcover_at_peak_bloom) 4.9259 1 0.02646 *
## scale(log(flowercover_at_peak_bloom)) 0.0196 1 0.88862
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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
M_wildbees2<- glm.nb(predictions[,8] ~
lifecycle +
scale(peak_bloom_timing) +
scale(groundcover_at_peak_bloom) +
scale(log(flowercover_at_peak_bloom)),
data = predictions)## Warning in dpois(y, mu, log = TRUE): non-integer x = 0.710281## Warning in dpois(y, mu, log = TRUE): non-integer x = 2.468564## Warning in dpois(y, mu, log = TRUE): non-integer x = 55.764439## Warning in dpois(y, mu, log = TRUE): non-integer x = 0.968623## Warning in dpois(y, mu, log = TRUE): non-integer x = 58.250825## Warning in dpois(y, mu, log = TRUE): non-integer x = 43.987527## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.314788## Warning in dpois(y, mu, log = TRUE): non-integer x = 2.256526## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.722237## Warning in dpois(y, mu, log = TRUE): non-integer x = 120.253788## Warning in dpois(y, mu, log = TRUE): non-integer x = 6.224524## Warning in dpois(y, mu, log = TRUE): non-integer x = 0.114970## Warning in dpois(y, mu, log = TRUE): non-integer x = 24.556783## Warning in dpois(y, mu, log = TRUE): non-integer x = 1.730733## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.636376## Warning in dpois(y, mu, log = TRUE): non-integer x = 3.137160## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.071975## Warning in dpois(y, mu, log = TRUE): non-integer x = 75.146545## Warning in dpois(y, mu, log = TRUE): non-integer x = 1.141354## Warning in dpois(y, mu, log = TRUE): non-integer x = 1.959463## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.866026## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.899245## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.378447## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.253272## Warning in dpois(y, mu, log = TRUE): non-integer x = 1.682292## Warning in dpois(y, mu, log = TRUE): non-integer x = 1.401364## Warning in dpois(y, mu, log = TRUE): non-integer x = 49.814926hist(residuals(M_wildbees2), breaks=20)
plot(simulateResiduals(M_wildbees2),quantreg = T) 
check_collinearity(M_wildbees2)## Warning: Using indexed data frames, such as `df[, 5]`, as model response can produce unexpected results. Specify your model using the literal name of
## the response variable instead.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## lifecycle 1.06 [1.00, 11.03] 1.03 0.94 [0.09, 1.00]
## scale(peak_bloom_timing) 1.22 [1.04, 2.36] 1.11 0.82 [0.42, 0.97]
## scale(groundcover_at_peak_bloom) 1.09 [1.00, 4.27] 1.05 0.92 [0.23, 1.00]
## scale(log(flowercover_at_peak_bloom)) 1.20 [1.03, 2.40] 1.10 0.83 [0.42, 0.97]Anova(M_wildbees2, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: predictions[, 8]
## LR Chisq Df Pr(>Chisq)
## lifecycle 3.1765 1 0.07471 .
## scale(peak_bloom_timing) 0.1333 1 0.71504
## scale(groundcover_at_peak_bloom) 2.0813 1 0.14912
## scale(log(flowercover_at_peak_bloom)) 0.5031 1 0.47815
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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
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)## Warning in dpois(y, mu, log = TRUE): non-integer x = 150.050819## Warning in dpois(y, mu, log = TRUE): non-integer x = 58.590033## Warning in dpois(y, mu, log = TRUE): non-integer x = 34.828775## Warning in dpois(y, mu, log = TRUE): non-integer x = 30.089985## Warning in dpois(y, mu, log = TRUE): non-integer x = 42.104132## Warning in dpois(y, mu, log = TRUE): non-integer x = 68.042479## Warning in dpois(y, mu, log = TRUE): non-integer x = 75.061162## Warning in dpois(y, mu, log = TRUE): non-integer x = 40.209113## Warning in dpois(y, mu, log = TRUE): non-integer x = 54.262799## Warning in dpois(y, mu, log = TRUE): non-integer x = 26.067625## Warning in dpois(y, mu, log = TRUE): non-integer x = 40.627870## Warning in dpois(y, mu, log = TRUE): non-integer x = 36.304345## Warning in dpois(y, mu, log = TRUE): non-integer x = 31.361279## Warning in dpois(y, mu, log = TRUE): non-integer x = 57.225146## Warning in dpois(y, mu, log = TRUE): non-integer x = 179.116701## Warning in dpois(y, mu, log = TRUE): non-integer x = 82.356216## Warning in dpois(y, mu, log = TRUE): non-integer x = 109.152164## Warning in dpois(y, mu, log = TRUE): non-integer x = 66.977117## Warning in dpois(y, mu, log = TRUE): non-integer x = 39.094985## Warning in dpois(y, mu, log = TRUE): non-integer x = 204.162806## Warning in dpois(y, mu, log = TRUE): non-integer x = 84.987082## Warning in dpois(y, mu, log = TRUE): non-integer x = 97.660414## Warning in dpois(y, mu, log = TRUE): non-integer x = 87.442664## Warning in dpois(y, mu, log = TRUE): non-integer x = 63.207544## Warning in dpois(y, mu, log = TRUE): non-integer x = 153.193681## Warning in dpois(y, mu, log = TRUE): non-integer x = 53.564880## Warning in dpois(y, mu, log = TRUE): non-integer x = 90.750837hist(residuals(M_herbivorous_arthropods2), breaks=20)
plot(simulateResiduals(M_herbivorous_arthropods2),quantreg = T) 
check_collinearity(M_herbivorous_arthropods2)## Warning: Using indexed data frames, such as `df[, 5]`, as model response can produce unexpected results. Specify your model using the literal name of
## the response variable instead.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## lifecycle 1.07 [1.00, 8.30] 1.03 0.94 [0.12, 1.00]
## scale(peak_bloom_timing) 1.24 [1.04, 2.34] 1.11 0.81 [0.43, 0.96]
## scale(groundcover_at_peak_bloom) 1.10 [1.00, 3.98] 1.05 0.91 [0.25, 1.00]
## scale(log(flowercover_at_peak_bloom)) 1.21 [1.03, 2.38] 1.10 0.83 [0.42, 0.97]Anova(M_herbivorous_arthropods2, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: predictions[, 7]
## LR Chisq Df Pr(>Chisq)
## lifecycle 0.47582 1 0.4903
## scale(peak_bloom_timing) 1.95776 1 0.1618
## scale(groundcover_at_peak_bloom) 1.37274 1 0.2413
## scale(log(flowercover_at_peak_bloom)) 0.01619 1 0.8988# 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
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)## Warning in dpois(y, mu, log = TRUE): non-integer x = 94.898270## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.234380## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.118495## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.680337## Warning in dpois(y, mu, log = TRUE): non-integer x = 27.403955## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.849220## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.004712## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.128274## Warning in dpois(y, mu, log = TRUE): non-integer x = 31.192957## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.647113## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.107540## Warning in dpois(y, mu, log = TRUE): non-integer x = 55.526716## Warning in dpois(y, mu, log = TRUE): non-integer x = 9.987496## Warning in dpois(y, mu, log = TRUE): non-integer x = 74.847005## Warning in dpois(y, mu, log = TRUE): non-integer x = 35.365396## Warning in dpois(y, mu, log = TRUE): non-integer x = 20.190194## Warning in dpois(y, mu, log = TRUE): non-integer x = 28.073080## Warning in dpois(y, mu, log = TRUE): non-integer x = 6.012383## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.368633## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.358770## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.886680## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.285458## Warning in dpois(y, mu, log = TRUE): non-integer x = 29.818688## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.245507## Warning in dpois(y, mu, log = TRUE): non-integer x = 37.592606## Warning in dpois(y, mu, log = TRUE): non-integer x = 27.491488## Warning in dpois(y, mu, log = TRUE): non-integer x = 33.273271hist(residuals(M_parasitic_wasps2), breaks=20)
plot(simulateResiduals(M_parasitic_wasps2),quantreg = T) 
check_collinearity(M_parasitic_wasps2)## Warning: Using indexed data frames, such as `df[, 5]`, as model response can produce unexpected results. Specify your model using the literal name of
## the response variable instead.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## lifecycle 1.06 [1.00, 11.69] 1.03 0.94 [0.09, 1.00]
## scale(peak_bloom_timing) 1.21 [1.03, 2.38] 1.10 0.82 [0.42, 0.97]
## scale(groundcover_at_peak_bloom) 1.09 [1.00, 4.55] 1.04 0.92 [0.22, 1.00]
## scale(log(flowercover_at_peak_bloom)) 1.20 [1.03, 2.40] 1.10 0.83 [0.42, 0.97]Anova(M_parasitic_wasps2, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: predictions[, 6]
## LR Chisq Df Pr(>Chisq)
## lifecycle 0.0021 1 0.96345
## scale(peak_bloom_timing) 2.8649 1 0.09053 .
## scale(groundcover_at_peak_bloom) 0.8764 1 0.34919
## scale(log(flowercover_at_peak_bloom)) 1.6544 1 0.19836
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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
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)## Warning in dpois(y, mu, log = TRUE): non-integer x = 47.589431## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.144657## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.248339## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.817193## Warning in dpois(y, mu, log = TRUE): non-integer x = 33.436492## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.912436## Warning in dpois(y, mu, log = TRUE): non-integer x = 14.089789## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.917700## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.091007## Warning in dpois(y, mu, log = TRUE): non-integer x = 34.791259## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.268756## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.027293## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.691970## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.168879## Warning in dpois(y, mu, log = TRUE): non-integer x = 23.203491## Warning in dpois(y, mu, log = TRUE): non-integer x = 27.547142## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.338976## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.033175## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.680784## Warning in dpois(y, mu, log = TRUE): non-integer x = 9.131877## Warning in dpois(y, mu, log = TRUE): non-integer x = 31.660473## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.744132## Warning in dpois(y, mu, log = TRUE): non-integer x = 26.340034## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.519068## Warning in dpois(y, mu, log = TRUE): non-integer x = 42.189787## Warning in dpois(y, mu, log = TRUE): non-integer x = 40.144902## Warning in dpois(y, mu, log = TRUE): non-integer x = 35.507587hist(residuals(M_predatory_arthropods2), breaks=20)
plot(simulateResiduals(M_predatory_arthropods2),quantreg = T) 
check_collinearity(M_predatory_arthropods2)## Warning: Using indexed data frames, such as `df[, 5]`, as model response can produce unexpected results. Specify your model using the literal name of
## the response variable instead.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## lifecycle 1.07 [1.00, 8.84] 1.03 0.94 [0.11, 1.00]
## scale(peak_bloom_timing) 1.22 [1.03, 2.37] 1.10 0.82 [0.42, 0.97]
## scale(groundcover_at_peak_bloom) 1.09 [1.00, 4.75] 1.04 0.92 [0.21, 1.00]
## scale(log(flowercover_at_peak_bloom)) 1.20 [1.03, 2.40] 1.10 0.83 [0.42, 0.97]Anova(M_predatory_arthropods2, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: predictions[, 5]
## LR Chisq Df Pr(>Chisq)
## lifecycle 3.02211 1 0.08214 .
## scale(peak_bloom_timing) 1.98947 1 0.15840
## scale(groundcover_at_peak_bloom) 0.66854 1 0.41356
## scale(log(flowercover_at_peak_bloom)) 0.19598 1 0.65798
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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
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)## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.770728## Warning in dpois(y, mu, log = TRUE): non-integer x = 5.708352## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.705039## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.518689## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.974054## Warning in dpois(y, mu, log = TRUE): non-integer x = 7.906214## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.724804## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.210605## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.999395## Warning in dpois(y, mu, log = TRUE): non-integer x = 53.329076## Warning in dpois(y, mu, log = TRUE): non-integer x = 46.362118## Warning in dpois(y, mu, log = TRUE): non-integer x = 20.856095## Warning in dpois(y, mu, log = TRUE): non-integer x = 14.698838## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.874159## Warning in dpois(y, mu, log = TRUE): non-integer x = 6.748630## Warning in dpois(y, mu, log = TRUE): non-integer x = 21.778979## Warning in dpois(y, mu, log = TRUE): non-integer x = 34.966540## Warning in dpois(y, mu, log = TRUE): non-integer x = 23.532834## Warning in dpois(y, mu, log = TRUE): non-integer x = 4.290037## Warning in dpois(y, mu, log = TRUE): non-integer x = 46.475547## Warning in dpois(y, mu, log = TRUE): non-integer x = 14.994856## Warning in dpois(y, mu, log = TRUE): non-integer x = 78.616383## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.153770## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.747787## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.807591## Warning in dpois(y, mu, log = TRUE): non-integer x = 21.937113## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.452387hist(residuals(M_herbivorous_nematodes2), breaks=20)
plot(simulateResiduals(M_herbivorous_nematodes2),quantreg = T) 
check_collinearity(M_herbivorous_nematodes2)## Warning: Using indexed data frames, such as `df[, 5]`, as model response can produce unexpected results. Specify your model using the literal name of
## the response variable instead.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## lifecycle 1.06 [1.00, 10.54] 1.03 0.94 [0.09, 1.00]
## scale(peak_bloom_timing) 1.25 [1.05, 2.33] 1.12 0.80 [0.43, 0.95]
## scale(groundcover_at_peak_bloom) 1.09 [1.00, 4.14] 1.05 0.91 [0.24, 1.00]
## scale(log(flowercover_at_peak_bloom)) 1.22 [1.04, 2.36] 1.11 0.82 [0.42, 0.96]Anova(M_herbivorous_nematodes2, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: predictions[, 2]
## LR Chisq Df Pr(>Chisq)
## lifecycle 0.10626 1 0.74444
## scale(peak_bloom_timing) 1.68117 1 0.19477
## scale(groundcover_at_peak_bloom) 1.96298 1 0.16120
## scale(log(flowercover_at_peak_bloom)) 2.72508 1 0.09878 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)##
## Call:
## glm.nb(formula = predictions[, 2] ~ lifecycle + scale(peak_bloom_timing) +
## scale(groundcover_at_peak_bloom) + scale(log(flowercover_at_peak_bloom)),
## data = predictions, init.theta = 2.822205276, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.00713 0.12485 24.087 <2e-16 ***
## lifecycle1 0.04119 0.12841 0.321 0.7484
## scale(peak_bloom_timing) -0.18281 0.14018 -1.304 0.1922
## scale(groundcover_at_peak_bloom) 0.19077 0.13033 1.464 0.1433
## scale(log(flowercover_at_peak_bloom)) -0.23731 0.13746 -1.726 0.0843 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.8222) family taken to be 1)
##
## Null deviance: 32.215 on 26 degrees of freedom
## Residual deviance: 27.733 on 22 degrees of freedom
## AIC: 219.45
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.822
## Std. Err.: 0.833
##
## 2 x log-likelihood: -207.453# 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
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)## Warning: Using indexed data frames, such as `df[, 5]`, as model response can produce unexpected results. Specify your model using the literal name of
## the response variable instead.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## lifecycle 1.18 [1.02, 2.48] 1.09 0.85 [0.40, 0.98]
## scale(peak_bloom_timing) 1.46 [1.15, 2.48] 1.21 0.68 [0.40, 0.87]
## scale(groundcover_at_peak_bloom) 1.20 [1.03, 2.39] 1.10 0.83 [0.42, 0.97]
## scale(log(flowercover_at_peak_bloom)) 1.32 [1.08, 2.33] 1.15 0.76 [0.43, 0.93]Anova(M_predatory_nematodes2, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: predictions[, 4]
## LR Chisq Df Pr(>Chisq)
## lifecycle 1.50549 1 0.2198
## scale(peak_bloom_timing) 1.79473 1 0.1804
## scale(groundcover_at_peak_bloom) 1.47817 1 0.2241
## scale(log(flowercover_at_peak_bloom)) 0.03276 1 0.8564# 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
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)## Warning in dpois(y, mu, log = TRUE): non-integer x = 248.459068## Warning in dpois(y, mu, log = TRUE): non-integer x = 157.834776## Warning in dpois(y, mu, log = TRUE): non-integer x = 180.502892## Warning in dpois(y, mu, log = TRUE): non-integer x = 91.421630## Warning in dpois(y, mu, log = TRUE): non-integer x = 287.284674## Warning in dpois(y, mu, log = TRUE): non-integer x = 141.648585## Warning in dpois(y, mu, log = TRUE): non-integer x = 169.313512## Warning in dpois(y, mu, log = TRUE): non-integer x = 178.889619## Warning in dpois(y, mu, log = TRUE): non-integer x = 137.713663## Warning in dpois(y, mu, log = TRUE): non-integer x = 382.609497## Warning in dpois(y, mu, log = TRUE): non-integer x = 249.647276## Warning in dpois(y, mu, log = TRUE): non-integer x = 197.158226## Warning in dpois(y, mu, log = TRUE): non-integer x = 170.902628## Warning in dpois(y, mu, log = TRUE): non-integer x = 128.228595## Warning in dpois(y, mu, log = TRUE): non-integer x = 85.107765## Warning in dpois(y, mu, log = TRUE): non-integer x = 93.861977## Warning in dpois(y, mu, log = TRUE): non-integer x = 201.390766## Warning in dpois(y, mu, log = TRUE): non-integer x = 195.661788## Warning in dpois(y, mu, log = TRUE): non-integer x = 154.354587## Warning in dpois(y, mu, log = TRUE): non-integer x = 258.642057## Warning in dpois(y, mu, log = TRUE): non-integer x = 236.455712## Warning in dpois(y, mu, log = TRUE): non-integer x = 182.696889## Warning in dpois(y, mu, log = TRUE): non-integer x = 114.645071## Warning in dpois(y, mu, log = TRUE): non-integer x = 118.903247## Warning in dpois(y, mu, log = TRUE): non-integer x = 152.079699## Warning in dpois(y, mu, log = TRUE): non-integer x = 184.673399## Warning in dpois(y, mu, log = TRUE): non-integer x = 275.170041hist(residuals(M_decomposer_nematodes2), breaks=20)
plot(simulateResiduals(M_decomposer_nematodes2),quantreg = T) 
check_collinearity(M_decomposer_nematodes2)## Warning: Using indexed data frames, such as `df[, 5]`, as model response can produce unexpected results. Specify your model using the literal name of
## the response variable instead.## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## lifecycle 1.06 [1.00, 9.63] 1.03 0.94 [0.10, 1.00]
## scale(peak_bloom_timing) 1.23 [1.04, 2.35] 1.11 0.81 [0.43, 0.96]
## scale(groundcover_at_peak_bloom) 1.09 [1.00, 4.23] 1.05 0.91 [0.24, 1.00]
## scale(log(flowercover_at_peak_bloom)) 1.21 [1.03, 2.38] 1.10 0.83 [0.42, 0.97]Anova(M_decomposer_nematodes2, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: predictions[, 3]
## LR Chisq Df Pr(>Chisq)
## lifecycle 15.8184 1 6.972e-05 ***
## scale(peak_bloom_timing) 0.1270 1 0.7215
## scale(groundcover_at_peak_bloom) 1.8895 1 0.1693
## scale(log(flowercover_at_peak_bloom)) 0.3738 1 0.5409
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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
coefficients_comb <- rbind(coefficients_soilpf,
coefficients_soildec,
coefficients_soilpred,
coefficients_leafpred,
coefficients_leafpara,
coefficients_leafherb,
coefficients_wildbee,
coefficients_hovfly)Figure 2
# 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()## png
## 2#===============================================================================================# # Corrleations
# load required package
library(metan) # 1.18.0## |=========================================================|## | Multi-Environment Trial Analysis (metan) v1.18.0 |## | Author: Tiago Olivoto |## | Type 'citation('metan')' to know how to cite metan |## | Type 'vignette('metan_start')' for a short tutorial |## | Visit 'https://bit.ly/pkgmetan' for a complete tutorial |## |=========================================================|##
## Attaching package: 'metan'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:tidyr':
##
## replace_na## The following object is masked from 'package:dplyr':
##
## recode_factorFigure 3
# 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")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.# export plot
tiff("Figure_3.tiff", res=600, width=230, height=170, units="mm", compression = "lzw")
Figure_3
dev.off()## png
## 2Figure S1
# 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")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.# export plot
tiff("Figure_S1.tiff", res=600, width=230, height=170, units="mm", compression = "lzw")
Figure_S1
dev.off()## png
## 2