## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = FALSE ) ## ----setup-------------------------------------------------------------------- # devtools::load_all() # library(dplyr) # library(tidyr) # library(mirai) ## ----------------------------------------------------------------------------- # # select data for a single chemical and species # my_data <- subset(cvt, # analyzed_chem_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 ## ----------------------------------------------------------------------------- # daemons(4L) # everywhere(devtools::load_all()) # system.time(my_pk <- do_fit(my_pk)) # dameons(0L) ## ----eval = FALSE------------------------------------------------------------- # my_pk <- pk(data = my_data) # # my_pk <- do_fit(my_pk) ## ----------------------------------------------------------------------------- # coef(my_pk) ## ----------------------------------------------------------------------------- # coef(my_pk) |> # slice_head(n = 1) |> # unnest_wider(coefs_vector) |> # glimpse() ## ----------------------------------------------------------------------------- # coef(my_pk) |> # slice_head(n = 5) |> # unnest_wider(coefs_vector) |> # 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( # Chemical = analyzed_chem_dtxsid, # Chemical_Name = analyzed_chem_name_original, # DTXSID = analyzed_chem_dtxsid, # CASRN = analyzed_chem_casrn, # Species = species, # Reference = fk_extraction_document_id, # Media = conc_medium_normalized, # Route = administration_route_normalized, # Dose = dose_level_normalized, # 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 = conc, # Value.Units = "mg/L", # Value_SD = conc_sd, # LOQ = 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(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() + # stat_nca_group() ## ----eval = FALSE------------------------------------------------------------- # my_pk <- pk(my_data) + # settings_preprocess() + # stat_nca_group() + # settings_optimx() ## ----------------------------------------------------------------------------- # formals(settings_optimx) ## ----eval = FALSE------------------------------------------------------------- # my_pk <- pk(my_data) + # settings_preprocess() + # stat_data_info() + # settings_optimx() + # scale_conc() ## ----------------------------------------------------------------------------- # formals(scale_conc) ## ----eval = FALSE------------------------------------------------------------- # my_pk <- pk(my_data) + # settings_preprocess() + # stat_nca_group() + # settings_optimx() + # scale_conc() + # scale_time() ## ----------------------------------------------------------------------------- # formals(scale_time) ## ----eval = FALSE------------------------------------------------------------- # my_pk <- pk(my_data) + # settings_preprocess() + # stat_nca_group() + # settings_optimx() + # scale_conc() + # scale_time() + # stat_model() ## ----------------------------------------------------------------------------- # formals(stat_model) ## ----eval = FALSE------------------------------------------------------------- # my_pk <- pk(my_data) + # settings_preprocess() + # stat_nca_group() + # 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_nca_group(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` ## ----model-setting------------------------------------------------------------ # new_2comp <- set_params_optimize(model_2comp, params = c("k12", "k21")) |> # set_params_starts(starts = list(V1 = 1)) # new_2comp <- adjust_model_name(new_2comp) # # new_pk <- my_pk + # stat_model(model = c("model_flat", "model_2comp", "new_2comp")) + # settings_optimx(method = "bobyqa") # # new_pk <- do_fit(new_pk) # # get_fit(new_pk) |> # filter(model == "new_2comp") |> # glimpse()