4 PCB ~ Exploratory

• Max LOD = combine all datasets, use PCBs only common to all datasets, and highest LOD = threshold for PCB detections

4.1 Libraries

library(ggplot2)
library(reshape2)
library(ggpubr)
library(gridExtra)
library(rstatix)
library(ggmosaic)
library(DescTools)
library(EnvStats)
library(car)
library(arm)
library(tidyverse)
library(vegan)

4.2 Data Import & Tidy

4.2.1 Data

meta <- 
  read.csv("Input Files/metadata_tidy.csv") 

maxlod <- 
  read.csv("Output Files/cleaned_pcb_maxLOD.csv") %>% 
  rename("Sample" = 1)

lod <- 
  read.csv("Output Files/lod_all.csv")

cyto <- 
  read.csv("../hg-cyto/Output Files/cleaned_cytokine_all.csv")

rna <- 
  read.csv("../hg-rna/Input Files/metadata.csv") %>% 
  filter(!sex=="male") %>% 
  mutate(sample=gsub("Hg","", sample),
         sample=gsub("a", "", sample),
         sample=gsub("b", "", sample),
         sample=gsub("11", "1", sample)) 

whole <- 
  read.csv("Output Files/cleaned_pcb2_wholeblood.csv")

4.2.2 Merge PCB data with metadata

• n = 127 samples

data <- 
  merge(maxlod, meta) %>% 
  mutate(year=as.numeric(year),
         molt.stage=factor(molt.stage, levels=c("II", "III", "IV", "V", "U")),
         iav=as.factor(iav),
         location=as.factor(location),
         sex=as.factor(sex),
         analysis.year=as.factor(analysis.year)) %>% 
  mutate_at(c(2:17), as.numeric)

4.3 Exploratory

4.3.1 Sample description

table(data$year)

## 
## 2014 2015 2016 2017 2018 2019 2020 
##    5    7    9   51   20   25   10

table(data$sex)

## 
##  F  M 
## 70 50

table(data$location)

## 
## Great Point     Monomoy    Muskeget 
##           6          81          40

table(data$molt.stage)

## 
##  II III  IV   V   U 
##   1  15  23  67   0

table(data$iav)

## 
## neg pos 
##  86  41

table(data$analysis.year)

## 
## 2016 2017 2019 2022 2023 
##   19   52   20   33    3

yearxanalysis<-
  table(data$year, data$analysis.year) %>%
  data.frame() %>% 
  pivot_wider(names_from = "Var2", values_from = "Freq") %>% 
  setNames(c("Year", "2016 Analysis", "2017 Analysis", "2019 Analysis", "2022 Analysis", "2023 Analysis"))

4.3.2 Supplemental table

Sample attributes

table <-
  data %>% 
  select(Sample, iav, location, year,  sex, molt.stage, analysis.year) %>% 
  arrange(year) %>% 
  mutate(batch = ifelse(analysis.year == "2016", "A", 
                   ifelse(analysis.year == "2017", "B", 
                     ifelse(analysis.year == "2019", "B", "C"))),
         PCBs = ifelse(Sample %in% whole$X, "x (P)", "x"),
         Cytokines= ifelse(Sample %in% cyto$Sample, "x", ""),
         RNA = ifelse(Sample %in% rna$sample, "x", "")) %>% 
  select(!analysis.year)

write.csv(table, "Output Files/supplemental_samplemetadata.csv")

4.3.3 Descriptive Stats

4.3.3.1 Mean Median, Range

data %>% 
  filter(sumpcb>0) %>% 
  summarise(mean=mean(sumpcb), 
            se=sd(sumpcb)/(sqrt(length(sumpcb[sumpcb>0]))),
            median=median(sumpcb),
            max=max(sumpcb),
            min=min(sumpcb))

##       mean       se median max  min
## 1 118.7485 17.66918   81.2 445 27.9

datalong <- 
  data %>% 
  pivot_longer(2:17, names_to = "pcb", values_to = "conc")

4.3.3.2 Congener Total

