Evaluating Tissue-to-Tissue Relationships with Mixtures Context

The first major assumption is that you have multiple PFAS measures paired per animal and per tissue.
i.e. 1633 results with dosed PFAS and other co-occuring PFAS.

The second major assumption is that you desire to relate tissue to tissue or dose to tissue while accounting or not accounting for mixture interactions. But, don’t have a “full factorial” study design, so can’t actually do the full work up.

We used a mixed effects model, maximum likelihood methods, so not Bayesian hierarchical approaches, but using this method, the toxicologists will listen. Importantly, this is still a hierarchical model–change the y-intercept based upon the PFAS and then estimate using the mean slope, etc. I also reported the indvidual-derived residuals as a 3rd level to highlight the amount of variability that is observed in the data.

(this is kind of an Andrew axe to grind, but individuals are more variable than the trends when we have laboratory study sample sizes, so folks tend to have likely inflated concerns about field data given the noise they observe in smaller sample sizes)

packages

library(knitr)
library(tidyverse)
library(lme4)
library(equatiomatic)
library(broom.mixed)

stacked data

Data should be “stacked” with sufficient ID columns and then a column for PFAS ID, a column for Tissue #1 ID, and a column for Tissue #2 ID.

Use log base e to clean things up, requires considering parameters as multiplicative per unit changes and some farting around to predict, but data are generally going to be log normal anyways, so I’m not concerned about it.

dataframe <- data.frame(
  animalID=c(1,1,1,1,1,2,2,2,2,2),
  treatment=c(rep("C68Low",10)),
  Sex=factor(c(rep("M",5),rep("F",5)), levels=c("F","M")),
  PFAS=factor(c("PFOS","PFHxS","PFOA","PFHxA","PFPeS","PFOS","PFHxS","PFOA","PFHxA","PFPeS")),
  Dose=c(0.428,0.3350,0.125,0.207,NA,0.428,0.3350,0.125,0.207,NA),
  Wholebody=c(3,2,1,2,1,3.1,2.1,1.1,2.1,1.1)
) # example vectors showing 2 animals.

dataframe$logDose <- log(dataframe$Dose)
dataframe$logWholebody <- log(dataframe$Wholebody)

full mixture model (i.e. normally distributed (by PFAS) slope and y-intercept)

mixturemodel <- lmer(logWholebody~logDose+Sex+(logDose|PFAS), data=dataframe, REML=T)  # this reprex doesn't work because too few data

# use below if you have many NAs in dataset:
# mixturemodel <- lmer(logWholebody~logDose+Sex+(logDose|PFAS), data=dataframe|>filter(!is.na(logWholebody)), REML=T) 

# this reprex works, but is ugly:
# mixturemodel <- lmer(logWholebody~logDose+Sex+(logDose|PFAS), data=rbind(dataframe,dataframe), REML=T)  # singular fit because overly simplistic reprex

summary(mixturemodel) # full model summary output  
# Pay attention to variance and stdev of y-intercept and slope by PFAS against residual

coef(mixturemodel) # actual parameter estimates
ranef(mixturemodel) # PFAS-specific adjustments to the mean model y-intercept and slope
fixef(mixturemodel) # overall mean model y-intercept and slope

In the mouse analysis, the variance and stdev associated with the PFAS-specific slope was substantially smaller than the variance and stdev associated with residuals. The interpretation is then that the effect of a PFAS-specific slope is smaller than the individual-to-individual variability–i.e. you could not detect a mixture interaction in the field, so it likely does not exist.

additive mixture model (i.e. normally distributed (by PFAS) y-intercept only)

nomixturemodel <- lmer(logWholebody~logDose+Sex+(1|PFAS), data=dataframe, REML=T)  # this reprex works because simpler model

# use below if you have many NAs in dataset:
# nomixturemodel <- lmer(logWholebody~logDose+Sex+(1|PFAS), data=dataframe|>filter(!is.na(logWholebody)), REML=T) 

summary(nomixturemodel) # full model summary output  

