# From NONMEM to Pumas - Bonus Case Study

## 1 From NONMEM to Pumas

In Case Study I through III (found here) we re-did the analyses from https://ctm.umaryland.edu/#/ms-pharma/model/bgt/nmc as they would have to be done in Pumas. Many new users of Pumas come from a background in NONMEM. A common question is then: “How do I translate my NONMEM model into Pumas?” Hopefully, this tutorial will help with the understand of the connection between parts of a NONMEM control stream and a Pumas model script for Case Study 1.

## 2 Case Study I

```
$PROBLEM INTRAVENOUS BOLUS STUDY; Units: Time=hr, Concentration=ug/ml, Amount=mg
; Age=yrs, Weight=kg, CLCR=ml/min,
; ISM = 1(is male) or 0(is female)
=C
$DATA CS1_IV1EST_PAR.CSV IGNORE; Specify name or path of data file. Data file should be comma delimited file
; IGNORE = # allows to comment out unwanted lines from the data file by prefixing
; them with "C"
=DV AMT DOSE MDV AGE WT SCR ISM CLCR
$INPUT ID TIME CONC; list of data items to be input to NONMEM. Items has
; DV = Dependant variable (concentration in plasma)
; AMT = Dose
; MDV = Missing dependant variable (MDV = 1 when DV = to match the datafile. 0 or missing(.))
$SUBROUTINE ADVAN1 TRANS2; $SUBROUTINE specifies which model from PREDPP(prediction for population pharmacokinetics)
; library is to be used to fit the data
; ADVAN1 = Model ( One compartment iv), TRANS2 = Parameterization in cl and v terms
; (translator will convert CL and V into Elimintaion rate constants(k))
$PK;$PK block assigns thetas to fixed effect parameters (eg.CL,V)
;Specifies random effect model(eg.additive,proportional or exponential) for
;structural model parameters
= THETA(1)*EXP(ETA(1))
CL ;THETA is the population value of clearance, exponential-between subject variability
;model is used
;ETA is the between subject variability parameter
= THETA(2)*EXP(ETA(2)) S1=V
V ;Sn refers to scaling factor to respective compartment volume
;(Eg.Dose is in mg, DV (Concentration) is in ug/ml,
;V will be in litres ie., no scaling factor required)
$ERROR= F
IPRED =F+F*ERR(1)+ERR(2)
Y;$ERROR block specifies a model for (residual) Within subject variability
;Combined error model is used ie., both additive and proportional error model
;F = Modeled value of the dependant variable or individual predicted concentrations
$THETA;$THETA specifies initial estimates and bounds for structural model parameters
0.1,1) ;THETACL ;Initial estimate for clearance (lower bound, estimate, upper bound)
(1,10,20) ;THETAVOL ;Initial estimate for volume (lower bound, estimate, upper bound)
(
$OMEGA;$OMEGA specifies initial estimates of the variance of between subject variability
0.09 ;ETA(1) ;Equivalent to 30% ((sqrt of 0.09) * 100) variability for clearance
0.09 ;ETA(2) ;Equivalent to 30% variability for volume
$SIGMA;$SIGMA specifies initial estimates of the variance of (residual) Within subject variability
0.09 ;ERR1 ;Proportional error of 30%
1 ;ERR2 ;Additive error = 0.32ug/ml-Normally the LOQ of the assay method
=0 MAXEVAL=9999 PRINT=5 POSTHOC
$ESTIMATION METH;$ESTIMATION sets condition for estimation of parameters
;METHOD = 0 refers to first order estimation method (ie.,etas set to zero during
;computation of objective function)
;MAXEVAL = 9999 refers to maximum allowable number of evaluations of the objective
;function during the estimation step
;PRINT = 5 refers to iteration summaries printed every 5th iteration including zeroth
;and last iteration
;POSTHOC refers to estimation of etas for each individual, ie., individual parameter
;estimates are obtained after estimation step
$COVARIANCE;$COVARIANCE provides standard error or precision of estimates
$TABLE ID TIME DV IPRED DOSE CL V ETA1 ETA2 AGE WT ISM CLCR=CS1_IV1ESTFPDF.fit
NOPRINTONEHEADERFILE;$TABLE specifies geneartion of tables for ID TIME DV Y
;NOPRINT specifies no printed table appears in the NONMEM output
;ONEHEADER specifies only one header for every 900 segment of the table output
;(used only with file option)
;file specifies the output table is written to given filename in ASCII format.
```

