Time-varying covariate implementation

Modeling
4

Hey Andras, I am also new to the stan and Torsten world but very excited to continue this path. I don’t have a perfect solution for you, but I think the key is to adjust the theta (or params) object and move it from a 1D array (real[]) to a 2D array (real[ , ]). If you have a look into the current user guide of Torsten (Torsten/docs/torsten_users_guide.pdf at master · metrumresearchgroup/Torsten · GitHub) under page 12 at Table 4.2:

Table 4.2. PMX model parameter overloadings. One can use 1d array real[] to indicate constants of all events, or 2d array real[ , ] so that the ith row of the array describes the model arguments for time interval (ti−1, ti), and the number of the rows equals to the size of time.

There it is described that the theta (or params) object can either be a 1d array real[] if your parameters stay constant for one individual and a 2d array real[ , ] if you have different parameters for one individual for a specific time interval. To my understanding this is exactly the case if we have time-varying covariates and/or inter-occasion variability. Then the params/thetas need to change within one individual. Currently you still use the 1d array real[] version in your code:

for(j in 1:n_subjects){
      real wt_adjustment_cl = (WT[j]/70)^0.75;
      array[4] real params = {CL[j]* wt_adjustment_cl, Q[j], VC[j], VP[j]};
      array[4] real params_tv = {TVCL* wt_adjustment_cl, TVQ, TVVC, TVVP};
      x_ipred[subj_start[j]:subj_end[j],] =
        pmx_solve_rk45(two_cmt_ir1_ode,
                       n_cmt,
                       time[subj_start[j]:subj_end[j]],
                       amt[subj_start[j]:subj_end[j]],
                       rate[subj_start[j]:subj_end[j]],
                       ii[subj_start[j]:subj_end[j]],
                       evid[subj_start[j]:subj_end[j]],
                       cmt[subj_start[j]:subj_end[j]],
                       addl[subj_start[j]:subj_end[j]],
                       ss[subj_start[j]:subj_end[j]],
                       params)';

(...)
}

This will tell Torsten that you assume constant parameters for one individual. If you want to add time-varying covariates, I think you need to add another loop within one individual and loop over each event (= 1 row in you NONMEM-style dataset) and comput the respective parameter value at that point. Here is how it looks in my code:

    // define the time-varying parameter matrix theta [nEvents, nTheta] (n rows = number of time points, per row one set of nTheta params)
    // see docs for details on 1d and 2d theta elements (2d needed for time-varying covariates)
    array[nEvents, nTheta] real theta;
    for (i in 1:nEvents) {
        // extract correct kappas 
        if (CYCLE[i] == 1) {
            KAPPA_CL = KAPPA_CL_C1;
            KAPPA_Q3 = KAPPA_Q3_C1;
        } else {
            KAPPA_CL = KAPPA_CL_C1;
            KAPPA_Q3 = KAPPA_Q3_C1;
        }

        // store thetas
        theta[i, 1] = TVCL * pow(EGFR[i] / 84, TH_EGFR_CL) * exp(ETA_CL + KAPPA_CL);
        theta[i, 2] = TVV1;
        theta[i, 3] = TVQ2;
        theta[i, 4] = TVV2;
        theta[i, 5] = TVQ3 * exp(ETA_Q3 + KAPPA_Q3);
        theta[i, 6] = TVV3;
    }

    // solve the ode system (x = drug amounts in each compartment over time)
    x = pmx_solve_rk45(
        ode_rhs,                          // ODE right hand side system
        nCmt,                             // number of compartments
        time,                             // time of events
        amt,                              // amount of drug administered
        rate,                             // rate of drug administered
        ii,                               // interdose interval
        evid,                             // event id
        cmt,                              // compartment number
        addl,                             // number of additional doses
        ss,                               // steady state flag
        theta,                            // individual parameters
        1e-5,                             // relative tolerance
        1e-8,                             // absolute tolerance    
        1e5                               // maximum stepsize
    );

In this case, EGFR[i] is the time-varying covariate (in mymodel I just have 1 indvidual, see below). Instead of defining

array[4] real params = {CL[j]* wt_adjustment_cl, Q[j], VC[j], VP[j]};

you would need to define per ith row of the respective individual something like this:

array[nEvents, nTheta] real params;
theta[i, 1] = CL[j] * (WT[i,j]/70)^0.75;      // CL
theta[i, 2] = Q[j];                           // Q
theta[i, 3] = VC[j];                          // VC
theta[i, 4] = VP[j];                          // VP

It is surely a bit tricky to get the indexing right and to find a smart way to code it. But this should be the inituiton behind it. The trick will be to pass all objects in form of xx[i,j] and i need to match the elements of the time array.

Please also have a look at the full model code I have posted in another question if you are interested: Simultaneously estimating IIV and IOV estimate for first cycle in pharmacokinetic NLME model using Torsten. It is a bit different since I do not fit a full NLME model there and rather try to replace the MAP estimation (mode of posterior) with a full Bayesian approach for a MIPD scenario. This means I only have 1 individual. But the concept for time-varying covariates and IOV remains the same. Let me know if that was helpful or rather confusing!