In RodolfoPelinson/PredatorIsolationStochasticity: Data Analyses for Pelinson et al 2021B
knitr::opts_chunk$set(echo = TRUE)
devtools::load_all()
This is the community variability analyses presented in the main paper.
These are the packages you will need to run this code:
library(lme4) # Version 1.1-23
library(car) # Version 3.0-7
library(emmeans) # Version 1.4.8
library(vegan) # Version 2.5-7
library(DHARMa) # Version 0.3.3.0
Community Variability
First Survey
Computing observed and expected distances to centroid, and beta-deviation.
beta_deviation_SS1 <- beta_deviation(com_SS1, strata = fish_isolation_SS1, times = 10000,
transform = NULL, dist = "bray", fixedmar="both",
shuffle = "both", method = "quasiswap", seed = 2, group = fish_isolation_SS1)
\
\
Observed Community Variability
Looking at residual plots for observed, expected distances to centroids and deviations.
fit_observed_SS1_G <- lm(beta_deviation_SS1$observed_distances~fish_SS1*isolation_SS1)
resid_observed <- simulateResiduals(fit_observed_SS1_G)
plot(resid_observed)
fit_expected_SS1_G <- lm(beta_deviation_SS1$expected_distances~fish_SS1*isolation_SS1)
resid_expected <- simulateResiduals(fit_expected_SS1_G)
plot(resid_expected)
fit_deviation_SS1_G <- lm(beta_deviation_SS1$deviation_distances~fish_SS1*isolation_SS1)
resid_deviation <- simulateResiduals(fit_deviation_SS1_G)
plot(resid_deviation)
Everything seems ok.
\
\
Running ANOVA table for observed distances to group centroids, or observed beta-diversity/community variability in each treatment.
fit_observed_SS1 <- lm(beta_deviation_SS1$observed_distances~fish_SS1*isolation_SS1)
round(Anova(fit_observed_SS1, test.statistic = "F"),3)
No effect of treatments.
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1))
boxplot(beta_deviation_SS1$observed_distances~isolation_SS1*fish_SS1, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,1), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(interaction(isolation_SS1, fish_SS1))
#levelProportions <- summary(fish_isolation_SS1)/length(beta_deviation_SS1$observed_distances)
col <- c(rep("sienna3",3), rep("dodgerblue3",3))
#bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- beta_deviation_SS1$observed_distances[interaction(isolation_SS1, fish_SS1)==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], cex = 1.5, lwd = 3)
}
boxplot(beta_deviation_SS1$observed_distances~isolation_SS1*fish_SS1, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Community variability", cex.lab = 1.3, line = 3)
title(ylab = "(Distance to Centroid)", cex.lab = 1.3, line = 1.75)
Expected Community Variability
Running ANOVA table for observed distances to group centroids, or observed beta-diversity/community variability in each treatment.
fit_expected_SS1 <- lm(beta_deviation_SS1$expected_distances~fish_SS1*isolation_SS1)
round(Anova(fit_expected_SS1, test.statistic = "F"),3)
emmeans(fit_expected_SS1, list(pairwise ~ isolation_SS1), adjust = "sidak")
There is an increase in expected distance to centroid from intermediate to high isolation.
\
\
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1))
boxplot(beta_deviation_SS1$expected_distances~isolation_SS1*fish_SS1, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,1), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(interaction(isolation_SS1, fish_SS1))
#levelProportions <- summary(fish_isolation_SS1)/length(beta_deviation_SS1$expected_distances)
col <- c(rep("sienna3",3), rep("dodgerblue3",3))
#bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- beta_deviation_SS1$expected_distances[interaction(isolation_SS1, fish_SS1)==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], cex = 1.5, lwd = 3)
}
boxplot(beta_deviation_SS1$expected_distances~isolation_SS1*fish_SS1, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Expected community variability", cex.lab = 1.3, line = 3)
title(ylab = "(Distance to Centroid)", cex.lab = 1.3, line = 1.75)
\
\
Beta-Deviation
Running ANOVA table for observed distances to group centroids, or observed beta-diversity/community variability in each treatment.
fit_deviation_SS1 <- lm(beta_deviation_SS1$deviation_distances~fish_SS1*isolation_SS1)
round(Anova(fit_deviation_SS1, test.statistic = "F"),3)
No significant effect of treatments.
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1))
boxplot(beta_deviation_SS1$deviation_distances~isolation_SS1*fish_SS1, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(interaction(isolation_SS1, fish_SS1))
#levelProportions <- summary(fish_isolation_SS1)/length(beta_deviation_SS1$deviation_distances)
col <- c(rep("sienna3",3), rep("dodgerblue3",3))
#bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- beta_deviation_SS1$deviation_distances[interaction(isolation_SS1, fish_SS1)==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], cex = 1.5, lwd = 3)
}
boxplot(beta_deviation_SS1$deviation_distances~isolation_SS1*fish_SS1, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n")
abline(h = 0, lty = 2, lwd = 2, col = "grey50")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Beta-deviation", cex.lab = 1.3, line = 2)
Whole community for the last two surveys.
First loading data
#data(com_SS2_SS3, All, fish_SS2_SS3, isolation_SS2_SS3, SS_SS2_SS3, ID_SS2_SS3,
# fish_isolation_SS2_SS3)
\
\
Computing observed and expected distances to centroid, and beta-deviation.
beta_deviation_SS2_SS3 <- beta_deviation(com_SS2_SS3, strata = All, times = 10000,
transform = NULL, dist = "bray", fixedmar="both",
shuffle = "both", method = "quasiswap", seed = 2,
group = All)
\
\
Looking at residual plots for observed, expected distances to centroids and deviations.
fit_observed_SS2_SS3_G <- lmer(beta_deviation_SS2_SS3$observed_distances~fish_SS2_SS3*
isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3))
resid_observed <- simulateResiduals(fit_observed_SS2_SS3_G)
plot(resid_observed)
fit_expected_SS2_SS3_G <- lmer(beta_deviation_SS2_SS3$expected_distances~fish_SS2_SS3*
isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3),
control = lmerControl(optimizer = c("bobyqa"), restart_edge = TRUE, boundary.tol = 1e-10))
resid_expected <- simulateResiduals(fit_expected_SS2_SS3_G)
plot(resid_expected)
fit_deviation_SS2_SS3_G <- lmer(beta_deviation_SS2_SS3$deviation_distances~fish_SS2_SS3*
isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3))
resid_deviation <- simulateResiduals(fit_deviation_SS2_SS3_G)
plot(resid_deviation)
Residual plots are not perfect, but they also don't seem too bad.
\
\
Observed Community Variability
Running ANOVA table for observed distances to group centroids, or observed beta-diversity/community variability in each treatment.
fit_observed_SS2_SS3 <- lmer(beta_deviation_SS2_SS3$observed_distances~fish_SS2_SS3*
isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3))
round(Anova(fit_observed_SS2_SS3, test.statistic = "Chisq"),3)
emmeans(fit_observed_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|fish_SS2_SS3),
adjust = "sidak")
emmeans(fit_observed_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|SS_SS2_SS3),
adjust = "sidak")
emmeans(fit_observed_SS2_SS3, list(pairwise ~ SS_SS2_SS3|isolation_SS2_SS3),
adjust = "sidak")
It seems that the effect of isolation is dependent on the presence or absence of fish. When fish is absent, there is no effect of isolation. When it is present, there is a negative effect of isolation.
\
\
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1))
boxplot(beta_deviation_SS2_SS3$observed_distances~isolation_SS2_SS3*fish_SS2_SS3,
outline = F, xlab = "", ylab = "",
at = c(1,2,3,5,6,7),ylim = c(0,1), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(All)
#levelProportions <- summary(All)/length(beta_deviation_SS2_SS3$observed_distances)
col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(22,21,24,22,21,24,15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- beta_deviation_SS2_SS3$observed_distances[All==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(beta_deviation_SS2_SS3$observed_distances~isolation_SS2_SS3*fish_SS2_SS3,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Community variability", cex.lab = 1.3, line = 3)
title(ylab = "(Distance to centroid)", cex.lab = 1.3, line = 1.75)
Expected Community Variability
Running ANOVA table for expected distances to group centroids, or expected beta-diversity/community variability in each treatment.
fit_expected_SS2_SS3 <- lmer(beta_deviation_SS2_SS3$expected_distances~fish_SS2_SS3*isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3), control = lmerControl(optimizer = "nlminbwrap"))
round(Anova(fit_expected_SS2_SS3, test.statistic = "Chisq"),3)
emmeans(fit_expected_SS2_SS3, list(pairwise ~ isolation_SS2_SS3),
adjust = "sidak")
emmeans(fit_expected_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|fish_SS2_SS3),
adjust = "sidak")
emmeans(fit_expected_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|SS_SS2_SS3),
adjust = "sidak")
emmeans(fit_expected_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|fish_SS2_SS3|SS_SS2_SS3),
adjust = "sidak")
Patterns are similar to those observed for the observed distances to centroid.
\
\
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1))
boxplot(beta_deviation_SS2_SS3$expected_distances~isolation_SS2_SS3*fish_SS2_SS3,
outline = F, xlab = "", ylab = "",
at = c(1,2,3,5,6,7),ylim = c(0,1), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(All)
#levelProportions <- summary(All)/length(beta_deviation_SS2_SS3$expected_distances)
col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(22,21,24,22,21,24,15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- beta_deviation_SS2_SS3$expected_distances[All==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(beta_deviation_SS2_SS3$expected_distances~isolation_SS2_SS3*fish_SS2_SS3,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Expected community variability", cex.lab = 1.3, line = 3)
title(ylab = "(Distance to centroid)", cex.lab = 1.3, line = 1.75)
\
\
Beta-Deviation
Running ANOVA table for the deviations of expected distances to group centroids from observed distances.
fit_deviation_SS2_SS3 <- lmer(beta_deviation_SS2_SS3$deviation_distances~fish_SS2_SS3*
isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3),
control = lmerControl(optimizer = "nlminbwrap"))
round(Anova(fit_deviation_SS2_SS3, test.statistic = "Chisq"),3)
emmeans(fit_deviation_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|fish_SS2_SS3),
adjust = "sidak")
Beta deviation seems to increase with isolation, but only in fishless ponds.
\
\
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1))
boxplot(beta_deviation_SS2_SS3$deviation_distances~isolation_SS2_SS3*fish_SS2_SS3,
outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),
ylim = c(-2,10), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(All)
#levelProportions <- summary(All)/length(beta_deviation_SS2_SS3$deviation_distances)
col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(22,21,24,22,21,24,15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- beta_deviation_SS2_SS3$deviation_distances[All==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(beta_deviation_SS2_SS3$deviation_distances~isolation_SS2_SS3*fish_SS2_SS3,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
abline(h = 0, lty = 2, lwd = 2, col = "grey50")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Beta-deviation", cex.lab = 1.3, line = 2)
\
\
\
Observed Environmental Variability
Multivariate effect of Environmental Variability
Getting the distances
env_data_SS2_SS3_st <- decostand(env_data_SS2_SS3, method = "stand", na.rm = TRUE)
env_data_SS2_SS3_dist <- vegdist(env_data_SS2_SS3_st, method = "euclidean", na.rm = TRUE)
betadisper_SS2_SS3 <- betadisper(env_data_SS2_SS3_dist, group = All)
\
\
Checking model fit
Environmental Variability as a function of treatments
fit_env_SS2_SS3 <- lmer(betadisper_SS2_SS3$distances~fish_SS2_SS3*
isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3))
resid_env <- simulateResiduals(fit_env_SS2_SS3)
plot(resid_env)
There seem to be some patterns in the residuals that are not being accounted for by this model. Anyway Even if there is something causing some pattern in community variability, it seem to not be related to our treatments.
\
\
Observed community variability as a function of environmental variability
fit_observed_SS2_SS3_env <- lmer(beta_deviation_SS2_SS3$observed_distances ~ betadisper_SS2_SS3$distances + (1|ID_SS2_SS3))
resid_observed_env <- simulateResiduals(fit_observed_SS2_SS3_env)
plot(resid_observed_env)
everything seems ok.
\
\
Beta deviation a function of environmental variability
fit_deviation_SS2_SS3_env <- lmer(beta_deviation_SS2_SS3$deviation_distances ~ betadisper_SS2_SS3$distances + (1|ID_SS2_SS3))
resid_deviation_env <- simulateResiduals(fit_deviation_SS2_SS3_env)
plot(resid_deviation_env)
everything seems ok.
\
\
running anovas for each of those models
anova_env <- round(Anova(fit_env_SS2_SS3, test.statistic = "Chisq"),3)
anova_env
anova_observed_env <- round(Anova(fit_observed_SS2_SS3_env, test.statistic = "Chisq"),3)
anova_observed_env
anova_deviation_env <- round(Anova(fit_deviation_SS2_SS3_env, test.statistic = "Chisq"),3)
anova_deviation_env
It seems that there is no evidence of any kind of effect of treatments on environmental variabiliy or of environmental variability on community variability or beta deviation
\
\
Ploting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1))
boxplot(betadisper_SS2_SS3$distances~isolation_SS2_SS3*fish_SS2_SS3,
outline = F, ylab = "", xlab = "",
at = c(1,2,3,5,6,7),ylim = c(0,5), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(All)
levelProportions <- summary(All)/length(betadisper_SS2_SS3$distances)
col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(22,21,24,22,21,24,15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- betadisper_SS2_SS3$distances[All==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(betadisper_SS2_SS3$distances~isolation_SS2_SS3*fish_SS2_SS3,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Environmental variability", cex.lab = 1.3, line = 3)
title(ylab = "(Distance to centroid)", cex.lab = 1.3, line = 1.75)
par(mfrow = c(1,2))
plot(beta_deviation_SS2_SS3$observed_distances ~ betadisper_SS2_SS3$distances, pch = 16, ylab = "Distance to centroid (Observed comm. variability)", xlab = "Distance to Centroid (Environmental variability)")
plot(beta_deviation_SS2_SS3$deviation_distances ~ betadisper_SS2_SS3$distances, pch = 16, ylab = "Beta Deviation", xlab = "Distance to centroid (environmental variability)")
abline(h = 0, lty = 2, col = "grey50", lwd = 2)
\
\
Univariate effect of Environmental Variability
First, we have to built a matrix with all univariate distances, that is, how much each observation of each variable differs from its treatment mean.
distances_env_uni <- data.frame(matrix(NA, ncol = ncol(env_data_SS2_SS3_st), nrow = nrow(env_data_SS2_SS3_st)))
for(i in 1:ncol(env_data_SS2_SS3_st)){
if(anyNA(env_data_SS2_SS3_st[,i])){
na_position <- match(NA, env_data_SS2_SS3_st[,i])
distances_env_uni_dist <- vegdist(env_data_SS2_SS3_st[-na_position,i], method = "euclidean", na.rm = FALSE)
betadisper_env_uni <- betadisper(distances_env_uni_dist, group = All[-na_position])
distances <- betadisper_env_uni$distances
distances <- append(distances, NA, after=na_position)
distances_env_uni[,i] <- distances
}else{
distances_env_uni_dist <- vegdist(env_data_SS2_SS3_st[,i], method = "euclidean", na.rm = FALSE)
betadisper_env_uni <- betadisper(distances_env_uni_dist, group = All)
distances_env_uni[,i] <- betadisper_env_uni$distances
}
}
colnames(distances_env_uni) <- colnames(env_data_SS2_SS3_st)
\
\
Now we can run ANOVAS for each environmental variable
uni_anova_env <- data.frame(matrix(NA, ncol = ncol(distances_env_uni), nrow = 7))
for(i in 1:ncol(distances_env_uni)){
fit <- lmer(distances_env_uni[,i]~fish_SS2_SS3*isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3))
uni_anova_env[,i] <- round(Anova(fit, test.statistic = "Chisq"),3)$`Pr(>Chisq)`
}
rownames(uni_anova_env) <- c("fish", "isolation", "survey", "fish:isolation", "fish:survey", "isolation:survey","fish:isolation:survey")
colnames(uni_anova_env) <- colnames(distances_env_uni)
uni_anova_env
\
\
Because we are blindly looking for an effect without previous hypothesis for specific variables, it is important to correct p values of each main and interactive effect for multiple comparisons.
uni_anova_env_adjusted_p <- uni_anova_env
for(i in 1:nrow(uni_anova_env)){
uni_anova_env_adjusted_p[i,] <- p.adjust(uni_anova_env_adjusted_p[i,], method = "fdr")
}
uni_anova_env_adjusted_p <- round(uni_anova_env_adjusted_p, 3)
uni_anova_env_adjusted_p
It seems like there is no important clear effects.
\
\
Plotting it.
par(mfrow = c(4,2))
for(j in 1:ncol(distances_env_uni)){
boxplot(distances_env_uni[,j]~isolation_SS2_SS3*fish_SS2_SS3,
outline = F, ylab = "Variability", xlab = "",
at = c(1,2,3,5,6,7), ylim = c(0,max(na.omit(distances_env_uni[,j]))*1.1), lwd = 1, col = "transparent", xaxt="n", main = paste(colnames(distances_env_uni)[j]))
mylevels <- levels(All)
levelProportions <- summary(All)/length(distances_env_uni[,j])
col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3))
c(22,21,24,22,21,24,15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- distances_env_uni[,j][All==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=levelProportions[i]/0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(distances_env_uni[,j]~isolation_SS2_SS3*fish_SS2_SS3,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1, xaxt="n")
axis(1,labels = c("30 m","120 m", "480 m","30 m","120 m", "480 m"), cex.axis = 0.8,
at =c(1,2,3,5,6,7), gap.axis = 0)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1, at =c(2,6), line = 1.5, tick = F )
}
\
\
Lets also check if any if the variability in any of those variables have an effect on observed community variability...
uni_anova_env_com <- data.frame(matrix(NA, nrow = ncol(distances_env_uni), ncol = 3))
for(i in 1:ncol(distances_env_uni)){
fit <- lmer(beta_deviation_SS2_SS3$observed_distances ~ distances_env_uni[,i] + (1|ID_SS2_SS3))
uni_anova_env_com[i,1] <- fit@beta[2]
uni_anova_env_com[i,2] <- round(Anova(fit, test.statistic = "Chisq"),3)$`Pr(>Chisq)`
}
uni_anova_env_com[,3] <- p.adjust(uni_anova_env_com[,2], method = "fdr")
rownames(uni_anova_env_com) <- colnames(env_data_SS2_SS3_st)
colnames(uni_anova_env_com) <- c("estimate","p","adjust.p")
uni_anova_env_com <- round(uni_anova_env_com, 3)
uni_anova_env_com
and beta deviation
uni_anova_env_dev <- data.frame(matrix(NA, nrow = ncol(distances_env_uni), ncol = 3))
for(i in 1:ncol(distances_env_uni)){
fit <- lmer(beta_deviation_SS2_SS3$deviation_distances ~ distances_env_uni[,i] + (1|ID_SS2_SS3))
uni_anova_env_dev[i,1] <- fit@beta[2]
uni_anova_env_dev[i,2] <- round(Anova(fit, test.statistic = "Chisq"),3)$`Pr(>Chisq)`
}
uni_anova_env_dev[,3] <- p.adjust(uni_anova_env_dev[,2], method = "fdr")
rownames(uni_anova_env_dev) <- colnames(env_data_SS2_SS3_st)
colnames(uni_anova_env_dev) <- c("estimate","p","adjust.p")
uni_anova_env_dev <-round(uni_anova_env_dev, 3)
uni_anova_env_dev
Again. It seems like there is no important clear effects.
\
\
We can plot it.
par(mfrow = c(4,2))
for(i in 1:ncol(distances_env_uni)){
plot(beta_deviation_SS2_SS3$observed_distances ~ distances_env_uni[,i],
pch = 16, ylab = "Observed Comm. variability",
xlab = "Environmental variability", main = paste(colnames(distances_env_uni)[i]))
}
par(mfrow = c(4,2))
for(i in 1:ncol(distances_env_uni)){
plot(beta_deviation_SS2_SS3$deviation_distances ~ distances_env_uni[,i],
pch = 16, ylab = "Beta-deviation",
xlab = "Environmental variability", main = paste(colnames(distances_env_uni)[i]))
}
RodolfoPelinson/PredatorIsolationStochasticity documentation built on Feb. 21, 2022, 12:03 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
devtools::load_all()
This is the community variability analyses presented in the main paper.
These are the packages you will need to run this code:
library(lme4) # Version 1.1-23 library(car) # Version 3.0-7 library(emmeans) # Version 1.4.8 library(vegan) # Version 2.5-7 library(DHARMa) # Version 0.3.3.0
Community Variability
First Survey
Computing observed and expected distances to centroid, and beta-deviation.
beta_deviation_SS1 <- beta_deviation(com_SS1, strata = fish_isolation_SS1, times = 10000, transform = NULL, dist = "bray", fixedmar="both", shuffle = "both", method = "quasiswap", seed = 2, group = fish_isolation_SS1)
\ \
Observed Community Variability
Looking at residual plots for observed, expected distances to centroids and deviations.
fit_observed_SS1_G <- lm(beta_deviation_SS1$observed_distances~fish_SS1*isolation_SS1) resid_observed <- simulateResiduals(fit_observed_SS1_G) plot(resid_observed) fit_expected_SS1_G <- lm(beta_deviation_SS1$expected_distances~fish_SS1*isolation_SS1) resid_expected <- simulateResiduals(fit_expected_SS1_G) plot(resid_expected) fit_deviation_SS1_G <- lm(beta_deviation_SS1$deviation_distances~fish_SS1*isolation_SS1) resid_deviation <- simulateResiduals(fit_deviation_SS1_G) plot(resid_deviation)
Everything seems ok. \ \
Running ANOVA table for observed distances to group centroids, or observed beta-diversity/community variability in each treatment.
fit_observed_SS1 <- lm(beta_deviation_SS1$observed_distances~fish_SS1*isolation_SS1) round(Anova(fit_observed_SS1, test.statistic = "F"),3)
No effect of treatments.
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1)) boxplot(beta_deviation_SS1$observed_distances~isolation_SS1*fish_SS1, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,1), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(interaction(isolation_SS1, fish_SS1)) #levelProportions <- summary(fish_isolation_SS1)/length(beta_deviation_SS1$observed_distances) col <- c(rep("sienna3",3), rep("dodgerblue3",3)) #bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) pch <- c(15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- beta_deviation_SS1$observed_distances[interaction(isolation_SS1, fish_SS1)==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], cex = 1.5, lwd = 3) } boxplot(beta_deviation_SS1$observed_distances~isolation_SS1*fish_SS1, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Community variability", cex.lab = 1.3, line = 3) title(ylab = "(Distance to Centroid)", cex.lab = 1.3, line = 1.75)
Expected Community Variability
Running ANOVA table for observed distances to group centroids, or observed beta-diversity/community variability in each treatment.
fit_expected_SS1 <- lm(beta_deviation_SS1$expected_distances~fish_SS1*isolation_SS1) round(Anova(fit_expected_SS1, test.statistic = "F"),3) emmeans(fit_expected_SS1, list(pairwise ~ isolation_SS1), adjust = "sidak")
There is an increase in expected distance to centroid from intermediate to high isolation. \ \
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1)) boxplot(beta_deviation_SS1$expected_distances~isolation_SS1*fish_SS1, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,1), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(interaction(isolation_SS1, fish_SS1)) #levelProportions <- summary(fish_isolation_SS1)/length(beta_deviation_SS1$expected_distances) col <- c(rep("sienna3",3), rep("dodgerblue3",3)) #bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) pch <- c(15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- beta_deviation_SS1$expected_distances[interaction(isolation_SS1, fish_SS1)==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], cex = 1.5, lwd = 3) } boxplot(beta_deviation_SS1$expected_distances~isolation_SS1*fish_SS1, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Expected community variability", cex.lab = 1.3, line = 3) title(ylab = "(Distance to Centroid)", cex.lab = 1.3, line = 1.75)
\ \
Beta-Deviation
Running ANOVA table for observed distances to group centroids, or observed beta-diversity/community variability in each treatment.
fit_deviation_SS1 <- lm(beta_deviation_SS1$deviation_distances~fish_SS1*isolation_SS1) round(Anova(fit_deviation_SS1, test.statistic = "F"),3)
No significant effect of treatments.
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1)) boxplot(beta_deviation_SS1$deviation_distances~isolation_SS1*fish_SS1, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(interaction(isolation_SS1, fish_SS1)) #levelProportions <- summary(fish_isolation_SS1)/length(beta_deviation_SS1$deviation_distances) col <- c(rep("sienna3",3), rep("dodgerblue3",3)) #bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) pch <- c(15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- beta_deviation_SS1$deviation_distances[interaction(isolation_SS1, fish_SS1)==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], cex = 1.5, lwd = 3) } boxplot(beta_deviation_SS1$deviation_distances~isolation_SS1*fish_SS1, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") abline(h = 0, lty = 2, lwd = 2, col = "grey50") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Beta-deviation", cex.lab = 1.3, line = 2)
Whole community for the last two surveys.
First loading data
#data(com_SS2_SS3, All, fish_SS2_SS3, isolation_SS2_SS3, SS_SS2_SS3, ID_SS2_SS3, # fish_isolation_SS2_SS3)
\ \
Computing observed and expected distances to centroid, and beta-deviation.
beta_deviation_SS2_SS3 <- beta_deviation(com_SS2_SS3, strata = All, times = 10000, transform = NULL, dist = "bray", fixedmar="both", shuffle = "both", method = "quasiswap", seed = 2, group = All)
\ \
Looking at residual plots for observed, expected distances to centroids and deviations.
fit_observed_SS2_SS3_G <- lmer(beta_deviation_SS2_SS3$observed_distances~fish_SS2_SS3* isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3)) resid_observed <- simulateResiduals(fit_observed_SS2_SS3_G) plot(resid_observed) fit_expected_SS2_SS3_G <- lmer(beta_deviation_SS2_SS3$expected_distances~fish_SS2_SS3* isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3), control = lmerControl(optimizer = c("bobyqa"), restart_edge = TRUE, boundary.tol = 1e-10)) resid_expected <- simulateResiduals(fit_expected_SS2_SS3_G) plot(resid_expected) fit_deviation_SS2_SS3_G <- lmer(beta_deviation_SS2_SS3$deviation_distances~fish_SS2_SS3* isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3)) resid_deviation <- simulateResiduals(fit_deviation_SS2_SS3_G) plot(resid_deviation)
Residual plots are not perfect, but they also don't seem too bad. \ \
Observed Community Variability
Running ANOVA table for observed distances to group centroids, or observed beta-diversity/community variability in each treatment.
fit_observed_SS2_SS3 <- lmer(beta_deviation_SS2_SS3$observed_distances~fish_SS2_SS3* isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3)) round(Anova(fit_observed_SS2_SS3, test.statistic = "Chisq"),3) emmeans(fit_observed_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|fish_SS2_SS3), adjust = "sidak") emmeans(fit_observed_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|SS_SS2_SS3), adjust = "sidak") emmeans(fit_observed_SS2_SS3, list(pairwise ~ SS_SS2_SS3|isolation_SS2_SS3), adjust = "sidak")
It seems that the effect of isolation is dependent on the presence or absence of fish. When fish is absent, there is no effect of isolation. When it is present, there is a negative effect of isolation. \ \
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1)) boxplot(beta_deviation_SS2_SS3$observed_distances~isolation_SS2_SS3*fish_SS2_SS3, outline = F, xlab = "", ylab = "", at = c(1,2,3,5,6,7),ylim = c(0,1), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(All) #levelProportions <- summary(All)/length(beta_deviation_SS2_SS3$observed_distances) col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3)) pch <- c(22,21,24,22,21,24,15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- beta_deviation_SS2_SS3$observed_distances[All==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(beta_deviation_SS2_SS3$observed_distances~isolation_SS2_SS3*fish_SS2_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Community variability", cex.lab = 1.3, line = 3) title(ylab = "(Distance to centroid)", cex.lab = 1.3, line = 1.75)
Expected Community Variability
Running ANOVA table for expected distances to group centroids, or expected beta-diversity/community variability in each treatment.
fit_expected_SS2_SS3 <- lmer(beta_deviation_SS2_SS3$expected_distances~fish_SS2_SS3*isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3), control = lmerControl(optimizer = "nlminbwrap")) round(Anova(fit_expected_SS2_SS3, test.statistic = "Chisq"),3) emmeans(fit_expected_SS2_SS3, list(pairwise ~ isolation_SS2_SS3), adjust = "sidak") emmeans(fit_expected_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|fish_SS2_SS3), adjust = "sidak") emmeans(fit_expected_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|SS_SS2_SS3), adjust = "sidak") emmeans(fit_expected_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|fish_SS2_SS3|SS_SS2_SS3), adjust = "sidak")
Patterns are similar to those observed for the observed distances to centroid. \ \
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1)) boxplot(beta_deviation_SS2_SS3$expected_distances~isolation_SS2_SS3*fish_SS2_SS3, outline = F, xlab = "", ylab = "", at = c(1,2,3,5,6,7),ylim = c(0,1), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(All) #levelProportions <- summary(All)/length(beta_deviation_SS2_SS3$expected_distances) col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3)) pch <- c(22,21,24,22,21,24,15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- beta_deviation_SS2_SS3$expected_distances[All==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(beta_deviation_SS2_SS3$expected_distances~isolation_SS2_SS3*fish_SS2_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Expected community variability", cex.lab = 1.3, line = 3) title(ylab = "(Distance to centroid)", cex.lab = 1.3, line = 1.75)
\ \
Beta-Deviation
Running ANOVA table for the deviations of expected distances to group centroids from observed distances.
fit_deviation_SS2_SS3 <- lmer(beta_deviation_SS2_SS3$deviation_distances~fish_SS2_SS3* isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3), control = lmerControl(optimizer = "nlminbwrap")) round(Anova(fit_deviation_SS2_SS3, test.statistic = "Chisq"),3) emmeans(fit_deviation_SS2_SS3, list(pairwise ~ isolation_SS2_SS3|fish_SS2_SS3), adjust = "sidak")
Beta deviation seems to increase with isolation, but only in fishless ponds. \ \
Plotting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1)) boxplot(beta_deviation_SS2_SS3$deviation_distances~isolation_SS2_SS3*fish_SS2_SS3, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7), ylim = c(-2,10), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(All) #levelProportions <- summary(All)/length(beta_deviation_SS2_SS3$deviation_distances) col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3)) pch <- c(22,21,24,22,21,24,15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- beta_deviation_SS2_SS3$deviation_distances[All==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(beta_deviation_SS2_SS3$deviation_distances~isolation_SS2_SS3*fish_SS2_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") abline(h = 0, lty = 2, lwd = 2, col = "grey50") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Beta-deviation", cex.lab = 1.3, line = 2)
\ \ \
Observed Environmental Variability
Multivariate effect of Environmental Variability
Getting the distances
env_data_SS2_SS3_st <- decostand(env_data_SS2_SS3, method = "stand", na.rm = TRUE) env_data_SS2_SS3_dist <- vegdist(env_data_SS2_SS3_st, method = "euclidean", na.rm = TRUE) betadisper_SS2_SS3 <- betadisper(env_data_SS2_SS3_dist, group = All)
\ \
Checking model fit
Environmental Variability as a function of treatments
fit_env_SS2_SS3 <- lmer(betadisper_SS2_SS3$distances~fish_SS2_SS3* isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3)) resid_env <- simulateResiduals(fit_env_SS2_SS3) plot(resid_env)
There seem to be some patterns in the residuals that are not being accounted for by this model. Anyway Even if there is something causing some pattern in community variability, it seem to not be related to our treatments. \ \
Observed community variability as a function of environmental variability
fit_observed_SS2_SS3_env <- lmer(beta_deviation_SS2_SS3$observed_distances ~ betadisper_SS2_SS3$distances + (1|ID_SS2_SS3)) resid_observed_env <- simulateResiduals(fit_observed_SS2_SS3_env) plot(resid_observed_env)
everything seems ok. \ \
Beta deviation a function of environmental variability
fit_deviation_SS2_SS3_env <- lmer(beta_deviation_SS2_SS3$deviation_distances ~ betadisper_SS2_SS3$distances + (1|ID_SS2_SS3)) resid_deviation_env <- simulateResiduals(fit_deviation_SS2_SS3_env) plot(resid_deviation_env)
everything seems ok. \ \
running anovas for each of those models
anova_env <- round(Anova(fit_env_SS2_SS3, test.statistic = "Chisq"),3) anova_env anova_observed_env <- round(Anova(fit_observed_SS2_SS3_env, test.statistic = "Chisq"),3) anova_observed_env anova_deviation_env <- round(Anova(fit_deviation_SS2_SS3_env, test.statistic = "Chisq"),3) anova_deviation_env
It seems that there is no evidence of any kind of effect of treatments on environmental variabiliy or of environmental variability on community variability or beta deviation \ \
Ploting it:
par(cex = 0.75, mar = c(4,4,0.1,0.1)) boxplot(betadisper_SS2_SS3$distances~isolation_SS2_SS3*fish_SS2_SS3, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,5), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(All) levelProportions <- summary(All)/length(betadisper_SS2_SS3$distances) col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3)) pch <- c(22,21,24,22,21,24,15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- betadisper_SS2_SS3$distances[All==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(betadisper_SS2_SS3$distances~isolation_SS2_SS3*fish_SS2_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Environmental variability", cex.lab = 1.3, line = 3) title(ylab = "(Distance to centroid)", cex.lab = 1.3, line = 1.75)
par(mfrow = c(1,2)) plot(beta_deviation_SS2_SS3$observed_distances ~ betadisper_SS2_SS3$distances, pch = 16, ylab = "Distance to centroid (Observed comm. variability)", xlab = "Distance to Centroid (Environmental variability)") plot(beta_deviation_SS2_SS3$deviation_distances ~ betadisper_SS2_SS3$distances, pch = 16, ylab = "Beta Deviation", xlab = "Distance to centroid (environmental variability)") abline(h = 0, lty = 2, col = "grey50", lwd = 2)
\ \
Univariate effect of Environmental Variability
First, we have to built a matrix with all univariate distances, that is, how much each observation of each variable differs from its treatment mean.
distances_env_uni <- data.frame(matrix(NA, ncol = ncol(env_data_SS2_SS3_st), nrow = nrow(env_data_SS2_SS3_st))) for(i in 1:ncol(env_data_SS2_SS3_st)){ if(anyNA(env_data_SS2_SS3_st[,i])){ na_position <- match(NA, env_data_SS2_SS3_st[,i]) distances_env_uni_dist <- vegdist(env_data_SS2_SS3_st[-na_position,i], method = "euclidean", na.rm = FALSE) betadisper_env_uni <- betadisper(distances_env_uni_dist, group = All[-na_position]) distances <- betadisper_env_uni$distances distances <- append(distances, NA, after=na_position) distances_env_uni[,i] <- distances }else{ distances_env_uni_dist <- vegdist(env_data_SS2_SS3_st[,i], method = "euclidean", na.rm = FALSE) betadisper_env_uni <- betadisper(distances_env_uni_dist, group = All) distances_env_uni[,i] <- betadisper_env_uni$distances } } colnames(distances_env_uni) <- colnames(env_data_SS2_SS3_st)
\ \
Now we can run ANOVAS for each environmental variable
uni_anova_env <- data.frame(matrix(NA, ncol = ncol(distances_env_uni), nrow = 7)) for(i in 1:ncol(distances_env_uni)){ fit <- lmer(distances_env_uni[,i]~fish_SS2_SS3*isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3)) uni_anova_env[,i] <- round(Anova(fit, test.statistic = "Chisq"),3)$`Pr(>Chisq)` } rownames(uni_anova_env) <- c("fish", "isolation", "survey", "fish:isolation", "fish:survey", "isolation:survey","fish:isolation:survey") colnames(uni_anova_env) <- colnames(distances_env_uni) uni_anova_env
\ \
Because we are blindly looking for an effect without previous hypothesis for specific variables, it is important to correct p values of each main and interactive effect for multiple comparisons.
uni_anova_env_adjusted_p <- uni_anova_env for(i in 1:nrow(uni_anova_env)){ uni_anova_env_adjusted_p[i,] <- p.adjust(uni_anova_env_adjusted_p[i,], method = "fdr") } uni_anova_env_adjusted_p <- round(uni_anova_env_adjusted_p, 3) uni_anova_env_adjusted_p
It seems like there is no important clear effects. \ \
Plotting it.
par(mfrow = c(4,2)) for(j in 1:ncol(distances_env_uni)){ boxplot(distances_env_uni[,j]~isolation_SS2_SS3*fish_SS2_SS3, outline = F, ylab = "Variability", xlab = "", at = c(1,2,3,5,6,7), ylim = c(0,max(na.omit(distances_env_uni[,j]))*1.1), lwd = 1, col = "transparent", xaxt="n", main = paste(colnames(distances_env_uni)[j])) mylevels <- levels(All) levelProportions <- summary(All)/length(distances_env_uni[,j]) col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3)) c(22,21,24,22,21,24,15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- distances_env_uni[,j][All==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=levelProportions[i]/0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(distances_env_uni[,j]~isolation_SS2_SS3*fish_SS2_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1, xaxt="n") axis(1,labels = c("30 m","120 m", "480 m","30 m","120 m", "480 m"), cex.axis = 0.8, at =c(1,2,3,5,6,7), gap.axis = 0) axis(1,labels = c("Fishless","Fish"), cex.axis = 1, at =c(2,6), line = 1.5, tick = F ) }
\ \
Lets also check if any if the variability in any of those variables have an effect on observed community variability...
uni_anova_env_com <- data.frame(matrix(NA, nrow = ncol(distances_env_uni), ncol = 3)) for(i in 1:ncol(distances_env_uni)){ fit <- lmer(beta_deviation_SS2_SS3$observed_distances ~ distances_env_uni[,i] + (1|ID_SS2_SS3)) uni_anova_env_com[i,1] <- fit@beta[2] uni_anova_env_com[i,2] <- round(Anova(fit, test.statistic = "Chisq"),3)$`Pr(>Chisq)` } uni_anova_env_com[,3] <- p.adjust(uni_anova_env_com[,2], method = "fdr") rownames(uni_anova_env_com) <- colnames(env_data_SS2_SS3_st) colnames(uni_anova_env_com) <- c("estimate","p","adjust.p") uni_anova_env_com <- round(uni_anova_env_com, 3) uni_anova_env_com
and beta deviation
uni_anova_env_dev <- data.frame(matrix(NA, nrow = ncol(distances_env_uni), ncol = 3)) for(i in 1:ncol(distances_env_uni)){ fit <- lmer(beta_deviation_SS2_SS3$deviation_distances ~ distances_env_uni[,i] + (1|ID_SS2_SS3)) uni_anova_env_dev[i,1] <- fit@beta[2] uni_anova_env_dev[i,2] <- round(Anova(fit, test.statistic = "Chisq"),3)$`Pr(>Chisq)` } uni_anova_env_dev[,3] <- p.adjust(uni_anova_env_dev[,2], method = "fdr") rownames(uni_anova_env_dev) <- colnames(env_data_SS2_SS3_st) colnames(uni_anova_env_dev) <- c("estimate","p","adjust.p") uni_anova_env_dev <-round(uni_anova_env_dev, 3) uni_anova_env_dev
Again. It seems like there is no important clear effects. \ \
We can plot it.
par(mfrow = c(4,2)) for(i in 1:ncol(distances_env_uni)){ plot(beta_deviation_SS2_SS3$observed_distances ~ distances_env_uni[,i], pch = 16, ylab = "Observed Comm. variability", xlab = "Environmental variability", main = paste(colnames(distances_env_uni)[i])) } par(mfrow = c(4,2)) for(i in 1:ncol(distances_env_uni)){ plot(beta_deviation_SS2_SS3$deviation_distances ~ distances_env_uni[,i], pch = 16, ylab = "Beta-deviation", xlab = "Environmental variability", main = paste(colnames(distances_env_uni)[i])) }
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.