The Pharmpy model#

At the heart of Pharmpy lies its non-linear mixed effects model abstraction. A model needs to follow the interface of the base model class pharmpy.model.Model. This means that any model format can in theory be supported by Pharmpy via subclasses that implement the same base interface. This makes any operation performed on a model be the same regardless of the underlying implementation of the model and it is one of the core design principles of Pharmpy. This allows the functions in pharmpy.modeling and tools in pharmpy.tools to be independent on the estimation software: it is only when writing and fitting a model that the implementation is estimation software specific.

Warning

This section is intended to give an overview of the Pharmpy model and its different components. For higher level manipulations, please check out the functions in pharmpy.modeling, which are described in further detail in the modeling user guide

Reading in a model#

Reading a model from a model specification file into a Pharmpy model object can be done by using the pharmpy.modeling.read_model() function.

from pharmpy.modeling import read_model

model = read_model(path / "pheno_real.mod")

model <- read_model(path / "pheno_real.mod")

Internally this will trigger a model type detection to select which model implementation to use, i.e. if it is an NM-TRAN control stream the Pharmpy NONMEM model subclass will be selected.

Writing a model#

A model object can be written to a file using its source format. By default the model file will be created in the current working directory using the name of the model.

from pharmpy.modeling import write_model
write_model(model)

write_model(model)

Optionally a path can be specified:

write_model(model, path='/home/user/mymodel.mod')

write_model(model, path='/home/user/mymodel.mod')

Inspecting the model attributes#

Name and description#

A model has a name property. After reading a model from a file the name is set to the filename without extension.

model.name

model$name

'pheno_real'

If the model has a description it can also be retrieved. For NONMEM models this is taken from $PROBLEM.

model.description

model$description

'PHENOBARB SIMPLE MODEL'

Parameters#

Model population parameters are scalar values that are used in the mathematical definition of the model and are estimated when a model is fit to a dataset. The parameters of a model are thus optimization parameters and can in turn be used as parameters in statistical distributions or as structural parameters. A parameter is represented by the pharmpy.model.Parameter class.

All parameters in a model are organized by the the pharmpy.model.Parameters class as an ordered set of pharmpy.model.Parameter. All parameters of a model can be accessed by using the parameters attribute:

params = model.parameters
params

params <- model$parameters
params

	value	lower	upper	fix
PTVCL	0.004693	0.00	∞	False
PTVV	1.009160	0.00	∞	False
THETA_3	0.100000	-0.99	∞	False
IVCL	0.030963	0.00	∞	False
IVV	0.031128	0.00	∞	False
SIGMA_1_1	0.013241	0.00	∞	False

Operations on multiple parameters are made easier using methods or properties on parameter sets. For example, to get all initial estimates as a dictionary:

params.inits

params$inits

{'PTVCL': 0.00469307,
 'PTVV': 1.00916,
 'THETA_3': 0.1,
 'IVCL': 0.0309626,
 'IVV': 0.031128,
 'SIGMA_1_1': 0.013241}

Each parameter can be retrieved using indexing:

par = params['PTVCL']

par <- params['PTVCL']

A model parameter must have a name and an initial value and can optionally be constrained to a lower and or upper bound. A parameter can also be fixed meaning that it will be set to its initial value. The parameter attributes can be read out via properties.

par.lower

par$lower

0.0

Random variables#

The random variables of a model are available through the random_variables property and are organized using the pharmpy.model.RandomVariables which is an ordered set of distributions of either pharmpy.model.NormalDistribution or pharmpy.model.JointNormalDistribution class. All random variables of a model can be accessed by using the random variables attribute:

rvs = model.random_variables
rvs

rvs <- model$random_variables
rvs

\[\begin{split}\begin{align*} \eta_{1} & \sim \mathcal{N} \left(0,IVCL\right) \\ \eta_{2} & \sim \mathcal{N} \left(0,IVV\right) \\ EPS_{1} & \sim \mathcal{N} \left(0,\sigma_{1 1}\right)\end{align*}\end{split}\]

The set of random variables can be split into subsets of random variables, for example IIVs:

rvs.iiv

rvs$iiv

\[\begin{split}\begin{align*} \eta_{1} & \sim \mathcal{N} \left(0,IVCL\right) \\ \eta_{2} & \sim \mathcal{N} \left(0,IVV\right)\end{align*}\end{split}\]

