models.qmd

---
title: "Models"
---

```{r}
#| include: false

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

```{r}
#| label: setup
#| echo: false
#| message: false
#| eval: true

library(Pmetrics)
library(glue)
library(knitr)
r_help <- function(pkg, name) {
    glue::glue("[`{name}`](https://rdrr.io/pkg/{pkg}/man/{name}.html)")
}

gh_help <- function(name) {

    glue::glue("[`{name}`](https://lapkb.github.io/Pmetrics/reference/{name}.html)")

}


pmetrics <- function(){
    knitr::asis_output("[Pmetrics]{style=\"color: #841010; font-family: 'Arial', Arial, sans-serif; font-weight: 900;\"}")
}

library(gt)

custom_table <- function(tab){
  #system.file("extData", tab, package = "Pmetrics") %>%
    read.csv(file.path("Data",tab), na.strings = ".") %>% gt() %>%
    tab_style(
      style = list(
        cell_fill(color = "black"),
        cell_text(color = "white", weight = "bold")
      ),
      locations = cells_column_labels(everything())
    )
}


```

## Introduction

:::{.callout-tip}
Make sure you have run the [tutorial setup code](index.qmd#tutorial-setup) in your R session before copying, pasting and running example code here.
:::

`PM_model` objects are one of two fundamental objects in Pmetrics, along with [PM_data](data.qmd) objects. Defining a `PM_model` allows for fitting it to the data via the `$fit()` method to conduct a population analysis, i.e. estimating the probability distribution of model equation parameter values in the population.

## Creating Models

There are three pathways to creating models in Pmetrics.

-   Use the new Pmetrics Model Builder app, coming soon. The Shiny app is retired.
-   Write a model list in R
-   Load an existing model text file. This is for legacy purposes only; we strongly recommend using model lists in R for all new modeling work.

Each of these use sections to define the model. In the app, the sections correspond to tabs. In the list, they are named elements. In the model file, they are code blocks delimited by "\#" and the name of the block. The sections are largely the same across the three model pathways, and we'll cover the details in this document.

<!-- ### Create a model with the app

Coming soon.

To launch the app, type the following into your console:

```{r}
#| eval: false
#| echo: false
#| label: mod_builder

build_model()
```

You can supply a `PM_data` and/or a `PM_model` as optional arguments to the function, e.g. `PM_model$new(dataEx)` or `PM_model$new(modEx)` or even `PM_model$new(dataEx, modEx)`. The order of data and/or model arguments doesn't matter. Pmetrics can figure out which is which. -->

### Model lists in R

You can write a Pmetrics model list in R. 
<!-- You can also generate a model in the app and copy it as a list to R which can be used in `PM_model$new(<copied_list>)`.  -->
Once in your R script or Quarto document, `PM_model` list objects can be updated simply by copying, pasting, and editing. This makes it very easy to see how your model evolves during development.

The list elements are the building blocks for the model. Building blocks are of two types.

#### Block types

-   **Lists** define *primary parameters*, *covariates*, and *error models*. These portions of the model are lists of specific creator functions and no additional R code is permissible. Note the commas separating the creator functions, `list(` to open the list and `)` to close the list. Names are case-insensitive and are converted to lowercase for Rust.

:::{code-copy='false'}
```{r}
#| label: mod-block-list
# example only; do not execute

# for primary parameters and covariates
block_name = list(     
  var1 = creator(), # var1 is the name of the parameter/covariate   
  var2 = creator(),
  ... # more variables
) # close the list

# for error block
err = list(     
    creator(),   
    creator(),
    ... # as many creators as output equations
) # close the list
```
:::



-   **Functions** define the other parts of the model, including *secondary (global) equations*, *model equations* (e.g. ODEs), *lag time*, *bioavailability*, *initial conditions*, and *outputs*. These parts of the model are defined as R functions without arguments, but whose body contains any permissible R code. Note the absence of arguments between the `()`, the opening curly brace `{` to start the function body and the closing curly brace "`}`" to end the body. Again, all R code will be converted to lowercase prior to translation into Rust.

:::{code-copy='false'}
```{r}
#| label: mod-block-func
# example only; do not execute
block_name = function() { 

      # any valid R code 
      # can use primary or secondary parameters and covariates
      # lines are not separated by commas
      
} # end function and block
```
:::



:::{.callout-note}
All models must have `pri`, `eqn`, `out`, and `err` blocks.
:::

### Model Files {#mb-files}

:::{.callout-warning}
As of `r pmetrics()` 3.0, you can no longer save models as text files. Model text files should be gradually phased out of your workflow, because it is cumbersome to copy and edit files between runs, and it makes it harder to go back and read the R script to understand how the modeling process evolved unless the script is extensively commented. Having the model code directly in the script is simple to edit and far more transparent. We **strongly** recommend using model lists in R for all modeling work.
:::

<!-- You can write the file by hand, or use the app to create a model and then save it to a file.  -->

The format of the .txt model file has been mostly unchanged since `r pmetrics()` v.0.4, and such files are created by hand. Once you have a file, it can be loaded in to R with `PM_model$new("file")`. To facilitate learning the preferred method of defining models directly in R rather than with files, when you create a model from a text file, the R code for the model is copied to the clipboard for easy pasting into your R script or Quarto document.

```{r}
#| label: load-model-file
mod1 <- PM_model$new("src/model.txt") 

