Error on real

Getting this error. Can you please suggest any solution?

// Define functions for the 2-compartment ODE model
functions {
  // Function defining the ODE system for amounts in compartments
  // y[1]: amount in depot (for SC)
  // y[2]: amount in central (Vc)
  // y[3]: amount in peripheral (Vp)

  // theta: parameters [CL, Vc, Vp, Q, ka]
  // x_r: real data (not used in this ODE function, but required by signature)
  // x_i: integer data (not used in this ODE function, but required by signature)
  vector two_cmt_ode(real t, vector y, real[] theta, real[] x_r, int[] x_i) {
    real current_CL = theta[1];
    real current_Vc = theta[2];
    real current_Vp = theta[3];
    real current_Q = theta[4];
    real current_ka = theta[5]; // Only used if depot has amount

    // Add a small floor to volumes to prevent division by near-zero during integration
    real safe_Vc = fmax(current_Vc, 1e-12); // Use 1e-12 or a similarly small positive number
    real safe_Vp = fmax(current_Vp, 1e-12);

    // Calculate concentrations in central and peripheral compartments
    real conc_c = y[2] / safe_Vc; // Amount in central / Safe Volume of central
    real conc_p = y[3] / safe_Vp; // Amount in peripheral / Safe Volume of peripheral

    // Calculate rates of transfer and elimination
    real rate_cl = current_CL * conc_c; // Elimination from central
    real rate_q_out = current_Q * conc_c; // Transfer rate from central to peripheral
    real rate_q_in = current_Q * conc_p;   // Transfer rate from peripheral to central
    real rate_abs = current_ka * y[1];     // Absorption rate from depot

    // Rates of change for [Depot, Central, Peripheral]
    vector[3] dydt = rep_vector(0.0, 3); // Correct syntax for declaring and initializing a sized vector

    // Rate of change for Depot (only applies if there's drug in the depot)
    dydt[1] = -rate_abs; // Drug leaves the depot by absorption

    // Rate of change for Central compartment
    dydt[2] = rate_abs     // Absorption into central (if SC)
              + rate_q_in  // Transfer from peripheral to central
              - rate_q_out // Transfer from central to peripheral
              - rate_cl;   // Elimination from central

    // Rate of change for Peripheral compartment
    dydt[3] = rate_q_out - rate_q_in; // Transfer between central and peripheral

    return dydt;
  }
}

data {
  int<lower=0> N; // Total number of observations (nIV + nSC)
  int<lower=1> N_subjects; // Number of subjects
  int<lower=1> N_drugs; // Number of drugs

  int<lower=1, upper=N_subjects> subject[N]; // Subject index for each observation
  int<lower=1, upper=N_drugs> drug[N]; // Drug index for each observation
  int<lower=0, upper=1> route[N]; // Route (0 for IV, 1 for SC)
  vector<lower=0>[N] time; // Time points for all observations
  vector<lower=0>[N] conc; // Concentrations for all observations
  vector<lower=0>[N] dose; // Dose for each observation (0 for non-dose points, actual dose for dose points)

  vector<lower=0>[N_subjects] subject_wt; // Subject weight
  real<lower=0> MEDWT; // <<-- Median weight (e.g., 70 kg)

  // Remove ODE solver control parameters from data block
  // They are moved to transformed data
}
transformed data {
  // Declare and assign ODE solver control parameters here as fixed constants
  real rel_tol = 1e-6;
  real abs_tol = 1e-6;
  int max_steps = 10000000; // Assign as integer literal (1e7)
}
parameters {
  // Typical population parameters (log scale) for standard 2-cmt model
  real log_CL; // Log of Typical elimination clearance
  real log_Vc; // Log of Typical central volume
  real log_Vp; // Log of Typical peripheral volume
  real log_Q; // Log of Typical intercompartmental clearance
  real log_ka; // Log of Absorption rate constant for SC

  real<lower=0,upper=1> F; // Bioavailability for SC

  real<lower=0,upper=100> sigma; // Residual error standard deviation (on log scale)

  // Random effects (on log scale)
  vector[N_drugs] eta_drug_CL;
  vector[N_subjects] eta_subject_CL;
  vector[N_drugs] eta_drug_Vc;
  vector[N_subjects] eta_subject_Vc;
  // Keep independent between-drug and between-subject Vp random effects
  vector[N_drugs] eta_drug_Vp;
  vector[N_subjects] eta_subject_Vp;
  vector[N_drugs] eta_drug_Q; // Between-drug Q rand effects
  vector[N_subjects] eta_subject_Q; // Between-subject Q rand effects
  vector[N_drugs] eta_drug_ka; // Between-drug ka rand effects
  vector[N_subjects] eta_subject_ka; // Between-subject ka rand effects


  // Correlation structure (Expanded to include CL, Vc, Vp)
  cholesky_factor_cov[3] L_drug_within_subject_cov; // Changed from [2] to [3] for CL, Vc, Vp
  matrix[N_drugs, 3] z_drug_subject[N_subjects]; // Changed from [N_drugs, 2] to [N_drugs, 3]

  // Covariate coefficients (on log scale) - Assuming weight on CL, Vc, and Vp
  real beta_CL_WT;
  real beta_Vc_WT;
  real beta_Vp_WT; // Uncommented and added beta for Vp
  // Add beta_Q_WT, beta_ka_WT priors if needed
}

