8 PCB Dataset A

2016 data set, 19 pups, file name = 160134

8.1 Libraries

library(tidyverse)
library(gridExtra)
library(ggpubr)
library(car)
library(arm)
library(ggmosaic)
library(EnvStats)

8.2 Data

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

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

data <- 
  merge (pcb, meta) %>% 
  mutate(year = as.numeric(year),
         molt.stage = as.factor(molt.stage),
         iav = as.factor(iav),
         location = as.factor(location),
         sex = as.factor(sex))

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

8.3 Exploratory

8.3.1 Sample Description

table(data$year)

## 
## 2014 2015 2016 
##    5    7    7

table(data$sex)

## 
##  F  M 
##  8 10

table(data$location)

## 
##  Monomoy Muskeget 
##       14        5

table(data$molt.stage)

## 
## III  IV   V 
##   6   5   7

table(data$iav)

## 
## neg pos 
##  14   5

table(data$pcbbin)

## 
##  0  1 
##  3 16

8.3.2 Descriptive Stats

8.3.2.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 229.7375 32.43486  199.4 450 27.5

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

8.3.2.2 Congener Total

sumcongener <- 
  data.frame(colSums(data[,2:36])) %>% 
  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_A.jpg", sumcongenerplot, height=5, width=10, units="in")

8.3.2.3 Congener Mean

meancongener <- 
  data %>% 
  dplyr::select(2:36) %>% 
  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_A.jpg", meancongenerplot, height=5, width=8, units="in")

8.3.2.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: 105 × 3
## # Groups:   year [3]
##     year pcb       sum
##    <dbl> <chr>   <dbl>
##  1  2014 pcb.105   0  
##  2  2014 pcb.110   0  
##  3  2014 pcb.118  84.7
##  4  2014 pcb.128   0  
##  5  2014 pcb.130   0  
##  6  2014 pcb.138   0  
##  7  2014 pcb.146   0  
##  8  2014 pcb.149   0  
##  9  2014 pcb.151   0  
## 10  2014 pcb.153   0  
## # ℹ 95 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_A.jpeg", yearsumplot, height=5, width=7, units="in")

8.3.2.5 Year Mean

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

## # A tibble: 3 × 3
##    year  mean    se
##   <dbl> <dbl> <dbl>
## 1  2014  68.3  12.6
## 2  2015 150.   49.5
## 3  2016  55.3  14.4

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_A.jpeg", yearmeanplot, height=5, width=7, units="in")

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

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

8.3.3 Table of Mean, SE, Ranges for each congener

Gives some errors because some PCBs have no detections

summetrics <- 
  data.frame(
    LOD = "", 
    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",
    . = "NA")

lod <- 
  lod %>% 
  rbind(sumpcb)

pcbtable <-
  data.frame(matrix(nrow=35, ncol=0)) %>% 
  mutate(LOD=lod$.[1:35],
         sampsize=sapply(data[,2:36], function(x) length(x[x>0])),
         sampperc=sapply(data[,2:36], function(x) (length(x[x>0])/127)*100),
         mean=sapply(data[,2:36], function(x) {mean(x[x>0])}),
         se=sapply(data[,2:36], function(x) sd(x)/(sqrt(length(x[x>0])))),
         median=sapply(data[,2:36], function(x) {median(x[x>0])}),
         min=sapply(data[,2:36], function(x) {min(x[x>0])}),
         max=sapply(data[,2:36], 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, ")"))

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

8.4 What influences PCBs?

8.4.1 Data

# Remove pcbs with 0 detections
data_tidy <- 
  pcb %>% 
  column_to_rownames(var = "Sample") %>% 
  dplyr::select(1:35) %>% 
  select_if(function (x) (sum(x) > 0)) %>% 
  rownames_to_column(var = "Sample") %>% 
  merge(., data[,c(1,37:49)], by = "Sample")

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

8.4.2 Year

8.4.2.1 Presence/Absence

Not significant, but looks like proportion of pups with PCBs increases from 2014 - 2016 (limited sample sizes)

# Logistic regression
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   3.7342  1    0.05331 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 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
data.frame(data_tidy, real_preds) %>% 
  ggplot(., 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))