A distribution can be extracted using the name of one of the etas:

dist = rvs['ETA_1']
dist

dist <- rvs['ETA_1']
dist

\[\eta_{1}\sim \mathcal{N} \left(0,IVCL\right)\]

Similarly to parameters, we can extract different attributes from the distribution:

dist.names

dist$names

('ETA_1',)

Statements#

The model statements represent the mathematical description of the model. All statements can be retrieved via the statements property as a pharmpy.model.Statements object, which is a list of model statements of either the class pharmpy.model.Assignment or pharmpy.model.CompartmentalSystem.

statements = model.statements
statements

statements <- model$statements
statements

\begin{align*}BTIME &= \begin{cases} TIME & \text{for}\: AMT > 0 \end{cases}\\TAD &= - BTIME + TIME\\TVCL &= PTVCL \cdot WGT\\TVV &= PTVV \cdot WGT\\TVV &= \begin{cases} TVV \cdot \left(\theta_{3} + 1\right) & \text{for}\: APGR < 5 \\TVV & \text{otherwise} \end{cases}\\CL &= TVCL \cdot e^{\eta_{1}}\\V &= TVV \cdot e^{\eta_{2}}\\S_{1} &= V\\\end{align*}

Bolus(AMT, admid=1) → CENTRAL 
┌───────┐       
│CENTRAL│──CL/V→
└───────┘

\begin{align*}F &= \frac{A_{CENTRAL}{\left(t \right)}}{S_{1}}\\W &= F\\Y &= EPS_{1} \cdot W + F\\IPRED &= F\\IRES &= DV - IPRED\\IWRES &= \frac{IRES}{W}\\\end{align*}

Changing the statements of a model can be done by setting the statements property. This way of manipulating a model is quite low level and flexible but cumbersome. For higher level model manipulation use the pharmpy.modeling module.

If the model has a system of ordinary differential equations this will be part of the statements. It can easily be retrieved from the statement object

statements.ode_system

statements$ode_system

Bolus(AMT, admid=1) → CENTRAL 
┌───────┐       
│CENTRAL│──CL/V→
└───────┘

Get the amounts vector:

statements.ode_system.amounts

statements$ode_system$amounts

[A_CENTRAL(t)]

Get the compartmental matrix:

statements.ode_system.compartmental_matrix

statements$ode_system$compartmental_matrix

⎡-CL ⎤
⎢────⎥
⎣ V  ⎦

Dataset and datainfo#

See Dataset.

Execution steps#

The pharmpy.model.ExecutionSteps object contains information on how to estimate or simulate the model.

ests = model.execution_steps
ests

ests <- model$execution_steps
ests

	method	interaction	parameter_uncertainty_method	evaluation	maximum_evaluations	laplace	isample	niter	auto	keep_every_nth_iter	tool_options
0	FOCE	True	SANDWICH	False	None	False	None	None	None	None	{}

Dependent variables#

A model can describe one or more dependent variables (output variables). Each dependent variable is defined in the dependent_variables attribute. This is a dictionary of each dependent variable symbol to the corresponding DVID. If there is only one dependent variable the DVID column in the dataset is not needed and its value in this definition is unimportant. The expressions of the dependent variables are all found in the statements.

model.dependent_variables

model$dependent_variables

{Y: 1}

Low level manipulations of the model object#

Creating and adding parameters#

It is possible to create a parameter of the pharmpy.model.Parameter class.

from pharmpy.model import Parameter

par = Parameter('THETA_1', 0.1, upper=2, fix=False)
par

par <- Parameter('THETA_1', 0.1, upper=2, fix=FALSE)
par

Parameter("THETA_1", 0.1, lower=-∞, upper=2, fix=False)

Substituting symbolic random variable distribution with numeric#

from pharmpy.tools import read_modelfit_results

frem_path = path / "frem" / "pheno" / "model_4.mod"
frem_model = read_model(frem_path)
frem_model_results = read_modelfit_results(frem_path)

rvs = frem_model.random_variables
rvs

frem_path <- path / "frem" / "pheno" / "model_4.mod"
frem_model <- read_model(frem_path)
frem_model_results <- read_modelfit_results(frem_path)

rvs <- frem_model$random_variables
rvs