transformed parameters {
  // Population parameters (original scale)
  real CL = exp(log_CL);
  real Vc = exp(log_Vc);
  real Vp = exp(log_Vp);
  real Q = exp(log_Q);
  real ka = exp(log_ka);

  // Effective parameters for each observation
  vector[N] CL_effective;
  vector[N] Vc_effective;
  vector[N] Vp_effective;
  vector[N] Q_effective;
  vector[N] ka_effective;

  // Combined predictions for all observations
  vector[N] pred_conc;

  // Transform within-subject/drug deviates for CL, Vc, and Vp
  vector[N_drugs] eta_drug_CL_within[N_subjects];
  vector[N_drugs] eta_drug_Vc_within[N_subjects];
  vector[N_drugs] eta_drug_Vp_within[N_subjects]; // Added within Vp


  for (s in 1:N_subjects) {
    for (d in 1:N_drugs) {
      // Joint deviation vector now has size 3
      vector[3] eta_drug_joint = L_drug_within_subject_cov * z_drug_subject[s, d]';
      eta_drug_CL_within[s, d] = eta_drug_joint[1];
      eta_drug_Vc_within[s, d] = eta_drug_joint[2];
      eta_drug_Vp_within[s, d] = eta_drug_joint[3]; // Extract Vp component
    }
  }

  // Calculate effective parameters for each observation
  for (i in 1:N) {
    real current_subject_wt = subject_wt[subject[i]];
    real log_wt_ratio = log(current_subject_wt) - log(MEDWT);

    // Effective CL (with all effects)
    CL_effective[i] = CL * exp(eta_drug_CL[drug[i]] + eta_subject_CL[subject[i]] + eta_drug_CL_within[subject[i], drug[i]] + beta_CL_WT * log_wt_ratio);

    // Effective Vc (with all effects)
    Vc_effective[i] = Vc * exp(eta_drug_Vc[drug[i]] + eta_subject_Vc[subject[i]] + eta_drug_Vc_within[subject[i], drug[i]] + beta_Vc_WT * log_wt_ratio);

    // Effective Vp (with all effects, including within Vp and weight)
    // Combines population, between-drug, between-subject, within-subject/drug, and covariate effects
    Vp_effective[i] = Vp * exp(eta_drug_Vp[drug[i]] + eta_subject_Vp[subject[i]] + eta_drug_Vp_within[subject[i], drug[i]] + beta_Vp_WT * log_wt_ratio);


    // Effective Q and ka (with independent random effects, assuming no weight effect for now)
    Q_effective[i] = Q * exp(eta_drug_Q[drug[i]] + eta_subject_Q[subject[i]]); // Add beta_Q_WT if needed
    ka_effective[i] = ka * exp(eta_drug_ka[drug[i]] + eta_subject_ka[subject[i]]); // Add beta_ka_WT if needed
  }

  // --- Implement 2-Compartment Prediction using ODE Solver ---

  // Solve ODEs for each observation point.
  // This approach solves from t=0 for each observation, handling dose at t=0
  // It's a simplification; a more accurate approach solves per subject with dose events at actual times.

  for (i in 1:N) {
     // Parameters for the ODE system: [CL, Vc, Vp, Q, ka]
     real theta[5] = {CL_effective[i], Vc_effective[i], Vp_effective[i], Q_effective[i], ka_effective[i]};

     // Integer data for ODE (not used by ODE function, but required by signature)
     int x_i[1] = {0};
      // Real data for ODE (not used by ODE function, but required by signature)
     real x_r[1] = {0.0};


     // Initial amounts at time 0: [Depot, Central, Peripheral]
     // Add dose instantly to the target compartment at t=0 based on route and F for SC
     real y0[3] = {0.0, 0.0, 0.0};
     if (route[i] == 0) { // IV Bolus
       y0[2] = dose[i]; // Add dose to Central (y[2])
     } else { // SC
       y0[1] = dose[i] * F; // Add dose*F to Depot (y[1])
     }


     // Solve ODE from time 0 up to observation time time[i]
     // The solver returns amounts at the requested output times (only time[i] in this case)
     real ts[1]; // Array of output times for the solver
     ts[1] = time[i];

     // amount_solution is a matrix: rows correspond to output times, columns to compartments
     matrix[rows(ts), cols(y0)] amount_solution =
        integrate_ode_bdf(two_cmt_ode, y0, 0.0, ts, theta, x_r, x_i,
                           rel_tol, abs_tol, max_steps);

     // Extract amount in Central compartment at time[i]
     // If time[i] == 0, amount_solution will contain y0
     real amount_central_at_obs_time = amount_solution[1, 2];

     // Calculate concentration (Mass/Volume)
     // Need to handle potential division by zero if Vc_effective is extremely small (should be positive)
     pred_conc[i] = amount_central_at_obs_time / Vc_effective[i];

     // Check for non-positive or NaN predictions from ODE solver output
     // This check helps identify problematic parameter values that lead to bad predictions
     // A check in the model block likelihood is also recommended.
     if (!is_finite(pred_conc[i])) {
        target += -1e18; // Penalize non-finite predictions
     } else if (pred_conc[i] < 0) {
        target += -1e18; // Penalize negative predictions
     }
  }
  // --- END OF ODE IMPLEMENTATION ---
}

