CAPA Annual Colony Loss Survey Data Analysis
CAPABEESURV
Brendan Daisley
Jan 25, 2025
Related article: Daisley et al. (2025).Antibiotic use,
air pollution, and seasonal temperature are strong predictors of honey
bee mortality in Canada . In draft
Related GitHub repository: https://github.com/bdaisley/CAPABEESURV
The R code provided below was used to generate the following figures and table presented in related article:
- Figure 1
- Figure 2
- Figure 3
- Table 1
- Supplementary Data 2
Load required packages
library(sf)
library(tidyr)
library(MuMIn)
library(dplyr)
library(purrr)
library(readxl)
library(tibble)
library(future)
library(eoffice)
library(ggplot2)
library(stringr)
library(cowplot)
library(ggrepel)
library(ggridges)
library(pdftools)
library(reshape2)
library(reactable)
library(tabulapdf)
library(packcircles)
library(rnaturalearth)
theme_set(theme_bw())
sessionInfo()## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_Canada.utf8 LC_CTYPE=English_Canada.utf8
## [3] LC_MONETARY=English_Canada.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_Canada.utf8
##
## time zone: America/Toronto
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rnaturalearth_1.0.1 packcircles_0.3.6 tabulapdf_1.0.5-5
## [4] reactable_0.4.4.9000 reshape2_1.4.4 pdftools_3.4.1
## [7] ggridges_0.5.6 ggrepel_0.9.6 cowplot_1.1.3
## [10] stringr_1.5.1 ggplot2_3.5.1 eoffice_0.2.2
## [13] future_1.34.0 tibble_3.2.1 readxl_1.4.3
## [16] purrr_1.0.2 dplyr_1.1.4 MuMIn_1.48.4
## [19] tidyr_1.3.1 sf_1.0-18
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.3 rlang_1.1.4 magrittr_2.0.3
## [4] e1071_1.7-16 compiler_4.4.1 png_0.1-8
## [7] systemfonts_1.1.0 vctrs_0.6.5 devEMF_4.5
## [10] pkgconfig_2.0.3 fastmap_1.2.0 backports_1.5.0
## [13] magick_2.8.5 utf8_1.2.4 rmarkdown_2.29
## [16] tzdb_0.4.0 ragg_1.3.3 xfun_0.49
## [19] cachem_1.1.0 jsonlite_1.8.9 uuid_1.2-1
## [22] terra_1.7-83 broom_1.0.7 parallel_4.4.1
## [25] R6_2.5.1 R.devices_2.17.2 bslib_0.8.0
## [28] stringi_1.8.4 parallelly_1.38.0 jquerylib_0.1.4
## [31] cellranger_1.1.0 Rcpp_1.0.12 knitr_1.48
## [34] base64enc_0.1-3 R.utils_2.12.3 readr_2.1.5
## [37] Matrix_1.7-1 tidyselect_1.2.1 rstudioapi_0.17.1
## [40] yaml_2.3.10 codetools_0.2-20 listenv_0.9.1
## [43] qpdf_1.3.4 lattice_0.22-6 plyr_1.8.9
## [46] withr_3.0.2 flextable_0.9.7 askpass_1.2.1
## [49] evaluate_1.0.1 gridGraphics_0.5-1 rJava_1.0-11
## [52] units_0.8-5 proxy_0.4-27 zip_2.3.1
## [55] xml2_1.3.6 pillar_1.9.0 KernSmooth_2.23-24
## [58] stats4_4.4.1 plotly_4.10.4 generics_0.1.3
## [61] hms_1.1.3 munsell_0.5.1 scales_1.3.0
## [64] globals_0.16.3 class_7.3-22 glue_1.7.0
## [67] gdtools_0.4.0 lazyeval_0.2.2 tools_4.4.1
## [70] data.table_1.16.2 fs_1.6.5 grid_4.4.1
## [73] colorspace_2.1-1 nlme_3.1-166 cli_3.6.3
## [76] textshaping_0.4.0 rvg_0.3.4 officer_0.6.7
## [79] fontBitstreamVera_0.1.1 fansi_1.0.6 viridisLite_0.4.2
## [82] gtable_0.3.6 R.methodsS3_1.8.2 yulab.utils_0.1.7
## [85] sass_0.4.9 digest_0.6.36 fontquiver_0.2.1
## [88] classInt_0.4-10 ggplotify_0.1.2 htmlwidgets_1.6.4
## [91] htmltools_0.5.8.1 R.oo_1.27.0 lifecycle_1.0.4
## [94] httr_1.4.7 fontLiberation_0.1.0 openssl_2.2.2Acquire datasets
Run the R scripts below to automatically download CAPA surveys, climate datasets, and air quality datasets
#source("https://bdaisley.github.io/CAPABEESURV/scripts/get_capa_survey_data.R")
#source("https://bdaisley.github.io/CAPABEESURV/scripts/get_air_quality_data.R")
#source("https://bdaisley.github.io/CAPABEESURV/scripts/get_climate_data.R")Import datasets
capa.df <- readxl::read_xlsx("tables/CAPA_survey_data.xlsx",
sheet="01 - CAPA Survey Data")
capa.df <- readxl::read_xlsx("tables/Supplementary_Data_1.xlsx",
sheet="S1 - CAPA Survey Data")
capa.df.filt <- capa.df %>%
mutate_at(vars(respondent_colonies_percent,
bac_spring_oxytet, bac_fall_oxytet,
nos_spring_fum,
nos_fall_fum,
var_spring_other,
var_fall_other), as.numeric) %>%
mutate(respondent_colonies_percent = ifelse(
respondent_colonies_percent > 100,
100,
respondent_colonies_percent)) %>%
mutate(oxytet_spring_fall = (bac_spring_oxytet +
bac_fall_oxytet)/2) %>%
mutate(fum_spring_fall = (nos_spring_fum +
nos_fall_fum)/2) %>%
mutate(var_spring_fall = (var_spring_other +
var_fall_other)/2) %>%
mutate(winterloss= winterloss_percent)
capa.melt <- melt(capa.df.filt , id.var=c("year",
"province",
"winterloss",
"raw_total",
"respondent_colonies")) %>%
mutate_at(vars("raw_total",
"respondent_colonies"),
as.numeric)
capa.melt.correlations <- melt(capa.df.filt,
id.var=c("year",
"province",
"winterloss_percent",
"raw_total",
"respondent_colonies")) %>%
filter(!is.na(winterloss_percent)) %>%
mutate_at(vars(winterloss_percent,
"raw_total",
"respondent_colonies"),
as.numeric) %>%
dplyr::rename("winterloss" = winterloss_percent)Figure 1
Figure 1A
capa.melt.g <- capa.melt %>%
filter(grepl("respondent_colonies_percent",
variable)) %>%
filter(!is.na(value)) %>%
filter(!grepl("TOTAL",
province)) %>%
mutate_at(vars(value),
as.numeric) %>%
arrange(province) %>%
group_by(province) %>%
mutate_at(vars(value),
funs(mean(.))) %>%
ungroup() %>%
droplevels() %>%
distinct(province,
.keep_all=TRUE)
#Get province polygons
provinces <- rnaturalearth::ne_states(country = "canada",
returnclass = "sf") %>%
filter(!grepl("Yukon|Nunavut|North",
name)) %>%
st_simplify(.,
dTolerance = 5000) %>% #simplify edges for plotting purposes
st_transform(., crs=3348) %>% #for 3348 info, see: https://epsg.io/3348
mutate(colonies_represented = dplyr::recode(!!!setNames(capa.melt.g$value,
paste0("CA-",
capa.melt.g$province)),
iso_3166_2))
#Get central points for each province
provinces.points <- provinces %>%
as.data.frame() %>%
st_as_sf(coords = c("longitude",
"latitude"),
crs=4326) %>%
st_transform(crs=3348)
#Plot
fig1.top.1 <- ggplot() +
geom_sf(data = provinces,
aes(fill = colonies_represented),
linewidth=0.05,
alpha=1,
colour="white") +
geom_sf(data = provinces.points,
colour="black",
size=2) +
ggtitle("CAPA Survey Overview (2015-2023)") +
theme_void() +
ggsflabel::geom_sf_label_repel(data= provinces.points,
aes(label=name),
nudge_x = c(-20, #BC
20, #AB
-20, #SK
15, #MB
-25, #ON
-50, #QC
-1.5, #NB
50, #NL
20, #NS
20)*10000,
nudge_y = c(-90, #BC
70, #AB
-75, #SK
75, #MB
-60, #ON
40, #QC
-70, #NB
50, #NL
-40, #NS
20)*10000) +
theme(panel.grid = element_blank(),
plot.title=element_text(size=24,
face="bold",
hjust = 0.5,
vjust=-2),
legend.position = "none")
fig1.top.2 <- ggplot() +
geom_sf(data = provinces,
aes(fill = colonies_represented),
linewidth=0.05,
alpha=1,
colour="white") +
geom_sf(data = provinces.points,
colour="black",
size=2) +
ggtitle("CAPA Survey Participation") +
scale_fill_gradientn(colours = c("#132B43",
"#56B1F7"),
limits = c(40, 100),
breaks = c(40, 60, 80, 100),
labels = c("40", "60", "80", "100"),
guide = guide_colourbar(direction = "horizontal",
title = "% participation",
title.position = "top")) +
theme_void() +
theme(panel.grid = element_blank(),
plot.title=element_text(size=24,
face="bold",
hjust = 0.5,
vjust=-2),
legend.position="bottom")
fig1.top.legend <- cowplot::plot_grid(NULL,
ggplotGrob(
fig1.top.2
)$grobs[[
which(sapply(ggplotGrob(fig1.top.2)$grobs,
function(x) x$name) == "guide-box")
]],
ncol=2,
rel_widths=c(0.95, 0.05))
fig1.a <- cowplot::plot_grid(fig1.top.legend +
theme(plot.margin = margin(4,3.25,-5,-3, "cm")),
fig1.top.1,
ncol=1,
rel_heights=c(0.2,0.8))
print(fig1.a)
Figure 1B
capa.melt.g <- capa.melt %>%
filter(grepl("respondent_colonies_percent",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
filter(grepl("TOTAL",
province)) %>%
droplevels() %>%
mutate(rel_respondent = (respondent_colonies/raw_total)*100)
capa.melt.g <- rbind(capa.melt.g %>%
mutate(new_respondents = respondent_colonies) %>%
mutate(rel_group="a_Colonies surveyed"),
capa.melt.g %>%
mutate(new_respondents = raw_total-respondent_colonies) %>%
mutate(rel_group="b_Colonies not surveyed"))
#All of Canada
bar.full <- ggplot(capa.melt.g,
aes(x=year,
y=new_respondents)) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
stat="identity",
position=position_dodge()) +
ylab("Total no. colonies") +
scale_x_continuous(breaks=unique(capa.melt.g$year)) +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50")) +
theme(axis.text.x = element_text(size = 8,
angle=45,
hjust=1),
axis.title = element_text(size=10),
axis.title.x = element_blank(),
axis.text.y = element_text(size=8),
legend.position = "none")
#---Pie graph 2015---
pie.2015 <- ggplot(capa.melt.g %>%
filter(year==2015),
aes(x=year,
y=(new_respondents))) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
linewidth=0.7,
stat="identity") +
coord_polar("y", start=0) +
theme_void() +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50")) +
theme(legend.position = "none")
#---Pie graph 2016---
pie.2016 <- ggplot(capa.melt.g %>%
filter(year==2016),
aes(x=year,
y=(new_respondents))) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
linewidth=0.7,
stat="identity") +
coord_polar("y", start=0) +
theme_void() +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50")) +
theme(legend.position = "none")
#---Pie graph 2017---
pie.2017 <- ggplot(capa.melt.g %>%
filter(year==2017),
aes(x=year,
y=(new_respondents))) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
linewidth=0.7,
stat="identity") +
coord_polar("y", start=0) +
theme_void() +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50")) +
theme(legend.position = "none")
#---Pie graph 2018---
pie.2018 <- ggplot(capa.melt.g %>%
filter(year==2018),
aes(x=year,
y=(new_respondents))) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
linewidth=0.7,
stat="identity") +
coord_polar("y", start=0) +
theme_void() +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50")) +
theme(legend.position = "none")
#---Pie graph 2019---
pie.2019 <- ggplot(capa.melt.g %>%
filter(year==2019),
aes(x=year,
y=(new_respondents))) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
linewidth=0.7,
stat="identity") +
coord_polar("y", start=0) +
theme_void() +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50")) +
theme(legend.position = "none")
#---Pie graph 2020---
pie.2020 <- ggplot(capa.melt.g %>%
filter(year==2020),
aes(x=year,
y=(new_respondents))) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
linewidth=0.7,
stat="identity") +
coord_polar("y", start=0) +
theme_void() +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50")) +
theme(legend.position = "none")
#---Pie graph 2021---
pie.2021 <- ggplot(capa.melt.g %>%
filter(year==2021),
aes(x=year,
y=(new_respondents))) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
linewidth=0.7,
stat="identity") +
coord_polar("y", start=0) +
theme_void() +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50")) +
theme(legend.position = "none")
#---Pie graph 2022---
pie.2022 <- ggplot(capa.melt.g %>%
filter(year==2022),
aes(x=year,
y=(new_respondents))) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
linewidth=0.7,
stat="identity") +
coord_polar("y", start=0) +
theme_void() +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50")) +
theme(legend.position = "none")
#---Pie graph 2023---
pie.2023 <- ggplot(capa.melt.g %>%
filter(year==2023),
aes(x=year,
y=(new_respondents))) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
linewidth=0.7,
stat="identity") +
coord_polar("y", start=0) +
theme_void() +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50")) +
theme(legend.position = "none")
#---Pie legend---
pie.legend <- ggplot(capa.melt.g %>%
filter(year==2015),
aes(x=year,
y=(new_respondents))) +
geom_bar(aes(group=rel_group,
fill= rel_group),
colour="white",
stat="identity") +
coord_polar("y", start=0) +
theme_void() +
scale_fill_manual(values= c("a_Colonies surveyed" = "#3875A7",
"b_Colonies not surveyed" = "grey50"),
labels = c("Colonies surveyed",
"Colonies not surveyed")) +
theme(legend.position = "bottom",
legend.title=element_blank())
pie.g <- cowplot::plot_grid(pie.2015,
pie.2016,
pie.2017,
pie.2018,
pie.2019,
pie.2020,
pie.2021,
pie.2022,
pie.2023,
ncol=9)
pie.g.legend <- cowplot::plot_grid(ggplotGrob(
pie.legend
)$grobs[[
which(sapply(ggplotGrob(pie.legend)$grobs,
function(x) x$name) == "guide-box")
]])
pie.full <- cowplot::plot_grid(pie.g.legend +
theme(plot.margin = margin(0,0,0,0, "cm")),
pie.g,
ncol=1,
rel_heights=c(0.1,0.9))
fig1.b <- cowplot::plot_grid(cowplot::plot_grid(NULL,
pie.full +
theme(plot.margin = margin(0,0,0,0, "cm")),
NULL,
ncol=3,
rel_widths=c(0.15, 0.8, 0.05)) +
theme(plot.margin = margin(0.2,0,-0.5,0, "cm")),
bar.full,
ncol=1,
rel_heights=c(0.25,0.75))
print(fig1.b)
Figure 1C
fig1.c <- readRDS("RDS/fig1c.RDS")
print(fig1.c)
Figure 1D
capa.melt.g <- capa.melt %>%
filter(grepl("^winterloss",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
mutate(value = ifelse(is.na(value),
0,
value)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(line_width = ifelse(grepl("TOTAL",
province),
2, 0.1)) %>%
mutate(line_cols = ifelse(grepl("TOTAL",
province),
"#640000",
"grey50")) %>%
arrange(year) %>% arrange(province)
#Plot
fig1.d <- ggplot(capa.melt.g,
aes(x=year,
y=value,
group=province,
colour=province)) +
geom_line(aes(linewidth=province)) +
scale_colour_manual(values=setNames(capa.melt.g$line_cols,
capa.melt.g$province)) +
scale_linewidth_manual(values = setNames(capa.melt.g$line_width,
capa.melt.g$province)) +
directlabels::geom_dl(aes(label=province),
method=list("last.points",
cex=0.75,
rot =0,
hjust = -.5)) +
scale_x_continuous(breaks=unique(capa.melt.g$year)) +
ylab(unique(capa.melt.g$variable)) +
xlab("Year") +
ylab("Overwintering loss (%)") +
theme_classic() +
theme(axis.text.x = element_text(size = 8,
angle=45,
hjust=1),
axis.title = element_text(size=10),
axis.title.x = element_blank(),
axis.text.y = element_text(size=8),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x.bottom = element_text(size=10),
legend.position = "none")
fig1.d
Figure 1E
capa.melt.g <- capa.melt %>%
filter(grepl("^winterloss",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
filter(!grepl("TOTAL",
province)) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL")))
#Per province
fig1.e <- ggplot(capa.melt.g,
aes(x=value,
y=province,
fill=stat(x))) +
geom_density_ridges_gradient(scale = 1.7,
panel_scaling = TRUE,
rel_min_height = 0.0005,
colour="black",
linewidth = 0.3) +
scale_fill_gradientn(colours=colorRampPalette(c("white",
"darkred"))(100)) +
geom_vline(xintercept=15,
linewidth = 0.8,
linetype = "dotted",
colour="red") +
annotate(geom="text",
x = 15,
y = 12.9,
label="Limit of \nsustainability",
colour="red") +
coord_cartesian(ylim = c(1,10.5),
clip = 'off') +
xlab("Overwintering loss (%)") +
ylab("Province") +
theme(plot.margin = unit(c(2.5,1,1,1),
"lines")) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size = 8),
axis.title = element_text(size=10),
axis.title.x = element_blank(),
axis.text.y = element_text(size=8),
legend.position = "none")
fig1.e
Figure 1 (Full panel)
fig1.bottom.left <- cowplot::plot_grid(fig1.b,
fig1.c,
ncol=1,
labels=c('B', "C"),
rel_heights=c(0.45, 0.55))
fig1.bottom.right <- cowplot::plot_grid(fig1.e,
cowplot::plot_grid(fig1.d,
NULL,
ncol=2,
rel_widths=c(0.9,0.1)),
ncol=1,
labels=c('D', "E"),
rel_heights=c(0.55, 0.45))
fig1.bottom <- cowplot::plot_grid(fig1.bottom.left,
fig1.bottom.right, ncol=2)
fig1.a <- cowplot::plot_grid(fig1.top.legend +
theme(plot.margin = margin(3,3,-5,-3, "cm")),
fig1.top.1,
ncol=1,
rel_heights=c(0.2,0.8))
fig1.full <- cowplot::plot_grid(fig1.a +
theme(plot.margin = unit(c(-0.04, -0.04, -0.04, -0.04),
"null")),
fig1.bottom,
ncol=1,
rel_heights=c(0.5, 0.5))
print(fig1.full)
Figure 2
Setup
capa.df.filt <- capa.df %>%
mutate_at(vars(respondent_colonies_percent,
bac_spring_oxytet,
bac_fall_oxytet,
nos_spring_fum,
nos_fall_fum,
var_spring_other,
var_fall_other),
as.numeric) %>%
mutate(respondent_colonies_percent = ifelse(
respondent_colonies_percent > 100, 100,
respondent_colonies_percent)) %>%
mutate(oxytet_spring_fall = (bac_spring_oxytet +
bac_fall_oxytet)/2) %>%
mutate(fum_spring_fall = (nos_spring_fum +
nos_fall_fum)/2) %>%
mutate(var_spring_fall = (var_spring_other +
var_fall_other)/2) %>%
mutate(winterloss= winterloss_percent)
capa.melt <- melt(capa.df.filt ,
id.var=c("year",
"province",
"winterloss",
"raw_total",
"respondent_colonies")) %>%
mutate_at(vars("raw_total",
"respondent_colonies"),
as.numeric)
capa.melt.correlations <- melt(capa.df.filt,
id.var=c("year",
"province",
"winterloss_percent",
"raw_total",
"respondent_colonies")) %>%
filter(!is.na(winterloss_percent)) %>%
mutate_at(vars(winterloss_percent, "raw_total",
"respondent_colonies"),
as.numeric) %>%
dplyr::rename("winterloss" = winterloss_percent)Figure 2A
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("bac_spring_oxytet|bac_fall_oxytet",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
arrange(year) %>%
arrange(province) %>%
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="oxytet_spring_fall") %>%
{. ->> tmp} %>%
#---Calculate weighted totals
filter(!grepl("TOTAL|NL", province)) %>%
group_by(year, variable) %>%
mutate(value_weighted1 = value * raw_total) %>%
mutate(value_weighted2 = sum(value_weighted1, na.rm=TRUE)/sum(raw_total, na.rm=TRUE)) %>%
mutate(value=value_weighted2) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>% dplyr::select(-value_weighted1, -value_weighted2) %>%
#---Combine totals + all provinces
rbind(., subset(tmp, province!="TOTAL")) %>%
mutate(line_width = ifelse(grepl("TOTAL",
province),
1.5,
0.1)) %>%
mutate(line_cols = ifelse(grepl("TOTAL",
province),
"#2980B9",
"grey50"))
fig2.a <- ggplot(capa.melt.g,
aes(x=year,
y=value,
group=province,
colour=province)) +
geom_line(aes(linewidth=province)) +
ylim(-2,86) +
scale_y_continuous(limits=c(-2,110),
breaks=c(0,25,50,75,100)) +
scale_colour_manual(values=setNames(capa.melt.g$line_cols,
capa.melt.g$province)) +
scale_linewidth_manual(values = setNames(capa.melt.g$line_width,
capa.melt.g$province)) +
directlabels::geom_dl(aes(label=province),
method=list("last.points",
cex=0.5,
rot =0,
hjust = 0)) +
ylab(unique(capa.melt.g$variable)[1]) +
xlab("Year") +
theme(legend.position = "none") +
ylab("Oxytetracycline use (%)") +
scale_x_continuous(breaks = seq(2015, 2023, 1)) +
theme(axis.text.x = element_text(angle=45,
hjust=1),
axis.text = element_text(size=8),
axis.title = element_text(size=10),
axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
fig2.a
Figure 2B
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("nos_spring_fum|nos_fall_fum",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
arrange(year) %>%
arrange(province) %>%
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="fum_spring_fall") %>%
{. ->> tmp} %>%
#---Calculate weighted totals
filter(!grepl("TOTAL|NL", province)) %>%
group_by(year, variable) %>%
mutate(value_weighted1 = value * raw_total) %>%
mutate(value_weighted2 = sum(value_weighted1, na.rm=TRUE)/sum(raw_total, na.rm=TRUE)) %>%
mutate(value=value_weighted2) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>% dplyr::select(-value_weighted1, -value_weighted2) %>%
#---Combine totals + all provinces
rbind(., subset(tmp, province!="TOTAL")) %>%
mutate(line_width = ifelse(grepl("TOTAL",
province),
1.5, 0.1)) %>%
mutate(line_cols = ifelse(grepl("TOTAL",
province),
"green4",
"grey50"))
fig2.b <-ggplot(capa.melt.g,
aes(x=year,
y=value,
group=province,
colour=province)) +
geom_line(aes(linewidth=province)) +
scale_y_continuous(limits=c(-2,110),
breaks=c(0,25,50,75,100)) +
scale_colour_manual(values=setNames(capa.melt.g$line_cols,
capa.melt.g$province)) +
scale_linewidth_manual(values = setNames(capa.melt.g$line_width,
capa.melt.g$province)) +
directlabels::geom_dl(aes(label=province),
method=list("last.points",
cex=0.5,
rot =0,
hjust = 0)) +
ylab(unique(capa.melt.g$variable)[1]) +
xlab("Year") + theme(legend.position = "none") +
ylab("Fumagillin use (%)") +
scale_x_continuous(breaks = seq(2015, 2023, 1)) +
theme(axis.text.x = element_text(angle=45,
hjust=1),
axis.text = element_text(size=8), axis.title = element_text(size=10),
axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
fig2.b
Figure 2C
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("var_spring_other|var_fall_other",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
arrange(year) %>%
arrange(province) %>%
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="var_spring_fall") %>%
{. ->> tmp} %>%
#---Calculate weighted totals
filter(!grepl("TOTAL|NL", province)) %>%
group_by(year, variable) %>%
mutate(value_weighted1 = value * raw_total) %>%
mutate(value_weighted2 = sum(value_weighted1, na.rm=TRUE)/sum(raw_total, na.rm=TRUE)) %>%
mutate(value=value_weighted2) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>% dplyr::select(-value_weighted1, -value_weighted2) %>%
#---Combine totals + all provinces
rbind(., subset(tmp, province!="TOTAL")) %>%
mutate(line_width = ifelse(grepl("TOTAL",
province),
1.5, 0.1)) %>%
mutate(line_cols = ifelse(grepl("TOTAL",
province),
"#AC8300",
"grey50"))
fig2.c <- ggplot(capa.melt.g,
aes(x=year,
y=value,
group=province,
colour=province)) +
geom_line(aes(linewidth=province)) +
scale_y_continuous(limits=c(-2,110),
breaks=c(0,25,50,75,100)) +
scale_colour_manual(values=setNames(capa.melt.g$line_cols,
capa.melt.g$province)) +
scale_linewidth_manual(values = setNames(capa.melt.g$line_width,
capa.melt.g$province)) +
directlabels::geom_dl(aes(label=province),
method=list("last.points",
cex=0.5,
rot =0,
hjust = 0)) +
ylab(unique(capa.melt.g$variable)[1]) +
xlab("Year") +
theme(legend.position = "none") +
ylab("Miticide use (%)") +
scale_x_continuous(breaks = seq(2015, 2023, 1)) +
theme(axis.text.x = element_text(angle=45, hjust=1),
axis.text = element_text(size=8),
axis.title = element_text(size=10),
axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
fig2.c
Figure 2D
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("bac_spring_oxytet|bac_fall_oxytet",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
mutate(before_after = ifelse(year<=2018,
"Before",
"After")) %>%
mutate(before_after = factor(before_after,
levels=c("Before",
"After"))) %>%
mutate(line_width = ifelse(grepl("TOTAL",
province),
2, 0.1)) %>%
mutate(line_cols = ifelse(grepl("TOTAL",
province),
"#2980B9",
"grey50")) %>%
arrange(year) %>%
arrange(province) %>%
#---Calculate national totals
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="oxytet_spring_fall") %>%
#---Calculate weighted totals
filter(!grepl("TOTAL", province)) %>%
group_by(year, variable) %>%
mutate(value_weighted1 = value * raw_total) %>%
mutate(value_weighted2 = sum(value_weighted1, na.rm=TRUE)/sum(raw_total, na.rm=TRUE)) %>%
mutate(value=value_weighted2) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>%
mutate(province= gsub("TOTAL",
"WEIGHTED TOTAL",
province))
#-------------STATS
stat.oxytet.before.after <- capa.melt.g %>%
group_by(before_after) %>%
summarise(mean = round(mean(value,
na.rm=TRUE), 2),
std_deviation = round(sd(value,
na.rm=TRUE), 2),
std_error = round((std_deviation / sqrt(n())), 2)) %>%
arrange(mean)
print(stat.oxytet.before.after)## # A tibble: 2 × 4
## before_after mean std_deviation std_error
## <fct> <dbl> <dbl> <dbl>
## 1 After 29.5 5.45 2.44
## 2 Before 55.5 5.88 2.94readr::write_tsv(stat.oxytet.before.after, "tables/stats_oxytet_before_after.tsv")
#-------------PLOT
norm.test <- capa.melt.g %>%
rstatix::shapiro_test(var=value)
pwc1 <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0, 0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
group_by(province) %>%
rstatix::wilcox_test(value ~ before_after,
p.adjust.method = "BH",
detailed=TRUE) %>%
rstatix::add_xy_position() %>%
mutate(p_str = format(round(p, 4),
scientific = FALSE))
pwc1.ft <- flextable::flextable(pwc1)
pwc1.ft <- flextable::theme_vanilla(pwc1.ft)
pwc1.ftprovince | estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative | y.position | groups | xmin | xmax | p_str |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
WEIGHTED TOTAL | 27.91386 | value | Before | After | 4 | 5 | 20 | 0.0159 | 16.26848 | 34.80699 | Wilcoxon | two.sided | 62.22496 | [[character]] | 1 | 2 | 0.0159 |
#Remove cases where there is no variance (i.e. all value equal zero or represent 'NA' observations)
summary.stats <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0, 0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
group_by(before_after) %>%
rstatix::get_summary_stats(value,
type = c("full"))
ft <- flextable::flextable(summary.stats)
ft <- flextable::theme_vanilla(ft)
ftbefore_after | variable | n | min | max | median | q1 | q3 | iqr | mad | mean | sd | se | ci |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Before | value | 4 | 46.910 | 59.638 | 57.780 | 54.213 | 59.095 | 4.882 | 2.217 | 55.527 | 5.884 | 2.942 | 9.363 |
After | value | 5 | 24.107 | 38.080 | 29.091 | 25.694 | 30.641 | 4.947 | 5.037 | 29.523 | 5.446 | 2.436 | 6.763 |
fig2.d <- ggplot(capa.melt.g,
aes(x=before_after,
y=value)) +
scale_y_continuous(limits=c(-2,110),
breaks=c(0,25,50,75,100)) +
facet_grid(~province) +
ylab("Oxytetracycline use (%)") +
geom_boxplot(aes(),
fill="#2980B9",
alpha=0.5,
outliers=FALSE) +
geom_point(position=position_jitter(width=0.05),
size=2,
alpha=0.5) +
ggpubr::stat_pvalue_manual(pwc1,
label="p_str",
hide.ns=FALSE,
size=2.5,
bracket.size=0.1,
vjust=-0.5,
step.increase=0.0000000000001) +
theme(axis.text.x = element_text(size = 10,
angle=90,
vjust=0.5),
axis.title = element_text(size=10),
axis.text.y = element_text(size=8),
axis.title.x = element_blank(),
strip.background = element_rect(fill=scales::alpha("#2980B9",
0.5)),
strip.text = element_text(colour = 'black',
size=8,
face="bold"))
fig2.d
Figure 2E
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("bac_spring_oxytet|bac_fall_oxytet",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
mutate(before_after = ifelse(year<=2018,
"Before",
"After")) %>%
mutate(before_after = factor(before_after,
levels=c("Before",
"After"))) %>%
mutate(line_width = ifelse(grepl("TOTAL",
province),
2, 0.1)) %>%
mutate(line_cols = ifelse(grepl("TOTAL",
province),
"#2980B9",
"grey50")) %>%
arrange(year) %>%
arrange(province) %>%
#---Calculate annual totals for each province
group_by(year, province) %>%
mutate(value=mean(value, na.rm=TRUE)) %>%
slice(1) %>% ungroup() %>%
mutate(variable="oxytet_spring_fall") %>%
filter(!grepl("TOTAL", province))
#-------------STATS
stat.oxytet.before.after <- capa.melt.g %>%
group_by(before_after) %>%
summarise(mean = round(mean(value,
na.rm=TRUE), 2),
std_deviation = round(sd(value,
na.rm=TRUE), 2),
std_error = round((std_deviation / sqrt(n())), 2)) %>%
arrange(mean)
stat.oxytet.before.after## # A tibble: 2 × 4
## before_after mean std_deviation std_error
## <fct> <dbl> <dbl> <dbl>
## 1 After 24.4 20.7 2.93
## 2 Before 42.6 28.6 4.53readr::write_tsv(stat.oxytet.before.after, "tables/stats_oxytet_before_after.tsv")
#-------------PLOT
norm.test <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0, 0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
rstatix::shapiro_test(var=value)
pwc1 <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0, 0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
group_by(province) %>%
rstatix::wilcox_test(value ~before_after,
p.adjust.method = "BH",
detailed=TRUE) %>%
rstatix::add_xy_position() %>%
mutate(p_str = format(round(p, 4),
scientific = FALSE))
pwc1.ft <- flextable::flextable(pwc1)
pwc1.ft <- flextable::theme_vanilla(pwc1.ft)
pwc1.ftprovince | estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative | y.position | groups | xmin | xmax | p_str |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
BC | 4.50000 | value | Before | After | 4 | 5 | 19 | 0.0317 | 0.5000000 | 17.500000 | Wilcoxon | two.sided | 33.22 | [[character]] | 1 | 2 | 0.0317 |
AB | 37.25000 | value | Before | After | 4 | 5 | 20 | 0.0159 | 11.0000000 | 54.500000 | Wilcoxon | two.sided | 79.22 | [[character]] | 1 | 2 | 0.0159 |
SK | 34.50000 | value | Before | After | 4 | 5 | 20 | 0.0159 | 14.0000000 | 47.500000 | Wilcoxon | two.sided | 95.72 | [[character]] | 1 | 2 | 0.0159 |
MB | 35.75000 | value | Before | After | 4 | 5 | 20 | 0.0159 | 15.5000000 | 48.900000 | Wilcoxon | two.sided | 83.22 | [[character]] | 1 | 2 | 0.0159 |
ON | 8.22500 | value | Before | After | 4 | 5 | 17 | 0.1110 | -2.5000000 | 17.500000 | Wilcoxon | two.sided | 81.67 | [[character]] | 1 | 2 | 0.1110 |
QC | 3.49991 | value | Before | After | 4 | 5 | 20 | 0.0195 | 2.0000088 | 4.500028 | Wilcoxon | two.sided | 17.72 | [[character]] | 1 | 2 | 0.0195 |
NB | 15.50835 | value | Before | After | 4 | 5 | 20 | 0.0195 | 4.0000420 | 24.000008 | Wilcoxon | two.sided | 68.22 | [[character]] | 1 | 2 | 0.0195 |
PE | 12.44553 | value | Before | After | 4 | 5 | 18 | 0.0651 | -0.7999496 | 31.499973 | Wilcoxon | two.sided | 43.72 | [[character]] | 1 | 2 | 0.0651 |
NS | 48.07500 | value | Before | After | 4 | 5 | 19 | 0.0317 | 10.1500000 | 56.500000 | Wilcoxon | two.sided | 79.22 | [[character]] | 1 | 2 | 0.0317 |
#Remove cases where there is no variance (i.e. all value equal zero or represent 'NA' observations)
summary.stats <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0, 0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
group_by(province,
before_after) %>%
rstatix::get_summary_stats(value,
type = c("full"))
ft <- flextable::flextable(summary.stats)
ft <- flextable::theme_vanilla(ft)
ftprovince | before_after | variable | n | min | max | median | q1 | q3 | iqr | mad | mean | sd | se | ci |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
BC | Before | value | 4 | 9.0 | 23.50 | 11.500 | 10.500 | 14.875 | 4.375 | 2.224 | 13.875 | 6.537 | 3.268 | 10.401 |
BC | After | value | 5 | 5.0 | 9.50 | 7.000 | 6.000 | 8.500 | 2.500 | 2.224 | 7.200 | 1.823 | 0.815 | 2.264 |
AB | Before | value | 4 | 36.5 | 69.50 | 58.000 | 48.875 | 64.625 | 15.750 | 12.231 | 55.500 | 14.370 | 7.185 | 22.866 |
AB | After | value | 5 | 14.0 | 32.50 | 24.000 | 15.000 | 25.500 | 10.500 | 12.602 | 22.200 | 7.735 | 3.459 | 9.604 |
SK | Before | value | 4 | 68.0 | 86.00 | 75.750 | 73.250 | 78.875 | 5.625 | 6.301 | 76.375 | 7.409 | 3.705 | 11.790 |
SK | After | value | 5 | 33.5 | 61.00 | 40.500 | 38.500 | 50.000 | 11.500 | 10.378 | 44.700 | 10.901 | 4.875 | 13.535 |
MB | Before | value | 4 | 60.0 | 73.50 | 70.450 | 66.375 | 72.675 | 6.300 | 3.706 | 68.600 | 6.122 | 3.061 | 9.741 |
MB | After | value | 5 | 23.5 | 53.00 | 31.500 | 29.500 | 38.500 | 9.000 | 10.378 | 35.200 | 11.300 | 5.054 | 14.031 |
ON | Before | value | 4 | 62.0 | 71.95 | 67.750 | 65.000 | 70.112 | 5.112 | 4.411 | 67.362 | 4.329 | 2.165 | 6.889 |
ON | After | value | 5 | 52.0 | 68.00 | 57.000 | 56.500 | 64.500 | 8.000 | 7.413 | 59.600 | 6.494 | 2.904 | 8.064 |
QC | Before | value | 4 | 6.5 | 8.00 | 7.500 | 7.250 | 7.625 | 0.375 | 0.371 | 7.375 | 0.629 | 0.315 | 1.001 |
QC | After | value | 5 | 3.0 | 5.00 | 4.000 | 3.500 | 4.500 | 1.000 | 0.741 | 4.000 | 0.791 | 0.354 | 0.982 |
NB | Before | value | 4 | 45.5 | 58.50 | 56.500 | 53.750 | 57.000 | 3.250 | 1.483 | 54.250 | 5.909 | 2.955 | 9.403 |
NB | After | value | 5 | 32.5 | 44.00 | 40.500 | 38.000 | 41.500 | 3.500 | 3.706 | 39.300 | 4.367 | 1.953 | 5.423 |
PE | Before | value | 4 | 10.0 | 34.00 | 18.100 | 13.525 | 24.625 | 11.100 | 8.525 | 20.050 | 10.430 | 5.215 | 16.596 |
PE | After | value | 5 | 0.0 | 15.50 | 9.000 | 2.500 | 9.000 | 6.500 | 9.637 | 7.200 | 6.109 | 2.732 | 7.586 |
NS | Before | value | 4 | 45.0 | 69.50 | 67.825 | 61.988 | 68.375 | 6.387 | 1.371 | 62.538 | 11.719 | 5.860 | 18.648 |
NS | After | value | 5 | 12.5 | 57.50 | 19.000 | 13.000 | 20.500 | 7.500 | 8.896 | 24.500 | 18.785 | 8.401 | 23.325 |
fig2.e <- ggplot(capa.melt.g,
aes(x=before_after,
y=value)) +
facet_grid(~province) +
scale_y_continuous(limits=c(-2,110),
breaks=c(0,25,50,75,100)) +
geom_boxplot(aes(),
fill="#2980B9",
alpha=0.5,
outliers=FALSE) +
geom_point(position=position_jitter(width=0.15),
size=1,
alpha=0.5) +
ggpubr::stat_pvalue_manual(pwc1,
label="p_str",
hide.ns=FALSE,
size=2.5,
bracket.size=0.1,
vjust=-0.5,
step.increase=0.0000000000001,
bracket.nudge.y=-7,
tip.length = 0.015) +
theme(axis.text.x = element_text(size = 10,
angle=90,
vjust = 0.5),
axis.title = element_text(size=10),
axis.text.y = element_text(size=8),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
strip.background = element_rect(fill=scales::alpha("#2980B9", 0.5)),
strip.text = element_text(colour = 'black', size=8, face="bold"))
fig2.e
Figure 2F
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("nos_spring_fum|nos_fall_fum",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province = factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
mutate(before_after = ifelse(year<=2018,
"Before",
"After")) %>%
mutate(before_after = factor(before_after,
levels=c("Before",
"After"))) %>%
mutate(line_width = ifelse(grepl("TOTAL",
province),
2, 0.1)) %>%
mutate(line_cols = ifelse(grepl("TOTAL",
province),
"green4",
"grey50")) %>%
arrange(year) %>%
arrange(province) %>%
#---Calculate national totals
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="fum_spring_fall") %>%
#---Calculate weighted totals
filter(!grepl("TOTAL", province)) %>%
group_by(year, variable) %>%
mutate(value_weighted1 = value * raw_total) %>%
mutate(value_weighted2 = sum(value_weighted1, na.rm=TRUE)/sum(raw_total, na.rm=TRUE)) %>%
mutate(value=value_weighted2) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>%
mutate(province= gsub("TOTAL",
"WEIGHTED TOTAL",
province))
#-------------STATS
stat.fum.before.after <- capa.melt.g %>%
group_by(before_after) %>%
summarise(mean = round(mean(value,
na.rm=TRUE), 2),
std_deviation = round(sd(value,
na.rm=TRUE), 2),
std_error = round((std_deviation / sqrt(n())), 2)) %>%
arrange(mean)
stat.fum.before.after## # A tibble: 2 × 4
## before_after mean std_deviation std_error
## <fct> <dbl> <dbl> <dbl>
## 1 After 25.3 3.69 1.65
## 2 Before 52.9 4.75 2.38readr::write_tsv(stat.fum.before.after, "tables/stats_oxytet_before_after.tsv")
#-------------PLOT
norm.test <- capa.melt.g %>%
rstatix::shapiro_test(var=value)
pwc1 <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0, 0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
group_by(province) %>%
rstatix::wilcox_test(value ~ before_after,
p.adjust.method = "BH",
detailed=TRUE) %>%
rstatix::add_xy_position() %>%
mutate(p_str = format(round(p, 4),
scientific = FALSE))
pwc1.ft <- flextable::flextable(pwc1)
pwc1.ft <- flextable::theme_vanilla(pwc1.ft)
pwc1.ftprovince | estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative | y.position | groups | xmin | xmax | p_str |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
WEIGHTED TOTAL | 26.50413 | value | Before | After | 4 | 5 | 20 | 0.0159 | 20.15936 | 34.31677 | Wilcoxon | two.sided | 62.7818 | [[character]] | 1 | 2 | 0.0159 |
#Remove cases where there is no variance (i.e. all value equal zero or represent 'NA' observations)
summary.stats <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0, 0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
group_by(before_after) %>%
rstatix::get_summary_stats(value,
type = c("full"))
ft <- flextable::flextable(summary.stats)
ft <- flextable::theme_vanilla(ft)
ftbefore_after | variable | n | min | max | median | q1 | q3 | iqr | mad | mean | sd | se | ci |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Before | value | 4 | 47.479 | 59.069 | 52.497 | 51.159 | 54.224 | 3.065 | 3.803 | 52.886 | 4.754 | 2.377 | 7.565 |
After | value | 5 | 19.109 | 28.129 | 27.233 | 24.753 | 27.320 | 2.567 | 1.328 | 25.309 | 3.690 | 1.650 | 4.581 |
fig2.f <- ggplot(capa.melt.g,
aes(x=before_after,
y=value)) +
scale_y_continuous(limits=c(-2,110),
breaks=c(0,25,50,75,100)) +
ylab("Fumagillin use (%)") +
facet_grid(~province) +
geom_boxplot(aes(),
fill="green4",
alpha=0.5,
outliers=FALSE) +
geom_point(position=position_jitter(width=0.05),
size=2,
alpha=0.5) +
ggpubr::stat_pvalue_manual(pwc1,
label="p_str",
hide.ns=FALSE,
size=2.5,
bracket.size=0.1,
vjust=-0.5,
step.increase=0.0000000000001) +
theme(axis.text.x = element_text(size = 10,
angle=90,
vjust=0.5),
axis.title = element_text(size=10),
axis.text.y = element_text(size=8),
axis.title.x = element_blank(),
strip.background =element_rect(fill=scales::alpha("green4", 0.5)),
strip.text = element_text(colour = 'black', size=8, face="bold"))
fig2.f
Figure 2G
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("nos_spring_fum|nos_fall_fum",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
mutate(before_after = ifelse(year<=2018,
"Before",
"After")) %>%
mutate(before_after = factor(before_after,
levels=c("Before",
"After"))) %>%
mutate(line_width = ifelse(grepl("TOTAL",
province),
2, 0.1)) %>%
mutate(line_cols = ifelse(grepl("TOTAL",
province),
"green4",
"grey50")) %>%
arrange(year) %>%
arrange(province) %>%
#---Calculate annual totals for each province
group_by(year, province) %>%
mutate(value=mean(value, na.rm=TRUE)) %>%
slice(1) %>% ungroup() %>%
mutate(variable="fum_spring_fall") %>%
filter(!grepl("TOTAL", province))
#-------------STATS
stat.fum.before.after <- capa.melt.g %>%
group_by(before_after,
province) %>%
summarise(mean = round(mean(value,
na.rm=TRUE), 2),
std_deviation = round(sd(value,
na.rm=TRUE), 2),
std_error = round((std_deviation / sqrt(n())), 2)) %>%
arrange(mean)
stat.fum.before.after## # A tibble: 20 × 5
## # Groups: before_after [2]
## before_after province mean std_deviation std_error
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 After NL 0 0 0
## 2 Before NL 1.25 2.5 1.25
## 3 After QC 2.6 1.56 0.7
## 4 After ON 7.6 2.88 1.29
## 5 After PE 9.2 7.39 3.3
## 6 Before QC 11.1 3.97 1.99
## 7 After MB 12.1 6.69 2.99
## 8 After BC 14.2 2.33 1.04
## 9 After NB 18.9 4.55 2.03
## 10 Before ON 20.1 3.65 1.82
## 11 Before PE 22.7 9.34 4.67
## 12 After NS 23.1 9.15 4.09
## 13 After SK 28.4 5.45 2.44
## 14 Before BC 32.9 14.5 7.26
## 15 Before MB 34.2 2.48 1.24
## 16 After AB 42.6 8.18 3.66
## 17 Before NB 46.8 11.0 5.5
## 18 Before SK 53.2 17.7 8.87
## 19 Before NS 54.8 6.16 3.08
## 20 Before AB 81.4 13.5 6.76readr::write_tsv(stat.fum.before.after, "tables/stats_oxytet_before_after.tsv")
#-------------PLOT
norm.test <- capa.melt.g %>%
rstatix::shapiro_test(var=value)
pwc1 <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0, 0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
group_by(province) %>%
rstatix::wilcox_test(value ~ before_after,
p.adjust.method = "BH",
detailed=TRUE) %>%
rstatix::add_xy_position() %>%
mutate(p_str = format(round(p, 4),
scientific = FALSE))
pwc1.ft <- flextable::flextable(pwc1)
pwc1.ft <- flextable::theme_vanilla(pwc1.ft)
pwc1.ftprovince | estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative | y.position | groups | xmin | xmax | p_str |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
BC | 14.85000 | value | Before | After | 4 | 5 | 20.0 | 0.0179 | 5.499977 | 38.50005 | Wilcoxon | two.sided | 64.22 | [[character]] | 1 | 2 | 0.0179 |
AB | 39.90000 | value | Before | After | 4 | 5 | 20.0 | 0.0159 | 14.000000 | 59.00000 | Wilcoxon | two.sided | 104.72 | [[character]] | 1 | 2 | 0.0159 |
SK | 24.75000 | value | Before | After | 4 | 5 | 19.0 | 0.0317 | 2.000000 | 48.50000 | Wilcoxon | two.sided | 86.22 | [[character]] | 1 | 2 | 0.0317 |
MB | 23.50000 | value | Before | After | 4 | 5 | 20.0 | 0.0159 | 10.650000 | 30.00000 | Wilcoxon | two.sided | 48.72 | [[character]] | 1 | 2 | 0.0159 |
ON | 12.50004 | value | Before | After | 4 | 5 | 20.0 | 0.0195 | 6.500019 | 19.50006 | Wilcoxon | two.sided | 36.72 | [[character]] | 1 | 2 | 0.0195 |
QC | 6.75000 | value | Before | After | 4 | 5 | 20.0 | 0.0159 | 5.000000 | 14.50000 | Wilcoxon | two.sided | 28.22 | [[character]] | 1 | 2 | 0.0159 |
NB | 30.18815 | value | Before | After | 4 | 5 | 20.0 | 0.0195 | 8.500054 | 40.49999 | Wilcoxon | two.sided | 65.72 | [[character]] | 1 | 2 | 0.0195 |
PE | 13.20000 | value | Before | After | 4 | 5 | 18.0 | 0.0635 | -1.100000 | 33.00000 | Wilcoxon | two.sided | 47.22 | [[character]] | 1 | 2 | 0.0635 |
NS | 32.20000 | value | Before | After | 4 | 5 | 20.0 | 0.0159 | 17.500000 | 44.00000 | Wilcoxon | two.sided | 73.22 | [[character]] | 1 | 2 | 0.0159 |
NL | 0.00000 | value | Before | After | 4 | 5 | 12.5 | 0.3710 | 0.000000 | 5.00000 | Wilcoxon | two.sided | 16.22 | [[character]] | 1 | 2 | 0.3710 |
summary.stats <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0, 0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
group_by(province,
before_after) %>%
rstatix::get_summary_stats(value,
type = c("full"))
ft <- flextable::flextable(summary.stats)
ft <- flextable::theme_vanilla(ft)
ftprovince | before_after | variable | n | min | max | median | q1 | q3 | iqr | mad | mean | sd | se | ci |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
BC | Before | value | 4 | 20.0 | 53.0 | 29.250 | 23.750 | 38.375 | 14.625 | 10.008 | 32.875 | 14.528 | 7.264 | 23.117 |
BC | After | value | 5 | 10.5 | 17.0 | 14.500 | 14.500 | 14.500 | 0.000 | 0.000 | 14.200 | 2.335 | 1.044 | 2.899 |
AB | Before | value | 4 | 62.5 | 93.5 | 84.900 | 76.975 | 89.375 | 12.400 | 8.673 | 81.450 | 13.507 | 6.754 | 21.493 |
AB | After | value | 5 | 29.0 | 49.0 | 45.000 | 41.500 | 48.500 | 7.000 | 5.189 | 42.600 | 8.181 | 3.659 | 10.158 |
SK | Before | value | 4 | 32.0 | 75.0 | 53.000 | 45.500 | 60.750 | 15.250 | 17.791 | 53.250 | 17.727 | 8.864 | 28.208 |
SK | After | value | 5 | 21.0 | 36.0 | 28.500 | 26.500 | 30.000 | 3.500 | 2.965 | 28.400 | 5.447 | 2.436 | 6.764 |
MB | Before | value | 4 | 31.5 | 37.5 | 33.825 | 33.112 | 34.875 | 1.763 | 1.853 | 34.163 | 2.484 | 1.242 | 3.953 |
MB | After | value | 5 | 6.0 | 23.0 | 11.000 | 7.500 | 13.000 | 5.500 | 5.189 | 12.100 | 6.693 | 2.993 | 8.311 |
ON | Before | value | 4 | 17.5 | 25.5 | 18.675 | 18.250 | 20.513 | 2.263 | 1.001 | 20.087 | 3.653 | 1.827 | 5.813 |
ON | After | value | 5 | 5.0 | 12.0 | 6.000 | 6.000 | 9.000 | 3.000 | 1.483 | 7.600 | 2.881 | 1.288 | 3.577 |
QC | Before | value | 4 | 8.5 | 17.0 | 9.500 | 8.875 | 11.750 | 2.875 | 1.112 | 11.125 | 3.966 | 1.983 | 6.311 |
QC | After | value | 5 | 0.0 | 4.0 | 3.000 | 2.500 | 3.500 | 1.000 | 0.741 | 2.600 | 1.557 | 0.696 | 1.934 |
NB | Before | value | 4 | 30.5 | 54.5 | 51.000 | 45.125 | 52.625 | 7.500 | 3.336 | 46.750 | 10.989 | 5.494 | 17.485 |
NB | After | value | 5 | 11.5 | 22.0 | 21.500 | 17.500 | 22.000 | 4.500 | 0.741 | 18.900 | 4.547 | 2.033 | 5.646 |
PE | Before | value | 4 | 16.0 | 36.0 | 19.450 | 16.300 | 25.875 | 9.575 | 4.818 | 22.725 | 9.336 | 4.668 | 14.856 |
PE | After | value | 5 | 0.0 | 17.5 | 12.000 | 3.000 | 13.500 | 10.500 | 8.154 | 9.200 | 7.387 | 3.304 | 9.173 |
NS | Before | value | 4 | 47.0 | 62.0 | 55.200 | 52.625 | 57.425 | 4.800 | 5.560 | 54.850 | 6.164 | 3.082 | 9.808 |
NS | After | value | 5 | 12.5 | 37.0 | 23.000 | 18.000 | 25.000 | 7.000 | 7.413 | 23.100 | 9.154 | 4.094 | 11.366 |
NL | Before | value | 4 | 0.0 | 5.0 | 0.000 | 0.000 | 1.250 | 1.250 | 0.000 | 1.250 | 2.500 | 1.250 | 3.978 |
NL | After | value | 5 | 0.0 | 0.0 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
fig2.g <- ggplot(capa.melt.g,
aes(x=before_after,
y=value)) +
facet_grid(~province) +
scale_y_continuous(limits=c(-2,110),
breaks=c(0,25,50,75,100)) +
geom_boxplot(aes(),
fill="green4",
alpha=0.5,
outliers=FALSE) +
geom_point(position=position_jitter(width=0.15),
size=1,
alpha=0.5) +
ggpubr::stat_pvalue_manual(pwc1,
label="p_str",
hide.ns=FALSE,
size=2.5, bracket.size=0.1,
vjust=-0.5,
step.increase=0.0000000000001,
bracket.nudge.y=-7,
tip.length = 0.015) +
theme(axis.text.x = element_text(size = 10,
angle=90,
vjust = 0.5),
axis.title = element_text(size=10),
axis.text.y = element_text(size=8),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
strip.background =element_rect(fill=scales::alpha("green4", 0.5)),
strip.text = element_text(colour = 'black',
size=8,
face="bold"))
fig2.g
Figure 2H
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("var_spring_other|var_fall_other",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
mutate(before_after = ifelse(year<=2018,
"Before",
"After")) %>%
mutate(before_after = factor(before_after,
levels=c("Before",
"After"))) %>%
mutate(line_width = ifelse(grepl("TOTAL",
province),
2, 0.1)) %>%
mutate(line_cols = ifelse(grepl("TOTAL",
province),
"#AC8300",
"grey50")) %>%
arrange(year) %>%
arrange(province) %>%
#---Calculate national totals
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="var_spring_fall") %>%
#---Calculate weighted totals
filter(!grepl("TOTAL", province)) %>%
group_by(year, variable) %>%
mutate(value_weighted1 = value * raw_total) %>%
mutate(value_weighted2 = sum(value_weighted1, na.rm=TRUE)/sum(raw_total, na.rm=TRUE)) %>%
mutate(value=value_weighted2) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>%
mutate(province= gsub("TOTAL",
"WEIGHTED TOTAL",
province))
#-------------STATS
stat.var.before.after <- capa.melt.g %>%
group_by(before_after) %>%
summarise(mean = round(mean(value,
na.rm=TRUE), 2),
std_deviation = round(sd(value,
na.rm=TRUE), 2),
std_error = round((std_deviation / sqrt(n())), 2)) %>%
arrange(mean)
stat.var.before.after## # A tibble: 2 × 4
## before_after mean std_deviation std_error
## <fct> <dbl> <dbl> <dbl>
## 1 Before 75.8 5.91 2.96
## 2 After 82.7 6.79 3.04readr::write_tsv(stat.var.before.after, "tables/stats_oxytet_before_after.tsv")
#-------------PLOT
norm.test <- capa.melt.g %>%
rstatix::shapiro_test(var=value)
pwc1 <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0,
0, 1)) %>% ungroup() %>%
filter(var!=0) %>%
group_by(province) %>%
rstatix::wilcox_test(value ~before_after,
p.adjust.method = "BH",
detailed=TRUE) %>%
rstatix::add_xy_position() %>%
mutate(p_str = format(round(p, 4),
scientific = FALSE))
pwc1.ft <- flextable::flextable(pwc1)
pwc1.ft <- flextable::theme_vanilla(pwc1.ft)
pwc1.ftprovince | estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative | y.position | groups | xmin | xmax | p_str |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
WEIGHTED TOTAL | -7.655427 | value | Before | After | 4 | 5 | 4 | 0.19 | -18.00971 | 5.745909 | Wilcoxon | two.sided | 90.89732 | [[character]] | 1 | 2 | 0.19 |
#Remove cases where there is no variance (i.e. all value equal zero or represent 'NA' observations)
summary.stats <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0,
0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
group_by(before_after) %>%
rstatix::get_summary_stats(value,
type = c("full"))
ft <- flextable::flextable(summary.stats)
ft <- flextable::theme_vanilla(ft)
ftbefore_after | variable | n | min | max | median | q1 | q3 | iqr | mad | mean | sd | se | ci |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Before | value | 4 | 68.774 | 82.765 | 75.781 | 72.654 | 78.903 | 6.249 | 6.537 | 75.775 | 5.905 | 2.953 | 9.396 |
After | value | 5 | 74.154 | 90.026 | 85.665 | 77.020 | 86.784 | 9.764 | 6.466 | 82.730 | 6.791 | 3.037 | 8.432 |
fig2.h <- ggplot(capa.melt.g,
aes(x=before_after,
y=value)) +
scale_y_continuous(limits=c(-2,110),
breaks=c(0,25,50,75,100)) +
ylab("Miticide use (%)") +
facet_grid(~province) +
geom_boxplot(aes(),
fill="#AC8300",
alpha=0.5,
outliers=FALSE) +
geom_point(position=position_jitter(width=0.05),
size=2, alpha=0.5) +
ggpubr::stat_pvalue_manual(pwc1,
label="p_str",
hide.ns=FALSE,
size=2.5,
bracket.size=0.1,
vjust=-0.5,
step.increase=0.0000000000001) +
theme(axis.text.x = element_text(size = 10,
angle=90,
vjust=0.5),
axis.title = element_text(size=10),
axis.text.y = element_text(size=8),
axis.title.x = element_blank(),
strip.background =element_rect(fill=scales::alpha("#AC8300", 0.5)),
strip.text = element_text(colour = 'black',
size=8,
face="bold"))
fig2.h
Figure 2I
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("var_spring_other|var_fall_other",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>% mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
mutate(before_after = ifelse(year<=2018,
"Before",
"After")) %>%
mutate(before_after = factor(before_after,
levels=c("Before",
"After"))) %>%
mutate(line_width = ifelse(grepl("TOTAL",
province),
2, 0.1)) %>%
mutate(line_cols = ifelse(grepl("TOTAL",
province),
"#AC8300",
"grey50")) %>%
arrange(year) %>%
arrange(province) %>%
#---Calculate annual totals for each province
group_by(year, province) %>%
mutate(value=mean(value, na.rm=TRUE)) %>%
slice(1) %>% ungroup() %>%
mutate(variable="var_spring_fall") %>%
filter(!grepl("TOTAL", province))
#-------------STATS
stat.var.before.after <- capa.melt.g %>%
group_by(before_after) %>%
summarise(mean = round(mean(value,
na.rm=TRUE), 2),
std_deviation = round(sd(value,
na.rm=TRUE), 2),
std_error = round((std_deviation / sqrt(n())), 2)) %>%
arrange(mean)
stat.var.before.after## # A tibble: 2 × 4
## before_after mean std_deviation std_error
## <fct> <dbl> <dbl> <dbl>
## 1 Before 67.5 25.8 4.07
## 2 After 72.6 26.8 3.79readr::write_tsv(stat.var.before.after, "tables/stats_oxytet_before_after.tsv")
#-------------PLOT
norm.test <- capa.melt.g %>%
rstatix::shapiro_test(var=value)
pwc1 <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0,
0, 1)) %>% ungroup() %>%
filter(var!=0) %>%
group_by(province) %>%
rstatix::wilcox_test(value ~before_after,
p.adjust.method = "BH",
detailed=TRUE) %>%
rstatix::add_xy_position() %>%
mutate(p_str = format(round(p, 4),
scientific = FALSE))
pwc1.ft <- flextable::flextable(pwc1)
pwc1.ft <- flextable::theme_vanilla(pwc1.ft)
pwc1.ftprovince | estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative | y.position | groups | xmin | xmax | p_str |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
BC | -5.500000 | value | Before | After | 4 | 5 | 5 | 0.2860 | -57.00000 | 28.0000000 | Wilcoxon | two.sided | 91.0 | [[character]] | 1 | 2 | 0.2860 |
AB | -9.200000 | value | Before | After | 4 | 5 | 5 | 0.2860 | -34.40000 | 13.5000000 | Wilcoxon | two.sided | 95.0 | [[character]] | 1 | 2 | 0.2860 |
SK | -10.250000 | value | Before | After | 4 | 5 | 6 | 0.4130 | -32.50000 | 6.5000000 | Wilcoxon | two.sided | 96.5 | [[character]] | 1 | 2 | 0.4130 |
MB | -6.519801 | value | Before | After | 4 | 5 | 8 | 0.7120 | -24.49992 | 39.3999779 | Wilcoxon | two.sided | 101.0 | [[character]] | 1 | 2 | 0.7120 |
ON | -5.250000 | value | Before | After | 4 | 5 | 4 | 0.1900 | -14.35000 | 4.0000000 | Wilcoxon | two.sided | 98.0 | [[character]] | 1 | 2 | 0.1900 |
QC | -6.250000 | value | Before | After | 4 | 5 | 4 | 0.1900 | -16.50000 | 11.5000000 | Wilcoxon | two.sided | 89.0 | [[character]] | 1 | 2 | 0.1900 |
NB | -1.250000 | value | Before | After | 4 | 5 | 9 | 0.9020 | -11.50000 | 8.4999578 | Wilcoxon | two.sided | 81.5 | [[character]] | 1 | 2 | 0.9020 |
PE | 2.075000 | value | Before | After | 4 | 5 | 12 | 0.7300 | -14.85000 | 22.0000000 | Wilcoxon | two.sided | 93.0 | [[character]] | 1 | 2 | 0.7300 |
NS | -8.499946 | value | Before | After | 4 | 5 | 1 | 0.0365 | -20.49997 | -0.5999206 | Wilcoxon | two.sided | 95.0 | [[character]] | 1 | 2 | 0.0365 |
#Remove cases where there is no variance (i.e. all value equal zero or represent 'NA' observations)
summary.stats <- capa.melt.g %>%
group_by(province) %>%
mutate(var=ifelse(var(value,
na.rm=TRUE)==0,
0, 1)) %>%
ungroup() %>%
filter(var!=0) %>%
group_by(province,
before_after) %>%
rstatix::get_summary_stats(value,
type = c("full"))
ft <- flextable::flextable(summary.stats)
ft <- flextable::theme_vanilla(ft)
ftprovince | before_after | variable | n | min | max | median | q1 | q3 | iqr | mad | mean | sd | se | ci |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
BC | Before | value | 4 | 24.00 | 77.0 | 74.25 | 60.375 | 76.250 | 15.875 | 3.336 | 62.375 | 25.656 | 12.828 | 40.824 |
BC | After | value | 5 | 48.00 | 88.0 | 78.50 | 73.000 | 81.000 | 8.000 | 8.154 | 73.700 | 15.344 | 6.862 | 19.053 |
AB | Before | value | 4 | 56.10 | 83.5 | 74.25 | 65.775 | 80.500 | 14.725 | 10.749 | 72.025 | 12.252 | 6.126 | 19.495 |
AB | After | value | 5 | 66.00 | 92.0 | 87.00 | 71.500 | 90.500 | 19.000 | 7.413 | 81.400 | 11.850 | 5.300 | 14.714 |
SK | Before | value | 4 | 57.00 | 93.0 | 78.75 | 66.750 | 88.875 | 22.125 | 17.050 | 76.875 | 16.484 | 8.242 | 26.230 |
SK | After | value | 5 | 84.50 | 93.5 | 88.00 | 86.500 | 89.500 | 3.000 | 2.224 | 88.400 | 3.399 | 1.520 | 4.220 |
MB | Before | value | 4 | 73.50 | 92.0 | 85.20 | 77.625 | 91.550 | 13.925 | 9.637 | 83.975 | 9.202 | 4.601 | 14.642 |
MB | After | value | 5 | 52.00 | 98.0 | 90.00 | 88.000 | 98.000 | 10.000 | 11.861 | 85.200 | 19.110 | 8.546 | 23.728 |
ON | Before | value | 4 | 80.15 | 89.5 | 86.50 | 83.037 | 89.125 | 6.088 | 4.077 | 85.662 | 4.435 | 2.218 | 7.058 |
ON | After | value | 5 | 85.00 | 95.0 | 91.50 | 87.500 | 94.500 | 7.000 | 5.189 | 90.700 | 4.367 | 1.953 | 5.423 |
QC | Before | value | 4 | 69.00 | 74.0 | 71.25 | 69.375 | 73.250 | 3.875 | 2.965 | 71.375 | 2.496 | 1.248 | 3.971 |
QC | After | value | 5 | 61.50 | 86.0 | 76.50 | 75.500 | 79.500 | 4.000 | 4.448 | 75.800 | 8.983 | 4.017 | 11.154 |
NB | Before | value | 4 | 66.00 | 73.0 | 69.25 | 66.750 | 71.875 | 5.125 | 4.077 | 69.375 | 3.400 | 1.700 | 5.411 |
NB | After | value | 5 | 63.00 | 78.5 | 69.00 | 69.000 | 72.000 | 3.000 | 4.448 | 70.300 | 5.630 | 2.518 | 6.991 |
PE | Before | value | 4 | 71.15 | 90.0 | 78.50 | 75.537 | 82.500 | 6.963 | 6.561 | 79.537 | 7.884 | 3.942 | 12.545 |
PE | After | value | 5 | 67.50 | 89.5 | 79.00 | 68.000 | 86.000 | 18.000 | 15.567 | 78.000 | 10.093 | 4.514 | 12.533 |
NS | Before | value | 4 | 69.50 | 79.4 | 74.00 | 71.000 | 77.225 | 6.225 | 5.189 | 74.225 | 4.535 | 2.268 | 7.217 |
NS | After | value | 5 | 79.00 | 92.0 | 80.00 | 80.000 | 84.000 | 4.000 | 1.483 | 83.000 | 5.385 | 2.408 | 6.687 |
fig2.i <- ggplot(capa.melt.g,
aes(x=before_after,
y=value)) +
facet_grid(~province) +
scale_y_continuous(limits=c(-2,110),
breaks=c(0,25,50,75,100)) +
geom_boxplot(aes(),
fill="#AC8300",
alpha=0.5,
outliers=FALSE) +
geom_point(position=position_jitter(width=0.15),
size=1,
alpha=0.5) +
ggpubr::stat_pvalue_manual(pwc1,
label="p_str",
hide.ns=FALSE,
size=2.5,
bracket.size=0.1,
vjust=-0.5,
step.increase=0.0000000000001,
bracket.nudge.y=-7,
tip.length = 0.015) +
theme(axis.text.x = element_text(size = 10,
angle=90,
vjust = 0.5),
axis.title = element_text(size=10),
axis.text.y = element_text(size=8),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
strip.background =element_rect(fill=scales::alpha("#AC8300",0.5)),
strip.text = element_text(colour = 'black',
size=8,
face="bold"))
fig2.i
Figure 2 (Full panel)
#::::::::::::::::::::
# Figure 2
#::::::::::::::::::::
fig2.top <-cowplot::plot_grid(fig2.a,
NULL,
fig2.b,
NULL,
fig2.c,
NULL,
ncol=6, labels=c('A',
'',
'B',
'',
'C',
''),
rel_widths=c(0.3,
0.03,
0.3,
0.03,
0.3,
0.03))
fig2.bottom.left <-cowplot::plot_grid(
fig2.d +
theme(axis.text.x = element_text(angle=90,
hjust=1),
axis.text = element_text(size=8),
axis.title = element_text(size=10)),
fig2.f +
theme(axis.text.x = element_text(angle=90,
hjust=1),
axis.text = element_text(size=8),
axis.title = element_text(size=10)),
fig2.h +
theme(axis.text.x = element_text(angle=90,
hjust=1),
axis.text = element_text(size=8),
axis.title = element_text(size=10)),
ncol=1,
labels=c('D',
'F',
'H'))
fig2.bottom.right <-cowplot::plot_grid(
fig2.e +
theme(
axis.text.x = element_text(angle=90,
hjust=1),
axis.text = element_text(size=8),
axis.title = element_text(size=10)),
fig2.g +
theme(axis.text.x = element_text(angle=90,
hjust=1),
axis.text = element_text(size=8),
axis.title = element_text(size=10)),
fig2.i +
theme(axis.text.x = element_text(angle=90,
hjust=1),
axis.text = element_text(size=8),
axis.title = element_text(size=10)),
ncol=1,
labels=c('E',
'G',
'I'))
fig2.bottom <- cowplot::plot_grid(
fig2.bottom.left,
fig2.bottom.right,
ncol=2,
rel_widths=c(0.25,
0.75))
fig2.full <-cowplot::plot_grid(fig2.top,
fig2.bottom,
ncol=1,
rel_heights=c(0.3,
0.7))
print(fig2.full)
Figure 3
Setup
capa.df.filt <- capa.df %>%
mutate_at(vars(respondent_colonies_percent,
bac_spring_oxytet,
bac_fall_oxytet,
nos_spring_fum,
nos_fall_fum,
var_spring_other,
var_fall_other),
as.numeric) %>%
mutate(respondent_colonies_percent = ifelse(
respondent_colonies_percent > 100,
100, respondent_colonies_percent)) %>%
mutate(oxytet_spring_fall = (bac_spring_oxytet +
bac_fall_oxytet)/2) %>%
mutate(fum_spring_fall = (nos_spring_fum +
nos_fall_fum)/2) %>%
mutate(var_spring_fall = (var_spring_other +
var_fall_other)/2) %>%
mutate(winterloss= winterloss_percent)
capa.melt <- melt(capa.df.filt ,
id.var=c("year",
"province",
"winterloss",
"raw_total",
"respondent_colonies")) %>%
mutate_at(vars("raw_total",
"respondent_colonies"),
as.numeric)
capa.melt.correlations <- melt(capa.df.filt,
id.var=c("year",
"province",
"winterloss_percent",
"raw_total",
"respondent_colonies")) %>%
filter(!is.na(winterloss_percent)) %>%
mutate_at(vars(winterloss_percent,
"raw_total",
"respondent_colonies"),
as.numeric) %>%
dplyr::rename("winterloss" = winterloss_percent)Figure 3E
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("bac_spring_oxytet|bac_fall_oxytet",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
arrange(year) %>% arrange(province) %>%
#---Calculate national totals
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="oxytet_spring_fall") %>%
filter(!grepl("TOTAL|NL", province)) %>%
group_by(year, variable) %>%
mutate(value=mean(value, na.rm=FALSE)) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>%
#---Add winterloss data
mutate(winterloss=(capa.melt.correlations %>% filter(grepl("^winterloss", variable) & grepl("TOTAL", province)))$winterloss)
#----------------------------------------------------
#Set functions for scaling secondary axis
#----------------------------------------------------
scale_fun <- function(x, scale, shift){
return ((x)*scale - shift)
}
scale_fun_inv <- function(x, scale, shift){
return ((x + shift)/scale)
}
#---------------------------------------------------- #Specify min/max of left (y1) and right (y2) y-axes
#Determined below values via inspection of min()/max() variables of interest
y1_min.left <- 10
y1_max.left <- 60
y2_min <- 0
y2_max <- 50
#----------------------------------------------------
scale_var.left = (y2_max - y2_min)/(y1_max.left - y1_min.left)
shift_var.left = y1_min.left - y2_min
#----------------------------------------------------
#Rename column for plotting
capa.melt.g.left <- capa.melt.g %>% dplyr::rename(usage = "value")
#Plot
fig3.a <- ggplot(capa.melt.g.left,
aes(x=year,
y=usage)) +
geom_smooth(method='loess',
span=1.5,
colour="#2980B9",
fill="#ADD3ED",
alpha=0.5,
size=0.7) +
geom_smooth(aes(y = scale_fun_inv(winterloss,
scale_var.left,
shift_var.left)),
method='loess',
span=1.5,
colour="black",
fill="grey50",
alpha=0.2,
size=0.7) +
geom_point(colour="#2980B9") +
geom_point(aes(y = scale_fun_inv(winterloss,
scale_var.left,
shift_var.left)),
colour="black") +
scale_x_continuous(breaks = seq(2015, 2023, 1)) +
scale_y_continuous(limits = c(y1_min.left,
y1_max.left),
sec.axis = sec_axis(~scale_fun(.,
scale_var.left,
shift_var.left),
name="Mortality (%)")) +
labs(x = "Year",
y = "Oxytetracycline use (%)",
color = "") +
theme(axis.title.y.left = element_text(size = 10,
colour = "#2980B9"),
axis.title.y.right = element_text(size = 10,
colour = "black"),
axis.title.x = element_text(size=10),
axis.text.x = element_text(size = 8,
angle=45,
hjust=1),
axis.text.y = element_text(size=8),
panel.background = element_blank(),
panel.grid.major = element_blank())
fig3.a
Figure 3B
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("nos_spring_fum|nos_fall_fum",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
arrange(year) %>% arrange(province) %>%
#---Calculate national totals
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="fum_spring_fall") %>%
filter(!grepl("TOTAL|NL", province)) %>%
group_by(year, variable) %>%
mutate(value=mean(value, na.rm=FALSE)) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>%
#---Add winterloss data
mutate(winterloss=(capa.melt.correlations %>% filter(grepl("^winterloss", variable) & grepl("TOTAL", province)))$winterloss)
#----------------------------------------------------
#Set functions for scaling secondary axis
#----------------------------------------------------
scale_fun <- function(x, scale, shift){
return ((x)*scale - shift)
}
scale_fun_inv <- function(x, scale, shift){
return ((x + shift)/scale)
}
#---------------------------------------------------- #Specify min/max of left (y1) and right (y2) y-axes
#Determined below values via inspection of min()/max() variables of interest
y1_min.mid <- 5
y1_max.mid <- 55
y2_min <- 0
y2_max <- 50
#----------------------------------------------------
scale_var.mid = (y2_max - y2_min)/(y1_max.mid - y1_min.mid)
shift_var.mid = y1_min.mid - y2_min
#----------------------------------------------------
#Rename column for plotting
capa.melt.g.mid <- capa.melt.g %>% dplyr::rename(usage = "value")
#Plot
fig3.b <- ggplot(capa.melt.g.mid,
aes(x=year,
y=usage)) +
geom_smooth(method='loess',
span=1.5,
colour="green4",
fill="#B2E8B3",
alpha=0.5,
size=0.7) +
geom_smooth(aes(y = scale_fun_inv(winterloss,
scale_var.mid,
shift_var.mid)),
method='loess',
span=1.5,
colour="black",
fill="grey50",
alpha=0.2,
size=0.7) +
geom_point(colour="green4") +
geom_point(aes(y = scale_fun_inv(winterloss,
scale_var.mid,
shift_var.mid)),
colour="black") +
scale_x_continuous(breaks = seq(2015, 2023, 1)) +
scale_y_continuous(limits = c(y1_min.mid,
y1_max.mid),
sec.axis = sec_axis(~scale_fun(.,
scale_var.mid,
shift_var.mid),
name="Mortality (%)")) +
labs(x = "Year",
y = "Fumagillin use (%)",
color = "green4") +
theme(axis.title.y.left = element_text(size = 10,
colour = "green4"),
axis.title.y.right = element_text(size = 10,
colour = "black"),
axis.title.x = element_text(size=10),
axis.text.x = element_text(size = 8,
angle=45,
hjust=1),
axis.text.y = element_text(size=8),
panel.background = element_blank(),
panel.grid.major = element_blank())
fig3.b
Figure 3C
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("var_spring_other|var_fall_other",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
arrange(year) %>% arrange(province) %>%
#---Calculate national totals
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="var_spring_fall") %>%
filter(!grepl("TOTAL|NL", province)) %>%
group_by(year, variable) %>%
mutate(value=mean(value, na.rm=FALSE)) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>%
#---Add winterloss data
mutate(winterloss=(capa.melt.correlations %>% filter(grepl("^winterloss", variable) & grepl("TOTAL", province)))$winterloss)
#----------------------------------------------------
#Set functions for scaling secondary axis
#----------------------------------------------------
scale_fun <- function(x, scale, shift){
return ((x)*scale - shift)
}
scale_fun_inv <- function(x, scale, shift){
return ((x + shift)/scale)
}
#---------------------------------------------------- #Specify min/max of left (y1) and right (y2) y-axes
#Determined below values via inspection of min()/max() variables of interest
y1_min.right <- 55 #Miticide min
y1_max.right <- 105 #Miticide max
y2_min <- 0 #Winterloss min
y2_max <- 50 #Winterloss max
#----------------------------------------------------
scale_var.right = (y2_max - y2_min)/(y1_max.right - y1_min.right)
shift_var.right = y1_min.right - y2_min
#----------------------------------------------------
#Rename column for plotting
capa.melt.g.right <- capa.melt.g %>% dplyr::rename(usage = "value")
#Plot
fig3.c <- ggplot(capa.melt.g.right,
aes(x=year,
y=usage)) +
geom_smooth(method='loess',
span=1.5,
colour="#AC8300",
fill="#FFE89F",
alpha=0.5,
size=0.7) +
geom_smooth(aes(y = scale_fun_inv(winterloss,
scale_var.right,
shift_var.right)),
method='loess',
span=1.5,
colour="black",
fill="grey50",
alpha=0.2,
size=0.7) +
geom_point(colour="#AC8300") +
geom_point(aes(y = scale_fun_inv(winterloss,
scale_var.right,
shift_var.right)),
colour="black") +
scale_x_continuous(breaks = seq(2015, 2023, 1)) +
scale_y_continuous(limits = c(y1_min.right,
y1_max.right),
sec.axis = sec_axis(~scale_fun(.,
scale_var.right,
shift_var.right),
name="Mortality (%)")) +
labs(x = "Year",
y = "Miticide use (%)",
color = "") +
theme(axis.title.y.left = element_text(size = 10,
colour = "#AC8300"),
axis.title.y.right = element_text(size = 10,
colour = "black"),
axis.title.x = element_text(size=10),
axis.text.x = element_text(size = 8,
angle=45,
hjust=1),
axis.text.y = element_text(size=8),
panel.background = element_blank(),
panel.grid.major = element_blank())
fig3.c
Figure 3D
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("bac_spring_oxytet|bac_fall_oxytet",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
arrange(year) %>% arrange(province) %>%
#---Calculate national totals
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="oxytet_spring_fall") %>%
filter(!grepl("TOTAL|NL", province)) %>%
group_by(year, variable) %>%
mutate(value=mean(value, na.rm=FALSE)) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>%
#---Add winterloss data
mutate(winterloss=(capa.melt.correlations %>% filter(grepl("^winterloss", variable) & grepl("TOTAL", province)))$winterloss)
#Plot
fig3.d <- ggplot(capa.melt.g,
aes(x=winterloss,
y=value)) +
smplot2::sm_statCorr(color = '#2980B9',
corr_method = 'pearson',
linetype="dashed",
text_size=3,
label_x = 20,
label_y = 15,
linewidth=0.7) +
xlab("Mortality (%)") +
ylab("Oxytetracycline use (%)") +
theme(axis.title.y.left = element_text(colour="#2980B9")) +
ylim(c(0,60)) +
geom_point() +
directlabels::geom_dl(aes(label=year),
method=list("bumpup",
cex=0.5,
rot =0,
hjust = -.5),
position=ggstance::position_dodgev(height=7)) +
theme(axis.title.y.left = element_text(size = 10,
colour = "#2980B9"),
axis.title.y.right = element_text(size = 10,
colour = "black"),
axis.title.x = element_text(size=10),
axis.text.x = element_text(size = 8),
axis.text.y = element_text(size=8),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank())
fig3.d
Figure 3E
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("nos_spring_fum|nos_fall_fum",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
arrange(year) %>% arrange(province) %>%
#---Calculate national totals
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="fum_spring_fall") %>%
filter(!grepl("TOTAL|NL", province)) %>%
group_by(year, variable) %>%
mutate(value=mean(value, na.rm=FALSE)) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>%
#---Add winterloss data
mutate(winterloss=(capa.melt.correlations %>% filter(grepl("^winterloss", variable) & grepl("TOTAL", province)))$winterloss)
#Plot
fig3.e <- ggplot(capa.melt.g,
aes(x=winterloss,
y=value)) +
smplot2::sm_statCorr(color = '#0f993d',
corr_method = 'pearson',
linetype="dashed",
text_size=3,
label_x = 20,
label_y = 15,
linewidth=0.7) +
xlab("Mortality (%)") +
ylab("Fumagillin use (%)") +
theme(axis.title.y.left = element_text(colour="#0f993d")) +
ylim(c(5,55)) +
geom_point() +
directlabels::geom_dl(aes(label=year),
method=list("bumpup",
cex=0.5,
rot =0,
hjust = -.5),
position=ggstance::position_dodgev(height=7)) +
theme(axis.title.y.left = element_text(size = 10,
colour = "#0f993d"),
axis.title.y.right = element_text(size = 10,
colour = "black"),
axis.title.x = element_text(size=10),
axis.text.x = element_text(size = 8),
axis.text.y = element_text(size=8),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank())
print(fig3.e)
Figure 3F
capa.melt.g <- capa.melt.correlations %>%
filter(grepl("var_spring_other|var_fall_other",
variable)) %>%
mutate_at(vars(value),
as.numeric) %>%
droplevels() %>%
mutate(province= factor(province,
levels=c("TOTAL",
"BC",
"AB",
"SK",
"MB",
"ON",
"QC",
"NB",
"PE",
"NS",
"NL"))) %>%
arrange(year) %>% arrange(province) %>%
#---Calculate national totals
group_by(year, province) %>% mutate(value=mean(value, na.rm=TRUE)) %>% slice(1) %>% ungroup() %>% mutate(variable="var_spring_fall") %>%
filter(!grepl("TOTAL|NL", province)) %>%
group_by(year, variable) %>%
mutate(value=mean(value, na.rm=FALSE)) %>%
slice(1) %>% mutate(province="TOTAL") %>% ungroup() %>%
#---Add winterloss data
mutate(winterloss=(capa.melt.correlations %>% filter(grepl("^winterloss", variable) & grepl("TOTAL", province)))$winterloss)
#Plot
fig3.f <- ggplot(capa.melt.g,
aes(x=winterloss,
y=value)) +
smplot2::sm_statCorr(color = '#AC8300',
corr_method = 'pearson',
linetype="dashed",
text_size=3,
label_x = 20,
label_y = 65,
linewidth=0.7) +
xlab("Mortality (%)") +
ylab("Miticide use (%)") +
theme(axis.title.y.left = element_text(colour="#AC8300")) +
ylim(c(60,90)) +
geom_point() +
directlabels::geom_dl(aes(label=year),
method=list("bumpup",
cex=0.5,
rot =0,
hjust = -.5)) +
theme(axis.title.y.left = element_text(size = 10,
colour = "#AC8300"),
axis.title.y.right = element_text(size = 10,
colour = "black"),
axis.title.x = element_text(size=10),
axis.text.x = element_text(size = 8),
axis.text.y = element_text(size=8),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank())
fig3.f
Table 1
capa.df.filt <- capa.df %>%
mutate_at(vars(respondent_colonies_percent,
bac_spring_oxytet,
bac_fall_oxytet,
nos_spring_fum,
nos_fall_fum,
var_spring_other,
var_fall_other),
as.numeric) %>%
mutate(respondent_colonies_percent = ifelse(
respondent_colonies_percent > 100,
100, respondent_colonies_percent)) %>%
mutate(oxytet_spring_fall = (bac_spring_combined +
bac_fall_combined)/2) %>%
mutate(fum_spring_fall = (nos_spring_combined +
nos_fall_combined)/2) %>%
mutate(var_spring_fall = (var_spring_other +
var_fall_other)/2) %>%
mutate(winterloss= winterloss_percent)
capa.melt <- melt(capa.df.filt ,
id.var=c("year",
"province",
"winterloss",
"raw_total",
"respondent_colonies")) %>%
mutate_at(vars("raw_total",
"respondent_colonies"),
as.numeric)
capa.melt.correlations <- melt(capa.df.filt,
id.var=c("year",
"province",
"winterloss_percent",
"raw_total",
"respondent_colonies")) %>%
filter(!is.na(winterloss_percent)) %>%
mutate_at(vars(winterloss_percent,
"raw_total",
"respondent_colonies"),
as.numeric) %>%
dplyr::rename("winterloss" = winterloss_percent)
#----------Air pollution and climate factors adjusted to 'bee year' format
air.vars <- cbind(readr::read_tsv("tables/air_quality_data_adj.tsv") %>%
filter(grepl("BC|AB|SK|MB|ON|QC|NB|PE|NS|NL",
Province)) %>%
mutate(year_prov = paste(Year_beeyear_lag,
"_",
Province,
sep="")) %>%
dcast(year_prov ~ Pollutant,
value.var = "mean_value") %>%
rowwise() %>%
mutate(NOX_custom = NO+NO2) %>%
ungroup(),
readr::read_tsv("tables/air_quality_data_adj.tsv") %>%
filter(grepl("BC|AB|SK|MB|ON|QC|NB|PE|NS|NL",
Province)) %>%
mutate(year_prov = paste(Year_beeyear_lag,
"_",
Province,
sep="")) %>%
dcast(year_prov ~ Pollutant,
value.var = "mean_winter") %>%
rowwise() %>%
mutate(NOX_custom = NO+NO2) %>%
ungroup() %>%
`colnames<-`(c("year_prov",
paste("winter",
colnames(.)[
2:length(colnames(.))
],
sep="_"))) %>%
select(-year_prov),
readr::read_tsv("tables/air_quality_data_adj.tsv") %>%
filter(grepl("BC|AB|SK|MB|ON|QC|NB|PE|NS|NL",
Province)) %>%
mutate(year_prov = paste(Year_beeyear_lag,
"_",
Province, sep="")) %>%
dcast(year_prov ~ Pollutant,
value.var = "mean_spring") %>%
rowwise() %>%
mutate(NOX_custom = NO+NO2) %>%
ungroup() %>%
`colnames<-`(c("year_prov",
paste("spring",
colnames(.)[
2:length(colnames(.))
],
sep="_"))) %>%
select(-year_prov),
readr::read_tsv("tables/air_quality_data_adj.tsv") %>%
filter(grepl("BC|AB|SK|MB|ON|QC|NB|PE|NS|NL",
Province)) %>%
mutate(year_prov = paste(Year_beeyear_lag,
"_",
Province,
sep="")) %>%
dcast(year_prov ~ Pollutant,
value.var = "mean_summer") %>%
rowwise() %>%
mutate(NOX_custom = NO+NO2) %>%
ungroup() %>%
`colnames<-`(c("year_prov",
paste("summer",
colnames(.)[
2:length(colnames(.))
],
sep="_"))) %>%
select(-year_prov),
readr::read_tsv("tables/air_quality_data_adj.tsv") %>%
filter(grepl("BC|AB|SK|MB|ON|QC|NB|PE|NS|NL",
Province)) %>%
mutate(year_prov = paste(Year_beeyear_lag,
"_",
Province,
sep="")) %>%
dcast(year_prov ~ Pollutant,
value.var = "mean_fall") %>%
rowwise() %>%
mutate(NOX_custom = NO+NO2) %>%
ungroup() %>%
`colnames<-`(c("year_prov",
paste("fall",
colnames(.)[
2:length(colnames(.))
], sep="_"))) %>%
select(-year_prov))
climate.vars <- readr::read_tsv("tables/climate_data_adj.tsv") %>%
filter(date_year_beeyear_lag != "2014") %>%
filter(grepl("BC|AB|SK|MB|ON|QC|NB|PE|NS|NL",
Prov_or_Ter)) %>%
mutate(year_prov = paste(date_year_beeyear_lag,
"_",
Prov_or_Ter, sep=""))
#---------------------------------------
capa.df.filt.mer <- merge(capa.df.filt %>%
mutate(year_prov = paste(year,
"_",
province,
sep="")),
climate.vars,
by="year_prov",
all=TRUE) %>%
merge(., air.vars,
by="year_prov",
all=TRUE)
capa.melt.correlations.mer <- melt(capa.df.filt.mer,
id.var=c("year",
"province",
"winterloss",
"raw_total",
"respondent_colonies",
"mean_Tm_spring",
"mean_Tm_summer",
"mean_Tm_fall",
"mean_Tm_winter",
"mean_P_spring",
"mean_P_summer",
"mean_P_fall",
"mean_P_winter",
"mean_S_spring",
"mean_S_summer",
"mean_S_fall",
"mean_S_winter",
"mean_Tm",
"mean_Tx",
"mean_Tn",
"mean_S",
"mean_P",
"mean_S_G",
"mean_Pd",
"mean_HDD",
"mean_CDD",
"CO",
"NO",
"NO2",
"NOX",
"O3",
"PM10",
"PM2.5",
"SO2",
"NOX_custom",
"winter_CO",
"winter_NO",
"winter_NO2",
"winter_NOX",
"winter_O3",
"winter_PM10",
"winter_PM2.5",
"winter_SO2",
"winter_NOX_custom",
"spring_CO",
"spring_NO",
"spring_NO2",
"spring_NOX",
"spring_O3",
"spring_PM10",
"spring_PM2.5",
"spring_SO2",
"spring_NOX_custom",
"summer_CO",
"summer_NO",
"summer_NO2",
"summer_NOX",
"summer_O3",
"summer_PM10",
"summer_PM2.5",
"summer_SO2",
"summer_NOX_custom",
"fall_CO",
"fall_NO",
"fall_NO2",
"fall_NOX",
"fall_O3",
"fall_PM10",
"fall_PM2.5",
"fall_SO2",
"fall_NOX_custom")) %>%
filter(!is.na(winterloss)) %>%
mutate_at(vars("winterloss",
"raw_total",
"respondent_colonies",
"mean_Tm_spring",
"mean_Tm_summer",
"mean_Tm_fall",
"mean_Tm_winter",
"mean_P_spring",
"mean_P_summer",
"mean_P_fall",
"mean_P_winter",
"mean_S_spring",
"mean_S_summer",
"mean_S_fall",
"mean_S_winter",
"mean_Tm",
"mean_S",
"CO",
"NO",
"NO2",
"NOX",
"O3",
"PM10",
"PM2.5",
"SO2",
"NOX_custom",
"mean_HDD",
"mean_CDD",
"winter_CO",
"winter_NO",
"winter_NO2",
"winter_NOX",
"winter_O3",
"winter_PM10",
"winter_PM2.5",
"winter_SO2",
"winter_NOX_custom",
"spring_CO",
"spring_NO",
"spring_NO2",
"spring_NOX",
"spring_O3",
"spring_PM10",
"spring_PM2.5",
"spring_SO2",
"spring_NOX_custom",
"summer_CO",
"summer_NO",
"summer_NO2",
"summer_NOX",
"summer_O3",
"summer_PM10",
"summer_PM2.5",
"summer_SO2",
"summer_NOX_custom",
"fall_CO",
"fall_NO",
"fall_NO2",
"fall_NOX",
"fall_O3",
"fall_PM10",
"fall_PM2.5",
"fall_SO2",
"fall_NOX_custom"),
as.numeric)
#::::::::::::::::::::::::::::
# Table 1
#::::::::::::::::::::::::::::
capa.melt.g <- capa.melt.correlations.mer %>%
filter(grepl("oxytet_spring_fall",
variable)) %>%
mutate_at(vars(value, winterloss),
as.numeric) %>%
filter(grepl("BC|AB|SK|MB|ON|QC|NB|PE|NS|NL",
province)) %>% droplevels() %>%
arrange(year) %>%
arrange(province) %>%
ungroup()
capa.melt.g2 <- capa.melt.correlations.mer %>%
filter(grepl("fum_spring_fall",
variable)) %>%
mutate_at(vars(value, winterloss),
as.numeric) %>%
filter(grepl("BC|AB|SK|MB|ON|QC|NB|PE|NS|NL",
province)) %>%
droplevels() %>%
arrange(year) %>%
arrange(province) %>%
ungroup()
capa.melt.g3 <- capa.melt.correlations.mer %>%
filter(grepl("var_spring_fall",
variable)) %>%
mutate_at(vars(value,
winterloss),
as.numeric) %>%
filter(grepl("BC|AB|SK|MB|ON|QC|NB|PE|NS|NL",
province)) %>%
droplevels() %>%
arrange(year) %>%
arrange(province) %>%
ungroup()
capa.melt.g$value2 <- capa.melt.g2$value
capa.melt.g$value3 <- capa.melt.g3$value
capa.melt.g <- capa.melt.g %>%
rowwise() %>%
mutate(value4 = mean(c(value,
value2),
na.rm=TRUE)) %>%
ungroup() %>%
mutate(year_prov = paste(year,
province,
sep="_")) %>%
mutate_at(vars(province),
as.factor) %>%
filter(!is.na(winterloss)) %>%
filter(!grepl("QC|NL",
province)) %>%
mutate_at(vars(value4,
value3,
mean_Tm,
mean_P,
mean_S,
mean_Tm_winter,
mean_Tm_spring,
mean_Tm_summer,
mean_Tm_fall,
mean_P_winter,
mean_P_spring,
mean_P_summer,
mean_P_fall,
mean_S_winter,
mean_S_spring,
mean_S_summer,
mean_S_fall,
NO,
NO2,
SO2,
PM2.5,
CO,
O3), funs(scale(.,
center=TRUE,
scale=TRUE)[,1])) %>%
mutate_at(vars(value3),
funs(randomForest::na.roughfix(.))) %>%
dplyr::rename(Antibiotics = value4,
Miticides=value3,
Tm = mean_Tm,
Tm_wi = mean_Tm_winter,
Tm_sp = mean_Tm_spring,
Tm_su = mean_Tm_summer,
Tm_fa = mean_Tm_fall,
P=mean_P,
P_wi = mean_P_winter,
P_sp = mean_P_spring,
P_su = mean_P_summer,
P_fa = mean_P_fall,
S = mean_S,
S_wi = mean_S_winter,
S_sp = mean_S_spring,
S_su = mean_S_summer,
S_fa = mean_S_fall)
#Plot correlogram to see individual correlations between predictors
#=================================================================================================================
GGally::ggpairs(capa.melt.g %>%
filter(!is.na(Miticides)) %>%
filter(!is.na(CO)) %>%
filter(!is.na(O3)) %>%
dplyr::select(winterloss,
Antibiotics,
Miticides,
Tm,
P,
S,
NO,
NO2,
CO,
O3,
SO2,
PM2.5))
#=================================================================================================================
model_1 <- lm(winterloss ~ Antibiotics,
data=capa.melt.g,
na.action=na.omit)
model_2 <- lme4::lmer(winterloss ~ Antibiotics +
(1 | province),
data=capa.melt.g,
REML=FALSE, na.action=na.omit)
model_3 <- lme4::lmer(winterloss ~ Antibiotics +
Miticides +
Tm +
P +
S +
(1 | province),
data=capa.melt.g %>%
filter(!is.na(Miticides)) %>%
filter(!is.na(CO)) %>%
filter(!is.na(O3)),
REML=FALSE, na.action=na.omit)
model_4 <- lme4::lmer(winterloss ~ Antibiotics +
Miticides +
NO +
NO2 +
CO +
O3 +
SO2 +
PM2.5 +
(1 | province),
data=capa.melt.g %>%
filter(!is.na(Miticides)) %>%
filter(!is.na(CO)) %>%
filter(!is.na(O3)),
REML=FALSE, na.action=na.omit)
model_5 <- lme4::lmer(winterloss ~ Antibiotics +
Miticides +
Tm +
P +
S +
NO +
NO2 +
CO +
O3 +
SO2 +
PM2.5 +
(1 | province),
data=capa.melt.g %>%
filter(!is.na(Miticides)) %>%
filter(!is.na(CO)) %>%
filter(!is.na(O3)),
REML=FALSE, na.action=na.omit)
model_6 <- lme4::lmer(winterloss ~ Antibiotics +
Miticides +
Tm_wi +
Tm_sp +
Tm_su +
Tm_fa +
P +
S +
NO +
NO2 +
CO +
O3 +
SO2 +
PM2.5 +
(1 | province),
data=capa.melt.g %>%
filter(!is.na(Miticides)) %>%
filter(!is.na(CO)) %>%
filter(!is.na(O3)),
REML=FALSE, na.action=na.omit)
model_7 <- lme4::lmer(winterloss ~ Antibiotics +
Miticides +
Tm +
S +
P_wi +
P_sp +
P_su +
P_fa +
NO +
NO2 +
CO +
O3 +
SO2 +
PM2.5 +
(1 | province),
data=capa.melt.g %>%
filter(!is.na(Miticides)) %>%
filter(!is.na(CO)) %>%
filter(!is.na(O3)),
REML=FALSE, na.action=na.omit)
model_8 <- lme4::lmer(winterloss ~
Antibiotics +
Miticides +
Tm +
P +
S_wi +
S_sp +
S_su +
S_fa +
NO +
NO2 +
CO +
O3 +
SO2 +
PM2.5 +
(1 | province),
data=capa.melt.g %>%
filter(!is.na(Miticides)) %>%
filter(!is.na(CO)) %>%
filter(!is.na(O3)),
REML=FALSE,
na.action=na.omit)
model_9 <- lme4::lmer(winterloss ~
Antibiotics +
Miticides +
Tm_wi +
Tm_sp +
Tm_su +
Tm_fa +
P_wi +
P_sp +
P_su +
P_fa +
S_wi +
S_sp +
S_su +
S_fa +
NO +
NO2 +
CO +
O3 +
SO2 +
PM2.5 +
(1 | province),
data=capa.melt.g %>%
filter(!is.na(Miticides)) %>%
filter(!is.na(CO)) %>%
filter(!is.na(O3)),
REML=FALSE,
na.action=na.omit)
sjPlot::tab_model(model_1,
model_2,
model_3,
model_4,
model_5,
model_6,
model_7,
model_8,
model_9,
show.ci=FALSE,
digits.p = 4,
show.aic = TRUE,
show.aicc = TRUE,
show.loglik = TRUE,
show.stat = TRUE,
show.re.var = TRUE,
show.icc=TRUE,
CSS = list(sjPlot::css_theme("regression"),
sjPlot::css_theme("cells"),
css.table = '+font-size: 8;'),
dv.labels = c("Model 1",
"Model 2",
"Model 3",
"Model 4",
"Model 5",
"Model 6",
"Model 7",
"Model 8",
"Model 9"))| Model 1 | Model 2 | Model 3 | Model 4 | Model 5 | Model 6 | Model 7 | Model 8 | Model 9 | |||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | Statistic | p | Estimates | Statistic | p | Estimates | Statistic | p | Estimates | Statistic | p | Estimates | Statistic | p | Estimates | Statistic | p | Estimates | Statistic | p | Estimates | Statistic | p | Estimates | Statistic | p |
| (Intercept) | 25.96 | 20.40 | <0.0001 | 25.96 | 13.56 | <0.0001 | 25.46 | 22.52 | <0.0001 | 24.96 | 21.57 | <0.0001 | 25.24 | 22.28 | <0.0001 | 24.99 | 30.94 | <0.0001 | 25.50 | 22.47 | <0.0001 | 25.13 | 21.71 | <0.0001 | 25.58 | 31.88 | <0.0001 |
| Antibiotics | -4.12 | -3.21 | 0.0020 | -4.65 | -3.32 | 0.0015 | -4.71 | -3.51 | 0.0009 | -4.27 | -3.67 | 0.0006 | -4.72 | -3.42 | 0.0013 | -5.29 | -5.18 | <0.0001 | -4.54 | -3.17 | 0.0028 | -4.11 | -2.71 | 0.0095 | -5.83 | -5.32 | <0.0001 |
| Miticides | 0.47 | 0.43 | 0.6676 | -0.82 | -0.78 | 0.4404 | -0.68 | -0.63 | 0.5344 | 1.32 | 1.60 | 0.1157 | -0.54 | -0.50 | 0.6225 | -0.58 | -0.54 | 0.5924 | 1.30 | 1.64 | 0.1095 | ||||||
| Tm | 5.18 | 2.42 | 0.0191 | 1.24 | 0.57 | 0.5685 | 0.86 | 0.34 | 0.7364 | 1.92 | 0.83 | 0.4123 | |||||||||||||||
| P | -10.29 | -3.69 | 0.0005 | -2.62 | -0.84 | 0.4057 | -3.24 | -1.40 | 0.1694 | -2.13 | -0.67 | 0.5050 | |||||||||||||||
| S | 5.08 | 2.60 | 0.0122 | 2.23 | 0.99 | 0.3293 | 5.44 | 3.26 | 0.0021 | 2.65 | 1.13 | 0.2641 | |||||||||||||||
| NO | -10.19 | -2.53 | 0.0146 | -8.74 | -1.87 | 0.0679 | -5.67 | -1.68 | 0.1002 | -9.99 | -2.04 | 0.0471 | -8.53 | -1.74 | 0.0893 | -8.53 | -2.41 | 0.0210 | |||||||||
| NO2 | 9.14 | 3.64 | 0.0006 | 8.33 | 2.47 | 0.0170 | 7.60 | 3.13 | 0.0031 | 8.53 | 2.43 | 0.0193 | 7.41 | 2.04 | 0.0474 | 8.85 | 3.44 | 0.0014 | |||||||||
| CO | 4.43 | 3.05 | 0.0036 | 4.10 | 2.65 | 0.0110 | 2.38 | 2.09 | 0.0424 | 4.40 | 2.77 | 0.0081 | 3.89 | 2.48 | 0.0170 | 2.63 | 2.32 | 0.0260 | |||||||||
| O3 | -0.44 | -0.17 | 0.8659 | -0.59 | -0.22 | 0.8251 | 1.13 | 0.60 | 0.5506 | -1.46 | -0.54 | 0.5896 | -0.44 | -0.16 | 0.8697 | -0.79 | -0.41 | 0.6810 | |||||||||
| SO2 | -2.39 | -1.79 | 0.0798 | -2.46 | -1.63 | 0.1101 | -2.26 | -2.09 | 0.0427 | -2.86 | -1.89 | 0.0650 | -2.31 | -1.52 | 0.1355 | -2.75 | -2.60 | 0.0132 | |||||||||
| PM2 5 | 2.14 | 1.50 | 0.1411 | 1.75 | 1.16 | 0.2500 | 0.32 | 0.29 | 0.7741 | 1.63 | 1.07 | 0.2884 | 1.88 | 1.24 | 0.2227 | 0.52 | 0.49 | 0.6280 | |||||||||
| Tm wi | -10.81 | -5.74 | <0.0001 | -12.14 | -6.18 | <0.0001 | |||||||||||||||||||||
| Tm sp | 8.15 | 6.33 | <0.0001 | 8.68 | 5.96 | <0.0001 | |||||||||||||||||||||
| Tm su | -4.30 | -4.49 | 0.0001 | -4.12 | -4.46 | 0.0001 | |||||||||||||||||||||
| Tm fa | 5.93 | 2.97 | 0.0049 | 7.88 | 3.20 | 0.0028 | |||||||||||||||||||||
| P wi | -3.09 | -0.97 | 0.3351 | -3.04 | -1.33 | 0.1898 | |||||||||||||||||||||
| P sp | 0.87 | 0.39 | 0.7010 | 1.06 | 0.66 | 0.5143 | |||||||||||||||||||||
| P su | -1.93 | -1.56 | 0.1249 | -1.38 | -1.46 | 0.1515 | |||||||||||||||||||||
| P fa | 0.46 | 0.22 | 0.8287 | -0.59 | -0.40 | 0.6907 | |||||||||||||||||||||
| S wi | 1.25 | 0.53 | 0.5995 | 1.78 | 0.97 | 0.3395 | |||||||||||||||||||||
| S sp | -0.27 | -0.12 | 0.9051 | 4.40 | 2.62 | 0.0125 | |||||||||||||||||||||
| S su | 0.05 | 0.04 | 0.9666 | 0.58 | 0.73 | 0.4674 | |||||||||||||||||||||
| S fa | 1.83 | 1.19 | 0.2407 | 1.76 | 1.35 | 0.1862 | |||||||||||||||||||||
| Random Effects | |||||||||||||||||||||||||||
| σ2 | 94.80 | 67.81 | 51.80 | 51.79 | 25.78 | 49.66 | 51.01 | 22.56 | |||||||||||||||||||
| τ00 | 18.80 province | 0.00 province | 0.88 province | 0.00 province | 0.00 province | 0.00 province | 0.00 province | 0.00 province | |||||||||||||||||||
| ICC | 0.17 | 0.02 | |||||||||||||||||||||||||
| N | 8 province | 7 province | 7 province | 7 province | 7 province | 7 province | 7 province | 7 province | |||||||||||||||||||
| Observations | 72 | 72 | 61 | 61 | 61 | 61 | 61 | 61 | 61 | ||||||||||||||||||
| R2 / R2 adjusted | 0.128 / 0.116 | 0.160 / 0.299 | 0.256 / NA | 0.425 / 0.435 | 0.433 / NA | 0.719 / NA | 0.457 / NA | 0.442 / NA | 0.754 / NA | ||||||||||||||||||
| AIC | 550.901 | 542.568 | 431.651 | 414.696 | 410.304 | 374.527 | 404.492 | 406.993 | 365.126 | ||||||||||||||||||
| AICc | 551.254 | 543.165 | 434.420 | 420.084 | 419.434 | 388.759 | 418.725 | 421.225 | 394.964 | ||||||||||||||||||
| log-Likelihood | -272.451 | -270.123 | -215.164 | -207.438 | -206.945 | -185.673 | -205.664 | -206.481 | -181.604 | ||||||||||||||||||
Supplementary Data 2
In additional to the model evaluation above, we used the dredge() function from the MuMIn package to systematically evaluate all possible subsets of predictor variables from the global model.
#Running with default functions takes ~17 hours to finish.
#dd <- dredge(model_9)
#Filter select top models
#dd.top <- as_tibble(as.data.frame(subset(dd, as.numeric(delta) < 4)))
#Save top models
#saveRDS(dd.top, "RDS/dredge_top_models.RDS")
#Alternatively, load the pre-computed results table
dd.top <- readRDS("RDS/dredge_top_models.RDS")
reactable::reactable(dd.top, width=1000, striped = TRUE, highlight = TRUE)Figure 3G
#Set seed for reproducible results from partR2 package
set.seed(123)
#Setup cores for multithreading
parallel::detectCores()## [1] 22future::plan(multisession, workers = 12)
#Run partR2 command get statistics for Model 6
R2_mod <- partR2::partR2(model_6,
partvars = c("Antibiotics",
"Miticides",
"Tm_wi",
"Tm_sp",
"Tm_su",
"Tm_fa",
"P",
"S",
"NO2",
"NO",
"CO",
"O3",
"SO2",
"PM2.5"),
R2_type = "conditional",
nboot = 100,
max_level=1,
parallel = TRUE)
summary(R2_mod)##
##
## R2 (conditional) and 95% CI for the full model:
## R2 CI_lower CI_upper ndf
## 0.7192 0.7028 0.8587 15
##
## ----------
##
## Part (semi-partial) R2:
## Predictor(s) R2 CI_lower CI_upper ndf
## Model 0.7192 0.7028 0.8587 15
## Antibiotics 0.1256 0.0000 0.3279 14
## Miticides 0.0121 0.0000 0.2354 14
## Tm_wi 0.1544 0.0043 0.3505 14
## Tm_sp 0.0000 0.0000 0.0866 14
## Tm_su 0.0893 0.0000 0.2970 14
## Tm_fa 0.0412 0.0000 0.2587 14
## P 0.0091 0.0000 0.2331 14
## S 0.1245 0.0000 0.3270 14
## NO2 0.1753 0.0169 0.3665 14
## NO 0.0132 0.0000 0.2364 14
## CO 0.0205 0.0000 0.2421 14
## O3 0.0017 0.0000 0.2272 14
## SO2 0.0204 0.0000 0.2421 14
## PM2.5 0.0004 0.0000 0.2262 14
##
## ----------
##
## Inclusive R2 (SC^2 * R2):
## Predictor IR2 CI_lower CI_upper
## Antibiotics 0.0503 0.0105 0.1398
## Miticides 0.0156 0.0003 0.0605
## Tm_wi 0.0278 0.0005 0.0938
## Tm_sp 0.0583 0.0173 0.1385
## Tm_su 0.0033 0.0000 0.0296
## Tm_fa 0.0039 0.0000 0.0458
## P 0.0350 0.0020 0.1167
## S 0.0068 0.0000 0.0466
## NO 0.0064 0.0000 0.0434
## NO2 0.0989 0.0459 0.1886
## CO 0.0879 0.0295 0.1549
## O3 0.0008 0.0000 0.0192
## SO2 0.0142 0.0000 0.0566
## PM2.5 0.0736 0.0247 0.1575
##
## ----------
##
## Structure coefficients r(Yhat,x):
## Predictor SC CI_lower CI_upper
## Antibiotics -0.2644 -0.4359 -0.1179
## Miticides 0.1471 -0.0020 0.2692
## Tm_wi -0.1966 -0.3440 -0.0016
## Tm_sp 0.2848 0.1494 0.4337
## Tm_su 0.0677 -0.0971 0.1988
## Tm_fa -0.0741 -0.2475 0.0972
## P -0.2207 -0.3850 -0.0467
## S -0.0972 -0.2399 0.0418
## NO 0.0944 -0.0223 0.2370
## NO2 0.3708 0.2417 0.4833
## CO 0.3495 0.1941 0.4444
## O3 0.0325 -0.1290 0.1607
## SO2 0.1406 -0.0038 0.2699
## PM2.5 0.3200 0.1819 0.4588
##
## ----------
##
## Beta weights (standardised estimates)
## Predictor BW CI_lower CI_upper
## Antibiotics -0.5204 -0.6901 -0.3448
## Miticides 0.1387 -0.0066 0.2795
## Tm_wi -1.1191 -1.4054 -0.7297
## Tm_sp 0.8790 0.6416 1.0781
## Tm_su -0.4586 -0.6183 -0.2708
## Tm_fa 0.6206 0.2591 0.9709
## P -0.3481 -0.8329 0.1554
## S 0.4884 0.2762 0.7259
## NO -0.5349 -1.2192 0.0047
## NO2 0.6688 0.3611 1.1194
## CO 0.2483 0.0552 0.4507
## O3 0.1186 -0.2304 0.4748
## SO2 -0.2321 -0.4648 -0.0145
## PM2.5 0.0307 -0.1611 0.2541
##
## ----------
##
## Parametric bootstrapping resulted in warnings or messages:
## Check r2obj$boot_warnings and r2obj$boot_messages.print(R2_mod$R2 %>% arrange(desc(estimate)) %>%
filter(grepl("Antibiotics|",
term)),
n=100)## # A tibble: 15 × 5
## term estimate CI_lower CI_upper ndf
## <chr> <dbl> <dbl> <dbl> <int>
## 1 Full 0.719 0.703 0.859 15
## 2 NO2 0.175 0.0169 0.367 14
## 3 Tm_wi 0.154 0.00433 0.350 14
## 4 Antibiotics 0.126 0 0.328 14
## 5 S 0.124 0 0.327 14
## 6 Tm_su 0.0893 0 0.297 14
## 7 Tm_fa 0.0412 0 0.259 14
## 8 CO 0.0205 0 0.242 14
## 9 SO2 0.0204 0 0.242 14
## 10 NO 0.0132 0 0.236 14
## 11 Miticides 0.0121 0 0.235 14
## 12 P 0.00914 0 0.233 14
## 13 O3 0.00169 0 0.227 14
## 14 PM2.5 0.000390 0 0.226 14
## 15 Tm_sp 0 0 0.0866 14print(R2_mod$SC %>% arrange(desc(estimate)) %>%
filter(grepl("Antibiotics|",
term)),
n=100)## # A tibble: 14 × 4
## term estimate CI_lower CI_upper
## <chr> <dbl> <dbl> <dbl>
## 1 NO2 0.371 0.242 0.483
## 2 CO 0.350 0.194 0.444
## 3 PM2.5 0.320 0.182 0.459
## 4 Tm_sp 0.285 0.149 0.434
## 5 Miticides 0.147 -0.00201 0.269
## 6 SO2 0.141 -0.00383 0.270
## 7 NO 0.0944 -0.0223 0.237
## 8 Tm_su 0.0677 -0.0971 0.199
## 9 O3 0.0325 -0.129 0.161
## 10 Tm_fa -0.0741 -0.248 0.0972
## 11 S -0.0972 -0.240 0.0418
## 12 Tm_wi -0.197 -0.344 -0.00161
## 13 P -0.221 -0.385 -0.0467
## 14 Antibiotics -0.264 -0.436 -0.118print(R2_mod$BW %>% arrange(desc(estimate)) %>%
filter(grepl("Antibiotics|",
term)),
n=100)## # A tibble: 14 × 4
## term estimate CI_lower CI_upper
## <chr> <dbl> <dbl> <dbl>
## 1 Tm_sp 0.879 0.642 1.08
## 2 NO2 0.669 0.361 1.12
## 3 Tm_fa 0.621 0.259 0.971
## 4 S 0.488 0.276 0.726
## 5 CO 0.248 0.0552 0.451
## 6 Miticides 0.139 -0.00658 0.280
## 7 O3 0.119 -0.230 0.475
## 8 PM2.5 0.0307 -0.161 0.254
## 9 SO2 -0.232 -0.465 -0.0145
## 10 P -0.348 -0.833 0.155
## 11 Tm_su -0.459 -0.618 -0.271
## 12 Antibiotics -0.520 -0.690 -0.345
## 13 NO -0.535 -1.22 0.00473
## 14 Tm_wi -1.12 -1.41 -0.730partR2.panel <- cowplot::plot_grid(partR2::forestplot(R2_mod,
type = "R2",
text_size = 10),
partR2::forestplot(R2_mod,
type = "IR2",
text_size = 10),
partR2::forestplot(R2_mod,
type = "SC",
text_size = 10),
partR2::forestplot(R2_mod,
type = "BW",
text_size = 10)
)
print(partR2.panel)
fig3.g <- partR2::forestplot(R2_mod,
type = "SC",
text_size = 10)
print(fig3.g)
Figure 3H
#proportion of variance in the dependent variable that can be explained by the independent variable.
R2_mod.df <- R2_mod$R2 %>%
mutate(relative_estimate = estimate) %>%
filter(term!="Full") %>%
arrange((estimate)) %>%
mutate(term = case_when(
term == "Temperature_winter" ~ "Tm_wi",
term == "Temperature_spring" ~ "Tm_sp",
term == "Temperature_summer" ~ "Tm_su",
term == "Temperature_fall" ~ "Tm_fa",
TRUE ~ term)) %>%
mutate(term = factor(term,
levels=rev(term)))
fig3.h <- ggplot(R2_mod.df,
aes(x=relative_estimate,
y=term)) +
coord_flip() +
xlab(bquote('Semi-partial R'^'2')) +
geom_bar(stat="identity",
color=NA,
width=0.9) +
theme(axis.text.x = element_text(size = 8,
angle=45,
hjust=1),
panel.background = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.y = element_text(size=8))
print(fig3.h)
Figure 3 (Full panel)
fig3.top <- cowplot::plot_grid(fig3.a,
fig3.b,
fig3.c,
ncol=3,
labels = c('A',
'B',
'C'))## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'fig3.mid <- cowplot::plot_grid(fig3.d,
NULL,
fig3.e,
NULL,
fig3.f,
NULL,
ncol=6,
labels = c('D',
'',
'E',
'',
'F',
''),
rel_widths=c(0.3,
0.033,
0.3,
0.033,
0.3,
0.033))## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'fig3.bottom <- cowplot::plot_grid(
NULL,
fig3.g,
fig3.h +
ggtitle("Variance uniquely explained by each predictor") +
theme(plot.title = element_text(size=7,
vjust=-1,
hjust=0.5)),
NULL,
labels = c('',
'G',
'H',
''),
ncol=4,
rel_widths=c(0.15,
0.3,
0.4,
0.05))
cowplot::plot_grid(fig3.top,
fig3.mid,
fig3.bottom,
ncol=1,
rel_heights=c(0.33,
0.33,
0.33))