Dynamically specifying which parameters to estimate in stan

Modeling
, , ,
1

Hello everyone,

I am wondering if it is possible to dynamically specify (e.g., from R) which parameters should be estimated in a Stan model.

Use case

I am developing an RShiny app that allows physicians/clinicians to input a patient’s clinical characteristics, drug administration details, and measured drug concentrations. The app then provides a posterior predictive distribution for the concentration-time profile based on a pharmacokinetic NLME model using Torsten.

In this model, parameters such as clearance are assumed to follow a log-normal distribution, and the measured concentrations inform the posterior distribution of these parameters. See also Simultaneously estimating IIV and IOV estimate for first cycle in pharmacokinetic NLME model using Torsten.

The challenge

Besides inter-individual variability (IIV), my model also accounts for inter-occasion variability (IOV), meaning that for each treatment cycle (i.e., the period between two drug administrations), additional variability is introduced for some model parameters.

Depending on how much data is available for n cycles (user input), we need to estimate n cycle-specific parameters. I want to avoid pre-compiling 20 different Stan models (one for each possible value of n) and would rather dynamically specify which parameters should be estimated directly from R.

Is this somehow possible?

Here is my code:

/////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////// PREAMBLE ///////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////

// stan representation of run5025
// 1. Based on:             run5013
// 2. Description:          train / log / 3cmt / iiv:CL/Q3 / iov:CL/Q3 / cov: eGFR on CL (pow) / M3 / no CWRES > 5 / LAPLACIAN
// 3. Structural model:     3cmt
// 4. Dispositon:           fo
// 5. IIV:                  CL/Q3
// 6. IOV:                  CL/Q3
// 7. Covariates:           eGFR on CL (pow)
// 8. RUV:                  exp
// 9. Estimation:           LAPLACIAN
// 10. Author:              Marian Klose

/////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////// FUNCTIONS BLOCK ////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////

// functions data block
functions {
    vector ode_rhs(
        real t,
        vector x,
        array[] real parms,
        array[] real x_r,
        array[] int x_i
    ){
        // extract parmeters from prms object in defined order, CL/V1/Q2/V2/Q3/V3
        real CL = parms[1];
        real V1 = parms[2];
        real Q2 = parms[3];
        real V2 = parms[4];
        real Q3 = parms[5];
        real V3 = parms[6];

        // define rate constants based on params
        real k10 = CL / V1;
        real k12 = Q2 / V1;
        real k21 = Q2 / V2;
        real k13 = Q3 / V1;
        real k31 = Q3 / V3;

        // define ODE system
        vector[3] y;
        y[1] = - k10 * x[1] - k12 * x[1] - k13 * x[1] + k21 * x[2] + k31 * x[3];
        y[2] = k12 * x[1] - k21 * x[2];
        y[3] = k13 * x[1] - k31 * x[3];

        // return the amount in the compartments
        return y;
    }
}

/////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////// DATA BLOCK /////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////

// input data block (passed from R)
data {
    // events
    int<lower = 1> nEvents;                 // number of events

    // observations
    int<lower = 1> nObs;                    // number of observations
    array[nObs] int<lower = 1> iObs;        // indices of events which are observations
    vector<lower = 0>[nObs] cObs;           // observed concentration (the dv) used to generate posterior

    // NONMEM standard event items
    array[nEvents] int<lower = 1> cmt;      // compartment number
    array[nEvents] int evid;                // event id (0=observation, 1=dose, 2=other)
    array[nEvents] int addl;                // number of additional identical doses given
    array[nEvents] int ss;                  // is it steady-state dosing? (0=false, 1=true)
    array[nEvents] real amt;                // dose amount administered at this time
    array[nEvents] real time;               // time of observation/administration
    array[nEvents] real rate;               // rate of drug infusion (0 for bolus administration)
    array[nEvents] real ii;                 // interdose interval: time between additional doses

    // typical parameter values
    real<lower = 0> TVCL;                   // typical value of central clearance
    real<lower = 0> TVV1;                   // typical value of central volume
    real<lower = 0> TVQ2;                   // typical value of shallow intercompartmental clearance
    real<lower = 0> TVV2;                   // typical value of shallow peripheral volume
    real<lower = 0> TVQ3;                   // typical value of deep intercompartmental clearance
    real<lower = 0> TVV3;                   // typical value of deep peripheral volume

    // fixed covariate effects
    real<lower = 0> TH_EGFR_CL;             // fixed effect of creatinine clearance on clearance

    // measured covariates 
    array[nEvents] real<lower = 0> EGFR;    // estimated glomerular filtration rate [ml/min]

    // cycle information
    array[nEvents] int<lower = 0> CYCLE;    // cyle number, needed for IOV

    // iiv standard deviations
    real<lower = 0> IIV_SD_CL;              // IIV (sd) on central clearance
    real<lower = 0> IIV_SD_Q3;              // IIV (sd) on intercompartimental clearance Q3

    // iov standard deviations
    real<lower = 0> IOV_SD_CL;              // IOV (sd) on central clearance
    real<lower = 0> IOV_SD_Q3;              // IIV (sd) on intercompartimental clearance Q3

    // ruv parameters
    real<lower = 0> EXP_RUV_SD;             // RUV (sd) for exponential error model
}