```

The above code will copy the following to your clipboard, allowing you to paste it into your script or Quarto document.

:::{code-copy='false'}
```{r}
#| label: pasted-model-code
#| echo: true
mod <- PM_model$new(

  pri = list(
    ka = ab(0.100, 0.900),
    ke = ab(0.001, 0.100),
    v = ab(30.000, 120.000),
    lag1 = ab(0.000, 4.000)
  ),
  cov = list(
    wt = interp(),
    africa = interp("none"),
    age = interp(),
    gender = interp("none"),
    height = interp()
  ),
  lag = function () 
  {
      lag[1] = lag1
  },
  eqn = function () 
  {
      two_comp_bolus
  },
  out = function () 
  {
      y[1] = x[2]/v
  },
  err = list(
    proportional(5, c(0.0, 0.1, -0.0, 0.0))
  )

)
```
:::

Saved models are only text files. You can write files directly yourself, although it is easier and more stable to use the model building app. We will review the format of the files in detail.

**Naming your model files.** The default model file name is "model.txt," but you can call them whatever you wish. 

**Structure of model files.** The model file is a text file with up to 9 blocks, each marked by "\#" followed by a header tag. For each header, only the capital letters are required for recognition by Pmetrics. The blocks can be in any order, and header names are case-insensitive (i.e. the capitalization here is just to show which letters are required). We include a [complete examples](#completeEx).

You can insert comments into your model text file by starting a line with a capital "C" followed by a space. These lines will be removed/ignored in the final Rust code.

## Model Blocks

<!-- *Note*: The model building app has a unique [Model library](#mb-lib) section not found in model files or lists. -->

Blocks common to all the model building pathways include the following. The capital letters are the shorthand way to designate the blocks in the R list format or model text files.

-   [PRImary variables](#mb-pri)
-   [COVariates](#mb-cov)
-   [SECcondary variables](#mb-sec)
-   [INItial conditions](#mb-ini)
-   [FA (bioavailability)](#mb-fa)
-   [LAG time](#mb-lag)
-   [EQuatioNs](#mb-eqn)
-   [OUTputs](#mb-out)
-   [ERRors](#mb-err)

### Model Library {#mb-lib}

<!-- 
![](Images/model_builder/front.png)

Here you can choose from pre-existing models, either which you have created yourself and load with the Previous Model dialogue, or from the Pmetrics Model Library. Use the filters and description search to help select the model you want. Matching models will update in the box at the bottom. If you select one, you'll see a model snapshot on the right.

In the snapshots, B is for bolus inputs, R is for Rate infusions (e.g. intravenous infusions), and Y is for observations. Arrows indicate the flow of drug. The grey compartment 0 is the environment.

Once you choose a model, hit the "Select" button at the bottom to populate remaining model components. Similarly, if you load a Previous Model, model components will be appropriately populated. -->

You can access the Pmetrics Model Library programmatically using the `r gh_help("model_lib")` function. This function returns a data frame of available models with their characteristics, including the model name to be used in the `EQN` block, number of compartments, number of inputs and outputs, and a brief description.

All of the models are pre-loaded as R6 objects. You can get details of the model by simply typing its primary name from the table into the console. For example:
```{r}
#| eval: false
#| label: model-lib-example
#| echo: true