model {
  // Priors for typical log parameters (Normal, centered at log(median))
  log_CL ~ normal(log(0.2), 1); // Prior for log of typical CL (elimination), median 0.2 L/time
  log_Vc ~ normal(log(4), 1); // Prior for log of typical Vc (central), median 4 L
  log_Vp ~ normal(log(10), 1); // Example Prior for log of typical Vp, median 10L
  log_Q ~ normal(log(5), 1); // Example Prior for log of typical Q, median 5 L/time
  log_ka ~ normal(log(0.4), 1); // Prior for log of typical ka, median 0.4 1/time

  // Prior for bioavailability
  F ~ beta(7, 3); // Example Beta prior centered around 0.7

  // Prior for residual error
  sigma ~ normal(0, 1); // Half-normal prior for sigma (due to lower=0 constraint)

  // Priors for covariance structure (Expanded to include CL, Vc, Vp correlation)
  L_drug_within_subject_cov ~ lkj_corr_cholesky(3); // Changed dimension
  // Priors for raw deviates (CL, Vc, Vp)
  for (s in 1:N_subjects) {
    to_vector(z_drug_subject[s]) ~ std_normal(); // Applies to the new [N_drugs, 3] matrix
  }


  // Priors for independent random effects (Q, ka)
  eta_drug_Q ~ normal(0, 0.5);
  eta_subject_Q ~ normal(0, 0.5);
  eta_drug_ka ~ normal(0, 0.5); // Prior for ka random effects
  eta_subject_ka ~ normal(0, 0.5); // Prior for ka random effects

  // Priors for between-drug and between-subject random effects (CL, Vc, Vp)
  eta_drug_CL ~ normal(0, 0.5);
  eta_subject_CL ~ normal(0, 0.5);
  eta_drug_Vc ~ normal(0, 0.5);
  eta_subject_Vc ~ normal(0, 0.5);
  eta_drug_Vp ~ normal(0, 0.5); // Prior for between-drug Vp
  eta_subject_Vp ~ normal(0, 0.5); // Prior for between-subject Vp


  // Priors for covariate coefficients (CL, Vc, Vp)
  beta_CL_WT ~ normal(0, 0.5);
  beta_Vc_WT ~ normal(0, 0.5);
  beta_Vp_WT ~ normal(0, 0.5); // Prior for weight effect on Vp


  // Likelihood for ALL data (using combined prediction vector)
  for (i in 1:N) {
    if (time[i] > 0) { // Only include observations after time 0 in likelihood
       // Check if prediction is positive before taking log
       // A robust ODE solver implementation should ideally not return non-positive for positive time.
       // This check penalizes problematic predictions during sampling.
       if (!is_finite(pred_conc[i]) || pred_conc[i] <= 0) {
         target += -1e18; // Add a large penalty for non-positive prediction
       } else {
         conc[i] ~ lognormal(log(pred_conc[i]), sigma);
       }
    }
  }
}