/////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////// TRANSFORMED DATA BLOCK /////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////

// transformed data block
transformed data {
    vector[nObs] logCObs = log(cObs);           // log-transformed observations
    int nTheta = 6;                             // number of parameters for ODE system
    int nCmt = 3;                               // number of compartments (1 = central, 2 = shallow peripheral, 3 = deep peripheral)
}

/////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////// PARAMETERS BLOCK ///////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////

// parameters block
parameters {
    // clearance
    real ETA_CL;              // ind. iiv estimate
    real KAPPA_CL_C1;         // ind. iov estimate (cycle 1)

    // intercompartmental clearance q3
    real ETA_Q3;              // individual iiv q3 estimates
    real KAPPA_Q3_C1;         // individual iov q3 estimates
}

/////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////// TRANSFORMED PARAMETERS BLOCK ///////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////

// transformed parameters block
transformed parameters {
    // predicted drug concentrations
    matrix<lower = 0>[nCmt, nEvents] x;  // drug amounts in each compartment over time
    row_vector<lower = 0>[nEvents] cHat; // row vector, one element per event
    vector<lower = 0>[nObs] cHatObs;     // column vector, one element per *observation*

    // retrieve correct KAPPA_CL and KAPPA_Q3 values
    real KAPPA_CL = 0;
    real KAPPA_Q3 = 0;

    // 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 
        // !!FIX ME: ADJUST KAPPA ACCORDING TO CYCLE; FOR NOW CYCLE 1 ONLY!!
      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
    );

    // compute drug concentrations
    for (i in 1:nEvents) {
        cHat[i] = x[1, i] / theta[i, 2];        // need to loop even if v1 remains constant
    }          
    cHatObs = to_vector(cHat[iObs]);            // transform to column vector & keep relevant cells only
    vector[nObs] log_cHatObs = log(cHatObs);    // log-transformed predictions
}

/////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////// MODEL BLOCK ////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////

// model block
model {
    // informative priors for clearance
    ETA_CL ~ normal(0, IIV_SD_CL);              // IIV (sd) on central clearance
    KAPPA_CL_C1 ~ normal(0, IOV_SD_CL);         // IIV (sd) on intercompartimental clearance Q3

    // informative priors for intercompartmental clearance Q3
    ETA_Q3 ~ normal(0, IIV_SD_Q3);              // IOV (sd) on central clearance
    KAPPA_Q3_C1 ~ normal(0, IOV_SD_Q3);         // IIV (sd) on intercompartimental clearance Q3

    // likelihood function for exponential error model (expects SD for normal distribution)
    logCObs ~ normal(log_cHatObs, EXP_RUV_SD);
}


/////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////// GENERATED QUANTITIES BLOCK /////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////

// generated quantities block
generated quantities {
    // sample prior etas
    real prior_ETA_CL = normal_rng(0, IIV_SD_CL);           // IIV (sd) on central clearance
    real prior_KAPPA_CL_C1 = normal_rng(0, IOV_SD_CL);      // IIV (sd) on intercompartimental clearance Q3
    real prior_ETA_Q3 = normal_rng(0, IIV_SD_Q3);           // IOV (sd) on central clearance
    real prior_KAPPA_Q3_C1 = normal_rng(0, IOV_SD_Q3);      // IIV (sd) on intercompartimental clearance Q3

    // sample posterior
    array[nObs] real logcObsPred;  // introduce observations for posterior
    for(i in 1:nObs) {
        logcObsPred[i] = normal_rng(log_cHatObs[i], EXP_RUV_SD);
    }
}

This is currently based on a single cycle/occasion, but I am in the process of extending it to multiple cycles. At the moment, we define the following parameters:

parameters {
    // clearance
    real ETA_CL;              // ind. iiv estimate
    real KAPPA_CL_C1;         // ind. iov estimate (cycle 1)

    // intercompartmental clearance q3
    real ETA_Q3;              // individual iiv q3 estimates
    real KAPPA_Q3_C1;         // individual iov q3 estimates (cycle 1)
}

(...)

transformed parameters {
    if (CYCLE[i] == 1) {
        KAPPA_CL = KAPPA_CL_C1;
        KAPPA_Q3 = KAPPA_Q3_C1;
    } 
}

If we have one cycle, this is sufficient. However, if we have three cycles, we would already need two additional sets of parameters:

