Plant PFAS mixture uptake and BCF

Data and package pre-work

library(tidyverse, quietly=T)
library(readxl)
library(broom)
library(knitr)
plantdata <- read_excel("Soil-Plant Stacked Data 8-22-25.xlsx", sheet="Comprehensive Data")
plantdata$PFAS <- factor(plantdata$PFAS, levels=c("PFBS","PFHxS","PFHpS","PFOS","PFNS","PFBA","PFPeA","PFHxA","PFOA","PFHpA","PFNA","PFDA","PFOSA","8:2 FTS"))

Data exploration

Kale concentration over time

ggplot() +
  geom_point(data=plantdata|>filter(Material=="Kale"), aes(x=`Age of plant (d)`, y=`concentration  µg/kg`, color=PFAS)) +
  geom_smooth(data=plantdata|>filter(Material=="Kale"), aes(x=`Age of plant (d)`, y=`concentration  µg/kg`, color=PFAS), method="lm", formula=y~x) +
  facet_wrap(~PFAS)+
  scale_y_log10()

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.

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

## Warning: Removed 59 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot() +
  geom_point(data=plantdata|>filter(Material=="Kale"), aes(x=`Age of plant (d)`, y=`concentration  µg/kg`, color=PFAS)) +
  geom_smooth(data=plantdata|>filter(Material=="Kale"), aes(x=`Age of plant (d)`, y=`concentration  µg/kg`, color=PFAS), method="lm", formula=y~x) +
  facet_wrap(~PFAS, scales="free")

## Warning: Removed 59 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 59 rows containing missing values or values outside the scale range
## (`geom_point()`).

Soil concentration over time

ggplot() +
  geom_point(data=plantdata|>filter(Material=="Soil"), aes(x=`Age of plant (d)`, y=`concentration  µg/kg`, color=PFAS)) +
  geom_smooth(data=plantdata|>filter(Material=="Soil"), aes(x=`Age of plant (d)`, y=`concentration  µg/kg`, color=PFAS), method="lm", formula=y~x) +
  facet_wrap(~PFAS)+
  scale_y_log10()

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.

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

## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot() +
  geom_point(data=plantdata|>filter(Material=="Soil"), aes(x=`Age of plant (d)`, y=`concentration  µg/kg`, color=PFAS)) +
  geom_smooth(data=plantdata|>filter(Material=="Soil"), aes(x=`Age of plant (d)`, y=`concentration  µg/kg`, color=PFAS), method="lm", formula=y~x) +
  facet_wrap(~PFAS, scales="free")

## Warning: Removed 5 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).

Data split

soilonly <- plantdata |>
  filter(Material=="Soil")

kaleonly <- plantdata |>
  filter(Material=="Kale")

combodata <- full_join(
  x=soilonly,
  y=kaleonly,
  by=c("Pairing ID","PFAS")
)

combocleanup <- data.frame(
  PFAS=combodata$PFAS,
  Pairing_ID=combodata$`Pairing ID`,
  time=combodata$`Age of plant (d).x`,
  soilconc=combodata$`concentration  µg/kg.x`,
  kaleconc=combodata$`concentration  µg/kg.y`
)

combocleanup$BCF <- combocleanup$kaleconc/combocleanup$soilconc

combocleanup$logkaleconc <- log(combocleanup$kaleconc) 
combocleanup$logsoilconc <- log(combocleanup$soilconc) 
combocleanup$logBCF <- log(combocleanup$BCF) 

BCF plot exploration

ggplot() +
  geom_point(data=combocleanup, aes(x=time, y=BCF, color=PFAS)) +
  geom_smooth(data=combocleanup, aes(x=time, y=BCF, color=PFAS), method="lm", formula=y~x) +
  facet_wrap(~PFAS)+
  scale_y_log10()

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.

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

## Warning: Removed 65 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot() +
  geom_point(data=combocleanup, aes(x=time, y=BCF, color=PFAS)) +
  geom_smooth(data=combocleanup, aes(x=time, y=BCF, color=PFAS), method="lm", formula=y~x) +
  facet_wrap(~PFAS, scales="free")

## Warning: Removed 66 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 65 rows containing missing values or values outside the scale range
## (`geom_point()`).

combocleanup |>
  group_by(PFAS) |>
  summarize(countBCF=sum(!is.na(BCF)),
            meanBCF=mean(BCF, na.rm=T),
            sdBCF=sd(BCF, na.rm=T),
            upperBCF=quantile(BCF, probs=0.975, names=F, na.rm=T),
            lowerBCF=quantile(BCF, probs=0.025, names=F, na.rm=T),
            semBCF=sdBCF/countBCF
            )

linear regression work

kale only

