---
title: 'Noncompartmental Analysis'
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Noncompartmental Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message=FALSE, eval=FALSE)
require(ggplot2)
require(rhandsontable)
require(flextable)
require(ruminate)
require(dplyr)
# Determining if ubiquity is installed
if(system.file(package="ubiquity") == ""){
  ubiquity_found = FALSE
} else {
  require(ubiquity)
  ubiquity_found = TRUE
}
if(system.file(package="gridExtra") == ""){
  gridExtra_found = FALSE
} else {
  require(gridExtra)
  gridExtra_found = TRUE
}

# The presim variable will contain presimulated data when eval is set to true
presim_loaded = FALSE
```

```{r echo=FALSE, results=FALSE}
presim= list()
if(file.exists("NCA_presim.RData")){
  file.remove("NCA_presim.RData")
}
```

```{r echo=FALSE, results=FALSE, eval=TRUE}
if(file.exists("NCA_presim.RData")){
  load("NCA_presim.RData")
  presim_loaded = TRUE
}

NCA_yaml = yaml::read_yaml(system.file(package="ruminate", "templates","NCA.yaml"))
```


# Introduction

Non-compartmental analysis (NCA) is a simple and quick method for evaluating the exposure of a drug. It allows you to evaluate things like linearity and in vivo exposure. To illustrate this consider an antibody given in a subcutaneous injection. The actual exposure profile a patient might experience is given in the solid black line in the left panel. But we don't yet have the ability to sample in a continuous fashion. What we might observer is given by the blue points.

```{r warning=FALSE, message=FALSE, echo=FALSE, error=FALSE, results="hide", fig.width=8, fig.height=4}
polydf     = NULL
som_smooth = NULL
som_sample = NULL



if(ubiquity_found){
  system_new(system_file="mab_pk", overwrite=TRUE, output_directory=tempdir())
  
  cfg = build_system(system_file         =file.path(tempdir(), "system.txt"),
                     output_directory    =file.path(tempdir(), "output"),
                     temporary_directory =file.path(tempdir(), "transient"))
  
  parameters = system_fetch_parameters(cfg)
  cfg = system_zero_inputs(cfg)
  cfg = system_set_bolus(cfg, state   ="At", 
                              times   = c(  0.0),  #  day
                              values  = c(400.0))  #  mg
  
  cfg=system_set_option(cfg, group  = "simulation",
                             option = "output_times",
                             linspace(0,30,100))
  
  som_smooth = run_simulation_ubiquity(parameters, cfg) 
  som_smooth$simout$time_C_ng_ml = som_smooth$simout$C_ng_ml*som_smooth$simout$ts.days
  
  cfg=system_set_option(cfg, group  = "simulation", 
                              option = "include_important_output_times", 
                              value  = "no")
  
  cfg=system_set_option(cfg, group  = "simulation",
                             option = "output_times",
                             c(0,.25, .5, 1, 2,7,14,21,28))
  som_sample = run_simulation_ubiquity(parameters, cfg) 
  som_sample$simout$time_C_ng_ml = som_sample$simout$C_ng_ml*som_sample$simout$ts.days
  
  
  for(tidx in 1:(nrow(som_sample$simout)-1)){
    xv = c(som_sample$simout$ts.days[tidx], som_sample$simout$ts.days[tidx+1],  som_sample$simout$ts.days[tidx+1], som_sample$simout$ts.days[tidx] )  
    yvC = c(som_sample$simout$C_ng_ml[tidx], som_sample$simout$C_ng_ml[tidx+1], 0, 0)  
    yvTC = c(som_sample$simout$time_C_ng_ml[tidx], som_sample$simout$time_C_ng_ml[tidx+1], 0, 0)  
    tmpdf = data.frame(xv = xv, yvC=yvC, yvTC=yvTC, sp=tidx)
    if(is.null(polydf)){
      polydf = tmpdf 
    } else {
      polydf = rbind(tmpdf, polydf)
    }
  }
}
```

```{r results="hide", warning=FALSE, echo=FALSE}
# When eval is set to TRUE we save the presimulated results
  presim$plots$som_smooth = som_smooth
  presim$plots$som_sample = som_sample
  presim$plots$polydf     = polydf
