## ----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)