generated quantities {
  // Declare vectors and arrays for output
  vector[N] pred_conc_subj_drug; // Combined posterior predictive concentrations
  int subject_output[N]; // Subject indices for combined output
  int drug_output[N];   // Drug indices for combined output
  int route_output[N]; // Route for combined output
  vector[N] time_output; // Time for combined output
  vector[N] dose_output; // Dose for combined output
  vector[N] conc_output; // Observed concentration for combined output

  // Effective parameters for output
  vector[N] effective_CL_output = CL_effective;
  vector[N] effective_Vc_output = Vc_effective;
  vector[N] effective_Vp_output = Vp_effective;
  vector[N] effective_Q_output = Q_effective;
  vector[N] effective_ka_output = ka_effective;


  // Typical parameters for output
  real typical_CL = exp(log_CL);
  real typical_Vc = exp(log_Vc);
  real typical_Vp = exp(log_Vp);
  real typical_Q = exp(log_Q);
  real typical_ka = exp(log_ka);


  // Weight-related quantities (combined size)
  vector[N] obs_wt;
  vector[N] obs_log_wt_ratio;
  vector[N] obs_CL_wt_factor;
  vector[N] obs_Vc_wt_factor;
  vector[N] obs_Vp_wt_factor; // Added weight factor for Vp


  // Loop to populate weight-related outputs
  for (i in 1:N) {
    real current_subject_wt = subject_wt[subject[i]];
    real log_wt_ratio = log(current_subject_wt) - log(MEDWT);

    obs_wt[i] = current_subject_wt;
    obs_log_wt_ratio[i] = log_wt_ratio;
    obs_CL_wt_factor[i] = exp(beta_CL_WT * log_wt_ratio);
    obs_Vc_wt_factor[i] = exp(beta_Vc_WT * log_wt_ratio);
    obs_Vp_wt_factor[i] = exp(beta_Vp_WT * log_wt_ratio); // Calculate Vp weight factor
  }

  // Generate posterior predictive samples and store original data/indices
  for (i in 1:N) {
    subject_output[i] = subject[i];
    drug_output[i] = drug[i];
    route_output[i] = route[i];
    time_output[i] = time[i];
    dose_output[i] = dose[i];
    conc_output[i] = conc[i]; // Store observed concentration

    if (time[i] > 0) {
      // Use the pred_conc calculated in transformed parameters
       if (!is_finite(pred_conc[i]) || pred_conc[i] <= 0) {
         pred_conc_subj_drug[i] = -99.0; // Indicate problematic prediction
       } else {
         pred_conc_subj_drug[i] = lognormal_rng(log(pred_conc[i]), sigma);
       }
    } else {
      // Indicator for time <= 0 points
      pred_conc_subj_drug[i] = -99.0;
    }
  }
}
```

> # Compile and fit the Stan model
> compiled_model  <- cmdstan_model(write_stan_file(stan_model_code))

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable "real" does not exist.
 error in '/tmp/Rtmp1s9jms/model-3fc52455e4d9.stan' at line 196, column 9
  -------------------------------------------------
   194:      // Solve ODE from time 0 up to observation time time[i]
   195:      // The solver returns amounts at the requested output times (only time[i] in this case)
   196:      real ts[1]; // Array of output times for the solver
                ^
   197:      ts[1] = time[i];
  -------------------------------------------------

There’s a couple places where you define reals with length 1. Reals don’t need length definitions, you can just write real variable;. Same thing with integer arrays — instead of int variable[N];, you should define them as array[N] int variable;.

real variable also giving error.

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable "real" does not exist.
 error in '/tmp/Rtmpo0u4B5/model-3b2d7198785b.stan' at line 196, column 9
  -------------------------------------------------
   194:      // Solve ODE from time 0 up to observation time time[i]
   195:      // The solver returns amounts at the requested output times (only time[i] in this case)
   196:      real ts; // Array of output times for the solver.  
                ^
   197:      ts[1] = time[i];
  -------------------------------------------------

did not work for array;;SYNTAX ERROR, MESSAGE(S) FROM PARSER:

 error in '/tmp/RtmpNu5DJA/model-71a071891d9.stan' at line 196, column 13
  -------------------------------------------------
   194:      // Solve ODE from time 0 up to observation time time[i]
   195:      // The solver returns amounts at the requested output times (only time[i] in this case)
   196:         array[1] real ts; //real ts[1];   
                    ^
   197:         ts[1] = time[i]; 
  -------------------------------------------------```

vector[1] ts; also failed.

I believe the error occurs because variable ts is being declared in the middle of your model block, surrounded by procedural code. You would have to move the variable declaration to the beginning of the model block, before any assignment or probability statements.

Moved to beginning. did not work. Here is an error on ode.