sumcongener <- 
  data.frame(colSums(data[,2:17])) %>% 
  mutate(congener=factor(row.names(.), levels=row.names(.))) %>% 
  rename(sum=c(1)) 

sumcongenerplot <- 
  ggplot(sumcongener, aes(x=congener, y=sum)) +
  geom_bar(stat="identity") + 
  labs(x="PCB Congener", y="Sum PCBs (ng/g wet weight)") +
  theme_classic() +
  theme(axis.text.x=element_text(angle=90, hjust=1), 
        text=element_text(size=12))
sumcongenerplot

ggsave("Figures/sumpcb.jpg", sumcongenerplot, height=5, width=7, units="in")

4.3.3.3 Congener Mean

meancongener <- 
  data %>% 
  dplyr::select(2:17) %>% 
  pivot_longer(everything(), names_to="pcb", values_to="conc") %>%
  filter(!conc==0) %>% 
  group_by(pcb) %>% 
  summarize(meanconc=mean(conc),
            se=sd(conc)/(sqrt(length(conc)))) 

meancongenerplot <-
  ggplot(meancongener, aes(x=pcb, y=meanconc)) +
  geom_bar(stat="identity") + 
  geom_errorbar(aes(ymin = meanconc - se, ymax = meanconc + se), width = 0.2) +
  labs(x="PCB Congener", y="Mean PCBs (ng/g wet weight)") +
  theme_classic() +
  theme(axis.text.x=element_text(angle=90, hjust=1), 
        text=element_text(size=12))
meancongenerplot

ggsave("Figures/meanpcb.jpg", meancongenerplot, height=5, width=7, units="in")

4.3.3.4 Year Total

yearsumcongener <-
  datalong %>% 
  group_by(year, pcb) %>% 
  summarize(sum = sum(conc))

## `summarise()` has grouped output by 'year'. You
## can override using the `.groups` argument.

yearsumcongener

## # A tibble: 112 × 3
## # Groups:   year [7]
##     year pcb      sum
##    <dbl> <chr>  <dbl>
##  1  2014 pcb105   0  
##  2  2014 pcb118  84.7
##  3  2014 pcb128   0  
##  4  2014 pcb138 133. 
##  5  2014 pcb153   0  
##  6  2014 pcb18    0  
##  7  2014 pcb180   0  
##  8  2014 pcb183   0  
##  9  2014 pcb184   0  
## 10  2014 pcb187   0  
## # ℹ 102 more rows

ggplot(yearsumcongener, aes(x=pcb, y=sum)) +
       geom_bar(stat="identity") + 
       labs(x="PCB Congener", y="sum PCBs (ng/g wet weight)") +
       facet_wrap(~year, ncol = 4) +
       theme(axis.text.x = element_text(angle=90))

yearsum <- 
  yearsumcongener %>% 
  group_by(year) %>% 
  summarise(sum=sum(sum)) 

yearsumplot <-
  ggplot(yearsum, aes(x=factor(year), y=sum)) +
  geom_bar(stat="identity") + 
  labs(x="Year", y="sum PCBs (ng/g wet weight)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90),
        panel.grid = element_blank())
yearsumplot

ggsave("Figures/yearsumpcb.jpeg", yearsumplot, height=5, width=7, units="in")

4.3.3.5 Year Mean

yearmean <-
  datalong %>% 
  filter(!conc==0) %>% 
  group_by(year) %>% 
  summarize(mean = mean(conc),
            se = sd(conc)/(sqrt(length(conc)))) 
yearmean

## # A tibble: 6 × 3
##    year  mean     se
##   <dbl> <dbl>  <dbl>
## 1  2014 109.   23.9 
## 2  2015 293.  109.  
## 3  2016  45.8   1.93
## 4  2017  74.7   5.97
## 5  2018  52    NA   
## 6  2019 103.    7.99

