.. _ddmore0061_conv_example:

DDMoRe0061 Conversion Example
####################################

Overview
==========

A real world model, based on a gentamicin |pk| study [Harling2015]_. 
It has a combined |bolus| + |infusion| dose and the original |pk| model uses 
ADVAN3 TRANS4 in |nonmem|. There are some covariate effects on the clearance 
and volume of distribution of the central compartment parameters. 
There are 210 individuals and 1949 data rows. See :ddmore_link:`0061`

Here we describe the manual conversion of the |nonmem| script file 
``c:\PoPy\validation\ddmore0061\Executable_gentamicin_pk.mod``
to create the |popy| :ref:`fit_script` 
``c:\PoPy\validation\ddmore0061\popy\fit_script.pyml``.

Data Source
------------

The |nonmem| script loads the data as follows:-

.. code-block:: none

    $DATA Simulated_gentamicin_pk.csv IGNORE=@
       
Whereas the |popy| script uses this syntax:-

.. code-block:: pyml

    FILE_PATHS: { input_data_file: Simulated_gentamicin_pk.csv }

Data Preprocessing
---------------------

Because PoPy expects data in a slightly different format to |nonmem|, we make 
some modifications to the incoming data using the PREPROCESS block of the 
input script.

The critical columns for the original |nonmem| data set looks something like
:numref:`table_nm_fields_ddmore0061`.

.. _table_nm_fields_ddmore0061:

.. csv-table:: Main data fields for |nonmem| (first six rows)
    :file: ddmore0061_nm_fields_first_6_rows.csv

and are missing a ``TYPE`` column along with a flag column that serves
essentially the same purpose as Nonmem's ``MDV`` column.

We add these missing columns using the :ref:`preprocess` section of the fit 
to modify every row of the incoming data file as it is read:

.. code-block:: pyml

    PREPROCESS: |
        # Call this helper function that (a) converts Nonmem's EVID column 
        # values into PoPy's TYPE column values and (b) adds a DV_FLAG 
        # column that is 1 for observation rows (EVID=0) and 0 everywhere else.
        nonmem2popy()
        
        # Augment PoPy's new TYPE column value for dose rows (where AMT>0) 
        # with information about whether the dose is an infusion (RATE>0)
        # or a bolus (RATE=0).
        if c[AMT] > 0:
            if c[RATE] > 0:
                c[TYPE] += ":DOSE_infrate"
            else:
                c[TYPE] += ":DOSE_bolus"

        # Make a copy of Nonmem's DV column to give it a more descriptive
        # name, DRUG_CONC, and update the corresponding flag.
        # (This is optional and could be left out after changing the 
        # PREDICTIONS block to use DV rather than DRUG_CONC.)
        c[DRUG_CONC] = c[DV]
        c[DRUG_CONC_FLAG] = c[DV_FLAG]

Technically the |nonmem| data set above has no CMT field, but CMT defaults 
to 1 in |nonmem| if it is not provided. |popy| does not assign a compartment
in the data file, preferring to assign the compartment in the model definition
(where compartments have meaning).

The (critical columns of the) preprocessed dataset used by |popy| will now 
look like that in :numref:`table_popy_fields_ddmore0061`.

.. _table_popy_fields_ddmore0061:

.. csv-table:: Extra data fields for |popy| (first six rows)
    :file: ddmore0061_popy_fields_first_6_rows.csv

Here the new ``DRUG_CONC`` field will be the target value for the new 
|popy| script. The ``DRUG_CONC_FLAG`` switches the observations on/off 
appropriately. Note that ``DRUG_CONC_FLAG`` is the inverse of |nonmem|'s'
``MDV`` flag. The new ``TYPE`` field is formatted to indicate either 
|bolus| or |infusion| doses in |popy|.

Fixed/Random Effects Specification
-----------------------------------

In the |nonmem| script, the theta/omega/sigma |fe| definitions are:-

.. code-block:: none

    $THETA  
        (0,4.1) ; CL
        (1,8.63) ; V1
        (0,1.21) ; Q
        (0,8.12) ; V2
        
        (-0.01,0.00982,0.016) ; BCLC
        (-0.409) ; ALB
        (0,0.00683) ; DCLC
        
    $OMEGA  
        0.0465 ; CL
        0.376 ; Q
        
    $SIGMA  
        0.0297 ; E1
        0.0486 ; E2
        

