10 PCB Dataset C
2022 & 2023 data sets, 36 pups, file name = 220264, 230045
10.2 Data
pcb <-
read.csv("Output Files/cleaned_pcb_C.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_C.csv")10.3 Exploratory
10.3.1 Sample Description
##
## 2016 2018 2019 2020
## 2 19 5 10##
## F M
## 24 12##
## Great Point Monomoy Muskeget
## 4 21 11##
## IV V
## 4 29##
## neg pos
## 26 10##
## 0 1
## 18 1810.3.2 Descriptive Stats
10.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 6.922222 1.019982 5.8 17.8 1.910.3.2.2 Congener Total
sumcongener <-
data.frame(colSums(data[,2:22])) %>%
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
10.3.2.3 Congener Mean
meancongener <-
data %>%
dplyr::select(2:22) %>%
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
10.3.2.4 Year Total
## `summarise()` has grouped output by 'year'. You
## can override using the `.groups` argument.## # A tibble: 84 × 3
## # Groups: year [4]
## year pcb sum
## <dbl> <chr> <dbl>
## 1 2016 pcb101 0
## 2 2016 pcb105 0
## 3 2016 pcb118 0
## 4 2016 pcb128 0
## 5 2016 pcb138 0
## 6 2016 pcb153 0
## 7 2016 pcb170 0
## 8 2016 pcb18 0
## 9 2016 pcb180 0
## 10 2016 pcb183 0
## # ℹ 74 more rowsggplot(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
10.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: 4 × 3
## year mean se
## <dbl> <dbl> <dbl>
## 1 2016 17.8 NA
## 2 2018 4.71 0.483
## 3 2019 4.18 0.749
## 4 2020 4.82 2.31yearmeanplot <-
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
10.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))
10.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))
10.3.3 Table of Mean, SE, Ranges for each congener
Gives some errors becauause some PCBs have no detections
summetrics <-
data.frame(
LOD = "NA",
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")
lod <-
lod %>%
rbind(sumpcb)
pcbtable <-
data.frame(matrix(nrow=21, ncol=0)) %>%
mutate(LOD=lod$X1[1:21],
sampsize=sapply(data[,2:22], function(x) length(x[x>0])),
sampperc=sapply(data[,2:22], function(x) (length(x[x>0])/127)*100),
mean=sapply(data[,2:22], function(x) {mean(x[x>0])}),
se=sapply(data[,2:22], function(x) sd(x)/(sqrt(length(x[x>0])))),
median=sapply(data[,2:22], function(x) {median(x[x>0])}),
min=sapply(data[,2:22], function(x) {min(x[x>0])}),
max=sapply(data[,2:22], 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(sampnum=paste0(sampsize, " (", sampperc, "%)"),
meanse=paste0(mean, " ± ", se),
medrange=paste0(median, " (", min, " - ", max, ")"))
write.csv(pcbtable, "Output Files/pcb_summarystats_C.csv", row.names=FALSE)10.4 What influences PCBs?
10.4.1 Data
# Remove pcbs with 0 detections
data_tidy <-
pcb %>%
column_to_rownames(var = "Sample") %>%
dplyr::select(1:21) %>%
select_if(function (x) (sum(x) > 0)) %>%
rownames_to_column(var = "Sample") %>%
merge(., data[,c(1,23:35)], by = "Sample")
# Subset only pups with PCBs
dataP <-
data_tidy %>%
filter(sumpcb > 0) %>%
mutate(logpcb=log(sumpcb))10.4.2 Year
10.4.2.1 Presence/Absence
Not significant - some decline from 2016 to 2020 but not steep
# 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 0.61853 1 0.4316# 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
10.4.2.2 Concentration
Not significant
##
## Pearson's product-moment correlation
##
## data: dataP$year and dataP$sumpcb
## t = -1.4477, df = 16, p-value = 0.167
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6965216 0.1504577
## sample estimates:
## cor
## -0.3403178# 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(2016, 2020, 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'
10.4.3 Sex
10.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))
})##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(data_tidy$sex, data_tidy$pcbbin)
## X-squared = 1.125, df = 1, p-value = 0.288810.4.3.2 Concentration
Significant - males lower sumpcbs than females, but low male sample size
##
## Welch Two Sample t-test
##
## data: sumpcb by sex
## t = 2.5407, df = 12.31, p-value = 0.02548
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## 0.5674322 7.2682821
## sample estimates:
## mean in group F mean in group M
## 7.792857 3.875000
10.4.4 Location
10.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.2903
## alternative hypothesis: two.sided10.4.5 Molt Stage
10.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 = 0.6012
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.2193726 179.5832915
## sample estimates:
## odds ratio
## 3.10924710.7 Data
pcb <-
read.csv("Output Files/cleaned_pcb_C.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")10.8 Sample Description
##
## 2016 2018 2019 2020
## 2 16 5 10##
## F M
## 21 12##
## Great Point Monomoy Muskeget
## 4 18 11##
## IV V
## 4 26##
## neg pos
## 23 10##
## PCB Absent PCB Present
## 18 15##
## 2022
## 3310.9 PCB - IAV interaction
Not significant
##
## PCB Absent PCB Present
## neg 13 10
## pos 5 5tryCatch({
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.2258772 7.4239096
## sample estimates:
## odds ratio
## 1.28962710.10 Cytokine Exploratory
10.10.1 Cytokine detections
detect <-
data %>%
pivot_longer(., cols=c(2:14), names_to="cytokine", values_to="conc") %>%
group_by(cytokine) %>%
summarize(n = sum(conc > 0)) %>%
arrange(desc(n))
detect %>%
kable() %>%
kable_styling("basic")| cytokine | n |
|---|---|
| IL.18 | 33 |
| IFNg | 22 |
| IL.2 | 15 |
| IL.7 | 15 |
| KC.like | 13 |
| IL.10 | 10 |
| IL.8 | 3 |
| IL.15 | 2 |
| IL.6 | 2 |
| GM.CSF | 0 |
| IP.10 | 0 |
| MCP.1 | 0 |
| TNFa | 0 |
10.10.2 Remove cytokines with low detection rates
Filtering out: IL-10, IL-6, IL-15, IP-10, GM-CSF, IL-8, MCP-1, TNFa
10.10.3 PCA
10.10.3.1 PCA & data
cyto_pca <-
prcomp(data[,2:8], scale = TRUE)
cyto_pca_data <-
data.frame(
x = cyto_pca$x[,1],
y = cyto_pca$x[,2],
pcbbin = factor(data$pcbbin),
analysis.year = factor(data$analysis.year.x),
iav = data$iav,
iavser = data$iavser,
location = data$location,
year = factor(data$year),
sex = data$sex,
molt.stage = data$molt.stage
)
cyto_pca_labelx <-
paste("PC1 (",
round(abs(summary(cyto_pca)$importance[2,1]*100), digits=0),
"%)",
sep="")
cyto_pca_labely <-
paste("PC2 (",
round(abs(summary(cyto_pca)$importance[2,2]*100), digits=0),
"%)",
sep="")10.10.3.2 PCBS
ggplot(cyto_pca_data,
aes(x=x,
y=y,
col=pcbbin)) +
geom_point() +
stat_ellipse() +
labs(x=cyto_pca_labelx,
y = cyto_pca_labely) +
scale_color_manual("PCBs",
values=c("black", "#008ba2")) +
theme_bw() +
theme(panel.grid=element_blank())
10.10.3.3 IAV
ggplot(cyto_pca_data,
aes(x=x,
y=y,
col=iav)) +
geom_point() +
stat_ellipse() +
labs(x=cyto_pca_labelx,
y = cyto_pca_labely) +
scale_color_manual("IAV",
values=c("#9d9d9d", "#008ba2")) +
theme_bw() +
theme(panel.grid=element_blank())
10.10.3.4 IAV Serology
ggplot(cyto_pca_data,
aes(x=x,
y=y,
col=iavser)) +
geom_point() +
stat_ellipse() +
labs(x=cyto_pca_labelx,
y = cyto_pca_labely) +
scale_color_manual("IAV Serology",
values=c("#9d9d9d", "#008ba2")) +
theme_bw() +
theme(panel.grid=element_blank())
10.10.3.5 Location
ggplot(cyto_pca_data,
aes(x=x,
y=y,
col=location)) +
geom_point() +
stat_ellipse() +
labs(x=cyto_pca_labelx,
y = cyto_pca_labely) +
scale_color_manual("Location",
values=c("#9d9d9d", "#008ba2", "black")) +
theme_bw() +
theme(panel.grid=element_blank())
10.10.3.6 Year
ggplot(cyto_pca_data,
aes(x=x,
y=y,
col=year)) +
geom_point() +
stat_ellipse() +
labs(x=cyto_pca_labelx,
y = cyto_pca_labely) +
scale_color_manual("Year",
values=c("#9d9d9d", "#008ba2", "orchid1", "darkseagreen")) +
theme_bw() +
theme(panel.grid=element_blank())## Too few points to calculate an ellipse## Warning: Removed 1 row containing missing values or
## values outside the scale range (`geom_path()`).
10.10.3.7 Sex
Maybe some separation here?
ggplot(cyto_pca_data,
aes(x=x,
y=y,
col=sex)) +
geom_point() +
stat_ellipse() +
labs(x=cyto_pca_labelx,
y = cyto_pca_labely) +
scale_color_manual("Sex",
values=c("#9d9d9d", "#008ba2")) +
theme_bw() +
theme(panel.grid=element_blank())
10.10.3.8 Molt Stage
ggplot(cyto_pca_data,
aes(x=x,
y=y,
col=molt.stage)) +
geom_point() +
stat_ellipse() +
labs(x=cyto_pca_labelx,
y = cyto_pca_labely) +
scale_color_manual("Molt Stage",
values=c("#9d9d9d", "#008ba2", "black")) +
theme_bw() +
theme(panel.grid=element_blank())## Too few points to calculate an ellipse## Warning: Removed 1 row containing missing values or
## values outside the scale range (`geom_path()`).
10.11 PERMANOVA
10.11.1 Cytokine Presence/Absence
10.11.1.1 Create similarity matrix (Sorensen)
Dissimilarity = 1 - Sorensen
# Dissimilarity summary stats
mean(1-databin_dist); min(1-databin_dist); max(1-databin_dist); median(1-databin_dist)## [1] 0.3914616## [1] 0## [1] 0.75## [1] 0.333333310.11.1.2 Run PERMANOVA
Not significant
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = (1 - databin_dist) ~ pcbbin, data = databin, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.11077 0.03635 1.1695 0.333
## Residual 31 2.93634 0.96365
## Total 32 3.04711 1.0000010.11.1.3 Homogeneity of group dispersions
Not significant
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.006434 0.0064341 1.188 0.2841
## Residuals 31 0.167890 0.0054158
10.11.1.4 Visualization
Principal coordinates analysis
# Use 1 - Sorensen's for dissimilarity
databin_pco <-
cmdscale((1-databin_dist), eig=TRUE)
plot(databin_pco$points, type="n",
cex.lab=1.5, cex.axis=1.1, cex.sub=1.1)
ordiellipse(ord = databin_pco,
groups = factor(data$pcbbin),
display = "sites",
col = c("grey", "darkseagreen"),
lwd = 2,
label = TRUE)
plot(databin_pco$points, type="n",
cex.lab=1.5, cex.axis=1.1, cex.sub=1.1)
ordispider(ord = databin_pco,
groups = factor(data$pcbbin),
display = "sites",
col = c("grey", "darkseagreen"),
lwd=2,
label=TRUE)
10.11.2 Cytokine Concentration
10.11.2.1 Run PERMANOVA
Not significant
data_dist <-
data %>%
select(2:8) %>%
vegdist(., method="bray")
adonis2(data_dist ~ pcbbin,
data=data,
permuations=4999)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = data_dist ~ pcbbin, data = data, permuations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.1743 0.03634 1.169 0.317
## Residual 31 4.6209 0.96366
## Total 32 4.7952 1.0000010.11.2.2 Homogeneity of group dispersions
Not significant
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.00554 0.0055358 0.2939 0.5916
## Residuals 31 0.58399 0.0188385
10.11.2.3 Visualization
# Use Bray-Curtis dissimilarity matrix
data_pco <-
cmdscale(data_dist, eig=TRUE)
plot(data_pco$points, type="n",
cex.lab=1.5, cex.axis=1.1, cex.sub=1.1)
ordiellipse(ord = data_pco,
groups = factor(data$pcbbin),
display = "sites",
col = c("grey", "darkseagreen"),
lwd = 2,
label = TRUE)
plot(data_pco$points, type="n",
cex.lab=1.5, cex.axis=1.1, cex.sub=1.1)
ordispider(ord = data_pco,
groups = factor(data$pcbbin),
display = "sites",
col = c("grey", "darkseagreen"),
lwd=2,
label=TRUE)
10.12 PLS-DA
10.12.3 Visualization
10.12.3.3 Plot
pcb_plsda_ggplot <-
ggplot(plsda_ggplot,
aes(x=x,
y=y,
col=pcb)) +
geom_point() +
stat_ellipse() +
labs(x = plsda_labelx,
y = plsda_labely) +
scale_color_manual("PCB Status",
values=c("black", "#9d9d9d")) +
theme_bw() +
theme(panel.grid = element_blank(),
legend.position = "bottom")
pcb_plsda_ggplot
10.12.4 Cytokine Contributions
plsda_model$loadings$X %>%
data.frame() %>%
select(1:2) %>%
arrange(comp1) %>%
kable() %>%
kable_styling("basic")| comp1 | comp2 | |
|---|---|---|
| IFNg | -0.1596069 | -0.3965643 |
| IL.10 | 0.0519775 | -0.3316747 |
| IL.7 | 0.2396897 | 0.8002565 |
| IL.2 | 0.3147829 | 0.0309254 |
| IL.18 | 0.3986756 | -0.1956570 |
| KC.like | 0.5624330 | -0.0760883 |
| IL.8 | 0.5831050 | -0.2174643 |
10.13 CART
10.13.2 Cross Validation
# Creating task and learner
task <- as_task_classif(pcbbin ~ .,
data = data_cart)
task <- task$set_col_roles(cols="pcbbin",
add_to="stratum")
# min pcb group = 15
learner <- lrn("classif.rpart",
predict_type = "prob",
maxdepth = to_tune(2, 5),
minbucket = to_tune(1, 15),
minsplit = to_tune(1, 30)
)
# Define tuning instance - info ~ tuning process
instance <- ti(task = task,
learner = learner,
resampling = rsmp("cv", folds = 10),
measures = msr("classif.ce"),
terminator = trm("none")
)
# Define how to tune the model
tuner <- tnr("grid_search",
batch_size = 10
)
# Trigger the tuning process
#tuner$optimize(instance)
# optimal values:
# maxdepth = 2
# minbucket = 6
# minsplit = 14
# classif.ce = 0.3510.13.3 Run model with optimized parameters
cart_model <- rpart(formula=pcbbin ~ .,
data=data_cart,
method="class",
maxdepth=2,
minbucket=6,
minsplit=14)
# plot optimized tree
rpart.plot(cart_model)
## IFNg IL.18 KC.like IL.10 IL.2 IL.7 IL.8
## 2.6654235 2.6577730 2.2125157 0.9915008 0.8476720 0.7757576 0.330500310.13.4 Plot
10.13.4.1 Variable Importance
varimp <-
data.frame(imp=cart_model$variable.importance) %>%
rownames_to_column() %>%
rename("variable" = rowname) %>%
arrange(imp) %>%
mutate(variable=gsub("\\.","-",variable),
variable = forcats::fct_inorder(variable))
varimpplot <-
ggplot(varimp) +
geom_segment(aes(x = variable,
y = 0,
xend = variable,
yend = imp),
linewidth = 0.5,
alpha = 0.7) +
geom_point(aes(x = variable,
y = imp),
size = 1,
show.legend = F,
col="black") +
labs(x="Cytokine",
y="Variable Importance") +
coord_flip() +
theme_classic() +
theme(text=element_text(size=10, color="black"))10.13.4.3 Proportions
##
## 0 1
## PCB Absent 22.22222 77.77778
## PCB Present 46.66667 53.33333##
## 1
## PCB Absent 100
## PCB Present 100ifngprop <-
ggplot(data=databin_plot) +
geom_mosaic(aes(x=product(IFNg, pcbbin), fill=IFNg)) +
scale_fill_manual(values=c("Present" = "gray60", "Absent" = "lightgray")) +
scale_y_continuous(labels = scales::percent) +
labs(title="IFNg") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position="bottom",
text = element_text(size=10),
plot.title = element_text(hjust = 0.5))
il18prop <-
ggplot(data=databin_plot) +
geom_mosaic(aes(x=product(IL.18, pcbbin), fill=IL.18)) +
scale_fill_manual(values=c("Present" = "gray60", "Absent" = "lightgray")) +
scale_y_continuous(labels = scales::percent) +
labs(title="IL-18") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position="bottom",
text = element_text(size=10),
plot.title = element_text(hjust = 0.5))10.13.4.4 Concentrations
ifng <-
data %>%
select(pcbbin, IFNg) %>%
filter(!IFNg==0) %>%
ggplot(., aes(x=pcbbin, y=log(IFNg))) +
geom_boxplot(color="black", fill="gray60", lwd=0.25,
outliers=FALSE) +
geom_point(size=0.75, position=position_jitter(width=0.2)) +
labs(x=NULL, y=NULL, title="IFNg") +
stat_n_text(size=3) +
theme_classic() +
theme(text = element_text(size=10),
panel.border=element_rect(color="black", fill=NA, linewidth=0.5))
il18 <-
data %>%
select(pcbbin, IL.18) %>%
filter(!IL.18==0) %>%
ggplot(., aes(x=pcbbin, y=log(IL.18))) +
geom_boxplot(color="black", fill="gray60", lwd=0.25,
outliers=FALSE) +
geom_point(size=0.75, position=position_jitter(width=0.2)) +
labs(x=NULL, y=NULL,title="IL-18") +
stat_n_text(size=3) +
theme_classic() +
theme(text = element_text(size=10),
panel.border=element_rect(color="black", fill=NA, linewidth=0.5))
10.14 Supplemental Figures
10.14.1 Cytokine Presence/Absence
ifng_prop <-
databin_plot %>%
ggplot(data=.) +
geom_mosaic(aes(x=product(IFNg, pcbbin), fill=IFNg)) +
scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray"),
name="Cytokine") +
scale_y_continuous(labels = scales::percent) +
labs(title="IFNg") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position="none",
text = element_text(size=10),
plot.title = element_text(hjust = 0.5))
il2_prop <-
databin_plot %>%
ggplot(data=.) +
geom_mosaic(aes(x=product(IL.2, pcbbin), fill=IL.2)) +
scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray")) +
scale_y_continuous(labels = scales::percent) +
labs(title="IL-2") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position="none",
text = element_text(size=10),
plot.title = element_text(hjust = 0.5))
il7_prop <-
databin_plot %>%
ggplot(data=.) +
geom_mosaic(aes(x=product(IL.7, pcbbin), fill=IL.7)) +
scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray")) +
scale_y_continuous(labels = scales::percent) +
labs(title="IL-7") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position="none",
text = element_text(size=10),
plot.title = element_text(hjust = 0.5))
kc_prop <-
databin_plot %>%
ggplot(data=.) +
geom_mosaic(aes(x=product(KC.like, pcbbin), fill=KC.like)) +
scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray")) +
scale_y_continuous(labels = scales::percent) +
labs(title="KC-like") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position="none",
text = element_text(size=10),
plot.title = element_text(hjust = 0.5))
il18_prop <-
databin_plot %>%
ggplot(data=.) +
geom_mosaic(aes(x=product(IL.18, pcbbin), fill=IL.18)) +
scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray")) +
scale_y_continuous(labels = scales::percent) +
labs(title="IL-18") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position="none",
text = element_text(size=10),
plot.title = element_text(hjust = 0.5))
cyto_prop_grid <-
grid.arrange(ifng_prop, il2_prop, il7_prop,
kc_prop, il18_prop,
ncol=3,
left=textGrob("Proportion of Samples", rot=90),
bottom=textGrob("Influenza A Virus Infection Status"))
10.14.2 Cytokine Concentration - Log
ifng <-
data %>%
select(pcbbin, IFNg) %>%
filter(!IFNg==0) %>%
ggplot(., aes(x=pcbbin, y=log(IFNg))) +
geom_boxplot(color="black", fill="gray60", lwd=0.25,
outliers=FALSE) +
geom_point(size=0.75, position=position_jitter(width=0.2)) +
labs(x=NULL, y=NULL, title = "IFNg") +
stat_n_text(size=3) +
theme_classic() +
theme(text = element_text(size=10),
panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
plot.title = element_text(hjust=0.5))
il2 <-
data %>%
select(pcbbin, IL.2) %>%
filter(!IL.2==0) %>%
ggplot(., aes(x=pcbbin, y=log(IL.2))) +
geom_boxplot(color="black", fill="gray60", lwd=0.25,
outliers=FALSE) +
geom_point(size=0.75, position=position_jitter(width=0.2)) +
labs(x=NULL, y=NULL, title = "IL-2") +
stat_n_text(size=3) +
theme_classic() +
theme(text = element_text(size=10),
panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
plot.title = element_text(hjust=0.5))
il7 <-
data %>%
select(pcbbin, IL.7) %>%
filter(!IL.7==0) %>%
ggplot(., aes(x=pcbbin, y=log(IL.7))) +
geom_boxplot(color="black", fill="gray60", lwd=0.25,
outliers=FALSE) +
geom_point(size=0.75, position=position_jitter(width=0.2)) +
labs(x=NULL, y=NULL, title = "IL-7") +
stat_n_text(size=3) +
theme_classic() +
theme(text = element_text(size=10),
panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
plot.title = element_text(hjust=0.5))
kc <-
data %>%
select(pcbbin, KC.like) %>%
filter(!KC.like==0) %>%
ggplot(., aes(x=pcbbin, y=log(KC.like))) +
geom_boxplot(color="black", fill="gray60", lwd=0.25,
outliers=FALSE) +
geom_point(size=0.75, position=position_jitter(width=0.2)) +
labs(x=NULL, y=NULL, title = "KC-like") +
stat_n_text(size=3) +
theme_classic() +
theme(text = element_text(size=10),
panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
plot.title = element_text(hjust=0.5))
il18 <-
data %>%
select(pcbbin, IL.18) %>%
filter(!IL.18==0) %>%
ggplot(., aes(x=pcbbin, y=log(IL.18))) +
geom_boxplot(color="black", fill="gray60", lwd=0.25,
outliers=FALSE) +
geom_point(size=0.75, position=position_jitter(width=0.2)) +
labs(x=NULL, y=NULL, title = "IL-18") +
stat_n_text(size=3) +
theme_classic() +
theme(text = element_text(size=10),
panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
plot.title = element_text(hjust=0.5))
cyto_conc_grid <-
grid.arrange(ifng, il2, il7,
kc, il18,
ncol=3,
left=textGrob("log(Cytokine Concentration (pg/mL))", rot=90),
bottom=textGrob("PCBs"))