OCN 683 - Homework 2
Sarah Pennington
Here we will analyze experimental data from Leibold et al. (2017),
in which they used mesocosm manipulations to test whether grazers
promote phytoplankton diversity in ponds. They factorially implemented a
number of experimental treatments: N:P ratio, light, nutrient level,
presence of grazers, and presence of a predator of the grazers. For the
purposes of this assignment we will only use the balanced N:P treatment,
and we will include both light levels but will not analyze the effect of
light (the authors found it has no effect).
The column NutrientLevel has two levels (L = low and H = high), and
the column FoodWeb has three levels (A = algae, G = algae + grazers, P =
algae + grazers + predator). The treatments were randomly assigned to
mesocosms within four blocks (column Block). The response variable we
will analyze is the column Phyto_Chao1, which is the Chao-1 estimator of
taxonomic richness.
getwd()
## [1] "C:/Users/sarah/Desktop/Masters program/Classes_spring_25"
div <- read_csv("leibold_mesocosm_data_subset.csv",
col_types = cols(NtoPratio = col_number(),
NutrientLevel = col_factor(levels = c("H",
"L")), ShadeCloth = col_factor(levels = c("N",
"Y")), FoodWeb = col_factor(levels = c("A",
"G", "P")), Zoop_Chao1 = col_number(),
`Phyto#taxa` = col_number(), Phyto_Chao1 = col_number(),
Block = col_number()))
names(div) #names of columns #Phyto_Chao1 is the response variable
## [1] "NtoPratio" "NutrientLevel" "ShadeCloth" "FoodWeb"
## [5] "Tank" "Zoop_Chao1" "Phyto#taxa" "Phyto_Chao1"
## [9] "Block"
str(div)
## spc_tbl_ [39 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ NtoPratio : num [1:39] 14 14 14 14 14 14 14 14 14 14 ...
## $ NutrientLevel: Factor w/ 2 levels "H","L": 1 1 1 1 1 1 1 1 1 1 ...
## $ ShadeCloth : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 2 ...
## $ FoodWeb : Factor w/ 3 levels "A","G","P": 1 1 2 2 2 2 3 3 3 1 ...
## $ Tank : num [1:39] 15 115 1 66 84 89 36 75 86 11 ...
## $ Zoop_Chao1 : num [1:39] 0 0 9.33 9.5 5 ...
## $ Phyto#taxa : num [1:39] 9 8 16 13 8 10 12 12 5 17 ...
## $ Phyto_Chao1 : num [1:39] 9 11 19 19 8 ...
## $ Block : num [1:39] 1 4 1 2 3 4 1 2 3 1 ...
## - attr(*, "spec")=
## .. cols(
## .. NtoPratio = col_number(),
## .. NutrientLevel = col_factor(levels = c("H", "L"), ordered = FALSE, include_na = FALSE),
## .. ShadeCloth = col_factor(levels = c("N", "Y"), ordered = FALSE, include_na = FALSE),
## .. FoodWeb = col_factor(levels = c("A", "G", "P"), ordered = FALSE, include_na = FALSE),
## .. Tank = col_double(),
## .. Zoop_Chao1 = col_number(),
## .. `Phyto#taxa` = col_number(),
## .. Phyto_Chao1 = col_number(),
## .. Block = col_number()
## .. )
## - attr(*, "problems")=<externalptr>
div$Block <-as.factor(div$Block)
str(div)
## spc_tbl_ [39 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ NtoPratio : num [1:39] 14 14 14 14 14 14 14 14 14 14 ...
## $ NutrientLevel: Factor w/ 2 levels "H","L": 1 1 1 1 1 1 1 1 1 1 ...
## $ ShadeCloth : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 2 ...
## $ FoodWeb : Factor w/ 3 levels "A","G","P": 1 1 2 2 2 2 3 3 3 1 ...
## $ Tank : num [1:39] 15 115 1 66 84 89 36 75 86 11 ...
## $ Zoop_Chao1 : num [1:39] 0 0 9.33 9.5 5 ...
## $ Phyto#taxa : num [1:39] 9 8 16 13 8 10 12 12 5 17 ...
## $ Phyto_Chao1 : num [1:39] 9 11 19 19 8 ...
## $ Block : Factor w/ 4 levels "1","2","3","4": 1 4 1 2 3 4 1 2 3 1 ...
## - attr(*, "spec")=
## .. cols(
## .. NtoPratio = col_number(),
## .. NutrientLevel = col_factor(levels = c("H", "L"), ordered = FALSE, include_na = FALSE),
## .. ShadeCloth = col_factor(levels = c("N", "Y"), ordered = FALSE, include_na = FALSE),
## .. FoodWeb = col_factor(levels = c("A", "G", "P"), ordered = FALSE, include_na = FALSE),
## .. Tank = col_double(),
## .. Zoop_Chao1 = col_number(),
## .. `Phyto#taxa` = col_number(),
## .. Phyto_Chao1 = col_number(),
## .. Block = col_number()
## .. )
## - attr(*, "problems")=<externalptr>
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. 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?
#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.
div.model = lm(Phyto_Chao1 ~ NutrientLevel * FoodWeb + Block, data = div)
summary(div.model)
##
## Call:
## lm(formula = Phyto_Chao1 ~ NutrientLevel * FoodWeb + Block, data = div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.9575 -3.1031 0.1463 2.2300 9.2563
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.7437 2.3171 4.637 6.48e-05 ***
## NutrientLevelL 2.3912 2.6524 0.902 0.3745
## FoodWebG 6.1100 2.6524 2.304 0.0284 *
## FoodWebP -0.3990 2.9133 -0.137 0.8920
## Block2 -0.1676 2.1969 -0.076 0.9397
## Block3 -0.9499 2.1048 -0.451 0.6550
## Block4 -1.6725 2.3149 -0.722 0.4756
## NutrientLevelL:FoodWebG 0.7125 3.7948 0.188 0.8523
## NutrientLevelL:FoodWebP 5.1877 3.9864 1.301 0.2030
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.885 on 30 degrees of freedom
## Multiple R-squared: 0.4046, Adjusted R-squared: 0.2458
## F-statistic: 2.548 on 8 and 30 DF, p-value: 0.03008
#Why is it important to account for the blocks?
###If block had a significant effect on richness(Phto_Chao1), then that would mean some other factor in individual blocks that was impacting our results.
#Analyze the results of the linear model as we did in Homework 1.
library(ggResidpanel)
resid_panel(div.model, plots = c('resid', 'qq', 'lev', 'hist'))
#Report F tests for terms of the model.
Anova(div.model)
## Anova Table (Type II tests)
##
## Response: Phyto_Chao1
## Sum Sq Df F value Pr(>F)
## NutrientLevel 163.07 1 6.8341 0.013853 *
## FoodWeb 314.91 2 6.5985 0.004217 **
## Block 14.80 3 0.2068 0.890886
## NutrientLevel:FoodWeb 45.82 2 0.9601 0.394315
## Residuals 715.86 30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##F values
##NutrientLevel 7.2289
##FoodWeb 7.0042
##NutrientLevel:FoodWeb 1.0368
#Create an effects plot displaying fitted effects
library(ggeffects)
plot(ggeffect(div.model, terms = c('FoodWeb', 'NutrientLevel')))
#What is your interpretation of the results so far?
###Richness is lower for higher nutrient levels across foodwebs. Foodwebs A and G appear to have the largest difference in richness based on the effect sizes model.
Now we will use contrasts to test some hypotheses, and we will
imagine that we formulated these hypotheses a priori (as opposed to
formulating them after we looked at the data). 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).
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 totest (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.
#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).
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
trt.control.contrast = emmeans(div.model, specs = trt.vs.ctrl ~ NutrientLevel * FoodWeb)
summary(trt.control.contrast)
## $emmeans
## NutrientLevel FoodWeb emmean SE df lower.CL upper.CL
## H A 10.05 2.01 30 5.94 14.2
## L A 12.44 1.73 30 8.91 16.0
## H G 16.16 1.73 30 12.63 19.7
## L G 19.26 2.05 30 15.07 23.5
## H P 9.65 2.06 30 5.44 13.9
## L P 17.23 2.26 30 12.61 21.8
##
## Results are averaged over the levels of: Block
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## L A - H A 2.391 2.65 30 0.902 0.7930
## H G - H A 6.110 2.65 30 2.304 0.1106
## L G - H A 9.214 2.85 30 3.237 0.0130
## H P - H A -0.399 2.91 30 -0.137 0.9985
## L P - H A 7.180 3.05 30 2.353 0.1000
##
## Results are averaged over the levels of: Block
## P value adjustment: dunnettx method for 5 tests
plot(trt.control.contrast$contrast)
#(1) whether G is different from A in the H treatment
levels(div$FoodWeb)
## [1] "A" "G" "P"
###[1] "A" "G" "P"
N.colimitation = contrast(emmeans(div.model, ~ FoodWeb * NutrientLevel, at = list(NutrientLevel = "H")), method
= list("G - A in H" = c(-1, 1, 0)))
print(N.colimitation)
## contrast estimate SE df t.ratio p.value
## G - A in H 6.11 2.65 30 2.304 0.0284
##
## Results are averaged over the levels of: Block
#(2) whether G is different from A in the L treatment
N.colimitation1 = contrast(emmeans(div.model, ~ FoodWeb * NutrientLevel, at = list(NutrientLevel = "L")), method
= list("G - A in L" = c(-1, 1, 0)))
print(N.colimitation1)
## contrast estimate SE df t.ratio p.value
## G - A in L 6.82 2.68 30 2.543 0.0164
##
## Results are averaged over the levels of: Block
#(3) whether G is different from A when averaging over the L and H treatments.
N.colimitation2 = contrast(emmeans(div.model, ~ FoodWeb), method
= list("G - A in H + L" = c(-1, 1, 0)))
## NOTE: Results may be misleading due to involvement in interactions
print(N.colimitation2)
## contrast estimate SE df t.ratio p.value
## G - A in H + L 6.47 1.88 30 3.448 0.0017
##
## Results are averaged over the levels of: NutrientLevel, Block
(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? Hint: for the interaction contrast you
are essentially taking contrast #2 and subtracting it from contrast #1,
in order to test whether the two contrasts are the same, which means
that the difference between them is zero.How do you interpret these
results?
interaction <- contrast(emmeans(div.model, ~ FoodWeb * NutrientLevel), method = list("G - A in H vs G - A in L" = c(-1, 1, 0, 1, -1, 0)))
print(interaction)
## contrast estimate SE df t.ratio p.value
## G - A in H vs G - A in L -0.712 3.79 30 -0.188 0.8523
##
## Results are averaged over the levels of: Block
#How do you interpret these results?
### Going all the way back to the first two contrasts, G-A in H or L, showed richness was significantly different between G and A when Nutrient levels were high and when they were low. The combined contrast with between G-A in H and L nutrient levels showed richness was also significantly different. Lastly, there was no significant difference in G and A difference between L and H treatments. The difference wasn't different between nutrient levels. This could imply perhaps nutrient level is the driver of the difference seen in the first contrasts (?)
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 G is different from P in the H treatment
levels(div$FoodWeb)
## [1] "A" "G" "P"
###[1] "A" "G" "P"
N.colimitation = contrast(emmeans(div.model, ~ FoodWeb * NutrientLevel, at = list(NutrientLevel = "H")), method
= list("G - P in H" = c(0,-1, 1)))
print(N.colimitation)
## contrast estimate SE df t.ratio p.value
## G - P in H -6.51 2.69 30 -2.421 0.0217
##
## Results are averaged over the levels of: Block
#(2) whether G is different from P in the L treatment
N.colimitation1 = contrast(emmeans(div.model, ~ FoodWeb * NutrientLevel, at = list(NutrientLevel = "L")), method
= list("G - P in L" = c(0,-1, 1)))
print(N.colimitation1)
## contrast estimate SE df t.ratio p.value
## G - P in L -2.03 3.11 30 -0.654 0.5182
##
## Results are averaged over the levels of: Block
#(3) whether G is different from P when averaging over the L and H treatments.
N.colimitation2 = contrast(emmeans(div.model, ~ FoodWeb), method
= list("G - P in H + L" = c(0,-1, 1)))
## NOTE: Results may be misleading due to involvement in interactions
print(N.colimitation2)
## contrast estimate SE df t.ratio p.value
## G - P in H + L -4.27 2.1 30 -2.033 0.0510
##
## Results are averaged over the levels of: NutrientLevel, Block
# (4) define an interaction contrast that tests whether the difference between G and P is itself different between L and H treatments.
interaction <- contrast(emmeans(div.model, ~ FoodWeb * NutrientLevel), method = list("G - A in H vs G - A in L" = c(0, -1, 1, 0, 1, -1)))
print(interaction)
## contrast estimate SE df t.ratio p.value
## G - A in H vs G - A in L -4.48 4.02 30 -1.114 0.2742
##
## Results are averaged over the levels of: Block
#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?
###Based on the results of these contrasts, we can see that food web does have an impact on richness. However, it appears the nutrient supply is a large driver of richness differnces also.The researchers probably shouldn't draw the conclusion that trophic cascade is more important than nutrient supply in determining species richness.