Select data to be fit

The first requirement is that the concentration vs. time data to be fit must be available in R in a data.frame format. You may import data from any source you like – external databases, Excel spreadsheets, CSV files, etc. – so long as it can be put into a data.frame format where each row is a single concentration vs. time observation.

invivoPKfit contains a built-in data set from the Concentration vs. Time Database (CvTdb) that is already in a compatible data.frame format: it is in the variable cvt that is loaded with the package. We will use this built-in data for our example.

For the sake of keeping run time reasonable, this example will include data on only two of the chemicals available in the full cvt dataset.

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

Initialize a pk object

Now, initialize a pk object: the basic unit of an invivoPKfit analysis. This object contains the data to be fit, along with a set of instructions for how to do the model fitting. (Note that none of the instructions are actually carried out at this stage: no fitting is done until you explicitly request it.)

my_pk <- pk(data = my_data)

When you call pk(), a default set of instructions will be supplied: how to divide the data into groups; how to pre-process the data; what PK models to fit to the data; what data transformations to apply; what error model to assume during data fitting; etc. More details on the instructions and how to change them will be provided below.

For now, let’s accept the default instructions. It will be easier to understand the instructions themselves if you first understand the overall steps of the model fitting workflow.

The model fitting workflow has four steps.

Step 1: Pre-processing data

The first is data pre-processing. Here, variable names are harmonized, and the data set is cleaned, filtered, and missing values are imputed when necessary.

Some common harmonized variable names that will be referenced in this document are:

• Chemical: the DTXSID for the analyte of a CvT experiment
  
• Species: the species of the subject in a CvT experiment
  
• Route: the route by which the Chemical was administered (iv or oral)
  
• Media: the medium which was collected for analysis (blood or plasma)

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

#preprocess data
my_pk <- do_preprocess(my_pk)

This step adds the pre-processed data to the pk object in a new element named data.

Step 2: Data information

Next, some summary information is gathered about each group of the data. These are things like the number of detects/non-detects, broken out by the different routes and media present in the data. This step also includes non-compartmental analysis of the data, to calculate things like the empirical Cmax (maximum concentration), empirical tmax (time of maximum concentration), and empirical AUC (area under the concentration-time curve). Instructions for this “data info” step are given in my_pk$settings_data_info.

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

This step adds the data summary info to the pk object in a new element named data_info.

Step 3: Pre-fitting

Next, a “pre-fitting” step is performed. This step looks at each group of the data, along with the PK models to be fitted, and determines the following:

• Which of the PK model parameters are to be fitted for each group. For example, if there is no oral data in a group, then parameters that define oral absorption rate and oral bioavailability will not be fitted.
• The lower and upper bounds for the PK model parameters to be fitted for each group, which may be constant or may be derived from the data in some way.
• A starting guess at the value for each PK model parameter to be fitted for each group, which may be constant or may be derived from the data in some way.

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

This step adds the parameter lower/upper bounds and starting guesses to the pk object.

Step 4: Model fitting

Finally, the model fitting can be performed. Each group of the data is fitted separately. By default, three different PK models are fitted to each group: a “flat” model, a 1-compartment model, and a 2-compartment model.

#model fitting
my_pk <- do_fit(my_pk)

This step adds the fitted model parameter values to the pk object. If multiple models were fitted, and/or if multiple optimization algorithms (“methods”) were used to do the fitting, then one set of fitted model parameter values is added for each model and each optimization algorithm (“method”).

Note that model fitting can be sped up substantially by using parallel processing. This is done simply by supplying argument n_cores to do_fit().

Time without parallel processing:

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

Time with parallel processing:

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

Fast-forward to model fitting

Note that you coan just skip straight to the last step. If you call do_fit() on a newly-initialized pk object, all of the preceding steps will be done automatically. The following code will give you exactly the same results as going through each step individually as above.

my_pk <- pk(data = my_data) 

my_pk <- do_fit(my_pk)

Post-fitting: Getting information about the fitted models

Once fitting has been completed, then you can use some familiar methods to extract information about the fitted model(s).

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)

Original data

The data_original element contains the data set you provided in argument data. This is the data set that will be fitted. In this case, it is a subset of CvTdb.

head(my_pk$data_original)
print(my_pk)

Variable mapping

