## ----configure_knitr, eval = TRUE---------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = '#>')
options(rmarkdown.html_vignette.check_title = FALSE)

## ----clear_memory, eval = TRUE------------------------------------------------
rm(list=ls()) 

## ----runchunks, eval = TRUE---------------------------------------------------
# Set whether or not the following chunks will be executed (run):
execute.vignette <- FALSE

## ----load_packages, eval = execute.vignette-----------------------------------
# #
# #
# # Setup
# #
# #
# #library(readxl)
# library(ggplot2)
# library(httk)
# library(scales)
# library(gridExtra)
# library(cowplot)

## ----load_code, eval = execute.vignette---------------------------------------
# scientific_10 <- function(x) {
#   out <- gsub("1e", "10^", scientific_format()(x))
#   out <- gsub("\\+","",out)
#   out <- gsub("10\\^01","10",out)
#   out <- parse(text=gsub("10\\^00","1",out))
# }
# 
# RMSE <- function(x)
# {
#   mean(x$residuals^2)^(1/2)
# }

## ----chem_numbers, eval=FALSE-------------------------------------------------
# length(get_cheminfo(model="fetal_pbtk",suppress.messages=TRUE))

## ----load_aylward_data, eval = execute.vignette-------------------------------
# #MFdata <- read_excel("Aylward-MatFet.xlsx")
# MFdata <- httk::aylward2014
# 
# cat(paste("summarized data from over 100 studies covering ",
#   length(unique(MFdata$DTXSID)[!(unique(MFdata$DTXSID)%in%c("","-"))]),
#   " unique chemicals structures\n",sep=""))
# 
# # We don't want to exclude the volatiles and PFAS at this point:
# MFdata.httk <- subset(MFdata,DTXSID %in% get_cheminfo(
#   info="DTXSID",
#   suppress.messages=TRUE))
# MFdata.httk[MFdata.httk$Chemical.Category=="bromodiphenylethers",
#   "Chemical.Category"] <- "Bromodiphenylethers"
# MFdata.httk[MFdata.httk$Chemical.Category=="organochlorine Pesticides",
#   "Chemical.Category"] <- "Organochlorine Pesticides"
#   MFdata.httk[MFdata.httk$Chemical.Category=="polyaromatic hydrocarbons",
#   "Chemical.Category"] <- "Polyaromatic Hydrocarbons"
# 
# colnames(MFdata.httk)[colnames(MFdata.httk) ==
#   "infant.maternal.conc...Central.tendency..calculate.j.k..or.report.paired.result."] <-
#   "MFratio"
# colnames(MFdata.httk)[colnames(MFdata.httk) ==
#   "PREFERRED_NAME"] <-
#   "Chemical"
# colnames(MFdata.httk)[colnames(MFdata.httk) ==
#   "details.on.matrix.comparison...e.g...cord.blood.lipid..maternal.serum.lipid..or.cord.blood.wet.weight..maternal.whole.blood.wet.weight"] <-
#   "Matrix"
# 
# # Format the columns:
# MFdata.httk$MFratio <- as.numeric(MFdata.httk$MFratio)
# MFdata.httk$Chemical <- as.factor(MFdata.httk$Chemical)
# MFdata.httk$Matrix <- as.factor(MFdata.httk$Chemical)
# MFdata.httk$Chemical.Category <- as.factor(MFdata.httk$Chemical.Category)
# 
# colnames(MFdata.httk)[15] <- "infant"
# colnames(MFdata.httk)[16] <- "maternal"
# colnames(MFdata.httk)[17] <- "obs.units"
# 
# MFdata.httk$infant <- suppressWarnings(as.numeric(MFdata.httk$infant))
# MFdata.httk$maternal <- suppressWarnings(as.numeric(MFdata.httk$maternal))
# MFdata.httk$AVERAGE_MASS <- as.numeric(MFdata.httk$AVERAGE_MASS)

## ----process_aylward_data, eval = execute.vignette----------------------------
# convert1.units <- c("ng/ml","ng/mL","ug/L","ug/l","ng/mL serum","ng/g",
#   "ng/g wet wt.","ppb")
# 
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"infant"] <-
#   MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"infant"] / # ng/ml = ug/L
#   MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"AVERAGE_MASS"]  # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"maternal"] <-
#   MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"maternal"] / # ng/ml = ug/L
#   MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"AVERAGE_MASS"]  # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"obs.units"] <- "uM"
# 
# convert2.units <- c("mg/L","ppm")
# 
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"infant"] <-
#   MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"infant"] * 1000 / # mg/L = ug/L
#   MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"AVERAGE_MASS"]  # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"maternal"] <-
#   MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"maternal"]* 1000 / # mg/L = ug/L
#   MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"AVERAGE_MASS"]  # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"obs.units"] <- "uM"

## ----make_aylward_preds, eval = execute.vignette------------------------------
# # Simulate starting from the 13th week of gestation until full term (40 weeks):
# times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))
# 
# # For each chemical with maternal-fetal ratio data and HTTK in vitro data:
# for (this.id in unique(MFdata.httk$DTXSID))
# {
#   print(this.id)
#   p <- parameterize_steadystate(dtxsid=this.id,
#                         suppress.messages=TRUE)
#   # skip chemicals where Fup is 0:
#   if (p$Funbound.plasma>1e-4)
#   {
# 
#     # Load the chemical-specifc paramaters:
#     p <- parameterize_fetal_pbtk(dtxsid=this.id,
#                                                   fetal_fup_adjustment =TRUE,
#                                                   suppress.messages=TRUE)
#     # Skip chemicals where the 95% credible interval for Fup spans more than 0.1 to
#     # 0.9 (that is, Fup is effectively unknown):
#     if (!is.na(p$Funbound.plasma.dist))
#     {
#       if (as.numeric(strsplit(p$Funbound.plasma.dist,",")[[1]][3])>0.9 &
#           as.numeric(strsplit(p$Funbound.plasma.dist,",")[[1]][2])<0.11)
#       {
#         skip <- TRUE
#       } else skip <- FALSE
#     } else skip <- FALSE
# 
#     if (!skip)
#     {
#       # Solve the PBTK equations for the full simulation time assuming 1 total daily
#       # dose of 1 mg/kg/day spread out over three evenly spaces exposures:
#       out <- solve_fetal_pbtk(
#         parameters=p,
#         dose=0,
#         times=times,
#         daily.dose=1,
#         doses.per.day=3,
#         output.units = "uM",
#         suppress.messages=TRUE)
#       # Identify the concentrations from the final (279th) day:
#       last.row <- which(out[,"time"]>279)
#       last.row <- last.row[!duplicated(out[last.row,"time"])]
#       # Average over the final day:
#       MFdata.httk[MFdata.httk$DTXSID==this.id,"Mat.pred"] <- mean(out[last.row,"Cplasma"])
#       MFdata.httk[MFdata.httk$DTXSID==this.id,"Fet.pred"] <- mean(out[last.row,"Cfplasma"])
#       MFdata.httk[MFdata.httk$DTXSID==this.id,"MFratio.pred"] <-
#         mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
# 
#       # Repeat the analysis without the adjustment to fetal Fup:
#       p <- parameterize_fetal_pbtk(dtxsid=this.id,
#                                                     fetal_fup_adjustment =FALSE,
#                                                     suppress.messages = TRUE)
#       out <- solve_fetal_pbtk(
#         parameters=p,
#         dose=0,
#         times=times,
#         daily.dose=1,
#         doses.per.day=3,
#         output.units = "uM",
#         maxsteps=1e7,
#         suppress.messages = TRUE)
# 
#       last.row <- which(out[,"time"]>279) # The whole final day
#       last.row <- last.row[!duplicated(out[last.row,"time"])]
#       MFdata.httk[MFdata.httk$DTXSID==this.id,"Mat.pred.nofup"] <- mean(out[last.row,"Cplasma"])
#       MFdata.httk[MFdata.httk$DTXSID==this.id,"Fet.pred.nofup"] <- mean(out[last.row,"Cfplasma"])
#       MFdata.httk[MFdata.httk$DTXSID==this.id,"MFratio.pred.nofup"] <-
#         mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
#     }
#   }
# }
# 
# # Something is wrong with cotinine:
# MFdata.httk <- subset(MFdata.httk,Chemical!="Cotinine")
# 
# max.chem <- MFdata.httk[which(MFdata.httk$MFratio==max(MFdata.httk$MFratio)),]
# min.chem <- MFdata.httk[which(MFdata.httk$MFratio==min(MFdata.httk$MFratio)),]
# cat(paste("The minimum observed ratio was ",
#   signif(min.chem[,"MFratio"],2),
#   " for ",
#   min.chem[,"Chemical"],
#   " and the maximum was ",
#   signif(max.chem[,"MFratio"],2),
#   " for ",
#   max.chem[,"Chemical"],
#   ".\n",sep=""))
# 
# 

## ----cleanup_repeated_aylward_figure, eval = execute.vignette-----------------
# MFdata.main <- NULL
# MFdata.outliers <- NULL
# for (this.id in unique(MFdata.httk$DTXSID))
# {
#   this.subset <- subset(MFdata.httk,DTXSID==this.id)
#   this.row <- this.subset[1,]
#   this.row$N.obs <- dim(this.subset)[1]
#   this.row$MFratio <- median(this.subset$MFratio)
#   this.row$MFratio.Q25 <- quantile(this.subset$MFratio,0.25)
#   this.row$MFratio.Q75 <- quantile(this.subset$MFratio,0.75)
#   MFdata.main <- rbind(MFdata.main,this.row)
#   this.subset <- subset(this.subset,
#     MFratio<this.row$MFratio.Q25 |
#     MFratio>this.row$MFratio.Q75)
#   MFdata.outliers <- rbind(MFdata.outliers,this.subset)
# }