kalereg <- combocleanup |>
  nest(data = -PFAS) |>
  mutate(
    fit = map(data, ~ lm(kaleconc ~ time, data = .x)),
    tidied = map(fit, tidy)
  ) |>
  unnest(tidied)

kalereg |>
  filter(term=="time") |>
  dplyr::select(c("PFAS","estimate","std.error","p.value")) |>
  kable(digits=3, caption="Slope of [kale] by time")

Slope of [kale] by time

PFAS	estimate	std.error	p.value
PFOS	0.364	0.119	0.004
8:2 FTS	-0.400	1.102	0.719
PFBA	1.875	5.358	0.728
PFOSA	-0.012	0.021	0.582
PFHpS	1.129	0.365	0.004
PFOA	2.376	0.691	0.001
PFHpA	2.407	0.817	0.005
PFBS	2.659	1.560	0.095
PFHxA	3.059	1.250	0.018
PFHxS	2.160	0.665	0.002
PFPeA	2.885	3.809	0.453
PFNA	0.062	0.038	0.135
PFDA	-0.140	0.119	0.262
PFNS	-0.002	0.007	0.797

logkalereg <- combocleanup |>
  filter(!logkaleconc%in%c(NA,"-Inf"))|>
  nest(data = -PFAS) |>
  mutate(
    fit = map(data, ~ lm(logkaleconc ~ time, data = .x)),
    tidied = map(fit, tidy)
  ) |>
  unnest(tidied)

logkalereg |>
  filter(term=="time") |>
  dplyr::select(c("PFAS","estimate","std.error","p.value")) |>
  kable(digits=3, caption="Slope of ln([kale]) by time")

Slope of ln([kale]) by time

PFAS	estimate	std.error	p.value
PFOS	0.060	0.015	0.000
8:2 FTS	0.049	0.024	0.045
PFBA	0.023	0.015	0.131
PFOSA	0.082	0.045	0.130
PFHpS	0.061	0.015	0.000
PFOA	0.066	0.016	0.000
PFHpA	0.045	0.014	0.002
PFBS	0.042	0.016	0.013
PFHxA	0.054	0.017	0.003
PFHxS	0.060	0.015	0.000
PFPeA	0.037	0.017	0.037
PFNA	0.067	0.039	0.117
PFDA	NA	NA	NA
PFNS	-0.004	0.008	0.662

Soil only

soilreg <- combocleanup |>
  nest(data = -PFAS) |>
  mutate(
    fit = map(data, ~ lm(soilconc ~ time, data = .x)),
    tidied = map(fit, tidy)
  ) |>
  unnest(tidied)

soilreg |>
  filter(term=="time") |>
  dplyr::select(c("PFAS","estimate","std.error","p.value")) |>
  kable(digits=3, caption="Slope of [soil] by time")

Slope of [soil] by time

PFAS	estimate	std.error	p.value
PFOS	0.052	0.271	0.847
8:2 FTS	-0.055	0.126	0.665
PFBA	-0.157	0.184	0.400
PFOSA	-0.189	0.232	0.421
PFHpS	-0.171	0.276	0.538
PFOA	-0.214	0.163	0.195
PFHpA	-0.195	0.155	0.215
PFBS	-0.191	0.316	0.549
PFHxA	-0.071	0.200	0.724
PFHxS	-0.057	0.185	0.762
PFPeA	-0.122	0.332	0.716
PFNA	0.001	0.006	0.842
PFDA	-0.001	0.001	0.296
PFNS	0.009	0.006	0.183

logsoilreg <- combocleanup |>
  filter(!logsoilconc%in%c(NA,"-Inf"))|>
  nest(data = -PFAS) |>
  mutate(
    fit = map(data, ~ lm(logsoilconc ~ time, data = .x)),
    tidied = map(fit, tidy)
  ) |>
  unnest(tidied)

logsoilreg |>
  filter(term=="time") |>
  dplyr::select(c("PFAS","estimate","std.error","p.value")) |>
  kable(digits=3, caption="Slope of ln([soil]) by time")

Slope of ln([soil]) by time

PFAS	estimate	std.error	p.value
PFOS	0.002	0.011	0.873
8:2 FTS	-0.010	0.009	0.287
PFBA	-0.030	0.017	0.087
PFOSA	-0.011	0.013	0.376
PFHpS	-0.015	0.014	0.290
PFOA	-0.027	0.015	0.084
PFHpA	-0.031	0.018	0.098
PFBS	-0.021	0.020	0.303
PFHxA	-0.008	0.019	0.677
PFHxS	-0.015	0.018	0.417
PFPeA	-0.018	0.021	0.399
PFNA	0.008	0.014	0.577
PFDA	-0.015	0.007	0.046
PFNS	0.011	0.017	0.542

