From NONMEM to Pumas - Bonus Case Study

Author

Patrick Kofod Mogensen

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)
$DATA CS1_IV1EST_PAR.CSV IGNORE=C
    ; 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"
$INPUT ID TIME CONC=DV AMT DOSE MDV AGE WT SCR ISM CLCR
    ; 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
    CL = THETA(1)*EXP(ETA(1))
        ;THETA is the population value of clearance, exponential-between subject variability
        ;model is used
        ;ETA is the between subject variability parameter
    V = THETA(2)*EXP(ETA(2)) S1=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
   IPRED = F
   Y=F+F*ERR(1)+ERR(2)
      ;$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

$ESTIMATION METH=0 MAXEVAL=9999 PRINT=5 POSTHOC
   ;$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
NOPRINTONEHEADERFILE=CS1_IV1ESTFPDF.fit
;$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
        desc = "INTRAVENOUS BOLUS STUDY"
        timeu = u"hr" # hour
    end
end

2.2 `$Data`

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

$DATA CS1_IV1EST_PAR.CSV IGNORE=C
    ; 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
pkdata = CSV.read("CS1_IV1EST_PAR.csv", DataFrame; header=4)
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.

$INPUT ID TIME CONC=DV AMT DOSE MDV AGE WT SCR ISM CLCR
    ; 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.

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

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.

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

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
        desc = "INTRAVENOUS BOLUS STUDY"
        timeu = u"hr" # hour
    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
    CL = THETA(1)*EXP(ETA(1))
        ;THETA is the population value of clearance, exponential-between subject variability
        ;model is used
        ;ETA is the between subject variability parameter
    V = THETA(2)*EXP(ETA(2)) S1=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
        desc = "INTRAVENOUS BOLUS STUDY"
       timeu = u"hr" # hour
    end

    @pre begin
         CL = θcl * exp(η[1])
         Vc = θvc * exp(η[2])
    end

    @dynamics Central1
end

2.6 `$ERROR`

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

$ERROR
   IPRED = F
   Y=F+F*ERR(1)+ERR(2)
      ;$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.

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

    @pre begin
         CL = θcl * exp(η[1])
         Vc = θvc * exp(η[2])
    end
    
    @dynamics Central1
    
    @derived begin
         conc_model := @. Central / Vc
         CONC ~ @. Normal(conc_model, sqrt(σ_add^2 + (conc_model*σ_prop)^2))
    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
    θcl ∈ RealDomain(lower=0.1, init=1.0)
    θvc ∈ RealDomain(lower=1.0, init=1.0, upper=20.0)
    Ω ∈ PDiagDomain(2)
    σ_add ∈ RealDomain(lower=1.0^2)
    σ_prop ∈ RealDomain(lower=0.09^2)
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 OMEGAs. 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.

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

    @param begin
         θcl ∈ RealDomain(lower=0.1, init = 1.0)
         θvc ∈ RealDomain(lower=1.0, upper=20.0, init=10.0)
         Ω ∈ PDiagDomain(2)
         σ_add ∈ RealDomain(lower=0.0, init = 1.0^2)
         σ_prop ∈ RealDomain(lower=0.0, init = 0.09^2)
    end
    @random begin
         η ~ MvNormal(Ω)
    end
    @pre begin
         CL = θcl * exp(η[1])
         Vc = θvc * exp(η[2])
    end
    @dynamics Central1
    @derived begin
         conc_model := @. Central / Vc
         CONC ~ @. Normal(conc_model, sqrt(σ_add^2 + (conc_model*σ_prop)^2))
    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.

$ESTIMATION METH=0 MAXEVAL=9999 PRINT=5 POSTHOC
   ;$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

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

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.

model_infer = infer(model_fit)
model_vcov = model_infer.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
NOPRINTONEHEADERFILE=CS1_IV1ESTFPDF.fit
;$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.

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

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.

Reuse

CC BY-SA 4.0

1 From NONMEM to Pumas

2 Case Study I

2.1 $Problem

2.2 $Data

2.3 $Input

2.4 $SUBROUTINE

2.5 $PK

2.6 $ERROR

2.7 $THETA, $OMEGA, $SIGMA

2.8 $ESTIMATION

2.9 $COVARIANCE

2.10 $TABLE