Prediction of Additive PFAS Exposure in Mice Example

An extensive amount of data management is obscured from this document, goal is to inform code writing, figure generation, and etc. not to evaluate these actual data.

library(knitr)
library(tidyverse)
library(readxl)
library(rio)
library(broom)
library(GGally)
library(ggridges)
library(lme4)
library(modelr)
library(nlme)
library(ggrepel)
library(broom.mixed)
library(ggrepel)
library(ggpubr)
library(drc)
library(mgcv)
library(purrr)
library(ggpmisc)
library(equatiomatic)
library(marginaleffects)
library(modelsummary)
library(sjPlot)
library(emmeans)
library(magic)
library(lemon)

###! Tread Carefully Here !##
options(dplyr.summarise.inform=F)
### This turns off the tidyverse summarize/summarise messages indicating default grouping patterns in resultant tibbles ##
### If you don't know what this means, check out https://stackoverflow.com/questions/62140483/how-to-interpret-dplyr-message-summarise-regrouping-output-by-x-override
### Implications of this being false are unknown errors and impacts in other activities in this R session!

summary plots

cbcolors <- c(
  "PFOS" = "#999999", 
  "PFOA" = "#E69F00", 
  "PFHxS" = "#56B4E9",
  "PFHxA" = "#999933",
  "PFNA" = "#F0E442",
  "PFBS" = "#0072B2",
  "6:2 FTS" = "#D55E00",
  "8:2 FTS" = "#CC79A7"
)

notsampleddf <- data.frame(
  treatment=c(rep("PFBS",3),rep("PFNA",3)),
  Tissue=c(rep(c("Brain","Kidney","Liver"),2)),
  COMPOUND=c("NA"),
  mean=c(rep(10,3),rep(10000,3)),
  label="NA",
  Sex="M"
)

reposition_legend(
  ggplot(data=bigplotstacker3|>filter(inmixture=="yes", singletoncompound=="yes"), aes(x=Tissue, y=Concentration, fill=Sex)) +
    geom_dotplot(binaxis="y",  stackdir="center")+
    facet_wrap(~fct_relevel(treatment, "PFOS","PFOA","PFHxS","PFNA","8:2 FTS","PFBS","PFHxA","6:2 FTS"), scales="free_y", axes="all") +
  geom_text(data=notsampleddf, aes(y=mean, label=label), color="black", show.legend=F)+
    scale_y_log10(labels=scales::label_comma())+
    theme_classic(base_size = 20) +
    labs(y="Wholebody, liver, brain, or kidney (ng/g); serum (ng/mL)", x=NULL)+
    theme(axis.text=element_text(color="black"), strip.background=element_rect(fill="white",color="white",linewidth=1), strip.text=element_text(face="bold"), legend.text=element_text(size=20), legend.title=element_text(size=20)) +
    guides(fill=guide_legend(override.aes=list(size=15))),
  position='center', panel='panel-3-3')

ggplot(data=bigplotstacker3|>filter(inmixture=="yes", singletoncompound=="no"), aes(x=Tissue, y=Concentration, fill=COMPOUND)) +
  geom_dotplot(binaxis="y",  stackdir="center", dotsize=1)+
  facet_wrap(~treatment, scales="free_y", axes="all") +
  scale_y_log10(labels=scales::label_comma())+
  theme_classic(base_size = 20) +
  labs(y="Wholebody, liver, brain, or kidney (ng/g); serum (ng/mL)", fill="", x=NULL) +
  guides(fill=guide_legend(title="PFAS"))+
  #scale_color_manual(values=cbcolors)+
  scale_fill_manual(values=cbcolors) +
  theme(axis.text=element_text(color="black"), strip.background=element_rect(fill="white",color="white",linewidth=1), strip.text=element_text(face="bold"), legend.text=element_text(size=14), legend.position="bottom") +
  guides(fill=guide_legend(override.aes=list(size=12), nrow=1))

