## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE
)


## ----setup--------------------------------------------------------------------
# devtools::load_all()
# library(dplyr)
# library(parallel)

## -----------------------------------------------------------------------------
# #select data for a single chemical and species
# my_data <- subset(cvt,
#                  analyte_dtxsid %in% c("DTXSID8031865",
#                                             "DTXSID0020442")
#                  )

## -----------------------------------------------------------------------------
# my_pk <- pk(data = my_data)

## -----------------------------------------------------------------------------
#   my_pk <- my_pk +
#   settings_preprocess(suppress.messages = FALSE,
#                       keep_data_original = TRUE)

## -----------------------------------------------------------------------------
# #preprocess data
# my_pk <- do_preprocess(my_pk)

## -----------------------------------------------------------------------------
# #get summary data info
# my_pk <- do_data_info(my_pk)

## -----------------------------------------------------------------------------
# #pre-fitting (get model parameter bounds & starting values)
# my_pk <- do_prefit(my_pk)

## -----------------------------------------------------------------------------
# #model fitting
# my_pk <- do_fit(my_pk)

## -----------------------------------------------------------------------------
# system.time(my_pk <- do_fit(my_pk)) #without parallel processing

## -----------------------------------------------------------------------------
# system.time(my_pk <- do_fit(my_pk,
#                 n_cores = parallel::detectCores()-1))

## ----eval = FALSE-------------------------------------------------------------
# my_pk <- pk(data = my_data)
# 
# my_pk <- do_fit(my_pk)

## -----------------------------------------------------------------------------
# coef(my_pk)

## -----------------------------------------------------------------------------
# coef(my_pk) %>% slice(1) %>%  mutate(as.data.frame(as.list(coefs_vector[[1]]))) %>%
#   glimpse()

## -----------------------------------------------------------------------------
# coef(my_pk) %>% slice(5) %>%  mutate(as.data.frame(as.list(coefs_vector[[1]]))) %>%
#   glimpse()

## -----------------------------------------------------------------------------
# my_resids <- residuals(my_pk)
# my_resids %>% glimpse()

## -----------------------------------------------------------------------------
# my_preds <- predict(my_pk)
# my_preds %>% glimpse()

## -----------------------------------------------------------------------------
# predict(my_pk, newdata = data.frame(
#   Chemical = "DTXSID8031865",
#   Species = "rat",
#   Time = seq(0, 5, by = 0.5),
#   Time.Units = "hours",
#   Conc.Units = "mg/kg",
#   Dose = 1,
#   Route = "iv",
#   Media = "plasma",
#   exclude = FALSE)) %>% glimpse()

## -----------------------------------------------------------------------------
# p <- plot(x = my_pk)
# print(p$final_plot)

## -----------------------------------------------------------------------------
# p2 <- plot(my_pk,
#      use_scale_conc = list(dose_norm = TRUE,
#                            log10_trans = FALSE)
#      )
# print(p2$final_plot)

## -----------------------------------------------------------------------------
# logLik(my_pk)

## -----------------------------------------------------------------------------
# AIC(my_pk)

## -----------------------------------------------------------------------------
# BIC(my_pk)

## -----------------------------------------------------------------------------
# rmse(my_pk)

## -----------------------------------------------------------------------------
# rsq(my_pk)

## -----------------------------------------------------------------------------
# AFE(my_pk)

## -----------------------------------------------------------------------------
# AAFE(my_pk)

## -----------------------------------------------------------------------------
# data_proc <- get_data(my_pk)
# data_proc %>% glimpse()

## -----------------------------------------------------------------------------
# my_data_info <- get_data_info(my_pk)
# names(my_data_info)
# 
# my_data_info$data_flags %>% glimpse()

## -----------------------------------------------------------------------------
# my_tkstats <- get_tkstats(my_pk, suppress.messages = TRUE)
# my_tkstats %>% glimpse()

## -----------------------------------------------------------------------------
# wm <- get_winning_model(my_pk)
# wm

## -----------------------------------------------------------------------------
# rsq_df <- rsq(my_pk)
# 
# rsq_win <- wm %>% dplyr::left_join(rsq_df)
# 
# rsq_win

## -----------------------------------------------------------------------------
# plot(my_pk, best_fit = TRUE, use_scale_conc = list(dose_norm = TRUE,
#                                                    log10_trans = FALSE)) %>%
#   pull(final_plot)