In order to understand how write this in Pumas we need to break it down. We will go by the order of the control stream and in the end we will write out the code in full.

### 2.1 `$Problem`

The first part of the control stream is the problem title

```
$PROBLEM INTRAVENOUS BOLUS STUDY; Units: Time=hr, Concentration=ug/ml, Amount=mg
; Age=yrs, Weight=kg, CLCR=ml/min,
; ISM = 1(is male) or 0(is female)
```

Pumas does not have a direct `$PROBLEM`

equivalent. In NONMEM, it is used to give a title to the project including model, estimation, inference, and more. In Pumas, each analysis is much more flexible such that one script can contain several calls to inference using `infer`

(including bootstraps). There is the possibility of adding a name to a model. Inside the the `@model`

block it is possible to add a description using the `@metadata`

block. Units are mentioned as comments after the `$PROBLEM`

statement. In Pumas, these can be added as `@metadata`

and as comments in the `@param`

block that we will see later.

```
@model begin
@metadata begin
= "INTRAVENOUS BOLUS STUDY"
desc = u"hr" # hour
timeu end
end
```

### 2.2 `$Data`

The next part of the control stream is the data input.

```
=C
$DATA CS1_IV1EST_PAR.CSV IGNORE; Specify name or path of data file. Data file should be comma delimited file
; IGNORE = # allows to comment out unwanted lines from the data file by prefixing
; them with "C"
```

Here, we point to the data set and rows to ignore. The specific line tells us to ignore all rows that start with C (for Comment). In Pumas, we also have to read in the data from a source such as `xpt`

, `sas7bdat`

, or `csv`

. In Case Study I, the equivalent line is

```
using CSV
= CSV.read("CS1_IV1EST_PAR.csv", DataFrame; header=4)
pkdata first(pkdata, 10)
```

`CSV.read`

expects a source first (the path to the file) and a sink (the type of table to construct, here a `DataFrame`

). The `header`

keyword tells us to ignore the first 4 lines (those that start with C),

### 2.3 `$Input`

The `$Input`

statement in NONMEM allows you to name the columns from the file in the `$DATA`

statement.

```
=DV AMT DOSE MDV AGE WT SCR ISM CLCR
$INPUT ID TIME CONC; list of data items to be input to NONMEM. Items has
; DV = Dependant variable (concentration in plasma)
; AMT = Dose
; MDV = Missing dependant variable (MDV = 1 when DV = to match the datafile. 0 or missing(.))
```

The `CONC=DV`

statements indicates that the third column should use the name CONC and that it is the reserved `DV`

keyword. In other words, the thirds column is the dependent variable and we want all output regarding this variable to have the `CONC`

label in the output. Then comes the two dose related columns `AMT`

and `DOSE`

, the `MDV`

column that indicates when an observation should be set to missing (and ignored) and finally five covariates.

The dataset used in the tutorial was ready for NONMEM, but not exactly in the format we expect in Pumas. The tutorial has the required data steps. We also renamed the covariate names to be easier to read. After those, the equivalent in Pumas is the following.

```
= read_pumas(
population
pkdata;= :ID,
id = :TIME,
time = [:CONC],
observations = :EVID,
evid = :AMT,
amt = :CMT,
cmt = [:WEIGHT, :CLCR, :AGE, :SEX], # WT, CLCR, AGE, ISM. SCR is not used in the model so we ignore it
covariates )
```

This looks extra verbose, but this is only because we are using “NONMEM” style column names. If all the column names had been lowercase instead, we can simply write the step as follows.

`= read_pumas(pkdata; observations = [:conc], covariates = [:WEIGHT, :CLCR, :AGE, :SEX]) population `

### 2.4 `$SUBROUTINE`

The next statement in the control stream is the `$SUBROUTINE`

. This specifies the ADVAN and the parameter format.

```
$SUBROUTINE ADVAN1 TRANS2; $SUBROUTINE specifies which model from PREDPP(prediction for population pharmacokinetics)
; library is to be used to fit the data
; ADVAN1 = Model ( One compartment iv), TRANS2 = Parameterization in cl and v terms
; (translator will convert CL and V into Elimintaion rate constants(k))
```

The equivalent code in Pumas is found in the `@model`

. The `@dynamics`

block specifies the equivalent to `ADVAN1 TRANS2`

below.

```
@model begin
@metadata begin
= "INTRAVENOUS BOLUS STUDY"
desc = u"hr" # hour
timeu end
@dynamics Central1
end
```

Where `Central1`

expects `Vc`

