v0.4.0

.Rbuildignore

@@ -5,7 +5,6 @@
 ^\_config\.yml$
 ^\.github$
 ^CONTRIBUTING\.md$
-^NEWS\.md$
 ^_pkgdown\.yml$
 pkgdown/
 .DS_Store

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: bmstate
 Type: Package
 Title: Bayesian multistate modeling
-Version: 0.3.3
+Version: 0.4.0
 Authors@R:
     c(person(given = "Juho",
            family = "Timonen",

NEWS.md

@@ -0,0 +1,85 @@
+# bmstate 0.4.0
+
+* **Correct the baseline-hazard prior scale.** The public model API stores
+  baseline-hazard reference values as natural-scale rates, but these values
+  were passed directly to a Stan parameter used as the location of
+  `log_w0`. For example, a reference rate of 0.001 therefore centered the
+  hazard near `exp(0.001) = 1.001`, rather than 0.001. Stan data creation now
+  applies `log()` exactly once at the R--Stan boundary. The default Stan data
+  field is consequently renamed from `mu_w0` to `mu_log_w0`; custom Stan
+  programs selected with `options(bmstate_stan_file = ...)` must adopt the new
+  name.
+  Fits produced by affected versions must be rerun.
+
+* **Make transition support and baseline-hazard calibration explicit.** The
+  previous cumulative-hazard summary silently omitted modeled transitions with
+  no events, which could leave empirical prior locations undefined and defer
+  failure until Stan data were validated. Fit-time validation now stops with
+  the affected transition names if any transition in the fitted matrix has no
+  observed occurrence or no positive at-risk time. This check is deliberately
+  not applied to out-of-sample prediction data. For valid transition systems,
+  reference rates are computed from event counts divided by exact at-risk
+  time, pooled within transition type. Documentation now calls the calibrated
+  quantity a natural-scale reference rate, not an arithmetic prior mean.
+
+* **Preserve joint posterior-draw alignment in prediction.** Replicating the
+  baseline-hazard arrays previously produced a different draw/subject ordering
+  from the hazard multipliers and prediction metadata, so parameters from one
+  posterior draw could be combined with another. The package now constructs an
+  explicit draw-major row index and uses the same ordering for spline weights,
+  log intercepts, subject multipliers, path links, and occupancy calculations.
+
+* **Retain zero state-visit probabilities.** `p_state_visit()` previously built
+  its output only from observed event rows, silently dropping requested
+  subjects or groups with no visit to an event state. It now completes the
+  state-by-group grid, joins event counts and denominators separately, and
+  returns explicit zero counts and probabilities.
+
+* **Keep simulated paths inside the requested horizon.** Applying
+  `min_t_step` could move a genuine event beyond `t_max`; the final-row cleanup
+  then relabeled that row as censoring without clearing its transition index.
+  Event-time adjustment now preserves a genuine pre-horizon event when the
+  adjustment would cross the horizon, every censor row has transition index
+  zero, invalid reversed horizons are rejected, and a zero-length horizon
+  returns a valid one-row censored path.
+
+* **Stabilize the oral one-compartment PK solution near equal rates.** Direct
+  evaluation divided by `ka - ke` and could produce cancellation, `0 / 0`, or
+  `0 * Inf` when absorption and elimination rates were equal or very different.
+  The transient solution now uses sign-aware `expm1()` factorizations, and the
+  transient and steady-state solutions use their continuous equal-rate limits
+  in a small numerical guard region. The limit uses the `ka` prefactor so that
+  both the value and first derivative match at equality to first order.
+
+* **Construct exactly the requested integration grid.** For a noninteger
+  maximum time, the old `seq(..., ceiling(t_max), by = delta)` expression could
+  generate more than `N_grid` points and make the declared Stan dimensions
+  disagree with the supplied arrays. The grid is now formed directly from
+  exactly `N_grid` midpoint indices.
+
+* **Honor delayed starts in transition-probability prediction.**
+  `solve_trans_prob_matrix()` previously initialized the forward equation at
+  the first output time, even when an earlier `t_start` was requested. It now
+  includes `t_start` in the integration times, returns only the requested
+  output rows, and validates times and hazard-parameter dimensions even on the
+  zero-duration identity path.
+
+* **Return the requested PK credible-interval width.** The plotting helper
+  divided the tail probability by two twice, so an 80% request used the 5th and
+  95th percentiles and returned a 90% interval. It now uses
+  `(1 - ci_alpha) / 2` and its complement directly.
+
+* **Reduce reverse-mode work for time-invariant hazard predictors.** The Stan
+  program previously expanded subject covariates to every interval and, before
+  hoisting, rebuilt the complete multiplier matrix inside each transition
+  loop. It now differentiates one predictor per subject and transition type,
+  then gathers those values for the relevant intervals. This is algebraically
+  identical and reduced gradient time by about 38% in a representative
+  benchmark with many intervals, with no measured regression when every
+  interval belonged to a different subject.
+
+* **Tighten edge-case validation and output construction.** The release adds
+  checks for positive PK rates, dosing intervals, concentration bounds, and
+  exposure-normalization scales; rejects malformed delayed-start requests;
+  handles base data frames without accidental dimension dropping; and builds
+  simulated path output in a list before binding it once.

R/DosingData.R

@@ -67,6 +67,9 @@ PSSDosingData <- R6::R6Class(
       self$set_sub_ids(subject_ids)
       N_sub <- self$num_subjects()
       checkmate::assert_number(tau_ss, lower = 0)
+      if (tau_ss <= 0) {
+        stop("tau_ss must be strictly positive")
+      }
       checkmate::assert_list(doses, len = N_sub)
       checkmate::assert_list(times, len = N_sub)
       self$doses <- doses
@@ -269,6 +272,9 @@ simulate_dosing <- function(df_subjects, tau = 24, p_miss = 0.2, t_jitter = 4) {
 pk_2cpt_pss <- function(t, dose_ss, times, doses, theta, tau, MAX_CONC) {
   ensure_exposed_stan_functions()
   checkmate::assert_number(tau, lower = 0)
+  if (tau <= 0) {
+    stop("tau must be strictly positive")
+  }
   checkmate::assert_number(MAX_CONC, lower = 0)
   checkmate::assert_numeric(dose_ss, lower = 0)
   N_sub <- length(dose_ss)

R/MultistateModel.R

@@ -110,7 +110,7 @@ MultistateModel <- R6::R6Class("MultistateModel",
 
       log_w0 <- matrix(rep(log_w0, N), N, S, byrow = TRUE)
       log_m <- private$simulate_log_hazard_multipliers(df_subjects, beta_haz)
-      paths <- self$system$simulate(w_all, log_w0, log_m, min_t_step = 0.1)
+      paths <- self$system$simulate(w_all, log_w0, log_m, min_t_step = 1e-6)
       as_tibble(paths)
     }
   ),
@@ -163,8 +163,11 @@ MultistateModel <- R6::R6Class("MultistateModel",
     #' @param loc Location
     #' @param scale Scale
     set_xpsr_normalizers = function(loc = 0, scale = 1) {
-      checkmate::assert_numeric(loc, lower = 0, len = 1)
-      checkmate::assert_numeric(scale, lower = 0, len = 1)
+      checkmate::assert_number(loc, finite = TRUE)
+      checkmate::assert_number(scale, finite = TRUE)
+      if (scale <= 0) {
+        stop("xpsr normalization scale must be strictly positive")
+      }
       message(
         "setting xpsr normalizers to loc = ",
         round(loc, 5), ", scale = ", round(scale, 5)
@@ -174,35 +177,43 @@ MultistateModel <- R6::R6Class("MultistateModel",
       invisible(NULL)
     },
 
-    #' @description Get assumed prior mean baseline hazard rates.
-    #' @return Numeric vector with length equal to number of transitions
+    #' @description Get baseline-hazard prior reference rates.
+    #' @return A positive numeric vector with length equal to the number of
+    #' transitions. These natural-scale rates are transformed to locations on
+    #' the log-hazard scale when Stan data are created; they are not arithmetic
+    #' prior means.
     get_prior_mean_h0 = function() {
       v <- private$prior_mean_h0
       if (is.null(v)) {
-        stop("prior mean h0 has not been set")
+        stop("prior reference h0 has not been set")
       }
       v
     },
 
-    #' @description Set assumed prior mean baseline hazard rates (side
-    #' effect).
-    #' @param mean_h0 Numeric vector with length equal to number of transitions
+    #' @description Set baseline-hazard prior reference rates (side effect).
+    #' @param mean_h0 A positive numeric vector of natural-scale reference rates
+    #' with length equal to the number of transitions. The argument name is
+    #' retained for compatibility; the values are not arithmetic prior means.
     set_prior_mean_h0 = function(mean_h0) {
       N_trans <- self$system$tm()$num_trans()
-      checkmate::assert_numeric(mean_h0, len = N_trans, lower = 0)
+      checkmate::assert_numeric(mean_h0, len = N_trans, finite = TRUE)
+      if (any(mean_h0 <= 0)) {
+        stop("prior baseline hazard rates must be strictly positive")
+      }
       private$prior_mean_h0 <- mean_h0
       invisible(NULL)
     },
 
-    #' @description Set assumed prior mean baseline hazard rates (side
-    #' effect) based on average hazards in data.
+    #' @description Set baseline-hazard prior reference rates (side effect)
+    #' using event counts and at-risk time pooled within transition type.
     #' @param data A \code{\link{JointData}} or \code{\link{PathData}} object.
     set_prior_mean_h0_data = function(data) {
       if (inherits(data, "PathData")) {
         data <- JointData$new(data, NULL)
       }
       checkmate::assert_class(data, "JointData")
-      df_ttype <- average_haz_per_ttype(data$paths) |> dplyr::arrange(.data$trans_idx)
+      df_ttype <- average_haz_per_ttype(data$paths) |>
+        dplyr::arrange(.data$trans_idx)
       self$set_prior_mean_h0(exp(df_ttype$log_h0_avg))
     },
 

R/MultistateModelFit.R

@@ -78,6 +78,9 @@ MultistateModelFit <- R6::R6Class("MultistateModelFit",
     #'
     #' @param name Param/quantity name
     get_draws = function(name = NULL) {
+      if (is.null(name)) {
+        return(private$draws)
+      }
       private$draws[[name]]
     },
 
@@ -261,7 +264,7 @@ MultistateModelFit <- R6::R6Class("MultistateModelFit",
     #' @description
     #' Full names of parameters that start with \code{log_z_}.
     log_z_pars = function() {
-      nams <- names(self$draws())
+      nams <- self$draws_names()
       match <- grepl(nams, pattern = "log_z_")
       nams[which(match)]
     },
@@ -508,40 +511,32 @@ msmfit_log_hazard_multipliers <- function(fit, oos = FALSE, data = NULL) {
     beta_xpsr <- array(0, dim = c(S, 1, 0, sd$N_trans_types))
   }
 
-  # Create x_haz_long (long version of hazard covariates vector)
-  N_sub <- sd$N_sub
-  first_indices <- sapply(seq_len(N_sub), function(x) which(sd$idx_sub == x)[1])
-  if (sd$nc_haz > 0) {
-    x_haz_long <- sd$x_haz[, sd$idx_sub, drop = FALSE]
-  } else {
-    x_haz_long <- array(0, dim = c(0, sd$N_int))
-  }
   an <- fit$model$get_xpsr_normalizers()
 
   # Call exposed Stan function for each draw (not optimal)
   out <- list()
   for (s in seq_len(S)) {
     if (sd$do_pk == 1) {
       ba <- list(beta_xpsr[s, 1, 1, ])
-      aa <- list(xpsr[[s]][sd$idx_sub])
+      aa <- list(xpsr[[s]])
     } else {
       ba <- NULL
       aa <- NULL
     }
     if (sd$nc_haz == 0 && sd$do_pk == 0) {
       r <- matrix(0, sd$N_sub, sd$N_trans)
     } else {
-      r <- compute_log_hazard_multiplier(
-        sd$N_int,
+      r_by_type <- compute_log_hazard_multiplier(
+        sd$N_sub,
+        sd$N_trans_types,
         mat2list(t(array(beta_oth[s, 1, , ], dim = c(sd$nc_haz, sd$N_trans_types)))),
         ba,
-        mat2list(t(x_haz_long)),
+        mat2list(t(sd$x_haz)),
         aa,
         an$loc,
-        an$scale,
-        sd$ttype
+        an$scale
       )
-      r <- r[first_indices, , drop = FALSE]
+      r <- r_by_type[, sd$ttype, drop = FALSE]
     }
     out[[s]] <- r
   }
@@ -581,6 +576,21 @@ msmfit_log_baseline_hazard <- function(fit, t = NULL) {
 }
 
 
+# Repeat each joint posterior draw for every subject while preserving draw-major
+# order. Keeping this in one helper makes the required parameter alignment
+# explicit and independently testable.
+repeat_hazard_parameter_draws <- function(w, log_w0, n_subjects) {
+  checkmate::assert_array(w, d = 3)
+  checkmate::assert_matrix(log_w0, nrows = dim(w)[1])
+  checkmate::assert_integerish(n_subjects, len = 1, lower = 1)
+  draw_rows <- rep(seq_len(dim(w)[1]), each = n_subjects)
+  list(
+    w = w[draw_rows, , , drop = FALSE],
+    log_w0 = log_w0[draw_rows, , drop = FALSE]
+  )
+}
+
+
 #' Extract and reshape draws of instant hazard related parameters
 #'
 #' @export
@@ -600,8 +610,9 @@ msmfit_inst_hazard_param_draws <- function(fit, oos = FALSE, data = NULL) {
   log_w0 <- fit$get_draws_of("log_w0")
   w <- access_one_dim(w, dim_index = 2, value = 1)
   log_w0 <- access_one_dim(log_w0, dim_index = 2, value = 1)
-  w_rep <- abind::abind(replicate(N, w, simplify = FALSE), along = 1)
-  log_w0_rep <- abind::abind(replicate(N, log_w0, simplify = FALSE), along = 1)
+  repeated <- repeat_hazard_parameter_draws(w, log_w0, N)
+  w_rep <- repeated$w
+  log_w0_rep <- repeated$log_w0
   log_m_reshaped <- do.call(rbind, log_m)
 
   # Subject-draw df
@@ -646,8 +657,6 @@ generate_paths <- function(fit, oos = FALSE, t_start = 0, t_max = NULL, n_rep =
   checkmate::assert_class(fit, "MultistateModelFit")
   checkmate::assert_integerish(n_rep, lower = 1, len = 1)
   sd <- msmfit_stan_data(fit, data)
-  log_m <- msmfit_log_hazard_multipliers(fit, oos, data)
-
   # Get and reshape draws and state vector at t_start
   sys <- fit$model$system
   S <- fit$num_draws()
@@ -729,9 +738,13 @@ p_state_occupancy <- function(fit, oos = FALSE, t_start = 0, t_out = NULL,
 
   sys <- fit$model$system
   if (is.null(t_out)) {
-    t_out <- seq(t_start, sys$get_tmax(), length.out = 30)
+    if (t_start == sys$get_tmax()) {
+      t_out <- t_start
+    } else {
+      t_out <- seq(t_start, sys$get_tmax(), length.out = 30)
+    }
   } else {
-    checkmate::assert_numeric(t_out, min.len = 2)
+    checkmate::assert_numeric(t_out, min.len = 1, finite = TRUE)
   }
 
   # Get and reshape draws