plot code
# libraries
library(ggplot2)
library(readr)
library(gridExtra)
library(ggpubr)
library(patchwork)
library(cowplot)
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:patchwork':
##
## align_plots
## The following object is masked from 'package:ggpubr':
##
## get_legend
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:gridExtra':
##
## combine
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(Rmisc)
## Loading required package: lattice
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following object is masked from 'package:ggpubr':
##
## mutate
# load data
nutrient_data2 <- read_csv("nutrient_data2.csv")
## New names:
## • `` -> `...1`
## Rows: 84 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): ...1, site
## dbl (3): group, Nitrite, Orthophosphate
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
water_quality_data3 <- read.csv("water_quality_data3.csv")
# nitrite nutrient data plot
#summary stats
nitritesumstat <- summarySE(data = nutrient_data2, measurevar = "Nitrite", groupvars = "site", na.rm = TRUE)
nitritesumstat
## site N Nitrite sd se ci
## 1 St_Johns_downstream 28 0.007142857 0.007069571 0.0013360234 0.002741294
## 2 St_Johns_upstream 27 0.007925926 0.008466795 0.0016294354 0.003349352
## 3 Townhouse_pond 27 0.003925926 0.003862302 0.0007433004 0.001527876
# making the plot
plotA <- ggplot(nitritesumstat, aes(x = site, y = Nitrite)) + geom_col(aes(fill = site)) +
#error bars
geom_errorbar(aes(ymin = Nitrite - se, ymax = Nitrite + se), width = 0.5, size = 1) +
# labels
labs(title = "A",
x = "Collection Site",
y = "Nitrite (mg/L NO2-N") +
#theme
theme_classic() +
theme(plot.title = element_text(size = 14),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1))
plotA
# orthophosphate plot
Orthophosphatesumstat <- summarySE(data = nutrient_data2, measurevar = "Orthophosphate", groupvars = "site", na.rm = TRUE)
Orthophosphatesumstat
## site N Orthophosphate sd se ci
## 1 St_Johns_downstream 26 0.1876538 0.3751585 0.07357464 0.15152981
## 2 St_Johns_upstream 27 0.1486667 0.1559349 0.03000969 0.06168579
## 3 Townhouse_pond 26 0.2655385 0.2454073 0.04812833 0.09912215
plotB <- ggplot(Orthophosphatesumstat, aes(x = site, y = Orthophosphate)) + geom_col(aes(fill = site)) +
#error bars
geom_errorbar(aes(ymin = Orthophosphate - se, ymax = Orthophosphate + se), width = 0.5, size = 1) +
# labels
labs(title = "B",
x = "Collection Site",
y = "Orthophosphate (mg/L PO4^3-)") +
#theme
theme_classic() +
theme(plot.title = element_text(size = 14),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1))
plot(plotB)
# multipanel plot
multipanel1 <- plot_grid(
plotA + theme(plot.title = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank()),
plotB + theme(plot.title = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank()),
ncol = 2,
Labels = c("A", "B"),
Label_size = 12,
legend.position = "right")
multipanel1
# different multipanel
multi2 <- grid.arrange(plotA, plotB,
ncol = 2)
# question 3
# tempature salinity do and tubridity
tempstat <- summarySE(data = water_quality_data3, measurevar = "Temperature", groupvars = "Site", na.rm = TRUE)
tempstat
## Site N Temperature sd se ci
## 1 St_Johns_downstream 69 10.31304 2.296859 0.2765094 0.5517659
## 2 St_Johns_upstream 69 11.53333 2.061458 0.2481703 0.4952162
## 3 Townhouse_pond 77 9.97013 1.976777 0.2252747 0.4486734
salinitystat <- summarySE(data = water_quality_data3, measurevar = "Salinity_refractometer", groupvars = "Site", na.rm = TRUE)
salinitystat
## Site N Salinity_refractometer sd se ci
## 1 St_Johns_downstream 69 4.050949 4.953582 0.5963412 1.1899803
## 2 St_Johns_upstream 69 1.000674 1.254104 0.1509764 0.3012687
## 3 Townhouse_pond 74 1.627182 2.703691 0.3142978 0.6263945
DOstat <- summarySE(data = water_quality_data3, measurevar = "DO_percent", groupvars = "Site", na.rm = TRUE)
DOstat
## Site N DO_percent sd se ci
## 1 St_Johns_downstream 69 86.91536 128.49013 15.468394 30.866699
## 2 St_Johns_upstream 69 66.80725 27.85591 3.353457 6.691719
## 3 Townhouse_pond 77 82.91273 31.20133 3.555722 7.081834
turbidstat <- summarySE(data = water_quality_data3, measurevar = "Secchi_depth_cm", groupvars = "Site", na.rm = TRUE)
turbidstat
## Site N Secchi_depth_cm sd se ci
## 1 St_Johns_downstream 16 32.65625 7.453676 1.863419 3.971784
## 2 St_Johns_upstream 26 12.80769 11.271270 2.210478 4.552565
## 3 Townhouse_pond 54 46.62759 22.990849 3.128658 6.275294
# temp graph:
plottemp <- ggplot(tempstat, aes(x = Site, y = Temperature)) + geom_col(aes(fill = Site)) +
#error bars
geom_errorbar(aes(ymin = Temperature - se, ymax = Temperature + se), width = 0.5, size = 1) +
# labels
labs(title = "A",
x = "Collection Site",
y = "Temperature C") +
#theme
theme_classic() +
theme(plot.title = element_text(size = 14),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1))
plottemp
# salinity graph
salinityplot <- ggplot(salinitystat, aes(x = Site, y = Salinity_refractometer)) + geom_col(aes(fill = Site)) +
#error bars
geom_errorbar(aes(ymin = Salinity_refractometer - se, ymax = Salinity_refractometer + se), width = 0.5, size = 1) +
# labels
labs(title = "B",
x = "Collection Site",
y = "salinity") +
#theme
theme_classic() +
theme(plot.title = element_text(size = 14),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1))
salinityplot
# DO plot
DOplot <- ggplot(DOstat, aes(x = Site, y = DO_percent)) + geom_col(aes(fill = Site)) +
#error bars
geom_errorbar(aes(ymin = DO_percent - se, ymax = DO_percent + se), width = 0.5, size = 1) +
# labels
labs(title = "C",
x = "Collection Site",
y = "DO%") +
#theme
theme_classic() +
theme(plot.title = element_text(size = 14),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1))
DOplot
#turbidity
turbplot <- ggplot(turbidstat, aes(x = Site, y = Secchi_depth_cm)) + geom_col(aes(fill = Site)) +
#error bars
geom_errorbar(aes(ymin = Secchi_depth_cm - se, ymax = Secchi_depth_cm + se), width = 0.5, size = 1) +
# labels
labs(title = "D",
x = "Collection Site",
y = "Turbidity") +
#theme
theme_classic() +
theme(plot.title = element_text(size = 14),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1))
turbplot
#2nd multipanel
multi3 <- grid.arrange(plottemp, salinityplot, DOplot, turbplot,
ncol = 2)
ggsave("multi3_plot.png", plot = multi3, width = 12, height = 8, dpi = 300)
# statistical tests
# ANOVA for Nitrite
anova_nitrite <- aov(Nitrite ~ site, data = nutrient_data2)
summary(anova_nitrite)
## Df Sum Sq Mean Sq F value Pr(>F)
## site 2 0.000243 1.217e-04 2.669 0.0756 .
## Residuals 79 0.003601 4.558e-05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2 observations deleted due to missingness
# ANOVA for Orthophosphate
anova_orthophosphate <- aov(Orthophosphate ~ site, data = nutrient_data2)
summary(anova_orthophosphate)
## Df Sum Sq Mean Sq F value Pr(>F)
## site 2 0.187 0.09339 1.255 0.291
## Residuals 76 5.656 0.07443
## 5 observations deleted due to missingness
# Post-hoc test: Tukey's HSD
TukeyHSD(anova_nitrite)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Nitrite ~ site, data = nutrient_data2)
##
## $site
## diff lwr upr
## St_Johns_upstream-St_Johns_downstream 0.0007830688 -0.003566873 0.005133011
## Townhouse_pond-St_Johns_downstream -0.0032169312 -0.007566873 0.001133011
## Townhouse_pond-St_Johns_upstream -0.0040000000 -0.008389309 0.000389309
## p adj
## St_Johns_upstream-St_Johns_downstream 0.9032362
## Townhouse_pond-St_Johns_downstream 0.1875384
## Townhouse_pond-St_Johns_upstream 0.0814717
TukeyHSD(anova_orthophosphate)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Orthophosphate ~ site, data = nutrient_data2)
##
## $site
## diff lwr upr
## St_Johns_upstream-St_Johns_downstream -0.03898718 -0.2181792 0.1402048
## Townhouse_pond-St_Johns_downstream 0.07788462 -0.1029900 0.2587592
## Townhouse_pond-St_Johns_upstream 0.11687179 -0.0623202 0.2960638
## p adj
## St_Johns_upstream-St_Johns_downstream 0.8617766
## Townhouse_pond-St_Johns_downstream 0.5607795
## Townhouse_pond-St_Johns_upstream 0.2696314