# see the library contents
model_lib()
# an example of a model in the library
two_comp_bolus
```


### PRImary {#mb-pri}
<!-- 
#### App

Coming soon. -->

<!-- ![](Images/model_builder/primary.png)

In this tab, you can choose the number and names of the primary parameters. Some variable names are [reserved](#reserved) for use by Pmetrics and cannot be used as primary variable names.

These are the parameters for which value probability distributions will be estimated. You can specify initial values as ranges or mean/SD. The mean is the mid point of a range and the SD is 1/6 of the range, i.e. 3 SD above and below the mean.

In Legacy Pmetrics it was possible to specify parameters as either fixed to a known value in the entire population (constant) or fixed to some unknown value in the population (fixed but unknown). To make parameters constant in Pmetrics ≥ 3.0, simply include the value in either the #SECondary block or in the equations directly. Because the Legacy optimization process was very suboptimal and expensive for fixed but unknown parameters, we have removed it in Pmetrics ≥ 3.0. To accomplish this, we suggest including such a parameter as a typical random parameter, and then fixing it to a constant value after inspecting the results of a satisfactory run when the parameter was random.

In the bottom you see the Model Diagram and the Model List. These dynamically build as you complete the Model Components. -->

#### List

Primary variables are defined in a `list()`. The list is named "pri", and is comprised of named elements which are the parameters, and the creator function to define each parameter, one parameter to a line.

The two creator functions for primary parameters are `ab` and `msd`. Choose one for each parameter. The first defines the absolute search space for that parameter for NPAG/NPOD, i.e., the range. `msd` is the companion function that specifies the range using a mean and standard deviation (SD), assuming the mean is the midpoint of the range, and the SD equals 1/6th of the range, i.3., there are 3 SD above and below the mean, covering 99.7% of the prior distribution. <!-- For IT2B/RPEM, it defines the mean of the prior as the midpoint  --> <!-- of the range, and the range covers 6 standard deviations, e.g. ±3 SD above and  --> <!-- below the mean, or 99.7% of the prior distribution. `msd` is the companion  --> <!-- function that specifies a mean and SD in IT2B and RPEM.  --> <!-- For NPAG/NPOD, it will be converted in to a range in the reverse fashion as --> <!-- described for `ab`.  -->

<!-- For both specifying functions, `gtz` is an argument to -->

<!-- force the parameter value to be positive, i.e. `gtz=T`, which is the default. -->

<!-- To allow negative parameters, set `gtz=F`. -->

:::{code-copy='false'}
```{r}
#| label: mod-block-pri
# this is only a partial model snippet for illustration; do not execute
mod <- PM_model$new(
  pri = list(
    Ke = ab(0, 5), # mean will be 2.5, SD = 5/6 = 0.83
    V = msd(100, 20) # range will be 100 ± 3*20, i.e. [40, 160]
  )
)
```
:::

#### File

Primary variables are a list of variable names, one name to a line. On each line, following the variable name, is the range for the parameter that defines the search space. These ranges behave slightly differently for NPAG and the simulator.

-   For all engines, the format of the limits is *min, max*.

-   For **NPAG/NPOD**, when *min, max* limits are specified, they are absolute, i.e. the algorithm will not search outside this range.

<!-- - For **IT2B**, the range defines the Bayesian prior distribution of -->

<!--     the parameter values for cycle 1. For each parameter, the mean of -->

<!--     the Bayesian prior distribution is taken as the middle of the range, -->

<!--     and the standard deviation is *xsig*\*range (see [[IT2B -->

<!--     runs]{.underline}](\l)). Adding a plus sign (+) to a line will -->

<!--     prevent that parameter from being assigned negative values. NPAG and -->

<!--     the simulator will ignore the pluses as the ranges are absolute for -->

<!--     these engines. Note that prior to version 1.5.0, this used to be an -->

<!--     exclamation point (!) but to be consistent throughout the model -->

<!--     file, the exclamation point is now used when fixed values are -->

<!--     desired. -->

-   The **simulator** will ignore the ranges with the default value of NULL for the `limits` argument to `r gh_help("PM_sim")`. If the simulator `limits` argument is set to NA, it will mean that these ranges will be used as the limits to truncate the [simulation](simulation.html).

Example:

::: script
#Pri\
KE, 0, 5\
V, 0.01, 100\
KA, 0, 5\
:::

### COVariates {#mb-cov}

Covariates are subject-specific data, such as body weight, contained in the data. The covariate names, which are the column names in the data file or `PM_data` object, must be declared if used in the model object. Once declared, they can be used in secondary variable and differential equations. The names should be the same as in the data file.

By default, missing covariate values are linearly interpolated between existing values, or carried forward if the first value is the only non-missing entry. However, you can suppress interpolation and carry forward the previous value in a piece-wise constant fashion, depending on the way you define the model as described below.

<!-- #### App

Coming soon. -->

<!-- ![](Images/model_builder/covariates.png)

If you launch the Model Builder with a `PM_data` object as an argument, the covariate tab will be pre-populated with the covariates in the data file, as above. The same is true if you include a `PM_model` object as an argument to `PM_model$new()`. If you include both, the covariates in the `PM_data` object will take precedence. You can select any of them to make them piece-wise constant, i.e. the value is held constant between measurements. If left unchecked, covariate values will be linearly interpolated between measurements. -->

#### List

The covariate element is a list whose names are the same as the covariates in the data file, and whose values are the `r gh_help("interp")` creator function to declare how each covariate is interpolated between entries in the data.

The default argument for `interp` is "lm" which means that values will be linearly interpolated between entries, like the R linear model function `r r_help("stats", "lm")`. The alternative is "none", which holds the covariate value the same as the previous entry until it changes, i.e., a carry-forward strategy.

:::{code-copy='false'}
```{r}
#| label:  mod-block-cov
# this is only a partial model snippet for illustration; do not execute
mod <- PM_model$new(
  pri = list(...),
  cov = list(
      wt = interp(),
      age = interp("none")
    )
)
```
:::

Here, `wt` will be linearly interpolated and `age` will be piece-wise constant, i.e. held constant until the next entry in the data file.

#### File

List the covariates by name, one per line in the #COV block. To make a covariate piece-wise constant, include an exclamation point (!) in the declaration line.

**Note** that any covariate relationship to any parameter may be described as the user wishes by mathematical equations and R code, allowing for exploration of complex, non-linear, time-dependent, and/or conditional relationships. These equations are included in any of the blocks which are functions (SEC, EQN, INI, FA, LAG, OUT).

Example:

::: script
#Cov\
wt\
cyp\
IC!
:::

Here, IC will be piece-wise constant and the other two will be linearly interpolated for missing values.

### Secondary {#mb-sec}

Secondary variables are those that are defined by equations that are combinations of primary, covariates, and other secondary variables. Pmetrics does *NOT* generate distributions for the values of secondary variables. They are only used internally to define the model. Variables declared in the secondary block are ***global***, i.e., they are available to all other blocks. This is in contrast to variables declared within other function blocks (INI, FA, LAG, OUT), which are only available within the block they are declared.

<!-- #### App

Coming soon. -->

<!-- ![](Images/model_builder/secondary.png)

If using secondary variables, define them first within this tab. Equation syntax for Secondary equations must be R. Each expression must be on a new line and contain only variables which have been previously defined in the Primary, Covariate, or Secondary blocks.

The image shows examples of secondary variable declarations without conditions. Here are two examples of conditional secondary variables chosen on the basis of sex. The primary variables are Vm, Vf, CLm, and CLf.

```         
V = Vm
if (sex == 1) V = Vf