and `CL`

to be defined in the following block.

### 2.5 `$PK`

The `$PK`

statement specifies the model parameters including random effects and covariate models in NONMEM.

```
$PK;$PK block assigns thetas to fixed effect parameters (eg.CL,V)
;Specifies random effect model(eg.additive,proportional or exponential) for
;structural model parameters
= THETA(1)*EXP(ETA(1))
CL ;THETA is the population value of clearance, exponential-between subject variability
;model is used
;ETA is the between subject variability parameter
= THETA(2)*EXP(ETA(2)) S1=V
V ;Sn refers to scaling factor to respective compartment volume
;(Eg.Dose is in mg, DV (Concentration) is in ug/ml,
;V will be in litres ie., no scaling factor required)
```

`ADVAN1 TRANS2`

understands the reserved keywords `CL`

, `V`

. There is no scale (`S1=V`

) equivalent in Pumas. Instead, you can scale the variables when needed. The other parameters have similar reserved parameter names in Pumas when using the `Central1`

model:

`CL`

in NONMEM is`CL`

in Pumas`V1`

in NONMEM is`Vc`

in Pumas (Volume of distribution for Central)

This results in the following `@pre`

block (the Pumas equivalent to `$PK`

)

```
@metadata begin
= "INTRAVENOUS BOLUS STUDY"
desc = u"hr" # hour
timeu end
@pre begin
= θcl * exp(η[1])
CL = θvc * exp(η[2])
Vc end
@dynamics Central1
end
```

### 2.6 `$ERROR`

The statistical model of the modeled concentrations are specified in the `$ERROR`

statement.

```
$ERROR= F
IPRED =F+F*ERR(1)+ERR(2)
Y;$ERROR block specifies a model for (residual) Within subject variability
;Combined error model is used ie., both additive and proportional error model
;F = Modeled value of the dependant variable or individual predicted concentrations
```

Here, `F`

specifies the prediction of the model given the parameters and solution. Since this statement is evaluated with `etas`

during estimation, `F`

is assigned to `IPRED`

here. Then, we build the error model with `ERR`

statements. In Pumas, we do not build the distribution of `Y`

like this. Rather we specify the distributional assumption directly.

```
= @model begin
iv1cmt @metadata begin
= "INTRAVENOUS BOLUS STUDY"
desc = u"hr" # hour
timeu end
@pre begin
= θcl * exp(η[1])
CL = θvc * exp(η[2])
Vc end
@dynamics Central1
@derived begin
:= @. Central / Vc
conc_model ~ @. Normal(conc_model, sqrt(σ_add^2 + (conc_model*σ_prop)^2))
CONC end
end
```

Here, `σ_add`

is the standard deviation associated with `ERR(2)`

in the NONMEM model and `σ_prop`

is the same for `ERR(1)`

.

### 2.7 `$THETA`

, `$OMEGA`

, `$SIGMA`

Finally, we get to the last part of the model code: the fixed effects and random effects specification. In Pumas, these are typically placed first instead of last. If we start with the fixed effects, these are specified as follows in the control stream.

```
$THETA;$THETA specifies initial estimates and bounds for structural model parameters
0.1,1) ;THETACL ;Initial estimate for clearance (lower bound, estimate)
(1,10,20) ;THETAVOL ;Initial estimate for volume (lower bound, estimate)
(
$OMEGA;$OMEGA specifies initial estimates of the variance of between subject variability
0.09 ;ETA(1) ;Equivalent to 30% ((sqrt of 0.09) * 100) variability for clearance
0.09 ;ETA(2) ;Equivalent to 30% variability for volume
$SIGMA;$SIGMA specifies initial estimates of the variance of (residual) Within subject variability
0.09 ;ERR1 ;Proportional error of 30%
1 ;ERR2 ;Additive error = 0.32ug/ml-Normally the LOQ of the assay method
```

The equivalent specification in Pumas is:

```
@param begin
∈ RealDomain(lower=0.1, init=1.0)
θcl ∈ RealDomain(lower=1.0, init=1.0, upper=20.0)
θvc ∈ PDiagDomain(2)
Ω ∈ RealDomain(lower=1.0^2)
σ_add ∈ RealDomain(lower=0.09^2)
σ_prop end
```

Since we specify the “SIGMA” parameters as standard deviations in this example in Pumas, we have to square the initial values. We put no restriction on names and where different parameters can be used in other blocks, but we do suggest the best practice of using `Ω`

for variance-covariance matrices for random effects, `ω`

