Week 2 Assignment
Load libraries and data
library(here)
library(tidyverse)
library(ggplot2)
library(car)
library(ggResidpanel)
library(ggeffects)
library(effects)
library(emmeans)
## Load data, manipulate to factors
mesocosmdata <- read_csv(here("Week_02", "Data", "leibold_mesocosm_data_subset.csv"))
mesocosmdata$NutrientLevel <- as.factor(mesocosmdata$NutrientLevel)
mesocosmdata$FoodWeb <- as.factor(mesocosmdata$FoodWeb)Linear model
Create a linear model that tests whether richness is explained by nutrient level and/or food web treatment, and whether the effect of food web treatment differs between nutrient levels, while also accounting for the blocked structure of the experiment.
richness.model = lm(Phyto_Chao1 ~ NutrientLevel #nutrient level is a predictor for richness
+ FoodWeb #foodweb as a predictor
+ FoodWeb:NutrientLevel #effect of foodweb changes between nutrient lvls
+ Block, #account for blocked structure
data = mesocosmdata) Analyze the results
#Report F-tests
Anova(richness.model)## Anova Table (Type II tests)
##
## Response: Phyto_Chao1
## Sum Sq Df F value Pr(>F)
## NutrientLevel 161.91 1 7.2289 0.011295 *
## FoodWeb 313.76 2 7.0042 0.002999 **
## Block 13.93 1 0.6217 0.436212
## NutrientLevel:FoodWeb 46.45 2 1.0368 0.366185
## Residuals 716.74 32
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Plots
resid_panel(richness.model, plots = c('resid', 'qq', 'lev', 'hist')) #residuals plot
plot(ggeffect(richness.model, terms = c('NutrientLevel', 'FoodWeb'))) #effects plot
Interpretation
Why is it important to account for the blocks? Analyze the results of the linear model as we did in Homework 1. What is your interpretation of these results so far?
Accounting for the blocks will help explain some residual variation if there are systematic differences between blocks. If the blocks had a statistically significant effect on species richness, there may be a different factor within the blocks themselves that impact species richness.
The nutrient level, whether high or low, is a predictor for species richness/phytoplankton diversity with a relatively large F-value of 7.23 and a p-value of 0.01. The food web structure is an even stronger predictor of phytoplankton diversity, with a lower p-value of 0.003 and a F-value of 7. The effect of food web treatment does not seem to significantly differ between nutrient levels, and the lack of interaction means that the effect of food web is consistent across nutrient levels. The blocks are also not a strong predictor of nutrient level. The Q-Q plot and histogram generally follow a normal distribution, which indicates a good fit for the model. The residual plot is homoscedastic, so residuals are randomly scattered with no clear pattern. The residual-leverage plot shows us there are no particular samples that have a large effect on the fitted relationship given the lack of outliers.
Estimated marginal means of each combo
Use emmeans to calculate the estimated marginal means of each combination of nutrient level and food web treatment (i.e., H + A, H + G, H + P, L + A, L + G, L + P). Now define contrasts to test (1) whether G is different from A in the H treatment, (2) whether G is different from A in the L treatment, and (3) whether G is different from A when averaging over the L and H treatments. (4) define an interaction contrast that tests whether the difference between G and A is itself different between L and H treatments. I.e., is (G – A) for the high nutrient treatments greater or smaller than (G – A) for the low nutrient treatments?
Emmeans of nutrient level & food web treatment
fw.nutrient.contrast = emmeans(richness.model,
specs = ~ NutrientLevel * FoodWeb)
fw.nutrient.contrast## NutrientLevel FoodWeb emmean SE df lower.CL upper.CL
## H A 10.08 1.93 32 6.14 14.0
## L A 12.52 1.68 32 9.10 15.9
## H G 16.23 1.68 32 12.82 19.6
## L G 19.25 1.94 32 15.29 23.2
## H P 9.77 1.95 32 5.81 13.7
## L P 17.36 2.15 32 12.97 21.7
##
## Confidence level used: 0.95G vs A treatment
#Now define contrasts to test (1) whether G is different from A in the H treatment
GvsAinH = contrast(emmeans(richness.model,
specs = ~ FoodWeb * NutrientLevel, #all combinations of factors
at = list(NutrientLevel = "H")), #fix for H treatment only
method = list("G-A in H" = c(-1,1,0))) #levels are A, G, P
GvsAinH## contrast estimate SE df t.ratio p.value
## G-A in H 6.16 2.56 32 2.409 0.0219#(2) whether G is different from A in the L treatment
GvsAinL = contrast(emmeans(richness.model,
specs = ~ FoodWeb * NutrientLevel, #all combinations of factors
at = list(NutrientLevel = "L")), #fix for L treatment only
method = list("G-A in L" = c(-1,1,0))) #levels are A, G, P
GvsAinL## contrast estimate SE df t.ratio p.value
## G-A in L 6.74 2.56 32 2.633 0.0129# (3) whether G is different from A when averaging over the L and H treatments
# (4) define an interaction contrast that tests whether the difference between G and A is itself different between L and H treatments. I.e., is (G – A) for the high nutrient treatments greater or smaller than (G – A) for the low nutrient treatments?
#levels are AH, GH, PH, AL, GL, PL
GAH = c(-1, 1, 0, 0, 0, 0) #G-A for H
GAL = c(0, 0, 0, -1, 1, 0) #G-A for L
Gaverage = (c(0, 1, 0, 0, 0, 0) + c(0, 0,0,0,1,0))/2 #(GH + GL)/2
Aaverage = (c(1, 0, 0, 0,0,0) + c(0, 0, 0, 1, 0, 0))/2 #(AH + AL)/2
GAcontrast = contrast(emmeans(richness.model,
specs = ~ FoodWeb * NutrientLevel), #all combinations of factors
method = list("G-A in HL" = Gaverage - Aaverage, #average over HL treatements
"G-A in H - G-A in L" = GAH-GAL)) #high nutrient vs low nutrient treatment
GAcontrast## contrast estimate SE df t.ratio p.value
## G-A in HL 6.447 1.81 32 3.565 0.0012
## G-A in H - G-A in L -0.582 3.62 32 -0.161 0.8732Interpretation
We would like to know specifically whether the grazer treatment (G) has greater richness than the algae treament (A), and whether the effect of grazers differs between high and low nutrient levels (we think the effect of grazers on coexistence may be greater at high nutrient loading). The grazer treatment has significantly greater phytoplankton richness than the algae treatment for both the high and the low nutrient levels with an estimate of 6.16 and p-value of 0.0219 for high nutrient levels. For the average over both nutrient levels, there is a highly significant effect that the grazer treatment has on phytoplankton richness with an estimate of 6.45 and p-value of 0.0012. The effect of grazers does not appear to differ between high and low nutrient levels - the interaction is not significant given the 0.8732 p-value.
P vs G treatment
Now, repeat the same set of 4 contrasts, but this time ask whether the P treatment is different from the G treatment. We are interested in these contrasts a priori because we think that adding a predator of the grazers may mean that grazers have a weaker effect on phytoplankton diversity, and that this trophic cascade may be more important under high nutrient supply. How do you interpret the results?
#(1) whether P is different from G in the H treatment
PvsGinH = contrast(emmeans(richness.model,
specs = ~ FoodWeb * NutrientLevel, #all combinations of factors
at = list(NutrientLevel = "H")), #fix for H treatment only
method = list("P-G in H" = c(0,-1,1))) #levels are A, G, P
PvsGinH## contrast estimate SE df t.ratio p.value
## P-G in H -6.46 2.58 32 -2.504 0.0176#(2) whether P is different from G in the L treatment
PvsGinL = contrast(emmeans(richness.model,
specs = ~ FoodWeb * NutrientLevel, #all combinations of factors
at = list(NutrientLevel = "L")), #fix for L treatment only
method = list("P-G in L" = c(0,-1,1))) #levels are A, G, P
PvsGinL## contrast estimate SE df t.ratio p.value
## P-G in L -1.9 2.93 32 -0.647 0.5222# (3) whether P is different from G when averaging over the L and H treatments
# (4) define an interaction contrast that tests whether the difference between G and A is itself different between L and H treatments. I.e., is (P – G) for the high nutrient treatments greater or smaller than (P – G) for the low nutrient treatments?
#levels are AH, GH, PH, AL, GL, PL
PGH = c(0, -1, 1, 0, 0, 0) #P-G for H
PGL = c(0, 0, 0, 0, -1, 1) #P-G for L
Paverage = (c(0, 0, 1, 0, 0, 0) + c(0, 0, 0, 0, 0, 1))/2 #(PH + PL)/2
Gaverage = (c(0, 1, 0, 0,0,0) + c(0, 0, 0, 0, 1, 0))/2 #(GH + GL)/2
PGcontrast = contrast(emmeans(richness.model,
specs = ~ FoodWeb * NutrientLevel), #all combinations of factors
method = list("P-G in HL" = Paverage - Gaverage, #average over HL treatements
"P-G in H - P-G in L" = PGH-PGL)) #high nutrient vs low nutrient
PGcontrast## contrast estimate SE df t.ratio p.value
## P-G in HL -4.18 1.98 32 -2.111 0.0426
## P-G in H - P-G in L -4.56 3.85 32 -1.186 0.2443Interpretation
We think that adding a predator of the grazers may mean that grazers have a weaker effect on phytoplankton diversity, and that this trophic cascade may be more important under high nutrient supply. How do you interpret the results?
Adding a predator of the grazers significantly reduces phytoplankton richness in a high nutrient level treatment with the estimate value of -6.46 and a p-value of 0.0176. However, the effect of predators on grazers at low nutrients is not significant given the 0.522 p-value. When the nutrient levels are averaged over high and low treatments, it appears that the predator reduces phytoplankton richness, but the effect is not as strong on grazers given the 0.0426 p-value and an estimate of -4.18 that is smaller in magnitude to the isolated high nutrient level treatment. The trophic cascade is likely stronger at higher nutrient level treatments, but there is no stastical significance comparing the interaction between high and low nutrient treatments.