V = Vm
CL = CLm
if(sex == 1) {
  V = Vf
  CL = CLf
}
``` -->

#### List

Specify each variable equation as a list of unnamed character elements in the same way as for the app. In the example below, V0 is the primary parameter which will be estimated, but internally, the model uses V as V0\*wt, unless age is \>18, in which case weight is capped at 75 kg. It's the same for CL0. Note that the conditional statement is not named.

:::{code-copy='false'}
```{r}
#| label: mod-block-sec
# this is only a partial model snippet for illustration; do not execute
mod <- PM_model$new(
    pri = list(
      CL0 = ab(0, 5),
      V0 = msd(10, 3)
    ),
    cov = list(
      wt = interp(),
      age = interp()
    ),
    sec = function(){
      V = V0*wt
      if (age >18) {V = V0 * 75}
      CL = CL0 * wt
    }
) # end of $new()
```
:::

#### File

Format is the same as for equations in the list.

Example:

::: script
#Sec\
CL = Ke \* V \* wt\*\*0.75\
if (cyp \> 1) {CL = CL \* cyp}
:::

### INItial conditions {#mb-ini}

By default, all model compartments have zero amounts at time 0. However, you can change the default initial condition of any compartment from 0 to something else. It can be an fixed value, primary or secondary variables, or covariate(s), or equations/conditionals based on these.
<!-- 
#### App

Coming soon. -->
<!-- We'll discuss more about the interface in the section on [lag times](#mb-lag). -->

#### List

Initial conditions can be changed as a function called "ini". Each line of the function specifies the compartment amount as "X\[.\] = expression", where "." is the compartment number. Primary and secondary variables and covariates may be used in the expression, as can conditional statements in R syntax. 

:::{code-copy='false'}
```{r}
#| label:  mod-block-ini
# this is only a partial model snippet for illustration; do not execute 
mod <- PM_model$new(
    pri = list(
      Ke = ab(0, 5),
      V = msd(100, 30),
      IC3 = ab(0, 1000)
    ),
    cov = list(
      wt = interp(),
      age = interp(),
      IC2 = interp("none")
    ),
    ini = function(){
      X[1] = IC2*V
      X[3] = IC3
    }
) # end of $new()
```
:::

In the example above, IC2 is a covariate with the measured trough concentration prior to an observed dose and IC3 is a fitted primary parameter specifying an initial amount in unobserved compartment 3.

In the first case, the initial condition for compartment 2 becomes the value of the IC2 covariate (defined in `cov` list) multiplied by the current estimate of V during each iteration. This is useful when a subject has been taking a drug as an outpatient, and comes in to the lab for PK sampling, with measurement of a concentration immediately prior to a witnessed dose, which is in turn followed by more sampling. In this case, IC2 or any other covariate can be set to the initial measured concentration, and if V is the volume of compartment 2, the initial condition (amount) in compartment 2 will now be set to the measured concentration of drug multiplied by the estimated volume for each iteration until convergence.

In the second case, the initial condition for compartment 3 becomes another variable, IC3 defined in the `pri` list, to fit in the model, given the observed data.

#### File

The same example as for the list above is shown below in the format for the file.

::: script
#Ini\
X\[2\] = IC\*V\
X\[3\] = IC3\
:::

### FA (bioavailability) {#mb-fa}

In this tab you can change the default bioavailability of any bolus input from 1 to something else. It can be an equation, primary or secondary variable, or covariate.
<!-- 
#### App

Coming soon. -->

<!-- We'll discuss more about the interface in the section on [lag times](#mb-lag). -->

#### List

Bioavailability for any bolus input can be changed from the default value of 1. Use a function called "fa", where each line is of the form "FA\[.\] = expression", and "." is the input number. Primary and secondary variables and covariates may be used in the expression, as can conditional statements.

:::{code-copy='false'}
```{r}
#| label:  mod-block-fa
# this is only a partial model snippet for illustration; do not execute
mod <- PM_model$new(
    pri = list(
      Ke = ab(0, 5),
      V = msd(100, 30),
      FA1 = ab(0, 1)
    ),
    fa = function(){
      FA[1] = FA1
    }
)
```
:::

#### File

The same example as for the list above is shown below in the format for the file.

::: script
#Fa\
FA\[1\] = FA1
:::

### LAG time {#mb-lag}

The lag time is a delay in absorption for a bolus input.

<!-- #### App

Coming soon. -->

<!-- ![](Images/model_builder/lag.png)

You can change the default delay in absorption of any bolus input from 0 to something else. It can be an equation, primary or secondary variable, or covariate. If you wish to use any of the latter three, select them from the drop down and the equation will pre-populate, as shown in the image. You can edit the equation or write your own equation. This is true for initial conditions and bioavailability in their tabs. -->

#### List

Use a function called "lag", each line of the form "LAG\[.\] = expression", where "." is the input number. Primary and secondary variables and covariates may be used in the expression, as can conditional statements.

:::{code-copy='false'}
```{r}
#| label:  mod-block-lag
# this is only a partial model snippet for illustration; do not execute 
mod <- PM_model$new(
    pri = list(
      Ke = ab(0, 5),
      V = msd(100, 30),
      lag1 = ab(0, 4)
    ),
    lag = function(){
      LAG[1] = lag1
    }
)
```
:::

#### File

The same example as for the list above is shown below in the format for the file.

::: script
#Lag\
LAG\[1\] = Tlag1
:::


### EQuatioNs {#mb-eqn}

These are the equations that define the structural model, i.e., the mathematical expressions that relate input (dose) to output (measurements). Pmetrics has a library of models with algebraic solutions and models can also use differential equations in a format compatible with R. Even algebraic models will contain differential equations for purposes of understanding and plotting the model, but algebraic models from the library include a token which tells Pmetrics how to choose the correct model. Details below.

<!-- #### App

Coming soon. -->

<!-- ![](Images/model_builder/eqn.png)

If you have not loaded a prior model file or selected a model from the library, you can write the differential equations here. The model diagram will update as information in the equations becomes available.

Use `dX[i]` for change in compartment amounts, where i is the compartment number, e.g. dX\[1\] or dX\[2\]. Compartment amounts are referred to as `X[i]`, e.g. X\[1\] or X\[2\]. Use `BOLUS[j]` or `B[j]` for bolus input j and `RATEIV[k]` or `R[k]` for infusion k. The indices j and k correspond to the `INPUT` column in the data file, which is omitted and assumed to be 1 for all doses for the most common scenario of modeling one drug input. The duration of the infusion and total dose is defined in the data. The `DUR` column in the data determines whether a dose is treated as a BOLUS (DUR = 0) or RATEIV (DUR \> 0). Any variable declared in `PRI`, `COV`, or `SEC` may be used in your equations.

**Note:** If you change the equations of an algebraic model loaded from the library, you will receive a warning in the app that you are forcing use of ODE solvers, and your model list (and file) will no longer have the algebraic token. Beware also of changing the parameter names, as the algebraic models must use specifically named parameters. However, it is perfectly acceptable to add parameters or even change primary parameters as long as the original appear in the secondary block. For example, if you wish to make volume dependent on weight, change `V` in the Primary section to something like `V0` with appropriate range or mean/SD, and add `V = V0 * wt` to your Secondary section, assuming you have a covariate named "wt". Because `V` is still defined, Pmetrics will know how to solve the algebraic model. -->

#### List

Specify a model as a function called "eqn", with each element an ordinary differential equation in R syntax. Additional equations can be included. Special elements are as follows:

* `dX[.]` is the change in amount in compartment ".", e.g. `dX[1]` for compartment 1. This is shorthand for $dX/dt$.
* `X[.]` is the amount in compartment ".", e.g. `X[2]` for compartment 2.
* `B[.]` or `BOLUS[.]` is the bolus input number ".", e.g. `B[1]` for bolus input 1. **Note** the the number is not the compartment number, but the number of the input, corresponding to the `INPUT` column in the data file. Which compartment the bolus enters is defined by which equation includes `B[.]`.
* `R[.]` or `RATEIV[.]` is the infusion input number ".", e.g. `R[1]` for infusion input 1, again corresponding to the data `INPUT` number, not the compartment number.

If your data file has dose inputs with `DUR = 0`, these are considered to be instantaneous bolus inputs, and you must use `B[.]` or `BOLUS[.]` to represent them in the equations. If your data file has dose inputs with `DUR > 0`, these are considered to be infusions, and you should use `R[.]` or `RATEIV[.]` to represent them in the equations.

:::{.callout-tip}
Local variables defined within the `eqn` function are not available outside the function. If you wish to use a variable in other blocks, define it in the `sec` function to make it global.
:::


:::{code-copy='false'}
```{r}
#| label:  mod-block-eqn
# this is only a partial model snippet for illustration; do not execute 
mod <- PM_model$new(
    pri = list(
      Ka = ab(0, 5),
      Ke = ab(0, 5),
      V = msd(100, 30),
      Kcp = ab(0, 5),
      Kpc = ab(0, 5)
    ),
    eqn = function(){
      ktotal = Ke + Kcp # local variable, not available outside eqn block
      dX[1] = B[1] - Ka * X[1]
      dX[2] = R[1] + Ka * X[1] - (ktotal) * X[2] + Kpc * X[3]
      dX[3] = Kcp * X[2] - Kpc * X[3]
    }
)
```
:::

To specify models from the library, first ensure the correct library model name by using the `r gh_help("model_lib")` function. Also, ensure that the variable names in the library are defined somewhere in the model (mostly in the `PRI` block), and similarly that `OUT` equations are correct.

```{r}
#| label:  mod-block-lib
# this is a full example of a model that can be compiled
mod <- PM_model$new(
    pri = list( # parameter names match those listed in model_lib(), except V0
      Ka = ab(0, 5),
      Ke = ab(0, 5),
      V0 = msd(10, 3),
      Kcp = ab(0, 5),
      Kpc = ab(0, 5)
    ),
    cov = list(
        wt = interp()
    ),
    sec = function(){
        V = V0 * wt # here we define the last parameter, V, found in model_lib() for this model
    },
    eqn = function(){
        "three_comp_bolus" # name acquired from model_lib()
    }, 
    out = function(){
        Y[1] = X[2]/V # correct compartment (2) identified as the output
    },
    err = list(
        proportional(1, c(1, 0.1, 0, 0))
    )
)
```

#### File

Use the same syntax as in the R function in a list specification.

Example:

::: script
#EQN\
ktotal = Ke + KCP\
dX\[1\] = B\[1\] -KA\*X\[1\]\
dX\[2\] = R\[1\] + KA\*X\[1\] - (ktotal)\*X\[2\] + KPC\*X\[3\]\
dX\[3\] = KCP\*X\[2\] - KPC\*X\[3\]
:::

***Note***: Pmetrics no longer recognizes the old Fortran format of XP(1), which is dX\[1\], and X(1), which is X\[1\].

### Outputs {#mb-out}

The outputs in Pmetrics are the values in the `OUT` column of the data that you are trying to predict. In the model objects, each output is referred to as `Y[.]`, where "." is the number of the output and the same as the corresponding `OUTEQ` in the data. Outputs can be defined as mathematical combinations of compartment amounts (`X[.]`), primary and/or secondary parameters, or covariates.

Output equations are in R format. If they include conditional statements, these should follow R's `if(){}` convention.
<!-- 
#### App

Coming soon. -->

<!-- ![](Images/model_builder/output.png)

On this tab, you can choose the number of output equations and define the equation for each output. You do not need to specify "Y\[.\] =" in the equations. Simply include the right side of the equation in the field.

Include the **Assay Error** polynomial coefficients for each output equation. The four coefficients estimate SD according to the folowing equation: $C0 + C1 * [obs] + C2 * [obs]^2 + C3 * [obs]^3$ and $[obs]$ is the observation. The values for the coefficients should ideally come from the analytic lab in the form of inter-run standard deviations or coefficients of variation at standard concentrations. You can use the Pmetrics function `makeErrorPoly()` to choose the best set of coefficients that fit the data from the laboratory. Alternatively, you can use a generic set of coefficients. We recommend that as a start, $C0$ be set to half of the lowest concentration in the dataset and $C1$ be set to 0.1. $C2$ and $C3$ can be 0.

If you check the *Use Always?* box, the model coefficients will be used regardless of what is in the data file. Leaving it unchecked, which is the default, means that coefficients in the data file will be used if present, and the ones here in the model definition will be used only if none are available for an observation in the data.

Choose your **Model Error** model from the drop down. The additive $\lambda$ and proportional $\gamma$ models both have fixed options. If not fixed, the parameter will be estimated, starting from the number in the *Value* field. If fixed, it will be held constant at the number in the *Value* field.

The proportional model weights each observation by $1/(\gamma * SD)^2$, where $\gamma$ is either fixed or estimated. The additive model weights each observation by $1/(\lambda + SD)^2$, where $\lambda$ is either fixed or estimated. The combination model uses $1/((\gamma * SD)^2 + (\lambda + SD)^2)$.

In the proportional model, $\gamma$ is a scalar on assay SD. In general, well-designed and executed studies and models with low mis-specification will have data with $\gamma$ values approaching 1. Values \<1 suggest over-inflated assay noise. Poor quality, noisy data will result in $\gamma$ of 5 or more. A good starting value for $\gamma$ is usually 5, and sometimes 10 early in model development if data are particularly complex or noisy.

In the additive model, $\lambda$ is additive to assay SD. In general, well-designed and executed studies and models with low mis-specification will have data with $\lambda$ values approaching 0. Values of 0 may suggest over inflated assay noise. Poor quality, noisy data will result in $\lambda$ of $5*C0$ or more. A good starting value for $\lambda$ is usually $3*C0$. Note, that $C0$ should generally not be 0, as it represents machine noise (e.g. HPLC or mass spectrometer) that is always present. -->

#### List

The `out` element of the model list is function. Any R code is valid within the function. Secondary variables declared here will only be available within the function. The output equations are of the form "Y[.] = expression", where "." is the output equation number in the data. Primary and secondary variables and covariates may be used in the expression, as can conditional statements in R code.

:::{code-copy='false'}
```{r}
#| label:  mod-block-out
# this is only a partial model snippet for illustration; do not execute
out = function(){
  Y[1] = X[1]/V
}
  