## ----aylward_figure_no_fup, eval = execute.vignette---------------------------
# FigCa  <- ggplot(data=MFdata.main) +
#   geom_segment(color="grey",aes(
#     x=MFratio.pred.nofup,
#     y=MFratio.Q25,
#     xend=MFratio.pred.nofup,
#     yend=MFratio.Q75))+
#   geom_point(aes(
#     x=MFratio.pred.nofup,
#     y=MFratio,
#     shape=Chemical.Category,
#     color=Chemical.Category),
#     size=4)   +
#   scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+
#   geom_point(data=MFdata.outliers,aes(
#     x=MFratio.pred.nofup,
#     y=MFratio,
#     shape=Chemical.Category,
#     color=Chemical.Category),
#     size=2)   +
#   xlim(0,2) +
#   ylim(0,3) +
# #   geom_text(aes(x=AUC,y=Critical.concentration,label=Compound.abbrev,color=Chemical)) +
# #   scale_y_log10(label=scientific_10) +
# #,limits=c(10^-7,100)) +
# #   scale_x_log10(label=scientific_10) +
# #   ,limits=c(10^-7,100)) +
# #    annotation_logticks() +
#   geom_abline(slope=1, intercept=0) +
# #    geom_abline(slope=1, intercept=1,linetype="dashed") +
# #    geom_abline(slope=1, intercept=-1,linetype="dashed") +
#   ylab(expression(paste(italic("In vivo")," Mat:Fet Plasma Ratio"))) +
#   xlab("Generic PBTK Predicted Ratio") +
# #    scale_color_brewer(palette="Set2") +
#   theme_bw()  +
#   theme(legend.position="bottom")+
#   theme( text  = element_text(size=14))+
#   theme(legend.text=element_text(size=10))+
#   guides(color=guide_legend(title="Class",nrow=3,byrow=TRUE))+
#   guides(shape=guide_legend(title="Class",nrow=3,byrow=TRUE))
# 
# print(FigCa)

## ----aylward_analysis_text, eval = execute.vignette---------------------------
# 
# cat(paste("In Figure 4 we compare predictions made with our high-throughput \
# human gestational PBTK model with experimental observations on a per chemical \
# basis wherever we had both in vitro HTTK data and in vivo observations (",
# length(unique(MFdata.main$DTXSID)),
# " chemicals).\n",sep=""))
# 
# 
# repeats <- subset(MFdata.main,N.obs>1)
# 
# cat(paste("Multiple observations were available for ",
#   dim(repeats)[1],
#   " of the chemicals,\n",sep=""))
# 
# 
# max.chem <- as.data.frame(repeats[which(repeats$MFratio==max(repeats$MFratio)),])
# min.chem <- as.data.frame(repeats[which(repeats$MFratio==min(repeats$MFratio)),])
# 
# cat(paste("However, among the chemicals with repeated observations, the median \
# observations only ranged from ",
#   signif(min.chem[,"MFratio"],2),
#   " for ",
#   min.chem[,"Chemical"],
#   " to ",
#   signif(max.chem[,"MFratio"],2),
#   " for ",
#   max.chem[,"Chemical"],
#   ".\n",sep=""))
# 
# max.chem <- as.data.frame(MFdata.main[which(MFdata.main$MFratio.pred==max(MFdata.main$MFratio.pred,na.rm=TRUE)),])
# min.chem <- as.data.frame(MFdata.main[which(MFdata.main$MFratio.pred==min(MFdata.main$MFratio.pred,na.rm=TRUE)),])
# 
# cat(paste("The predictions for all chemicals ranged from ",
#   signif(min.chem[,"MFratio.pred"],2),
#   " for ",
#   min.chem[,"Chemical"],
#   " to ",
#   signif(max.chem[,"MFratio.pred"],2),
#   " for ",
#   max.chem[,"Chemical"],
#   ".\n",sep=""))
# 
# 
# 
# fit1 <- lm(data=MFdata.main,MFratio~MFratio.pred.nofup)
# summary(fit1)
# RMSE(fit1)
# 
# fit1sub <- lm(data=subset(MFdata.main,
#   !(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons"))),
#   MFratio~MFratio.pred.nofup)
# summary(fit1sub)
# RMSE(fit1sub)
# 

## ----aylward_figure_with_fup, eval = execute.vignette-------------------------
# FigCb  <- ggplot(data=MFdata.main) +
#   geom_segment(color="grey",aes(
#     x=MFratio.pred,
#     y=MFratio.Q25,
#     xend=MFratio.pred,
#     yend=MFratio.Q75))+
#   geom_point(aes(
#     x=MFratio.pred,
#     y=MFratio,
#     shape=Chemical.Category,
#     color=Chemical.Category),
#     size=3)   +
#   scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+
#   geom_point(data=MFdata.outliers,aes(
#     x=MFratio.pred,
#     y=MFratio,
#     shape=Chemical.Category,
#     color=Chemical.Category),
#     size=1)   +
#   xlim(0,2) +
#   ylim(0,3) +
# #   geom_text(aes(x=AUC,y=Critical.concentration,label=Compound.abbrev,color=Chemical)) +
# #   scale_y_log10(label=scientific_10) +
# #,limits=c(10^-7,100)) +
# #   scale_x_log10(label=scientific_10) +
# #   ,limits=c(10^-7,100)) +
# #    annotation_logticks() +
#   geom_abline(slope=1, intercept=0) +
# #    geom_abline(slope=1, intercept=1,linetype="dashed") +
# #    geom_abline(slope=1, intercept=-1,linetype="dashed") +
#   ylab(expression(paste(italic("In vivo")," Mat:Fet Plasma Ratio"))) +
#   xlab(expression(paste(italic("In vitro")," Predicted Ratio"))) +
# #    scale_color_brewer(palette="Set2") +
#   theme_bw()  +
#   theme(legend.position="bottom")+
#   theme( text  = element_text(size=14))+
#   theme(legend.text=element_text(size=10))+
#   guides(color=guide_legend(title="Class",nrow=3,byrow=TRUE))+
#   guides(shape=guide_legend(title="Class",nrow=3,byrow=TRUE))
# 
# print(FigCb)
# 

## ----voc_analysis, eval = execute.vignette------------------------------------
# # Mean logHenry's law constant for PAH's:
# mean(subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFdata.main,Chemical.Category=="Polyaromatic Hydrocarbons")$DTXSID)$logHenry)
# 
# nonvols <- subset(chem.physical_and_invitro.data,logHenry < -4.5)$DTXSID
# fluoros <- chem.physical_and_invitro.data$DTXSID[regexpr("fluoro",tolower(chem.physical_and_invitro.data$Compound))!=-1]
# 
# fit2 <- lm(data=MFdata.main,MFratio~MFratio.pred)
# summary(fit2)
# RMSE(fit2)
# 
# fit2sub <- lm(data=subset(MFdata.main,
#   !(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons"))),
#   MFratio~MFratio.pred)
# summary(fit2sub)
# 
# RMSE(fit2sub)
# 
# cat(paste("When volatile and fluorinated chemicals are omitted only ",
#   dim(subset(MFdata.main,
#   !(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons"))))[1],
#   " evaluation chemicals remain, but the R2 is ",
#   signif(summary(fit1sub)$adj.r.squared,2),
#   " and the RMSE is ",
#   signif(RMSE(fit1sub),2),
#   " without the correction. When the fetal plasma fraction unbound correction is used, the predictivity decreases, slightly: R2 is ",
#   signif(summary(fit2sub)$adj.r.squared,2),
#   " and the RMSE is ",
#   signif(RMSE(fit2sub),2),
#   " for the non-volatile, non-fluorinated chemicals.\n",
#   sep=""))
# 
# 
# 
# cat(paste("We compare the RMSE for our predictions to the standard deviation \
# of the observations ",
#   signif(sd(MFdata.main$MFratio)[1],2),
#   " (",
#   signif(sd(subset(MFdata.main,
#   !(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons")))$MFratio),2),
#   " for non PAH or fluorinated compounds).\n",sep=""))
# 
# cat(paste("The average standard deviation for chemicals with repeated observations was ",
#   signif(sd(subset(MFdata.main,N.obs>1)$MFratio)[1],2),
#   " (",
#   signif(sd(subset(MFdata.main,
#   N.obs > 1 &
#   !(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons")))$MFratio),2),
#   " for non PAH or fluorinated compounds).\n",sep=""))
# 
# fit3 <- lm(data=repeats,MFratio~MFratio.pred.nofup)
# summary(fit3)
# RMSE(fit3)
# 
# fit3sub <- lm(data=subset(MFdata.main, N.obs > 1 &
#   !(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons"))),
#   MFratio~MFratio.pred.nofup)
# summary(fit3sub)
# 
# 
# fit4 <- lm(data=subset(MFdata.main,N.obs > 1),MFratio~MFratio.pred)
# summary(fit4)
# RMSE(fit4)
# 
# fit4sub <- lm(data=subset(MFdata.main, N.obs > 1 &
#   !(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons"))),
#   MFratio~MFratio.pred)
# summary(fit4sub)
# 
# cat(paste("Overall, we observed relatively poor correlation (R2 = ",
#   signif(summary(fit1)$adj.r.squared,2),
#   ", RMSE = ",
#   signif(RMSE(fit1),2),
#   ") without our fetal fup correction that was unchanged with that assumption (R2 = ",
#   signif(summary(fit2)$adj.r.squared,2),
#   ", RMSE = ",
#   signif(RMSE(fit2),2),
#   ").\n",sep=""))
# 
# repeats <-subset(MFdata.main,N.obs > 1)
# cat(paste("The RMSE of the predictions for the ",
#   dim(subset(repeats,!(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons"))))[1],
#   " non-PAH and non-fluorinated compounds with repeated observations is ",
#   signif(RMSE(fit4sub),2),
#   " with the fup correction and ",
#   signif(RMSE(fit3sub),2),
#   " without.\n",sep=""))
# 
# cat(paste("These values are close to the standard deviation of the mean but the p-value for the chemicals with repeated observations is ",
#   signif(pf(
#     summary(fit4sub)$fstatistic[1],
#     summary(fit4sub)$fstatistic[2],
#     summary(fit4sub)$fstatistic[3]),2),
# " indicating some value for the predictive model, albeit for only ",
# dim(subset(repeats,!(Chemical.Category %in% c(
#   "Fluorinated compounds",
#   "Polyaromatic Hydrocarbons"))))[1],
#  " chemicals",sep=""))
# 
# 
# 
# Fig4.table <- data.frame()
# Fig4.table["Number of Chemicals","All Fig 4"] <- length(unique(MFdata.httk$DTXSID))
# Fig4.table["Number of Observations","All Fig 4"] <- dim(MFdata.httk)[1]
# Fig4.table["Observed Mean (Min - Max)","All Fig 4"] <- paste(
#   signif(mean(MFdata.httk$MFratio),2)," (",
#   signif(min(MFdata.httk$MFratio),2)," - ",
#   signif(max(MFdata.httk$MFratio),2),")",sep="")
# Fig4.table["Observed Standard Deviation","All Fig 4"] <- signif(sd(MFdata.httk$MFratio),2)
# Fig4.table["Predicted Mean (Min - Max)","All Fig 4"] <-  paste(
#   signif(mean(MFdata.main$MFratio.pred,na.rm=TRUE),2)," (",
#   signif(min(MFdata.main$MFratio.pred,na.rm=TRUE),2)," - ",
#   signif(max(MFdata.main$MFratio.pred,na.rm=TRUE),2),")",sep="")
# Fig4.table["RMSE","All Fig 4"] <- signif(RMSE(fit2),2)
# Fig4.table["R-squared","All Fig 4"] <- signif(summary(fit2)[["adj.r.squared"]],2)
# Fig4.table["p-value","All Fig 4"] <- signif(summary(fit2)[["coefficients"]]["MFratio.pred",4],2)
# 
# MFdata.sub1 <- subset(MFdata.httk,
#   !(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons")))
# 
# Fig4.table["Number of Chemicals","No PFAS/PAH"] <- length(unique(MFdata.sub1$DTXSID))
# Fig4.table["Number of Observations","No PFAS/PAH"] <- dim(MFdata.sub1)[1]
# Fig4.table["Observed Mean (Min - Max)","No PFAS/PAH"] <- paste(
#   signif(mean(MFdata.sub1$MFratio,na.rm=TRUE),2)," (",
#   signif(min(MFdata.sub1$MFratio,na.rm=TRUE),2)," - ",
#   signif(max(MFdata.sub1$MFratio,na.rm=TRUE),2),")",sep="")
# Fig4.table["Observed Standard Deviation","No PFAS/PAH"] <- signif(sd(MFdata.sub1$MFratio),2)
# 
# MFdata.sub2 <- subset(MFdata.main,
#   !(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons")))
# 
# Fig4.table["Predicted Mean (Min - Max)","No PFAS/PAH"] <-  paste(
#   signif(mean(MFdata.sub2$MFratio.pred,na.rm=TRUE),2)," (",
#   signif(min(MFdata.sub2$MFratio.pred,na.rm=TRUE),2)," - ",
#   signif(max(MFdata.sub2$MFratio.pred,na.rm=TRUE),2),")",sep="")
# Fig4.table["RMSE","No PFAS/PAH"] <- signif(RMSE(fit2sub),2)
# Fig4.table["R-squared","No PFAS/PAH"] <- signif(summary(fit2sub)[["adj.r.squared"]],2)
# Fig4.table["p-value","No PFAS/PAH"] <- signif(summary(fit2sub)[["coefficients"]]["MFratio.pred",4],2)
# 
# 
# MFdata.sub3 <- subset(MFdata.main,N.obs > 1 &
#   !(Chemical.Category %in% c(
#     "Fluorinated compounds",
#     "Polyaromatic Hydrocarbons")))
# 
# Fig4.table["Number of Chemicals","Replicates Only"] <- length(unique(MFdata.sub3$DTXSID))
# Fig4.table["Number of Observations","Replicates Only"] <- dim(MFdata.sub3)[1]
# Fig4.table["Observed Mean (Min - Max)","Replicates Only"] <- paste(
#   signif(mean(MFdata.sub3$MFratio),2)," (",
#   signif(min(MFdata.sub3$MFratio),2)," - ",
#   signif(max(MFdata.sub3$MFratio),2),")",sep="")
# Fig4.table["Observed Standard Deviation","Replicates Only"] <- signif(sd(MFdata.sub3$MFratio),2)
# 
# Fig4.table["Predicted Mean (Min - Max)","Replicates Only"] <-  paste(
#   signif(mean(MFdata.sub3$MFratio.pred,na.rm=TRUE),2)," (",
#   signif(min(MFdata.sub3$MFratio.pred,na.rm=TRUE),2)," - ",
#   signif(max(MFdata.sub3$MFratio.pred,na.rm=TRUE),2),")",sep="")
# Fig4.table["RMSE","Replicates Only"] <- signif(RMSE(fit4sub),2)
# Fig4.table["R-squared","Replicates Only"] <- signif(summary(fit4sub)[["adj.r.squared"]],2)
# Fig4.table["p-value","Replicates Only"] <- signif(summary(fit4sub)[["coefficients"]]["MFratio.pred",4],2)
# 
# write.csv(Fig4.table[1:6,],file="Fig4Table.csv")

