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

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

A real world model, based on a gentamicin |pk| study [Germovsek2017]_. Population PK, 205 individuals. 2788 total data rows. Infusion dosing, IOV + custom derivative equations to implement PK. See :ddmore_link:`0238`

Data Conversion
==================

We describe the conversion of this data set in |nonmem| format:-

.. code-block:: console

    c:\PoPy\validation\ddmore0238\Simulated_simdataDDM.csv
    
To this data set in |popy| format:-

.. code-block:: console

    c:\PoPy\validation\ddmore0238\popy\output_popy_data.csv

Using the :ref:`n2pdat_script` located here:-

.. code-block:: console

    c:\PoPy\validation\ddmore0238\popy\n2p_script.pyml

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 8 rows) 
    :file: ddmore0238_nm_fields_first_8_rows.csv
    
The final data set looks something like :numref:`table_popy_fields_ddmore0238`.

.. _table_popy_fields_ddmore0238:

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

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'.
    
Script Conversion
====================

We describe the manual conversion of the |nonmem| script file:-

.. code-block:: console

    c:\PoPy\validation\ddmore0238\Executable_run35b_ddm2.ctl
    
To create the |popy| :ref:`fit_script` here:-

.. code-block:: console

    c:\PoPy\validation\ddmore0238\popy\fit_script.pyml

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: output_popy_data.csv
        
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

    LEVEL_PARAMS:
        GLOBAL:
            params: |
                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
    
    LEVEL_PARAMS:
        GLOBAL:
            params: |
                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 'INDIV' and 'OCCASION' sections:-

.. code-block:: pyml
    
    LEVEL_PARAMS:
        INDIV:
            split_field: ID
            split_dict: "*"
            params: |
                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])
        OCCASION:
            split_field: OCC
            split_dict: "*"
            params: |
                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
    
    LEVEL_PARAMS:
        GLOBAL:
            params: |
                f[PNOISE_var] ~ P 0.036033
                f[ANOISE_var] ~ P 0.0164023
                
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.

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|.  

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 scipy Odeint LSODA 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:
        SCIPY_ODEINT:
            # 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

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) 
        
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: [JOE: {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.
    
        