.. _ddmore0238_conv_example:
    
DDMoRe0238 Conversion Example
####################################

Overview
=============

A real world model, based on a gentamicin |pk| study [Germovsek2017]_. 

A Population PK model with IOV, infusion dosing and custom derivative 
equations to implement the PK. The dataset consists of 205 individuals 
and 2788 rows. See :ddmore_link:`0238`

Here we describe the manual conversion of the |nonmem| script file,
``c:\PoPy\validation\ddmore0238\Executable_run35b_ddm2.ctl``    
to create the |popy| :ref:`fit_script` 
``c:\PoPy\validation\ddmore0238\popy\fit_script.pyml``.

Here the |popy| 'TYPE' field is a more explicit version of the |nonmem| 
'EVID' field. The |popy| 'DV_CENTRAL' field is a copy of the |nonmem| 
'DV' field, with an additional 'DV_CENTRAL_FLAG' field to make the true 
observed values more obvious. Note here |popy| could also use the |nonmem| 
'DV' field because it is always defined for all rows with :pyml:`c[EVID]` =0
or equivalently :pyml:`c[TYPE]` ='obs'.
    
Data Source
------------

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

.. code-block:: none

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

.. code-block:: pyml
      
    FILE_PATHS:
        # path to input comma separated value file in popy data format
        input_data_file: Simulated_simdataDDM.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_ddmore0238`.

.. _table_nm_fields_ddmore0238:

.. csv-table:: Main data fields for |nonmem| (first eight rows)
    :file: ddmore0238_nm_fields_first_8_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 indicating that the dose is an infusion (RATE>0).
        if c[AMT] > 0 and c[RATE] > 0:
            c[TYPE] += ":DOSE_infrate"
                    
        # Make a copy of Nonmem's DV column to give it a more descriptive
        # name, DV_CENTRAL, and update the corresponding flag.
        # (This is optional and could be left out after changing the 
        # PREDICTIONS block to use DV rather than DV_CENTRAL.)
        c[DV_CENTRAL] = c[DV]
        c[DV_CENTRAL_FLAG] = c[DV_FLAG]

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

.. _table_popy_fields_ddmore0238:

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

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

In the |nonmem| script, the theta variable definitions are:-

.. code-block:: none

    $THETA  (0,6.20684) ; 1. TVCL (lower bound,initial estimate)
    $THETA  (0,26.5004) ; 2. TVV1  (lower bound,initial estimate)
    $THETA  (0,2.15099) ; 3. TVQ
    $THETA  (0,21.151) ; 4. TVV2
    $THETA  (0,0.270697) ; 5. TVQ2
    $THETA  (0,147.893) ; 6. TVV3
    $THETA  55.4 FIX ; 7. T50
    $THETA  3.33 FIX ; 8. Hill
    $THETA  -0.129934 ; 9. power exponent on creatinine
    $THETA  (0,1.70302) ; 10. PNA50
    
In the |popy| script the equivalent |fx| variables are:-

.. code-block:: pyml

    EFFECTS:
        POP: |
            f[TVCL] ~ P 6.20684
            f[TVV1] ~ P 26.5004
            f[TVQ] ~ P 2.15099
            f[TVV2] ~ P 21.151
            f[TVQ2] ~ P 0.270697
            f[TVV3] ~ P 147.893
            f[T50] = 55.4
            f[HILL] = 3.33
            f[CRPWR] = -0.129934
            f[P50] ~ P 1.70302

Notice that |popy| uses the natural equal sign for |fx| that are fixed and the '~' for |fx| that need to be estimated.
    
In the |nonmem| script, the omega variable definitions are:-

.. code-block:: none
    
    $OMEGA  BLOCK(2)
     0.175278  ; variance for ETA(1), initial estimate
     0.115896 0.112362  ; COvariance ETA(1)-ETA(2), var for ETA(2), initial estimate
    $OMEGA  0  FIX
    $OMEGA  0.131759
    $OMEGA  0  FIX
    $OMEGA  0.177214
    $OMEGA  BLOCK(1)
     0.0140684  ;  7. IOV_CL
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME
    $OMEGA  BLOCK(1) SAME

The ugly code above defines a 2x2 matrix of omega |fes|, followed by 4 scalar omega |fes|. The 
'$OMEGA  BLOCK(1)' + '$OMEGA  BLOCK(1) SAME' lines defines a single |iov| omega |fes| to be estimated over up to 22 occasions. |nonmem| creates an ETA normal variable for each individual corresponding to each $OMEGA above. Hence, each individual will have 4 ETAs (with 2 set to zero) and an additional 22 occasion ETAS.

You just have to do all the indexing in your head to use the ETAs in the subsequent PK block. In contrast the |popy| script defines the equivalent var |fx| variables as:-

.. code-block:: pyml
    
    EFFECTS:
        POP: |
            f[TCVL_isv, TVV1_isv] ~ spd_matrix() [
                [ 0.175278 ],
                [ 0.115896, 0.112362 ]
            ]
            f[TVQ_isv] = 0
            f[TVV2_isv] ~ P 0.131759
            f[TVQ2_isv] = 0
            f[TVV3_isv] ~ P 0.177214
            f[TVCL_iov] ~ P 0.0140684
                
Here the inter-subject variability (isv) variables are just more |fx| variables with convenient labels. Note in contrast to |nonmem| the :pyml:`f[TVCL_iov]` occasion variance is just defined once here (it being one |fe|). 

In |popy| the mechanism for creating the appropriate number of |rx| variables given the var |fx| is as shown in the 'ID' and 'OCC' sections:-

.. code-block:: pyml
    
    EFFECTS:
        ID: |
            r[TVCL_isv, TVV1] ~ mnorm([0,0], f[TCVL_isv, TVV1_isv])
            r[TVQ] ~ norm(0, f[TVQ_isv])
            r[TVV2] ~ norm(0, f[TVV2_isv])
            r[TVQ2] ~ norm(0, f[TVQ2_isv])
            r[TVV3] ~ norm(0, f[TVV3_isv])
        OCC: |
            r[TVCL_iov] ~ norm(0, f[TVCL_iov])

This syntax means that each individual has one :pyml:`r[TVCL_isv, TVV1]` and one :pyml:`r[TVQ]` instance etc. However the number of instances of :pyml:`r[TVCL_iov]` are dependent on the number of values of the 'OCC' field for each individual in the |data_file|. |ie| Each individual will have a different sample of :pyml:`r[TVCL_iov]` for each occasion. You do not need to know the highest occasion number here (unlike in |nonmem|).

In the |nonmem| script, the sigma variable definitions are:-

.. code-block:: none

    $SIGMA  0.036033  ; variance PROP res error, initial estimate
    $SIGMA  0.0164023 ; additional res error, initial estimate
    
In the |popy| script the equivalent noise |fx| variables are:-

.. code-block:: pyml
    
    EFFECTS:
        POP: |
            f[PNOISE_var] ~ P 0.036033
            f[ANOISE_var] ~ P 0.0164023
                
PK Model Parameters Specification
----------------------------------

The |nonmem| 'PK' section is defined as follows:-

.. code-block:: none

    $PK
        ; Three-comp model
        IF(NEWIND.NE.2)OTIM1=0
        IF(NEWIND.NE.2)OCOV1=0
        IF(NEWIND.NE.2)OTIM2=0
        IF(NEWIND.NE.2)OCOV2=0
        ;
        STUDY=0
        IF(ID.LT.2000) STUDY=1                     ;Glasgow, Thomson1988
        IF(ID.GE.2000.AND.ID.LT.3000) STUDY=2      ;Uppsala, Nielsen2009
        IF(ID.GE.3000) STUDY=3                     ;Estonia, unpublished
        ;
        WTKG = WT/1000
        ;
        T50  = THETA(7)
        HILL = THETA(8)
        MF   = PMA**HILL/(PMA**HILL+T50**HILL)
        ;
        CREAT2 = CREAT
        IF(CREAT.LT.0) CREAT2 = TCREA       ; when SCr is NA=-99, it is the typical SCr
        ;OF = (CREAT2/TCREA)**(THETA(9))
        ;
        P50  = THETA(10)
        ;PNAF = PNA/(P50+PNA)
        ;
        CRPWR = THETA(9)
        ;IOV code
        BOVC = 0
        IF(OCC.EQ.1)  BOVC = ETA(7)
        IF(OCC.EQ.2)  BOVC = ETA(8)
        IF(OCC.EQ.3)  BOVC = ETA(9)
        IF(OCC.EQ.4)  BOVC = ETA(10)
        IF(OCC.EQ.5)  BOVC = ETA(11)
        IF(OCC.EQ.6)  BOVC = ETA(12)
        IF(OCC.EQ.7)  BOVC = ETA(13)
        IF(OCC.EQ.8)  BOVC = ETA(14)
        IF(OCC.EQ.9)  BOVC = ETA(15)
        IF(OCC.EQ.10) BOVC = ETA(16)
        IF(OCC.EQ.11) BOVC = ETA(17)
        IF(OCC.EQ.12) BOVC = ETA(18)
        IF(OCC.EQ.13) BOVC = ETA(19)
        IF(OCC.EQ.14) BOVC = ETA(20)
        IF(OCC.EQ.15) BOVC = ETA(21)
        IF(OCC.EQ.16) BOVC = ETA(22)
        IF(OCC.EQ.17) BOVC = ETA(23)
        IF(OCC.EQ.18) BOVC = ETA(24)
        IF(OCC.EQ.19) BOVC = ETA(25)
        IF(OCC.EQ.20) BOVC = ETA(26)
        IF(OCC.EQ.21) BOVC = ETA(27)
        IF(OCC.EQ.22) BOVC = ETA(28)
        ;
        TVCL = THETA(1)*MF*(WTKG/70)**(0.632)   ; typical value of CL
        TVV1 = THETA(2)*(WTKG/70)                       ; typical value of V1
        TVQ  = THETA(3)*(WTKG/70)**(0.75)               ; ty. value of intercompartmental CL
        TVV2 = THETA(4)*(WTKG/70)                       ; ty. value of V2
        TVQ2 = THETA(5)*(WTKG/70)**(0.75)               ; ty value of CL3
        TVV3 = THETA(6)*(WTKG/70)                       ; ty value of V3
        ;
        CL   = TVCL*EXP(ETA(1)+BOVC)  ; individual value of CL
        V1   = TVV1*EXP(ETA(2))
        Q    = TVQ*EXP(ETA(3))
        V2   = TVV2*EXP(ETA(4))
        Q2   = TVQ2*EXP(ETA(5))
        V3   = TVV3*EXP(ETA(6))
        ;
        K    = CL/V1
        K12  = Q/V1
        K21  = Q/V2
        K13  = Q2/V1
        K31  = Q2/V3
        ;
        IF(EVID.EQ.1) TM=TIME
        IF(EVID.EQ.1) TAD=0
        IF(EVID.NE.1) TAD=TIME-TM
        ;
        SL1 = 0
        IF(TIME.GT.OTIM1) SL1 = (PNA-OCOV1)/(TIME-OTIM1)
        A_0(4) = PNA
        ;
        SL2 = 0
        IF(TIME.GT.OTIM2) SL2 = (CREAT2-OCOV2)/(TIME-OTIM2)
        A_0(5) = CREAT2
        
The equivalent in |popy| is the |model_params| section:-

.. code-block:: pyml

    MODEL_PARAMS: |
       
        WTKG = c[WT] / 1000
        MF = c[PMA]**f[HILL]/(c[PMA]**f[HILL] + f[T50]**f[HILL])
        TVCL = f[TVCL] * (WTKG/70)**(0.632) * MF
        TVV1 = f[TVV1] * (WTKG/70)
        TVQ  = f[TVQ]  * (WTKG/70)**(0.75)
        TVV2 = f[TVV2] * (WTKG/70)
        TVQ2 = f[TVQ2] * (WTKG/70)**(0.75)
        TVV3 = f[TVV3] * (WTKG/70)
        BOVC = r[TVCL_iov]
        CL = TVCL * exp(r[TVCL_isv] + BOVC)
        V1 = TVV1 * exp(r[TVV1])
        Q  = TVQ * exp(r[TVQ])
        V2 = TVV2 * exp(r[TVV2])
        Q2 = TVQ2 * exp(r[TVQ2])
        V3 = TVV3 * exp(r[TVV3])
        m[K]   = CL/V1
        m[K12] = Q/V1
        m[K21] = Q/V2
        m[K13] = Q2/V1
        m[K31] = Q2/V3
        m[CRPWR] = f[CRPWR]
        m[P50] = f[P50]   
        m[SL1] = 0
        if c[CREAT] < 0:
            m[CREAT2] = c[TCREA]
        else:
            m[CREAT2] = c[CREAT]  
        m[SL2] = 0
        m[V1] = V1
        m[PNOISE_var] = f[PNOISE_var]
        m[ANOISE_var] = f[ANOISE_var]    
        
The |popy| version is shorter, mainly because it does |not| have to write an if statement for all values of the OCC variable |ie| 'IF(OCC.EQ.1) BOVC = ETA(7)' etc.

Initial State Specification
------------------------------

Note that the |nonmem| $PK section contains the lines:-

.. code-block:: none

    A_0(4) = PNA
    A_0(5) = CREAT2

These are necessary to provide initial values for the $DES section. The equivalent in |popy| is the |states| section:-

.. code-block:: pyml

    STATES: |
        s[COVCMT1] = c[PNA]
        s[COVCMT2] = m[CREAT2] 
        
which provides |sx| values at t=0.0 for the |popy| |derivatives| section. Note by default in both |popy| |sx| =0.0, unless the |states| section says otherwise. A similar convention is used in |nonmem|.  

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

The |nonmem| $DES section is as follows:-

.. code-block:: none

    $DES
        TCOV1 = A(4)
        TCOV2 = A(5)
        PNAF = TCOV1/(P50+TCOV1)
        OF = (TCOV2/TCREA)**CRPWR

        DADT(1) = A(3)*K31+A(2)*K21-A(1)*(K*PNAF*OF+K12+K13)
        DADT(2) = A(1)*K12-A(2)*K21
        DADT(3) = A(1)*K13-A(3)*K31
        DADT(4)= SL1
        DADT(5)= SL2
        
This is very similar to the |popy| |derivatives| section:-

.. code-block:: pyml

    DERIVATIVES: |
        TCOV1 = s[COVCMT1]
        TCOV2 = s[COVCMT2]
        PNAF = TCOV1 / (m[P50]+TCOV1)
        OF = (TCOV2/c[TCREA])**m[CRPWR]
        dose[DOSE_infrate] = @inf_rate{ amt: c[AMT], rate: c[RATE], lag: 0 }
        d[CENTRAL] = dose[DOSE_infrate] + m[K31]*s[PERIPH2] + m[K21]*s[PERIPH1] - (m[K]*PNAF*OF+m[K12]+m[K13])*s[CENTRAL]
        d[PERIPH1] = - m[K21]*s[PERIPH1] + m[K12] *s[CENTRAL]
        d[PERIPH2] = - m[K31]*s[PERIPH2] + m[K13] *s[CENTRAL]
        d[COVCMT1] =   m[SL1]
        d[COVCMT2] =   m[SL2]

except that in |popy| the :pyml:`dose[DOSE_infrate]` is defined and explicitly placed in the :pyml:`s[CENTRAL]` compartment. |nonmem| puts all boluses/infusions in compartment one by convention, unless a CMT field is defined in the data set (it isn't in this case). In |nonmem| you have to examine the data file to determine the type of dose, however in the |popy| syntax above it's pretty clear it's an infusion defined by the :pyml:`c[AMT]` and :pyml:`c[RATE]` data entries in each dosing data row.

The |ode| solver parameters for |nonmem| are set here:-

.. code-block:: none
    
    $SUBROUTINE ADVAN6 TOL=6
    
Here 'ADVAN6' is the |nonmem| ode solver for non-stiff systems. With a relative tolerance of 1e-6. Note this solver is different from the CPPLSODA method employed by |popy| which is similar to selecting 'ADVAN13' in |nonmem|. In this case the |ode| solver method makes little difference, as both |popy| and |nonmem| return similar objective values for this model.
    
In |popy| the |ode| solver parameters are defined as follows:-

.. code-block:: pyml

    ODE_SOLVER:
        CPPLSODA:
            # Absolute tolerance of ode solver.
            atol: 1e-12

            # Relative tolerance of ode solver.
            rtol: 1e-12

            # Maximum number of steps allowed in ode solver.
            max_nsteps: 10000000

Likelihood Error Specification
--------------------------------

The $ERROR block in the |nonmem| script defines a proportional + additional noise model:-

.. code-block:: none

    $ERROR
        IPRED  = A(1)/V1
        Y      = IPRED*(1+EPS(1)) + EPS(2)
        
The |predictions| section in the |popy| script defines the same error model:-

.. code-block:: pyml

    PREDICTIONS: |
        p[DV_CENTRAL] = s[CENTRAL]/m[V1]
        var = p[DV_CENTRAL]**2 * m[PNOISE_var] + m[ANOISE_var]
        c[DV_CENTRAL] ~ norm(p[DV_CENTRAL], var) 
        
Parameter Fitting Methods Specification
------------------------------------------

The |nonmem| script defines the FOCE fitting method (METHOD=1) here:-

.. code-block:: none

      $ESTIMATION METHOD=1 INTER MAXEVAL=0 PRINT=1 
      
The equivalent in the |popy| script is here:-

.. code-block:: pyml

    FIT_METHODS: [ND: {max_n_main_iterations: 0}]

Note both |nonmem| and |popy| only optimize the |rx| in the fitting run as the number of main iterations is set to zero.
    
        