cbcolors2 <- c(
  "PFOS" = "#999999", 
  "PFOA" = "#E69F00", 
  "PFHxS" = "#56B4E9",
  "PFHxA" = "#999933",
  "PFNA" = "#F0E442",
  "PFBS" = "#0072B2",
  "6:2 FTS" = "#D55E00",
  "8:2 FTS" = "#CC79A7",
  
  "PFPeA"="#AA0DFE",
  "PFPeS"="#3283FE",
  "PFTeDA"="#85660D",
  "PFTrDA"="#782AB6",
  "PFUnA"="#565656",
  "PFDS"="#1C8356",
  "PFDoA"="#16FF32",
  "PFDoS"="#F7E1A0",
  "PFHpA"="#E2E2E2",
  "PFHpS"="#1CBE4F",
  "PFMBA"="#C4451C",
  "PFNS"="#DEA0FD",
  "7:3 FTCA"="#FE00FA",
  "4:2 FTS"="#009e73",
  "EtFOSAA"="#325A9B",
  "MeFOSAA"="#FEAF16",
  "N-EtFOSE"="#F8A19F",
  "N-MeFOSE"="#90AD1C",
  "PFBA"="#F6222E",
  "PFDA"="#1CFFCE"
)



bigplotstacker3 <- bigplotstacker3 |>
  mutate(
    logConcentration=log(Concentration),
    logDose=log(dose)
  )


sumplot2 <- ggplot(data=bigplotstacker3, aes(x=dose, y=Concentration, color=COMPOUND)) +
  geom_point(size=3, alpha=0.25) +
  geom_smooth(method="lm", formula=y~x, se=F, linewidth=2) +
  facet_wrap(~Tissue) +
  scale_y_log10(labels=scales::label_comma())+
  scale_x_log10(labels=scales::label_comma())+
  theme_classic(base_size = 22) +
  labs(y="Wholebody, liver, brain, or kidney (ng/g); serum (ng/mL)", caption="", x="Dose (mg/kg-d)") +
  scale_color_manual(values=cbcolors2) +
  theme(strip.background=element_rect(fill="white",color="white",linewidth=1)) #+
  #scale_fill_manual(values=cbcolors2)

reposition_legend(sumplot2, 'center', panel='panel-3-2')

reposition_legend(
ggplot(data=bigplotstacker3, aes(x=dose, y=Concentration, color=COMPOUND)) +
  geom_point(size=3, alpha=0.25) +
  geom_smooth(method="lm", formula=y~x, se=F, linewidth=2) +
  facet_wrap(~Tissue*Sex) +
  scale_y_log10(labels=scales::label_comma())+
  scale_x_log10(labels=scales::label_comma())+
  theme_classic(base_size = 22) +
  labs(y="Wholebody, liver, brain, or kidney (ng/g); serum (ng/mL)", caption=NULL, color=NULL,x="Dose (mg/kg-d)") +
  scale_color_manual(values=cbcolors2) +
  theme(strip.background=element_rect(fill="white",color="white",linewidth=1)),
 'center', panel='panel-4-3')

all tissues all mixture using broom

fullmixedmodelsdosebased <- bigplotstacker3|>filter(!is.na(logDose)) |>
  nest(data=-Tissue) |>
  mutate(
    fit=map(data, ~ lmer(logConcentration ~logDose+Sex+(logDose|COMPOUND), data=.x)),
    tidied=map(fit, tidy),
    glanced=map(fit, glance),
    augmented=map(fit, augment)
  )

## boundary (singular) fit: see help('isSingular')

fullmixedmodelsdosebased |>
  unnest(tidied) 

fullmixedmodelsdosebased |>
  unnest(tidied) |>
  dplyr::select(Tissue, effect, term, estimate, std.error)|>
  group_by(Tissue, effect, term) |>
  summarize(`Mean, SE`=paste(round(estimate,3),", ",round(std.error,3))) |>
  pivot_wider(names_from=c(effect,term), values_from=`Mean, SE`) 

fullmixedmodelsdosebased |>
  unnest(glanced) 

#fullmixedmodelsdosebased|>
#  unnest(augmented)