yearmeanplot <- 
  ggplot(yearmean, aes(x=factor(year), y=mean)) +
  geom_bar(stat="identity") + 
  geom_errorbar(aes(ymin = mean-se, ymax = mean+se), 
                width = 0.2) +
  labs(x="PCB Congener", y="Mean PCBs (ng/g wet weight)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90),
        panel.grid = element_blank())
yearmeanplot

ggsave("Figures/yearmeanpcb.jpeg", yearmeanplot, height=5, width=7, units="in")

4.3.3.6 What congeners are present in each animal (separated by year/iav status)?

Additive/synergistic effects - is # congeners associated with IAV?

ggplot(data=datalong, aes(x=Sample, y=conc, fill=pcb)) +
  geom_col() +
  labs(x="Seal ID", y="PCB Concentration") +
  facet_wrap(~iav, drop=TRUE, scales="free", ncol=1) +
  theme_bw() +
  theme(panel.grid=element_blank(), 
        axis.text.x=element_text(angle=90, vjust=0.5, hjust=1))

4.3.3.7 How many seals were each congener detected in?

datalong %>%
  group_by(pcb) %>%
  summarise(count=sum(conc>0)) %>% 
  ggplot(data=., aes(x=pcb, y=count)) +
  geom_col() +
  labs(x="Congener", y="Number of Seals") + 
  geom_text(aes(label=count, vjust=-0.75)) +
  theme_bw() +
  theme(panel.grid=element_blank(), axis.text.x=element_text(angle=90, vjust=0.5, hjust=1))

4.3.4 Table of Mean, SE, Ranges for each congener

Gives some errors because some PCBs have no detections

summetrics <- 
  data.frame(
    maxLOD = "", 
    sampsize = table(data$pcbbin)[2],
    sampperc = (table(data$pcbbin)[2]/nrow(data))*100,
    mean = mean(data$sumpcb[data$sumpcb>0]),
    se = sd(data$sumpcb)/(sqrt(length(data$sumpcb[data$sumpcb>0]))),
    median = median(data$sumpcb[data$sumpcb>0]),
    min = min(data$sumpcb[data$sumpcb>0]),
    max = max(data$sumpcb[data$sumpcb>0])
  )

sumpcb <- 
  data.frame(
    X = "sumPCB",
    X1 = "NA",
    X2 = "NA", 
    X3 = "NA", 
    X4 = "NA", 
    X5 = "NA", 
    max = "NA")

lod <- 
  lod %>% 
  rbind(sumpcb)
  
pcbtable <-
  data.frame(matrix(nrow=16, ncol=0)) %>% 
  mutate(maxLOD=lod$max[1:16],
         sampsize=sapply(data[,2:17], function(x) length(x[x>0])),
         sampperc=sapply(data[,2:17], function(x) (length(x[x>0])/127)*100),
         mean=sapply(data[,2:17], function(x) {mean(x[x>0])}),
         se=sapply(data[,2:17], function(x) sd(x)/(sqrt(length(x[x>0])))),
         median=sapply(data[,2:17], function(x) {median(x[x>0])}),
         min=sapply(data[,2:17], function(x) {min(x[x>0])}),
         max=sapply(data[,2:17], max)) %>% 
  rbind(summetrics) %>%
  mutate(across(everything(), as.numeric)) %>% 
  sapply(., function(x) round(x, digits=2)) %>% 
  data.frame() %>% 
  add_column(., pcb=lod$X, .before = 1) %>% 
  mutate(pcb=gsub("pcb", "PCB ", pcb),
         sampnum=paste0(sampsize, " (", sampperc, "%)"),
         meanse=paste0(mean, " ± ", se),
         medrange=paste0(median, " (", min, " - ", max, ")")) %>% 
  rbind()

write.csv(pcbtable, "Output Files/pcb_summarystats.csv", row.names=FALSE)

4.4 Tidy data

Remove PCBs with 0 detections across all samples

