.. _ddmore0093_conv_example:
   
DDMoRe0093 Conversion Example
#################################################

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

A real world model, based on a QT study [Cheung2015]_. 
The model incorporates a circadian function with no compartment derivatives. 
There are 59 individuals and 5790 total data rows. See :ddmore_link:`0093`

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

Data Source & Preprocessing
-------------------------------

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

.. code-block:: none

    $DATA Simulated_dataset.csv
       IGNORE = #
       IGNORE = (EVID == 1)

Note the 'IGNORE' syntax removes the dosing rows from the data set, effectively
preprocessing the dataset at the point of loading.
       
We achieve the same thing in |popy| using the :ref:`preprocess` block:

.. code-block:: pyml

    FILE_PATHS: { input_data_file: Simulated_dataset.csv }
    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()

        # Suppress warning about having doses we don't use
        # by ignoring dose rows (EVID==1) altogether.
        if c[TYPE] == 'dose': return

which loads the |popy| data, converts ``EVID`` values to ``TYPE`` values 
and also excludes the dosing rows.
    
Fixed/Random Effects Specification
-----------------------------------

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

.. code-block:: none

    $THETA (0, 372.6)      ; 1 Baseline QTcF [ms]
    $THETA (0, 0.01844)    ; 2 Amplitude  24h
    $THETA (0, 3.62)       ; 3 Peak shift 24h
    $THETA (0, 0.01392)    ; 4 Amplitude  12h circadian rhythm
    $THETA (0, 1.301)      ; 5 Peak shift 12h circadian rhythm
    $THETA (0, 0.003441)   ; 6 Slope (linear effect) parameter

    $OMEGA  0.0392         ; 1 Baseline QTcF - (BSV)
    $OMEGA  0.3523         ; 2 Amplitude 24h  - BSV
    $OMEGA  2.007          ; 3 Peak shift 24h - BSV
    $OMEGA  0.3848         ; 4 Amplitude 12h  - BSV
    $OMEGA  0.993          ; 5 Peak shift 12h - BSV
    $OMEGA  0.0001         ; 6 Slope          - BSV

    $SIGMA  0.01802        ; 1 Proportional residual error [ms]
    
The equivalent |fx| in the |popy| script are:-

.. code-block:: pyml

    EFFECTS:
        POP: |
            # THETAs
            f[BASE] ~ P 372.6
            f[AMP24] ~ P 0.01844
            f[SHFT24] ~ P 3.62
            f[AMP12] ~ P 0.01392
            f[SHFT12] ~ P 1.301
            f[SLOPE] ~ P 0.003441

            # OMEGAs
            f[BASE_bsv,AMP24_bsv,SHFT24_bsv,AMP12_bsv,SHFT12_bsv,SLOPE_bsv] ~ diag_matrix() [
                [ 0.0392, 0.3523, 2.007, 0.3848, 0.993, 0.0001 ]
            ]

            # SIGMAs
            f[NOISEVAR] ~ P 0.01802
        ...

Additionally |popy| has the 'ID' section of the |effects|:-

.. code-block:: pyml

    EFFECTS:
        ...
        ID: |
            r[ BASE, AMP24, SHFT24, AMP12, SHFT12, SLOPE ] ~ mnorm(
                [0,0,0,0,0,0],
                f[BASE_bsv,AMP24_bsv,SHFT24_bsv,AMP12_bsv,SHFT12_bsv,SLOPE_bsv]
            )
                
This section defines :pyml:`r[BASE]` and :pyml:`r[AMP24]` etc. 
These |rx| definitions are a more explicit version of |nonmem| implicitly 
creating ETA(X) variables for each of the omega variables. With |nonmem| you 
have to remember what each of the ETA indices represent. Hence the heavily 
commented 'PRED' block below, which defined the variables for each individual
in the |nonmem| script.
                
PK Model Parameters Specification
----------------------------------