\[\begin{split}\begin{align*} \left[\begin{matrix}\eta_{1}\\\eta_{2}\\\eta_{3}\\\eta_{4}\end{matrix}\right] & \sim \mathcal{N} \left(\left[\begin{matrix}0\\0\\0\\0\end{matrix}\right],\left[\begin{matrix}IVCL & \omega_{2 1} & \omega_{3 1} & \omega_{4 1}\\\omega_{2 1} & IVV & \omega_{3 2} & \omega_{4 2}\\\omega_{3 1} & \omega_{3 2} & BSV_{WGT} & \omega_{4 3}\\\omega_{4 1} & \omega_{4 2} & \omega_{4 3} & BSV_{APGR}\end{matrix}\right]\right) \\ EPS_{1} & \sim \mathcal{N} \left(0,\sigma_{1 1}\right) \\ EPS_{2} & \sim \mathcal{N} \left(0,EPSCOV\right)\end{align*}\end{split}\]

Starting by extracting the variance:

omega = rvs['ETA_1'].variance
omega

omega <- rvs['ETA_1'].variance
omega

⎡  IVCL    OMEGA₂ ₁  OMEGA₃ ₁  OMEGA₄ ₁⎤
⎢                                      ⎥
⎢OMEGA₂ ₁    IVV     OMEGA₃ ₂  OMEGA₄ ₂⎥
⎢                                      ⎥
⎢OMEGA₃ ₁  OMEGA₃ ₂  BSV_WGT   OMEGA₄ ₃⎥
⎢                                      ⎥
⎣OMEGA₄ ₁  OMEGA₄ ₂  OMEGA₄ ₃  BSV_APGR⎦

Substitution of numerical values can be done directly from initial values

omega.subs(frem_model.parameters.inits)

omega$subs(frem_model$parameters$inits)

⎡ 0.0293508   0.000286193   0.0256033   -0.0676481⎤
⎢                                                 ⎥
⎢0.000286193   0.027906    -0.00161838  0.0235094 ⎥
⎢                                                 ⎥
⎢ 0.0256033   -0.00161838      1.0       0.244579 ⎥
⎢                                                 ⎥
⎣-0.0676481    0.0235094    0.244579       1.0    ⎦

or from estimated values

omega_est = omega.subs(dict(frem_model_results.parameter_estimates))
omega_est

omega_est <- omega$subs(list(frem_model_results$parameter_estimates))
omega_est

⎡0.0393451    0.0168364    0.0426897   -0.0441918⎤
⎢                                                ⎥
⎢0.0168364    0.0259368   -0.00168931  0.0654487 ⎥
⎢                                                ⎥
⎢0.0426897   -0.00168931   0.983048     0.240426 ⎥
⎢                                                ⎥
⎣-0.0441918   0.0654487    0.240426     0.983046 ⎦

Operations on this parameter matrix can be done either by using SymPy

omega_est.cholesky()

omega_est$cholesky()

⎡0.198355993103309           0                   0                  0        ⎤
⎢                                                                            ⎥
⎢0.0848797141774847  0.136865752184937           0                  0        ⎥
⎢                                                                            ⎥
⎢0.215217596061068   -0.145813819169769  0.956800771574114          0        ⎥
⎢                                                                            ⎥
⎣-0.222790344312832  0.616363695081209   0.395326548580974  0.630256438169308⎦

or in a pure numerical setting in NumPy

import numpy as np

a = np.array(omega_est).astype(np.float64)
a

a <- np$array(omega_est).astype(np$numeric64)
a

array([[ 0.0393451 ,  0.0168364 ,  0.0426897 , -0.0441918 ],
       [ 0.0168364 ,  0.0259368 , -0.00168931,  0.0654487 ],
       [ 0.0426897 , -0.00168931,  0.983048  ,  0.240426  ],
       [-0.0441918 ,  0.0654487 ,  0.240426  ,  0.983046  ]])

np.linalg.cholesky(a)

np$linalg$cholesky(a)

array([[ 0.19835599,  0.        ,  0.        ,  0.        ],
       [ 0.08487971,  0.13686575,  0.        ,  0.        ],
       [ 0.2152176 , -0.14581382,  0.95680077,  0.        ],
       [-0.22279034,  0.6163637 ,  0.39532655,  0.63025644]])