stan_model_code <- "
// Define functions for the 2-compartment ODE model
functions {
  // Function defining the ODE system for amounts in compartments
  // y[1]: amount in depot (for SC)
  // y[2]: amount in central (Vc)
  // y[3]: amount in peripheral (Vp)

  // theta: parameters [CL, Vc, Vp, Q, ka]
  // x_r: real data (not used in this ODE function, but required by signature)
  // x_i: integer data (not used in this ODE function, but required by signature)
  // NOTE: x_r and x_i must be declared as 'data' in the function signature for ODE solvers
  // NOTE: Changing 'vector y' to 'real[] y' and return type to 'real[]' based on parser error message
  real[] two_cmt_ode(real t, real[] y, real[] theta, data real[] x_r, data int[] x_i) {
    real current_CL = theta[1];
    real current_Vc = theta[2];
    real current_Vp = theta[3];
    real current_Q = theta[4];
    real current_ka = theta[5]; // Only used if depot has amount

    // Add a small floor to volumes to prevent division by near-zero during integration
    real safe_Vc = fmax(current_Vc, 1e-12); // Use 1e-12 or a similarly small positive number
    real safe_Vp = fmax(current_Vp, 1e-12);

    // Calculate concentrations in central and peripheral compartments
    real conc_c = y[2] / safe_Vc; // Amount in central / Safe Volume of central
    real conc_p = y[3] / safe_Vp; // Amount in peripheral / Safe Volume of peripheral

    // Calculate rates of transfer and elimination
    real rate_cl = current_CL * conc_c; // Elimination from central
    real rate_q_out = current_Q * conc_c; // Transfer rate from central to peripheral
    real rate_q_in = current_Q * conc_p;    // Transfer rate from peripheral to central
    real rate_abs = current_ka * y[1];     // Absorption rate from depot

    // Rates of change for [Depot, Central, Peripheral]
    // Changed declaration to real[]
    real dydt[3]; // Correct syntax for declaring a sized real array

    // Rate of change for Depot (only applies if there's drug in the depot)
    dydt[1] = -rate_abs; // Drug leaves the depot by absorption

    // Rate of change for Central compartment
    dydt[2] = rate_abs       // Absorption into central (if SC)
              + rate_q_in    // Transfer from peripheral to central
              - rate_q_out // Transfer from central to peripheral
              - rate_cl;    // Elimination from central

    // Rate of change for Peripheral compartment
    dydt[3] = rate_q_out - rate_q_in; // Transfer between central and peripheral

    return dydt; // Return as real[]
  }
}

data {
  int<lower=0> N; // Total number of observations (nIV + nSC)
  int<lower=1> N_subjects; // Number of subjects
  int<lower=1> N_drugs; // Number of drugs

  int<lower=1, upper=N_subjects> subject[N]; // Subject index for each observation
  int<lower=1, upper=N_drugs> drug[N]; // Drug index for each observation
  int<lower=0, upper=1> route[N]; // Route (0 for IV, 1 for SC)
  vector<lower=0>[N] time; // Time points for all observations
  vector<lower=0>[N] conc; // Concentrations for all observations
  vector<lower=0>[N] dose; // Dose for each observation (0 for non-dose points, actual dose for dose points)

  vector<lower=0>[N_subjects] subject_wt; // Subject weight
  real<lower=0> MEDWT; // <<-- Median weight (e.g., 70 kg)

  // Dummy variables required by ODE solver signature - MUST be declared in data block
  // These must be provided as input data when fitting the model.
  real x_r[1];
  int x_i[1];

  // Remove ODE solver control parameters from data block
  // They are moved to transformed data
}
transformed data {
  // Declare and assign ODE solver control parameters here as fixed constants
  real rel_tol = 1e-6;
  real abs_tol = 1e-6;
  int max_steps = 10000000; // Assign as integer literal (1e7)

  // Removed assignments for x_r and x_i.
  // These must be provided in the data list when calling the model from R/Python.
}
parameters {
  // Typical population parameters (log scale) for standard 2-cmt model
  real log_CL; // Log of Typical elimination clearance
  real log_Vc; // Log of Typical central volume
  real log_Vp; // Log of Typical peripheral volume
  real log_Q; // Log of Typical intercompartmental clearance
  real log_ka; // Log of Absorption rate constant for SC

  real<lower=0,upper=1> F; // Bioavailability for SC

  real<lower=0,upper=100> sigma; // Residual error standard deviation (on log scale)

  // Random effects (on log scale)
  vector[N_drugs] eta_drug_CL;
  vector[N_subjects] eta_subject_CL;
  vector[N_drugs] eta_drug_Vc;
  vector[N_subjects] eta_subject_Vc;
  // Keep independent between-drug and between-subject Vp random effects
  vector[N_drugs] eta_drug_Vp;
  vector[N_subjects] eta_subject_Vp;
  vector[N_drugs] eta_drug_Q; // Between-drug Q rand effects
  vector[N_subjects] eta_subject_Q; // Between-subject Q rand effects
  vector[N_drugs] eta_drug_ka; // Between-drug ka rand effects
  vector[N_subjects] eta_subject_ka; // Between-subject ka rand effects


  // Correlation structure (Expanded to include CL, Vc, Vp correlation)
  cholesky_factor_cov[3] L_drug_within_subject_cov; // Changed from [2] to [3] for CL, Vc, Vp
  matrix[N_drugs, 3] z_drug_subject[N_subjects]; // Changed from [N_drugs, 2] to [N_drugs, 3]

  // Covariate coefficients (on log scale) - Assuming weight on CL, Vc, and Vp
  real beta_CL_WT;
  real beta_Vc_WT;
  real beta_Vp_WT; // Uncommented and added beta for Vp
  // Add beta_Q_WT, beta_ka_WT priors if needed
}

