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.
+ }
+ }
+ }
+ }