Plant Community Composition Results
dat_try2019=
read.csv("cover19.csv")%>%
mutate(DM_Cover=as.numeric(DM_Cover)) %>%
unite("code",c(Plot,Subplot), sep="_") %>%
as_tibble() %>%
select(code,Species,DM_Cover) %>%
pivot_wider(names_from = "Species", values_from = "DM_Cover", values_fill=0) %>%
select (-"Oxalis stricta",-"Berberis trifoliolata", -"Vine A") %>%
filter(code !="41_3")
trt_try=
read.csv("trt.csv") %>%
unite("code",c(Plot,Subplot), sep="_") %>%
rename(herb=Herbivory.Tr, fire= Fire.Energy.Tr)
trt2019=
dat_try2019 %>%
inner_join(trt_try) %>%
select(code, herb, fire)
dat2019ord=
read.csv("cover19.csv")%>%
mutate(DM_Cover=as.numeric(DM_Cover)) %>%
unite("code",c(Plot,Subplot), sep="_") %>%
as_tibble() %>%
select(code,Species,DM_Cover) %>%
pivot_wider(names_from = "Species", values_from = "DM_Cover", values_fill=0) %>%
select (-"Oxalis stricta",-"Berberis trifoliolata", -"Vine A") %>%
filter(code !="41_3")%>%
select(-code)1 2019 Community data
Ordination - NMDS (using metaMDS from the vegan package in R)
Assessment of homogeneity of variances among groups (fire treatments and herbivory treatments) using betadisper from vegan package
PERMANOVA to compare community composition - model includes fire treatment, herbivory treatment, and fire*herbivory - use ADONIS function in vegan
To cite package ‘vegan’ in publications use:
Jari Oksanen, F. Guillaume Blanchet, Michael Friendly, Roeland Kindt, Pierre Legendre, Dan McGlinn, Peter R. Minchin, R. B. O’Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens, Eduard Szoecs and Helene Wagner (2020). vegan: Community Ecology Package. R package version 2.5-7. https://CRAN.R-project.org/package=vegan
set.seed(101)
dat19_NMDS=metaMDS(dat2019ord,k=4,trymax =1000)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1301534
## Run 1 stress 0.1300628
## ... New best solution
## ... Procrustes: rmse 0.01661139 max resid 0.1074907
## Run 2 stress 0.1300739
## ... Procrustes: rmse 0.009724374 max resid 0.0568385
## Run 3 stress 0.1343117
## Run 4 stress 0.1305389
## ... Procrustes: rmse 0.01829743 max resid 0.07550609
## Run 5 stress 0.1308982
## Run 6 stress 0.1301583
## ... Procrustes: rmse 0.01335899 max resid 0.08099544
## Run 7 stress 0.1300954
## ... Procrustes: rmse 0.003342052 max resid 0.02742332
## Run 8 stress 0.13008
## ... Procrustes: rmse 0.0137003 max resid 0.0698675
## Run 9 stress 0.1310645
## Run 10 stress 0.1328464
## Run 11 stress 0.13005
## ... New best solution
## ... Procrustes: rmse 0.001849439 max resid 0.01147041
## Run 12 stress 0.1301131
## ... Procrustes: rmse 0.00640712 max resid 0.03489886
## Run 13 stress 0.1312104
## Run 14 stress 0.1300683
## ... Procrustes: rmse 0.002835654 max resid 0.01250772
## Run 15 stress 0.1301051
## ... Procrustes: rmse 0.009172482 max resid 0.06723056
## Run 16 stress 0.1300595
## ... Procrustes: rmse 0.001450472 max resid 0.009475022
## ... Similar to previous best
## Run 17 stress 0.1324438
## Run 18 stress 0.1319723
## Run 19 stress 0.1300665
## ... Procrustes: rmse 0.009345564 max resid 0.0561176
## Run 20 stress 0.130069
## ... Procrustes: rmse 0.01099618 max resid 0.05606138
## *** Solution reacheddat19_NMDS##
## Call:
## metaMDS(comm = dat2019ord, k = 4, trymax = 1000)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(dat2019ord))
## Distance: bray
##
## Dimensions: 4
## Stress: 0.13005
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(dat2019ord))'# First create a data frame of the scores from the individual sites.
# This data frame will contain x and y values for where sites are located.
data_scores <- as.data.frame(scores(dat19_NMDS))
# Now add the extra aquaticSiteType column
data_scores <- cbind(data_scores, trt2019$fire)
colnames(data_scores)[5] <- "FireTreatment"
# Next, we can add the scores for species data
species_scores <- as.data.frame(scores(dat19_NMDS, "species"))
# Add a column equivalent to the row name to create species labels
species_scores$species <- rownames(species_scores)
# Now we can build the plot!
ggplot() +
geom_point(data = data_scores, aes(x = NMDS1, y = NMDS2,
color = FireTreatment), size = 3) +
scale_color_manual(values = c("black", "orangered","darkorchid1"),
name = "Fire Treatment") +
annotate(geom = "label", x = -1, y = 1.25, size = 6,
label = paste("Stress: ", round(dat19_NMDS$stress, digits = 3))) +
theme_minimal() +
theme(legend.position = "bottom",
text = element_text(size = 15))
# First create a data frame of the scores from the individual sites.
# This data frame will contain x and y values for where sites are located.
data_scores <- as.data.frame(scores(dat19_NMDS))
# Now add the extra aquaticSiteType column
data_scores <- cbind(data_scores, trt2019$herb)
colnames(data_scores)[5] <- "HerbTreatment"
# Next, we can add the scores for species data
species_scores <- as.data.frame(scores(dat19_NMDS, "species"))
# Add a column equivalent to the row name to create species labels
species_scores$species <- rownames(species_scores)
# Now we can build the plot!
ggplot() +
geom_point(data = data_scores, aes(x = NMDS1, y = NMDS2,
color = HerbTreatment), size = 3) +
scale_color_manual(values = c("darkolivegreen", "goldenrod2"),
name = "Cage") +
annotate(geom = "label", x = -1, y = 1.25, size = 6,
label = paste("Stress: ", round(dat19_NMDS$stress, digits = 3))) +
theme_minimal() +
theme(legend.position = "bottom",
text = element_text(size = 15))
groupsFire=trt2019$fire
groupsHerb=trt2019$herb
disFire <- vegdist(dat2019ord,method="bray")
modFire <- betadisper(disFire, groupsFire)
anova(modFire)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.64128 0.32064 19.553 3.337e-08 ***
## Residuals 138 2.26296 0.01640
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(modFire)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## High Energy-Control 0.06566638 0.002711641 0.12862111 0.0387513
## Low Energy-Control -0.09896321 -0.160894221 -0.03703219 0.0006628
## Low Energy-High Energy -0.16462958 -0.227584320 -0.10167485 0.0000000ordiplot(dat19_NMDS,type="n")
ordispider(dat19_NMDS,groups=trt2019$fire,col=c("black", "orangered", "darkorchid1"),label=T)
ordiplot(dat19_NMDS,type="text")
ordispider(dat19_NMDS,groups=trt2019$fire,col=c("black", "orangered", "darkorchid1"),label=T)
So you can identify influential species…
boxplot(modFire, xlab="Fire Treatment", notch=TRUE, col=c( "black", "orangered", "darkorchid1"))
groupsHerb=trt2019$herb
disHerb <- vegdist(dat2019ord,method="bray")
modHerb <- betadisper(disHerb, groupsHerb)
anova(modHerb)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.00759 0.0075944 0.3672 0.5455
## Residuals 139 2.87440 0.0206791TukeyHSD(modHerb)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## Control-Cage -0.01467836 -0.06256819 0.03321148 0.5454964ordiplot(dat19_NMDS,type="n")
ordispider(dat19_NMDS,groups=trt2019$herb,col=c("darkolivegreen","goldenrod2"),label=T)
boxplot(modHerb, xlab="Cage", notch=TRUE, col=c( "darkolivegreen","goldenrod"))
adon.results<-adonis(dat2019ord ~ trt2019$fire*trt2019$herb, method="bray",perm=999)
print(adon.results)##
## Call:
## adonis(formula = dat2019ord ~ trt2019$fire * trt2019$herb, permutations = 999, method = "bray")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## trt2019$fire 2 2.890 1.44497 4.4060 0.06040 0.001 ***
## trt2019$herb 1 0.343 0.34332 1.0469 0.00718 0.377
## trt2019$fire:trt2019$herb 2 0.336 0.16784 0.5118 0.00702 0.988
## Residuals 135 44.274 0.32796 0.92540
## Total 140 47.843 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12 2020 Community Data Analysis and Visualization
Ordination - NMDS (using metaMDS from the vegan package in R)
Assessment of homogeneity of variances among groups (fire treatments and herbivory treatments) using betadisper from vegan package
PERMANOVA to compare community composition - model includes fire treatment, herbivory treatment, and fire*herbivory - use ADONIS function in vegan
To cite package ‘vegan’ in publications use:
Jari Oksanen, F. Guillaume Blanchet, Michael Friendly, Roeland Kindt, Pierre Legendre, Dan McGlinn, Peter R. Minchin, R. B. O’Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens, Eduard Szoecs and Helene Wagner (2020). vegan: Community Ecology Package. R package version 2.5-7. https://CRAN.R-project.org/package=vegan
dat_try2020=
read.csv("cover20.csv")%>%
mutate(DM_Cover=as.numeric(DM_Cover)) %>%
unite("code",c(Plot,Subplot), sep="_") %>%
as_tibble() %>%
select(code,Species,DM_Cover) %>%
pivot_wider(names_from = "Species", values_from = "DM_Cover", values_fill=0)
trt_try2020=
read.csv("trt.csv")%>%
unite("code",c(Plot,Subplot), sep="_") %>%
rename(herb=Herbivory.Tr, fire= Fire.Energy.Tr)
trt2020=
dat_try2020 %>%
inner_join(trt_try2020) %>%
select(code, herb, fire)
dat2020ord=
read.csv("cover20.csv")%>%
mutate(DM_Cover=as.numeric(DM_Cover)) %>%
unite("code",c(Plot,Subplot), sep="_") %>%
as_tibble() %>%
select(code,Species,DM_Cover) %>%
pivot_wider(names_from = "Species", values_from = "DM_Cover", values_fill=0) %>%
select(-code)set.seed(106)
dat20_NMDS=metaMDS(dat2020ord,k=4,trymax =1000)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.120954
## Run 1 stress 0.1207188
## ... New best solution
## ... Procrustes: rmse 0.01458404 max resid 0.1211542
## Run 2 stress 0.1208106
## ... Procrustes: rmse 0.005498439 max resid 0.05866822
## Run 3 stress 0.1208457
## ... Procrustes: rmse 0.007596885 max resid 0.07897238
## Run 4 stress 0.1208674
## ... Procrustes: rmse 0.01138681 max resid 0.1129392
## Run 5 stress 0.1206619
## ... New best solution
## ... Procrustes: rmse 0.004852536 max resid 0.05272744
## Run 6 stress 0.1251471
## Run 7 stress 0.1214091
## Run 8 stress 0.1209534
## ... Procrustes: rmse 0.0136288 max resid 0.1195356
## Run 9 stress 0.1207866
## ... Procrustes: rmse 0.007955399 max resid 0.08541724
## Run 10 stress 0.1206621
## ... Procrustes: rmse 0.0008873284 max resid 0.006413112
## ... Similar to previous best
## Run 11 stress 0.122984
## Run 12 stress 0.1214209
## Run 13 stress 0.1206616
## ... New best solution
## ... Procrustes: rmse 0.0005542569 max resid 0.004356755
## ... Similar to previous best
## Run 14 stress 0.120662
## ... Procrustes: rmse 0.0006469272 max resid 0.005365695
## ... Similar to previous best
## Run 15 stress 0.121414
## Run 16 stress 0.1208393
## ... Procrustes: rmse 0.01074132 max resid 0.116876
## Run 17 stress 0.1206616
## ... New best solution
## ... Procrustes: rmse 0.0003450684 max resid 0.002325966
## ... Similar to previous best
## Run 18 stress 0.1207643
## ... Procrustes: rmse 0.004425409 max resid 0.02047469
## Run 19 stress 0.1207885
## ... Procrustes: rmse 0.008054406 max resid 0.08784118
## Run 20 stress 0.1206613
## ... New best solution
## ... Procrustes: rmse 0.000342827 max resid 0.002218599
## ... Similar to previous best
## *** Solution reacheddat20_NMDS##
## Call:
## metaMDS(comm = dat2020ord, k = 4, trymax = 1000)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(dat2020ord))
## Distance: bray
##
## Dimensions: 4
## Stress: 0.1206613
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(dat2020ord))'# First create a data frame of the scores from the individual sites.
# This data frame will contain x and y values for where sites are located.
data_scores <- as.data.frame(scores(dat20_NMDS))
# Now add the extra aquaticSiteType column
data_scores <- cbind(data_scores, trt2020$fire)
colnames(data_scores)[5] <- "FireTreatment"
# Next, we can add the scores for species data
species_scores <- as.data.frame(scores(dat20_NMDS, "species"))
# Add a column equivalent to the row name to create species labels
species_scores$species <- rownames(species_scores)
# Now we can build the plot!
ggplot() +
geom_point(data = data_scores, aes(x = NMDS1, y = NMDS2,
color = FireTreatment), size = 3) +
scale_color_manual(values = c("black", "orangered","darkorchid1"),
name = "Fire Treatment") +
annotate(geom = "label", x = -1.5, y = 1.25, size = 6,
label = paste("Stress: ", round(dat20_NMDS$stress, digits = 3))) +
theme_minimal() +
theme(legend.position = "bottom",
text = element_text(size = 15))
# First create a data frame of the scores from the individual sites.
# This data frame will contain x and y values for where sites are located.
data_scores <- as.data.frame(scores(dat20_NMDS))
# Now add the extra aquaticSiteType column
data_scores <- cbind(data_scores, trt2020$herb)
colnames(data_scores)[5] <- "HerbTreatment"
# Next, we can add the scores for species data
species_scores <- as.data.frame(scores(dat20_NMDS, "species"))
# Add a column equivalent to the row name to create species labels
species_scores$species <- rownames(species_scores)
# Now we can build the plot!
ggplot() +
geom_point(data = data_scores, aes(x = NMDS1, y = NMDS2,
color = HerbTreatment), size = 3) +
scale_color_manual(values = c("darkolivegreen", "goldenrod2"),
name = "Cage") +
annotate(geom = "label", x = -1.5, y = 1.25, size = 6,
label = paste("Stress: ", round(dat20_NMDS$stress, digits = 3))) +
theme_minimal() +
theme(legend.position = "bottom",
text = element_text(size = 15))
groupsFire=trt2020$fire
disFire20 <- vegdist(dat2020ord,method="bray")
modFire20 <- betadisper(disFire20, groupsFire)
anova(modFire20)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.1724 0.086203 2.2872 0.1054
## Residuals 139 5.2388 0.037690TukeyHSD(modFire20)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## High Energy-Control 0.08268997 -0.01220775 0.17758769 0.1011038
## Low Energy-Control 0.06149055 -0.03340718 0.15638827 0.2777385
## Low Energy-High Energy -0.02119942 -0.11508216 0.07268332 0.8542733ordiplot(dat20_NMDS,type="n")
ordispider(dat20_NMDS,groups=trt2020$fire,col=c("black", "orangered", "darkorchid1"),label=T)
ordiplot(dat20_NMDS,type="text")
ordispider(dat20_NMDS,groups=trt2020$fire,col=c("black", "orangered", "darkorchid1"),label=T)
So you can identify influential species…
boxplot(modFire20, xlab="Fire Treatment", notch=TRUE, col=c( "black", "orangered", "darkorchid1"))
groupsHerb=trt2020$herb
disHerb20 <- vegdist(dat2020ord,method="bray")
modHerb20 <- betadisper(disHerb20, groupsHerb)
anova(modHerb20)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.0342 0.034248 1.1908 0.277
## Residuals 140 4.0264 0.028760TukeyHSD(modHerb20)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## Control-Cage 0.03106338 -0.02521519 0.08734196 0.2770375ordiplot(dat20_NMDS,type="n")
ordispider(dat20_NMDS,groups=trt2020$herb,col=c("darkolivegreen","goldenrod2"),label=T)
boxplot(modHerb20, xlab="Cage", notch=TRUE, col=c( "darkolivegreen","goldenrod"))
adon.results<-adonis(dat2020ord ~ trt2020$fire*trt2020$herb, method="bray",perm=999)
print(adon.results)##
## Call:
## adonis(formula = dat2020ord ~ trt2020$fire * trt2020$herb, permutations = 999, method = "bray")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## trt2020$fire 2 2.906 1.45287 5.5214 0.07349 0.001 ***
## trt2020$herb 1 0.424 0.42403 1.6115 0.01072 0.120
## trt2020$fire:trt2020$herb 2 0.421 0.21058 0.8003 0.01065 0.610
## Residuals 136 35.786 0.26313 0.90513
## Total 141 39.537 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1