yearprop <-
  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(legend.position = "inside",
        legend.position.inside = c(0.83, 0.15),
        axis.title.x=element_blank(),
        text=element_text(size=13),
        axis.text.x=element_text(angle=45, hjust=1))
yearprop

8.4.2.2 Concentration

Not significant

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

## 
##  Pearson's product-moment correlation
## 
## data:  dataP$year and dataP$sumpcb
## t = 0.37733, df = 14, p-value = 0.7116
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4160632  0.5678005
## sample estimates:
##       cor 
## 0.1003359

# Plot
yearconc <-
  ggplot(dataP, aes(x=year, y=sumpcb)) +
  geom_point(size=0.5, position=position_jitter(width=0.2)) +
  geom_smooth(method = "lm", colour="black") +
  labs(y="log(ΣPCB ng/g wet weight)") +
  scale_x_continuous(breaks = seq(2014,2016, by=1)) +
  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))
yearconc

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

8.4.2.3 Combine

yearplots <- grid.arrange(yearprop, yearconc, ncol=2)

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

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

8.4.3 Sex

8.4.3.1 Presence/Absence

Not significant

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

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(data_tidy$sex, data_tidy$pcbbin)
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.008498059 13.719161956
## sample estimates:
## odds ratio 
##  0.5889506

8.4.3.2 Concentration

Not significant

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

## 
##  Welch Two Sample t-test
## 
## data:  sumpcb by sex
## t = 0.67796, df = 12.935, p-value = 0.5097
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -99.98219 191.36433
## sample estimates:
## mean in group F mean in group M 
##        262.6286        216.9375

8.4.4 Location

8.4.4.1 Presence/Absence

Not significant

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.1548
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.001829541 3.338779848
## sample estimates:
## odds ratio 
##  0.1345122

8.4.4.2 Concentration

Not significant

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

## 
##  Welch Two Sample t-test
## 
## data:  sumpcb by location
## t = 1.3336, df = 4.2969, p-value = 0.2486
## alternative hypothesis: true difference in means between group Monomoy and group Muskeget is not equal to 0
## 95 percent confidence interval:
##  -88.77244 261.74680
## sample estimates:
##  mean in group Monomoy mean in group Muskeget 
##               245.9538               159.4667

8.4.5 Molt Stage

8.4.5.1 Presence/Absence

Not significant

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 = 1
## alternative hypothesis: two.sided

8.4.5.2 Concentration

Significant - molt stage III has higher sumPCBs

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

##             Df Sum Sq Mean Sq F value Pr(>F)  
## molt.stage   2  95667   47834   3.965 0.0452 *
## Residuals   13 156817   12063                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova_molt <- aov(sumpcb ~ molt.stage, data = dataP)
TukeyHSD(anova_molt, "molt.stage")

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = sumpcb ~ molt.stage, data = dataP)
## 
## $molt.stage
##           diff       lwr        upr     p adj
## IV-III -177.52 -360.9335   5.893504 0.0582486
## V-III  -155.86 -331.4651  19.745060 0.0847512
## V-IV     21.66 -153.9451 197.265060 0.9434580

# Plot
ggplot(dataP, aes(x=molt.stage, y=sumpcb)) +
  geom_boxplot() +
  geom_point(size=0.5, position=position_jitter(width=0.2)) +
  labs(y="Sum PCBs") +
  stat_n_text() +
  theme_classic() +
  theme(axis.title.x=element_blank(), 
        text=element_text(size=13))

8.5 PCB Dataset A - Cytokines

n = 17 pups, 2016 analysis year

8.6 Libraries

library(mixOmics)
library(tidyverse)
library(ggfortify)
library(vegan)
library(car)
library(EnvStats)
library(gridExtra)
library(ggpubr)
library(mlr3)
library(mlr3tuning)
library(rpart)
library(rpart.plot)
library(grid)
library(ggmosaic)
library(kableExtra)

8.7 Data

pcb <- 
  read.csv("Output Files/cleaned_pcb_A.csv") %>% 
  rename("Sample" = c(1)) %>% 
  mutate(pcbbin = ifelse(pcbbin=="0", "PCB Absent", "PCB Present"))

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