data_tidy <- 
  data %>% 
  dplyr::select(2:17) %>% 
  select_if(function (x) (sum(x) > 0)) %>% 
  mutate(Sample=data$Sample,
         sumpcb=data$sumpcb, 
         pcbbin=data$pcbbin,
         year=data$year,
         sex=data$sex,
         location=data$location,
         molt.stage=data$molt.stage,
         iav=data$iav,
         iavser=data$iavser,
         bc_bin=data$bc_bin,
         analysis.year=data$analysis.year)

4.5 Batch Effects?

4.5.1 Data

pcb_pca <- 
  data_tidy %>% 
  mutate(across(.cols=c(1:13), function(x) log(x+0.01)))

pcb_dist <- 
  pcb_pca %>% 
  column_to_rownames("Sample") %>% 
  select(1:13) %>% 
  vegdist(., method="euclidean")

pcb_remove0 <- 
  data_tidy %>% 
  filter(!sumpcb==0)

4.5.2 PCA

pca <- 
  prcomp(pcb_pca[,1:13], center=TRUE, scale=TRUE)

pca_df <- 
  data.frame(
    x=pca$x[,1],
    y=pca$x[,2],
    molt.stage=pcb_pca$molt.stage,
    year=factor(pcb_pca$year),
    location=pcb_pca$location,
    sex=pcb_pca$sex,
    iav=pcb_pca$iav,
    iavser=pcb_pca$iavser,
    analysis.year=pcb_pca$analysis.year
  )

ggplot(pca_df, 
       aes(x = x, 
           y = y,
           col = analysis.year)) +
  geom_point() +
  stat_ellipse() +
  theme_bw() +
  theme(panel.grid=element_blank())

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning: Computation failed in `stat_ellipse()`.
## Caused by error in `if (all(abs(w - w0) < tol)) ...`:
## ! missing value where TRUE/FALSE needed

4.5.3 PCO

Dissimilarity matrix to help deal with zero-inflated data

mds <- cmdscale(pcb_dist, eig=TRUE)

mds_df <-
  mds$points %>% 
  data.frame() %>% 
  rownames_to_column("Sample") %>% 
  merge(., data_tidy, by="Sample")

ggplot(mds_df, 
       aes(x=X1, 
           y=X2,
           col=analysis.year)) +
  geom_point() +
  stat_ellipse() +
  theme_bw() +
  theme(panel.grid=element_blank())

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Too few points to calculate an ellipse

## Warning: Removed 1 row containing missing values or
## values outside the scale range (`geom_path()`).

4.5.4 PCA without 0’s

pca2 <- prcomp(pcb_remove0[,1:13], center=TRUE, scale=TRUE)

pca_df2 <- 
  data.frame(
    x=pca2$x[,1],
    y=pca2$x[,2],
    molt.stage=pcb_remove0$molt.stage,
    year=factor(pcb_remove0$year),
    location=pcb_remove0$location,
    sex=pcb_remove0$sex,
    iav=pcb_remove0$iav,
    iavser=pcb_remove0$iavser,
    analysis.year=pcb_remove0$analysis.year
  )

ggplot(pca_df2, 
       aes(x = x, 
           y = y,
           col = analysis.year)) +
  geom_point() +
  stat_ellipse() +
  theme_bw() +
  theme(panel.grid=element_blank())

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Too few points to calculate an ellipse

## Warning: Removed 1 row containing missing values or
## values outside the scale range (`geom_path()`).

4.6 What influences PCBs?

Use Hurdle Model framework - first presence/absence of PCBs and then use sum of non-zero PCBs.

4.6.1 Data

# Subset pups with PCBs Present
dataP <- 
  data_tidy %>% 
  filter(sumpcb > 0) %>% 
  mutate(logpcb=log(sumpcb))

4.6.1.1 Are the non-zero data normal?

• Raw data
  • not normally distributed
• Log-transformed data
  • not significantly different than a normal distribution!

# Raw data
hist <- ggplot(dataP, aes(x=sumpcb)) + 
 geom_histogram()
dens <- ggdensity(dataP$sumpcb)
qq <- ggqqplot(dataP$sumpcb)
grid.arrange(hist, dens, qq, ncol=3)

shapiro.test(dataP$sumpcb)

