Error with Parser

[krinaj]

Hello,

I am learning stan for PK PD modeling and running my model that I wrote using your neutropenia example. In my model, I just need to estimate Ka and Kdeg. All other parameters are supplied in the dataset. Below is the code. I am getting error message: SYNTAX ERROR, MESSAGE(S) FROM PARSER:

ERROR at line 111

109: V2[start[i]:end[i]], 110: Q[start[i]:end[i]], 111: parms, 1e-6, 1e-6, 1e-6); ^ 112:
PARSER EXPECTED: ","

I have spent a lot of time trying to debug this, but no luck. Any help/suggestions will be appreciated.

Thanks, Krina

functions{

  real[] twoCptModelODE(real t,
                            real[] x,
                            real[] parms,
                            real[] CL,
                            real[] V1,
                            real[] V2,
                            real[] Q){
    real ka = parms[1];
    real Kdeg = parms[2];
    real k10;
    real k12;
    real k21;

    real dxdt[3];

    dxdt[1] = -ka * x[1] - Kdeg * x[1];
    dxdt[2] = ka * x[1] - (k10 + k12) * x[2] + k21 * x[3];
    dxdt[3] = k12 * x[2] - k21 * x[3];

    return dxdt;
  }

}

data{
  int<lower = 1> nt;
  int<lower = 1> nObsPK;
  int<lower = 1> iObsPK[nObsPK];
  real<lower = 0> AMT[nt];
  int CMT[nt];
  int<lower = 0> EVID[nt];
  real<lower = 0> TIME[nt];
  real RATE[nt];
  vector<lower = 0>[nObsPK] cObs;

  ## data for population model
  int<lower = 1> nSubjects;
  int<lower = 1> start[nSubjects];
  int<lower = 1> end[nSubjects];
  real<lower = 0> CL[nSubjects];
  real<lower = 0> V1[nSubjects];
  real<lower = 0> V2[nSubjects];
  real<lower = 0> Q[nSubjects];

  ## data for priors
  real<lower = 0> kaHatPrior;
  real<lower = 0> kaHatPriorCV;
  real<lower = 0> KdegHatPrior;
  real<lower = 0> KdegHatPriorCV;
}

transformed data{
  vector[nObsPK] logCObs = log(cObs);
  int nTheta = 2;  ## number of ODE parameters
  int nIIV = 2;  ## parameters with IIV
  int nCmt = 3;  ## number of compartments
}

parameters{
  real<lower = 0> kaHat;
  real<lower = 0> KdegHat;
  real<lower = 0> sigma;

  ## IIV parameters
  cholesky_factor_corr[nIIV] L;
  vector<lower = 0>[nIIV] omega;
  matrix[nIIV, nSubjects] etaStd;
}

transformed parameters{
  vector[nt] cHat;
  vector[nObsPK] cHatObs;
  matrix[nt, nCmt] x;
  real<lower = 0> parms[nTheta]; 
  vector[nObsPK] k10;
  vector[nObsPK] k12;
  vector[nObsPK] k21;

  ## variables for Matt's trick
  vector<lower = 0>[nIIV] thetaHat;
  matrix<lower = 0>[nSubjects, nIIV] thetaM; 

  ## Matt's trick to use unit scale
  thetaHat[1] = kaHat; 
  thetaHat[2] = KdegHat;

  thetaM = (rep_matrix(thetaHat, nSubjects) .* 
              exp(diag_pre_multiply(omega, L * etaStd)))';

  for(i in 1:nSubjects) {

  parms[1] = thetaM[i, 1]; # ka
  parms[2] = thetaM[i, 2]; # Kdeg
  k10[i] = CL[i]/V1[i];
  k12[i] = Q[i]/V1[i];
  k21[i] = Q[i]/V2[i];

  x[start[i]:end[i]] = generalOdeModel_rk45(twoCptModelODE, nCmt,
  TIME[start[i]:end[i]], 
  AMT[start[i]:end[i]], 
  RATE[start[i]:end[i]], 
  EVID[start[i]:end[i]], 
  CMT[start[i]:end[i]], 
  CL[start[i]:end[i]], 
  V1[start[i]:end[i]], 
  V2[start[i]:end[i]], 
  Q[start[i]:end[i]], 
  parms, 1e-6, 1e-6, 1e-6);

  cHat[start[i]:end[i]] = x[start[i]:end[i], 2] / V1[i];  ## divide by V1
  }

  cHatObs = cHat[iObsPK];
}