The mapping element specifies how the variable names in the original data should be translated to the harmonized/standardized variable names that are used internally to invivopkfit. This is analogous to providing an aes mapping in ggplot2 – and in fact, it uses the same aes() function as ggplot2!

In ggplot2, you would specify which variables in your data map to the standardized “aesthetic” variables x, y, color, shape, etc., using a command like this:

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

In invivopkfit, you specify which variables in your data map to the standardized variables Chemical, Species, Time, Dose, etc., using a command like this:

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

The above command shows the default mapping that will be applied if you don’t specify a mapping argument when you call pk. It’s a series of “standardized = original” pairs, where the standardized invivoPKfit internal variable name is on the left-hand side, and the original data variable name is on the right-hand side. Here, the original data variable names are all from CvTdb. These are variables appearing in the invivoPKfit built-in data object cvt, which contains a subset of CvTdb data.

names(cvt)

At minimum, you will need to specify mappings for Chemical, Species, Reference, Media, Route, Dose, Time, Value, and LOQ. Most of these names are self-explanatory. Media refers to the biological media in which concentration is measured (e.g. blood, plasma, liver tissue). Route refers to the route of administration (e.g. intravenous or oral). Time refers to the time point at which concentration observation was made. Value refers to the measured concentration. LOQ refers to the limit of quantification of the measured concentration.

Reference = as.character(
  ifelse(
    is.na(
      documents_reference.id
    ),
    documents_extraction.id,
    documents_reference.id
  )
)

Value.Units = "mg/L"

Current status

my_pk$status is an integer code representing the current status of the pk object in the overall pk workflow. The workflow steps are numbered like this:

1. Initialize the pk object in pk()
2. Preprocess data in do_preprocess()
3. Get data summary info (including non-compartmental analysis) in do_data_info()
4. Do pre-fitting (get parameter bounds & starting values) in do_prefit()
5. Do model fitting in `do_fit()

You can get the status of a pk object at any time using

get_status(my_pk)

Data grouping

The data may include the results of multiple experiments for multiple chemicals, species, routes of administration, biological media, etc. You control how to split up the data into batches for fitting. Do you want to fit a PK model to all data for a given chemical and species, even if that includes more than one experiment? Do you want to fit separate PK models to each experiment? It’s your choice.

Data grouping instructions are given using a function called facet_data() (named by analogy with ggplot2). You specify data grouping using a call to dplyr::vars(). You specify data grouping using the names of the new, standardized variables that you created using the mapping argument to pk(). Whatever variables you specify, the data will be grouped by unique combinations of those variables.

The default data grouping is the result of calling the function facet_data() with its default arguments. In other words,

my_pk <- pk(my_data)

is the same as

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

This specification is the default grouping: if you just call pk() without adding a + facet_data() specification, the data will be grouped by unique combinations of Chemical and Species.

Settings for data pre-processing

my_pk$settings_preprocess contains the instructions for pre-processing the data. Its default value is the result of calling the function settings_preprocess() with its default arguments. In other words,

my_pk <- pk(my_data)

is the same as

my_pk <- pk(my_data) + settings_preprocess()

Here are the default arguments to settings_preprocess().

formals(settings_preprocess)

Explanations of each argument follow, to show how you can control the data pre-processing settings.

Data information settings

my_pk$settings_data_info provides control settings for getting summary data information, including non-compartmental analysis. Its default value is the result of calling the function settings_data_info() with its default arguments.

In other words,

my_pk <- pk(my_data)

is the same as

my_pk <- pk(my_data) +
  settings_preprocess() +
  settings_data_info()

Here are the default arguments to settings_data_info():

formals(settings_data_info)

Optimizer settings

my_pk$optimx_settings provides control settings for the optimizer used to fit the model(s) to the data, which is optimx::optimx(). It is the result of a call to settings_optimx() with default arguments. In other words, just calling pk() by itself is the same as

my_pk <- pk(my_data) +
  settings_preprocess() +
  settings_data_info() +
  settings_optimx()

Here are the default arguments to settings_optimx():

formals(settings_optimx)

These are simply the arguments for optimx::optimx(): method, itnmax, hessian, control, and ....

In general, you’ll need to choose a value for method that names one of the optimization algorithms that can handle parameters bounds: either "bobyqa" or "L-BFGS-B".

Data transformations (scalings)

my_pk$scales is itself a list with elements conc and time, providing instructions on how to transform concentration and time variables before fitting any models.

my_pk <- pk(my_data) +
  settings_preprocess() +
  settings_data_info() +
  settings_optimx() +
  scale_conc()

formals(scale_conc)

my_pk <- pk(my_data) +
  settings_preprocess() +
  settings_data_info() +
  settings_optimx() +
  scale_conc() +
  scale_time()

formals(scale_time)

Models to be fitted

my_pk$stat_model is a named list, with one item named for each model to be fitted. It is the result of a call to stat_model() with its default arguments. In other words, calling pk() by itself is the same as calling

my_pk <- pk(my_data) +
  settings_preprocess() +
  settings_data_info() +
  settings_optimx() +
  scale_conc() +
  scale_time() +
  stat_model()

Here are the default arguments to stat_model:

formals(stat_model)

The model argument is a character vector (or scalar), listing the names of one or more models to fit to the data. The default value lists the three models that are already implemented and built in to invivopkfit.

Error model to apply during fitting

my_pk$stat_error_model lists the error model to apply during fitting. This is the result of a call to stat_error_model() with its default arguments.

In other words, calling pk() by itself is the same as calling

my_pk <-pk(my_data) +
  settings_preprocess() +
  settings_data_info() +
  settings_optimx() +
  scale_conc() +
  scale_time() +
  stat_model() +
  stat_error_model()

Here are the default arguments to stat_error_model():

formals(stat_error_model)

There is only one default argument: error_group. It is specified, just like data grouping in facet_data(), using a call to dplyr::vars().

This specifies the error model as, essentially, a fixed-effects model. Within each data group, the data are further split into subgroups according to unique combinations of the variables specified in error_group. Residual errors within each error subgroup are assumed to be independent and identically distributed according to a zero-mean normal distribution with a standard deviation that is unique to that error subgroup.

The error standard deviations for each error subgroup are hyperparameters of the model. They will be estimated from the data simultaneously with the model parameters for each data group. The argument error_group defines the number of different error subhgroups, and therefore the numnber of error SD hyperparameters that need to be estimated.

If you specify error_group to be the same as the data group, then that simply means there is only one error group within each data group, and therefore only one error standard deviation for that data group.

Here are three examples:

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

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

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

Check original data

get_data_original(my_pk)

Check mapping

get_mapping(my_pk)

Check status

get_status(my_pk)

Check data grouping

get_data_group(my_pk)

Check settings_preprocess

get_settings_preprocess(my_pk)

Check settings_data_info

get_settings_data_info(my_pk)

Check settings_optimx

get_settings_optimx(my_pk)

Check scale_conc

get_scale_conc(my_pk)

Check scale_time

get_scale_time(my_pk)

Check stat_model

get_stat_model(my_pk)

Check stat_error_model

get_stat_error_model(my_pk)

Replacing instructions with new ones

You can then overwrite any instructions by adding new/different ones. Let’s say you made a mistake – you actually wanted to dose-normalize concentrations, but not log-transform them.

You can change the concentration-scaling instructions by simply adding new ones:

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

Now the new instructions have replaced the old ones.

get_scale_conc(my_pk)

And let’s say you only want to fit the one-compartment model – you don’t want the default behavior, which fits the flat and two-compartment models as well as the one-compartment model. You can do this:

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

Now, “model_flat” and “model_2comp” models have been removed from the list of models to fit.

get_stat_model(my_pk)

Initialization

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

The only items in the pk object are those that contain the original data and instructions.

Pre-processing data

my_pk <- do_preprocess(my_pk)

Re-check the names of the pk object:

names(my_pk)

Notice that a new element data has been added to the end. This contains the pre-processed data. You can access it using the function get_data().

get_data(my_pk) %>% glimpse()

Notice that the variables Value and LOQ have been translated into Conc and Detect. If Value < LOQ then Detect = FALSE; otherwise Detect = TRUE. If Detect = FALSE, then Conc = LOQ; otherwise, Conc = Value. This is just a different way of presenting the same information, but it’s a bit easier to handle with some of the internal calculations of invivoPKfit.

Also, notice that the requested time and concentration transformations have been applied – compare Time and Time.Units vs. Time_trans and Time_trans.Units; compare Conc and Conc.Units vs. Conc_trans and Conc_trans.Units.

Notice that missing LOQ values have been imputed; compare LOQ.orig to LOQ. Similarly, missing Value_SD values have been imputed; compare Value_SD.orig and Value_SD. Note that the standard deviations are put in a new variable Conc_SD to match with the Conc variable. Also, note that the concentration transformations have been applied to the Conc_SD variable as well – see Conc_SD_trans. That way, the transformation concentration standard deviations are in the same units as the transformed concentrations.

(A warning: If you do scale_conc(log10_trans = TRUE), then the log10 transformation will be applied to the concentration SD as well as the concentration. This means you can no longer take Conc_trans + Conc_SD or Conc_trans - Conc_SD as upper and lower bounds. This math is handled correctly internal to invivopkfit, but if you extract the pre-processed data and use it for another purpose, be aware of that caveat!)

Data summary info

my_pk <- do_data_info(my_pk)

names(my_pk)

A new element data_info has been added. This element is itself a named list containing elements data_summary, data_flags, dose_norm_check, and nca. You can extract these elements using the following functions:

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

Pre-fitting

my_pk <- do_prefit(my_pk)

names(my_pk)

An new list element prefit has been added. It can be extracted using the function get_prefit(). It is itself a named list, with one element for every PK model to be fitted to the data.

get_prefit(my_pk) %>% names()

Element prefit$stat_error_model has one element sigma_DF, which gives the starting values and lower/upper bounds for the error-group-specific standard deviations.

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

Element prefit$par_DF gives the starting values and lower/upper bounds for each parameter of each model to be fit.

my_pk$prefit$par_DF %>% glimpse()

Element fit_check gives information about whether each fit should be attempted (i.e., whether there are enough observations).

my_pk$prefit$fit_check %>% glimpse()

Keeping track of workflow status as you change instructions

Let’s say you go through all the steps of the workflow for a pk object, all the way to fitting. And then, you change one of the instructions. For example, you change the concentration scaling, which will change the results of everything from pre-processing the data (step 2) onwards. Now, the fit results stored in the fitted pk object don’t match the concentration scaling stored in the same object. You need to re-do all steps of the workflow.

However, if you had changed the error model instead, you wouldn’t need to re-do all steps. You’d only need to re-do the pre-fitting and fitting steps. The error model doesn’t affect data pre-processing or data summary information steps.

How can you keep track of which steps you need to re-do after changing one or more of the instructions?

The answer is that invivopkfit will automatically check the object status every time you issue a new instruction, and will re-set it back to the point where you need to re-do the steps of the workflow.

Let’s say I’ve created a pk object. I don’t apply any scalings/transformations to concentration or time.

The status of the newly-created pk object is 1.

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

Now I do all the steps up to step 5. (If you just call do_fit(), it will fast-forward through all the steps.)

my_pk <- do_fit(my_pk)
get_status(my_pk)

I can now extract the AIC for my fitted models.

AIC(my_pk, suppress.messages = TRUE) 

But now I realize that I should have dose-normalized the concentration. I do this by adding a scale_conc(). invivoPKfit throws a warning to tell me that the status will be reset back to 1.

my_pk <- my_pk + scale_conc(dose_norm = TRUE)

And in fact, the status is now reset to 1.

get_status(my_pk)

This means that I can’t do anything that requires a fully fitted model – not until I re-do all the workflow steps.

coef(my_pk) #throws an error

In effect, invivopkfit throws out any fit results that were made before I changed the concentration scaling.

But after I re-fit the new model (setting the status back to 5), I can now extract the new AIC values.

my_pk <- do_fit(my_pk)
get_status(my_pk)

AIC(my_pk, suppress.messages = TRUE) 

If I wanted to keep my original fit results, and also create a new pk object that was the same except for the concentration scaling, then I should copy my pk object into a new variable and modify the new one.

#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

Compare fits made using different error models

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)

Many more and detailed examples can be seen in the R Markdown document that gives the code to produce the figures and tables for the manuscript.

Model function requirements

Here are the requirements for each of those functions (which can also be found in help(pk_model)):

Caveats

Currently, some functions in invivopkfit may assume particular units for the data, and therefore particular units for the fitted model parameters. For example, time is generally assumed to be in units of hours. We plan to relax these assumptions in future releases of invivoPKfit.