## ----dallmann2018_data, eval = execute.vignette-------------------------------
# #TKstats <- as.data.frame(read_excel("MaternalFetalAUCData.xlsx"))
# TKstats <- httk::pregnonpregaucs
# 
# 
# ids <- unique(TKstats$DTXSID)
# cat(paste(sum(ids %in% get_cheminfo(
#   model="fetal_pbtk",
#   info="dtxsid",
#   suppress.messages=TRUE)),
#           "chemicals for which the model fetal_pbtk can run that are in the Dallmann data set."))
# 
# 
# TKstats.Cmax <- subset(TKstats,DTXSID!="" & Parameter=="Cmax")
# TKstats <- subset(TKstats,DTXSID!="" & Parameter %in% c("AUCinf","AUClast"))
# 
# ids <- unique(TKstats$DTXSID)
# cat(paste(sum(ids %in% get_cheminfo(
#   model="fetal_pbtk",
#   info="dtxsid",
#   suppress.messages=TRUE)),
#           "chemicals for which the model fetal_pbtk can run that have AUCs in the Dallmann data set."))
# 
# 
# # Assumed body weight (kg) from Kapraun 2019
# BW.nonpreg <- 61.103
# 
# #TKstats$Dose <- TKstats$Dose/BW
# #TKstats$Dose.Units <- "mg/kg"
# colnames(TKstats)[colnames(TKstats)=="Observed Pregnant"] <- "Observed.Pregnant"
# colnames(TKstats)[colnames(TKstats)=="Observed Pregnant5"] <- "Observed.Pregnant5"
# colnames(TKstats)[colnames(TKstats)=="Observed Non-Pregnant"] <- "Observed.Non.Pregnant"
# colnames(TKstats)[colnames(TKstats)=="Observed Non-Pregnant2"] <- "Observed.Non.Pregnant2"
# colnames(TKstats)[colnames(TKstats)=="Dose Route"] <- "Dose.Route"

## ----dallmann2018_preds, eval = execute.vignette------------------------------
# for (this.id in unique(TKstats$DTXSID))
# {
#   if (any(regexpr("ng",TKstats[TKstats$DTXSID==this.id,"Dose Units"])!=-1))
#   {
#   }
#   if (this.id %in% get_cheminfo(
#     info="DTXSID",
#     model="fetal_pbtk",
#     suppress.messages=TRUE))
#   {
#     this.subset <- subset(TKstats,DTXSID==this.id)
#     p <- parameterize_pbtk(
#       dtxsid=this.id,
#       suppress.messages=TRUE)
#     p$hematocrit <- 0.39412 # Kapraun 2019 (unitless)
#     p$Rblood2plasma <- calc_rblood2plasma(
#       parameters=p,
#       suppress.messages=TRUE)
#     p$BW <- BW.nonpreg
#     p$Qcardiacc <- 301.78 / p$BW^(3/4) # Kapraun 2019 (L/h/kg^3/4)
#     for (this.route in unique(this.subset$Dose.Route))
#     {
#       this.route.subset <- subset(this.subset, Dose.Route==this.route)
#       if (this.route.subset[1,"Gestational.Age.Weeks"] > 12)
#       {
#         this.dose <- this.route.subset$Dose
#         # Non-pregnant PBPK:
#         out.nonpreg <- solve_pbtk(
#           parameters=p,
#           times = seq(0, this.route.subset[1,"NonPreg.Duration.Days"],
#                       this.route.subset[1,"NonPreg.Duration.Days"]/100),
#           dose=this.dose/BW.nonpreg, # mg/kg
#           daily.dose=NULL,
#           iv.dose=(this.route=="iv"),
#           output.units="uM",
#           suppress.messages=TRUE)
#         # Pregnant PBPK:
#         BW.preg <- suppressWarnings(calc_maternal_bw(
#           week=this.route.subset[1,"Gestational.Age.Weeks"]))
#         out.preg <- solve_fetal_pbtk(
#           dtxsid=this.id,
#           times = seq(
#             this.route.subset[1,"Gestational.Age.Weeks"]*7,
#             this.route.subset[1,"Gestational.Age.Weeks"]*7 +
#               this.route.subset[1,"Preg.Duration.Days"],
#             this.route.subset[1,"Preg.Duration.Days"]/100),
#           dose=this.dose/BW.preg, # mg/kg
#           iv.dose=(this.route=="iv"),
#           daily.dose=NULL,
#           output.units="uM",
#           suppress.messages=TRUE)
#         if (any(regexpr("AUC",this.subset$Parameter)!=-1))
#         {
#           TKstats[TKstats$DTXSID==this.id &
#                     TKstats$Dose.Route == this.route &
#                     regexpr("AUC",TKstats$Parameter)!=-1,
#                   "Predicted.Non.Pregnant.httk"] <- max(out.nonpreg[,"AUC"])*24 #uMdays->uMh
#           TKstats[TKstats$DTXSID==this.id &
#                     TKstats$Dose.Route == this.route &
#                     regexpr("AUC",TKstats$Parameter)!=-1,
#                   "Predicted.Pregnant.httk"] <- max(out.preg[,"AUC"])*24 # uM days -> uM h
#         }
#         if (any(regexpr("Cmax",this.subset$Parameter)!=-1))
#         {
#           TKstats[TKstats$DTXSID==this.id &
#                     TKstats$Dose.Route == this.route &
#                     regexpr("Cmax",TKstats$Parameter)!=-1,
#                   "Predicted.Non.Pregnant.httk"] <- max(out.nonpreg[,"Cplasma"])
#           TKstats[TKstats$DTXSID==this.id &
#                     TKstats$Dose.Route == this.route &
#                     regexpr("Cmax",TKstats$Parameter)!=-1,
#                   "Predicted.Pregnant.httk"] <- max(out.preg[,"Cfplasma"])
#         }
#       }
#     }
#   }
# }
# 
# TKstats$Ratio.obs <- TKstats$Observed.Pregnant / TKstats$Observed.Non.Pregnant
# TKstats$Ratio.httk <- TKstats$Predicted.Pregnant.httk / TKstats$Predicted.Non.Pregnant.httk
# TKstats <- subset(TKstats, !is.na(TKstats$Ratio.httk))
# 
# write.csv(TKstats,file="Table-Dallmann2018.csv",row.names=FALSE)