## -----------------------------------------------------------------------------
# my_pk <- pk(data = my_data)

## -----------------------------------------------------------------------------
# names(my_pk)

## -----------------------------------------------------------------------------
# head(my_pk$data_original)
# print(my_pk)

## ----eval = FALSE-------------------------------------------------------------
# ggplot(data = my_data,
#        aes(x = time_hr,
#            y = conc,
#            color = dose_level_normalized_corrected,
#            shape = as.factor(document_id)
#        )
# )

## ----eval = FALSE-------------------------------------------------------------
# pk(data = my_data,
#    mapping = ggplot2::aes(
#                   mapping = ggplot2::aes(
#                    Chemical = analyte_dtxsid,
#                    Chemical_Name = analyte_name_original,
#                    DTXSID = analyte_dtxsid,
#                    CASRN = analyte_casrn,
#                    Species = species,
#                    Reference = fk_extraction_document_id,
#                    Media = conc_medium_normalized,
#                    Route = administration_route_normalized,
#                    Dose = invivPK_dose_level,
#                    Dose.Units = "mg/kg",
#                    Subject_ID = fk_subject_id,
#                    Series_ID = fk_series_id,
#                    Study_ID = fk_study_id,
#                    ConcTime_ID = conc_time_id,
#                    N_Subjects = n_subjects_normalized,
#                    Weight = weight_kg,
#                    Weight.Units = "kg",
#                    Time = time_hr,
#                    Time.Units = "hours",
#                    Value = invivPK_conc,
#                    Value.Units = "mg/L",
#                    Value_SD = invivPK_conc_sd,
#                    LOQ = invivPK_loq
#                  ))

## -----------------------------------------------------------------------------
# names(cvt)

## ----eval = FALSE-------------------------------------------------------------
# Reference = as.character(
#   ifelse(
#     is.na(
#       documents_reference.id
#     ),
#     documents_extraction.id,
#     documents_reference.id
#   )
# )
# 

## ----eval=FALSE---------------------------------------------------------------
# Value.Units = "mg/L"

## -----------------------------------------------------------------------------
# get_status(my_pk)

## -----------------------------------------------------------------------------
# my_pk <- pk(my_data)

## ----eval = FALSE-------------------------------------------------------------
# my_pk <-   pk(my_data) +
#   facet_data(facets = dplyr::vars(Chemical, Species))

## ----eval = FALSE-------------------------------------------------------------
# my_pk <- pk(my_data)

## ----eval = FALSE-------------------------------------------------------------
# my_pk <- pk(my_data) + settings_preprocess()

## -----------------------------------------------------------------------------
# formals(settings_preprocess)

## ----eval = FALSE-------------------------------------------------------------
# my_pk <- pk(my_data)

## ----eval = FALSE-------------------------------------------------------------
# my_pk <- pk(my_data) +
#   settings_preprocess() +
#   settings_data_info()

## -----------------------------------------------------------------------------
# formals(settings_data_info)

## ----eval = FALSE-------------------------------------------------------------
# my_pk <- pk(my_data) +
#   settings_preprocess() +
#   settings_data_info() +
#   settings_optimx()

## -----------------------------------------------------------------------------
# formals(settings_optimx)

## ----eval = FALSE-------------------------------------------------------------
# my_pk <- pk(my_data) +
#   settings_preprocess() +
#   settings_data_info() +
#   settings_optimx() +
#   scale_conc()

## -----------------------------------------------------------------------------
# formals(scale_conc)

## ----eval = FALSE-------------------------------------------------------------
# my_pk <- pk(my_data) +
#   settings_preprocess() +
#   settings_data_info() +
#   settings_optimx() +
#   scale_conc() +
#   scale_time()

## -----------------------------------------------------------------------------
# formals(scale_time)

## ----eval = FALSE-------------------------------------------------------------
# my_pk <- pk(my_data) +
#   settings_preprocess() +
#   settings_data_info() +
#   settings_optimx() +
#   scale_conc() +
#   scale_time() +
#   stat_model()

## -----------------------------------------------------------------------------
# formals(stat_model)

## ----eval = FALSE-------------------------------------------------------------
# my_pk <-pk(my_data) +
#   settings_preprocess() +
#   settings_data_info() +
#   settings_optimx() +
#   scale_conc() +
#   scale_time() +
#   stat_model() +
#   stat_error_model()