ggplot(data=fullmixedmodelsdosebased|>unnest(augmented) , aes(x=logConcentration, y=.fitted, color=COMPOUND, shape=Sex))+
  geom_point() +
  geom_abline(slope=1)+
  geom_abline(slope=1, intercept=-1, linetype="dashed")+
  geom_abline(slope=1, intercept=1, linetype="dashed")+
  facet_wrap(~Tissue) +
  theme_classic(base_size=18) +
  labs(y="Predicted ln(Tissue (ng/g) or serum (ng/mL))", x="Observed ln(Tissue (ng/g) or serum (ng/mL))")+
  scale_color_manual(values=cbcolors2)

ggplot(data=fullmixedmodelsdosebased|>unnest(augmented) , aes(x=logConcentration, y=.fitted, color=COMPOUND, shape=Sex))+
  geom_point() +
  geom_abline(slope=1)+
  geom_abline(slope=1, intercept=-1, linetype="dashed")+
  geom_abline(slope=1, intercept=1, linetype="dashed")+
  facet_wrap(~Tissue*Sex) +
  theme_classic(base_size=18) +
  labs(y="Predicted ln(Tissue (ng/g) or serum (ng/mL))", x="Observed ln(Tissue (ng/g) or serum (ng/mL))")+
  scale_color_manual(values=cbcolors2)

all tissues all mixed to get random effect param values

Whole Body

mixallWB <- lmer(logConcentration ~logDose+Sex+(logDose|COMPOUND), data = bigplotstacker3|>filter(Tissue=="Wholebody",!is.na(logConcentration)), REML=T)
summary(mixallWB)

## Linear mixed model fit by REML ['lmerMod']
## Formula: logConcentration ~ logDose + Sex + (logDose | COMPOUND)
##    Data: 
## filter(bigplotstacker3, Tissue == "Wholebody", !is.na(logConcentration))
## 
## REML criterion at convergence: 1493.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.9700 -0.3582  0.0128  0.4269  4.0496 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  COMPOUND (Intercept) 7.99624  2.8278       
##           logDose     0.04017  0.2004   0.88
##  Residual             0.22966  0.4792       
## Number of obs: 1010, groups:  COMPOUND, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  7.67585    0.86576   8.866
## logDose      0.83113    0.06435  12.916
## SexM         0.31552    0.03031  10.409
## 
## Correlation of Fixed Effects:
##         (Intr) logDos
## logDose  0.870       
## SexM    -0.015  0.011

fixef(mixallWB)

## (Intercept)     logDose        SexM 
##   7.6758520   0.8311269   0.3155191

ranef(mixallWB)

## $COMPOUND
##         (Intercept)     logDose
## 6:2 FTS  -0.5043785 -0.19609479
## 8:2 FTS   1.4884508  0.05694876
## PFBS     -3.6437169 -0.30698430
## PFHpS     2.4612902  0.26647197
## PFHxA    -6.0778594 -0.27626678
## PFHxS     1.6868963  0.13368940
## PFNA      1.3178685  0.01177662
## PFNS      1.5408796  0.06594370
## PFOA      1.7498237  0.13916327
## PFOS      2.0978518  0.16504180
## PFPeS    -2.1171061 -0.05968965
## 
## with conditional variances for "COMPOUND"

coef(mixallWB) 

## $COMPOUND
##         (Intercept)   logDose      SexM
## 6:2 FTS    7.171474 0.6350321 0.3155191
## 8:2 FTS    9.164303 0.8880757 0.3155191
## PFBS       4.032135 0.5241426 0.3155191
## PFHpS     10.137142 1.0975989 0.3155191
## PFHxA      1.597993 0.5548602 0.3155191
## PFHxS      9.362748 0.9648163 0.3155191
## PFNA       8.993721 0.8429035 0.3155191
## PFNS       9.216732 0.8970706 0.3155191
## PFOA       9.425676 0.9702902 0.3155191
## PFOS       9.773704 0.9961687 0.3155191
## PFPeS      5.558746 0.7714373 0.3155191
## 
## attr(,"class")
## [1] "coef.mer"