## ----dallmann2018_figure, eval = execute.vignette-----------------------------
# FigEa  <- ggplot(data=TKstats) +
#   geom_point(aes(
#     y=Observed.Non.Pregnant2,
#     x=Predicted.Non.Pregnant.httk,
#     shape=Dose.Route,
#     alpha=Drug,
#     fill=Drug),
#     size=5)   +
#   geom_abline(slope=1, intercept=0) +
#   geom_abline(slope=1, intercept=1, linetype=3) +
#   geom_abline(slope=1, intercept=-1, linetype=3) +
#   scale_shape_manual(values=c(22,24))+
#   xlab("httk Predicted (uM*h)") +
#   ylab("Observed") +
#   scale_x_log10(#limits=c(10^-3,10^3),
#     label=scientific_10) +
#   scale_y_log10(#limits=c(10^-3,10^3),
#     label=scientific_10) +
#   ggtitle("Non-Pregnant") +
#   theme_bw()  +
#   theme(legend.position="none")+
#   theme( text  = element_text(size=12)) +
#   annotate("text", x = 0.1, y = 1000, label = "A", fontface =2)
# 
# #print(FigEa)
# cat(paste(
#   "For the Dallman et al. (2018) data we observe an average-fold error (AFE)\n\ of",
#   signif(mean(TKstats$Predicted.Non.Pregnant.httk/TKstats$Observed.Non.Pregnant2),2),
#   "and a RMSE (log10-scale) of",
#   signif((mean((log10(TKstats$Predicted.Non.Pregnant.httk)-log10(TKstats$Observed.Non.Pregnant2))^2))^(1/2),2),
#   "for non-pregnant woman.\n"))
# 
# FigEb  <- ggplot(data=TKstats) +
#   geom_point(aes(
#     y=Observed.Pregnant5,
#     x=Predicted.Pregnant.httk,
#     shape=Dose.Route,
#     alpha=Drug,
#     fill=Drug),
#     size=5)   +
#       geom_abline(slope=1, intercept=0) +
#   geom_abline(slope=1, intercept=1, linetype=3) +
#   geom_abline(slope=1, intercept=-1, linetype=3) +
#   scale_shape_manual(values=c(22,24))+
#   xlab("httk Predicted (uM*h)") +
#   ylab("Observed") +
#   scale_x_log10(#limits=c(10^-5,10^3),
#                 label=scientific_10) +
#   scale_y_log10(#limits=c(10^-5,10^3),
#                 label=scientific_10) +
#   ggtitle("Pregnant")+
#   theme_bw()  +
#   theme(legend.position="none")+
#   theme( text  = element_text(size=12)) +
#   annotate("text", x = 0.1, y = 300, label = "B", fontface =2)
# 
# #print(FigEb)
# cat(paste(
#   "We observe an AFE of",
#   signif(mean(TKstats$Predicted.Pregnant.httk/TKstats$Observed.Pregnant5),2),
#   "and a RMSE (log10-scale) of",
#   signif((mean((log10(TKstats$Predicted.Pregnant.httk)-log10(TKstats$Observed.Pregnant5))^2))^(1/2),2),
#   "for pregnant woman.\n"))
# 
# 
# 
# FigEc  <- ggplot(data=TKstats) +
#   geom_point(aes(
#     x=Predicted.Non.Pregnant.httk,
#     y=Predicted.Pregnant.httk,
#     shape=Dose.Route,
#     alpha=Drug,
#     fill=Drug),
#     size=5)   +
#       geom_abline(slope=1, intercept=0) +
#   geom_abline(slope=1, intercept=1, linetype=3) +
#   geom_abline(slope=1, intercept=-1, linetype=3) +
#   scale_shape_manual(values=c(22,24),name="Route")+
#   ylab("httk Pregnant (uM*h)") +
#   xlab("httk Non-Pregnant (uM*h)") +
#   scale_x_log10(#limits=c(10^-1,10^2),
#                 label=scientific_10) +
#   scale_y_log10(#limits=c(10^-1,10^2),
#                 label=scientific_10) +
#   ggtitle("Model Comparison") +
#   theme_bw()  +
#   theme(legend.position="left")+
#   theme( text  = element_text(size=12))+
#   theme(legend.text=element_text(size=10))+
#   guides(fill=guide_legend(ncol=2,override.aes=list(shape=21)),alpha=guide_legend(ncol=2),shape=guide_legend(ncol=2))
#  print(FigEc)
# FigEc <- get_legend(FigEc)
# 
# 
# FigEd  <- ggplot(data=subset(TKstats,!is.na(Ratio.httk))) +
#   geom_point(aes(
#     y=Ratio.obs,
#     x=Ratio.httk,
#     shape=Dose.Route,
#     alpha=Drug,
#     fill=Drug),
#     size=5)   +
#   scale_shape_manual(values=c(22,24))+
#   xlab("httk Predicted") +
#   ylab("Observed") +
#   scale_x_continuous(limits=c(0.25,3)) +
#   scale_y_continuous(limits=c(0.25,3)) +
#   geom_abline(slope=1, intercept=0) +
#  geom_abline(slope=1, intercept=1, linetype=3) +
#   geom_abline(slope=1, intercept=-1, linetype=3) +
#   ggtitle("Ratio") +
#   theme_bw()  +
#   theme(legend.position="none")+
#   theme( text  = element_text(size=12)) +
#   annotate("text", x = 0.5, y = 2.75, label = "C", fontface =2)
# 
# dev.new()
# grid.arrange(FigEa,FigEb,FigEc,FigEd,nrow=2)
# 
# write.csv(subset(TKstats,
#   !is.na(Ratio.httk)),
#   file="DallmanTable.txt")
# 

## ----curley1969_data, eval = execute.vignette---------------------------------
# Curley <- as.data.frame(read_excel("Curley1969.xlsx"))
# dim(Curley)
# Curley.compounds <- Curley[1,4:13]
# Curley <- Curley[4:47,]
# colnames(Curley)[1] <- "Tissue"
# colnames(Curley)[2] <- "N"
# colnames(Curley)[3] <- "Stat"

## ----curley1969_calcpcs, eval = execute.vignette------------------------------
# Curley.pcs <- NULL
# cord.blood <- subset(Curley, Tissue == "Cord Blood")
# suppressWarnings(
# for (this.tissue in unique(Curley$Tissue))
#   if (this.tissue != "Cord Blood")
#   {
#     this.subset <- subset(Curley, Tissue == this.tissue)
#     for (this.chemical in colnames(Curley)[4:13])
#     {
#       if (!is.na(as.numeric(subset(this.subset,Stat=="Mean")[,this.chemical])) &
#         !is.na(as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])))
#       {
#         this.row <- data.frame(
#           Compound = Curley.compounds[,this.chemical],
#           DTXSID = this.chemical,
#           Tissue = this.tissue,
#           PC = as.numeric(subset(this.subset,Stat=="Mean")[,this.chemical]) /
#             as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])
#           )
#         Curley.pcs <- rbind(Curley.pcs,this.row)
#       } else if (!is.na((as.numeric(subset(this.subset,Stat=="Range")[,this.chemical]))) &
#         !is.na((as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical]))))
#       {
#         this.row <- data.frame(
#           Compound = Curley.compounds[,this.chemical],
#           DTXSID = this.chemical,
#           Tissue = this.tissue,
#           PC = as.numeric(subset(this.subset,Stat=="Range")[,this.chemical]) /
#             as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])
#           )
#         Curley.pcs <- rbind(Curley.pcs,this.row)
#       }
#     }
#   }
# )
# Curley.pcs[Curley.pcs$Tissue=="Lungs","Tissue"] <- "Lung"
# 
# pearce2017 <- read_excel("PC_Data.xlsx")
# pearce2017 <- subset(pearce2017,DTXSID%in%Curley.pcs$DTXSID)
# any(pearce2017$DTXSID%in%Curley.pcs$DTXSID)
# print("None of the Curley chems were in the Pearce 2017 PC predictor evaluation.")

## ----csanady2002_pcs, eval = execute.vignette---------------------------------
# csanadybpa <- read_excel("Csanady2002.xlsx",sheet="Table 2")
# csanadybpa$Compound <- "Bisphenol A"
# csanadybpa$DTXSID <- "DTXSID7020182"
# csanadybpa <- csanadybpa[,c("Compound","DTXSID","Tissue","Mean")]
# colnames(csanadybpa) <- colnames(Curley.pcs)
# 
# csanadydaid <- read_excel("Csanady2002.xlsx",sheet="Table 3",skip=1)
# csanadydaid$Compound <- "Daidzein"
# csanadydaid$DTXSID <- "DTXSID9022310"
# csanadydaid <- csanadydaid[,c("Compound","DTXSID","Tissue","Mean...6")]
# colnames(csanadydaid) <- colnames(Curley.pcs)
# 
# Curley.pcs <- rbind(Curley.pcs,csanadybpa,csanadydaid)
# Curley.pcs <- subset(Curley.pcs,Tissue!="Mammary gland")

## ----weijs2013_loadPCs, eval = execute.vignette-------------------------------
# weijstab3 <- read_excel("Weijs2013.xlsx",sheet="Table3",skip=1)
# weijstab3 <- subset(weijstab3, !is.na(Compound) & !is.na(Tissue))
# weijstab4 <- read_excel("Weijs2013.xlsx",sheet="Table4",skip=1)
# weijstab4 <- subset(weijstab4, !is.na(Compound) & !is.na(Tissue))
# 
# for (this.compound in unique(weijstab3$Compound))
#   for (this.tissue in unique(weijstab3$Tissue))
#   {
#     Curley.pcs[
#       Curley.pcs$DTXSID==this.compound & Curley.pcs$Tissue==this.tissue,
#       "Weijs2013"] <- weijstab3[weijstab3$Compound==this.compound &
#                                   weijstab3$Tissue==this.tissue,"value"]
# 
#   }
# 
# for (this.compound in unique(weijstab4$Compound))
#   for (this.tissue in unique(weijstab4$Tissue))
#   {
#     Curley.pcs[
#       Curley.pcs$DTXSID==this.compound & Curley.pcs$Tissue==this.tissue,
#       "Weijs2013"] <- weijstab4[weijstab4$Compound==this.compound &
#                                   weijstab4$Tissue==this.tissue,"value"]
# 
#   }
# 
# write.csv(Curley.pcs,"PCs-table.csv",row.names=FALSE)
# 

