Kona Dilution Quaterly Analysis
Sofia Ferreira
Data
library(readr)
library(ggplot2)
library(ggpubr)
library(dplyr)
library(dunn.test)
library(car)
library(tidyr)
library(infer)
# UH
ALL_UH <- read_csv("ALL_UH.csv")
ALL_UH$month <- factor(ALL_UH$month, levels = c("June", "September", "December", "March"))
ALL_UH$Station <- factor(ALL_UH$Station, levels = c("NS", "MID", "NP", "O"))
# NELHA
ALL_NELHA <- read_csv("ALL_NELHA.csv")
ALL_NELHA$month <- factor(ALL_NELHA$month, levels = c("June", "September", "December", "March"))
ALL_NELHA$Station <- factor(ALL_NELHA$Station, levels = c("NS", "MID", "NP", "O"))
# UH & NELHA
ALL_KD <- read_csv("ALL_KD.csv")
ALL_KD$month <- factor(ALL_KD$month, levels = c("June", "September", "December", "March"))
ALL_KD$Station <- factor(ALL_KD$Station, levels = c("NS", "MID", "NP", "O"))DIFFERENCES AMONG SITES
NO3NO2
########## NO3NO2 ############
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(NO3NO2 ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.00000258UH_NO3NO2_SiteDUNN <- ALL_UH %>% rstatix::dunn_test(NO3NO2 ~ Site) %>%
rstatix::add_xy_position(x = "Site")
UH_NO3NO2_Site <- ggboxplot(ALL_UH, x = "Site", y = "NO3NO2", color = "Site") +
stat_pvalue_manual(UH_NO3NO2_SiteDUNN, hide.ns = TRUE,y.position = 100,step.increase=0.05) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(NO3NO2 ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.00589NE_NO3NO2_SiteDUNN <- ALL_NELHA %>% rstatix::dunn_test(NO3NO2 ~ Site) %>%
rstatix::add_xy_position(x = "Site")
NE_NO3NO2_Site <- ggboxplot(ALL_NELHA, x = "Site", y = "NO3NO2", color = "Site") +
stat_pvalue_manual(NE_NO3NO2_SiteDUNN, hide.ns = TRUE,y.position = 110) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ggarrange(UH_NO3NO2_Site, NE_NO3NO2_Site, ncol = 2)
PO4
########## PO4 ############
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(PO4 ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.0393UH_PO4_SiteDUNN <- ALL_UH %>% rstatix::dunn_test(PO4 ~ Site) %>%
rstatix::add_xy_position(x = "Site")
UH_PO4_Site <- ggboxplot(ALL_UH, x = "Site", y = "PO4", color = "Site") +
stat_pvalue_manual(UH_PO4_SiteDUNN, hide.ns = TRUE,y.position = 62) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(PO4 ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.00830NE_PO4_SiteDUNN <- ALL_NELHA %>% rstatix::dunn_test(PO4 ~ Site) %>%
rstatix::add_xy_position(x = "Site")
NE_PO4_Site <- ggboxplot(ALL_NELHA, x = "Site", y = "PO4", color = "Site") +
stat_pvalue_manual(NE_PO4_SiteDUNN, hide.ns = TRUE,y.position = 62) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ggarrange(UH_PO4_Site, NE_PO4_Site, ncol = 2)
H4SiO4 & Si
########## H4SiO4 ############
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(H4SiO4 ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.00191UH_Si_SiteDUNN <- ALL_UH %>% rstatix::dunn_test(H4SiO4 ~ Site) %>%
rstatix::add_xy_position(x = "Site")
UH_Si_Site <- ggboxplot(ALL_UH, x = "Site", y = "H4SiO4", color = "Site") +
stat_pvalue_manual(UH_Si_SiteDUNN, hide.ns = TRUE,y.position = 1000,step.increase=0.05) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(Si ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.0000288NE_Si_SiteDUNN <- ALL_NELHA %>% rstatix::dunn_test(Si ~ Site) %>%
rstatix::add_xy_position(x = "Site")
NE_Si_Site <-ggboxplot(ALL_NELHA, x = "Site", y = "Si", color = "Site") +
stat_pvalue_manual(NE_Si_SiteDUNN, hide.ns = TRUE,y.position = 1000,step.increase=0.05) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ggarrange(UH_Si_Site, NE_Si_Site, ncol = 2)
NH4 & NH3NH4
########## NH4 ############
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(NH4 ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.0217UH_NH4_SiteDUNN <- ALL_UH %>% rstatix::dunn_test(NH4 ~ Site) %>%
rstatix::add_xy_position(x = "Site")
UH_NH4_Site <- ggboxplot(ALL_UH, x = "Site", y = "NH4", color = "Site") +
stat_pvalue_manual(UH_NH4_SiteDUNN, hide.ns = TRUE,y.position = 200,step.increase=0.05) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(NH3NH4 ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 3.15e-19NE_NH4_SiteDUNN <- ALL_NELHA %>% rstatix::dunn_test(NH3NH4 ~ Site) %>%
rstatix::add_xy_position(x = "Site")
NE_NH4_Site <- ggboxplot(ALL_NELHA, x = "Site", y = "NH3NH4", color = "Site") +
stat_pvalue_manual(NE_NH4_SiteDUNN, hide.ns = TRUE,y.position =30,step.increase=0.07) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ggarrange(UH_NH4_Site, NE_NH4_Site, ncol = 2)
TN
########## TN ############
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TN ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.00322UH_TN_SiteDUNN <- ALL_UH %>% rstatix::dunn_test(TN ~ Site) %>%
rstatix::add_xy_position(x = "Site")
UH_TN_Site <- ggboxplot(ALL_UH, x = "Site", y = "TN", color = "Site") +
stat_pvalue_manual(UH_TN_SiteDUNN, hide.ns = TRUE,y.position = 100,step.increase=0.05) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TN ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.000000295NE_TN_SiteDUNN <- ALL_NELHA %>% rstatix::dunn_test(TN ~ Site) %>%
rstatix::add_xy_position(x = "Site")
NE_TN_Site <- ggboxplot(ALL_NELHA, x = "Site", y = "TN", color = "Site") +
stat_pvalue_manual(NE_TN_SiteDUNN, hide.ns = TRUE,y.position = 150,step.increase=0.05) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ggarrange(UH_TN_Site, NE_TN_Site, ncol = 2)
TP
########## TP ############
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TP ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.000415UH_TP_SiteDUNN <- ALL_UH %>% rstatix::dunn_test(TP ~ Site) %>%
rstatix::add_xy_position(x = "Site")
UH_TP_Site <- ggboxplot(ALL_UH, x = "Site", y = "TP", color = "Site") +
stat_pvalue_manual(UH_TP_SiteDUNN, hide.ns = TRUE,y.position = 110,step.increase=0.05) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TP ~ Site)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.000537NE_TP_SiteDUNN <- ALL_NELHA %>% rstatix::dunn_test(TP ~ Site) %>%
rstatix::add_xy_position(x = "Site")
NE_TP_Site <- ggboxplot(ALL_NELHA, x = "Site", y = "TP", color = "Site") +
stat_pvalue_manual(NE_TP_SiteDUNN, hide.ns = TRUE,y.position = 70,step.increase=0.05) +
scale_color_brewer(palette = "Paired") + theme(legend.position = "none")
ggarrange(UH_TP_Site, NE_TP_Site, ncol = 2)
Factorial ANOVA – (Parameter ~ Site * Station)
########### PO4 ##############
PO4_twoway <- aov(PO4 ~ Site * Station, data = ALL_UH)
summary(PO4_twoway)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 11 967 87.89 5.924 4.00e-09 ***
## Station 3 464 154.68 10.426 1.13e-06 ***
## Site:Station 33 2098 63.58 4.285 2.87e-13 ***
## Residuals 528 7834 14.84
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1########### NO3NO2 ##############
NO3NO2_twoway <- aov(NO3NO2 ~ Site * Station, data = ALL_UH)
summary(NO3NO2_twoway)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 11 3000 272.7 5.033 1.72e-07 ***
## Station 3 3332 1110.6 20.496 1.41e-12 ***
## Site:Station 33 7476 226.5 4.181 8.32e-13 ***
## Residuals 528 28610 54.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1########### H4SiO4 ##############
H4SiO4_twoway <- aov(H4SiO4 ~ Site * Station, data = ALL_UH)
summary(H4SiO4_twoway)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 11 933399 84854 5.427 3.27e-08 ***
## Station 3 2285979 761993 48.739 < 2e-16 ***
## Site:Station 33 2174678 65899 4.215 5.87e-13 ***
## Residuals 528 8254869 15634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1########### NH4 ##############
NH4_twoway <- aov(NH4 ~ Site * Station, data = ALL_UH)
summary(NH4_twoway)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 11 2826 256.9 0.968 0.474
## Station 3 547 182.5 0.688 0.560
## Site:Station 33 8537 258.7 0.975 0.510
## Residuals 528 140113 265.4########### TN ##############
TN_twoway <- aov(TN ~ Site * Station, data = ALL_UH)
summary(TN_twoway)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 11 6359 578.1 3.552 7.80e-05 ***
## Station 3 2628 876.1 5.383 0.00118 **
## Site:Station 33 14620 443.0 2.722 1.61e-06 ***
## Residuals 528 85935 162.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1########### TP ##############
TP_twoway <- aov(TP ~ Site * Station, data = ALL_UH)
summary(TP_twoway)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 11 9224 838.6 3.347 0.000176 ***
## Station 3 817 272.2 1.087 0.354248
## Site:Station 33 8184 248.0 0.990 0.486032
## Residuals 528 132273 250.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1DIFFERENCES AMONG STATIONS – GROUPED BY SITE
PO4
########### PO4 ##############
ALL_UH %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(PO4 ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 1 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 HA 0.000935UH_PO4_Station <- ggplot(ALL_UH, aes(x = Station, y = PO4, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 40,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ALL_NELHA %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(PO4 ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 2 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 HA 0.00491
## 2 KO 0.0266NE_PO4_Station <- ggplot(ALL_NELHA, aes(x = Station, y = PO4, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 40,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggarrange(UH_PO4_Station, NE_PO4_Station, ncol = 2)
NH4 & NH3NH4
ALL_UH %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(NH4 ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 0 × 2
## # ℹ 2 variables: Site <chr>, Kruskal_Wallis_p_value <dbl>UH_NH3_Station <- ggplot(ALL_UH, aes(x = Station, y = NH4, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 40,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ALL_NELHA %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(NH3NH4 ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 2 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 HA 0.0181
## 2 W 0.0211NE_NH3NH4_Station <- ggplot(ALL_NELHA, aes(x = Station, y = NH3NH4, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 40,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggarrange(UH_NH3_Station, NE_NH3NH4_Station, ncol = 2)
NO2NO3
########### NO3NO2 ##############
ALL_UH %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(NO3NO2 ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 7 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 C 0.00414
## 2 HA 0.00000890
## 3 HO 0.0303
## 4 KP 0.0245
## 5 NP 0.0134
## 6 W 0.000000323
## 7 Z 0.000969UH_NO3NO2_Station <-ggplot(ALL_UH, aes(x = Station, y = NO3NO2, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 60,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ALL_NELHA %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(NO3NO2 ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 11 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 C 0.0000161
## 2 HA 0.00000615
## 3 HO 0.0000969
## 4 KI 0.00000918
## 5 KK 0.0172
## 6 KO 0.000115
## 7 KP 0.00408
## 8 NP 0.0000189
## 9 PT 0.000161
## 10 W 0.0000107
## 11 Z 0.0000110NE_NO3NO2_Station <-ggplot(ALL_NELHA, aes(x = Station, y = NO3NO2, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 100,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggarrange(UH_NO3NO2_Station, NE_NO3NO2_Station, ncol = 2)
H4SiO4
########### H4SiO4 ##############
ALL_UH %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(H4SiO4 ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 7 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 HA 0.00000337
## 2 HO 0.0126
## 3 KI 0.00567
## 4 KO 0.000263
## 5 NP 0.00113
## 6 W 0.000384
## 7 Z 0.00905UH_Si_Station <-ggplot(ALL_UH, aes(x = Station, y = H4SiO4, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 1000,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ALL_NELHA %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(Si ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 12 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 C 0.00411
## 2 HA 0.00000810
## 3 HO 0.00498
## 4 KI 0.00318
## 5 KK 0.0321
## 6 KO 0.0000551
## 7 KP 0.000468
## 8 M 0.00395
## 9 NP 0.0000545
## 10 PT 0.0000183
## 11 W 0.0000116
## 12 Z 0.0157NE_Si_Station <-ggplot(ALL_NELHA, aes(x = Station, y = Si, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 1000,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggarrange(UH_Si_Station, NE_Si_Station, ncol = 2)
########### TN ##############
ALL_UH %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(TN ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 1 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 HA 0.00113UH_TN_Station <- ggplot(ALL_UH, aes(x = Station, y = TN, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 200,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ALL_NELHA %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(TN ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 8 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 HA 0.000194
## 2 KK 0.000496
## 3 KO 0.00418
## 4 KP 0.00155
## 5 M 0.00470
## 6 PT 0.000868
## 7 W 0.000861
## 8 Z 0.0181NE_TN_Station <- ggplot(ALL_NELHA, aes(x = Station, y = TN, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 200,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggarrange(UH_TN_Station, NE_TN_Station, ncol = 2)
########### TP ##############
ALL_UH %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(TP ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 1 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 HO 0.0293UH_TP_Station <- ggplot(ALL_UH, aes(x = Station, y = TP, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 75,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ALL_NELHA %>% group_by(Site) %>%
summarize(Kruskal_Wallis_p_value = kruskal.test(TP ~ Station)$p.value) %>%
filter(Kruskal_Wallis_p_value < 0.05)## # A tibble: 1 × 2
## Site Kruskal_Wallis_p_value
## <chr> <dbl>
## 1 HA 0.000706NE_TP_Station <- ggplot(ALL_NELHA, aes(x = Station, y = TP, fill = Station)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 45,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggarrange(UH_TP_Station, NE_TP_Station, ncol = 2)
Extra Figures – Parameter ~ Station * Depth / Site
###### graphs #########
## UH
UH_TP <- ggplot(ALL_UH, aes(x = Station, y = TP, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
UH_TN <-ggplot(ALL_UH, aes(x = Station, y = TN, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
UH_NH4 <-ggplot(ALL_UH[-375,], aes(x = Station, y = NH4, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
UH_NO3NO2 <- ggplot(ALL_UH, aes(x = Station, y = NO3NO2, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
UH_PO4 <- ggplot(ALL_UH, aes(x = Station, y = PO4, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
UH_H4SiO4 <-ggplot(ALL_UH, aes(x = Station, y = H4SiO4, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
UH_Turb <- ggplot(ALL_UH, aes(x = Station, y = Turbidity, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
## NELHA
NE_TP <-ggplot(ALL_NELHA, aes(x = Station, y = TP, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
NE_TN <- ggplot(ALL_NELHA, aes(x = Station, y = TN, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
NE_NH3NH4 <- ggplot(ALL_NELHA, aes(x = Station, y = NH3NH4, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
NE_NO3NO2 <- ggplot(ALL_NELHA, aes(x = Station, y = NO3NO2, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
NE_PO4 <- ggplot(ALL_NELHA, aes(x = Station, y = PO4, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
NE_Si <- ggplot(ALL_NELHA, aes(x = Station, y = Si, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
NE_Turb <-ggplot(ALL_NELHA, aes(x = Station, y = Turbidity, fill = Sample_ID)) +
geom_boxplot(lwd = 0.2,outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() + scale_fill_brewer(palette = "Set2")
## arrange to compare labs ##
ggpubr::ggarrange(UH_TP, NE_TP, ncol = 2)
ggarrange(UH_TN, NE_TN, ncol = 2)
ggarrange(UH_NO3NO2, NE_NO3NO2, ncol = 2)
ggarrange(UH_NH4, NE_NH3NH4, ncol = 2)
ggarrange(UH_PO4, NE_PO4, ncol = 2)
ggarrange(UH_Turb, NE_Turb, ncol = 2)
ggarrange(UH_H4SiO4, NE_Si, ncol = 2)
DIFFERENCES AMONG MONTHS
Turbidity
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(Turbidity ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 3.86e-19#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_Turb_Month <- ggboxplot(ALL_UH, x = "month", y = "Turbidity", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(Turbidity ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 3.01e-60#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_Turb_Month <- ggboxplot(ALL_NELHA, x = "month", y = "Turbidity", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_Turb_Month, NE_Turb_Month, ncol = 2)
### PO4
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(PO4 ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 3.55e-93#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_PO4_Month <- ggboxplot(ALL_UH, x = "month", y = "PO4", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(PO4 ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 3.67e-59#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_PO4_Month <- ggboxplot(ALL_NELHA, x = "month", y = "PO4", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_PO4_Month, NE_PO4_Month, ncol = 2)
NH4 & NH3NH4
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(NH4 ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.0000398#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_NH4_Month <- ggboxplot(ALL_UH, x = "month", y = "NH4", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 200,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(NH3NH4 ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 1.01e-21#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_NH4_Month <- ggboxplot(ALL_NELHA, x = "month", y = "NH3NH4", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 30,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_NH4_Month, NE_NH4_Month, ncol = 2)
H4SiO4 & Si
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(H4SiO4 ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 2.62e-27#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_Si_Month <- ggboxplot(ALL_UH, x = "month", y = "H4SiO4", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 1000,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(Si ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 1.45e-27#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_Si_Month <- ggboxplot(ALL_NELHA, x = "month", y = "Si", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 1000,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_Si_Month, NE_Si_Month, ncol = 2)
TN
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TN ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 1.34e-49#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_TN_Month <- ggboxplot(ALL_UH, x = "month", y = "TN", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 230,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TN ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 1.67e-31#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_TN_Month <- ggboxplot(ALL_NELHA, x = "month", y = "TN", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 230,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_TN_Month, NE_TN_Month, ncol = 2)
TP
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TP ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 3.64e-36#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_TP_Month <- ggboxplot(ALL_UH, x = "month", y = "TP", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 80,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TP ~ month)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 3.27e-88#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_TP_Month <- ggboxplot(ALL_NELHA, x = "month", y = "TP", fill = "month") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 70,step.increase = 0.06) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_TP_Month, NE_TP_Month, ncol = 2)
DIFFERENCES AMONG MONTHS - GROUPED BY SITE
ggplot(ALL_UH, aes(x = month, y = Turbidity, fill = month)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 1,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggplot(ALL_UH, aes(x = month, y = PO4, fill = month)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 15,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggplot(ALL_UH, aes(x = month, y = NH4, fill = month)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 20,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggplot(ALL_UH, aes(x = month, y = NO3NO2, fill = month)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 45,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggplot(ALL_UH, aes(x = month, y = H4SiO4, fill = month)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 700,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggplot(ALL_UH, aes(x = month, y = TP, fill = month)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 75,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
ggplot(ALL_UH, aes(x = month, y = TN, fill = month)) +
geom_boxplot(lwd = 0.2, outlier.size = 0.3) + facet_wrap(~Site) + theme_bw() +
scale_fill_brewer(palette = "Set2") + geom_pwc(method = "dunn.test",
label = "p.signif", hide.ns = TRUE, y.position = 110,
label.size = 3, tip.length = 0.02,step.increase = 0.08)
Extra Figures – Parameter ~ Month * Depth / Site
## UH
UH_TP2 <- ggplot(ALL_UH, aes(x = as.numeric(month), y = TP, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
UH_PO42 <- ggplot(ALL_UH, aes(x = as.numeric(month), y = PO4, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
UH_TN2 <- ggplot(ALL_UH, aes(x = as.numeric(month), y = TN, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
UH_H4SiO42 <- ggplot(ALL_UH, aes(x = as.numeric(month), y = H4SiO4, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
UH_NO3NO22 <- ggplot(ALL_UH, aes(x = as.numeric(month), y = NO3NO2, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
UH_NH42 <- ggplot(ALL_UH, aes(x = as.numeric(month), y = NH4, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
## NELHA
NE_TP2 <- ggplot(ALL_NELHA, aes(x = as.numeric(month), y = TP, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
NE_PO42 <- ggplot(ALL_NELHA, aes(x = as.numeric(month), y = PO4, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
NE_TN2 <- ggplot(ALL_NELHA, aes(x = as.numeric(month), y = TN, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
NE_Si2 <- ggplot(ALL_NELHA, aes(x = as.numeric(month), y = Si, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
NE_NO3NO22 <- ggplot(ALL_NELHA, aes(x = as.numeric(month), y = NO3NO2, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
NE_NH3NH42 <- ggplot(ALL_NELHA, aes(x = as.numeric(month), y = NH3NH4, color = Station)) +
geom_point(size=0.7) + geom_smooth(method = "lm", se = FALSE, size=0.5) +
facet_wrap(~Site) + scale_colour_brewer(palette = "Set2") +
scale_x_continuous(breaks = 1:4, labels = c("June", "September", "December", "March")) +
theme_bw() + labs(x = NULL) + theme(axis.text.x = element_text(angle = 20, hjust = 1, size = 8))
ggarrange(UH_TP2, NE_TP2, ncol = 2)
ggarrange(UH_TN2, NE_TN2, ncol = 2)
ggarrange(UH_NO3NO22, NE_NO3NO22, ncol = 2)
ggarrange(UH_NH42, NE_NH3NH42, ncol = 2)
ggarrange(UH_PO42, NE_PO42, ncol = 2)
ggarrange(UH_H4SiO42, NE_Si2, ncol = 2)
DIFFERENCES AMONG LABS
ggboxplot(ALL_KD, x = "Lab", y = "Turbidity", fill = "Lab") +
stat_compare_means(method = 'kruskal.test') +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggboxplot(ALL_KD, x = "Lab", y = "PO4", fill = "Lab") +
stat_compare_means(method = 'kruskal.test') +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggboxplot(ALL_KD, x = "Lab", y = "TP", fill = "Lab") +
stat_compare_means(method = 'kruskal.test') +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggboxplot(ALL_KD, x = "Lab", y = "TN", fill = "Lab") +
stat_compare_means(method = 'kruskal.test') +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggboxplot(ALL_KD, x = "Lab", y = "NH4", fill = "Lab") +
stat_compare_means(method = 'kruskal.test') +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggboxplot(ALL_KD, x = "Lab", y = "NO3NO2", fill = "Lab") +
stat_compare_means(method = 'kruskal.test') +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
Extra Figures – Parameter ~ Station * Lab / Site
ggplot(ALL_KD, aes(x = Station, y = TP, fill = Lab)) +
geom_boxplot(lwd = 0.3, outlier.size = 0.5) + facet_wrap(~ Site) +
scale_fill_brewer(palette = "Set2") + theme_bw()
ggplot(ALL_KD, aes(x = Station, y = PO4, fill = Lab)) +
geom_boxplot(lwd = 0.3, outlier.size = 0.5) + facet_wrap(~ Site) +
scale_fill_brewer(palette = "Set2") + theme_bw()
ggplot(ALL_KD, aes(x = Station, y = TN, fill = Lab)) +
geom_boxplot(lwd = 0.3, outlier.size = 0.5) + facet_wrap(~ Site) +
scale_fill_brewer(palette = "Set2") + theme_bw()
ggplot(ALL_KD, aes(x = Station, y = NO3NO2, fill = Lab)) +
geom_boxplot(lwd = 0.3, outlier.size = 0.5) + facet_wrap(~ Site) +
scale_fill_brewer(palette = "Set2") + theme_bw()
ggplot(ALL_KD, aes(x = Station, y = Turbidity, fill = Lab)) +
geom_boxplot(lwd = 0.3, outlier.size = 0.5) + facet_wrap(~ Site) +
scale_fill_brewer(palette = "Set2") + theme_bw()
ggplot(ALL_KD, aes(x = Station, y = NH4, fill = Lab)) +
geom_boxplot(lwd = 0.3, outlier.size = 0.5) + facet_wrap(~ Site, scales = "free_y") +
scale_fill_brewer(palette = "Set2") + theme_bw()
DIFFERENCES AMONG REGIONS
ALL_UH <- ALL_UH %>%
mutate(Region = case_when(
Site %in% c("KK", "KI", "M") ~ "North",
Site %in% c("C", "NP", "Z") ~ "Net Pen",
Site %in% c("W", "KP", "HO") ~ "Keahole Point",
Site %in% c("HA", "KO", "PT") ~ "South",
TRUE ~ "Other" ))
ALL_NELHA <- ALL_NELHA %>%
mutate(Region = case_when(
Site %in% c("KK", "KI", "M") ~ "North",
Site %in% c("C", "NP", "Z") ~ "Net Pen",
Site %in% c("W", "KP", "HO") ~ "Keahole Point",
Site %in% c("HA", "KO", "PT") ~ "South",
TRUE ~ "Other" ))PO4
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(PO4 ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.0122#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_PO4_Reg <- ggboxplot(ALL_UH, x = "Region", y = "PO4", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 40,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(PO4 ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.00157#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_PO4_Reg <- ggboxplot(ALL_NELHA, x = "Region", y = "PO4", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 40,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_PO4_Reg, NE_PO4_Reg, ncol = 2)
NO3NO2
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(NO3NO2 ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.0000273#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_NO3NO2_Reg <- ggboxplot(ALL_UH, x = "Region", y = "NO3NO2", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 80,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(NO3NO2 ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.0114#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_NO3NO2_Reg <- ggboxplot(ALL_NELHA, x = "Region", y = "NO3NO2", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 80,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_NO3NO2_Reg, NE_NO3NO2_Reg, ncol = 2)
H4SiO4 & Si
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(H4SiO4 ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.0585#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_H4SiO4_Reg <- ggboxplot(ALL_UH, x = "Region", y = "H4SiO4", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 1000,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(Si ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.000720#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_Si_Reg <- ggboxplot(ALL_NELHA, x = "Region", y = "Si", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 1000,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_H4SiO4_Reg, NE_Si_Reg, ncol = 2)
NH4 & NH3NH4
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(NH4 ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.00530#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_NH4_Reg <- ggboxplot(ALL_UH, x = "Region", y = "NH4", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 100,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(NH3NH4 ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 1.00e-17#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_NH3NH4_Reg <- ggboxplot(ALL_NELHA, x = "Region", y = "NH3NH4", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 45,step.increase = 0.07) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_NH4_Reg, NE_NH3NH4_Reg, ncol = 2)
TP
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TP ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.0000676#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_TP_Reg <- ggboxplot(ALL_UH, x = "Region", y = "TP", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 100,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TP ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.000427#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_TP_Reg <- ggboxplot(ALL_NELHA, x = "Region", y = "TP", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 45,step.increase = 0.07) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_TP_Reg, NE_TP_Reg, ncol = 2)
### TN
ALL_UH %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TN ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.158#dunn.test(ALL_UH$Turbidity, g = ALL_UH$month, method = "bonferroni")
UH_TN_Reg <- ggboxplot(ALL_UH, x = "Region", y = "TN", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 150,step.increase = 0.05) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ALL_NELHA %>% summarize(Kruskal_Wallis_p_value = kruskal.test(TN ~ Region)$p.value) ## # A tibble: 1 × 1
## Kruskal_Wallis_p_value
## <dbl>
## 1 0.000379#dunn.test(ALL_NELHA$Turbidity, g = ALL_NELHA$month, method = "bonferroni")
NE_TN_Reg <- ggboxplot(ALL_NELHA, x = "Region", y = "TN", fill = "Region") +
geom_pwc(method = "dunn.test", label = "p.signif", hide.ns = TRUE, y.position = 150,step.increase = 0.07) +
scale_fill_brewer(palette = "Pastel2") + theme(legend.position = "none")
ggarrange(UH_TN_Reg, NE_TN_Reg, ncol = 2)
07/02
Bootstraping
# - NH3NH4, NO2NO3 - all data
# - Each region (north, etc) — but do it for surface only & NS and NP only
UH_subset <- ALL_UH %>% select(Region, Site, Station,Sample_ID, Lab, NH4, NO3NO2) %>%
filter(Station %in% c("NS", "NP")) %>% filter(Sample_ID == "S")
NELHA_subset <- ALL_NELHA %>% select(Region, Site, Station, Sample_ID,Lab, NH3NH4, NO3NO2) %>%
filter(Station %in% c("NS", "NP")) %>% filter(Sample_ID == "S")
# explore initial distributions
ggplot(UH_subset %>% filter(Station == "NS"), aes(x = NO3NO2, fill = Region)) +
geom_density(alpha = 0.4) + # Density plot for NO3NO2 distribution
facet_wrap(~ Region, scales = "free") 
## bootstrap
NELHA_boot <- NELHA_subset %>%
group_by(Region, Station) %>%
group_modify(~ {NO3NO2_bootstrap <- .x %>%
specify(response = NO3NO2) %>%generate(reps = 1000, type = "bootstrap")
tibble( NO3NO2_bootstrap = list(NO3NO2_bootstrap))
})
NELHA_boot2 <- NELHA_subset %>%
group_by(Region, Station) %>%
group_modify(~ {NH3NH4_bootstrap <- .x %>%
specify(response = NH3NH4) %>%generate(reps = 1000, type = "bootstrap")
tibble( NH3NH4_bootstrap = list(NH3NH4_bootstrap))
})
# Compute confidence intervals for each bootstrapped dataset
NELHA_boot_ci <- NELHA_boot %>%
rowwise() %>%mutate( NO3NO2_summary = list(NO3NO2_bootstrap %>%summarise(stat = mean(NO3NO2))),
NO3NO2_ci = list( NO3NO2_summary %>%
get_confidence_interval( level = 0.95, type = "percentile") ) ) %>%
unnest(cols = c(NO3NO2_summary, NO3NO2_ci))
NELHA_boot_ci2 <- NELHA_boot2 %>%
rowwise() %>%mutate( NH3NH4_summary = list(NH3NH4_bootstrap %>%summarise(stat = mean(NH3NH4))),
NO3NO2_ci = list( NH3NH4_summary %>%
get_confidence_interval( level = 0.95, type = "percentile") ) ) %>%
unnest(cols = c(NH3NH4_summary, NO3NO2_ci))
# View the confidence intervals
NELHA_boot_ci## # A tibble: 8,000 × 7
## # Groups: Region, Station [8]
## Region Station NO3NO2_bootstrap replicate stat lower_ci upper_ci
## <chr> <fct> <list> <int> <dbl> <dbl> <dbl>
## 1 Keahole Point NS <infer [12,000 × 2]> 1 28.2 12.4 35.1
## 2 Keahole Point NS <infer [12,000 × 2]> 2 30.4 12.4 35.1
## 3 Keahole Point NS <infer [12,000 × 2]> 3 31.1 12.4 35.1
## 4 Keahole Point NS <infer [12,000 × 2]> 4 23.2 12.4 35.1
## 5 Keahole Point NS <infer [12,000 × 2]> 5 7.09 12.4 35.1
## 6 Keahole Point NS <infer [12,000 × 2]> 6 26.3 12.4 35.1
## 7 Keahole Point NS <infer [12,000 × 2]> 7 17.7 12.4 35.1
## 8 Keahole Point NS <infer [12,000 × 2]> 8 20.5 12.4 35.1
## 9 Keahole Point NS <infer [12,000 × 2]> 9 24.4 12.4 35.1
## 10 Keahole Point NS <infer [12,000 × 2]> 10 25.8 12.4 35.1
## # ℹ 7,990 more rowsNELHA_boot_ci2## # A tibble: 8,000 × 7
## # Groups: Region, Station [8]
## Region Station NH3NH4_bootstrap replicate stat lower_ci upper_ci
## <chr> <fct> <list> <int> <dbl> <dbl> <dbl>
## 1 Keahole Point NS <infer [12,000 × 2]> 1 4.14 2.20 4.32
## 2 Keahole Point NS <infer [12,000 × 2]> 2 2.39 2.20 4.32
## 3 Keahole Point NS <infer [12,000 × 2]> 3 2.76 2.20 4.32
## 4 Keahole Point NS <infer [12,000 × 2]> 4 3.67 2.20 4.32
## 5 Keahole Point NS <infer [12,000 × 2]> 5 2.94 2.20 4.32
## 6 Keahole Point NS <infer [12,000 × 2]> 6 4.64 2.20 4.32
## 7 Keahole Point NS <infer [12,000 × 2]> 7 3.37 2.20 4.32
## 8 Keahole Point NS <infer [12,000 × 2]> 8 3.15 2.20 4.32
## 9 Keahole Point NS <infer [12,000 × 2]> 9 2.17 2.20 4.32
## 10 Keahole Point NS <infer [12,000 × 2]> 10 2.56 2.20 4.32
## # ℹ 7,990 more rows#### filter & plot #####
# NO3NO2
NELHA_boot_NO3NO2_NS <- NELHA_boot_ci %>% filter(Station == "NS") %>%
ggplot(aes(x = stat, fill = Region)) + geom_density(alpha = 0.4) +
geom_vline(aes(xintercept = mean), data = . %>% group_by(Region) %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) + labs(tag = "NS") +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
scale_fill_brewer(palette = "Set2") + facet_wrap(~ Region, scales = "free") + theme_minimal() +
theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
theme(plot.tag = element_text(face = "bold", size=12))
NELHA_boot_NO3NO2_NP <- NELHA_boot_ci %>% filter(Station == "NP") %>%
ggplot(aes(x = stat, fill = Region)) + geom_density(alpha = 0.4) +
geom_vline(aes(xintercept = mean), data = . %>% group_by(Region) %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) + labs(tag = "NP") +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
scale_fill_brewer(palette = "Set2") + facet_wrap(~ Region, scales = "free") + theme_minimal() +
theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
theme(plot.tag = element_text(face = "bold", size=12))
NE_NO3NO2_bootr <- ggarrange(NELHA_boot_NO3NO2_NS, NELHA_boot_NO3NO2_NP, ncol = 2,common.legend = TRUE, legend = "bottom")
annotate_figure(NE_NO3NO2_bootr, top = text_grob("NELHA -- NO3NO2 -- NearShore / NetPen",
color = "black", face = "bold", size = 12))
# NH3NH4
NELHA_boot_NH3NH4_NS <- NELHA_boot_ci2 %>% filter(Station == "NS") %>%
ggplot(aes(x = stat, fill = Region)) + geom_density(alpha = 0.4) +
geom_vline(aes(xintercept = mean), data = . %>% group_by(Region) %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) + labs(tag = "NS") +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
scale_fill_brewer(palette = "Set2") + facet_wrap(~ Region, scales = "free") + theme_minimal() +
theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
theme(plot.tag = element_text(face = "bold", size=12))
NELHA_boot_NH3NH4_NP <- NELHA_boot_ci2 %>% filter(Station == "NP") %>%
ggplot(aes(x = stat, fill = Region)) + geom_density(alpha = 0.4) +
geom_vline(aes(xintercept = mean), data = . %>% group_by(Region) %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35)+ labs(tag = "NP") +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
scale_fill_brewer(palette = "Set2") + facet_wrap(~ Region, scales = "free") + theme_minimal() +
theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
theme(plot.tag = element_text(face = "bold", size=12))
NE_NH3NH4_bootr <- ggarrange(NELHA_boot_NH3NH4_NS, NELHA_boot_NH3NH4_NP, ncol = 2,common.legend = TRUE, legend = "bottom")
annotate_figure(NE_NH3NH4_bootr, top = text_grob("NELHA -- NH3NH4 -- NearShore / NetPen",
color = "black", face = "bold", size = 12))
###### combined regions #######
## bootstrap
NELHA_boot_reg <- NELHA_subset %>%
group_by(Station) %>%
group_modify(~ {NO3NO2_bootstrap <- .x %>%
specify(response = NO3NO2) %>%generate(reps = 1000, type = "bootstrap")
tibble( NO3NO2_bootstrap = list(NO3NO2_bootstrap))
})
NELHA_boot_reg2 <- NELHA_subset %>%
group_by(Station) %>%
group_modify(~ {NH3NH4_bootstrap <- .x %>%
specify(response = NH3NH4) %>%generate(reps = 1000, type = "bootstrap")
tibble( NH3NH4_bootstrap = list(NH3NH4_bootstrap))
})
## ci
NELHA_boot_regci <- NELHA_boot_reg %>%
rowwise() %>%mutate( NO3NO2_summary = list(NO3NO2_bootstrap %>%summarise(stat = mean(NO3NO2))),
NO3NO2_ci = list( NO3NO2_summary %>%
get_confidence_interval( level = 0.95, type = "percentile") ) ) %>%
unnest(cols = c(NO3NO2_summary, NO3NO2_ci))
NELHA_boot_regci2 <- NELHA_boot_reg2 %>%
rowwise() %>%mutate( NH3NH4_summary = list(NH3NH4_bootstrap %>%summarise(stat = mean(NH3NH4))),
NO3NO2_ci = list( NH3NH4_summary %>%
get_confidence_interval( level = 0.95, type = "percentile") ) ) %>%
unnest(cols = c(NH3NH4_summary, NO3NO2_ci))
## figs
NELHA_boot_regci %>%
ggplot(aes(x = stat, fill = Station)) + geom_density(alpha = 0.4) +
geom_vline(aes(xintercept = mean), data = . %>% group_by(Station) %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
scale_fill_brewer(palette = "Set1") + facet_wrap(~ Station, scales = "free") + theme_minimal() +
theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
labs(title = "NO3NO2 at NearShore (NS) and NetPen (NP)") +
theme(plot.title = element_text(hjust=0.5))
NELHA_boot_regci2 %>%
ggplot(aes(x = stat, fill = Station)) + geom_density(alpha = 0.4) +
geom_vline(aes(xintercept = mean), data = . %>% group_by(Station) %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
scale_fill_brewer(palette = "Set1") + facet_wrap(~ Station, scales = "free") + theme_minimal() +
theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
labs(title = "NH3NH4 at NearShore (NS) and NetPen (NP)") +
theme(plot.title = element_text(hjust=0.5))
###### combined regions -- UH #######
## bootstrap
UH_boot_reg <- UH_subset %>%
group_by(Station) %>%
group_modify(~ {NO3NO2_bootstrap <- .x %>%
specify(response = NO3NO2) %>%generate(reps = 1000, type = "bootstrap")
tibble( NO3NO2_bootstrap = list(NO3NO2_bootstrap))
})
UH_boot_reg2 <- UH_subset %>%
group_by(Station) %>%
group_modify(~ {NH4_bootstrap <- .x %>%
specify(response = NH4) %>%generate(reps = 1000, type = "bootstrap")
tibble( NH4_bootstrap = list(NH4_bootstrap))
})
## ci
UH_boot_regci <- UH_boot_reg %>%
rowwise() %>%mutate( NO3NO2_summary = list(NO3NO2_bootstrap %>%summarise(stat = mean(NO3NO2))),
NO3NO2_ci = list( NO3NO2_summary %>%
get_confidence_interval( level = 0.95, type = "percentile") ) ) %>%
unnest(cols = c(NO3NO2_summary, NO3NO2_ci))
UH_boot_regci2 <- UH_boot_reg2 %>%
rowwise() %>%mutate( NH4_summary = list(NH4_bootstrap %>%summarise(stat = mean(NH4))),
NH4_ci = list( NH4_summary %>%
get_confidence_interval( level = 0.95, type = "percentile") ) ) %>%
unnest(cols = c(NH4_summary, NH4_ci))
## figs
UH_boot_regci %>%
ggplot(aes(x = stat, fill = Station)) + geom_density(alpha = 0.4) +
geom_vline(aes(xintercept = mean), data = . %>% group_by(Station) %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
scale_fill_brewer(palette = "Set1") + facet_wrap(~ Station, scales = "free") + theme_minimal() +
theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
labs(title = "NO3NO2 at NearShore (NS) and NetPen (NP)") +
theme(plot.title = element_text(hjust=0.5))
UH_boot_regci2 %>%
ggplot(aes(x = stat, fill = Station)) + geom_density(alpha = 0.4) +
geom_vline(aes(xintercept = mean), data = . %>% group_by(Station) %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
scale_fill_brewer(palette = "Set1") + facet_wrap(~ Station, scales = "free") + theme_minimal() +
theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
labs(title = "NH4 at NearShore (NS) and NetPen (NP)") +
theme(plot.title = element_text(hjust=0.5))
### ALL UH ####
ALL_UH_boot <- UH_subset %>%
group_modify(~ {NO3NO2_bootstrap <- .x %>%
specify(response = NO3NO2) %>%generate(reps = 1000, type = "bootstrap")
tibble( NO3NO2_bootstrap = list(NO3NO2_bootstrap))
})
ALL_UH_boot2 <- UH_subset %>%
group_modify(~ {NH4_bootstrap <- .x %>%
specify(response = NH4) %>%generate(reps = 1000, type = "bootstrap")
tibble( NH4_bootstrap = list(NH4_bootstrap))
})
ALL_UH_bootci <- ALL_UH_boot %>%
rowwise() %>%mutate( NO3NO2_summary = list(NO3NO2_bootstrap %>%summarise(stat = mean(NO3NO2))),
NO3NO2_ci = list( NO3NO2_summary %>%
get_confidence_interval( level = 0.95, type = "percentile") ) ) %>%
unnest(cols = c(NO3NO2_summary, NO3NO2_ci))
ALL_UH_bootci2 <- ALL_UH_boot2 %>%
rowwise() %>%mutate( NH4_summary = list(NH4_bootstrap %>%summarise(stat = mean(NH4))),
NH4_ci = list( NH4_summary %>%
get_confidence_interval( level = 0.95, type = "percentile") ) ) %>%
unnest(cols = c(NH4_summary, NH4_ci))
UH_NO3_boot <- ALL_UH_bootci %>%
ggplot(aes(x = stat)) + geom_density(fill="#6FBFD0",alpha = 0.6) +
geom_vline(aes(xintercept = mean), data = . %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
theme_minimal() + theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
labs(title = "UH NO3NO2") + theme(plot.title = element_text(hjust=0.5))
UH_NH4_boot <- ALL_UH_bootci2 %>%
ggplot(aes(x = stat)) + geom_density(fill="#6FBFD0",alpha = 0.6) +
geom_vline(aes(xintercept = mean), data = . %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
theme_minimal() + theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
labs(title = "UH NH4") + theme(plot.title = element_text(hjust=0.5))
### ALL NELHA ####
ALL_NELHA_boot <- NELHA_subset %>%
group_modify(~ {NO3NO2_bootstrap <- .x %>%
specify(response = NO3NO2) %>%generate(reps = 1000, type = "bootstrap")
tibble( NO3NO2_bootstrap = list(NO3NO2_bootstrap))
})
ALL_NELHA_boot2 <- NELHA_subset %>%
group_modify(~ {NH3NH4_bootstrap <- .x %>%
specify(response = NH3NH4) %>%generate(reps = 1000, type = "bootstrap")
tibble( NH3NH4_bootstrap = list(NH3NH4_bootstrap))
})
ALL_NELHA_bootci <- ALL_NELHA_boot %>%
rowwise() %>%mutate( NO3NO2_summary = list(NO3NO2_bootstrap %>%summarise(stat = mean(NO3NO2))),
NO3NO2_ci = list( NO3NO2_summary %>%
get_confidence_interval( level = 0.95, type = "percentile") ) ) %>%
unnest(cols = c(NO3NO2_summary, NO3NO2_ci))
ALL_NELHA_bootci2 <- ALL_NELHA_boot2 %>%
rowwise() %>%mutate( NH3NH4_summary = list(NH3NH4_bootstrap %>%summarise(stat = mean(NH3NH4))),
NH3NH4_ci = list( NH3NH4_summary %>%
get_confidence_interval( level = 0.95, type = "percentile") ) ) %>%
unnest(cols = c(NH3NH4_summary, NH3NH4_ci))
NELHA_NO3_boot <- ALL_NELHA_bootci %>%
ggplot(aes(x = stat)) + geom_density(fill="#E5DC5A",alpha = 0.6) +
geom_vline(aes(xintercept = mean), data = . %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
theme_minimal() + theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
labs(title = "NELHA NO3NO2") + theme(plot.title = element_text(hjust=0.5))
NELHA_NH4_boot <- ALL_NELHA_bootci2 %>%
ggplot(aes(x = stat)) + geom_density(fill="#E5DC5A",alpha = 0.6) +
geom_vline(aes(xintercept = mean), data = . %>% summarise(mean = mean(stat)),
color = "black", linetype = "solid", size = 0.35) +
geom_vline(aes(xintercept = lower_ci), color = "black", linetype = "dashed", size = 0.3) +
geom_vline(aes(xintercept = upper_ci), color = "black", linetype = "dashed", size = 0.3) +
theme_minimal() + theme(panel.border = element_rect(color = "black", linewidth = 0.5, fill = NA)) +
labs(title = "NELHA NH3NH4") + theme(plot.title = element_text(hjust=0.5))
####
ggarrange(UH_NO3_boot, NELHA_NO3_boot, ncol = 2)
ggarrange(UH_NH4_boot, NELHA_NH4_boot, ncol = 2)