BCF only

BCFreg <- combocleanup |>
  filter(!BCF%in%c(NA,"Inf","NaN"))|>
  nest(data = -PFAS) |>
  mutate(
    fit = map(data, ~ lm(BCF ~ time, data = .x)),
    tidied = map(fit, tidy)
  ) |>
  unnest(tidied)

BCFreg |>
  filter(term=="time") |>
  dplyr::select(c("PFAS","estimate","std.error","p.value")) |>
  kable(digits=3, caption="Slope of [BCF] by time")

Slope of [BCF] by time

PFAS	estimate	std.error	p.value
PFOS	0.021	0.006	0.001
8:2 FTS	-0.062	0.163	0.705
PFBA	4.268	7.666	0.581
PFOSA	-0.001	0.002	0.695
PFHpS	0.175	0.066	0.011
PFOA	1.233	0.542	0.028
PFHpA	3.110	1.566	0.053
PFBS	3.304	3.894	0.401
PFHxA	1.478	1.616	0.365
PFHxS	1.485	1.061	0.169
PFPeA	0.730	8.180	0.929
PFNA	0.226	0.128	0.108
PFDA	-1.079	0.873	0.242
PFNS	0.061	0.105	0.569

BCFmedianpredictions <- combocleanup |>
  filter(!BCF%in%c(NA,"Inf","NaN"))|>
  nest(data = -PFAS) |>
  mutate(
    fit = map(data, ~ lm(BCF ~ time, data = .x))
  ) |>
  mutate(
    prediction=map2(fit,data,~augment(.x, newdata=.y))
  ) |>
  unnest(prediction)

BCFmedianpredictionssumm <- BCFmedianpredictions |>
  filter(time==61)|>
  group_by(PFAS) |>
  summarize(estimate=mean(.fitted),
            sd=sd(.resid),
            count=sum(!is.na(.fitted)),
            sem=sd/sqrt(count))





logBCFreg <- combocleanup |>
  filter(!logBCF%in%c(NA,"-Inf","Inf","NaN"))|>
  nest(data = -PFAS) |>
  mutate(
    fit = map(data, ~ lm(logBCF ~ time, data = .x)),
    tidied = map(fit, tidy)
  ) |>
  unnest(tidied)


logBCFreg |>
  filter(term=="time") |>
  dplyr::select(c("PFAS","estimate","std.error","p.value")) |>
  kable(digits=3, caption="Slope of ln([BCF]) by time")

Slope of ln([BCF]) by time

PFAS	estimate	std.error	p.value
PFOS	0.057	0.017	0.002
8:2 FTS	0.058	0.024	0.022
PFBA	0.053	0.023	0.030
PFOSA	0.064	0.053	0.282
PFHpS	0.073	0.020	0.001
PFOA	0.093	0.023	0.000
PFHpA	0.076	0.025	0.003
PFBS	0.063	0.025	0.017
PFHxA	0.062	0.024	0.012
PFHxS	0.070	0.024	0.006
PFPeA	0.055	0.026	0.037
PFNA	0.066	0.030	0.052
PFDA	NA	NA	NA
PFNS	-0.016	0.024	0.510

boxplot of BCFs

boxplot(logBCF~PFAS, data=combocleanup, las=2)

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outliers (Inf, -Inf) in boxplot 12 are not drawn

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (-Inf) in boxplot 14 is not drawn

boxplot(BCF~PFAS, data=combocleanup, ylim=c(0,10))

## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (Inf) in boxplot 12 is not drawn

Marginal means for BCF at time while accounting for time

ggplot()+
  geom_point(data=BCFmedianpredictionssumm, aes(x=PFAS,y=estimate, color=PFAS), show.legend=F, size=5) +
  geom_linerange(data=BCFmedianpredictionssumm, aes(x=PFAS, ymax=estimate+1.96*sd, ymin=estimate-1.96*sd, color=PFAS), show.legend=F, size=1.1)+
  scale_y_log10()+
  labs(y="estimate plus CI")+
  theme_classic()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning in transformation$transform(x): NaNs produced

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.

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

ggplot()+
  geom_point(data=BCFmedianpredictionssumm, aes(x=PFAS,y=log10(estimate), color=PFAS), show.legend=F, size=5) +
  geom_linerange(data=BCFmedianpredictionssumm, aes(x=PFAS, ymax=log10(estimate+1.96*sd), ymin=log10(estimate-1.96*sd), color=PFAS), show.legend=F, size=1.1)+
  #scale_y_log10()+
  labs(y="log10(estimate) plus CI")+
  theme_classic()

## Warning in FUN(X[[i]], ...): NaNs produced
## Warning in FUN(X[[i]], ...): Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).