```
:::

#### File

Outputs are of the form Y\[.\] = expression, where "." is the output equation number. Primary and secondary variables and covariates may be used in the expression, as can conditional statements in Fortran code. An "&" continuation prefix is not necessary in this block for any statement, although if present, will be ignored. **There can be a maximum of 6 outputs.** They are referred to as Y\[1\], Y\[2\], etc. These equations may also define a model explicitly as a function of primary and secondary variables and covariates.


Examples:

::: script
#Out\
Y\[1\] = X\[2\]/V
:::


Here's an example of an explicit model defined in the output and not the equation block.

::: script
#OUT\

RES = BOLUS\[1\] \* KIN/(KIN-KOUT) \* (EXP(-KOUT\*TPD)-EXP(-KIN\*TPD))\
Y\[1\] = RES/VD
:::

### Error {#mb-err}

All models predict observations with error. Part of this error is due to measurement or assay error. The assay error in Pmetrics is explicitly accounted for, and it may be inflated by an additive model error component $\lambda$ or a proportional component $\gamma$ that reflect additional noise and model misspecification. The remaining error is random and not optimized. 

When fitting, the error in a prediction determines how closely the model must match the observation. Pmetrics weights predictions by $1/error^2$. The standard deviation of an prediction is calculated from the assay error model as $SD_{pred} = C0 + C1*pred + C2*pred^2 + C3*pred^3$. Error is then calculated according to the selected model error structure: 

-   **Proportional**: $error = \gamma * SD_{pred}$
-   **Additive**: $error = \lambda + SD_{pred}$

**Note:** Each output equation must have an error model. In Pmetrics versions prior to 3.0.0, each output equation had its own assay error polynomial, but there was only one model error for all outputs. In Pmetrics ≥ 3.0.0, each output equation has its own assay error polynomial and its own model error, allowing for greater flexibility in modeling complex data.

<!-- #### App

Coming soon. -->

#### List

The error list element is a list of creator functions. The two possible creator functions are `proportional` or `additive`. Each of these functions has the same arguments: `initial`, `coeff`, and `fixed`.

* The `initial` argument is the starting value for $\gamma$ in the `proportional()` creator function or $\lambda$ for the `additive()` creator function. 
* For both creators, `coeff` is the vector of 4 numbers that defines a polynomial equation to permit calculation of the standard deviation of an observation, based on the noise of the assay (see `r gh_help("makeErrorPoly")`). 
* The final argument `fixed` is `FALSE` by default, which means that $\gamma$ or $\lambda$ will be estimated for each output equation as measures of additional noise in the data and/or model misspecification. If `fixed` is set to `TRUE`, then the value in `initial` will be held constant during estimation.

In the code snippet below, there are two output equations so there are two error models. The first is a $\gamma$ proportional model with starting value 1, which will be updated based on the model fit of the data, and assay polynomial coefficients of 0.15, 0.1, 0, and 0.  The second is an additive $\lambda$ model with starting value 0.1, assay polynomial coefficients of 0.1, 0, 0, and 0, and because `fixed = TRUE`, it will not be estimated and will be held constant at 0.1.

:::{code-copy='false'}
```{r}
#| label:  mod-block-err
# this is only a partial model snippet for illustration; do not execute
mod <- PM_model$new(
  out = function(){
    Y[1] = X[1]/V1
    Y[2] = X[2]/V2
  },
  err = list(
    proportional(1, c(0.15, 0.1, 0, 0)),
    additive(0.1, c(0.1, 0, 0, 0), fixed = TRUE)
  )
)
```
:::

#### File

Unlike the R6 model builder app or model list, the error block in the file is separated from the output block. This is to maintain backward compatibility with Legacy Pmetrics.

To specify the model in this block, each output requires two lines. The first line is either L=number or G=number for a $\lambda$ or $\gamma$ error model. The number is the starting value for $\lambda$ or $\gamma$. If you include an exclamation point (!) in the declaration, then $\lambda$ or $\gamma$ will be fixed and not estimated. 

The second line contains the values for the assay error polynomial coefficients: $C0$, $C1$, $C2$, and $C3$, separated by commas. `r pmetrics()` will use these values for the coefficients unless they are included in the data. 

Example 1: estimated $\lambda$, starting at 0.4, one output, use 0.1,0.1,0,0 for the coefficients unless present in the data.

::: script
#Err\
L=0.4\
0.1,0.1,0,0\
:::

Example 2: two outputs with fixed $\gamma$ of 2 for output 1, and $\lambda$ with starting value of 0.1 but estimated for output 2. Coefficients of 0.1,0.1,0,0 for the first output, and  0.3, 0.1, 0, 0 for output 2.

::: script
#Err\
G=2!\
0.1,0.1,0,0\
L=0.1\
0.3,0.1,0,0\
:::


<!-- ### Extra {#extra}

#### App, List

This block is not currently implemented in either the app or list forms.

#### File

This block will be for advanced Rust programmers only who wish to include additional code in the model. It is not yet implemented. -->

## Complete Model Examples

### List Examples {#completeEx}

Here are two complete examples of model definitions using the list format.
```{r}
#| eval: true
#| label:  mod-block-ex1
# this example of a one-compartment IV model will compile 
mod1 <- PM_model$new(
  pri = list(
    Ke = ab(0, 5),
    V = msd(100, 30)
  ),
  eqn = function(){
    one_comp_iv # in model_lib()
  },
  out = function(){
    Y[1] = X[1]/V
  },
  err = list(
    proportional(5, c(0.05, 0.1, 0, 0))
  )
)
```

When you execute the above code in your own script or Quarto document, you will see a message indicating that the model has compiled.

*Notes*:

* The #Pri block contains the primary parameters and their initial ranges, which are the only parameters that will be estimated.
  - Ke is the elimination rate constant defined by the creator function `ab()` with lower and upper bounds of 0 and 5, respectively.
  - V is the volume of distribution defined by the creator function `msd()` with a
* The value of "one_comp_iv" in the #Eqn block tells Pmetrics to use the algebraic solution for the model in the library. 
* The #Out block contains the output equation. The relevant compartment number that contains the amount of drug, 1, is obtained by looking at `model_lib()` to see which compartment is the central, so that it can be divided by `V` to obtain concentration. 
* The #Err block contains the model error structure (proportional) with a starting $\gamma$ value of 5, and assay error coefficients of 0.05, 0.1, 0, and 0.
* The comments in the code above will be ignored.

You can see a summary of the model by printing the model object, and you can plot the model structure by calling the `plot()` method.

```{r}
#| eval: true
#| label: example-mod1_print
mod1 # prints the model object
```

```{r}
#| eval: true
#| message: false
#| label: example-mod1-plot

