Simple estimation example#

Here is a simple example of how to create a model and performing an estimation.

Create a model#

We start with creating a model. A model can either be imported using read_model:

from pharmpy.modeling import read_model

model = read_model('path/to/model')

model <- read_model('path/to/model')

or it can be created from a dataset with create_basic_pk_model:

from pharmpy.modeling import create_basic_pk_model

model = create_basic_pk_model('path/to/dataset')

model <- create_basic_pk_model('path/to/dataset')

Lets create a basic IV PK model:

from pharmpy.modeling import *

dataset_path = 'tests/testdata/nonmem/pheno.dta'
model = create_basic_pk_model(administration='iv',
                                    dataset_path=dataset_path,
                                    cl_init=0.01,
                                    vc_init=1.0,
                                    mat_init=0.1)

dataset_path <- 'tests/testdata/nonmem/pheno.dta'
model <- create_basic_pk_model(administration='iv',
                                    dataset_path=dataset_path,
                                    cl_init=0.01,
                                    vc_init=1.0,
                                    mat_init=0.1)

We can now examine the model code:

model.statements

model$statements

\begin{align*}CL &= POP_{CL} \cdot e^{\eta_{CL}}\\VC &= POP_{VC} \cdot e^{\eta_{VC}}\\\end{align*}

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

\begin{align*}IPRED &= \frac{A_{CENTRAL}{\left(t \right)}}{VC}\\IPREDADJ &= \begin{cases} 2.225 \cdot 10^{-16} & \text{for}\: IPRED &= 0 \\IPRED & \text{otherwise} \end{cases}\\Y &= IPRED + IPREDADJ \cdot \epsilon_{p}\\\end{align*}

When creating a model with create_basic_pk_model the model will be a pharmpy model. To convert it to a NONMEM model use:

model = convert_model(model, 'nonmem')

model <- convert_model(model, 'nonmem')

We can then examine the NONMEM model code:

print_model_code(model)

print_model_code(model)

$PROBLEM Start model
$INPUT ID TIME AMT WGT APGR DV FA1 FA2
$DATA file.csv IGNORE=@
$SUBROUTINES ADVAN1 TRANS2
$ABBR REPLACE ETA_CL=ETA(1)
$ABBR REPLACE ETA_VC=ETA(2)
$PK
CL = THETA(1)*EXP(ETA_CL)
VC = THETA(2)*EXP(ETA_VC)
V = VC
$ERROR
IPRED = A(1)/VC
IF (IPRED.EQ.0) THEN
    IPREDADJ = 2.22500000000000E-16
ELSE
    IPREDADJ = IPRED
END IF
Y = IPRED + EPS(1)*IPREDADJ
$THETA  (0,0.01) ; POP_CL
$THETA  (0,1.0) ; POP_VC
$OMEGA BLOCK(2)
0.1	; IIV_CL
0.01	; IIV_CL_IIV_VC
0.1	; IIV_VC
$SIGMA  0.09 ; sigma
$TABLE ID TIME DV CIPREDI PRED CWRES FILE=mytab ONEHEADER NOAPPEND NOPRINT
$ESTIMATION METHOD=COND INTER MAXEVAL=99999

Modify model#

Now the model can be manipulated. We can for example add a peripheral compartment:

model = add_peripheral_compartment(model)
model.statements

model <- add_peripheral_compartment(model)
model$statements

\begin{align*}VP_{1} &= POP_{VP1}\\QP_{1} &= POP_{QP1}\\CL &= POP_{CL} \cdot e^{\eta_{CL}}\\VC &= POP_{VC} \cdot e^{\eta_{VC}}\\V_{1} &= VC\\Q &= QP_{1}\\V_{2} &= VP_{1}\\\end{align*}

Bolus(AMT, admid=1) → CENTRAL 
┌───────────┐        
│PERIPHERAL1│        
└───────────┘        
  ↑       │          
Q/V1     Q/V2        
  │       ↓          
┌───────────┐        
│  CENTRAL  │──CL/V1→
└───────────┘

We can now see that a compartment has been added to the model.

We can also remove IIV from a model parameter:

model = remove_iiv(model, "CL")
model.statements.find_assignment("CL")

model <- remove_iiv(model, "CL")
model$statements$find_assignment("CL")

\[CL = POP_{CL}\]

or add IIV to a parameter:

model = add_iiv(model, "CL", "exp")
model.statements.find_assignment("CL")

model <- add_iiv(model, "CL", "exp")
model$statements$find_assignment("CL")

\[CL = POP_{CL} \cdot e^{\eta_{CL}}\]

For more information about what transformations can be applied to a model see here.

Estimate model#

When we have created our model we can perform an estimation using fit:

from pharmpy.tools import fit

res = fit(model)

res <- fit(model)

Note

In order to esimate using any of the supported softwares (NONMEM, nlmixr2, rxode2) you need to have a configuration file set up with a path to NONMEM, instructions can be found here.

Analyze the results#

We can now analyze the fit results.

Let’s look at the ofv:

res.ofv

res$ofv

684.9291087566957

and the parameter estimates:

res.parameter_estimates

res$parameter_estimates

POP_CL           0.006721
POP_VC           0.778863
POP_QP1          4.332420
POP_VP1          0.563051
IIV_CL           0.239591
IIV_CL_IIV_VC    0.314295
IIV_VC           0.451086
sigma            0.013053
Name: estimates, dtype: float64

For more information about model estimation see here.