## -----------------------------------------------------------------------------
# formals(stat_error_model)

## ----eval=FALSE---------------------------------------------------------------
# hier_pk <- my_pk +
#   facet_data(ggplot2::vars(Chemical, Species)) +
#   stat_error_model(error_group = vars(Chemical, Species, Reference))

## ----eval=FALSE---------------------------------------------------------------
# pooled_pk <- my_pk +
#   facet_data(dplyr::vars(Chemical, Species)) +
#   stat_error_model(error_group = vars(Chemical, Species))

## ----eval=FALSE---------------------------------------------------------------
# separate_pk <- my_pk +
#   facet_data(dplyr::vars(Chemical, Species, Reference)) +
#   stat_error_model(error_group = vars(Chemical, Species, Reference))

## ----eval = FALSE-------------------------------------------------------------
# my_pk <- pk(data = my_data)
# my_pk <- do_fit(my_pk)

## -----------------------------------------------------------------------------
# my_pk <- pk(my_data) +
#     #instructions to use an error model that puts all observations in the same group
#   stat_error_model(error_group = ggplot2::vars(Chemical, Species)) +
#   #instructions for concentration scaling/transformation
#   scale_conc(dose_norm = TRUE,
#              log10_trans = TRUE) +
#   #instructions for time rescaling
#   scale_time(new_units = "auto") +
#   #instructions to use only one method for optimx::optimx()
#   settings_optimx(method = "L-BFGS-B") +
#   #instructions to impute missing LOQs slightly differently
#   settings_preprocess(calc_loq_factor = 0.5)

## -----------------------------------------------------------------------------
# get_data_original(my_pk)

## -----------------------------------------------------------------------------
# get_mapping(my_pk)

## -----------------------------------------------------------------------------
# get_status(my_pk)

## -----------------------------------------------------------------------------
# get_data_group(my_pk)

## -----------------------------------------------------------------------------
# get_settings_preprocess(my_pk)

## -----------------------------------------------------------------------------
# get_settings_data_info(my_pk)

## -----------------------------------------------------------------------------
# get_settings_optimx(my_pk)

## -----------------------------------------------------------------------------
# get_scale_conc(my_pk)

## -----------------------------------------------------------------------------
# get_scale_time(my_pk)

## -----------------------------------------------------------------------------
# get_stat_model(my_pk)

## -----------------------------------------------------------------------------
# get_stat_error_model(my_pk)

## -----------------------------------------------------------------------------
# my_pk <- my_pk +
#   scale_conc(dose_norm = TRUE, log10_trans = FALSE)

## -----------------------------------------------------------------------------
# get_scale_conc(my_pk)

## -----------------------------------------------------------------------------
# my_pk <- my_pk + stat_model(model = "model_1comp")

## -----------------------------------------------------------------------------
# get_stat_model(my_pk)

## -----------------------------------------------------------------------------
# my_pk <- pk(data = my_data)
# names(my_pk)

## -----------------------------------------------------------------------------
# my_pk <- do_preprocess(my_pk)

## -----------------------------------------------------------------------------
# names(my_pk)

## -----------------------------------------------------------------------------
# get_data(my_pk) %>% glimpse()

## -----------------------------------------------------------------------------
# my_pk <- do_data_info(my_pk)

## -----------------------------------------------------------------------------
# names(my_pk)

## -----------------------------------------------------------------------------
# my_data_info <- get_data_info(my_pk) #extracts the `data_info` element as a named list
# names(my_data_info)

## -----------------------------------------------------------------------------
# my_nca <- get_nca(my_pk) #extracts the `data_info$nca` element specifically

## -----------------------------------------------------------------------------
# my_pk <- do_prefit(my_pk)

## -----------------------------------------------------------------------------
# names(my_pk)

## -----------------------------------------------------------------------------
# get_prefit(my_pk) %>% names()

## -----------------------------------------------------------------------------
# my_pk$prefit$stat_error_model$sigma_DF %>% glimpse()

## -----------------------------------------------------------------------------
# my_pk$prefit$par_DF %>% glimpse()

## -----------------------------------------------------------------------------
# my_pk$prefit$fit_check %>% glimpse()

## -----------------------------------------------------------------------------
# my_pk <- do_fit(my_pk)

## -----------------------------------------------------------------------------
# names(my_pk)

## -----------------------------------------------------------------------------
# get_fit(my_pk) %>% glimpse()