.. code-block:: none
       
    $PRED
       BASE   = THETA(1) * EXP(ETA(1))                       ; Baseline QTcF with between subject varaibility (BSV)
       AMP24  = THETA(2) * EXP(ETA(2))                       ; The amplitude of the 24h circadian rhythm
       SHFT24 = THETA(3) + ETA(3)                            ; The peak shift of the 24h circadian rhythm
       AMP12  = THETA(4) * EXP(ETA(4))                       ; The amplitude of the 12h circadian rhythm
       SHFT12 = THETA(5) + ETA(5)                            ; The peak shift of the 12h circadian rhythm

       CIRC24 = AMP24 * COS(2 * 3.14 * (TIME - SHFT24)/24)   ; 24 hour circadian rhythm
       CIRC12 = AMP12 * COS(2 * 3.14 * (TIME - SHFT12)/12)   ; 12 hour circadian rhythm

       RYTM   = BASE * (1 + CIRC24 + CIRC12)                 ; Change in baseline QTcF over the day due to circadian rhythm

       SLOPE  = THETA(6) + ETA(6)                            ; Linear effect
       EFF    = SLOPE * CPP                                  ; Linear effect

       ...
    
The |model_params| section of the |popy| script is very similar:-

.. code-block:: pyml

    MODEL_PARAMS: |
       
        BASE   = f[BASE] * exp(r[BASE])
        AMP24  = f[AMP24] * exp(r[AMP24])
        SHFT24 = f[SHFT24] + r[SHFT24]
        AMP12  = f[AMP12] * exp(r[AMP12])
        SHFT12 = f[SHFT12] + r[SHFT12]
        
        CIRC24 = AMP24 * COS(2 * 3.14 * (c[TIME] - SHFT24)/24)
        CIRC12 = AMP12 * COS(2 * 3.14 * (c[TIME] - SHFT12)/12)
        m[RYTM] = BASE * (1 + CIRC24 + CIRC12)
        
        SLOPE = f[SLOPE] + r[SLOPE]
        m[EFF] = SLOPE * c[CPP]
        m[NOISEVAR] = f[NOISEVAR]
        
The above is almost identical apart from the line 'm[NOISEVAR] = f[NOISEVAR]'.
Note in |model_params| it is necessary to define |mx| variables that you wish
to use in subsequent sections, e.g. |derivatives| or |predictions|.

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

This model does not require or use differential equations so the ``$DES``
section is absent in the |nonmem| file and the :ref:`derivatives` section
is absent from the |popy| file.

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

The |nonmem| error model (also here part of the 'PRED' section) 
is defined as follows:-

.. code-block:: none

   IPRED  = RYTM  + EFF        ; Linear direct effect model
   W      = IPRED * SIGMA(1,1)
   IWRES  = (QTCF - IPRED)/W

   Y      = IPRED + IPRED*EPS(1)
    
The equivalent in |popy| is the |predictions| block:-

.. code-block:: pyml

    PREDICTIONS: |
        p[CONC_PLASMA] = m[RYTM] + m[EFF]
        var = m[NOISEVAR] * p[CONC_PLASMA]**2
        c[QTCF] ~ norm(p[CONC_PLASMA], var)

Note |popy| explicitly defines the likelihood by comparing :pyml:`c[QTCF]` 
from the data file with the predicition :pyml:`p[CONC_PLASMA]`, using a 
normal distribution with proportional variance.

This is more explicit than the |nonmem| line 'Y      = IPRED + IPRED*EPS(1)', 
which implicitly uses the 'DV' |nonmem| magic variable for Y.
The likelihood is expressed as a function of EPS normal distributions.

Parameter Fitting Methods Specification
------------------------------------------

The |nonmem| control file specifies using FOCE fitting (METHOD=1) below:-
        
.. code-block:: none
        
   $ESTIMATION METHOD=1 MAXEVALS=99999 INTER NOABORT PRINT=5
   
The |popy| equivalent, specifying the |nd| fitting method is:-

.. code-block:: pyml

    FIT_METHODS: [ND: {max_n_main_iterations: 100}]