## Linear mixed model fit by REML ['lmerMod']
## Formula: logWholebody ~ logDose + Sex + (1 | PFAS)
##    Data: dataframe
## 
## REML criterion at convergence: -11.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0398 -0.3319  0.0047  0.3450  0.9970 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  PFAS     (Intercept) 0.0383413 0.19581 
##  Residual             0.0003645 0.01909 
## Number of obs: 8, groups:  PFAS, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  1.73790    0.30725   5.656
## logDose      0.75776    0.20803   3.643
## SexM        -0.05642    0.01350  -4.179
## 
## Correlation of Fixed Effects:
##         (Intr) logDos
## logDose  0.947       
## SexM    -0.022  0.000

# Pay attention to variance and stdev of y-intercept by PFAS against residual

coef(nomixturemodel) # actual parameter estimates

## $PFAS
##       (Intercept)   logDose        SexM
## PFHxA    1.938303 0.7577625 -0.05642008
## PFHxS    1.575233 0.7577625 -0.05642008
## PFOA     1.651996 0.7577625 -0.05642008
## PFOS     1.786050 0.7577625 -0.05642008
## 
## attr(,"class")
## [1] "coef.mer"

ranef(nomixturemodel) # PFAS-specific adjustments to the mean model y-intercept 

## $PFAS
##       (Intercept)
## PFHxA  0.20040747
## PFHxS -0.16266227
## PFOA  -0.08589941
## PFOS   0.04815421
## 
## with conditional variances for "PFAS"

fixef(nomixturemodel) # overall mean model y-intercept 

## (Intercept)     logDose        SexM 
##  1.73789565  0.75776247 -0.05642008

# see the broom and broom.mixed packages for increased efficiency and labeling cleanup.
nomixturemodel |>
  broom.mixed::tidy() |>
  #dplyr::select(effect, term, estimate, std.error)|>
  group_by(effect, term) |>
  summarize(`Mean, SE`=paste(round(estimate,3),", ",round(std.error,3))) |>
  pivot_wider(names_from=c(effect,term), values_from=`Mean, SE`) 

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

#equatiomatic just makes nice html math from the lmer functions.

extract_eq(nomixturemodel, 
           use_coefs=T,
           mean_separate=T,
           swap_subscript_names=c("M"="Male"),
           swap_var_names=c("logWholebody"="ln(Wholebody, ng/g)",
                            "logDose"="ln(Dose, mg/kg-d)")
           )

\[ \begin{aligned} \operatorname{\widehat{ln(Wholebody,\ ng/g)}}_{i} &\sim N \left(\mu, \sigma^2 \right) \\ \mu &=1.74_{\alpha_{j[i]}} - 0.06_{\beta_{1}}(\operatorname{Sex}_{\operatorname{Male}}) \\ \alpha_{j} &\sim N \left(0.76_{\gamma_{1}^{\alpha}}(\operatorname{ln(Dose,\ mg/kg-d)}), 0.2 \right) \text{, for PFAS j = 1,} \dots \text{,J} \end{aligned} \]

The below is from a mix of hand coding and R tricks to get nice math.

\[ \widehat{ln(Whole\ Body,\ ng/g)} \sim N (\mu, 0.02) \\ \mu =\begin{bmatrix} PFHxA & 1.94 \\ PFHxS & 1.58 \\ PFOA & 1.65 \\ PFOS & 1.79 \\ \end{bmatrix} + 0.76(ln(Dose, mg/kg))+ -0.06(Sex_{M=1;F=0}) \]

Comes from:

$$ 
\widehat{ln(Whole\ Body,\ ng/g)}  \sim N (\mu, `r round(as.data.frame(VarCorr(nomixturemodel))[2,5],2)`) \\
\mu =\begin{bmatrix}
  PFHxA & 1.94 \\ 
  PFHxS & 1.58 \\ 
  PFOA & 1.65 \\ 
  PFOS & 1.79 \\ 
   \end{bmatrix} + `r round(fixef(nomixturemodel),2)[[2]]`(ln(Dose, mg/kg))+ `r round(fixef(nomixturemodel),2)[[3]]`(Sex_{M=1;F=0})
$$

written in line in the .Rmd file.