mod1$plot() #plots the model structure
```
      

Here's an example for a parent and metabolite showing differential equations and multiple outputs . 

```{r}
#| label:  mod-block-ex2
# this example of a parent and metabolite model will compile
mod2_ode <- PM_model$new(
  pri = list(
    Ka = ab(0, 5),
    Ke = ab(0, 5),
    V = msd(100, 30),
    fm = ab(0, 1),
    Ke_m = ab(0, 5)
  ),
  eqn = function(){
    # same as two_comp_bolus in model_lib()
    dX[1] = B[1] - Ka * X[1] # bolus input here with B[1]
    dX[2] = R[1] + Ka * X[1] - Ke * X[2] # can also accept infusion input here with R[1]
    dX[3] = fm * X[2] - Ke_m * X[3] # a metabolite compartment
  },
  out = function(){
    Vm = V * 0.5 # fix metabolite volume at 50% of parent
    Y[1] = X[1]/V
    Y[2] = X[3]/Vm
  },
  err = list( # one error model for each output
    proportional(5, c(0.05, 0.1, 0, 0)),
    additive(0.1, c(0.1, 0, 0, 0))
  )
)
```


<!-- ```{r}
#| eval: true
#| message: false
#| label: example-mod-bateman

mod_bateman <- PM_model$new(
  pri = list(
    kin = ab(0, 5),
    kout = ab(0, 5),
    tpd = ab(0, 5),
    V = msd(100, 30)
  ),
  eqn = function(){
    X[1] = (X[1] + b[1]) * kin/(kin-kout) * (exp(-kout*tpd)-exp(-kin*tpd))
  },
  out = function(){
    Y[1] = X[1]/V
  },
  err = list(
    additive(0.4, c(0.3, 0.15, 0, 0))
  )
)
```

This last example is known as the Bateman equation for a model with linear absorption (KIN) into and elimination (KOUT) from a central compartment, and a time post-dose (TPD) or lag time. Here B[1] is the oral bolus dosing vector for drug 1, and V is the volume of the central compartment. -->

### File Example {#completeFileEx}


Here is a complete example of a model file, as of Pmetrics version 3.0 and higher:

::: script
#Pri\
KE, 0, 5\
V0, 0.1, 100\
KA, 0, 5\
Lag1, 0, 3\

#Cov\
wt\
C this weight is in kg\

#Sec\
V = V0\*wt\

#Lag\
LAG\[1\] = Tlag1

#Eqn\
two_comp_bolus\

#Out\
Y\[1\] = X\[2\]/V\

#Err\
L=0.4\
0.1, 0.1, 0, 0\
:::

*Notes*:

* The #Pri block contains the primary parameters and their initial ranges, which are the only parameters that will be estimated. 
* The #Cov block contains the covariate names. 
* The #Sec block contains the global secondary equations used to incorporate the covariate "wt" in the model definition. 
* The #Lag block contains the lag time for the bolus input. 
* The value of "two_comp_bolus" in the #Eqn block tells Pmetrics to use the algebraic solution for the model. 
* The #Out block contains the output equation. The relevant compartment number that contains the amount of drug, 2, is obtained by looking at `model_lib()` to see which compartment is the central, so that it can be divided by `V` to obtain concentration. 
* The #Err block contains the model error structure (additive) with a starting $\lambda$ value of 0.4, and assay error coefficients of 0.1, 0.1, 0, and 0.
* The comment line "C this weight is in kg" will be ignored.


Create and save such a file in any text editor, giving it a name ending in ".txt".

`PM_model$new()` also accepts the path to a model file to create the same model using the file. **Remember** that you should phase model files out of your workflow in favor of direct definition in R. 

```{r}
#| echo: false
#| eval: true
mod1b <- PM_model$new("Data/src/model.txt")
mod1b # look at it
```

```{r}
#| label: example-mod1b-def
# if not using the Rscript/Learn.R template created by PM_tutorial(), 
# modify the path as needed
mod1b <- PM_model$new("src/model.txt")
mod1b # look at it
```


<!-- If you build the same algebraic model in the app by selecting it from the library, and save the file, Pmetrics will write the #EQN blocks for you. -->



<!-- ## Saving and copying models from the app

![](Images/model_builder/output2.png)

The Model List tab in the bottom window of the app shows the List form of the model. You can see it develop as you complete components. Once the model is finished, you can use one or both of the buttons at the bottom to copy or save the model.

Choosing *Copy* will put the code in you clipboard and you can paste it into R. It will include `PM_model$new(...)`, so you should add code such as `mod1 <-` to your R script, and then paste. This will entirely define `mod1` as a `PM_model` equivalent to what you built in the app.

Choosing *Save* in the app will write a file named "model.txt" to your current working directory in R that contains all the code necessary to define the same model as you build in the app. Be aware that if there is another file named "model.txt" in your working directory, it will be overwritten.

While incomplete models can be copied, a model can only be saved when it is complete. It is necessary for you to at least open relevant components to ensure that all fields are finalized in the Model List before you can save it. -->
## Reserved Names {#reserved}

The following cannot be used as primary, covariate, or secondary variable names. They can be used in equations, however.

```{r}
#| echo: false
#| eval: true
#| label: example-mod1-reserved
custom_table("reserved.csv")
```

## Updating models {#update-model}

<!-- ### App

You can update an existing model in one of two ways in the app:

-   Pass the `PM_model` object to the app as an argument to `build_model()`, for example, `build_model(NPex$model)`. Then, edit the appropriate fields in the app. Save the model as a new model and/or copy it into R. Note that the app saves all models as "model.txt" so be sure not to overwrite your old model if you wish to keep it.
-   Load a previously saved model file using the dialog button in the app and proceed as for the first option above. -->

### List

To update a model, simply copy the text from the model list in R, and modify the relevant fields. Then reassign the modified model to a new object.

```{r}
#| label:  mod-block-update-example
# these examples will compile and can be executed 

