Map_rect for censored data

Modeling
1

Hi, I’m new to using map_rect and have been facing issues with trying to code the function to calculate the log-posterior of each shard with censored data (due to the limit of detection). I’m using a simulated dataset of 3 patients’ viral load with 25 data points each while in the process of coding this. So each patient is meant to be a shard. The data vLoad is the viral load on each day (1-25) for patient 1 to 3, appended one after another. The LOD looks similar, except that for each patient’s daily viral load, they have a corresponding LOD value (which is the same throughout the 25 days. E.g. 1500 for patient 1, 600 for patient 2, etc)

This is what I have done so far for the code with map_rect

functions {
  real[] OneSys(real t, real[] y, real[] params, real[] x_r, int[] x_i) {
    real dydt[5];
  
    dydt[1] = -params[1]*y[1]*y[4];
    dydt[2] = params[1]*y[1]*y[4] - params[2]*y[2];
    dydt[3] = params[2]*y[2] - params[3]*y[3];
    dydt[4] = 15*y[3] - 20*y[4] - params[4]*y[5]*y[4];
    dydt[5] = params[5]*(y[2]+y[3])*y[5];
  
    return dydt;
  }
  
// function lp calculates log posterior on one shard (patient) at a time
  vector lp(vector pred, real par, real[] lod, real[] xr) {
    real sigma = par;
    real ll;
    int N = 25;
    real measure = xr[1:25];
    real limit = lod[1:25];
    real ypred = pred[1:25];
    
    for (i in 1:N) {
      if (measure[i] <= log(limit[i])) {
        ll = normal_lcdf(log(limit[i]) | ypred[i], sigma);
      }
      if (measure[i] > log(limit[i])) {
        ll = normal_lpdf(measure[i] | ypred[i], sigma);
      }
      
    }
    return [ll]';
  }
}


data {
  int<lower=0> L; //rows in dataset
  int<lower=0> M; //no. of patients
  vector[L] LOD;
  real vLoad[L,1];
  real y0[5];
  real t0;
  real ts[L];
  int patID[L];
  int patCount[M];
  int pos_pat[M];
}

transformed data {
  real x_r[0];
  int x_i[0];
  real xr[M,25] = {vLoad[1:25,1], vLoad[26:50,1], vLoad[51:75,1]};
  real limit[M,25] = {LOD[1:25], LOD[26:50], LOD[51:75]};
}

parameters {
  real<lower = 10^-10, upper = 10^-6> beta_par;
  real<lower = 0, upper = 10> k_par;
  real<lower = 0, upper = 1> delta_par;
  real<lower = 0, upper = 1> gamma_par[M];
  real<lower = 0, upper = 1> gamma_mu;
  real<lower = 0, upper = 0.1> gamma_sigma;
  real<lower = 0, upper = 1> omega_par[M];
  real<lower = 0, upper = 1> omega_mu;
  real<lower = 0, upper = 0.1> omega_sigma;
  real<lower = 0, upper = 20> sigma;
}

transformed parameters {
  real predVal[L];
  real params[5];
  real logpred[M,25];
  
  params[1] = beta_par;
  params[2] = k_par;
  params[3] = delta_par;
  
  for (m in 1:M) {
    real value_m[patCount[m],5];
    params[4] = gamma_par[m]; //patient-specific
    params[5] = omega_par[m]; //patient-specific
    
    value_m = integrate_ode_bdf(OneSys, y0, t0, segment(ts,pos_pat[m],patCount[m]), params, x_r, x_i);
    // store predicted viral load values
    for (j in 1:patCount[m]) {
      predVal[pos_pat[m] + j - 1] = value_m[j,4];
    }
  }
  logpred = {log(predVal[1:25]), log(predVal[26:50]), log(predVal[51:75])};
}

model {
  gamma_par ~ normal(gamma_mu,gamma_sigma);
  omega_par ~ normal(omega_mu,omega_sigma);
  beta_sigma ~ cauchy(0,10^-7);
  gamma_sigma ~ cauchy(0,0.01);
  omega_sigma ~ cauchy(0,10^-7);
  sigma ~ cauchy(0,10);
  inc_per ~ normal(5,2);
  
  target += sum(map_rect(lp, logpred, sigma, limit, xr));
}

I get the following error:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable definition base type mismatch, variable declared as base type real variable definition has base type real[  ]Variable definition dimensions mismatch, definition specifies 0, declaration specifies 1 error in 'model310c6c4b6013_20220718_maprect' at line 20, column 27
  -------------------------------------------------
    18:     real ll;
    19:     int N = 25;
    20:     real measure = xr[1:25];
                                  ^

I understand that this is likely a trivial problem but I’m just really stumped – Any help on this would be much appreciated! Thank you!

2

Not exactly the cleanest error message, I’m afraid. The two first lines are key. What they’re saying is that the variable you’re trying to assign to, measure, is of type real, whereas the value defined on the right-hand side is an array of real. Assuming you meant 1:N instead of hard coding 1:25, that’d be:

array[N] real measure = xr[1:N];

Then, of course, xr is now an N-dimensional array and needs to be used as such. Generally, I tend to define vectors rather than arrays, but xr itself was defined as an array.

3

I see! Thanks for the explanation and troubleshoot! I replaced the line for measure with yours but get this error

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable "array" does not exist.

Which I’m assuming is because Rstan is using an older version of Stan? I’m currently unable to use Cmdstan on the HPC so I unfortunately can’t switch for now. Is it possible to work around this by not having xr as an array? I wasn’t sure if it was possible since the Stan guide for map_rect stated that map function needs data arrays?