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'))
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 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 species richness, 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 block is 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.
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 = trt.vs.ctrl ~ NutrientLevel * FoodWeb)
fw.nutrient.contrast## $emmeans
## 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.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## L A - H A 2.438 2.56 32 0.954 0.7647
## H G - H A 6.156 2.56 32 2.409 0.0877
## L G - H A 9.175 2.73 32 3.355 0.0092
## H P - H A -0.304 2.75 32 -0.110 0.9991
## L P - H A 7.280 2.91 32 2.504 0.0714
##
## P value adjustment: dunnettx method for 5 testsplot(fw.nutrient.contrast$contrasts)
##Contrasts
#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" = c(-1,1,0))) #levels are A, G, P
GvsAinH## contrast estimate SE df t.ratio p.value
## G-A 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" = c(-1,1,0))) #levels are A, G, P
GvsAinL## contrast estimate SE df t.ratio p.value
## G-A 6.74 2.56 32 2.633 0.0129#(3)whether G is different from A when averaging over the L and H treatments
GvsAaverage = contrast(emmeans(richness.model,
specs = ~ FoodWeb, #L and H treatments averaged
method = list("G-A" = c(-1,1,0)))) #levels are A, G, P## NOTE: Results may be misleading due to involvement in interactionsGvsAaverage## contrast estimate SE df t.ratio p.value
## A effect -2.905 1.07 32 -2.709 0.0161
## G effect 3.542 1.08 32 3.278 0.0076
## P effect -0.636 1.17 32 -0.544 0.5899
##
## Results are averaged over the levels of: NutrientLevel
## P value adjustment: fdr method for 3 tests# (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?
GvsAinteraction = contrast(emmeans(richness.model,
specs = ~ FoodWeb:NutrientLevel, #interaction
method = list("G-A in H vs G- A in L" = c(-1, 1, 0, 1, -1, 0)))) #levels are AH, GH, PH, AL, GL, PL
GvsAinteraction## contrast estimate SE df t.ratio p.value
## A H effect -4.12 1.76 32 -2.345 0.0508
## G H effect 2.03 1.57 32 1.292 0.2466
## P H effect -4.43 1.77 32 -2.502 0.0508
## A L effect -1.69 1.57 32 -1.073 0.2914
## G L effect 5.05 1.77 32 2.852 0.0453
## P L effect 3.16 1.93 32 1.637 0.1670
##
## P value adjustment: fdr method for 6 testsInterpretation
Marginal means
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" = c(0,-1,1))) #levels are A, G, P
PvsGinH## contrast estimate SE df t.ratio p.value
## P-G -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" = c(0,-1,1))) #levels are A, G, P
PvsGinL## contrast estimate SE df t.ratio p.value
## P-G -1.9 2.93 32 -0.647 0.5222#(3)whether P is different from G when averaging over the L and H treatments
PvsGaverage = contrast(emmeans(richness.model,
specs = ~ FoodWeb, #L and H treatments averaged
method = list("P-G" = c(0,-1,1)))) #levels are A, G, P## NOTE: Results may be misleading due to involvement in interactionsPvsGaverage## contrast estimate SE df t.ratio p.value
## A effect -2.905 1.07 32 -2.709 0.0161
## G effect 3.542 1.08 32 3.278 0.0076
## P effect -0.636 1.17 32 -0.544 0.5899
##
## Results are averaged over the levels of: NutrientLevel
## P value adjustment: fdr method for 3 tests# (4) define an interaction contrast that tests whether the difference between P and G 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?
PvsGinteraction = contrast(emmeans(richness.model,
specs = ~ FoodWeb:NutrientLevel, #interaction
method = list("P-G in H vs P-G in L" = c(0, -1, 1, 0, -1, 1)))) #levels are AH, GH, PH, AL, GL, PL
PvsGinteraction## contrast estimate SE df t.ratio p.value
## A H effect -4.12 1.76 32 -2.345 0.0508
## G H effect 2.03 1.57 32 1.292 0.2466
## P H effect -4.43 1.77 32 -2.502 0.0508
## A L effect -1.69 1.57 32 -1.073 0.2914
## G L effect 5.05 1.77 32 2.852 0.0453
## P L effect 3.16 1.93 32 1.637 0.1670
##
## P value adjustment: fdr method for 6 tests