for individual standard deviations used in scalar random effect specifications, and `σ`

for standard deviations used in error models. If you prefer to specify variances over standard deviations, we suggest using `ω²`

and `σ²`

respectively.

In NONMEM, the random effect specification was implicit based on the `OMEGA`

s. In Pumas, we have more flexibility with respect to the specification of named random effects and the distributions of them (say, a `Beta`

distributed random effect for a bioavailability). In the current case study, the random effects specification is simple.

```
@random begin
~ MvNormal(Ω)
η end
```

Then, the final model looks as follows.

```
= @model begin
iv1cmt @metadata begin
= "INTRAVENOUS BOLUS STUDY"
desc = u"hr" # hour
timeu end
@param begin
∈ RealDomain(lower=0.1, init = 1.0)
θcl ∈ RealDomain(lower=1.0, upper=20.0, init=10.0)
θvc ∈ PDiagDomain(2)
Ω ∈ RealDomain(lower=0.0, init = 1.0^2)
σ_add ∈ RealDomain(lower=0.0, init = 0.09^2)
σ_prop end
@random begin
~ MvNormal(Ω)
η end
@pre begin
= θcl * exp(η[1])
CL = θvc * exp(η[2])
Vc end
@dynamics Central1
@derived begin
:= @. Central / Vc
conc_model ~ @. Normal(conc_model, sqrt(σ_add^2 + (conc_model*σ_prop)^2))
CONC end
end
```

As you may have noticed, we write the parameter, random effects, pk parameters, dynamical system specification and error model in a different order than what was in the original control stream. This is because we find it more clear when reading the code to define things before they are used.

### 2.8 `$ESTIMATION`

Just as we had for the `$DATA`

and `$INPUT`

statements, we have a difference between the two systems when it comes to `$ESTIMATION`

. In Pumas, we do not include this information in the “model” because the workflow is typically more interactive.

```
=0 MAXEVAL=9999 PRINT=5 POSTHOC
$ESTIMATION METH;$ESTIMATION sets condition for estimation of parameters
;METHOD = 0 refers to first order estimation method (ie.,etas set to zero during
;computation of objective function)
;MAXEVAL = 9999 refers to maximum allowable number of evaluations of the objective
;function during the estimation step
```

The equivalent to the above is something like

```
= init_params(model)
param0 = fit(model, population, param0, FO(); optim_options = (iterations = 9999)) model_fit
```

Here, `METH=0`

is equivalent to specifying `FO()`

in Pumas, `MEXAEVAL`

is the same as the `iterations`

key in `optim_options`

, `PRINT`

has no direct equivalent in Pumas. `POSTHOC`

can be obtained by called the `empirical_bayes`

function later in the script. The specification in NONMEM is necessary because FO does not need to estimate the empirical bayes estimates (EBEs) during fitting. POSTHOC forces the calculation of the EBEs at the end.

### 2.9 `$COVARIANCE`

In Pumas, the `$COVARIANCE`

step is equivalent to the `infer`

function.

```
$COVARIANCE;$COVARIANCE provides standard error or precision of estimates
```

To calculate the asymptotic variance-covariance matrix, use the `infer`

function on the fit output and grab the `vcov`

field.

```
= infer(model_fit)
model_infer = model_infer.vcov model_vcov
```

### 2.10 `$TABLE`

Since each NONMEM run is invoked based on the control stream it is necessary to specify what to output when the fitting and inference steps have completed. In Pumas, it is possible interactively work with the objects and save what is needed when the users wishes.

```
$TABLE ID TIME DV IPRED DOSE CL V ETA1 ETA2 AGE WT ISM CLCR=CS1_IV1ESTFPDF.fit
NOPRINTONEHEADERFILE;$TABLE specifies geneartion of tables for ID TIME DV Y
;NOPRINT specifies no printed table appears in the NONMEM output
;ONEHEADER specifies only one header for every 900 segment of the table output
;(used only with file option)
;file specifies the output table is written to given filename in ASCII format.
```

The above information can also be saved by computing the `inspect`

quantities, convert them to a `DataFame`

and save them to a file.

```
= inspect(model_fit)
model_inspect = DataFrame(model_inspect)
model_inspect_df write("model2_inspect.csv", model_inspect_df) CSV.
```

## 3 Conclusion

In this tutorial, we saw how to translate the case study I from NONMEM to Pumas. Hopefully, this helps new users connect the dots and understand how one section of a NONMEM control stream relates to a model block or function call in Pumas.