transformed parameters {
  // Population parameters (original scale)
  real CL = exp(log_CL);
  real Vc = exp(log_Vc);
  real Vp = exp(log_Vp);
  real Q = exp(log_Q);
  real ka = exp(log_ka);
  real ts[1]; // Declare ts outside the loop
  // Changed amount_solution from matrix to real array to match ODE solver return type
  real amount_solution[1, 3];
  // Declare amount_central_at_obs_time outside the loop
  real amount_central_at_obs_time;


  // Effective parameters for each observation
  vector[N] CL_effective;
  vector[N] Vc_effective;
  vector[N] Vp_effective;
  vector[N] Q_effective;
  vector[N] ka_effective;

  // Combined predictions for all observations
  vector[N] pred_conc;

  // Transform within-subject/drug deviates for CL, Vc, and Vp
  vector[N_drugs] eta_drug_CL_within[N_subjects];
  vector[N_drugs] eta_drug_Vc_within[N_subjects];
  vector[N_drugs] eta_drug_Vp_within[N_subjects]; // Added within Vp


  for (s in 1:N_subjects) {
    for (d in 1:N_drugs) {
      // Joint deviation vector now has size 3
      vector[3] eta_drug_joint = L_drug_within_subject_cov * z_drug_subject[s, d]';
      eta_drug_CL_within[s, d] = eta_drug_joint[1];
      eta_drug_Vc_within[s, d] = eta_drug_joint[2];
      eta_drug_Vp_within[s, d] = eta_drug_joint[3]; // Extract Vp component
    }
  }

  // Calculate effective parameters for each observation
  for (i in 1:N) {
    real current_subject_wt = subject_wt[subject[i]];
    real log_wt_ratio = log(current_subject_wt) - log(MEDWT);

    // Effective CL (with all effects)
    CL_effective[i] = CL * exp(eta_drug_CL[drug[i]] + eta_subject_CL[subject[i]] + eta_drug_CL_within[subject[i], drug[i]] + beta_CL_WT * log_wt_ratio);

    // Effective Vc (with all effects)
    Vc_effective[i] = Vc * exp(eta_drug_Vc[drug[i]] + eta_subject_Vc[subject[i]] + eta_drug_Vc_within[subject[i], drug[i]] + beta_Vc_WT * log_wt_ratio);

    // Effective Vp (with all effects, including within Vp and weight)
    // Combines population, between-drug, between-subject, within-subject/drug, and covariate effects
    Vp_effective[i] = Vp * exp(eta_drug_Vp[drug[i]] + eta_subject_Vp[subject[i]] + eta_drug_Vp_within[subject[i], drug[i]] + beta_Vp_WT * log_wt_ratio);


    // Effective Q and ka (with independent random effects, assuming no weight effect for now)
    Q_effective[i] = Q * exp(eta_drug_Q[drug[i]] + eta_subject_Q[subject[i]]); // Add beta_Q_WT if needed
    ka_effective[i] = ka * exp(eta_drug_ka[drug[i]] + eta_subject_ka[subject[i]]); // Add beta_ka_WT if needed
  }

  // --- Implement 2-Compartment Prediction using ODE Solver ---

  // Solve ODEs for each observation point.
  // This approach solves from t=0 for each observation, handling dose at t=0
  // It's a simplification; a more accurate approach solves per subject with dose events at actual times.

  for (i in 1:N) {
    // Parameters for the ODE system: [CL, Vc, Vp, Q, ka]
    // Explicitly declare theta as real[5]
    real theta[5] = {CL_effective[i], Vc_effective[i], Vp_effective[i], Q_effective[i], ka_effective[i]};

    // Initial amounts at time 0: [Depot, Central, Peripheral]
    // y0 remains real[] as required by the ODE solver's second argument type
    real y0[3] = {0.0, 0.0, 0.0};
    if (route[i] == 0) { // IV Bolus
      y0[2] = dose[i]; // Add dose to Central (y[2])
    } else { // SC
      y0[1] = dose[i] * F; // Add dose*F to Depot (y[1])
    }

    // Solve ODE from time 0 up to observation time time[i]
    // The solver returns amounts at the requested output times (only time[i] in this case)
    ts[1] = time[i];

    // amount_solution is a real array: rows correspond to output times, columns to compartments
    amount_solution =
        integrate_ode_bdf(two_cmt_ode, y0, 0.0, ts, theta, x_r, x_i,
                           rel_tol, abs_tol, max_steps);

    // Extract amount in Central compartment at time[i]
    // If time[i] == 0, amount_solution will contain y0
    // Separated declaration and assignment
    amount_central_at_obs_time = amount_solution[1, 2];

    // Calculate concentration (Mass/Volume)
    // Need to handle potential division by zero if Vc_effective is extremely small (should be positive)
    pred_conc[i] = amount_central_at_obs_time / Vc_effective[i];

    // Check for non-positive or NaN predictions from ODE solver output
    // This check helps identify problematic parameter values that lead to bad predictions
    // A check in the model block likelihood is also recommended.
    if (!is_finite(pred_conc[i])) {
        target += -1e18; // Penalize non-finite predictions
    } else if (pred_conc[i] < 0) {
        target += -1e18; // Penalize negative predictions
    }
  }
  // --- END OF ODE IMPLEMENTATION ---
}