## ----curley1969_makepreds, eval = execute.vignette----------------------------
# Curley.pcs <- httk::fetalpcs
# dtxsid.list <- get_cheminfo(
#   info="DTXSID",
#   model="fetal_pbtk",
#   suppress.messages=TRUE)
# 
# suppressWarnings(load_sipes2017())
# for (this.chemical in unique(Curley.pcs$DTXSID))
#   if (this.chemical %in% dtxsid.list)
#   {
#     this.subset <- subset(Curley.pcs,DTXSID==this.chemical)
#     p <- parameterize_fetal_pbtk(dtxsid=this.chemical,
#       fetal_fup_adjustment = FALSE,
#       suppress.messages=TRUE,
#       minimum.Funbound.plasma = 1e-16)
#     fetal.blood.pH <- 7.26
#     Fup <- p$Fraction_unbound_plasma_fetus
#     fetal_schmitt_parms <- parameterize_schmitt(
#       dtxsid=this.chemical,
#       suppress.messages=TRUE,
#       minimum.Funbound.plasma = 1e-16)
#     fetal_schmitt_parms$plasma.pH <- fetal.blood.pH
#     fetal_schmitt_parms$Funbound.plasma <- Fup
#     Curley.pcs[Curley.pcs$DTXSID==this.chemical,"Fup"] <- Fup
#     # httk gives tissue:fup (unbound chemical in plasma) PC's:
#     fetal_pcs <- predict_partitioning_schmitt(
#       parameters = fetal_schmitt_parms,
#       suppress.messages=TRUE,
#       model="fetal_pbtk",
#       minimum.Funbound.plasma = 1e-16)
#     fetal_pcs.nocal <- predict_partitioning_schmitt(
#       parameters = fetal_schmitt_parms,
#       regression=FALSE,
#       suppress.messages=TRUE,
#       model="fetal_pbtk",
#       minimum.Funbound.plasma = 1e-16)
#     out <- solve_fetal_pbtk(
#       dtxsid = this.chemical,
#       fetal_fup_adjustment =FALSE,
#       suppress.messages=TRUE,
#       minimum.Funbound.plasma = 1e-16)
#     Rb2p <- out[dim(out)[1],"Rfblood2plasma"]
#     Curley.pcs[Curley.pcs$DTXSID==this.chemical,"Rb2p"] <- Rb2p
#     # Convert to tissue:blood PC's:
#     for (this.tissue in this.subset$Tissue)
#       if (tolower(this.tissue) %in%
#         unique(subset(tissue.data,Species=="Human")$Tissue))
#       {
#         Curley.pcs[Curley.pcs$DTXSID==this.chemical &
#           Curley.pcs$Tissue == this.tissue, "HTTK.pred.t2up"] <-
#           fetal_pcs[[paste("K",tolower(this.tissue),"2pu",sep="")]]
#         Curley.pcs[Curley.pcs$DTXSID==this.chemical &
#           Curley.pcs$Tissue == this.tissue, "HTTK.pred.nocal.t2up"] <-
#           fetal_pcs.nocal[[paste("K",tolower(this.tissue),"2pu",sep="")]]
#         Curley.pcs[Curley.pcs$DTXSID==this.chemical &
#           Curley.pcs$Tissue == this.tissue, "HTTK.pred"] <-
#           fetal_pcs[[paste("K",tolower(this.tissue),"2pu",sep="")]]*Fup/Rb2p
#         Curley.pcs[Curley.pcs$DTXSID==this.chemical &
#           Curley.pcs$Tissue == this.tissue, "HTTK.pred.nocal"] <-
#           fetal_pcs.nocal[[paste("K",tolower(this.tissue),"2pu",sep="")]]*Fup/Rb2p
#       } else {
#        print(this.tissue)
#       }
#   } else print(this.chemical)
# reset_httk()
# 
# cat(paste(
#   "For the two placental observations (",
#   signif(subset(Curley.pcs,Compound=="Bisphenol A" & Tissue=="Placenta")[,"PC"],2),
#   " for Bisphenol A and ",
#   signif(subset(Curley.pcs,Compound=="Daidzein" & Tissue=="Placenta")[,"PC"],2),
#   " for Diadzen) the predictions were ",
#   signif(subset(Curley.pcs,Compound=="Bisphenol A" & Tissue=="Placenta")[,"HTTK.pred"],2),
#   " and ",
#   signif(subset(Curley.pcs,Compound=="Daidzein" & Tissue=="Placenta")[,"HTTK.pred"],2),
#   ", respectively.\n",sep=""))
# 
# 

## ----curley1969_figure_compounds, eval = execute.vignette---------------------
# FigFa  <- ggplot(data=subset(Curley.pcs,!is.na(HTTK.pred.nocal))) +
#   geom_point(size=2,aes(
#     y=Weijs2013,
#     x=HTTK.pred,
#     shape=Compound,
#     color=Compound)) +
#   geom_point(size=5,aes(
#     y=PC,
#     x=HTTK.pred,
#     shape=Compound,
#     color=Compound)) +
#   geom_abline(slope=1, intercept=0, size=2) +
#   geom_abline(slope=1, intercept=1, linetype=3, size=2) +
#   geom_abline(slope=1, intercept=-1, linetype=3, size=2) +
#   scale_shape_manual(values=rep(c(15,16,17,18),4))+
#   xlab("Predicted Tissue:Blood Partition Coefificent") +
#   ylab("Observed") +
#   scale_x_log10(label=scientific_10,limits=c(10^-1,10^2)) +
#   scale_y_log10(label=scientific_10,limits=c(10^-1,10^4)) +
#   theme_bw()  +
#   theme(legend.position="bottom")+
#   theme( text  = element_text(size=28))+
#   theme(legend.text=element_text(size=12))+
#    theme(legend.title = element_blank())+
#   guides(shape=guide_legend(nrow=3,byrow=TRUE))
# 
# print(FigFa)
# 
# fitFa <- lm(data=Curley.pcs,log10(PC)~log10(HTTK.pred))
# RMSE(fitFa)
# fitFb <- lm(data=Curley.pcs,log10(PC)~log10(HTTK.pred.nocal))
# RMSE(fitFb)

## ----curley1969_figure_tissues, eval = execute.vignette-----------------------
# FigFb  <- ggplot(data=subset(Curley.pcs,!is.na(HTTK.pred))) +
#   geom_point(size=2,aes(
#     y=Weijs2013,
#     x=HTTK.pred,
#     shape=Tissue,
#     color=Tissue)) +
#   geom_point(size=5, aes(
#     y=PC,
#     x=HTTK.pred,
#     shape=Tissue,
#     color=Tissue))   +
#   geom_abline(slope=1, intercept=0, size=2) +
#   geom_abline(slope=1, intercept=1, linetype=3, size=2) +
#   geom_abline(slope=1, intercept=-1, linetype=3, size=2) +
#   scale_shape_manual(values=rep(c(15,16,17,18),4))+
#   xlab("Predicted Tissue:Blood Partition Coefificent") +
#   ylab("Observed") +
#   scale_x_log10(label=scientific_10,limits=c(10^-1,10^2)) +
#   scale_y_log10(label=scientific_10,limits=c(10^-1,10^4)) +
#   theme_bw()  +
#   theme(legend.position="bottom")+
#   theme( text  = element_text(size=28))+
#   theme(legend.text=element_text(size=20))+
#    theme(legend.title = element_blank())+
#   guides(shape=guide_legend(nrow=3,byrow=TRUE))
# 
# print(FigFb)

## ----read_pksim_pcs, eval = execute.vignette----------------------------------
# pksim.pcs <- as.data.frame(read_excel("PartitionCoefficients.xlsx"))
# dim(pksim.pcs)
# for (this.id in unique(pksim.pcs$DTXSID))
# {
#   httk.PCs <- predict_partitioning_schmitt(dtxsid=this.id,
#                                            suppress.messages = TRUE)
#   p <-
#     parameterize_fetal_pbtk(dtxsid=this.id,
#     suppress.messages = TRUE)
#   httk.PCs[["Kplacenta2pu"]] <- p[["Kplacenta2pu"]]
#   httk.PCs[["Fup"]] <- p[["Funbound.plasma"]]
#   lapply(httk.PCs,as.numeric)
#   this.subset <- subset(pksim.pcs,DTXSID==this.id)
#   for (this.tissue in unique(this.subset$Tissue))
#   {
#     if (any(regexpr(tolower(this.tissue),names(httk.PCs))!=-1))
#     {
#       pksim.pcs[pksim.pcs$DTXSID==this.id & pksim.pcs$Tissue==this.tissue,
#         "HTTK.pred"] <-
#         httk.PCs[regexpr(tolower(this.tissue),names(httk.PCs))!=-1][[1]]*
#         httk.PCs[["Fup"]][[1]]
#     }
#   }
# }
# colnames(pksim.pcs)[colnames(pksim.pcs) ==
#                       "Tissue-to-plasma partition coefficient"] <- "PKSim.pred"
# colnames(pksim.pcs)[colnames(pksim.pcs) ==
#         "Method for calculating tissue-to-plasma partition coefficients"] <-
#   "Method"
# pksim.pcs[,"Method"] <- as.factor(pksim.pcs[,"Method"])
# pksim.pcs[,"Drug"] <- as.factor(pksim.pcs[,"Drug"])
# pksim.pcs[,"Tissue"] <- as.factor(pksim.pcs[,"Tissue"])
# 
# 
# pksim.pcs[,"PKSim.pred"] <- as.numeric(
#   pksim.pcs[,"PKSim.pred"])
# 

## ----compare_pksim_pcs, eval = execute.vignette-------------------------------
# pksim.pcs <- httk::pksim.pcs
# 
# pksim.fit1 <- lm(data=pksim.pcs,
#                 log10(PKSim.pred)~log10(HTTK.pred))
# summary(pksim.fit1)
# 
# pksim.fit2 <- lm(data=pksim.pcs,
#                 log10(PKSim.pred)~log10(HTTK.pred)+
#                   Drug+Tissue+Method)
# summary(pksim.fit2)
# 

## ----wang2018_loaddata, eval = execute.vignette-------------------------------
# #Wangchems <- read_excel("Wang2018.xlsx",sheet=3,skip=2)
# Wangchems <- httk::wang2018
# maternal.list <- Wangchems$CASRN[Wangchems$CASRN %in%
#   get_cheminfo(model="fetal_pbtk",
#   suppress.messages=TRUE)]
# nonvols <- subset(chem.physical_and_invitro.data,logHenry < -4.5)$CAS
# nonfluoros <- chem.physical_and_invitro.data$CAS[
#   regexpr("fluoro",tolower(chem.physical_and_invitro.data$Compound))==-1]
# maternal.list <- maternal.list[maternal.list %in% intersect(nonvols,nonfluoros)]

