# DDMoRe0093 Conversion Example

## Overview

A real world model, based on a QT study [Cheung2015]. Circadian function with covariates, 59 individuals. 5790 total data rows. No compartment derivatives etc. See DDMoRe: 0093

## Data Conversion

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

    c:\PoPy\validation\ddmore0093\Simulated_dataset.csv

To this data set in PoPy format:-

    c:\PoPy\validation\ddmore0093\popy\popy_data.csv

Using the N2PDat Script located here:-

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

The final data set looks something like Table 50.

Table 50 Main data fields for Nonmem and PoPy
ID | TIME | AMT | EVID | QTCF | MDV | CPP | TYPE
1001 | 0 | 0 | 2 | 0 | 1 | 0 | pred
1001 | 31.2575 | 0 | 0 | 351.7 | 0 | 0 | obs
1001 | 31.2655556 | 0 | 0 | 343.7 | 0 | 0 | obs
1001 | 31.2736111 | 0 | 0 | 354.1 | 0 | 0 | obs
1001 | 31.4072222 | 0 | 0 | 351.5 | 0 | 0 | obs
1001 | 31.4152778 | 0 | 0 | 359 | 0 | 0 | obs
1001 | 33.6672222 | 600 | 1 | 0 | 1 | 0 | dose:dose_bolus
1001 | 34.1683333 | 0 | 2 | 0 | 1 | 1083 | pred
1001 | 34.5841667 | 0 | 0 | 353.6 | 0 | 1447 | obs

Note the main output of the data conversion is the ‘TYPE’ field which has the correct entries for PoPy corresponding to the Nonmem ‘EVID’ field. as follows:-

    2 -> pred
    0 -> obs
    1 -> dose:dose_bolus

Note however that in this model the dosing is not actually used in the Nonmem control file.

## Script Conversion

We describe the manual conversion of the Nonmem script file:-

    c:\PoPy\validation\ddmore0093\Executable_run7b.ctl

To create the PoPy Fit Script here:-

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

The Nonmem script loads the data as follows:-

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

Note the ‘IGNORE’ syntax removes the dosing rows from the data set.

Whereas the PoPy script uses this syntax:-

    FILE_PATHS: {input_data_file: output_popy_data.csv}
    PREPROCESS: |
        if c[TYPE] == 'dose:dose_bolus': return

Which loads the PoPy data and also excludes the dosing rows.

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

    $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 `f[X]` in the PoPy script are:-

    LEVEL_PARAMS:
        GLOBAL:
            params: |
                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
                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 ]
                ]
                f[NOISEVAR] ~ P 0.01802

Additionally PoPy has the ‘INDIV’ section of the LEVEL_PARAMS:-

    LEVEL_PARAMS:
        INDIV:
            params: |
                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]
                )
            split_field: ID
            split_dict: "*"

This section defines `r[BASE]` and `r[AMP24]` etc. These `r[X]` 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:-

    $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:-

    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 `m[X]` variables that you wish to use in subsequent sections, e.g. DERIVATIVES or PREDICTIONS.

The Nonmem error model (also here part of the ‘PRED’ section) is defined as follows:-

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

    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 `c[QTCF]` from the data file with the predicition `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.

The Nonmem control file specifies using FOCE fitting (METHOD=1) below:-

    $ESTIMATION METHOD=1 MAXEVALS=99999 INTER NOABORT PRINT=5

The PoPy equivalent, specifying the JOE fitting method (roughly equivalent to FOCE) is:-

    FIT_METHODS: [JOE: {max_n_main_iterations: 100}]