12 PCB ~ Exploratory Batch 1

2016 - 2019 Analysis Years (excluding 2022 & 2023)

12.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)

12.2 Data Import & Tidy

12.2.1 Data

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

maxlod <- 
  read.csv("Output Files/cleaned_pcb_maxLOD.csv") %>% 
  rename("Sample" = 1) %>% 
  filter(analysis.year == "2016" |
         analysis.year == "2017" |
         analysis.year == "2019")

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

12.2.2 Merge PCB data with metadata

• n = 127 samples

data <- 
  merge(maxlod,meta) %>% 
  mutate(year=as.numeric(year),
         molt.stage=as.factor(molt.stage),
         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)

12.2.3 Tidy data

Remove PCBs with 0 detections across all samples

data_tidy <- 
  data %>%
  column_to_rownames("Sample") %>% 
  dplyr::select(1:16) %>% 
  select_if(function (x) (sum(x) > 0)) %>% 
  rownames_to_column("Sample") %>% 
  merge(., data[,c(1, 18:36)], by = "Sample")

12.3 What influences PCBs?

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

12.3.1 Data

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

12.3.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

12.3.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   2.6059  1     0.1065

# Make predictions from model
logit_preds <- data.frame(stats::predict(yearmodel, type="link", se.fit=TRUE))
logit_preds$lwr <- logit_preds$fit + 1.96 * logit_preds$se.fit
logit_preds$upr <- logit_preds$fit - 1.96 * logit_preds$se.fit
real_preds <- apply(logit_preds, 2, invlogit)

# Combine with original data & plot
datapreds <- data.frame(data_tidy, real_preds)

yearlogplot <- 
  ggplot(datapreds, aes(x = year, y = fit)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.25) +
  geom_line() +
  labs(x="Year", y="Probality of PCBs") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        text=element_text(size=13))
yearlogplot

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"="gray60", "Absent"="lightgray"), name="PCB Status") +
  scale_y_continuous(labels = scales::percent) +
  labs(y="Proportion of Samples") +
  theme_classic() +
  theme(axis.title.x=element_blank(),
        text=element_text(size=13),
        axis.text.x=element_text(angle=45, hjust=1))
yearpresplot