## ----wang2018_makepreds, eval = execute.vignette------------------------------
# 
# pred.table1 <- subset(get_cheminfo(
#   info=c("Compound","CAS","DTXSID","logP","pka_accept","pka_donor"),
#   model="fetal_pbtk"),
#   CAS %in% maternal.list,
#   suppress.messages=TRUE)
# pred.table1$Compound <- gsub("\"","",pred.table1$Compound)
# 
# for (this.cas in maternal.list)
# {
#   Fup <- subset(get_cheminfo(
#     info="all",
#     suppress.messages=TRUE),
#     CAS==this.cas)$Human.Funbound.plasma
#   # Check if Fup is provided as a distribution:
#   if (regexpr(",",Fup)!=-1)
#   {
#     if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
#       (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
#       as.numeric(strsplit(Fup,",")[[1]][2])<0.11))
#     {
#       skip <- TRUE
#     } else skip <- FALSE
#     Fup.upper <- as.numeric(strsplit(Fup,",")[[1]][3])
#     Fup <- as.numeric(strsplit(Fup,",")[[1]][1])
#     # These results are actually from a Bayesian posterior, but we can approximate that
#     # they are normally distributed about the median to estimate a standard deviation:
#     Fup.sigma <- (Fup.upper - Fup)/1.96
#   # If it's not a distribution use defaults (see Wambaugh et al. 2019)
#   } else if (as.numeric(Fup)== 0)
#   {
#     skip <- TRUE
#   } else {
#     skip <- FALSE
#     Fup <- as.numeric(Fup)
#     Fup.sigma <- Fup*0.4
#     Fup.upper  <- min(1,(1+1.96*0.4)*Fup)
#   }
# 
#   # Only run the simulation if we have the necessary parameters:
#   if (!skip)
#   {
#     print(this.cas)
#     out <- solve_fetal_pbtk(
#       chem.cas=this.cas,
#       dose=0,
#       daily.dose=1,
#       doses.per.day=3,
#       fetal_fup_adjustenment=TRUE,
#       suppress.messages=TRUE)
#     last.row <- which(out[,"time"]>279)
#     last.row <- last.row[!duplicated(out[last.row,"time"])]
#     pred.table1[pred.table1$CAS==this.cas,"fup"] <- signif(Fup,3)
#     pred.table1[pred.table1$CAS==this.cas,"fup.sigma"] <- signif(Fup.sigma,3)
#     pred.table1[pred.table1$CAS==this.cas,"fup.upper"] <- signif(Fup.upper,3)
#     pred.table1[pred.table1$CAS==this.cas,"Ratio.pred"] <-
#       signif(mean(out[last.row,"Cfplasma"])/mean(out[last.row,"Cplasma"]),3)
#   }
# }

## ----brain_uncertainty_prop_table, eval = execute.vignette--------------------
# PC.table <- pred.table1
# colnames(PC.table)[colnames(PC.table)=="Ratio.pred"] <- "R.plasma.FtoM"
# pred.table1$Uncertainty <- "Predicted F:M Plasma Ratio"

## ----Rfetmat_uncertainty, eval = execute.vignette-----------------------------
# pred.table3 <- pred.table1
# pred.table3$Uncertainty <- "Plasma Error (Fig. 4)"
# empirical.error <- RMSE(fit2sub)
# for (this.cas in maternal.list)
# {
# 
#   pred.table3[pred.table3$CAS==this.cas,"Ratio.pred"] <- signif(
#     1/((1/pred.table1[pred.table3$CAS==this.cas,"Ratio.pred"])*
#          (1-1.96*empirical.error)), 3)
# }
# # Update final table for paper:
# PC.table$RMSE <- signif(RMSE(fit2sub),3)
# PC.table$R.plasma.FtoM.upper <- pred.table3$Ratio.pred

## ----Rfetbrain, eval = execute.vignette---------------------------------------
# pred.table4 <- pred.table1
# pred.table4$Uncertainty <- "Fetal Brain Partitioning"
# for (this.cas in maternal.list)
# {
#   print(this.cas)
#   p <- parameterize_fetal_pbtk(chem.cas=this.cas,
#       fetal_fup_adjustment=TRUE,
#       suppress.messages=TRUE)
#   Kbrain2pu <- p$Kfbrain2pu
#   fup <- pred.table4[pred.table4$CAS==this.cas,"fup"]
#   pred.table4[pred.table4$CAS==this.cas,"Ratio.pred"] <- signif(
#      PC.table[PC.table$CAS==this.cas,"R.plasma.FtoM"]*
#     Kbrain2pu * fup, 3)
#   PC.table[PC.table$CAS==this.cas,"Kbrain2pu"] <- Kbrain2pu
#   PC.table[PC.table$CAS==this.cas,"R.brain.FtoM"] <-
#     pred.table4[pred.table4$CAS==this.cas,"Ratio.pred"]
# }

## ----Rfetbrain_uncertainty, eval = execute.vignette---------------------------
# 
# pred.table5 <- pred.table1
# pred.table5$Uncertainty <- "Brain Partitioning Error (Pearce et al., 2017)"
# for (this.cas in maternal.list)
# {
#   p <- parameterize_fetal_pbtk(chem.cas=this.cas,
#       fetal_fup_adjustment=TRUE,
#       suppress.messages=TRUE)
#   Kbrain2pu <- p$Kfbrain2pu
# 
#   fup <- pred.table5[pred.table5$CAS==this.cas,"fup"]
#   sigma.fup <- pred.table5[pred.table5$CAS==this.cas,"fup.sigma"]
#   Rmatfet <- 1/PC.table[PC.table$CAS==this.cas,"R.plasma.FtoM"]
#   Rbrain2matblood <- Kbrain2pu * fup / Rmatfet
# 
# # From Pearce et al. (2017) PC paper:
#   sigma.logKbrain <- 0.647
#   Kbrain2pu.upper <- signif(10^(log10(Kbrain2pu)+1.96*sigma.logKbrain),3)
# 
#   error.Rmatfet <- PC.table[PC.table$CAS==this.cas,"RMSE"]
#   sigma.Rbrain2matblood <- Rbrain2matblood *
#     (log(10)^2*sigma.logKbrain^2 +
#     sigma.fup^2/fup^2 +
#     error.Rmatfet^2/Rmatfet^2)^(1/2)
#   Rbrain2matblood.upper <- Rbrain2matblood *
#     (1 + 1.96*(log(10)^2*sigma.logKbrain^2 +
#     sigma.fup^2/fup^2 +
#     error.Rmatfet^2/Rmatfet^2)^(1/2))
#   pred.table5[pred.table5$CAS==this.cas,"Ratio.pred"] <-
#     signif(Rbrain2matblood.upper,3)
#   PC.table[PC.table$CAS==this.cas,"sigma.logKbrain"] <-
#     signif(sigma.logKbrain, 3)
#   PC.table[PC.table$CAS==this.cas,"Kbrain2pu.upper"] <-
#     signif(Kbrain2pu.upper, 3)
#   PC.table[PC.table$CAS==this.cas,"sigma.Rbrain2matblood"] <-
#     signif(sigma.Rbrain2matblood, 3)
#   PC.table[PC.table$CAS==this.cas,"R.brain.FtoM.upper"] <-
#     pred.table5[pred.table5$CAS==this.cas,"Ratio.pred"]
# }

## ----Wang_Uncertainty_Propagation, eval = execute.vignette--------------------
# 
# pred.levels <- pred.table5$Compound[order(pred.table5$Ratio.pred)]
# 
# pred.table <- rbind(
#   pred.table1,
# #  pred.table2,
#   pred.table3,
#   pred.table4,
#   pred.table5)
# pred.table$Compound <- factor(pred.table$Compound,
#   levels = pred.levels)
# 
# pred.table$Uncertainty <- factor(pred.table$Uncertainty,
#   levels = c(pred.table1[1,"Uncertainty"],
# #    pred.table2[1,"Uncertainty"],
#     pred.table3[1,"Uncertainty"],
#     pred.table4[1,"Uncertainty"],
#     pred.table5[1,"Uncertainty"]))

## ----wang2018_figure, eval = execute.vignette---------------------------------
# #Wang 2018 confirmed 6 chemical identities:
# confirmed.chemicals <- c(
#   "2,4-Di-tert-butylphenol",
#   "2,4-Dinitrophenol",
#   "Pyrocatechol",
#   "2'-Hydroxyacetophenone",
#   "3,5-Di-tert-butylsalicylic acid",
#   "4-Hydroxycoumarin"
#   )
# confirmed.chemicals <- c(
#   "96-76-4",
#   "19715-19-6",
#   "51-28-5",
#   "120-80-9",
#   "118-93-4",
#   "1076-38-6")
# 
# 
# FigG  <- ggplot(data=pred.table) +
#   geom_point(aes(
#     x=Ratio.pred,
#     y=Compound,
#     color = Uncertainty,
#     shape = Uncertainty),
#     size=3)   +
#     scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+
#   scale_x_log10(limits=c(10^-2,10^3),label=scientific_10)+
#   ylab(expression(paste(
#     "Chemicals Found in Maternal Plasma by Wang   ",italic("et al.")," (2018)"))) +
#   xlab("Predicted Ratio to Maternal Plasma") +
#   theme_bw()  +
# #  theme(legend.position="bottom")+
#   theme( text  = element_text(size=14))+
#   theme(legend.text=element_text(size=10))#+
#  # guides(color=guide_legend(nrow=4,byrow=TRUE))+
#   #guides(shape=guide_legend(nrow=4,byrow=TRUE))
#   #+
#  # theme(legend.justification = c(0, 0), legend.position = c(0, 0))
# 
# print(FigG)

## ----wang2018_details, eval = execute.vignette--------------------------------
# 
# for (this.col in 7:14)
#   PC.table[,this.col] <- signif(PC.table[,this.col],3)
# 
# PC.table <- PC.table[order(PC.table$R.brain.FtoM.upper,decreasing=TRUE),]
# 
# for (this.row in 1:dim(PC.table)[1])
# {
#   out <- calc_ionization(
#     pH=7.26,
#     pKa_Donor=PC.table[this.row,"pKa_Donor"],
#     pKa_Accept=PC.table[this.row,"pKa_Accept"])
#   if (out$fraction_neutral>0.9) PC.table[this.row,"Charge_726"] <- "Neutral"
#   else if (out$fraction_positive>0.1) PC.table[this.row,"Charge_726"] <-
#     paste(signif(out$fraction_positive*100,2),"% Positive",sep="")
#   else if (out$fraction_negative>0.1) PC.table[this.row,"Charge_726"] <-
#     paste(signif(out$fraction_negative*100,2),"% Negative",sep="")
#   else if (out$fraction_zwitter>0.1) PC.table[this.row,"Charge_726"] <-
#     paste(signif(out$fraction_zwitter*100,2),"% Zwitterion",sep="")
# }
# 
# PC.table <- PC.table[,c(
#   "Compound",
#   "CAS",
#   "DTXSID",
#   "logP",
#   "Charge_726",
#   "R.plasma.FtoM",
#   "RMSE",
#   "R.plasma.FtoM.upper",
#   "Kbrain2pu",
#   "fup",
#   "R.brain.FtoM",
#   "Kbrain2pu.upper",
#   "R.brain.FtoM.upper")]
# 
# write.csv(PC.table,
#   file="WangTable.txt",
#   row.names=FALSE)

