Combined Density RA
Frank
2022-05-30
Introduction
These are the codes I used to run basic reciprocal averaging for my combined density data (part a) using the function decorana() from the vegan package (Oksanen et al. 2022).
Test run
First, I want to try out this function by using the same dataset that I ran in PC-ORD version 7.09 (McCune and Mefford, 2016).
Load data
library(readxl)
count_partb <- read_excel("~/thesis-codes/Part b (09-19)/RA (count)/Count Part b.xlsx",
sheet = "Count")Make a matrix from count data
m_count_partb <- as.matrix(count_partb[,4:ncol(count_partb)]) #exclude the first column that contains SU informationRun CA/RA
library(vegan)
ra_count_partb <- decorana(m_count_partb, iweigh=0, iresc=4, ira=1, short=0, before=NULL, after=NULL) #no rare spp downweighting;no rescaling; basic RA; Add columns of Site, Year, Subwatershed
Again, I add columns to the existing dataset by using pipes from the dplyr package (Wickham et al., 2022).
library(dplyr) #load the package
data.scores <- as.data.frame(ra_count_partb$rproj) #this is to extract site scores
data.scores.dplyr <- data.scores %>% mutate(Site = count_partb$Site,Year=count_partb$Year,
Subwatershed=count_partb$Subwatershed) %>% relocate(Site, Year,
Subwatershed)
head(data.scores.dplyr)## Site Year Subwatershed RA1 RA2 RA3 RA4
## 1 B1 2015 B1 0.24777670 0.35397883 -0.2565207 -0.2462255
## 2 B5 2015 B1 -0.04005602 -0.04121631 -0.2019329 -0.4469857
## 3 B7 2015 B7 2.56290863 -1.33742912 -0.1844083 -0.1540035
## 4 B10 2015 B9A 0.67047304 1.39308762 -1.1459170 1.2122555
## 5 REF 2015 B1 -0.65603562 -0.51285260 -0.4728799 0.2561581
## 6 B1 2017 B1 -0.09198814 0.39513611 0.5919194 -0.3436318Plot the graph using ggplot2
library(ggplot2) #load the package
count_partb_graph <- ggplot(data.scores.dplyr, aes(x = RA1, y = RA2)) +
geom_point(size = 4,aes( shape = as.character(Year), colour = Subwatershed))+ # geom_polygon(data=data.scores.dplyr,aes(x = RA1, y = RA2,group=Subwatershed),alpha=0.30); this is to add a little convex polygon to visualize the clusters better. You can try to see how it looks.
theme(axis.text.y = element_text(colour = "black", size = 12, face = "bold"),
axis.text.x = element_text(colour = "black", face = "bold", size = 12),
legend.text = element_text(size = 12, face ="bold", colour ="black"),
legend.position = "right", axis.title.y = element_text(face = "bold", size = 14),
axis.title.x = element_text(face = "bold", size = 14, colour = "black"),
legend.title = element_text(size = 14, colour = "black", face = "bold"),
panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA, size = 1.2),
legend.key=element_blank(),
plot.title = element_text(color = "black", size = 30, face = "bold", hjust = 0.5)) +
labs(
title = "RA graph of count data (2015-2019)") +
theme(axis.title.x = element_text(margin=margin(t=10)), #add margin to x-axis title
axis.title.y = element_text(margin=margin(r=10)))+
labs(x = "RA1", colour = "Subwatershed", y = "RA2", shape = "Year") +
scale_colour_manual(values = c("#b10026", "#feb14c","#fc4d2a")) +
geom_text(aes(label=Site),hjust=-0.15, vjust=1.2)
count_partb_graph
ggsave("Part B Count RA graph.jpeg",count_partb_graph,width=15,height=8)Real data
Load data
Nice, now that I have confirmed that my codes for RA have worked and produced the same results generated by PC-ORD version 7.09 (McCune and Mefford, 2016), I am going to use the same codes for my other dataset.
library(readxl)
density_combined <- read_excel("~/thesis-codes/Part a (87-09)/RA (combined density)/density_combined.xlsx",
sheet = "Final")Make a matrix from abundance data
m_density_combined <- as.matrix(density_combined[,4:ncol(density_combined)]) #exclude the first column that contains SU informationRun CA/RA
library(vegan)
ra_density_combined <- decorana(m_density_combined, iweigh=0, iresc=4, ira=1, short=0, before=NULL, after=NULL) #no rare spp downweighting;no rescaling; basic RA; Add columns of Site, Year, Subwatershed
Again, I add columns to the existing dataset by using pipes from the dplyr package (Wickham et al., 2022).
library(dplyr) #load the package
data.scores <- as.data.frame(ra_density_combined$rproj) #this is to extract site scores
data.scores.dplyr <- data.scores %>% mutate(Site = density_combined$Site,Year=density_combined$Year,
Subwatershed=density_combined$Subwatershed) %>% relocate(Site, Year,
Subwatershed)
head(data.scores.dplyr)## Site Year Subwatershed RA1 RA2 RA3 RA4
## 1 B1 1987-1988 B1 0.9775763 -0.7942114 -0.13570811 0.03498867
## 2 B2 1987-1988 B1 0.7368090 -0.6738703 0.11327041 -0.02727801
## 3 B3 1987-1988 B1 0.8233325 -0.5929189 0.01619114 -0.29177196
## 4 B4 1987-1988 B1 0.6901656 -0.2893593 0.19171759 -0.32415977
## 5 B5 1987-1988 B1 1.0273260 0.4074791 -0.03716040 -0.82147778
## 6 B7 1987-1988 B7 1.0282250 1.3610880 -0.22472660 0.49493718Plot the graph using ggplot2
library(ggplot2) #load the package
combined_density_ra_graph <- ggplot(data.scores.dplyr, aes(x = RA1, y = RA2)) +
geom_point(size = 4,aes( shape = Year, colour = Subwatershed))+ # geom_polygon(data=data.scores.dplyr,aes(x = RA1, y = RA2,group=Subwatershed),alpha=0.30); this is to add a little convex polygon to visualize the clusters better. You can try to see how it looks.
theme(axis.text.y = element_text(colour = "black", size = 12, face = "bold"),
axis.text.x = element_text(colour = "black", face = "bold", size = 12),
legend.text = element_text(size = 12, face ="bold", colour ="black"),
legend.position = "right", axis.title.y = element_text(face = "bold", size = 14),
axis.title.x = element_text(face = "bold", size = 14, colour = "black"),
legend.title = element_text(size = 14, colour = "black", face = "bold"),
panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA, size = 1.2),
legend.key=element_blank(),
plot.title = element_text(color = "black", size = 30, face = "bold", hjust = 0.5)) +
labs(
title = "RA graph of abundance data (1987 and 2009)") +
theme(axis.title.x = element_text(margin=margin(t=10)), #add margin to x-axis title
axis.title.y = element_text(margin=margin(r=10)))+
labs(x = "RA1", colour = "Subwatershed", y = "RA2", shape = "Year") +
scale_colour_manual(values = c("#b10026", "#feb14c","#ffeea1","#fc4d2a")) +
geom_text(aes(label=Site),hjust=-0.15, vjust=1.2)
combined_density_ra_graph
ggsave("Part A Abundance RA graph.jpeg",combined_density_ra_graph,width=15,height=8)References
- Jari Oksanen, Gavin L. Simpson, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, Peter R. Minchin, R.B. O’Hara, Peter Solymos, M. Henry H. Stevens, Eduard Szoecs, Helene Wagner, Matt Barbour, Michael Bedward, Ben Bolker, Daniel Borcard, Gustavo Carvalho, Michael Chirico, Miquel De Caceres, Sebastien Durand, Heloisa Beatriz Antoniazi Evangelista, Rich FitzJohn, Michael Friendly, Brendan Furneaux, Geoffrey Hannigan, Mark O. Hill, Leo Lahti, Dan McGlinn, Marie-Helene Ouellette, Eduardo Ribeiro Cunha, Tyler Smith, Adrian Stier, Cajo J.F. Ter Braak and James Weedon (2022). vegan: Community Ecology Package. R package version 2.6-2. https://CRAN.R-project.org/package=vegan
- Hadley Wickham, Romain François, Lionel Henry and Kirill Müller (2022). dplyr: A Grammar of Data Manipulation. R package version 1.0.8. https://CRAN.R-project.org/package=dplyr
- Hadley Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.
- McCune, B., & Mefford, M. J. (2016). PC-ORD. Multivariate Analysis of Ecological Data. Version 7. MjM Software Design, Gleneden Beach, Oregon, U.S.A.