mixallWB |>
  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`) 

extract_eq(mixallWB, use_coefs=F, mean_separate=T,
           intercept="beta",
           swap_var_names=c(
             "logDose"="ln(Dose, mg/kg-d)",
             "logConcentration"="ln(Whole Body, ng/g)",
             "Sex"="Sex"),
           swap_subscript_names=c(
             "M"="Male"
           ))

\[ \begin{aligned} \operatorname{ln(Whole\ Body,\ ng/g)}_{i} &\sim N \left(\mu, \sigma^2 \right) \\ \mu &=\alpha_{j[i]} + \beta_{1j[i]}(\operatorname{ln(Dose,\ mg/kg-d)}) + \beta_{2}(\operatorname{Sex}_{\operatorname{Male}}) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{1j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &\mu_{\alpha_{j}} \\ &\mu_{\beta_{1j}} \end{aligned} \end{array} \right) , \left( \begin{array}{cc} \sigma^2_{\alpha_{j}} & \rho_{\alpha_{j}\beta_{1j}} \\ \rho_{\beta_{1j}\alpha_{j}} & \sigma^2_{\beta_{1j}} \end{array} \right) \right) \text{, for COMPOUND j = 1,} \dots \text{,J} \end{aligned} \]

extract_eq(mixallWB, use_coefs=T, mean_separate=T, index_factors=T,
           intercept="beta",
           swap_var_names=c(
             "logDose"="ln(Dose, mg/kg-d)",
             "logConcentration"="ln(Whole Body, ng/g)",
             "Sex"="Sex"),
           swap_subscript_names=c(
             "M"="Male"
           ))

\[ \begin{aligned} \operatorname{\widehat{ln(Whole\ Body,\ ng/g)}}_{i} &\sim N \left(\mu, \sigma^2 \right) \\ \mu &=7.68_{\alpha_{j[i]}} + 0.83_{\beta_{1j[i]}}(\operatorname{ln(Dose,\ mg/kg-d)}) + 0.32_{\beta_{2}}(\operatorname{Sex}_{\operatorname{Male}}) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{1j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &0 \\ &0 \end{aligned} \end{array} \right) , \left( \begin{array}{cc} 2.83 & 0.88 \\ 0.88 & 0.2 \end{array} \right) \right) \text{, for COMPOUND j = 1,} \dots \text{,J} \end{aligned} \]

tab_model(mixallWB)

	logConcentration
Predictors	Estimates	CI	p
(Intercept)	7.68	5.98 – 9.37	<0.001
logDose	0.83	0.70 – 0.96	<0.001
Sex [M]	0.32	0.26 – 0.38	<0.001
Random Effects
σ2	0.23
τ00 COMPOUND	8.00
τ11 COMPOUND.logDose	0.04
ρ01 COMPOUND	0.88
ICC	0.96
N COMPOUND	11
Observations	1010
Marginal R2 / Conditional R2	0.460 / 0.979

modelsummary(dvnames(mixallWB), 
             estimate="{estimate} [{conf.low}, {conf.high}]", 
             statistic=NULL,
             coef_rename=coef_rename)

	logConcentration
(Intercept)	7.676 [5.977, 9.375]
logDose	0.831 [0.705, 0.957]
SexM	0.316 [0.256, 0.375]
SD (Intercept COMPOUND)	2.828
SD (logDose COMPOUND)	0.200
Cor (Intercept~logDose COMPOUND)	0.881
SD (Observations)	0.479
Num.Obs.	1010
R2 Marg.	0.460
R2 Cond.	0.979
AIC	1507.6
BIC	1542.0
ICC	1.0
RMSE	0.47

plot_predictions(mixallWB, condition=list("logDose","COMPOUND")) +
  theme_classic() +
  labs(y="Estimated logConcentration")

plot_predictions(mixallWB, condition=list("logDose","COMPOUND")) +
  facet_wrap(~COMPOUND)+
  theme_classic() +
  labs(y="Estimated logConcentration")

all tissues no mixture using broom

nonemixedmodelsdosebased <- bigplotstacker3|>filter(!is.na(logDose)) |>
  nest(data=-Tissue) |>
  mutate(
    fit=map(data, ~ lmer(logConcentration ~logDose+Sex+(1|COMPOUND), data=.x)),
    tidied=map(fit, tidy),
    glanced=map(fit, glance),
    augmented=map(fit, augment)
  )

nonemixedmodelsdosebased |>
  unnest(tidied) 

nonemixedmodelsdosebased |>
  unnest(tidied) |>
  dplyr::select(Tissue, effect, term, estimate, std.error)|>
  group_by(Tissue, effect, term) |>
  summarize(`Mean, SE`=paste(round(estimate,3),", ",round(std.error,3))) |>
  pivot_wider(names_from=c(effect,term), values_from=`Mean, SE`) 

nonemixedmodelsdosebased |>
  unnest(glanced) 

#nonemixedmodelsdosebased|>
#  unnest(augmented)


ggplot(data=nonemixedmodelsdosebased|>unnest(augmented) , aes(x=logConcentration, y=.fitted, color=COMPOUND, shape=Sex))+
  geom_point() +
  geom_abline(slope=1)+
  geom_abline(slope=1, intercept=-1, linetype="dashed")+
  geom_abline(slope=1, intercept=1, linetype="dashed")+
  facet_wrap(~Tissue) +
  theme_classic(base_size=18) +
  labs(y="Predicted ln(Tissue (ng/g) or serum (ng/mL))", x="Observed ln(Tissue (ng/g) or serum (ng/mL))")+
  scale_color_manual(values=cbcolors2)

ggplot(data=nonemixedmodelsdosebased|>unnest(augmented) , aes(x=logConcentration, y=.fitted, color=COMPOUND, shape=Sex))+
  geom_point() +
  geom_abline(slope=1)+
  geom_abline(slope=1, intercept=-1, linetype="dashed")+
  geom_abline(slope=1, intercept=1, linetype="dashed")+
  facet_wrap(~Tissue*Sex) +
  theme_classic(base_size=18) +
  labs(y="Predicted ln(Tissue (ng/g) or serum (ng/mL))", x="Observed ln(Tissue (ng/g) or serum (ng/mL))")+
  scale_color_manual(values=cbcolors2)

all tissues no mixed to get random effect param values

Whole Body

nomixallWB <- lmer(logConcentration ~logDose+Sex+(1|COMPOUND), data = bigplotstacker3|>filter(Tissue=="Wholebody",!is.na(logConcentration)), REML=T)
summary(nomixallWB)

## Linear mixed model fit by REML ['lmerMod']
## Formula: logConcentration ~ logDose + Sex + (1 | COMPOUND)
##    Data: 
## filter(bigplotstacker3, Tissue == "Wholebody", !is.na(logConcentration))
## 
## REML criterion at convergence: 1763.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.5698 -0.3434  0.0374  0.4403  5.7816 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  COMPOUND (Intercept) 6.5691   2.5630  
##  Residual             0.3071   0.5541  
## Number of obs: 1010, groups:  COMPOUND, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 7.716252   0.773705   9.973
## logDose     0.908548   0.008793 103.323
## SexM        0.326414   0.035007   9.324
## 
## Correlation of Fixed Effects:
##         (Intr) logDos
## logDose  0.033       
## SexM    -0.022  0.037

fixef(nomixallWB)

## (Intercept)     logDose        SexM 
##   7.7162521   0.9085479   0.3264141

ranef(nomixallWB)

## $COMPOUND
##         (Intercept)
## 6:2 FTS  -0.1624136
## 8:2 FTS   1.4996219
## PFBS     -3.7808041
## PFHpS     1.5836937
## PFHxA    -5.6084134
## PFHxS     1.5314868
## PFNA      1.4245553
## PFNS      1.5908916
## PFOA      1.4490210
## PFOS      1.9662436
## PFPeS    -1.4938828
## 
## with conditional variances for "COMPOUND"

coef(nomixallWB) 

## $COMPOUND
##         (Intercept)   logDose      SexM
## 6:2 FTS    7.553839 0.9085479 0.3264141
## 8:2 FTS    9.215874 0.9085479 0.3264141
## PFBS       3.935448 0.9085479 0.3264141
## PFHpS      9.299946 0.9085479 0.3264141
## PFHxA      2.107839 0.9085479 0.3264141
## PFHxS      9.247739 0.9085479 0.3264141
## PFNA       9.140807 0.9085479 0.3264141
## PFNS       9.307144 0.9085479 0.3264141
## PFOA       9.165273 0.9085479 0.3264141
## PFOS       9.682496 0.9085479 0.3264141
## PFPeS      6.222369 0.9085479 0.3264141
## 
## attr(,"class")
## [1] "coef.mer"

nomixallWB |> glance()

nomixallWB |>
  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`) 