## -----------------------------------------------------------------------------
# my_pk <- pk(data = my_data) + #initialize a `pk` object
#   stat_model(model = c("model_flat",
#                        "model_1comp",
#                        "model_2comp")) + #add PK models to fit
#   settings_optimx(method = "L-BFGS-B") #use only this optimx::optimx() algorithm
# 
# get_status(my_pk) #status is 1

## -----------------------------------------------------------------------------
# my_pk <- do_fit(my_pk)
# get_status(my_pk)

## -----------------------------------------------------------------------------
# AIC(my_pk, suppress.messages = TRUE)

## -----------------------------------------------------------------------------
# my_pk <- my_pk + scale_conc(dose_norm = TRUE)

## -----------------------------------------------------------------------------
# get_status(my_pk)

## ----error = TRUE-------------------------------------------------------------
try({
# coef(my_pk) #throws an error
})

## -----------------------------------------------------------------------------
# my_pk <- do_fit(my_pk)
# get_status(my_pk)

## -----------------------------------------------------------------------------
# AIC(my_pk, suppress.messages = TRUE)

## -----------------------------------------------------------------------------
# #fit a pk object
# my_pk <- pk(data = my_data) +
#   settings_preprocess(suppress.messages = TRUE)
# 
# my_pk <- do_fit(my_pk)
# 
# get_status(my_pk) #status is now 5
# 
# #copy it to a new variable
# my_pk_new <- my_pk
# 
# #and modify scale_conc() on the new copy
# my_pk_new <- my_pk + scale_conc(dose_norm = TRUE)
# 
# get_status(my_pk_new) #status has been reset to 1 for the new copy
# get_status(my_pk) #but status of the original is still 5

## -----------------------------------------------------------------------------
# #suppress messages
# my_pk <- my_pk + settings_preprocess(suppress.messages = TRUE)
# 
# #dose normalization, no log transformation
# pk1 <- my_pk + scale_conc(dose_norm = TRUE, log10_trans = FALSE)
# 
# #log transformation, no dose normalization
# pk2 <- my_pk + scale_conc(dose_norm = FALSE, log10_trans = TRUE)
# 
# #both normalization and log transformation
# pk3 <- my_pk + scale_conc(dose_norm = TRUE, log10_trans = TRUE)
# 
# #do fits
# pk1 <- do_fit(pk1, n_cores = parallel::detectCores()-1)
# pk2 <- do_fit(pk2, n_cores = parallel::detectCores()-1)
# pk3 <- do_fit(pk3, n_cores = parallel::detectCores()-1)

## -----------------------------------------------------------------------------
# plot(pk1) %>% pull(final_plot)

## -----------------------------------------------------------------------------
# plot(pk2) %>% pull(final_plot)

## -----------------------------------------------------------------------------
# plot(pk3) %>% pull(final_plot)

## -----------------------------------------------------------------------------
# plot(pk1, use_scale_conc = list(dose_norm = TRUE,
#                                 log10_trans = FALSE)) %>% pull(final_plot)

## -----------------------------------------------------------------------------
# plot(pk2, use_scale_conc = list(dose_norm = TRUE,
#                                 log10_trans = FALSE)) %>% pull(final_plot)

## -----------------------------------------------------------------------------
# plot(pk3, use_scale_conc = list(dose_norm = TRUE,
#                                 log10_trans = FALSE)) %>% pull(final_plot)

## -----------------------------------------------------------------------------
# pk_hier <- my_pk +
#   stat_error_model(error_group = vars(Chemical, Species, Reference))
# 
# pk_pool <- my_pk +
#   stat_error_model(error_group = vars(Chemical, Species))
# 
# 
# pk_hier <- do_fit(pk_hier, n_cores = parallel::detectCores()-1)
# pk_pool <- do_fit(pk_pool, n_cores = parallel::detectCores()-1)

## -----------------------------------------------------------------------------
# plot(pk_hier, use_scale_conc = list(dose_norm = TRUE,
#                                 log10_trans = FALSE)) %>% pull(final_plot)

## -----------------------------------------------------------------------------
# plot(pk_pool, use_scale_conc = list(dose_norm = TRUE,
#                                 log10_trans = FALSE)) %>% pull(final_plot)

## -----------------------------------------------------------------------------
# `model_1comp`

## -----------------------------------------------------------------------------
# `model_2comp`