## 
##  Shapiro-Wilk normality test
## 
## data:  dataP$sumpcb
## W = 0.7037, p-value = 7.433e-07

# Log-transformed
hist <- ggplot(dataP, aes(x=logpcb)) + 
 geom_histogram()
dens <- ggdensity(dataP$logpcb)
qq <- ggqqplot(dataP$logpcb)
grid.arrange(hist, dens, qq, ncol=3)

shapiro.test(dataP$logpcb)

## 
##  Shapiro-Wilk normality test
## 
## data:  dataP$logpcb
## W = 0.9387, p-value = 0.06229

4.6.1.2 Year - Presence/Absence

yearmodel <- glm(pcbbin ~ year, data=data_tidy, family=binomial) 
Anova(yearmodel)

## Analysis of Deviance Table (Type II tests)
## 
## Response: pcbbin
##      LR Chisq Df Pr(>Chisq)    
## year   11.769  1  0.0006021 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(yearmodel)

## 
## Call:
## glm(formula = pcbbin ~ year, family = binomial, data = data_tidy)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) 1021.8055   318.2875   3.210  0.00133 **
## year          -0.5070     0.1578  -3.213  0.00131 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 145.51  on 126  degrees of freedom
## Residual deviance: 133.75  on 125  degrees of freedom
## AIC: 137.75
## 
## Number of Fisher Scoring iterations: 4

yearpresplot <-
  data_tidy %>% 
  mutate(pcbbin=ifelse(pcbbin=="1", "Present", "Absent"),
         pcbbin=factor(pcbbin, levels=c("Present", "Absent"))) %>% 
  ggplot(data=.) +
  geom_mosaic(aes(x=product(pcbbin, year), fill=pcbbin)) +
  scale_fill_manual(values=c("Present"="#7eaaac", "Absent"="#b4b4b4"), name="PCB Status") +
  scale_y_continuous(labels = scales::percent) +
  labs(y="Proportion of Samples") +
  theme_classic() +
  theme(legend.position=c(0.83, 0.86), 
        axis.title.x=element_blank(),
        text=element_text(size=13),
        axis.text.x=element_text(angle=45, hjust=1))
yearpresplot

4.6.1.3 Year - Concentration

# Pearson's Correlation
cor.test(dataP$year, dataP$logpcb, method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  dataP$year and dataP$logpcb
## t = -0.28919, df = 31, p-value = 0.7744
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3882654  0.2967224
## sample estimates:
##         cor 
## -0.05187053

yearconcplot <- 
  ggplot(dataP, aes(x=year, y=log(sumpcb))) +
  geom_point(size=2, position=position_jitter(width=0.2)) +
  geom_smooth(method = "lm", colour="#7eaaac") +
  labs(y="log(ΣPCB) (ng/g wet weight)") +
  stat_n_text() +
  theme_classic() +
  theme(axis.title.x=element_blank(), 
        text=element_text(size=13),
        axis.text.x=element_text(angle=45, hjust=1))
yearconcplot

## `geom_smooth()` using formula = 'y ~ x'

yearplots <- grid.arrange(yearpresplot, yearconcplot, ncol=2)

## `geom_smooth()` using formula = 'y ~ x'

ggsave("Figures/year-pcbbin-sumpcb.jpeg", yearplots, width=12, height=6, units="in")

4.6.1.3.1 Year Presentation Plots

year_pres_pres <- 
  data_tidy %>% 
  group_by(year) %>% 
  summarize(sum = sum(pcbbin)) %>% 
  mutate(yeartotal = table(data_tidy$year),
         prop = sum/yeartotal,
         year = factor(year))

prop_pres <-
  ggplot(year_pres_pres, aes(x = year, y = prop, fill = year)) +
  geom_col(position = "identity") +
  scale_fill_manual(values = c("2014" = "#75B6E5", 
                               "2015" = "#B5CDD3",
                               "2016" = "#AFBABB",
                               "2017" = "#ACD7CA",
                               "2018" = "#9BD3D9",
                               "2019" = "#7FC1DB",
                               "2020" = "#7F7B9A")) +
  scale_y_continuous(labels = scales::percent) +
  labs(x = NULL, y = "Proportion of Samples") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        text = element_text(size = 16),
        legend.position = "none")