extract_eq(nomixallWB, use_coefs=F, mean_separate=T,
           intercept="beta",
           swap_var_names=c(
             "logDose"="ln(Dose, mg/kg-d)",
             "logConcentration"="ln(Whole Body, ng/g)",
             "Sex"="Sex"),
           swap_subscript_names=c(
             "M"="Male"
           ))

\[ \begin{aligned} \operatorname{ln(Whole\ Body,\ ng/g)}_{i} &\sim N \left(\mu, \sigma^2 \right) \\ \mu &=\alpha_{j[i]} + \beta_{1}(\operatorname{ln(Dose,\ mg/kg-d)}) + \beta_{2}(\operatorname{Sex}_{\operatorname{Male}}) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for COMPOUND j = 1,} \dots \text{,J} \end{aligned} \]

extract_eq(nomixallWB, use_coefs=T, mean_separate=T, index_factors=T,
           intercept="beta",
           swap_var_names=c(
             "logDose"="ln(Dose, mg/kg-d)",
             "logConcentration"="ln(Whole Body, ng/g)",
             "Sex"="Sex"),
           swap_subscript_names=c(
             "M"="Male"
           ))

\[ \begin{aligned} \operatorname{\widehat{ln(Whole\ Body,\ ng/g)}}_{i} &\sim N \left(\mu, \sigma^2 \right) \\ \mu &=7.72_{\alpha_{j[i]}} + 0.91_{\beta_{1}}(\operatorname{ln(Dose,\ mg/kg-d)}) + 0.33_{\beta_{2}}(\operatorname{Sex}_{\operatorname{Male}}) \\ \alpha_{j} &\sim N \left(0, 2.56 \right) \text{, for COMPOUND j = 1,} \dots \text{,J} \end{aligned} \]