```

```{r results="hide", warning=FALSE, echo=FALSE, eval=TRUE}
if(presim_loaded){
  som_smooth = presim$plots$som_smooth 
  som_sample = presim$plots$som_sample 
  polydf     = presim$plots$polydf
}
```

```{r warning=FALSE, message=FALSE, echo=FALSE, error=FALSE, eval=TRUE, results="hide", fig.width=8, fig.height=4}
if(ubiquity_found){
  p_C = ggplot()
  p_C = p_C +  geom_line(data=som_smooth$simout, aes(x=ts.days, y=C_ng_ml))
  p_C = p_C + geom_point(data=som_sample$simout, aes(x=ts.days, y=C_ng_ml), color="blue")
  p_C = p_C + xlab("Time") + ylab("Concentration") + ggtitle("AUC")
  p_C = p_C +  geom_polygon(data=polydf, aes(x=xv, y=yvC, group=sp), color="blue", linetype='dashed', fill="lightblue")
  p_C = p_C + theme(plot.title = element_text(hjust = 0.5))
  p_C = prepare_figure(fo=p_C, purpose="shiny")
  
  
  p_TC = ggplot()
  p_TC = p_TC +  geom_line(data=som_smooth$simout, aes(x=ts.days, y=time_C_ng_ml))
  p_TC = p_TC + geom_point(data=som_sample$simout, aes(x=ts.days, y=time_C_ng_ml), color="blue")
  p_TC = p_TC + xlab("Time") + ylab("Time x Concentration") + ggtitle("AUMC")
  p_TC = p_TC +  geom_polygon(data=polydf, aes(x=xv, y=yvTC, group=sp), color="blue", linetype='dashed', fill="lightblue")
  p_TC = p_TC + theme(plot.title = element_text(hjust = 0.5))
  p_TC = prepare_figure(fo=p_TC, purpose="shiny")
  
  
  p_AUC = p_C
  p_AUMC = p_TC
} else {
  p_AUC = ggplot2::ggplot()
  p_AUMC = ggplot2::ggplot()
}
```

```{r warning=FALSE, message=FALSE, echo=FALSE, error=FALSE, eval=TRUE, fig.width=8, fig.height=3.5, fig.alt="Two figures: left panel showing PK profile and right showing integration of the PK profile."}
if(gridExtra_found & ubiquity_found){
  gridExtra::grid.arrange(p_AUC, p_AUMC, ncol=2)
}
```

Generally NCA will determine the following directly from the data:

* ``Cmax`` - Maximum observed concentration (units=concentration)
* ``Tmax`` - The time where the maximum concentration was observed (units=time)
* ``AUC`` - The area under the curve  ($units=time \times concentration$)
* ``AUMC`` - The area under the first moment curve  ($units=time^2 \times concentration$)

These properties are all based on observational data. So the ``Cmax`` and ``Tmax`` will most certainly not be at the actual maximum concentration but as long as we sample judiciously it will give us a good approximation. Similarly, the calculated ``AUC`` and ``AUMC`` will be different than the actual values. To calculate the areas you need to dig back into your calculus text books to the trapezoid method. Basically each sampling interval is a trapezoid and the area of each is calculated and added up for all of the ``n`` samples:

$$
\begin{eqnarray}
AUC  &=& \int_0^{t_f} Cdt &\approx  \sum_{i=1}^{n-1}{\frac{C_i+C_{i+1}}{2}\times (t_{i+1}-t_{i})}  \\ 
AUMC &=& \int_0^{t_f} t\times Cdt &\approx  \sum_{i=1}^{n-1}{\frac{t_iC_i+t_{i+1}C_{i+1}}{2}\times (t_{i+1}-t_{i})}
\end{eqnarray}
$$

This can be done in Excel pretty easily. Depending on the data and the analysis other properties can be calculated. For example we can calculate the clearance, mean residence time, steady-state volume of distribution and terminal half-life:

* Clearance: $CL = \frac{Dose}{AUC}$
* Mean residence time: $MRT = \frac{AUMC}{AUC}$
* Steady state volume of distribution: $V_{ss} = MRT \times CL$
* Half-life: Terminal slope of the natural log of the data

Properties like AUC and AUMC can also be be calculated using extrapolation from the last time point to infinity to account for data beyond the observations at hand. The subsequent values of clearance, volumes of distribution, etc can also change with extrapolation. 

There is a lot of nuance associated with these calculations, and it is good to rely on software that focuses on this type of analysis. The `{PKNCA}` package has been developed with this in mind. The ruminate app provides an interface to functionality found in `{PKNCA}` and will generate code you can use to get started with `{PKNCA}`.

# Data format

To run NCA you need your data formatted in a specific way. The names of the columns can be defined by the user. Some columns are required, others are optional, and others are only required depending on how the user wants to specify dosing. See below for a general description of the columns followed by more information about how dosing is inferred: 

* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_id"]]`:**        `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_id"]][["tooltip"]]`
* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_time"]]`:**      `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_time"]][["tooltip"]]`
* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_ntime"]]`:**     `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_ntime"]][["tooltip"]]`
* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_dose"]]`:**      `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_dose"]][["tooltip"]]`
* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_conc"]]`:**      `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_conc"]][["tooltip"]]`
* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_route"]]`:**     `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_route"]][["tooltip"]]`
* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_group"]]`:**     `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_group"]][["tooltip"]]`
* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_dur"]]`:**       `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_dur"]][["tooltip"]]`
* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_cycle"]]`:**     `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_cycle"]][["tooltip"]]`
* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_evid"]]`:**      `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_evid"]][["tooltip"]]`
* **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_analyte"]]`:**   `r NCA_yaml[["MC"]][["formatting"]][["select_ana_col_analyte"]][["tooltip"]]`


