Torsten - isize

Checked, all groups in the data passed to Stan have at least one event, ; Run still fails with this message; Tried with & without initial value assignment. ;;;;; Any help will be appreciated.

Chain 1 Rejecting initial value:
Chain 1   Error evaluating the log probability at the initial value.
Chain 1 Exception: EventHistory::EventHistory: isize is 0, but must be greater than or equal to 1.000000 (in '/tmp/RtmpDeAgU3/model-3be6d2c5a47.stan', line 296, column 1 to line 316, column 2)
functions {
+  // ODE system for a 2-compartment model with a depot
+ // y[1]: amount in depot
+  // y[2]: amount in central
+  // y[3]: amount in peripheral
+  // theta[1]: CL
+  // theta[2]: Vc
+  // theta[3]: Vp
+ // theta[4]: Q
+  // theta[5]: ka
+  // Note: x_r and x_i are dummy variables required by Stan's ODE solver signature
+  vector two_cmt_ode_torsten(real t, vector y, array[] real theta, array[] real x_r, array[] int x_i) {
+  real CL = theta[1];
+ real Vc = theta[2];
+  real Vp = theta[3];
+  real Q = theta[4];
+  real ka = theta[5];
+ 
+  // Ensure volumes are positive and not near zero for division
+  real safe_Vc = fmax(Vc, 1e-12);
+  real safe_Vp = fmax(Vp, 1e-12);
+ 
+  // Micro-rate constants
+  real k10 = CL / safe_Vc;
+  real k12 = Q / safe_Vc;
+  real k21 = Q / safe_Vp;
+ 
+  // Rates of change for [Depot, Central, Peripheral]
+  vector[3] dydt;
+ 
+  dydt[1] = -ka * y[1]; // Absorption from depot
+ 
+ dydt[2] = ka * y[1] // Absorption into central
+  + k21 * y[3] // From peripheral
+  - k12 * y[2] // To peripheral
+  - k10 * y[2]; // Elimination
+ 
+  dydt[3] = k12 * y[2] // From central
+  - k21 * y[3]; // To central
+ 
+  return dydt;
+  }
+ }
+ data {
+  // Declare N_subjects and N_drugs FIRST, as they are used in subsequent declarations
+  int<lower=1> N_subjects; // Number of subjects
+  int<lower=1> N_drugs; // Number of drugs
+  vector<lower=0>[N_subjects] subject_wt; // Subject weight
+  real<lower=0> MEDWT; // Median weight (e.g., 70 kg)
+ 
+  // Consolidated Event Data Structure for Torsten
+  int<lower=0> N_events; // Total number of events (doses + observations)
+  array[N_events] real event_time; // Corrected array declaration
+  array[N_events] real event_amt; // Corrected array declaration
+  // Use standard array declarations: type variable_name[dimension];
+  array[N_events] int event_cmt; // Corrected array declaration // Compartment for dose input (1=Depot, 2=Central) or observation (e.g., 2=Central)
+  array[N_events] int<lower=0, upper=1> event_evid; // Corrected array declaration // Event ID (1=Dose, 0=Observation)
+  // The event_id passed to pmx_solve_group MUST be the mapped group ID (1 to N_subjects*N_drugs)
+  array[N_events] int event_id; // Corrected array declaration // Constraint removed based on previous error
+ 
+  // event_drug_map and event_subject_map are needed to apply drug/subject random effects and covariates
+  // These indices should map back to the rows of eta_drug_*, eta_subject_*, z_drug_subject, and subject_wt
+  // Let's pass the mapped subject and drug IDs corresponding to each event - Use standard array declarations
+  array[N_events] int<lower=1, upper=N_subjects> event_subject_map; // Corrected array declaration // Mapped Subject ID (1 to N_subjects) for each event
+  array[N_events] int<lower=1, upper=N_drugs> event_drug_map; // Corrected array declaration // Mapped Drug ID (1 to N_drugs) for each event
+ 
+  array[N_events] real event_rate; // Corrected array declaration // Rate (0 for bolus/observation, 1 for SC )
+ 
+  int<lower=0> N_obs; // Total number of observations (EVID=0)
+  // conc_obs now contains the original observed values (LLOQ/2 if BLQ, or measured value if >LLOQ)
+  vector<lower=0>[N_obs] conc_obs;
+  array[N_obs] int<lower=1, upper=N_events> obs_event_idx; // Corrected array declaration // Index in event_time/... for each observation
+ 
+  // Additional data variables passed from R - Use standard array declarations
+  array[N_events] int event_subject_id_orig; // Corrected array declaration // Original Subject ID for each event (for output)
+  array[N_events] int event_drug_id_orig; // Corrected array declaration // Original Drug ID for each event (for output)
+  // event_blq flag is now used in the likelihood
+  array[N_events] int event_blq; // Corrected array declaration // BLQ flag for each event
+  array[N_events] int event_mdv; // Corrected array declaration // MDV flag for each event (already filtered in R prep)
+ 
+  // LLOQ value for censored likelihood # ADD THIS
+  real<lower=0> LLOQ;
+ 
+  // Dummy variables required by ODE solver signature (dim=0 if not used by ODE system function theta)
+  // CORRECTED: Declare x_r and x_i as 2D arrays of size 0x0
+  array[0, 0] real x_r; // Dummy real data needed for ODE solver signature
+  array[0, 0] int x_i; // Dummy integer data needed for ODE solver signature
+ 
+  // Torsten ODE solver controls (can also be in transformed data)
+  real rel_tol; // Relative tolerance for ODE solver
+  real abs_tol; // Absolute tolerance for ODE solver
+  int max_steps; // Maximum steps for ODE solver
+ 
+  // Priors for log_sigma (define here or directly in model block)
+  real mu_log_sigma; // Mean for the normal prior on log_sigma (e.g., log(typical_sigma))
+  real sigma_log_sigma; // Std dev for the normal prior on log_sigma (e.5 or 1)
+ }
+ 
+ 
+ transformed data {
+   int nGroups = N_subjects * N_drugs;
+  int nTheta = 5; // number of parameters
+  int nCmt = 3; // number of compartments
+  array[nGroups] int len = rep_array(0, nGroups); // Initialize len with zeros
+  array[nGroups] int len_theta;
+ array[nGroups] int len_biovar;
+  array[nGroups] int len_tlag;
+  array[N_events] int ii;
+  array[N_events] int addl;
+  array[N_events] int ss;
+ 
+  // Robust calculation of len by counting events per group ID
+  for (i in 1:N_events) {
+  if (event_id[i] >= 1 && event_id[i] <= nGroups) {
+  len[event_id[i]] += 1;
+  } else {
+  // Handle invalid event_id if necessary (e.g., print warning or error)
+  reject(\"Invalid event_id encountered: \", event_id[i], \" at event index \", i);
+  }
+  }
+ 
+  // len_theta, len_biovar, len_tlag should match len for pmx_solve_group_bdf
+  for (i in 1:nGroups) {
+  len_theta[i] = len[i];
+  len_biovar[i] = len[i];
+ len_tlag[i] = len[i];
+  }
+ 
+  // Debugging: Print the sum of the len array
+  int sum_len = 0;
+  for (i in 1:nGroups) {
+  sum_len += len[i];
+  }
+  print(\"Sum of len: \", sum_len); // Escaped quotes
+  print(\"N_events: \", N_events); // Escaped quotes
+ 
+  // Check if sum of len matches N_events
+  if (sum_len != N_events) {
+  reject(\"Sum of len (\", sum_len, \") does not match N_events (\", N_events, \").\");
+  }
+ 
+ 
+  for (i in 1:N_events) {
+  ii[i] = 0;// No inter-dose interval for a single dose
+  addl[i] = 0;// No additional doses for a single dose
+  ss[i] = 0;// Not steady state for a single dose
+  }
+ 
+ }
+ 
+ parameters {
+  // Residual error standard deviation on the log scale
+  real log_sigma;
+ 
+  // 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
+ 
+  // Population typical bioavailability (logit scale)
+  real logit_F_pop;
+ 
+  // Random effects (on log scale, except for F which is logit 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;
+  vector[N_drugs] eta_drug_Vp;
+  vector[N_subjects] eta_subject_Vp;
+  vector[N_drugs] eta_drug_Q;
+  vector[N_subjects] eta_subject_Q;
+  vector[N_drugs] eta_drug_ka;
+  vector[N_subjects] eta_subject_ka;
+  vector[N_drugs] eta_drug_F; // Drug random effects for F (logit scale)
+ 
+  // Standard deviations for random effects (Half-Normal or similar)
+  real<lower=0> omega_drug_CL;
+  real<lower=0> omega_subject_CL;
+  real<lower=0> omega_drug_Vc;
+  real<lower=0> omega_subject_Vc;
+  real<lower=0> omega_drug_Vp;
+  real<lower=0> omega_subject_Vp;
+  real<lower=0> omega_drug_Q;
+  real<lower=0> omega_subject_Q;
+  real<lower=0> omega_drug_ka;
+  real<lower=0> omega_subject_ka;
+  real<lower=0> omega_drug_F; // Std dev for eta_drug_F
+ 
+  // Correlation structure for CL, Vc, Vp correlation)
+  cholesky_factor_cov[3] L_drug_within_subject_cov;
+  array[N_subjects, N_drugs] matrix[1, 3] z_drug_subject; // Corrected array of matrices declaration - Note: matrix[1,3] is likely what was intended as raw_deviates_sd is row_vector[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;
+  // Add beta_Q_WT, beta_ka_WT, beta_F_WT if needed
+ }
+ transformed parameters {
+  // Exponentiate log_sigma to get sigma on the original scale
+  real sigma = exp(log_sigma);
+ 
+  // Population parameters on original scale
+  real CL_pop = exp(log_CL);
+  real Vc_pop = exp(log_Vc);
+  real Vp_pop = exp(log_Vp);
+  real Q_pop = exp(log_Q);
+  real ka_pop = exp(log_ka);
+  real pop_F = inv_logit(logit_F_pop); // Population typical F on original scale
+ 
+  // Transform within-subject/drug deviates for CL, Vc, and Vp
+  array[N_subjects, N_drugs] real eta_drug_CL_within; // Corrected array declaration
+  array[N_subjects, N_drugs] real eta_drug_Vc_within; // Corrected array declaration
+  array[N_subjects, N_drugs] real eta_drug_Vp_within; // Corrected array declaration
+ 
+  for (s in 1:N_subjects) {
+  for (d in 1:N_drugs) {
+ row_vector[3] raw_deviates_sd = z_drug_subject[s, d][1]; // Access the 1x3 matrix
+  vector[3] eta_drug_joint = L_drug_within_subject_cov * raw_deviates_sd';
+  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];
+  }
+  }
+ 
+  // Calculate drug-specific F values on original scale
+  array[N_drugs] real F_drug; // Corrected array declaration
+  for (d in 1:N_drugs) {
+  // Combine population logit_F with drug-specific random effect
+  real logit_F_drug_d = logit_F_pop + eta_drug_F[d];
+  F_drug[d] = inv_logit(logit_F_drug_d);
+  }
+ 
+  // Individual parameters for each subject-drug combination (used by pmx_solve_group)
+  array[nGroups, nTheta] real theta; // Corrected array declaration
+ 
+  // Calculate individual parameters for each subject-drug group
+  for (s in 1:N_subjects) {
+  for (d in 1:N_drugs) {
+  real current_subject_wt = subject_wt[s];
+  real log_wt_ratio = log(current_subject_wt) - log(MEDWT);
+ 
+  // Calculate individual log parameters including all effects
+  real ind_log_CL = log(CL_pop) + eta_drug_CL[d] + eta_subject_CL[s] + eta_drug_CL_within[s, d] + beta_CL_WT * log_wt_ratio;
+ real ind_log_Vc = log(Vc_pop) + eta_drug_Vc[d] + eta_subject_Vc[s] + eta_drug_Vc_within[s, d] + beta_Vc_WT * log_wt_ratio;
+  real ind_log_Vp = log(Vp_pop) + eta_drug_Vp[d] + eta_subject_Vp[s] + eta_drug_Vp_within[s, d] + beta_Vp_WT * log_wt_ratio;
+  real ind_log_Q = log(Q_pop) + eta_drug_Q[d] + eta_subject_Q[s]; // Add beta_Q_WT if needed
+  real ind_log_ka = log(ka_pop) + eta_drug_ka[d] + eta_subject_ka[s]; // Add beta_ka_WT if needed
+ 
+ // Assign exponentiated individual parameters to the theta matrix for the group
+  theta[(s - 1) * N_drugs + d, 1] = exp(ind_log_CL); // CL
+  theta[(s - 1) * N_drugs + d, 2] = exp(ind_log_Vc); // Vc
+  theta[(s - 1) * N_drugs + d, 3] = exp(ind_log_Vp); // Vp
+  theta[(s - 1) * N_drugs + d, 4] = exp(ind_log_Q); // Q
+  theta[(s - 1) * N_drugs + d, 5] = exp(ind_log_ka); // ka
+  }
+  }
+ 
+ // Bioavailability matrix for each group and compartment
+  array[nGroups, nCmt] real biovar_group; // Renamed to avoid confusion
+  for (s in 1:N_subjects) {
+  for (d in 1:N_drugs) {
+  // Use the drug-specific F_drug value for the Depot (CMT 1)
+  biovar_group[(s - 1) * N_drugs + d, 1] = F_drug[d];
+  biovar_group[(s - 1) * N_drugs + d, 2] = 1.0; // Bioavailability for Central (CMT 2) is 1 (for IV)
+  biovar_group[(s - 1) * N_drugs + d, 3] = 1.0; // Bioavailability for Peripheral (CMT 3)
+  }
+  }
+  array[nGroups, nCmt] real tlag_group; // Renamed to avoid confusion
+  for (i in 1:nGroups){
+ for (j in 1:nCmt){
+ tlag_group[i,j] = 0;
+ }
+  }
+ 
+  // Create event-indexed bioavailability and tlag matrices for the ODE solver
+  array[N_events, nCmt] real biovar_event;
+  array[N_events, nCmt] real tlag_event;
+ 
+  for (i in 1:N_events) {
+  int current_group_id = event_id[i]; // Get the group ID for this event
+  for (j in 1:nCmt) {
+  biovar_event[i, j] = biovar_group[current_group_id, j]; // Assign group biovar to event
+  tlag_event[i, j] = tlag_group[current_group_id, j]; // Assign group tlag to event
+  }
+  }
+ 
+ 
+ 
+  // Solve ODE system for all events and all groups using pmx_solve_group_bdf
+  matrix[N_events, 3] amount_solution;
+ 
+  amount_solution = pmx_solve_group_bdf(
+  two_cmt_ode_torsten, // ODE system function
+  nCmt, // Number of compartments
+  len, // Length of event arrays for each group - CORRECTED
+  event_time, // Time of events
+   event_amt, // Amount of drug for events
+  event_rate, // Rate of infusion (if any)
+  ii, // Inter-dose interval
+  event_evid, // Event type (dose or observation)
+  event_cmt, // Compartment number
+  addl, // Number of additional doses
+ ss, // Steady-state indicator
+  theta, // Individual parameters
+  biovar_group, // Bioavailability matrix - CORRECTED to group-indexed
+  tlag_group, // Tlag matrix - CORRECTED to group-indexed
+  x_r, // Real data for ODE function (dummy - CORRECTED TO 2D)
+  x_i, // Integer data for ODE function (dummy - CORRECTED TO 2D)
+  rel_tol, // Relative tolerance for ODE solver
+  abs_tol, // Absolute tolerance for ODE solver
+ max_steps // Maximum steps for ODE solver
+ );
+ 
+ 
+  // Extract predicted concentrations at observation times (Central compartment, CMT=2)
+ array[N_obs] real pred_conc; // Corrected array declaration
+  array[nGroups] real Vc_effective_group; // Corrected array declaration // Store Vc for each group
+ 
+  // Store Vc for each group to use for concentration calculation
+  for (g in 1:nGroups) {
+  Vc_effective_group[g] = theta[g, 2]; // Vc is the 2nd parameter in theta matrix
+  }
+ 
+  for (i in 1:N_obs) {
+  // Get the event index corresponding to this observation
+  int current_event_idx = obs_event_idx[i];
+  // Need to determine the group ID for this observation event
+  int current_group_id = event_id[current_event_idx]; // Assuming event_id maps to the group
+ 
+  // Extract the amount in the central compartment at this observation event time
+  real amount_central_at_obs = amount_solution[current_event_idx, 2]; // Central is the 2nd compartment (y[2])
+ 
+  // Calculate concentration (Amount / Volume)
+  if (Vc_effective_group[current_group_id] > 1e-12) {
+  pred_conc[i] = amount_central_at_obs / Vc_effective_group[current_group_id];
+  } else {
+  pred_conc[i] = -1.0; // Indicate a problem if Vc is zero or negative
+  }
+  }
+ }
+ 
+ 
+ 
+ 
+ 
+ model {
+  // Priors for typical log parameters (Normal, centered at log(median))
+  log_CL ~ normal(log(0.2), 1); // Prior for log of typical CL
+  log_Vc ~ normal(log(3), 1); // Prior for log of typical Vc
+  log_Vp ~ normal(log(4), 1); // Prior for log of typical Vp
+  log_Q ~ normal(log(0.5), 1); // Prior for log of typical Q
+  log_ka ~ normal(log(0.4), 1); // Prior for log of typical ka
+ 
+  // Prior for population typical bioavailability (logit scale) - Center around a reasonable value like logit(0.7)
+  logit_F_pop ~ normal(logit(0.7), 1); // Prior for logit of typical F
+ 
+  // Prior for log_sigma (Normal prior on the log-scaled residual error standard deviation)
+  log_sigma ~ normal(mu_log_sigma, sigma_log_sigma); // Use data-defined prior parameters
+ 
+  // Priors for standard deviations for random effects (Half-Normal or similar)
+  omega_drug_CL ~ cauchy(0, 0.3); // Example Half-Cauchy prior
+  omega_subject_CL ~ cauchy(0, 0.3);
+  omega_drug_Vc ~ cauchy(0, 0.3);
+  omega_subject_Vc ~ cauchy(0, 0.3);
+  omega_drug_Vp ~ cauchy(0, 0.3);
+  omega_subject_Vp ~ cauchy(0, 0.3);
+  omega_drug_Q ~ cauchy(0, 0.3);
+  omega_subject_Q ~ cauchy(0, 0.3);
+  omega_drug_ka ~ cauchy(0, 0.3);
+  omega_subject_ka ~ cauchy(0, 0.3);
+  omega_drug_F ~ cauchy(0, 0.3); // Prior for std dev of eta_drug_F
+ 
+  // Priors for random effects (Normal, centered at 0, scaled by corresponding omega)
+  eta_drug_CL ~ normal(0, omega_drug_CL);
+  eta_subject_CL ~ normal(0, omega_subject_CL);
+  eta_drug_Vc ~ normal(0, omega_drug_Vc);
+  eta_subject_Vc ~ normal(0, omega_subject_Vc);
+  eta_drug_Vp ~ normal(0, omega_drug_Vp);
+  eta_subject_Vp ~ normal(0, omega_subject_Vp); // Corrected: Use omega_subject_Vp
+  eta_drug_Q ~ normal(0, omega_drug_Q);
+  eta_subject_Q ~ normal(0, omega_subject_Q);
+  eta_drug_ka ~ normal(0, omega_drug_ka);
+  eta_subject_ka ~ normal(0, omega_subject_ka);
+  eta_drug_F ~ normal(0, omega_drug_F);
+ 
+  // Priors for covariance structure (CL, Vc, Vp correlation)
+  L_drug_within_subject_cov ~ lkj_corr_cholesky(3); // Dimension 3 for CL, Vc, Vp
+ 
+  // Priors for raw deviates (CL, Vc, Vp) - applied per subject and drug
+  // Corrected: Iterate through subjects and drugs and apply prior to each matrix[1, 3]
+  for (s in 1:N_subjects) {
+  for (d in 1:N_drugs) {
+  to_vector(z_drug_subject[s, d]) ~ std_normal();
+  }
+  }
+ 
+  // Priors for covariate coefficients (CL, Vc, Vp)
+  beta_CL_WT ~ normal(0, 0.3);
+  beta_Vc_WT ~ normal(0, 0.3);
+  beta_Vp_WT ~ normal(0, 0.3);
+  // Add beta_Q_WT, beta_ka_WT, beta_F_WT if needed
+ 
+  // Likelihood for ALL observations (EVID = 0 and not MDV)
+ // Handles BLQ data using censored likelihood
+  for (i in 1:N_obs) {
+  // Get the event index corresponding to this observation
+  int current_event_idx = obs_event_idx[i];
+  // Get the BLQ status and original observed value for this observation row
+  int is_blq = event_blq[current_event_idx];
+  real observed_value = conc_obs[i]; // This holds the DVOR value (LLOQ/2 if BLQ, or measured value if >LLOQ)
+ 
+  real predicted_concentration = pred_conc[i]; // Prediction for this observation
+ 
+  // Ensure prediction is positive before using lognormal/censoring functions
+ if (predicted_concentration <= 0) {
+ // Apply large penalty for non-positive prediction - indicates model failure
+  target += -1e18;
+ } else {
+  if (is_blq == 0) { // Not a BLQ observation (original value > LLOQ)
+  // Standard lognormal likelihood for observed value > 0
+  // We use the original observed_value here (should be > 0)
+  if (observed_value > 0) {
+  observed_value ~ lognormal(log(predicted_concentration), sigma);
+ } else {
+  // This case suggests a data issue where a non-BLQ observation is <= 0.
+ // Apply a penalty for unexpected observed value.
+  target += -1e18;
+  }
+ 
+  } else { // is_blq == 1, it's a BLQ observation (original value was <= LLOQ, imputed as LLOQ/2)
+  // Censored lognormal likelihood: Probability that the true value is <= LLOQ
+  // P(Y <= LLOQ) where Y ~ Lognormal(log(pred), sigma)
+ // On log scale: P(log(Y) <= log(LLOQ)) where log(Y) ~ Normal(log(pred), sigma)
+  // Use the Normal CDF: Phi((log(LLOQ) - mu_log) / sigma_log)
+  // Stan's normal_lcdf computes log(P(X <= y)) for X ~ Normal(mu, sigma)
+  target += normal_lcdf(log(LLOQ) | log(predicted_concentration), sigma);
+  // Note: We do NOT use the imputed 'observed_value' (LLOq/2) here for thelikelihood.
+  // The imputed value in conc_obs for BLQ rows is not used in this branch.
+  }
+  }
+  }
+ }

Torsten’s not part of the Stan project. It’s owned and distributed by Metrum. Three folks on this list were involved in developing Torsten, and they may be able to help: @yizhang, @charlesm93, and @billg.

Thank you. Will contact them.

@somkamalin I can’t pick out exactly what the issue is without the data, but first I’d try to print() things. From looking at the C++ code in the Torsten repo, it prints that error when something in your len is 0 (which I think you’ve figured out from your comment that you checked that all groups have at least one event), but I’d still bet something in len is 0 where it has to be \geq 1. So do a print(len); on line 121 just after you fill it with values. You can also probably more easily create this vector in R and pass it in as data. Either way is fine, though. It’s worth doing it in R as a self-check, regardless.

Also, I’m not sure why you’re wanting to use an ODE for a two-compartment model other than I guess you want to be able to use Torsten’s group solver for within-chain parallelization, but to be honest, that’s probably still slower than no within-chain parallelization while using the analytical (or linear ODE) solution - the extra cost in solving the ODE vs. the analytical solution is more than the savings of within-chain parallelization. That being said, you probably still want to use within-chain parallelization, so as a shameless plug, check out this paper and this repo. It’ll help you use within-chain parallelization using multi-threading with reduce_sum(), which can use any of the analytical, linear ODE, or general ODE solutions (you want to use the analytical solution for your two-compartment depot model) and will provide pretty big speedups. You should be able to integrate your model with this model without too much trouble.

Casey

Thank you. Tried both anal.sol. and ODE and references suggested above. Still have issues, may be due to data structure.(example data structure is attached). Any insight would be appreciated.

nm_sim.csv (913 Bytes)

I was able to run your data with one of my .stan files without any problems. Can you attach an R file and your .stan file with the data (the data you already attached is fine if it recreates the problem) that produces the error?

Hi cbdavis33, Thank you, This stan code runs. I am not using stanpmx-lib functions this time. But it takes long time, 3 days, still running for data around (8000) obs. Using parallel_chains = 4. Any suggestion for reducing run time or by using stanpmx-lib?. Above Stan and R code is attached.
nm_sim.csv (1.4 KB)

testmod.R (15.7 KB)

Ok, that’s a bit too messy for me to go through (try attaching, not copy-and-paste. The comments in an R file create sections in Markdown), but a couple preliminary thoughts:

  1. In the data, no need for an observation record at t = 0. No drug yet = no concentration. Removing that will help something, not sure what, in your first code, though (the one with Torsten),
  2. In the other one where you’re writing the pre-Torsten version of code, it doesn’t have a problem with that observation record at t = 0. I wouldn’t really bother spendiing time speeding it up. I’d recommend getting the isize issue sorted out using Torsten and then parallelizing that.

Can you attach (or at least put in a code block) an R file and your .stan file that produces the isize error? I need something that reproduces the isize error to debug it.

Thanks, Attached above. isize was resolved by changing data structure.