In |popy| these |fe| are defined in |effects|:-

.. code-block:: pyml

    EFFECTS:
        POP: |
            f[CL] ~ P 4.1
            f[V1] ~ unif(1, +inf) 8.63
            f[Q] ~ P 1.21
            f[V2] ~ P 8.12
            f[BCLC_eff] ~ unif(-0.01, 0.016) 0.00982
            f[ALB_eff] = -0.409
            f[DCLC_eff] ~ P 0.00683
            
            f[CL_isv, Q_isv] ~ diag_matrix() [ [ 0.0465, 0.376 ] ]
            
            f[PNOISE_var] ~ P 0.0297
            f[ANOISE_var] ~ P 0.0486
                
Note the ranges and initial values of the |fx| variables are defined using 
different syntax from the $THETA variables. 
Note also that the |fx| have explicit names, whereas the |nonmem| list of 
thetas has been annotated with comments to make it human readable.

In |popy| you can specify a :ref:`diagonal matrix<pos_diag_mat_dist>`, which 
can be explicitly used as a covariance matrix input to a :ref:`mnorm_dist` 
distributed |rx| vector, as follows:-

.. code-block:: pyml

    EFFECTS:
        ID: |
            r[CL_isv, Q_isv] ~ mnorm([0,0], f[CL_isv, Q_isv])
            
There :pyml:`r[CL_isv]` is the equivalent of 'ETA(1)' and :pyml:`r[Q_isv]` 
is the equivalent of 'ETA(2)'. In |nonmem| you just have to remember the 
index values, whereas |popy| uses actual real life variable names.

In |popy| the ``SIGMA``s are treated as |fes| and also defined in |effects|.
                 
PK Model Parameters Specification
----------------------------------

The |nonmem| PK section is as follows:-
       
.. code-block:: none
  
    $PK       
        BSA = 71.84*(WT**0.425)*(HT**0.725)
        BSA = BSA/10000
        HT1 = HT/100
        BMI = WT/HT1**2
        
        IF(NEWIND.NE.2) BCLC = CLC
        DCLC = CLC-BCLC
        IF(NEWIND.NE.2) BBSA = BSA
        
        TVCL = THETA(1) * (1 + THETA(5)*(BCLC-81) + THETA(7)*DCLC)
        TVV1 = THETA(2) * BBSA*(ALB/34)**THETA(6)
        TVQ  = THETA(3)
        TVV2 = THETA(4)
        
        CL = TVCL * EXP(ETA(1))
        V1 = TVV1
        Q  = TVQ  * EXP(ETA(2))
        V2 = TVV2
        
        S1 = V1
    
This maps almost exactly to |popy|'s' |model_params| section:-

.. code-block:: pyml

    MODEL_PARAMS: |
        BSA = 71.84 * (c[WT]**0.425) * (c[HT]**0.725)
        BSA = BSA / 10000
        HT1 = c[HT] / 100
        BMI = c[WT] / HT1**2
   
        BCLC = c[CLC]
        DCLC = c[CLC] - BCLC
        BBSA = BSA
        
        TVCL = f[CL] * (1 + f[BCLC_eff]*(BCLC-81) + f[DCLC_eff]*DCLC)
        TVV1 = f[V1] * BBSA * (c[ALB]/34)**f[ALB_eff]
        TVQ  = f[Q]
        TVV2 = f[V2]
        
        m[CL] = TVCL * EXP(r[CL_isv])
        m[V1] = TVV1
        m[Q]  = TVQ  * EXP(r[Q_isv])
        m[V2] = TVV2
        
        m[PNOISE_var] = f[PNOISE_var]
        m[ANOISE_var] = f[ANOISE_var]
    
Note in the above |nonmem| specifies 'S1 = V1', to scale the amounts in 
compartment one. This is not necessary in |popy|, as the compartment amounts 
are converted to concentrations in the |predictions| section explicitly.

Also note that |popy| requires the lines:-

.. code-block:: pyml

    m[PNOISE_var] = f[PNOISE_var]
    m[ANOISE_var] = f[ANOISE_var]
    