Note that  **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_dose"]]`**  and  **`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_conc"]]`** need to have the same mass units.


## Dosing information

Dosing will be extracted from the datasets, and it can be incorporated in two ways. One is that the dosing information is stored in columns of the dataset. This is the simplest and most straight forward way of storing the dosing information. The other method is that dosing is specified as individual records (rows) in the dataset (if your data is formatted for NONMEM, you can use this). 

### Determining dosing from columns

Beyond the required columns to infer dosing from columns you must also have the following columns defined:
**`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_dose"]]`**,
**`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_cycle"]]`**, and
**`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_ntime"]]`**.


### Determining dosing from rows (NONMEM format)

If your dataset is formatted for NONMEM you can extract the dosing from those rows that contain the dosing records. For this case it is expected you will have the following columns in your dataset:
**`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_dose"]]`** (AMT in NONMEM) and
**`r NCA_yaml[["MC"]][["labels"]][["select_ana_col_evid"]]`** (EVID in NONMEM).

Some notes on dosing information. For each interval where you want to calculate dose-dependent parameters, you will need to provide information about dosing in those intervals. In a multiple dose setting where you have intervals where you do not have enough PK data to run NCA it, is not necessary to provide dosing information for those intervals. For example if you have four weekly doses but only have intensive sampling on weeks 1 and 4 but peak and trough on weeks 2 and 3, you do not have to provide dosing information for weeks 2 and 3. However there are certain parameters that `{PKNCA}` will provide if that dosing and intermediate sampling information is provided (e.g. time to steady-state). So if all you are interested in is parameters for intervals with intensive sampling, then you don't need to worry about having dosing information for the other intervals. 

# Manually flagging data

It is possible to manually flag data for the following purposes:

* Censoring: This will remove data from the analysis completely (flag: `censor`)
* Specifying half-life: For a given interval/analyte combination only these flagged points will be used to determine the half-life (flag: `hlin`)
* Exclusion from half-life: Records with this flag (`hlex`) will not be used in the half-life calculation but will otherwise be included in the analysis

Note that if you have an interval/analyte combination with points specified to be included in the half-life (`hlin`), any flags for exclusion (`hlex`) in this interval/analyte combination will be ignored. At the scripting level the function `flag_nca_ds()` is used to apply flags to the dataset. See the examples section of that function documentation for an example of how flags are applied. 

In the Shiny App manual flags can be applied by:

* First running the analysis with no flags
* Navigate to the figures tab 
* Select the individual profiles
* Select the interactive view
* Select the manual flag to apply
* Add any notes you want to associated with the flagged records
* Select the points to apply that flag

Once manual flags have been applied it is necessary to rerun the analysis. These points will be carried through to the reported tables.

# NCA details

## NCA analysis intervals

If you select an interval for analysis, it is necessary that you have observations at the beginning and end of that interval. If you do not have observations at dosing times (e.g. you dose at time zero but your first observation is at 15 minutes post-dose) and you want to use the dosing time as the beginning of your interval, you can have an observation that is BQL or missing (`NA`) at the dosing time (zero).

```{r warning=FALSE, message=FALSE, echo=FALSE}
save(presim, file="NCA_presim.RData")
```


## NCA parameters 

The following shows the different NCA parameters that can be calculated. The PKNCA parameter is value used in PKNCA, and the App Parameter Name is the value shown in the app. 

```{r warning=FALSE, message=FALSE, echo=FALSE, eval=TRUE}
NCA_meta = NCA_fetch_PKNCA_meta()
NCA_meta = NCA_meta[["parameters"]] |>
  dplyr::select(-data_type) |>
  dplyr::rename("PKNCA Parameter"    = parameter)    |>
  dplyr::rename("App Parameter Name" = pretty_name)  |>
  dplyr::rename("Type of Units"      = unit_type)    |>
  dplyr::rename("Description"        = desc)      

DT::datatable(NCA_meta)
```


# Configuration File

```{r echo=FALSE, message=FALSE, warning=FALSE, eval=TRUE,  style="max-height: 100px;", comment=""}
#yaml= file.path(system.file(package="ruminate"), "templates", "NCA.yaml")
cat(readLines(file.path(system.file(package="ruminate"), "templates", "NCA.yaml")), sep="\n")
```