model {
  // Priors for typical log parameters (Normal, centered at log(median))
  log_CL ~ normal(log(0.2), 1); // Prior for log of typical CL (elimination), median 0.2 L/time
  log_Vc ~ normal(log(4), 1); // Prior for log of typical Vc (central), median 4 L
  log_Vp ~ normal(log(10), 1); // Example Prior for log of typical Vp, median 10L
  log_Q ~ normal(log(5), 1); // Example Prior for log of typical Q, median 5 L/time
  log_ka ~ normal(log(0.4), 1); // Prior for log of typical ka, median 0.4 1/time

  // Prior for bioavailability
  F ~ beta(7, 3); // Example Beta prior centered around 0.7

  // Prior for residual error
  sigma ~ normal(0, 1); // Half-normal prior for sigma (due to lower=0 constraint)

  // Priors for covariance structure (Expanded to include CL, Vc, Vp correlation)
  L_drug_within_subject_cov ~ lkj_corr_cholesky(3); // Changed dimension
  // Priors for raw deviates (CL, Vc, Vp)
  for (s in 1:N_subjects) {
    to_vector(z_drug_subject[s]) ~ std_normal(); // Applies to the new [N_drugs, 3] matrix
  }


  // Priors for independent random effects (Q, ka)
  eta_drug_Q ~ normal(0, 0.5);
  eta_subject_Q ~ normal(0, 0.5);
  eta_drug_ka ~ normal(0, 0.5); // Prior for ka random effects
  eta_subject_ka ~ normal(0, 0.5); // Prior for ka random effects

  // Priors for between-drug and between-subject random effects (CL, Vc, Vp)
  eta_drug_CL ~ normal(0, 0.5);
  eta_subject_CL ~ normal(0, 0.5);
  eta_drug_Vc ~ normal(0, 0.5);
  eta_subject_Vc ~ normal(0, 0.5);
  eta_drug_Vp ~ normal(0, 0.5); // Prior for between-drug Vp
  eta_subject_Vp ~ normal(0, 0.5); // Prior for between-subject Vp


  // Priors for covariate coefficients (CL, Vc, Vp)
  beta_CL_WT ~ normal(0, 0.5);
  beta_Vc_WT ~ normal(0, 0.5);
  beta_Vp_WT ~ normal(0, 0.5); // Prior for weight effect on Vp


  // Likelihood for ALL data (using combined prediction vector)
  for (i in 1:N) {
    if (time[i] > 0) { // Only include observations after time 0 in likelihood
      // Check if prediction is positive before taking log
      // A robust ODE solver implementation should ideally not return non-positive for positive time.
      // This check penalizes problematic predictions during sampling.
      if (!is_finite(pred_conc[i]) || pred_conc[i] <= 0) {
        target += -1e18; // Add a large penalty for non-positive prediction
      } else {
        conc[i] ~ lognormal(log(pred_conc[i]), sigma);
      }
    }
  }
}