To propagate the |fx| noise |fes| which are similar to the |nonmem| builtin 
sigma variables. 
   
Note in this example the |nonmem| 'NEWIND' keyword is superfluous, as the 
:pyml:`c[CLC]` covariate and BSA variable are constant within each individual.

ODEs Specification
--------------------

The |nonmem| control file prescribes the compartment model on this line:-

.. code-block:: none

    $SUBROUTINE ADVAN3 TRANS4
    
That uses the inbuilt 'ADVAN3' model with 'TRANS4' parametrisations. 
This uses the analytic solution for a two compartment model that expects 
the CL, V1, Q, V2 parameters to be defined in the |nonmem| PK section.
    
The |popy| control file here specifies the compartment model equations 
explicitly in the |derivatives| section:-

.. code-block:: pyml

    DERIVATIVES: |
        dose[DOSE_bolus] = @bolus{amt:c[AMT], lag:0.0}
        dose[DOSE_infrate] = @inf_rate{amt:c[AMT], lag:0.0, rate:c[RATE]}
        d[CENTRAL] = dose[DOSE_bolus] + dose[DOSE_infrate] - s[CENTRAL]*m[Q]/m[V1] + s[PERI]*m[Q]/m[V2] - s[CENTRAL]*m[CL]/m[V1]
        d[PERI]    =                                         s[CENTRAL]*m[Q]/m[V1] - s[PERI]*m[Q]/m[V2]

        
Note that |popy| has an equivalent builtin version of 'ADVAN3 TRANS4' as follows:-

.. code-block:: pyml

    DERIVATIVES: |
        s[CENTRAL, PERI] = @iv_two_cmp_cl{ dose: @bolus{amt:c[AMT]}, CL: m[CL], V1: m[V1], Q: m[Q], V2: m[V2] }
       
However this function is unable to process two different types of dose 
(bolus + infusion), unlike the long hand version above. Note that |popy| 
will be using an |ode| solver here (instead of a built in analytic solution), 
so we also need to specify this |popy| field:-

.. code-block:: pyml

    ODE_SOLVER: {CPPLSODA: {}}
    
Which tells |popy| to use the c++ LSODA solver, the same method as 
|nonmem| Advan13.
       
Likelihood Error Specification
--------------------------------

The |nonmem| ERROR section is as follows:-
       
.. code-block:: none

    $ERROR
        Y = F*EXP(ERR(1)) + ERR(2)
        IPRED=F
        IRES=DV-IPRED

The |popy| equivalent is shown below in the |predictions| section:-
        
.. code-block:: pyml

    PREDICTIONS: |
        p[DV_CENTRAL] = s[CENTRAL]/m[V1]
        var = m[PNOISE_var] * p[DV_CENTRAL]**2 + m[ANOISE_var]
        c[DRUG_CONC] ~ norm(p[DV_CENTRAL], var)
        
Note here the line 'Y = F*EXP(ERR(1)) + ERR(2)' defines an proportional 
exponential normal noise model with an additional additive normal noise 
component (even though it is not clear - to the |popy| developers, at least - 
what the distribution of a log normal + a normal distribution is). 
|popy| requires a known distribution to calculate a likelihood. 
If you approximate the exponential (assuming ERR(1) is small) and get 
'Y = F*(1 + ERR(1)) + ERR(2)' then this is the standard proportional + 
additive noise model. In this case, the |popy| :ref:`fit_script` return a 
very similar objective function to the |ddmore| |nonmem| objective function, 
indicating that |nonmem| also makes the same approximation.

Parameter Fitting Methods Specification
------------------------------------------
              
The |nonmem| control file specifies using FOCE fitting (METHOD=1) below:-
        
.. code-block:: none
        
   $ESTIMATION MAXEVALS=0 SIG=3 PRINT=10 METHOD=1 INTER
   
The |popy| equivalent, specifying the |nd| fitting method is:-

.. code-block:: pyml

    FIT_METHODS: [ND: {max_n_main_iterations: 0}]
      
Setting 'max_n_main_iterations' to zero, means that |popy| will optimise the |rx| parameters and leave the |fx| unchanged.