prop_pres

conc_pres <- 
  data_tidy %>% 
  mutate(year = factor(year)) %>% 
  ggplot(., aes(x = year, y = log(sumpcb), fill = year)) +
  geom_boxplot() +
  scale_fill_manual(values = c("2014" = "#75B6E5", 
                               "2015" = "#B5CDD3",
                               "2016" = "#AFBABB",
                               "2017" = "#ACD7CA",
                               "2018" = "#9BD3D9",
                               "2019" = "#7FC1DB",
                               "2020" = "#7F7B9A")) +
  labs(x = NULL, y = "log(Sum PCB Concentration)") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        text = element_text(size = 16),
        legend.position = "none")
conc_pres

## Warning: Removed 94 rows containing non-finite outside
## the scale range (`stat_boxplot()`).

year_pres <- 
  grid.arrange(prop_pres, conc_pres, ncol = 2)

## Warning: Removed 94 rows containing non-finite outside
## the scale range (`stat_boxplot()`).

ggsave("Figures/year-pcbbin-sumpcb_pres.jpeg", year_pres, 
       width = 10, height = 5, units = "in")

4.6.1.3.2 Year Production Plots

yearpresplot_prod <-
  yearpresplot +
  theme(text = element_text(size = 9),
        legend.position = c(0.81, 0.85),
        legend.key.size = unit(0.25, "cm"))

yearconcplot_prod <-
  ggplot(dataP, aes(x=year, y=log(sumpcb))) +
  geom_point(size=1, position=position_jitter(width=0.2)) +
  geom_smooth(method = "lm", colour="#7eaaac") +
  labs(y="log(ΣPCB) (ng/g wet weight)") +
  stat_n_text(size = 2) +
  theme_classic() +
  theme(axis.title.x=element_blank(), 
        text=element_text(size=9),
        axis.text.x=element_text(angle=45, hjust=1))

production <- grid.arrange(yearpresplot_prod, yearconcplot_prod, ncol=2)

## `geom_smooth()` using formula = 'y ~ x'

ggsave("Figures/Figure_2.jpeg", production, 
        width = 7480, height = 3500, units = "px", dpi = 1000)

4.6.1.4 Sex - Presence/Absence

tryCatch({
  chisq.test(table(data_tidy$sex, data_tidy$pcbbin))
}, warning = function(w) {
  fisher.test(table(data_tidy$sex, data_tidy$pcbbin))
})

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(data_tidy$sex, data_tidy$pcbbin)
## X-squared = 0.1331, df = 1, p-value = 0.7152

4.6.1.5 Sex - Concentration

# Concentration
t.test(logpcb ~ sex, data=dataP)

## 
##  Welch Two Sample t-test
## 
## data:  logpcb by sex
## t = -1.3245, df = 19.709, p-value = 0.2005
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -0.9152528  0.2047721
## sample estimates:
## mean in group F mean in group M 
##        4.393187        4.748427

ggplot(dataP, aes(x = sex, y = sumpcb)) +
  geom_boxplot() +
  stat_n_text()

4.6.1.6 Location - Presence/Absence

tryCatch({
  chisq.test(table(data_tidy$location, data_tidy$pcbbin))
}, warning = function(w) {
  fisher.test(table(data_tidy$location, data_tidy$pcbbin))
})

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(data_tidy$location, data_tidy$pcbbin)
## p-value = 0.4352
## alternative hypothesis: two.sided

4.6.1.7 Location - Concentration

t.test(logpcb ~ location, data=dataP)

## 
##  Welch Two Sample t-test
## 
## data:  logpcb by location
## t = -0.87241, df = 23.415, p-value = 0.3918
## alternative hypothesis: true difference in means between group Monomoy and group Muskeget is not equal to 0
## 95 percent confidence interval:
##  -0.6442122  0.2617620
## sample estimates:
##  mean in group Monomoy mean in group Muskeget 
##               4.479841               4.671066