generated quantities {
  // Declare vectors and arrays for output
  vector[N] pred_conc_subj_drug; // Combined posterior predictive concentrations
  int subject_output[N]; // Subject indices for combined output
  int drug_output[N];    // Drug indices for combined output
  int route_output[N]; // Route for combined output
  vector[N] time_output; // Time for combined output
  vector[N] dose_output; // Dose for combined output
  vector[N] conc_output; // Observed concentration for combined output

  // Effective parameters for output
  vector[N] effective_CL_output = CL_effective;
  vector[N] effective_Vc_output = Vc_effective;
  vector[N] effective_Vp_output = Vp_effective;
  vector[N] effective_Q_output = Q_effective;
  vector[N] effective_ka_output = ka_effective;


  // Typical parameters for output
  real typical_CL = exp(log_CL);
  real typical_Vc = exp(log_Vc);
  real typical_Vp = exp(log_Vp);
  real typical_Q = exp(log_Q);
  real typical_ka = exp(log_ka);


  // Weight-related quantities (combined size)
  vector[N] obs_wt;
  vector[N] obs_log_wt_ratio;
  vector[N] obs_CL_wt_factor;
  vector[N] obs_Vc_wt_factor;
  vector[N] obs_Vp_wt_factor; // Added weight factor for Vp


  // Loop to populate weight-related outputs
  for (i in 1:N) {
    real current_subject_wt = subject_wt[subject[i]];
    real log_wt_ratio = log(current_subject_wt) - log(MEDWT);

    obs_wt[i] = current_subject_wt;
    obs_log_wt_ratio[i] = log_wt_ratio;
    obs_CL_wt_factor[i] = exp(beta_CL_WT * log_wt_ratio);
    obs_Vc_wt_factor[i] = exp(beta_Vc_WT * log_wt_ratio);
    obs_Vp_wt_factor[i] = exp(beta_Vp_WT * log_wt_ratio); // Calculate Vp weight factor
  }

  // Generate posterior predictive samples and store original data/indices
  for (i in 1:N) {
    subject_output[i] = subject[i];
    drug_output[i] = drug[i];
    route_output[i] = route[i];
    time_output[i] = time[i];
    dose_output[i] = dose[i];
    conc_output[i] = conc[i]; // Store observed concentration

    if (time[i] > 0) {
      // Use the pred_conc calculated in transformed parameters
      if (!is_finite(pred_conc[i]) || pred_conc[i] <= 0) {
        pred_conc_subj_drug[i] = -99.0; // Indicate problematic prediction
      } else {
        pred_conc_subj_drug[i] = lognormal_rng(log(pred_conc[i]), sigma);
      }
    } else {
      // Indicator for time <= 0 points
      pred_conc_subj_drug[i] = -99.0;
    }
  }
}


"


Dimension mismatch in assignment; variable name = amount_solution, type = matrix; right-hand side type = real[ , ].
Illegal statement beginning with non-void expression parsed as
  amount_solution
Not a legal assignment, sampling, or function statement.  Note that
  * Assignment statements only allow variables (with optional indexes) on the left;
  * Sampling statements allow arbitrary value-denoting expressions on the left.
  * Functions used as statements must be declared to have void returns

 error in '/tmp/Rtmpn4bC92/model-401c77fde23b.stan' at line 207, column 4
  -------------------------------------------------
   205: 
   206:     // amount_solution is a matrix: rows correspond to output times, columns to compartments
   207:     amount_solution =
           ^
   208:         integrate_ode_bdf(two_cmt_ode, y0, 0.0, ts, theta, x_r, x_i,
  -------------------------------------------------

PARSER EXPECTED: "}"


make: *** [/tmp/Rtmpn4bC92/model-401c77fde23b.hpp] Error 253

Error: An error occured during compilation! See the message above for more information.
Execution halted```

It looks like you’re using an older version of CmdStan in cmdstanr, because we deprecated for more than a year and then finally removed the array syntax you are using. Our current interfaces won’t parse it at all. I would suggest updating your versions of CmdStan and cmdstanr so that they match all of our current doc.

Everywhere you have something like

real a[5];

you need to replace with

array[5] real a;

In function arguments, that goes from real[,] to array[,] real.

We did this so that our types were contiguous, which is required by the tuple syntax.

Thank you for the suggestion . I am using cmdstan-2.21.0, cmdstanr version 0.6.1 which are older versions, causing these issues. I will update them. Thanks.