model{
## Priors
kaHat ~ lognormal(log(kaHatPrior), kaHatPriorCV);
KdegHat ~ lognormal(log(KdegHatPrior), KdegHatPriorCV);
sigma ~ cauchy(0, 1);

## Parameters for Matt's trick
L ~ lkj_corr_cholesky(1);
to_vector(etaStd) ~ normal(0, 1);

## observed data likelihood
logCObs ~ normal(log(cHatObs), sigma);
}

generated quantities{
  matrix[nt, nCmt] xPred;
  real<lower = 0> parmsPred[nTheta];
  vector[nt] cHatPred;
  vector<lower = 0>[nObsPK] cHatObsCond;
  vector<lower = 0>[nObsPK] cHatObsPred;

  ## Variables for IIV  
  matrix[nIIV, nSubjects] etaStdPred;
  matrix<lower = 0>[nSubjects, nIIV] thetaPredM;
  corr_matrix[nIIV] rho;

  rho = L * L';
  for(i in 1:nSubjects) {
  for(j in 1:nIIV) {
  etaStdPred[j, i] = normal_rng(0, 1);
  }
  }
  thetaPredM = (rep_matrix(thetaHat, nSubjects) .* 
  exp(diag_pre_multiply(omega, L * etaStdPred)))';

  for(i in 1:nSubjects) {
    parmsPred[1] = kaHat; # ka
    parmsPred[2] = KdegHat; # ka

    xPred[start[i]:end[i]] = generalOdeModel_rk45(twoCptModelODE, nCmt,
                                                  TIME[start[i]:end[i]], 
                                                  AMT[start[i]:end[i]],
                                                  RATE[start[i]:end[i]],
                                                  EVID[start[i]:end[i]],
                                                  CMT[start[i]:end[i]],
                                                  CL[start[i]:end[i]], 
                                                  V1[start[i]:end[i]], 
                                                  V2[start[i]:end[i]], 
                                                  Q[start[i]:end[i]], 
                                                  parmsPred, 
                                                  1e-6, 1e-6, 1e-6);

    cHatPred[start[i]:end[i]] = xPred[start[i]:end[i], 2] / V1[i]; # divide by V1
   }

  ## predictions for observed data records
  cHatObsPred = cHatPred[iObsPK];

  for(i in 1:nObsPK) {
    cHatObsCond[i] = exp(normal_rng(log(fmax(machine_precision(), cHatObs[i])), sigma));
    cHatObsPred[i] = exp(normal_rng(log(fmax(machine_precision(), cHatObsPred[i])), sigma));
  }
}

[billgillespie]

Hi Krina,

The problem appears to be a misunderstanding about the correct signature for generalOdeModel_rk45. The correct signature is:

generalOdeModel_rk45(ODE system, nCmt,
                     time, amt, rate, ii, evid, cmt, addl, ss,
                     theta, biovar, tlag,
                     rel tol, abs tol, max step)

This is for Torsten version 0.83. All of the model parameters need to be passed in theta (parms in your example) even if they are fixed values. Your arguments CL, V1, V2 and Q are not recognized by the compiler. Also keep in mind that Stan only understands the order and types of function arguments. Names are not recognized.

I should also mention that the active github repo has moved to https://github.com/metrumresearchgroup/example-models.

Regards, Bill

[krinaj]

thanks, Bill. Very helpful. I am now able to compile the stan model.