4.6.1.8 Molt Stage - Presence/Absence

tryCatch({
  chisq.test(table(data_tidy$molt.stage, data_tidy$pcbbin))
}, warning = function(w) {
  fisher.test(table(data_tidy$molt.stage, data_tidy$pcbbin))
})

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(data_tidy$molt.stage, data_tidy$pcbbin)
## p-value = 0.03023
## alternative hypothesis: two.sided

pairwise_fisher_test(table(data_tidy$molt.stage, data_tidy$pcbbin))

## # A tibble: 10 × 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 0.292  1     ns          
##  3 II     V         68 0.191  1     ns          
##  4 II     U          1 1      1     ns          
##  5 III    IV        38 0.295  1     ns          
##  6 III    V         82 0.0367 0.367 ns          
##  7 III    U         15 1      1     ns          
##  8 IV     V         90 0.384  1     ns          
##  9 IV     U         23 1      1     ns          
## 10 V      U         67 1      1     ns

4.6.1.9 Molt Stage - Concnetration

summary(aov(logpcb ~ molt.stage, data=dataP))

##             Df Sum Sq Mean Sq F value Pr(>F)
## molt.stage   3  1.848  0.6159   1.446  0.257
## Residuals   22  9.370  0.4259               
## 7 observations deleted due to missingness

4.6.2 Batch - Presence/Absence

Sig: 2022 vs 2016 & 2017

tryCatch({
  chisq.test(table(data_tidy$analysis.year, data_tidy$pcbbin))
}, warning = function(w) {
  fisher.test(table(data_tidy$analysis.year, data_tidy$pcbbin))
})

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(data_tidy$analysis.year, data_tidy$pcbbin)
## p-value = 5.841e-06
## alternative hypothesis: two.sided

pairwise_fisher_test(table(data_tidy$analysis.year, data_tidy$pcbbin))

## # A tibble: 10 × 6
##    group1 group2     n          p     p.adj p.adj.signif
##  * <chr>  <chr>  <int>      <dbl>     <dbl> <chr>       
##  1 2016   2017      71 1          1         ns          
##  2 2016   2019      39 0.0824     0.494     ns          
##  3 2016   2022      52 0.0001     0.0009    ***         
##  4 2016   2023      22 0.273      1         ns          
##  5 2017   2019      72 0.0511     0.39      ns          
##  6 2017   2022      85 0.00000229 0.0000229 ****        
##  7 2017   2023      55 0.267      1         ns          
##  8 2019   2022      53 0.0487     0.39      ns          
##  9 2019   2023      23 1          1         ns          
## 10 2022   2023      36 1          1         ns

table(data_tidy$pcbbin, data_tidy$analysis.year)

##    
##     2016 2017 2019 2022 2023
##   0   11   30   17   33    3
##   1    8   22    3    0    0

batch_pres_plot <-
  data_tidy %>% 
  mutate(pcbbin=ifelse(pcbbin=="1", "Present", "Absent"),
         pcbbin=factor(pcbbin, levels=c("Present", "Absent"))) %>% 
  ggplot(data=.) +
  geom_mosaic(aes(x=product(pcbbin, analysis.year), fill=pcbbin)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray"), 
                    name="PCB Status") +
  scale_y_continuous(labels = scales::percent) +
  labs(y="Proportion of Samples") +
  theme_classic() +
  theme(legend.position=c(0.83, 0.86), 
        axis.title.x=element_blank(),
        text=element_text(size=13),
        axis.text.x=element_text(angle=45, hjust=1))
batch_pres_plot

4.6.3 Batch - Concentration

Not significant

summary(aov(logpcb ~ analysis.year, data=dataP))

##               Df Sum Sq Mean Sq F value Pr(>F)
## analysis.year  2  1.366  0.6828   1.674  0.204
## Residuals     30 12.234  0.4078