tab_model(nomixallWB)

	logConcentration
Predictors	Estimates	CI	p
(Intercept)	7.72	6.20 – 9.23	<0.001
logDose	0.91	0.89 – 0.93	<0.001
Sex [M]	0.33	0.26 – 0.40	<0.001
Random Effects
σ2	0.31
τ00 COMPOUND	6.57
ICC	0.96
N COMPOUND	11
Observations	1010
Marginal R2 / Conditional R2	0.464 / 0.976

modelsummary(dvnames(nomixallWB), 
             estimate="{estimate} [{conf.low}, {conf.high}]", 
             statistic=NULL,
             coef_rename=coef_rename)

	logConcentration
(Intercept)	7.716 [6.198, 9.235]
logDose	0.909 [0.891, 0.926]
SexM	0.326 [0.258, 0.395]
SD (Intercept COMPOUND)	2.563
SD (Observations)	0.554
Num.Obs.	1010
R2 Marg.	0.464
R2 Cond.	0.976
AIC	1773.2
BIC	1797.8
ICC	1.0
RMSE	0.55

plot_predictions(nomixallWB, condition=list("logDose","COMPOUND")) +
  theme_classic() +
  labs(y="Estimated logConcentration")

plot_predictions(nomixallWB, condition=list("logDose","COMPOUND")) +
  facet_wrap(~COMPOUND)+
  theme_classic() +
  labs(y="Estimated logConcentration")