## ----impact_of_model, eval = execute.vignette---------------------------------
# MFbrainfit <- lm(R.brain.FtoM.upper~Kbrain2pu.upper,data=PC.table)
# summary(MFbrainfit)
# 
# cat(paste("As expected, the predicted fetal brain to maternal blood ratio was strongly correlated with the predicted brain partition coefficient (R2 = ",
#   signif(summary(MFbrainfit)$adj.r.square,2),
#   ", p-value ",
#   signif(summary(MFbrainfit)$coefficients[2,4],2),
#   "). However, the PBTK simulation impacted the plasma and tissue concentrations such that ",
#   dim(subset(PC.table, rank(R.brain.FtoM.upper) <  rank(Kbrain2pu.upper)))[1],
#   " chemicals were ranked higher than they would have been using partitioning alone.",
#   sep=""))

## ----wang2018_makepreds_nofup, eval = execute.vignette------------------------
# 
# pred.table1.nofup <- subset(get_cheminfo(
#   info=c("Compound","CAS","DTXSID","logP","pka_accept","pka_donor"),
#   model="fetal_pbtk",
#   suppress.messages=TRUE),
#   CAS %in% maternal.list)
# pred.table1.nofup$Compound <- gsub("\"","",pred.table1.nofup$Compound)
# 
# for (this.cas in maternal.list)
# {
#   Fup <- subset(get_cheminfo(
#     info="all",
#     suppress.messages=TRUE),
#     CAS==this.cas)$Human.Funbound.plasma
#   # Check if Fup is provided as a distribution:
#   if (regexpr(",",Fup)!=-1)
#   {
#     if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
#       (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
#       as.numeric(strsplit(Fup,",")[[1]][2])<0.11))
#     {
#       skip <- TRUE
#     } else skip <- FALSE
#     Fup.upper <- as.numeric(strsplit(Fup,",")[[1]][3])
#     Fup <- as.numeric(strsplit(Fup,",")[[1]][1])
#     # These results are actually from a Bayesian posterior, but we can approximate that
#     # they are normally distributed about the median to estimate a standard deviation:
#     Fup.sigma <- (Fup.upper - Fup)/1.96
#   # If it's not a distribution use defaults (see Wambaugh et al. 2019)
#   } else if (as.numeric(Fup)== 0)
#   {
#     skip <- TRUE
#   } else {
#     skip <- FALSE
#     Fup <- as.numeric(Fup)
#     Fup.sigma <- Fup*0.4
#     Fup.upper  <- min(1,(1+1.96*0.4)*Fup)
#   }
# 
#   # Only run the simulation if we have the necessary parameters:
#   if (!skip)
#   {
#     out <- solve_fetal_pbtk(
#       chem.cas=this.cas,
#       dose=0,
#       daily.dose=1,
#       doses.per.day=3,
#       fetal_fup_adjustenment=FALSE,
#       suppress.messages=TRUE)
#     last.row <- which(out[,"time"]>279)
#     last.row <- last.row[!duplicated(out[last.row,"time"])]
#     pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup"] <- signif(Fup,3)
#     pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup.sigma"] <- signif(Fup.sigma,3)
#     pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup.upper"] <- signif(Fup.upper,3)
#     pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"Ratio.pred"] <-
#       signif(mean(out[last.row,"Cfplasma"])/mean(out[last.row,"Cplasma"]),3)
#   }
# }

## ----brain_uncertainty_prop_table_nofup, eval = execute.vignette--------------
# PC.table.nofup <- pred.table1.nofup
# colnames(PC.table.nofup)[colnames(PC.table.nofup)=="Ratio.pred"] <- "R.plasma.FtoM"
# pred.table1.nofup$Uncertainty <- "Predicted F:M Plasma Ratio"

## ----Rfetmat_uncertainty_nofup, eval = execute.vignette-----------------------
# pred.table3.nofup <- pred.table1.nofup
# pred.table3.nofup$Uncertainty <- "Plasma Error (Fig. 4)"
# empirical.error <- RMSE(fit2sub)
# for (this.cas in maternal.list)
# {
# 
#   pred.table3.nofup[pred.table3.nofup$CAS==this.cas,"Ratio.pred"] <- signif(
#     1/((1/pred.table1.nofup[pred.table3.nofup$CAS==this.cas,"Ratio.pred"])*
#          (1-1.96*empirical.error)), 3)
# }
# # Update final table for paper:
# PC.table.nofup$RMSE <- signif(RMSE(fit2sub),3)
# PC.table.nofup$R.plasma.FtoM.upper <- pred.table3.nofup$Ratio.pred

## ----Rfetbrain_nofup, eval = execute.vignette---------------------------------
# pred.table4.nofup <- pred.table1.nofup
# pred.table4.nofup$Uncertainty <- "Fetal Brain Partitioning"
# for (this.cas in maternal.list)
# {
#   p <- parameterize_fetal_pbtk(chem.cas=this.cas,
#       fetal_fup_adjustment=FALSE,
#       suppress.messages=TRUE)
#   Kbrain2pu <- p$Kfbrain2pu
#   fup <- pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"fup"]
#   pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"Ratio.pred"] <- signif(
#      PC.table.nofup[PC.table.nofup$CAS==this.cas,"R.plasma.FtoM"]*
#     Kbrain2pu * fup, 3)
#   PC.table.nofup[PC.table.nofup$CAS==this.cas,"Kbrain2pu"] <- Kbrain2pu
#   PC.table.nofup[PC.table.nofup$CAS==this.cas,"R.brain.FtoM"] <-
#     pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"Ratio.pred"]
# }

## ----Rfetbrain_uncertainty_nofup, eval = execute.vignette---------------------
# 
# pred.table5.nofup <- pred.table1.nofup
# pred.table5.nofup$Uncertainty <- "Brain Partitioning Error (Pearce et al., 2017)"
# for (this.cas in maternal.list)
# {
#   p <- parameterize_fetal_pbtk(chem.cas=this.cas,
#       fetal_fup_adjustment=FALSE,
#       suppress.messages=TRUE)
#   Kbrain2pu <- p$Kfbrain2pu
# 
#   fup <- pred.table5.nofup[pred.table5.nofup$CAS==this.cas,"fup"]
#   sigma.fup <- pred.table5.nofup[pred.table5.nofup$CAS==this.cas,"fup.sigma"]
#   Rmatfet <- 1/PC.table[PC.table.nofup$CAS==this.cas,"R.plasma.FtoM"]
#   Rbrain2matblood <- Kbrain2pu * fup / Rmatfet
# 
# # From Pearce et al. (2017) PC paper:
#   sigma.logKbrain <- 0.647
#   Kbrain2pu.upper <- signif(10^(log10(Kbrain2pu)+1.96*sigma.logKbrain),3)
# 
#   error.Rmatfet <- PC.table.nofup [PC.table.nofup $CAS==this.cas,"RMSE"]
#   sigma.Rbrain2matblood <- Rbrain2matblood *
#     (log(10)^2*sigma.logKbrain^2 +
#     sigma.fup^2/fup^2 +
#     error.Rmatfet^2/Rmatfet^2)^(1/2)
#   Rbrain2matblood.upper <- Rbrain2matblood *
#     (1 + 1.96*(log(10)^2*sigma.logKbrain^2 +
#     sigma.fup^2/fup^2 +
#     error.Rmatfet^2/Rmatfet^2)^(1/2))
#   pred.table5.nofup [pred.table5.nofup $CAS==this.cas,"Ratio.pred"] <-
#     signif(Rbrain2matblood.upper,3)
#   PC.table.nofup [PC.table.nofup $CAS==this.cas,"sigma.logKbrain"] <-
#     signif(sigma.logKbrain, 3)
#   PC.table.nofup [PC.table.nofup $CAS==this.cas,"Kbrain2pu.upper"] <-
#     signif(Kbrain2pu.upper, 3)
#   PC.table.nofup [PC.table.nofup $CAS==this.cas,"sigma.Rbrain2matblood"] <-
#     signif(sigma.Rbrain2matblood, 3)
#   PC.table.nofup [PC.table.nofup $CAS==this.cas,"R.brain.FtoM.upper"] <-
#     pred.table5.nofup [pred.table5.nofup $CAS==this.cas,"Ratio.pred"]
# }

## ----compare_fup_correction_case_study,eval=FALSE-----------------------------
# case.study.fup.correct.table <- merge(PC.table,PC.table.nofup,by=c("Compound","CAS","DTXSID"))
# 
# MFbrainfit.fup <- lm(R.brain.FtoM.upper.x~R.brain.FtoM.upper.y,
#                      data=case.study.fup.correct.table)
# summary(MFbrainfit.fup)
# 
# cat(paste("The predictions for fetal brain to maternal blood ratio with or without the fetal plasma binding correction (Eqn. 1) were strongly correlated (R2 = ",
#   signif(summary(MFbrainfit.fup)$adj.r.square,2),
#   "). There were ",
#   dim(subset(case.study.fup.correct.table,
#              rank(R.brain.FtoM.upper.x) > rank(R.brain.FtoM.upper.y)))[1],
#   " chemicals that were ranked higher with the correction than without, with an average increase of ",
#   signif(mean(
#     case.study.fup.correct.table$R.brain.FtoM.upper.y /
#       case.study.fup.correct.table$R.brain.FtoM.upper.x),2),
#   " times when the plasma binding change was included.\n",
#   sep=""))
# 