parameters {
    // clearance
    real ETA_CL;              // ind. iiv estimate
    real KAPPA_CL_C1;         // ind. iov estimate (cycle 1)
    real KAPPA_CL_C2;         // ind. iov estimate (cycle 2)
    real KAPPA_CL_C3;         // ind. iov estimate (cycle 3)

    // intercompartmental clearance q3
    real ETA_Q3;              // individual iiv q3 estimates
    real KAPPA_Q3_C1;         // individual iov q3 estimates (cycle 1)
    real KAPPA_Q3_C2;         // individual iov q3 estimates (cycle 2)
    real KAPPA_Q3_C3;         // individual iov q3 estimates (cycle 3)
}

(...)

transformed parameters {
    if (CYCLE[i] == 1) {
        KAPPA_CL = KAPPA_CL_C1;
        KAPPA_Q3 = KAPPA_Q3_C2;
    } else if (CYCLE[i] == 2){
        KAPPA_CL = KAPPA_CL_C2;
        KAPPA_Q3 = KAPPA_Q3_C2;
    } else if (CYCLE[i] == 3){
        KAPPA_CL = KAPPA_CL_C3;
        KAPPA_Q3 = KAPPA_Q3_C3;
    }
}

My current idea is to define one single model with parameters for 20 cycles/occasions to accommodate even extreme scenarios. I would like to always use this model, regardless of whether data for only 1 or 20 cycles is provided.

Question

Is there a way to dynamically specify which parameters should be estimated and which should be ignored when fewer cycles are provided?

Let’s assume a scenario where we have data for 2 cycles, but KAPPA_CL_C3 and KAPPA_Q3_C3 are still defined in the model. Since CYCLE[i] will never take the value 3, the data contains no information for KAPPA_CL_C3 and KAPPA_Q3_C3, making these cycle 3 related parameters complely insensitive to the data.

  • Will this destabilize the model, or will the posterior for these parameters simply remain equal to the prior without additional convergence issues?
  • Does it make sense to assign a very tight prior to parameters that do not need to be estimated? Alternatively, could I even set a prior density of 0 for these parameters?
  • Is there any built-in way to disable parameters in the model when they are not needed?

Thanks a lot for your support!

2

I think you could vectorize your parameters, with the vector length changing depending on number of cycles. In the data block, I would add

int ncycles;

Then, when you declare parameters, use lines like

vector[ncycles] KAPPA_CL;

Then, instead of the if-else statements in transformed parameters, you could choose a specific parameter via

KAPPA_CL[ CYCLE[i] ];
1 Like
3

Thanks a lot for getting back to me, @edm. Sounds like a handy solution, it would be neat if it works as suggested. I will give it a try soon and report back if it works!

4

The latter—the posterior will be equal to the prior so it shouldn’t cause any convergence issues unless your priors are very wide or if you have very high dimensionality of uninformative parameters, at which point they might slow things down a little.

No. This will have problems if you make them too tight. Stan has to figure out the scale of parameters and if you make them very tight (e.g., normal(0, 1e-10), then random inits will be problematic because they’re so far into the tail.

That’s even worse—this will grind HMC to a halt.

Not directly—we usually just do what @edm suggested.

1 Like
5

Thanks a lot, this solution worked for me. Exactly what I needed! Here are the updated blocks:

data {
    // cycles
    int<lower = 1> nCycles;                 // number of cycles
}

parameters {
    // clearance
    real ETA_CL;                                 // ind. iiv estimate
    vector[nCycles] KAPPA_CL;                    // ind. iov estimates (vectorized)

    // intercompartmental clearance q3
    real ETA_Q3;                                 // ind. iiv q3 estimates
    vector[nCycles] KAPPA_Q3;                    // ind. iov q3 estimates (vectorized)
}

transformed parameters {
    array[nEvents, nTheta] real theta;
    for (i in 1:nEvents) {

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

model {
    // informative priors for clearance
    ETA_CL ~ normal(0, IIV_SD_CL);              // IIV (sd) on central clearance
    KAPPA_CL ~ normal(0, IOV_SD_CL);            // IIV (sd) on central clearance (vectorized)

    // informative priors for intercompartmental clearance Q3
    ETA_Q3 ~ normal(0, IIV_SD_Q3);              // IOV (sd) on intercompartimental clearance Q3
    KAPPA_Q3 ~ normal(0, IOV_SD_Q3);            // IIV (sd) on intercompartimental clearance Q3 (vectorized)

    // likelihood function for exponential error model (expects SD for normal distribution)
    logCObs ~ normal(log_cHatObs, EXP_RUV_SD);
}

It is good to know that this vectorized approach works, makes my life much easier!

Thanks a lot for the additional explanations, makes all sense now!

1 Like