data <- 
  read.csv("../hg-cyto/Output Files/cleaned_cytokine_all.csv") %>% 
  filter(Sample %in% pcb$Sample) %>% 
  merge(., pcb, by="Sample") %>% 
  merge(., meta, by = "Sample")

databin <- 
  read.csv("../hg-cyto/Output Files/cleaned_cytokine_bin_all.csv") %>% 
  filter(Sample %in% pcb$Sample) %>%  
  merge(., pcb, by="Sample") %>% 
  merge(., meta, by = "Sample")

8.7.1 Sample Description

table(data$year)

## 
## 2014 2015 2016 
##    3    7    7

table(data$sex)

## 
## F M 
## 8 8

table(data$location)

## 
##  Monomoy Muskeget 
##       14        3

table(data$molt.stage)

## 
## III  IV   V 
##   5   5   7

table(data$iav)

## 
## neg pos 
##  13   4

table(data$pcbbin)

## 
##  PCB Absent PCB Present 
##           1          16

table(data$analysis.year.x)

## 
## 2016 
##   17

Only 1 pup with PCBs - not continuing with this analysis.

8.8 PCB Dataset A- IAV

n = 19 pups

8.9 Load required libraries

library(tidyverse)
library(car)
library(ggmosaic)
library(gridExtra)
library(EnvStats)
library(rstatix)

8.10 Data

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

pcb <- 
  read.csv("Output Files/cleaned_pcb_A.csv") %>% 
  dplyr::rename("Sample" = 1)

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% pcb$Sample) %>% 
  rbind(., viro_two)

8.11 Merge data & configure variables

data <- 
  merge(meta, pcb, by="Sample") %>% 
  select(Sample, iav, sumpcb, pcbbin, year, location, sex, molt.stage) %>% 
  mutate(iavbin = as.numeric(ifelse(iav=="neg", "0", "1")))

swabdata <- 
  merge(viro_tidy, pcb, by="Sample")

8.12 IAV ~ Variables

8.12.1 Year

Not significant

yeariav <- glm(iavbin ~ year, data=data, family=binomial) 
Anova(yeariav)

## Analysis of Deviance Table (Type II tests)
## 
## Response: iavbin
##      LR Chisq Df Pr(>Chisq)
## year 0.098841  1     0.7532

8.12.2 Sex

Not significant

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

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(data$sex, data$iavbin)
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.1042587 20.1438280
## sample estimates:
## odds ratio 
##   1.267987

8.12.3 Location

Not significant

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 = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.0101872 9.9071291
## sample estimates:
## odds ratio 
##  0.6397607

8.12.4 Molt Stage

Not significant

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.424
## alternative hypothesis: two.sided

8.13 IAV ~ PCB p/a

Not significant

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

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(data$iav, data$pcbbin)
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.02763137 49.05655211
## sample estimates:
## odds ratio 
##  0.6820379

8.14 IAV ~ PCB conc

Significant

# Fit model
model <- glm(iavbin ~ sumpcb, data=data, family=binomial)
Anova(model)

## Analysis of Deviance Table (Type II tests)
## 
## Response: iavbin
##        LR Chisq Df Pr(>Chisq)  
## sumpcb   3.9712  1    0.04628 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.15 Figures

8.15.1 Presence/Absence

IAV+ pups less likely to have PCBs present (barely)

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"="lightgray", "Present"="gray60"), 
                    name = "PCB Status") +
  scale_y_continuous(labels = scales::percent) +
  labs(y = "Proportion of Samples") +  
  theme_classic() +
  theme(legend.position=c(0.83, 0.15),
        axis.title.x = element_blank(),
        text = element_text(size=13),
        axis.text = element_text(size=13))

8.15.2 Concentration

.. but have greater concentrations of PCBs (barely) Only 4 pups with IAV in this dataset.

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 () +
  geom_point(position = position_jitter(width = 0.2)) +
  labs(y="log(ΣPCB Concentration)") +
  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))

8.15.3 Combine

pcbiav <- grid.arrange(iavmosaic, iavbox, ncol=2)

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

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

ggsave("Figures/pcb-iav_A.jpeg", pcbiav, width=10, height=5, units="in")