## ----fup_table, eval = execute.vignette---------------------------------------
# fup.table <- NULL
# all.chems <- get_cheminfo(
#   model="fetal_pbtk",
#   info="all",
#   suppress.messages=TRUE)
# # Get rid of median fup 0:
# all.chems <- subset(all.chems,
#   as.numeric(unlist(lapply(strsplit(
#     all.chems$Human.Funbound.plasma,","),function(x) x[[1]])))!=0)
# for (this.chem in all.chems[,"CAS"])
# {
#   temp <- parameterize_fetal_pbtk(chem.cas=this.chem,suppress.messages = TRUE)
#   state <- calc_ionization(
#       pH=7.26,
#       pKa_Donor=temp$pKa_Donor,
#       pKa_Accept=temp$pKa_Accept)
#   if (state$fraction_positive > 0.5) this.charge <- "Positive"
#   else if (state$fraction_negative > 0.5) this.charge <- "Negative"
#   else this.charge <- "Neutral"
#   this.row <- data.frame(DTXSID=all.chems[all.chems[,"CAS"]==this.chem,"DTXSID"],
#     Compound=all.chems[all.chems[,"CAS"]==this.chem,"Compound"],
#     CAS=this.chem,
#     Fup.Mat.Pred = temp$Funbound.plasma,
#     Fup.Neo.Pred = temp$Fraction_unbound_plasma_fetus,
#     Charge = this.charge
#     )
#   fup.table <- rbind(fup.table,this.row)
# }
# 
# fup.table[fup.table$Charge=="Positive","Charge"] <- paste("Positive (n=",
#   sum(fup.table$Charge=="Positive"),
#   ")",sep="")
# fup.table[fup.table$Charge=="Negative","Charge"] <- paste("Negative (n=",
#   sum(fup.table$Charge=="Negative"),
#   ")",sep="")
# fup.table[fup.table$Charge=="Neutral","Charge"] <- paste("Neutral (n=",
#   sum(fup.table$Charge=="Neutral"),
#   ")",sep="")

## ----fup_figure, eval = execute.vignette--------------------------------------
# FigA  <- ggplot(data=fup.table) +
#   geom_point(alpha=0.25, aes(
#     x=Fup.Mat.Pred,
#     y=Fup.Neo.Pred,
#     shape=Charge,
#     color=Charge),
#     size=3)   +
#   geom_abline(slope=1, intercept=0) +
#   ylab(expression(paste("Adjusted Neonate ",f[up]))) +
#   xlab(expression(paste(italic("In vitro")," Measured Adult ",f[up]))) +
#    scale_x_log10(label=scientific_10) +
#    scale_y_log10(label=scientific_10) +
#   theme_bw()  +
#   theme( text  = element_text(size=14))
# 
# print(FigA)

## ----MFratio_allhttk_preds, eval = execute.vignette---------------------------
# times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))
# 
# MFratio.pred <- NULL
# all.chems <- get_cheminfo(
#   model="fetal_pbtk",
#   info=c("Compound","DTXSID","CAS","Funbound.plasma"),
#   suppress.messages=TRUE)
# for (this.cas in all.chems$CAS)
#   if ((this.cas %in% nonvols) &
#     !(this.cas %in% fluoros))
# {
#   this.id <- all.chems[all.chems$CAS==this.cas,"DTXSID"]
#   Fup <- subset(all.chems,DTXSID==this.id)$Human.Funbound.plasma
#   if (regexpr(",",Fup)!=-1)
#   {
#     if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
#       (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
#       as.numeric(strsplit(Fup,",")[[1]][2])<0.1))
#     {
#       skip <- TRUE
#     } else skip <- FALSE
#   } else if (Fup== 0)
#   {
#     skip <- TRUE
#   } else skip <- FALSE
# 
#   if (!skip)
#   {
#     p <- parameterize_fetal_pbtk(dtxsid=this.id,
#     fetal_fup_adjustment =TRUE,
#     suppress.messages=TRUE)
#     out <- solve_fetal_pbtk(
#       parameters=p,
#       dose=0,
#       times=times,
#       daily.dose=1,
#       doses.per.day=3,
#       output.units = "uM",
#       suppress.messages=TRUE)
#     last.row <- which(out[,"time"]>279)
#     last.row <- last.row[!duplicated(out[last.row,"time"])]
#     new.row <- data.frame(
#       Chemical = all.chems[all.chems$DTXSID==this.id,"Compound"],
#       DTXSID = this.id,
#       Mat.pred = mean(out[last.row,"Cplasma"]),
#       Fet.pred = mean(out[last.row,"Cfplasma"]),
#       MFratio.pred = mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
#       )
#     MFratio.pred <- rbind(MFratio.pred,new.row)
#   }
#   }

## ----MFratio_allhttk_figure, eval = execute.vignette--------------------------
# FigD <- ggplot(data=MFratio.pred)+
#    geom_histogram(binwidth = 0.05,fill="Red",aes(MFratio.pred))+
#   xlab("Maternal:Fetal Plasma Concentration Ratio") +
#   ylab("Number of chemicals")+
#     theme_bw()+
#     theme( text  = element_text(size=14))
# 
# 
# print(FigD)

## ----MFratio_allhttk_stats, eval = execute.vignette---------------------------
# max.chem <- MFratio.pred[which(
#   MFratio.pred$MFratio.pred==max(MFratio.pred$MFratio.pred,na.rm=TRUE)),]
# min.chem <- MFratio.pred[which(
#   MFratio.pred$MFratio.pred==min(MFratio.pred$MFratio.pred,na.rm=TRUE)),]
# cat(paste("In Figure X we examine the ratios predicted for the ",
#   dim(MFratio.pred)[1],
#   " appropriate (non-volatile or PFAS) chemicals with measured HTTK data.\n",
#   sep=""))
# 
# 
# cat(paste("We observe a median ratio of ",
#   signif(median(MFratio.pred$MFratio.pred,na.rm=TRUE),3),
#   ", with values ranging from ",
#   signif(min.chem[,"MFratio.pred"],3),
#   " for ",
#   min.chem[,"DTXSID"],
#   " to ",
#   signif(max.chem[,"MFratio.pred"],3),
#   " for ",
#   max.chem[,"DTXSID"],
#   ".\n",sep=""))
# 
# # Check out phys-chem > 1.6, < 1:
# highratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred,MFratio.pred>1.6)$DTXSID)
# 
# cat(paste("There are",
#   dim(highratio)[1],
#   "chemicals with ratios greater than 1.6, indicating a tendency to have higher concentrations"))
# 
# # all highly bound
# highratio$Compound
# suppressWarnings(apply(highratio,2,function(x) mean(as.numeric(x),na.rm=TRUE)))
# 
# 
# lowratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred,MFratio.pred<0.9)$DTXSID)
# # No obvious pattern

## ----MFratio_allhttk_preds_nofup, eval = execute.vignette---------------------
# times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))
# 
# MFratio.pred.nofup <- NULL
# all.chems <- get_cheminfo(
#   model="fetal_pbtk",
#   info=c("Compound","DTXSID","CAS","Funbound.plasma"),
#   suppress.messages=TRUE)
# for (this.cas in all.chems$CAS)
#   if ((this.cas %in% nonvols) &
#     !(this.cas %in% fluoros))
# {
#   this.id <- all.chems[all.chems$CAS==this.cas,"DTXSID"]
#   Fup <- subset(all.chems,DTXSID==this.id)$Human.Funbound.plasma
#   if (regexpr(",",Fup)!=-1)
#   {
#     if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
#       (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
#       as.numeric(strsplit(Fup,",")[[1]][2])<0.1))
#     {
#       skip <- TRUE
#     } else skip <- FALSE
#   } else if (Fup== 0)
#   {
#     skip <- TRUE
#   } else skip <- FALSE
# 
#   if (!skip)
#   {
#     p <- parameterize_fetal_pbtk(dtxsid=this.id,
#     fetal_fup_adjustment =FALSE,
#     suppress.messages=TRUE))
#     out <- suppressWarnings(solve_fetal_pbtk(
#       parameters=p,
#       dose=0,
#       times=times,
#       daily.dose=1,
#       doses.per.day=3,
#       output.units = "uM",
#       suppress.messages=TRUE)
#     last.row <- which(out[,"time"]>279)
#     last.row <- last.row[!duplicated(out[last.row,"time"])]
#     new.row <- data.frame(
#       Chemical = all.chems[all.chems$DTXSID==this.id,"Compound"],
#       DTXSID = this.id,
#       Mat.pred = mean(out[last.row,"Cplasma"]),
#       Fet.pred = mean(out[last.row,"Cfplasma"]),
#       MFratio.pred = mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
#       )
#     MFratio.pred.nofup <- rbind(MFratio.pred.nofup,new.row)
#   }
#   }

## ----MFratio_allhttk_figure_nofup, eval = execute.vignette--------------------
# FigSupD <- ggplot(data=MFratio.pred.nofup)+
#    geom_histogram(binwidth = 0.05,fill="Red",aes(MFratio.pred))+
#   xlab("Maternal:Fetal Plasma Concentration Ratio (No Fup Corr.)") +
#   ylab("Number of chemicals")+
#     theme_bw()+
#     theme( text  = element_text(size=14))
# 
# 
# print(FigSupD)

## ----MFratio_allhttk_stats_nofup, eval = execute.vignette---------------------
# max.chem <- MFratio.pred.nofup[which(
#   MFratio.pred.nofup$MFratio.pred==max(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE)),]
# min.chem <- MFratio.pred.nofup[which(
#   MFratio.pred.nofup$MFratio.pred==min(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE)),]
# cat(paste("In Figure X we examine the ratios predicted for the ",
#   dim(MFratio.pred)[1],
#   " appropriate (non-volatile or PFAS) chemicals with measured HTTK data.\n",
#   sep=""))
# 
# 
# cat(paste("We observe a median ratio of ",
#   signif(median(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE),3),
#   ", with values ranging from ",
#   signif(min.chem[,"MFratio.pred"],3),
#   " for ",
#   min.chem[,"DTXSID"],
#   " to ",
#   signif(max.chem[,"MFratio.pred"],3),
#   " for ",
#   max.chem[,"DTXSID"],
#   ".\n",sep=""))
# 
# # Check out phys-chem > 1.6, < 1:
# highratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred.nofup,MFratio.pred>1.6)$DTXSID)
# 
# cat(paste("There are",
#   dim(highratio)[1],
#   "chemicals with ratios greater than 1.6, indicating a tendency to have higher concentrations"))
# 
# # all highly bound
# highratio$Compound
# suppressWarnings(apply(highratio,2,function(x) mean(as.numeric(x),na.rm=2)))
# 
# 
# lowratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred.nofup,MFratio.pred<0.9)$DTXSID)
# # No obvious pattern

## ----create_distributable_data, eval = execute.vignette-----------------------
# aylward2014 <-MFdata
# pregnonpregaucs <- TKstats
# fetalpcs <- Curley.pcs
# wang2018 <- Wangchems
# 
# save(aylward2014,pregnonpregaucs,fetalpcs,wang2018,pksim.pcs,
#      file="Kapraun2022Vignette.RData",version=2)