matrix(c(rownames(as.data.frame(coef(nomixallWB)$COMPOUND)),
         rep(" & ", NROW(coef(nomixallWB)$COMPOUND)),
         round(as.data.frame(coef(nomixallWB)$COMPOUND),2)[,1],
         rep(" \\", NROW(coef(nomixallWB)$COMPOUND))), ncol=4)

##       [,1]      [,2]  [,3]   [,4] 
##  [1,] "6:2 FTS" " & " "7.55" " \\"
##  [2,] "8:2 FTS" " & " "9.22" " \\"
##  [3,] "PFBS"    " & " "3.94" " \\"
##  [4,] "PFHpS"   " & " "9.3"  " \\"
##  [5,] "PFHxA"   " & " "2.11" " \\"
##  [6,] "PFHxS"   " & " "9.25" " \\"
##  [7,] "PFNA"    " & " "9.14" " \\"
##  [8,] "PFNS"    " & " "9.31" " \\"
##  [9,] "PFOA"    " & " "9.17" " \\"
## [10,] "PFOS"    " & " "9.68" " \\"
## [11,] "PFPeS"   " & " "6.22" " \\"

bmatrix = function(x, digits=NULL, ...) {
  library(xtable)
  default_args = list(include.colnames=FALSE, only.contents=TRUE,
                      include.rownames=FALSE, hline.after=NULL, comment=FALSE,
                      print.results=FALSE)
  passed_args = list(...)
  calling_args = c(list(x=xtable(x, digits=digits)),
                   c(passed_args,
                     default_args[setdiff(names(default_args), names(passed_args))]))
  cat("\\begin{bmatrix}\n",
      do.call(print.xtable, calling_args),
      "\\end{bmatrix}\n")
}

bmatrix(matrix(c(rownames(as.data.frame(coef(nomixallWB)$COMPOUND)),
         round(as.data.frame(coef(nomixallWB)$COMPOUND),2)[,1]), ncol=2))

## \begin{bmatrix}
##   6:2 FTS & 7.55 \\ 
##   8:2 FTS & 9.22 \\ 
##   PFBS & 3.94 \\ 
##   PFHpS & 9.3 \\ 
##   PFHxA & 2.11 \\ 
##   PFHxS & 9.25 \\ 
##   PFNA & 9.14 \\ 
##   PFNS & 9.31 \\ 
##   PFOA & 9.17 \\ 
##   PFOS & 9.68 \\ 
##   PFPeS & 6.22 \\ 
##    \end{bmatrix}

\[ \widehat{ln(Whole\ Body,\ ng/g)} \sim N (\mu, 0.55) \\ \mu =\begin{bmatrix} 6:2 FTS & 7.55 \\ 8:2 FTS & 9.22 \\ PFBS & 3.94 \\ PFHpS & 9.3 \\ PFHxA & 2.11 \\ PFHxS & 9.25 \\ PFNA & 9.14 \\ PFNS & 9.31 \\ PFOA & 9.17 \\ PFOS & 9.68 \\ PFPeS & 6.22 \\ \end{bmatrix} + 0.91(ln(Dose, mg/kg))+ 0.33(Sex_{M=1;F=0}) \]

Comes from:

$$ 
\widehat{ln(Whole\ Body,\ ng/g)}  \sim N (\mu, `r round(as.data.frame(VarCorr(nomixallWB))[2,5],2)`) \\
\mu =\begin{bmatrix}
  6:2 FTS & 7.55 \\ 
  8:2 FTS & 9.22 \\ 
  PFBS & 3.94 \\ 
  PFHpS & 9.3 \\ 
  PFHxA & 2.11 \\ 
  PFHxS & 9.25 \\ 
  PFNA & 9.14 \\ 
  PFNS & 9.31 \\ 
  PFOA & 9.17 \\ 
  PFOS & 9.68 \\ 
  PFPeS & 6.22 \\ 
   \end{bmatrix} + `r round(fixef(nomixallWB),2)[[2]]`(ln(Dose, mg/kg))+ `r round(fixef(nomixallWB),2)[[3]]`(Sex_{M=1;F=0})
$$

written “in line” in the .Rmd file.