mod1 <- PM_model$new(
  pri = list(
    Ka = ab(0.1, 0.9),
    Ke = ab(0.001, 0.1),
    V = ab(30, 120),
    lag1 = ab(0, 4)
  ),
  cov = list(
    wt = interp(),
    africa = interp("none"),
    age = interp(),
    gender = interp("none"),
    height = interp()
  ),
  eqn = function() {
    two_comp_bolus
  },
  lag = function() {
    lag[1] <- lag1
  },
  out = function() {
    Y[1] <- X[2] / V
  },
  err = list(
    proportional(5, c(0.02, 0.05, -0.0002, 0))
  )
)

mod2 <- PM_model$new(
  pri = list(
    Ka = ab(0.1, 0.9),
    Ke = ab(0.001, 0.5), # changed range
    V0 = ab(30, 120), # note the new name
    lag1 = ab(0, 4)
  ),
  cov = list(
    wt = interp(),
    africa = interp("none"),
    age = interp(),
    gender = interp("none"),
    height = interp()
  ),
  sec = function() { # new secondary block
    V = V0 * wt # the library model expects a parameter named "V" defined somewhere
  },
  eqn = function() {
    two_comp_bolus
  },
  lag = function() {
    lag[1] <- lag1
  },
  out = function() {
    Y[1] <- X[2] / V
  },
  err = list(
    proportional(5, c(0.02, 0.05, -0.0002, 0))
  )
)
```

This example changes the range of *Ke*, deletes *V*, adds *V0*, and defines a new secondary relationship for *V = V0 \* wt*.

### File

The only way to update a model file is to edit the text file and load it with `PM_model$new()`. The disadvantage is that you have to keep editing, copying and pasting files, and unless your documentation in R or elsewhere is good, you have no idea how you changed the model from run to run. In contrast, by using the direct model definition in R, you can see a written record in your script of how the model evolved over the project.