7 PCB - IAV
7.2 Data
meta <-
read.csv("Input Files/metadata_tidy.csv")
maxlod <-
read.csv("Output Files/cleaned_pcb_maxLOD.csv") %>%
dplyr::rename("Sample" = 1) %>%
select(!pcb180 &
!pcb183 &
!pcb187)
viro <-
read.csv("Input Files/Hg_virology_2023.csv", sep=',', strip.white=TRUE) %>%
select(Tag.ID., year, IAV, IAV.Nasal, IAV.Conj, IAV.Rectal) %>%
setNames(c("Sample", "year", "IAV", "Nasal", "Conj", "Rectal"))
viro_two <-
viro %>%
filter(Sample == "204" |
Sample == "268")
viro_tidy <-
viro %>%
filter(Sample %in% maxlod$Sample) %>%
rbind(., viro_two)7.4 IAV ~ Variables
## Analysis of Deviance Table (Type II tests)
##
## Response: iavbin
## LR Chisq Df Pr(>Chisq)
## year 1.1267 1 0.2885# Sex
tryCatch({
chisq.test(table(data$sex, data$iavbin))
}, warning = function(w) {
fisher.test(table(data$sex, data$iavbin))
})##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(data$sex, data$iavbin)
## X-squared = 0.10714, df = 1, p-value = 0.7434# Location
tryCatch({
chisq.test(table(data$location, data$iavbin))
}, warning = function(w) {
fisher.test(table(data$location, data$iavbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(data$location, data$iavbin)
## p-value = 0.04527
## alternative hypothesis: two.sided## # A tibble: 3 × 6
## group1 group2 n p p.adj p.adj.signif
## * <chr> <chr> <int> <dbl> <dbl> <chr>
## 1 Great Point Monomoy 87 0.186 0.372 ns
## 2 Great Point Muskeget 46 1 1 ns
## 3 Monomoy Muskeget 121 0.0365 0.109 ns# Molt Stage
tryCatch({
chisq.test(table(data$molt.stage, data$iavbin))
}, warning = function(w) {
fisher.test(table(data$molt.stage, data$iavbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(data$molt.stage, data$iavbin)
## p-value = 0.03975
## alternative hypothesis: two.sided## # A tibble: 6 × 6
## group1 group2 n p p.adj p.adj.signif
## * <chr> <chr> <int> <dbl> <dbl> <chr>
## 1 II III 16 1 1 ns
## 2 II IV 24 1 1 ns
## 3 II V 68 1 1 ns
## 4 III IV 38 0.326 1 ns
## 5 III V 82 0.506 1 ns
## 6 IV V 90 0.0159 0.0954 ns7.5 IAV ~ PCB p/a
tryCatch({
chisq.test(table(data$iav, data$pcbbin))
}, warning = function(w) {
fisher.test(table(data$iav, data$pcbbin))
})##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(data$iav, data$pcbbin)
## X-squared = 0.2492, df = 1, p-value = 0.61767.6 IAV ~ PCB conc
## Analysis of Deviance Table (Type II tests)
##
## Response: iavbin
## LR Chisq Df Pr(>Chisq)
## sumpcb 0.03147 1 0.85927.7 Figures
7.7.1 Presence/Absence
iavmosaic <-
data %>%
mutate(iav=ifelse(iav=="pos", "IAV+", "IAV-"),
iav=factor(iav, levels=c("IAV+", "IAV-")),
pcbbin=ifelse(pcbbin=="1", "Present", "Absent"),
pcbbin=factor(pcbbin, levels=c("Present", "Absent"))) %>%
ggplot(data=.) +
geom_mosaic(aes(x=product(pcbbin, iav), fill=pcbbin)) +
scale_fill_manual(values=c("Absent"="#b4b4b4", "Present"="#7eaaac"),
name="PCB Status") +
scale_y_continuous(labels = scales::percent) +
labs(y="Proportion of Samples") +
theme_classic() +
theme(legend.position=c(0.81, 0.85),
axis.title.x=element_blank(),
text=element_text(size=13),
axis.text = element_text(size=13))
iavmosaic
7.7.2 Concentration
iavbox <-
data %>%
mutate(iav=ifelse(iav=="pos", "IAV+", "IAV-"),
iav=factor(iav, levels=c("IAV+", "IAV-"))) %>%
ggplot(., aes(x=iav, y=log(sumpcb))) +
geom_boxplot(fill = "#7eaaac", color = "black") +
labs(y="log(ΣPCB) (ng/g wet weight)") +
stat_n_text(y.pos=3) +
theme_classic() +
theme(axis.title.x=element_blank(),
text=element_text(size=13),
axis.text = element_text(size=13))
iavbox## Warning: Removed 94 rows containing non-finite outside
## the scale range (`stat_boxplot()`).## Warning: Removed 94 rows containing non-finite outside
## the scale range (`stat_n_text()`).
7.8 Figures - Presentation
prop_pres <-
data %>%
group_by(iav, pcbbin) %>%
summarize(sum = n()) %>%
mutate(prop = ifelse(iav == "neg", sum/(62+24), sum/(32+9))) %>%
filter(pcbbin == "1") %>%
mutate(iav = ifelse(iav == "neg", "IAV-", "IAV+")) %>%
ggplot(., aes(x = iav, y = prop, fill = iav)) +
geom_col(position = "identity") +
scale_fill_manual(values = c("IAV+" = "#ACD7CA", "IAV-" = "#7FC1DB"),
name = "IAV Status") +
scale_y_continuous(labels = scales::percent) +
labs(x = NULL, y = "Proportion of Samples with PCBs") +
theme_bw() +
theme(panel.grid = element_blank(),
text = element_text(size = 20),
legend.position = "none")## `summarise()` has grouped output by 'iav'. You
## can override using the `.groups` argument.
conc_pres <-
data %>%
mutate(iav=ifelse(iav=="pos", "IAV+", "IAV-")) %>%
ggplot(., aes(x=iav, y=log(sumpcb), fill = iav)) +
geom_boxplot () +
geom_point(position = position_jitter(width = 0.2)) +
labs(y="log(ΣPCB Concentration)") +
scale_fill_manual(values = c("IAV+" = "#ACD7CA", "IAV-" = "#7FC1DB")) +
theme_classic() +
theme(axis.title.x=element_blank(),
text=element_text(size=20),
legend.position = "none")
conc_pres## Warning: Removed 94 rows containing non-finite outside
## the scale range (`stat_boxplot()`).
## Warning: Removed 94 rows containing non-finite outside
## the scale range (`stat_boxplot()`).
7.9 Figures - Production
iav_mosaic_prod <-
iavmosaic +
theme(text = element_text(size = 8),
legend.position = c(0.82, 0.84),
axis.text = element_text(size=8),
legend.key.size = unit(0.25, "cm"))
iav_mosaic_prod
iavbox_prod <-
data %>%
mutate(iav=ifelse(iav=="pos", "IAV+", "IAV-"),
iav=factor(iav, levels=c("IAV+", "IAV-"))) %>%
ggplot(., aes(x=iav, y=log(sumpcb))) +
geom_boxplot(fill = "#7eaaac", color = "black") +
labs(y="log(ΣPCB) (ng/g wet weight)") +
stat_n_text(y.pos=3, size = 2.5) +
theme_classic() +
theme(axis.title.x=element_blank(),
text = element_text(size = 8),
axis.text = element_text(size=8))
iavbox_prod## Warning: Removed 94 rows containing non-finite outside
## the scale range (`stat_boxplot()`).## Warning: Removed 94 rows containing non-finite outside
## the scale range (`stat_n_text()`).
## Warning: Removed 94 rows containing non-finite outside
## the scale range (`stat_boxplot()`).
## Removed 94 rows containing non-finite outside
## the scale range (`stat_n_text()`).