From c960f39c2bcaec8782e7c4255995bc602637fbce Mon Sep 17 00:00:00 2001 From: jtimonen Date: Mon, 8 Dec 2025 10:52:29 +0200 Subject: [PATCH 1/7] Increment version number to 0.3.3.9000 --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index e837edc..6fc326c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: bmstate Type: Package Title: Bayesian multistate modeling -Version: 0.3.3 +Version: 0.3.3.9000 Authors@R: c(person(given = "Juho", family = "Timonen", From 6a46ede7117c57098639e4c0905bb4ad1b43c48f Mon Sep 17 00:00:00 2001 From: jtimonen Date: Tue, 28 Jul 2026 11:11:25 +0300 Subject: [PATCH 2/7] bugfixes --- DESCRIPTION | 2 +- NEWS.md | 14 + R/DosingData.R | 6 + R/MultistateModel.R | 17 +- R/MultistateModelFit.R | 27 +- R/MultistateSystem.R | 100 +++++-- R/PKModel.R | 22 +- R/PathData.R | 62 +++-- R/stan-data.R | 82 ++++-- R/stan.R | 15 +- R/utils.R | 14 +- inst/stan/msm.stan | 97 ++++--- man/MultistateSystem.Rd | 3 +- tests/testthat/test-correctness-regressions.R | 257 ++++++++++++++++++ tests/testthat/test-stan.R | 13 + 15 files changed, 608 insertions(+), 123 deletions(-) create mode 100644 NEWS.md create mode 100644 tests/testthat/test-correctness-regressions.R diff --git a/DESCRIPTION b/DESCRIPTION index 6fc326c..ffffec2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: bmstate Type: Package Title: Bayesian multistate modeling -Version: 0.3.3.9000 +Version: 0.3.3.9001 Authors@R: c(person(given = "Juho", family = "Timonen", diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..afdd900 --- /dev/null +++ b/NEWS.md @@ -0,0 +1,14 @@ +# bmstate 0.3.3.9001 + +* Correct the empirical baseline-hazard prior scale and use robust event-time + exposure estimates for sparse transitions. +* Preserve joint posterior-draw alignment in path and occupancy prediction. +* Return explicit zero state-visit probabilities for all requested groups. +* Use a guaranteed B-spline thinning bound and remove rare-event truncation. +* Stabilize the oral one-compartment PK solution when absorption and elimination + rates are equal or nearly equal. +* Fix noninteger integration grids, delayed transition-probability start times, + requested PK credible-interval widths, and zero-event Stan data. +* Center spline weights to remove their shift redundancy with the baseline + log-hazard intercept, and avoid recomputing invariant hazard multipliers in + the Stan transition loop. diff --git a/R/DosingData.R b/R/DosingData.R index 5a76b73..789fb57 100644 --- a/R/DosingData.R +++ b/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) diff --git a/R/MultistateModel.R b/R/MultistateModel.R index 020b4e1..eb21545 100644 --- a/R/MultistateModel.R +++ b/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) @@ -189,7 +192,10 @@ MultistateModel <- R6::R6Class("MultistateModel", #' @param mean_h0 Numeric vector with length equal to number of transitions 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) }, @@ -202,7 +208,8 @@ MultistateModel <- R6::R6Class("MultistateModel", 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)) }, diff --git a/R/MultistateModelFit.R b/R/MultistateModelFit.R index c16d01c..3f8f4f7 100644 --- a/R/MultistateModelFit.R +++ b/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)] }, @@ -581,6 +584,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 +618,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 +665,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() diff --git a/R/MultistateSystem.R b/R/MultistateSystem.R index 3078bf0..8f8dfbc 100644 --- a/R/MultistateSystem.R +++ b/R/MultistateSystem.R @@ -68,19 +68,23 @@ MultistateSystem <- R6::R6Class("MultistateSystem", # Generate transition given state and transition intensity functions generate_transition = function(state, t_init, t_max, w, log_w0, log_m) { - tol <- 1.05 possible <- private$transmat$possible_transitions_from(state) - J <- private$transmat$num_trans() - UB <- rep(0, J) - for (j in seq_len(J)) { - UB[j] <- tol * self$max_inst_hazard(t_init, t_max, w[j, ], log_w0[j], log_m[j]) - } - possible <- possible[which(UB[possible] > 1e-9)] # so rare that not possible if (length(possible) == 0) { # Absorbing state, no transitions possible return(list(t = t_max, new_state = 0, idx = 0)) } - UB <- UB[possible] + UB <- vapply(possible, function(j) { + self$max_inst_hazard(t_init, t_max, w[j, ], log_w0[j], log_m[j]) + }, numeric(1)) + if (any(!is.finite(UB))) { + stop("non-finite transition-hazard upper bound") + } + keep <- UB > 0 + possible <- possible[keep] + UB <- UB[keep] + if (length(possible) == 0) { + return(list(t = t_max, new_state = 0, idx = 0)) + } w <- w[possible, , drop = FALSE] log_w0 <- log_w0[possible] log_m <- log_m[possible] @@ -276,7 +280,11 @@ MultistateSystem <- R6::R6Class("MultistateSystem", knots <- self$get_knots() L <- length(knots) BK <- knots[c(1, L)] - knots <- knots[2:(L - 1)] + if (L > 2) { + knots <- knots[2:(L - 1)] + } else { + knots <- numeric(0) + } bspline_basis(t, private$spline_k, knots, BK) }, @@ -284,7 +292,8 @@ MultistateSystem <- R6::R6Class("MultistateSystem", #' #' @param t Time point(s). Not used if \code{SBF} is given. #' @param w Spline basis function weights (vector) - #' @param log_w0 Intercept (log) + #' @param log_w0 Intercept (log), or \code{-Inf} for an identically zero + #' baseline hazard #' @param log_m Hazard multiplier (log) #' @param SBF Pre-computed basis function matrix at \code{t}. log_inst_hazard = function(t, w, log_w0, log_m, SBF = NULL) { @@ -336,12 +345,51 @@ MultistateSystem <- R6::R6Class("MultistateSystem", #' @param t1 Start time point #' @param t2 End time point #' @param w Spline basis function weights (vector) - #' @param log_w0 Intercept (log) + #' @param log_w0 Intercept (log), or \code{-Inf} for an identically zero + #' baseline hazard #' @param log_m Hazard multiplier (log) max_inst_hazard = function(t1, t2, w, log_w0, log_m) { - ttt <- seq(t1, t2, length.out = 100) - log_haz <- log_m + self$log_baseline_hazard(ttt, log_w0, w) - return(exp(max(log_haz))) + checkmate::assert_number(t1, lower = min(self$get_knots())) + checkmate::assert_number(t2, lower = t1, upper = self$get_tmax()) + checkmate::assert_numeric(w, len = self$num_weights(), finite = TRUE) + checkmate::assert_number(log_w0, upper = .Machine$double.xmax) + checkmate::assert_number(log_m, finite = TRUE) + if (is.infinite(log_w0)) { + return(0) + } + + # An intercept-inclusive B-spline basis is nonnegative and partitions + # unity inside its boundary knots. On a given interval, its weighted sum + # is therefore bounded above by the largest coefficient whose compact + # support intersects that interval. + knots <- self$get_knots() + order <- private$spline_k + internal_knots <- if (length(knots) > 2) { + knots[2:(length(knots) - 1)] + } else { + numeric(0) + } + augmented_knots <- c( + rep(knots[1], order), + internal_knots, + rep(knots[length(knots)], order) + ) + weight_idx <- seq_len(self$num_weights()) + support_left <- augmented_knots[weight_idx] + support_right <- augmented_knots[weight_idx + order] + active <- support_right > t1 & support_left < t2 + + # Include basis functions active exactly at either closed endpoint. This + # also handles a zero-length interval without broadening to global support. + endpoint_basis <- self$basisfun_matrix(unique(c(t1, t2))) + active <- active | colSums(endpoint_basis > 0) > 0 + if (!any(active)) { + stop("could not identify an active B-spline basis function") + } + # A negligible upward margin protects thinning from floating-point + # overshoot when the same hazard is subsequently evaluated numerically. + log_margin <- log1p(sqrt(.Machine$double.eps)) + exp(log_w0 + log_m + max(w[active]) + log_margin) }, #' @description Generate paths @@ -379,7 +427,7 @@ MultistateSystem <- R6::R6Class("MultistateSystem", pb <- progress::progress_bar$new(total = n_paths) # Set max time - out <- NULL + out <- vector("list", n_paths) if (is.null(t_max)) { t_max <- self$get_tmax() } @@ -397,10 +445,10 @@ MultistateSystem <- R6::R6Class("MultistateSystem", min_t_step ) p <- cbind(p, rep(cnt, nrow(p))) - out <- rbind(out, p) + out[[cnt]] <- p } } - df <- data.frame(out) + df <- data.frame(do.call(rbind, out)) colnames(df)[ncol(df)] <- "path_id" df } @@ -469,9 +517,15 @@ solve_trans_prob_matrix <- function(system, t_out, log_w0, w = NULL, log_m = NULL, t_start = 0, ...) { checkmate::assert_class(system, "MultistateSystem") - checkmate::assert_numeric(t_out) + checkmate::assert_numeric(t_out, min.len = 1, finite = TRUE) K <- length(t_out) checkmate::assert_number(t_start, lower = 0) + if (any(t_out < t_start)) { + stop("all t_out values must be greater than or equal to t_start") + } + if (is.unsorted(t_out, strictly = TRUE)) { + stop("t_out must be strictly increasing") + } H <- system$num_trans() W <- system$num_weights() S <- system$num_states() @@ -481,10 +535,16 @@ solve_trans_prob_matrix <- function(system, t_out, log_w0, w = NULL, if (is.null(log_m)) { log_m <- rep(0, H) } - kfe <- solve_time_evolution(system, t_out, log_w0, w, log_m, ...) + solve_times <- unique(c(t_start, t_out)) + if (length(solve_times) == 1) { + P <- array(diag(1, S, S), dim = c(1, S, S)) + return(P) + } + kfe <- solve_time_evolution(system, solve_times, log_w0, w, log_m, ...) + output_rows <- match(t_out, kfe[, 1]) P <- array(0, dim = c(K, S, S)) for (k in seq_len(K)) { - P[k, , ] <- matrix(kfe[k, 2:ncol(kfe)], S, S) + P[k, , ] <- matrix(kfe[output_rows[k], 2:ncol(kfe)], S, S) } P } diff --git a/R/PKModel.R b/R/PKModel.R index dac6723..5369fda 100644 --- a/R/PKModel.R +++ b/R/PKModel.R @@ -23,6 +23,9 @@ PKModel <- R6::R6Class("PKModel", #' @param value Upper bound for concentration, to avoid numerical issues. set_max_conc = function(value) { checkmate::assert_number(value, lower = 0) + if (value <= 0) { + stop("concentration upper bound must be strictly positive") + } message("setting max conc = ", round(value, 5)) private$MAX_CONC <- value invisible(NULL) @@ -90,16 +93,23 @@ PKModel <- R6::R6Class("PKModel", checkmate::assert_number(theta$V2, lower = 0) checkmate::assert_number(dose, lower = 0) checkmate::assert_number(tau, lower = 0) + if (theta$ka <= 0 || theta$CL <= 0 || theta$V2 <= 0 || tau <= 0) { + stop("ka, CL, V2, and tau must be strictly positive") + } N <- length(t) ka <- theta$ka ke <- theta$CL / theta$V2 - A <- (dose / theta$V2) * (ka / (ka - ke)) - conc <- rep(0, N) - ma <- A / (-expm1(-ka * tau)) - me <- A / (-expm1(-ke * tau)) tt <- t %% tau - conc <- me * exp(-ke * tt) - ma * exp(-ka * tt) - conc + d <- ka - ke + rel_tol <- 1e-8 * max(abs(ka), abs(ke)) + if (abs(d) <= rel_tol) { + k <- 0.5 * (ka + ke) + f <- exp(-k * tt) / (-expm1(-k * tau)) + return((dose / theta$V2) * k * f * (tt + tau / expm1(k * tau))) + } + f_ke <- exp(-ke * tt) / (-expm1(-ke * tau)) + f_ka <- exp(-ka * tt) / (-expm1(-ka * tau)) + (dose / theta$V2) * ka / d * (f_ke - f_ka) }, #' @description Compute exposure, which is the steady-state diff --git a/R/PathData.R b/R/PathData.R index 36d52e5..68bf9db 100644 --- a/R/PathData.R +++ b/R/PathData.R @@ -85,7 +85,7 @@ PathData <- R6::R6Class( checkmate::assert_integerish(path_df$path_id) # As tibbles - subject_df <- as_tibble(subject_df[, cols1]) + subject_df <- as_tibble(subject_df[, cols1, drop = FALSE]) path_df <- as_tibble(path_df[, cols2]) link_df <- as_tibble(link_df[, cols3]) @@ -444,8 +444,8 @@ msfit_average_hazard <- function(msfit) { msfit$Haz |> dplyr::group_by(.data$trans) |> summarise( - avg_haz = (dplyr::last(.data$Haz) - dplyr::first(.data$Haz)) / - (dplyr::last(.data$time) - dplyr::first(.data$time)) + avg_haz = dplyr::last(.data$Haz) / dplyr::last(.data$time), + .groups = "drop" ) } @@ -650,29 +650,50 @@ p_state_visit <- function(pd, t = NULL, by = NULL) { } checkmate::assert_numeric(t, lower = 0) - if (!is.null(by)) { - checkmate::assert_character(by, len = 1) - c <- pd$as_data_frame(covariates = by) |> - dplyr::group_by(.data$state, .data[[by]]) - } else { - c <- pd$as_data_frame() |> - dplyr::group_by(.data$state) - } estates <- which(pd$state_names() %in% pd$get_event_state_names()) - df <- count_paths_with_event(c, t, S) |> dplyr::filter(.data$state %in% estates) + events <- pd$as_data_frame(covariates = by) |> + dplyr::filter( + .data$time <= t, + .data$state %in% estates + ) + if (!is.null(by)) { - df_all <- c |> dplyr::ungroup() - df_all <- df_all |> dplyr::group_by(.data[[by]]) - df_all <- df_all |> - dplyr::distinct(.data$path_id) |> - dplyr::count() - df <- df |> dplyr::left_join(df_all, by = by) + checkmate::assert_character(by, len = 1) + links <- pd$full_link(covariates = by) |> + dplyr::distinct(.data$path_id, .data[[by]]) + groups <- links |> + dplyr::distinct(.data[[by]]) + grid <- tidyr::expand_grid( + state_idx = estates, + group_value = groups[[by]] + ) + names(grid)[2] <- by + counts <- events |> + dplyr::transmute( + path_id = .data$path_id, + state_idx = .data$state, + group_value = .data[[by]] + ) |> + dplyr::distinct() |> + dplyr::count(.data$state_idx, .data$group_value, name = "n_event") + names(counts)[2] <- by + totals <- links |> + dplyr::count(.data[[by]], name = "n") + df <- grid |> + dplyr::left_join(counts, by = c("state_idx", by)) |> + dplyr::left_join(totals, by = by) + df$n_event[is.na(df$n_event)] <- 0L df$prob <- df$n_event / df$n } else { + counts <- events |> + dplyr::transmute(path_id = .data$path_id, state_idx = .data$state) |> + dplyr::distinct() |> + dplyr::count(.data$state_idx, name = "n_event") + df <- data.frame(state_idx = estates) |> + dplyr::left_join(counts, by = "state_idx") + df$n_event[is.na(df$n_event)] <- 0L df$prob <- df$n_event / pd$n_paths() } - df$state_idx <- df$state - df$state <- NULL sdf <- pd$transmat$states_df() df |> dplyr::left_join(sdf, by = "state_idx") } @@ -783,6 +804,7 @@ df_to_paths_df_part2 <- function(pdf, tm) { #' @return A \code{\link{PathData}} object df_to_pathdata <- function(df, tm, covs = NULL, validate = TRUE) { check_columns(df, c("state", "time", "subject_id", "is_transition")) + df <- tibble::as_tibble(df) checkmate::assert_integerish(df$state) checkmate::assert_numeric(df$time) checkmate::assert_character(df$subject_id) diff --git a/R/stan-data.R b/R/stan-data.R index b67bd64..7794fc0 100644 --- a/R/stan-data.R +++ b/R/stan-data.R @@ -31,7 +31,9 @@ create_stan_data_model <- function(model) { N_trans = tm$num_trans(), N_trans_types = tm$num_trans_types(), ttype = tm$trans_df()$trans_type, - mu_w0 = model$get_prior_mean_h0(), + # The public API specifies prior locations as hazard rates. Stan models + # log hazards, so transform exactly once at this boundary. + mu_log_w0 = log(model$get_prior_mean_h0()), nc_haz = length(model$data_covs("haz")), N_sbf = model$system$num_weights() ) @@ -85,7 +87,7 @@ create_stan_data_timegrid <- function(model) { if (delta_grid > 0.25 * t_max) { stop("delta_grid is very large compared to t_max") } - t_grid <- seq(delta_grid / 2, ceiling(t_max), by = delta_grid) # midpoints + t_grid <- (seq_len(G) - 0.5) * delta_grid # exactly G midpoints # Return list( @@ -227,11 +229,17 @@ standata_scaled_covariates <- function(pd, model, name) { for (cn in covs) { j <- j + 1 xx <- sub_df[[cn]] + if (!is.numeric(xx)) { + stop("covariate '", cn, "' must be numeric before normalization") + } xj_loc <- norms$locations[[cn]] xj_scale <- norms$scales[[cn]] if (is.null(xj_loc) || is.null(xj_scale)) { stop("no normalizers set for ", cn) } + if (!is.finite(xj_loc) || !is.finite(xj_scale) || xj_scale <= 0) { + stop("normalization scale for '", cn, "' must be finite and positive") + } x_norm <- (xx - xj_loc) / xj_scale check_normalized_covariate(x_norm, cn) x[[j]] <- x_norm @@ -295,29 +303,69 @@ create_stan_data_intervalidx <- function(t_start, t_end, t_grid, delta_grid) { ) } -# Average hazard per transition type, ignoring transitions that did not -# occur +# Estimate a constant hazard per transition type by pooling event counts and +# at-risk time over all transitions of that type. If an exposed type has no +# events, one Jeffreys-style half-event avoids an undefined empirical prior. +# Entirely unexposed types borrow the global pooled transition rate. average_haz_per_ttype <- function(pd) { - msfit <- pd$fit_mstate() - h0 <- msfit_average_hazard(msfit) |> dplyr::arrange(.data$trans) - df_trans <- pd$transmat$trans_df() - df_ttype <- df_trans |> dplyr::select("trans_idx", "trans_type") - df_ttype$trans <- df_ttype$trans_idx - h0 <- h0 |> left_join(df_ttype, by = "trans") - df_mean_log_h0 <- h0 |> - dplyr::filter(.data$avg_haz > 0) |> + checkmate::assert_class(pd, "PathData") + dat <- pd$as_transitions(truncate = TRUE) + trans <- pd$transmat$trans_df() + exposure <- numeric(nrow(trans)) + n_event <- integer(nrow(trans)) + + for (h in seq_len(nrow(trans))) { + at_risk <- dat$from == trans$prev_state[h] + exposure[h] <- sum(dat$time[at_risk] - dat$time_prev[at_risk]) + n_event[h] <- sum(dat$trans_idx == trans$trans_idx[h]) + if (exposure[h] <= 0 && n_event[h] > 0) { + stop("transition ", h, " has events but no positive at-risk time") + } + } + + total_exposure <- sum(exposure) + total_events <- sum(n_event) + if (total_exposure <= 0) { + stop("cannot calibrate baseline-hazard priors: no transition has positive at-risk time") + } + global_rate <- if (total_events > 0) { + total_events / total_exposure + } else { + 0.5 / total_exposure + } + + rate_df <- data.frame( + trans_idx = trans$trans_idx, + trans_type = trans$trans_type, + exposure = exposure, + n_event = n_event + ) + type_df <- rate_df |> dplyr::group_by(.data$trans_type) |> - mutate(log_haz = log(.data$avg_haz)) |> - summarize(log_h0_avg = mean(.data$log_haz)) |> - dplyr::arrange(.data$trans_type) - df_ttype |> left_join(df_mean_log_h0, by = "trans_type") + dplyr::summarise( + exposure = sum(.data$exposure), + n_event = sum(.data$n_event), + .groups = "drop" + ) |> + dplyr::mutate( + rate = ifelse( + .data$exposure > 0, + ifelse(.data$n_event > 0, .data$n_event, 0.5) / .data$exposure, + global_rate + ), + log_h0_avg = log(.data$rate) + ) |> + dplyr::select("trans_type", "log_h0_avg") + rate_df |> + dplyr::select("trans_idx", "trans_type") |> + dplyr::left_join(type_df, by = "trans_type") } # Edit Stan data, sort of computes which() for the binary vectors of each # transition which_format_for_stan <- function(x, name) { H <- nrow(x) - M <- max(rowSums(x)) + M <- max(1L, max(rowSums(x))) N_sum <- rep(0, H) a <- matrix(0, H, M) for (h in seq_len(H)) { diff --git a/R/stan.R b/R/stan.R index 0403066..526ff3e 100644 --- a/R/stan.R +++ b/R/stan.R @@ -100,11 +100,24 @@ fit_stan <- function(model, data, # Set covariate normalizing locations and scales (side effect) if (set_normalizers) { model$set_normalizers(data) - if (!is.null(data$dosing)) { + if (model$has_pk()) { + if (is.null(data$dosing)) { + stop("a model with a PK component requires dosing data") + } log_mu_CL <- -2 # should match msm.stan + if (any(!is.finite(data$dosing$dose_ss)) || any(data$dosing$dose_ss <= 0)) { + stop("steady-state doses must be finite and positive for log-exposure modeling") + } aaa <- log(data$dosing$dose_ss) - log_mu_CL loc <- mean(aaa) sca <- stats::sd(aaa) + if (!is.finite(sca) || sca <= sqrt(.Machine$double.eps)) { + sca <- 0.5 + warning( + "steady-state dose has no usable between-subject variation; ", + "using the default xpsr normalization scale of 0.5" + ) + } model$set_xpsr_normalizers(loc, sca) } } diff --git a/R/utils.R b/R/utils.R index 3c97474..469186e 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1,11 +1,11 @@ # Helper pksim_to_quantiles <- function(sim, ci_alpha) { checkmate::assert_number(ci_alpha, lower = 0, upper = 1) - av <- (1 - ci_alpha) / 2 + tail_prob <- (1 - ci_alpha) / 2 sim <- sim |> dplyr::group_by(.data$subject_id, .data$time) |> dplyr::summarise(q = list( - stats::quantile(.data$val, probs = c(av / 2, 0.5, 1 - av / 2)) + stats::quantile(.data$val, probs = c(tail_prob, 0.5, 1 - tail_prob)) ), .groups = "drop") |> tidyr::unnest_wider(q, names_sep = "_") colnames(sim)[3:5] <- c("lower", "val", "upper") @@ -117,10 +117,14 @@ sim_subject_ids <- function(N) { # Truncate path data frame truncate_after_terminal_events <- function(df, term_state_inds) { - term_events <- df |> - dplyr::filter(.data$state %in% term_state_inds, .data$trans_idx > 0) |> + term_rows <- df |> + dplyr::filter(.data$state %in% term_state_inds, .data$trans_idx > 0) + if (nrow(term_rows) == 0) { + return(df) + } + term_events <- term_rows |> dplyr::group_by(.data$path_id) |> - summarise(term_time = min(.data$time, na.rm = T)) |> + summarise(term_time = min(.data$time, na.rm = TRUE)) |> dplyr::ungroup() no_terms <- df |> dplyr::anti_join(term_events, by = "path_id") diff --git a/inst/stan/msm.stan b/inst/stan/msm.stan index 389dd88..80ed5e0 100644 --- a/inst/stan/msm.stan +++ b/inst/stan/msm.stan @@ -86,7 +86,20 @@ functions { reject("t = ", t,", should be non-negative"); } real ke = CL / V2; - return(A0 * ka/(ka-ke) * (exp(-ke*t) - exp(-ka*t)) + C0 * exp(-ke*t)); + real d = ka - ke; + real scale = fmax(abs(ka), abs(ke)); + real absorbed; + if(abs(d) <= 1e-8 * scale){ + real k = 0.5 * (ka + ke); + absorbed = A0 * k * t * exp(-k*t); + } else if(d > 0){ + absorbed = A0 * ka * exp(-ke*t) * (-expm1(-d*t)) / d; + } else { + // When ka < ke, factoring by exp(-ke*t) can create 0 * Inf. + // Factor by the more slowly decaying exp(-ka*t) instead. + absorbed = A0 * ka * exp(-ka*t) * expm1(d*t) / d; + } + return(absorbed + C0 * exp(-ke*t)); } // Analytic solution with general initial condition A0 @@ -102,18 +115,23 @@ functions { real two_cpt_central_ss(real t, real tau, real dose, real ka, real CL, real V2) { real ke = CL / V2; - real A = dose * (ka / (ka-ke)); - real tt = fmod(t, tau); - real ma = A * inv(-expm1(-ka*tau)); // 1/(1-exp(-ka*tau)) - real me = A * inv(-expm1(-ke*tau)); - return(me * exp(-ke*tt) - ma * exp(-ka * tt)); + real d = ka - ke; + real tt = t - floor(t / tau) * tau; + real scale = fmax(abs(ka), abs(ke)); + if(abs(d) <= 1e-8 * scale){ + real k = 0.5 * (ka + ke); + real f = exp(-k*tt) * inv(-expm1(-k*tau)); + return(dose * k * f * (tt + tau / expm1(k*tau))); + } + return(dose * ka / d * ( + exp(-ke*tt) * inv(-expm1(-ke*tau)) + - exp(-ka*tt) * inv(-expm1(-ka*tau)) + )); } // Two-cpt PK model (steady state at trough) real two_cpt_central_ss0(real tau, real dose, real ka, real CL, real V2) { - real ke = CL / V2; - real A = dose * (ka / (ka-ke)); - return(A * (inv(-expm1(-ke*tau)) - inv(-expm1(-ka*tau)))); + return(two_cpt_central_ss(0.0, tau, dose, ka, CL, V2)); } // Two-cpt PK model (steady state at trough) @@ -294,7 +312,7 @@ data { int omit_lik_haz; // flag int omit_lik_pk; //flag int nc_haz; // number of hazard covariates - vector[N_trans] mu_w0; // Assumed mean h0 + vector[N_trans] mu_log_w0; // Prior location for log baseline hazard int N_grid; // number of integration grid points real delta_grid; // grid step size matrix[N_grid, N_sbf] SBF_grid; // basis functions evaluated at t_grid @@ -318,7 +336,7 @@ data { int D_trans; // max num of occurrences for a transition int D_risk; // max num of at-risk intervals for a transition array[N_trans] int sum_trans; // total number of occurred - array[N_trans] int sum_risk; // total number at risk + array[N_trans] int sum_risk; // total number at risk array[N_trans, D_trans] int which_trans; array[N_trans, D_risk] int which_risk; @@ -406,9 +424,12 @@ transformed parameters { // Baseline hazard if(do_haz == 1) { for(j in 1:N_trans){ - weights[1,j] = mu_weights[1][:,ttype[j]] + + vector[N_sbf] weights_raw = mu_weights[1][:,ttype[j]] + sig_weights[1][:,ttype[j]] .* z_weights[1][:,j]; - log_w0[1,j] = mu_w0[j] + sig_w0[1] * z_w0[1][j]; + // The intercept-inclusive B-spline basis partitions unity. Centering + // removes its otherwise exact shift redundancy with log_w0. + weights[1,j] = weights_raw - mean(weights_raw); + log_w0[1,j] = mu_log_w0[j] + sig_w0[1] * z_w0[1][j]; } } @@ -482,36 +503,28 @@ model { // Hazard model likelihood if (omit_lik_haz == 0 && do_haz == 1) { + // Invariant in h: compute once, rather than once per transition. + matrix[N_int, N_trans] log_C_haz = compute_log_hazard_multiplier( + N_int, beta_oth[1], beta_xpsr[1], x_haz_long, x_xpsr_long, xpsr_loc, xpsr_scale, + ttype + ); for(h in 1:N_trans){ - - // log of hazard multiplier on each interval - matrix[N_int, N_trans] log_C_haz = compute_log_hazard_multiplier( - N_int, beta_oth[1], beta_xpsr[1], x_haz_long, x_xpsr_long, xpsr_loc, xpsr_scale, - ttype - ); - - // Ragged array access - array[sum_risk[h]] int idx_atr = which_risk[h, 1:sum_risk[h]]; - array[sum_trans[h]] int idx_occ = which_trans[h, 1:sum_trans[h]]; - - // Occurred transitions (log hazard at interval end time) - target += log_hazard( - log_C_haz[idx_occ, h], SBF_end[idx_occ,:], weights[1,h], log_w0[1,h] - ); - - // Evaluate baseline hazard at grid, prepad with zero - vector[N_grid+1] h0_grid = rep_vector(0.0, N_grid+1); - h0_grid[2:(N_grid+1)] = exp(log_basehaz(SBF_grid, weights[1,h], log_w0[1,h])); - - // Integrated baseline hazard at grid - vector[N_grid+1] h0_int = cumulative_sum(h0_grid) * delta_grid; - - // Transitions that were at risk (- integrated hazard over interval) - // Taking into account that indices are off by one because of - // The prepadded zero - target += - exp(log_C_haz[idx_atr, h]) .* - (h0_int[t_end_idx[idx_atr]] - h0_int[t_start_idx_m1[idx_atr]]) .* - correction_multiplier[idx_atr]; + if(sum_trans[h] > 0){ + array[sum_trans[h]] int idx_occ = which_trans[h, 1:sum_trans[h]]; + target += log_hazard( + log_C_haz[idx_occ, h], SBF_end[idx_occ,:], weights[1,h], log_w0[1,h] + ); + } + if(sum_risk[h] > 0){ + array[sum_risk[h]] int idx_atr = which_risk[h, 1:sum_risk[h]]; + vector[N_grid+1] h0_grid = rep_vector(0.0, N_grid+1); + vector[N_grid+1] h0_int; + h0_grid[2:(N_grid+1)] = exp(log_basehaz(SBF_grid, weights[1,h], log_w0[1,h])); + h0_int = cumulative_sum(h0_grid) * delta_grid; + target += - exp(log_C_haz[idx_atr, h]) .* + (h0_int[t_end_idx[idx_atr]] - h0_int[t_start_idx_m1[idx_atr]]) .* + correction_multiplier[idx_atr]; + } } } diff --git a/man/MultistateSystem.Rd b/man/MultistateSystem.Rd index 9e550af..c2ab6c8 100644 --- a/man/MultistateSystem.Rd +++ b/man/MultistateSystem.Rd @@ -309,7 +309,8 @@ Max instant hazard on interval (t1, t2) \item{\code{w}}{Spline basis function weights (vector)} -\item{\code{log_w0}}{Intercept (log)} +\item{\code{log_w0}}{Intercept (log), or \code{-Inf} for an identically zero +baseline hazard} \item{\code{log_m}}{Hazard multiplier (log)} } diff --git a/tests/testthat/test-correctness-regressions.R b/tests/testthat/test-correctness-regressions.R new file mode 100644 index 0000000..a84bd42 --- /dev/null +++ b/tests/testthat/test-correctness-regressions.R @@ -0,0 +1,257 @@ +make_sparse_survival_paths <- function(event = TRUE) { + tm <- transmat_survival(c("Healthy", "Dead")) + subject_df <- data.frame(subject_id = c("a", "b")) + path_df <- data.frame( + path_id = c(1L, 1L, 2L, 2L), + state = c(1L, if (event) 2L else 1L, 1L, 1L), + time = c(0, 1, 0, 2), + trans_idx = c(0L, if (event) 1L else 0L, 0L, 0L) + ) + link_df <- data.frame( + path_id = 1:2, + subject_id = c("a", "b"), + draw_idx = 1L, + rep_idx = 1L + ) + PathData$new(subject_df, path_df, link_df, tm) +} + +test_that("baseline-hazard prior uses log scale in Stan data", { + pd <- make_sparse_survival_paths() + mod <- create_msm(pd$transmat, t_max = 2, n_grid = 10) + mod$set_prior_mean_h0_data(pd) + + expect_equal(mod$get_prior_mean_h0(), 1 / 3) + expect_equal(create_stan_data_model(mod)$mu_log_w0, log(1 / 3)) + expect_error(mod$set_prior_mean_h0(0), "strictly positive") +}) + +test_that("all-censored data get finite prior and padded Stan indices", { + pd <- make_sparse_survival_paths(event = FALSE) + mod <- create_msm(pd$transmat, t_max = 2, n_grid = 10) + mod$set_prior_mean_h0_data(pd) + sd <- create_stan_data(mod, pd) + + expect_equal(mod$get_prior_mean_h0(), 0.5 / 3) + expect_equal(sd$D_trans, 1) + expect_equal(sd$sum_trans, 0) + expect_equal(sd$sum_risk, 2) + + empty <- which_format_for_stan(matrix(0, nrow = 2, ncol = 3), "risk") + expect_equal(empty$D_risk, 1) + expect_equal(empty$sum_risk, c(0, 0)) +}) + +test_that("empirical hazard calibration pools within transition type", { + tm <- transmat_illnessdeath(c("Healthy", "Ill", "Dead")) + subject_df <- data.frame(subject_id = c("a", "b")) + path_df <- data.frame( + path_id = c(1L, 1L, 1L, 2L, 2L), + state = c(1L, 2L, 2L, 1L, 1L), + time = c(0, 1, 3, 0, 2), + trans_idx = c(0L, 1L, 0L, 0L, 0L) + ) + link_df <- data.frame( + path_id = 1:2, + subject_id = c("a", "b"), + draw_idx = 1L, + rep_idx = 1L + ) + pd <- PathData$new(subject_df, path_df, link_df, tm) + rates <- average_haz_per_ttype(pd) |> + dplyr::arrange(.data$trans_idx) + + expect_equal(exp(rates$log_h0_avg), c(1 / 3, 0.5 / 5, 0.5 / 5)) +}) + +test_that("an unexposed transition type borrows the global pooled rate", { + tm <- TransitionMatrix$new( + matrix(c( + 0, 1, 0, 0, + 0, 0, 0, 0, + 0, 0, 0, 1, + 0, 0, 0, 0 + ), nrow = 4, byrow = TRUE), + c("A", "B", "C", "D") + ) + subject_df <- data.frame(subject_id = c("a", "b")) + path_df <- data.frame( + path_id = c(1L, 1L, 2L, 2L), + state = c(1L, 2L, 1L, 1L), + time = c(0, 1, 0, 2), + trans_idx = c(0L, 1L, 0L, 0L) + ) + link_df <- data.frame( + path_id = 1:2, + subject_id = c("a", "b"), + draw_idx = 1L, + rep_idx = 1L + ) + pd <- PathData$new(subject_df, path_df, link_df, tm) + rates <- average_haz_per_ttype(pd) |> + dplyr::arrange(.data$trans_idx) + + expect_equal(exp(rates$log_h0_avg), c(1 / 3, 1 / 3)) + + mod <- create_msm(tm, t_max = 2, n_grid = 10) + mod$set_prior_mean_h0_data(pd) + sd <- create_stan_data(mod, pd) + expect_equal(sd$sum_trans, c(1, 0)) + expect_equal(sd$sum_risk, c(2, 0)) + expect_equal(sd$D_trans, 1) +}) + +test_that("state-visit probabilities retain zero-event subjects", { + pd <- make_sparse_survival_paths() + p <- p_state_visit(pd, by = "subject_id") |> + dplyr::arrange(.data$subject_id) + + expect_equal(p$subject_id, c("a", "b")) + expect_equal(p$prob, c(1, 0)) + expect_equal(p$n_event, c(1, 0)) +}) + +test_that("posterior hazard draws remain draw-major", { + w <- array(1:3, dim = c(3, 1, 1)) + log_w0 <- matrix(11:13, nrow = 3) + out <- repeat_hazard_parameter_draws(w, log_w0, n_subjects = 2) + + expect_equal(as.vector(out$w), c(1, 1, 2, 2, 3, 3)) + expect_equal(as.vector(out$log_w0), c(11, 11, 12, 12, 13, 13)) +}) + +test_that("noninteger time grids contain exactly the requested midpoints", { + mod <- create_msm(transmat_survival(), t_max = 365.25, n_grid = 1000) + grid <- create_stan_data_timegrid(mod) + + expect_length(grid$t_grid, 1000) + expect_equal(grid$N_grid, 1000) + expect_true(all(grid$t_grid > 0 & grid$t_grid < mod$get_tmax())) + expect_equal(diff(range(grid$t_grid)), 999 * grid$delta_grid) +}) + +test_that("transition probabilities honor a delayed start time", { + sys <- MultistateSystem$new(transmat_survival()) + sys$set_knots(c(0, 10)) + lambda <- 0.2 + p <- solve_trans_prob_matrix( + sys, t_out = c(6, 7), log_w0 = log(lambda), t_start = 5 + ) + + expect_equal(p[, 1, 1], exp(-lambda * c(1, 2)), tolerance = 1e-6) + expect_equal(p[, 1, 2], 1 - exp(-lambda * c(1, 2)), tolerance = 1e-6) + expect_equal( + solve_trans_prob_matrix(sys, t_out = 5, log_w0 = log(lambda), t_start = 5)[1, , ], + diag(2) + ) + expect_error( + solve_trans_prob_matrix(sys, t_out = c(7, 6), log_w0 = log(lambda), t_start = 5), + "strictly increasing" + ) +}) + +test_that("B-spline thinning bound dominates the evaluated hazard", { + sys <- MultistateSystem$new(transmat_survival()) + sys$set_knots(c(0, 0.4999, 0.5001, 1)) + w <- c(-3, 3, -3, 3, -3) + bound <- sys$max_inst_hazard(0, 1, w, log_w0 = -2, log_m = 0.4) + dense <- exp(sys$log_inst_hazard(seq(0, 1, length.out = 10001), w, -2, 0.4)) + + expect_gte(bound, max(dense)) + expect_gt(bound, exp(-2 + 0.4 + max(w))) +}) + +test_that("B-spline thinning bound is local to the requested interval", { + sys <- MultistateSystem$new(transmat_survival()) + sys$set_knots(c(0, 0.4999, 0.5001, 1)) + w <- c(-10, 10, -10, -10, -10) + times <- seq(0.8, 1, length.out = 10001) + + local_bound <- sys$max_inst_hazard( + 0.8, 1, w, log_w0 = 0, log_m = 0 + ) + global_bound <- exp(max(w)) + dense <- exp(sys$log_inst_hazard(times, w, 0, 0)) + + expect_gte(local_bound, max(dense)) + expect_lt(local_bound, global_bound / 1e6) +}) + +test_that("zero baseline hazard has zero bound and generates no event", { + sys <- MultistateSystem$new(transmat_survival()) + sys$set_knots(c(0, 1)) + w <- array(0, dim = c(1, 1, sys$num_weights())) + + expect_equal(sys$max_inst_hazard(0, 1, w[1, 1, ], -Inf, 0), 0) + + path <- suppressMessages(sys$simulate( + w = w, + log_w0 = matrix(-Inf, nrow = 1), + log_m = matrix(0, nrow = 1), + t_max = 1 + )) + expect_equal(path$time, c(0, 1)) + expect_false(any(path$is_event == 1)) +}) + +test_that("steady-state PK solution is finite when ka equals ke", { + pk <- PKModel$new(list(ka = NULL, CL = NULL, V2 = NULL)) + times <- c(0, 1, 6, 12) + theta_equal <- list(ka = 0.2, CL = 0.2, V2 = 1) + theta_near <- list(ka = 0.2, CL = 0.2 * (1 + 1e-10), V2 = 1) + x <- pk$simulate_ss(times, theta_equal, dose = 30, tau = 12) + y <- pk$simulate_ss(times, theta_near, dose = 30, tau = 12) + + expect_true(all(is.finite(x) & x >= 0)) + expect_equal(x, y, tolerance = 1e-7) +}) + +test_that("small but distinct PK rates are not treated as equal", { + pk <- PKModel$new(list(ka = NULL, CL = NULL, V2 = NULL)) + times <- c(0, 1, 6, 12) + dose <- 30 + tau <- 12 + ka <- 1e-9 + ke <- 2e-9 + tt <- times %% tau + f_ke <- exp(-ke * tt) / (-expm1(-ke * tau)) + f_ka <- exp(-ka * tt) / (-expm1(-ka * tau)) + expected <- dose * ka / (ka - ke) * (f_ke - f_ka) + actual <- pk$simulate_ss( + times, + list(ka = ka, CL = ke, V2 = 1), + dose = dose, + tau = tau + ) + + expect_equal(actual, expected, tolerance = 1e-8) +}) + +test_that("requested PK interval width uses the correct quantiles", { + sim <- data.frame( + subject_id = "a", + time = 1, + val = seq_len(100) + ) + out <- pksim_to_quantiles(sim, ci_alpha = 0.8) + expect_equal(out$lower, unname(stats::quantile(sim$val, 0.1))) + expect_equal(out$upper, unname(stats::quantile(sim$val, 0.9))) +}) + +test_that("base data frames work without covariates", { + df <- data.frame( + subject_id = rep(c("a", "b"), each = 2), + time = c(0, 1, 0, 2), + state = c(1L, 2L, 1L, 1L), + is_transition = c(FALSE, TRUE, FALSE, FALSE) + ) + pd <- df_to_pathdata(df, transmat_survival()) + expect_s3_class(pd, "PathData") + expect_equal(pd$n_paths(), 2) +}) + +test_that("exposure normalization scale must be positive", { + mod <- create_msm(transmat_survival()) + expect_error(mod$set_xpsr_normalizers(0, 0), "strictly positive") + expect_message(mod$set_xpsr_normalizers(-1, 0.5), "loc = -1") +}) diff --git a/tests/testthat/test-stan.R b/tests/testthat/test-stan.R index 569e603..8b5d39a 100644 --- a/tests/testthat/test-stan.R +++ b/tests/testthat/test-stan.R @@ -69,3 +69,16 @@ test_that("fitting with Stan works (multi-transition)", { pp <- plot_state_occupancy(r) expect_true(is_ggplot(pp)) }) + +test_that("transient PK solution is stable when elimination is much faster", { + ensure_exposed_stan_functions() + ka <- 0.1 + ke <- 100 + t <- 10 + A0 <- 1 + expected <- A0 * ka / (ka - ke) * (exp(-ke * t) - exp(-ka * t)) + actual <- two_cpt_central(t, ka, ke, 1, A0, 0) + + expect_true(is.finite(actual)) + expect_equal(actual, expected, tolerance = 1e-12) +}) From 267209a43ca2e42919e52fb9b9ce8abaf2f68dba Mon Sep 17 00:00:00 2001 From: jtimonen Date: Tue, 28 Jul 2026 12:04:37 +0300 Subject: [PATCH 3/7] v0.4.0 rc --- .Rbuildignore | 1 - DESCRIPTION | 2 +- NEWS.md | 118 +++++++++++++++--- R/MultistateModel.R | 20 +-- R/MultistateModelFit.R | 30 ++--- R/MultistateSystem.R | 51 ++++++-- R/PKModel.R | 2 +- R/stan.R | 8 +- inst/stan/msm.stan | 82 ++++++------ man/MultiStateModel.Rd | 18 +-- man/MultistateSystem.Rd | 6 +- man/fit_stan.Rd | 6 +- tests/testthat/test-correctness-regressions.R | 92 ++++++++++++++ tests/testthat/test-stan.R | 58 +++++++++ vignettes/math.Rmd | 20 +-- 15 files changed, 400 insertions(+), 114 deletions(-) diff --git a/.Rbuildignore b/.Rbuildignore index 60657a7..808480e 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,7 +5,6 @@ ^\_config\.yml$ ^\.github$ ^CONTRIBUTING\.md$ -^NEWS\.md$ ^_pkgdown\.yml$ pkgdown/ .DS_Store diff --git a/DESCRIPTION b/DESCRIPTION index ffffec2..6fb62ef 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: bmstate Type: Package Title: Bayesian multistate modeling -Version: 0.3.3.9001 +Version: 0.4.0 Authors@R: c(person(given = "Juho", family = "Timonen", diff --git a/NEWS.md b/NEWS.md index afdd900..ab7eb4c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,14 +1,104 @@ -# bmstate 0.3.3.9001 - -* Correct the empirical baseline-hazard prior scale and use robust event-time - exposure estimates for sparse transitions. -* Preserve joint posterior-draw alignment in path and occupancy prediction. -* Return explicit zero state-visit probabilities for all requested groups. -* Use a guaranteed B-spline thinning bound and remove rare-event truncation. -* Stabilize the oral one-compartment PK solution when absorption and elimination - rates are equal or nearly equal. -* Fix noninteger integration grids, delayed transition-probability start times, - requested PK credible-interval widths, and zero-event Stan data. -* Center spline weights to remove their shift redundancy with the baseline - log-hazard intercept, and avoid recomputing invariant hazard multipliers in - the Stan transition loop. +# 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 supplied to `fit_stan(filepath = ...)` must adopt the new name. + Fits produced by affected versions must be rerun. + +* **Make sparse-transition prior calibration finite and explicit.** The + previous cumulative-hazard summary omitted transitions with no events and + could leave their empirical prior locations undefined. Reference rates are + now computed from event counts divided by exact at-risk time, pooled within + transition type. An exposed type with no events uses a half-event + regularization, and an entirely unexposed type borrows the global pooled + rate. Documentation now calls this 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. + +* **Make event-time thinning valid for flexible spline hazards.** The former + upper bound was the largest hazard on a 100-point grid, which could miss a + narrow spline peak and invalidate thinning. It also discarded transitions + whose bound was at most `1e-9`, introducing rare-event truncation. The new + bound uses non-negativity, partition of unity, and local support of the + B-spline basis to dominate the hazard over the complete requested interval; + all transitions with a positive bound remain eligible. + +* **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. + +* **Support transitions with zero events or zero at-risk intervals in Stan + data.** Zero counts previously produced zero-width ragged arrays while the + Stan declarations required positive dimensions and at least one risk row. + R now pads storage to width one while retaining a count of zero; Stan accepts + zero counts and skips the corresponding likelihood contribution. + +* **Separate baseline level from spline shape.** Because the + intercept-inclusive B-spline basis partitions unity, adding a constant to + every spline coefficient and subtracting it from `log_w0` left the likelihood + unchanged. Evaluated spline weights are now centered to sum to zero, so the + baseline intercept alone carries the constant level. The full-dimensional + raw hierarchy still contains prior-only common-shift directions; this change + removes the likelihood confounding, not every latent nuisance direction. + +* **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. diff --git a/R/MultistateModel.R b/R/MultistateModel.R index eb21545..dc3f158 100644 --- a/R/MultistateModel.R +++ b/R/MultistateModel.R @@ -177,19 +177,23 @@ 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, finite = TRUE) @@ -200,8 +204,8 @@ MultistateModel <- R6::R6Class("MultistateModel", 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")) { diff --git a/R/MultistateModelFit.R b/R/MultistateModelFit.R index 3f8f4f7..bb3387a 100644 --- a/R/MultistateModelFit.R +++ b/R/MultistateModelFit.R @@ -511,14 +511,6 @@ 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) @@ -526,7 +518,7 @@ msmfit_log_hazard_multipliers <- function(fit, oos = FALSE, data = NULL) { 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 @@ -534,17 +526,17 @@ msmfit_log_hazard_multipliers <- function(fit, oos = FALSE, data = 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 } @@ -746,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 diff --git a/R/MultistateSystem.R b/R/MultistateSystem.R index 8f8dfbc..295111a 100644 --- a/R/MultistateSystem.R +++ b/R/MultistateSystem.R @@ -16,6 +16,15 @@ MultistateSystem <- R6::R6Class("MultistateSystem", w <- matrix(w, nrow = 1) } checkmate::assert_array(w, d = 2) + if (t_start == t_max) { + return(cbind( + time = t_start, + state = init_state, + is_event = 0, + is_censor = 1, + trans_idx = 0 + )) + } # Setup states <- init_state @@ -33,8 +42,27 @@ MultistateSystem <- R6::R6Class("MultistateSystem", states[j], t_cur, t_max, w, log_w0, log_m ) t_next <- trans$t - if (t_next - t_cur < min_t_step) { - t_next <- t_cur + min_t_step + if (trans$idx > 0) { + original_t_next <- t_next + if (t_next <= t_cur) { + # Ensure representable forward progress even when an exponential + # waiting time rounds to zero and min_t_step is zero. + numerical_step <- 2 * .Machine$double.eps * max(1, abs(t_cur)) + t_next <- t_cur + max(min_t_step, numerical_step) + } else if (t_next - t_cur < min_t_step) { + t_next <- t_cur + min_t_step + } + if (t_next >= t_max) { + if (original_t_next > t_cur && original_t_next < t_max) { + # Do not move a genuine pre-horizon event to or beyond censoring. + t_next <- original_t_next + } else { + # No representable interior event time remains. + t_next <- t_max + trans$idx <- 0 + trans$new_state <- 0 + } + } } # Update time and state @@ -57,6 +85,7 @@ MultistateSystem <- R6::R6Class("MultistateSystem", is_event[1] <- 0 # initial state is never an event is_event[L] <- 0 # last state is never an event (it is censoring time) states[L] <- states[L - 1] + tidx[L] <- 0 cbind( time = times, state = states, @@ -111,12 +140,6 @@ MultistateSystem <- R6::R6Class("MultistateSystem", new_state <- 0 } - # Safeguard against infinite loop - dt_min <- 1e-9 - if (t_min_found - t_init < dt_min) { - t_min_found <- t_init + dt_min - } - # Return list( t = t_min_found, new_state = new_state, idx = trans_idx @@ -134,7 +157,7 @@ MultistateSystem <- R6::R6Class("MultistateSystem", n_tries <- n_tries + 1 u1 <- stats::runif(1) t <- t - 1 / lambda_ub * log(u1) - if (t > t_max) { + if (t >= t_max) { return(list(t = t_max, n_tries = n_tries, censor = TRUE)) } lambda_t <- exp(self$log_inst_hazard(t, w, log_w0, log_m)) @@ -404,7 +427,8 @@ MultistateSystem <- R6::R6Class("MultistateSystem", #' @param t_max Max time. If \code{NULL}, the max #' time of the model is used. #' @param n_rep Number of repetitions to do for each draw. - #' @param min_t_step Minimal time step. + #' @param min_t_step Minimum separation applied to consecutive simulated + #' event times when this does not move an event to or beyond \code{t_max}. #' @return A data frame with \code{n_draws} x \code{n_rep} paths. simulate = function(w, log_w0, log_m, init_state = 1, t_start = 0, t_max = NULL, n_rep = 1, min_t_step = 1e-6) { @@ -431,7 +455,9 @@ MultistateSystem <- R6::R6Class("MultistateSystem", if (is.null(t_max)) { t_max <- self$get_tmax() } - checkmate::assert_number(t_max, lower = 0) + checkmate::assert_number( + t_max, lower = t_start, upper = self$get_tmax(), finite = TRUE + ) # Should not be done in parallel as such because can mess order in link df cnt <- 0 @@ -535,6 +561,9 @@ solve_trans_prob_matrix <- function(system, t_out, log_w0, w = NULL, if (is.null(log_m)) { log_m <- rep(0, H) } + checkmate::assert_matrix(w, ncols = W, nrows = H) + checkmate::assert_numeric(log_w0, len = H) + checkmate::assert_numeric(log_m, len = H) solve_times <- unique(c(t_start, t_out)) if (length(solve_times) == 1) { P <- array(diag(1, S, S), dim = c(1, S, S)) diff --git a/R/PKModel.R b/R/PKModel.R index 5369fda..290070e 100644 --- a/R/PKModel.R +++ b/R/PKModel.R @@ -105,7 +105,7 @@ PKModel <- R6::R6Class("PKModel", if (abs(d) <= rel_tol) { k <- 0.5 * (ka + ke) f <- exp(-k * tt) / (-expm1(-k * tau)) - return((dose / theta$V2) * k * f * (tt + tau / expm1(k * tau))) + return((dose / theta$V2) * ka * f * (tt + tau / expm1(k * tau))) } f_ke <- exp(-ke * tt) / (-expm1(-ke * tau)) f_ka <- exp(-ka * tt) / (-expm1(-ka * tau)) diff --git a/R/stan.R b/R/stan.R index 526ff3e..7720669 100644 --- a/R/stan.R +++ b/R/stan.R @@ -50,7 +50,7 @@ ensure_exposed_stan_functions <- function(...) { #' #' @description #' \emph{NOTE:} This function has a side effect of setting covariate -#' normalizers, prior assumed mean baseline hazard, and concentration upper +#' normalizers, baseline-hazard prior reference rates, and concentration upper #' bound (PK) based on data. #' #' @export @@ -59,8 +59,8 @@ ensure_exposed_stan_functions <- function(...) { #' @param return_stanfit Return also the raw 'Stan' fit object? #' @inheritParams create_stan_model #' @param set_normalizers Set covariate normalization automatically? -#' @param set_prior_h0 Set prior mean average baseline hazard levels based -#' on data? +#' @param set_prior_h0 Set baseline-hazard prior reference rates from event +#' counts and at-risk time in the data? #' @param max_conc_factor Factor to multiply observed max concentration by #' to get concentration upper bound. #' @param method Must be one of \code{"sample"} (default), @@ -92,7 +92,7 @@ fit_stan <- function(model, data, # Get Stan model object stan_model <- create_stan_model(filepath = filepath) - # Set prior mean baseline hazard rates (side effect) + # Set baseline-hazard prior reference rates (side effect) if (set_prior_h0) { model$set_prior_mean_h0_data(data) } diff --git a/inst/stan/msm.stan b/inst/stan/msm.stan index 80ed5e0..eb11a1b 100644 --- a/inst/stan/msm.stan +++ b/inst/stan/msm.stan @@ -7,29 +7,27 @@ functions { // log of hazard multiplier matrix compute_log_hazard_multiplier( - data int N_int, + data int N_row, + data int N_trans_types, array[] vector beta_oth, array[] vector beta_xpsr, data array[] vector x_haz, array[] vector x_xpsr, real xpsr_loc, - real xpsr_scale, - data array[] int ttype + real xpsr_scale ) { - int N_trans = size(ttype); int nc_haz = size(x_haz); - matrix[N_int, N_trans] log_C_haz = rep_matrix(0.0, N_int, N_trans); - for(j in 1:N_trans){ - int h = ttype[j]; - if(nc_haz > 0){ - for(k in 1:nc_haz){ - log_C_haz[,j] += beta_oth[k][h] * x_haz[k]; - } - } - if(size(beta_xpsr)==1){ - log_C_haz[,j] += beta_xpsr[1][h] * ((x_xpsr[1] - xpsr_loc) / xpsr_scale); + matrix[N_row, N_trans_types] log_C_haz = + rep_matrix(0.0, N_row, N_trans_types); + if(nc_haz > 0){ + for(k in 1:nc_haz){ + log_C_haz += x_haz[k] * beta_oth[k]'; } } + if(size(beta_xpsr)==1){ + log_C_haz += ((x_xpsr[1] - xpsr_loc) / xpsr_scale) * + beta_xpsr[1]'; + } return(log_C_haz); } @@ -89,9 +87,11 @@ functions { real d = ka - ke; real scale = fmax(abs(ka), abs(ke)); real absorbed; + // Parameter-dependent numerical guard. The branch approximation is + // first-order matched, so its value and gradient agree at ka = ke. if(abs(d) <= 1e-8 * scale){ real k = 0.5 * (ka + ke); - absorbed = A0 * k * t * exp(-k*t); + absorbed = A0 * ka * t * exp(-k*t); } else if(d > 0){ absorbed = A0 * ka * exp(-ke*t) * (-expm1(-d*t)) / d; } else { @@ -118,10 +118,11 @@ functions { real d = ka - ke; real tt = t - floor(t / tau) * tau; real scale = fmax(abs(ka), abs(ke)); + // As above, this guard is first-order matched at the equal-rate limit. if(abs(d) <= 1e-8 * scale){ real k = 0.5 * (ka + ke); real f = exp(-k*tt) * inv(-expm1(-k*tau)); - return(dose * k * f * (tt + tau / expm1(k*tau))); + return(dose * ka * f * (tt + tau / expm1(k*tau))); } return(dose * ka / d * ( exp(-ke*tt) * inv(-expm1(-ke*tau)) @@ -280,7 +281,6 @@ functions { // Find drug amounts in both compartments at dose_times int N_t = num_elements(t[1]); - vector[N_t] max_conc = rep_vector(MAX_CONC, N_t); int N_sub = num_elements(dose_ss); int D = num_elements(dose_times[1]); array[2, N_sub] vector[D] amounts = pop_2cpt_partly_ss_stage1( @@ -373,18 +373,28 @@ data { transformed data { array[N_sub] vector[1] t0_ss; - for(n in 1:N_sub){ - t0_ss[n] = rep_vector(0.0, 1); + if(I_xpsr != do_pk){ + reject("I_xpsr must equal do_pk"); } - // Set x corresponding to each interval - array[nc_haz] vector[N_int] x_haz_long; - if(nc_haz > 0){ - for(j in 1:nc_haz){ - for(n in 1:N_int){ - x_haz_long[j][n] = x_haz[j][idx_sub[n]]; + if(I_xpsr == 1 && xpsr_scale <= 0){ + reject("xpsr_scale must be strictly positive when exposure is modeled"); + } + if(do_pk == 1){ + if(tau_ss <= 0){ + reject("tau_ss must be strictly positive when PK is modeled"); + } + if(MAX_CONC <= 0){ + reject("MAX_CONC must be strictly positive when PK is modeled"); + } + for(n in 1:N_sub){ + if(dose_ss[n] <= 0){ + reject("dose_ss must be strictly positive when PK is modeled"); } } } + for(n in 1:N_sub){ + t0_ss[n] = rep_vector(0.0, 1); + } } parameters { @@ -437,7 +447,6 @@ transformed parameters { array[do_pk] matrix[N_sub, 3] log_theta_pk; array[do_pk, N_sub] vector[2] conc_mu_pk; array[do_pk] vector[N_sub] ss_xpsr; // xpsr for each subject - array[do_pk] vector[N_int] x_xpsr_long; // xpsr for each interval if(do_pk == 1){ @@ -455,11 +464,6 @@ transformed parameters { // Concentration xpsr at steady state ss_xpsr[1] = log_ss_area_under_conc(dose_ss, log_theta_pk[1]); - - // Set xpsr corresponding to each interval - for(n in 1:N_int){ - x_xpsr_long[1][n] = ss_xpsr[1][idx_sub[n]]; - } } } @@ -489,7 +493,7 @@ model { // PK PARAM PRIOR if(do_pk == 1){ - for(n in 1:size(log_z_pk[1])){ + for(n in 1:N_sub){ log_z_pk[1, n] ~ normal(0, 1); } log_mu_pk[1] ~ normal(0, 2); @@ -503,16 +507,18 @@ model { // Hazard model likelihood if (omit_lik_haz == 0 && do_haz == 1) { - // Invariant in h: compute once, rather than once per transition. - matrix[N_int, N_trans] log_C_haz = compute_log_hazard_multiplier( - N_int, beta_oth[1], beta_xpsr[1], x_haz_long, x_xpsr_long, xpsr_loc, xpsr_scale, - ttype + // Covariates are time invariant, so differentiate one predictor per + // subject and transition type, then gather it for each interval. + matrix[N_sub, N_trans_types] log_C_haz = compute_log_hazard_multiplier( + N_sub, N_trans_types, beta_oth[1], beta_xpsr[1], x_haz, ss_xpsr, + xpsr_loc, xpsr_scale ); for(h in 1:N_trans){ if(sum_trans[h] > 0){ array[sum_trans[h]] int idx_occ = which_trans[h, 1:sum_trans[h]]; target += log_hazard( - log_C_haz[idx_occ, h], SBF_end[idx_occ,:], weights[1,h], log_w0[1,h] + log_C_haz[idx_sub[idx_occ], ttype[h]], SBF_end[idx_occ,:], + weights[1,h], log_w0[1,h] ); } if(sum_risk[h] > 0){ @@ -521,7 +527,7 @@ model { vector[N_grid+1] h0_int; h0_grid[2:(N_grid+1)] = exp(log_basehaz(SBF_grid, weights[1,h], log_w0[1,h])); h0_int = cumulative_sum(h0_grid) * delta_grid; - target += - exp(log_C_haz[idx_atr, h]) .* + target += - exp(log_C_haz[idx_sub[idx_atr], ttype[h]]) .* (h0_int[t_end_idx[idx_atr]] - h0_int[t_start_idx_m1[idx_atr]]) .* correction_multiplier[idx_atr]; } diff --git a/man/MultiStateModel.Rd b/man/MultiStateModel.Rd index 5eaa27e..2ae9c68 100644 --- a/man/MultiStateModel.Rd +++ b/man/MultiStateModel.Rd @@ -126,21 +126,23 @@ Set normalization constants for exposure (side effect) \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-MultistateModel-get_prior_mean_h0}{}}} \subsection{Method \code{get_prior_mean_h0()}}{ -Get assumed prior mean baseline hazard rates. +Get baseline-hazard prior reference rates. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{MultistateModel$get_prior_mean_h0()}\if{html}{\out{
}} } \subsection{Returns}{ -Numeric vector with length equal to number of transitions +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. } } \if{html}{\out{
}} \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-MultistateModel-set_prior_mean_h0}{}}} \subsection{Method \code{set_prior_mean_h0()}}{ -Set assumed prior mean baseline hazard rates (side -effect). +Set baseline-hazard prior reference rates (side effect). \subsection{Usage}{ \if{html}{\out{
}}\preformatted{MultistateModel$set_prior_mean_h0(mean_h0)}\if{html}{\out{
}} } @@ -148,7 +150,9 @@ effect). \subsection{Arguments}{ \if{html}{\out{
}} \describe{ -\item{\code{mean_h0}}{Numeric vector with length equal to number of transitions} +\item{\code{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.} } \if{html}{\out{
}} } @@ -157,8 +161,8 @@ effect). \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-MultistateModel-set_prior_mean_h0_data}{}}} \subsection{Method \code{set_prior_mean_h0_data()}}{ -Set assumed prior mean baseline hazard rates (side -effect) based on average hazards in data. +Set baseline-hazard prior reference rates (side effect) +using event counts and at-risk time pooled within transition type. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{MultistateModel$set_prior_mean_h0_data(data)}\if{html}{\out{
}} } diff --git a/man/MultistateSystem.Rd b/man/MultistateSystem.Rd index c2ab6c8..47a4327 100644 --- a/man/MultistateSystem.Rd +++ b/man/MultistateSystem.Rd @@ -243,7 +243,8 @@ Evaluate log instant hazard \item{\code{w}}{Spline basis function weights (vector)} -\item{\code{log_w0}}{Intercept (log)} +\item{\code{log_w0}}{Intercept (log), or \code{-Inf} for an identically zero +baseline hazard} \item{\code{log_m}}{Hazard multiplier (log)} @@ -355,7 +356,8 @@ time of the model is used.} \item{\code{n_rep}}{Number of repetitions to do for each draw.} -\item{\code{min_t_step}}{Minimal time step.} +\item{\code{min_t_step}}{Minimum separation applied to consecutive simulated +event times when this does not move an event to or beyond \code{t_max}.} } \if{html}{\out{}} } diff --git a/man/fit_stan.Rd b/man/fit_stan.Rd index 8ba7968..1545611 100644 --- a/man/fit_stan.Rd +++ b/man/fit_stan.Rd @@ -23,8 +23,8 @@ fit_stan( \item{set_normalizers}{Set covariate normalization automatically?} -\item{set_prior_h0}{Set prior mean average baseline hazard levels based -on data?} +\item{set_prior_h0}{Set baseline-hazard prior reference rates from event +counts and at-risk time in the data?} \item{filepath}{Deprecated.} @@ -45,7 +45,7 @@ A \code{\link{MultistateModelFit}} object. } \description{ \emph{NOTE:} This function has a side effect of setting covariate -normalizers, prior assumed mean baseline hazard, and concentration upper +normalizers, baseline-hazard prior reference rates, and concentration upper bound (PK) based on data. } \seealso{ diff --git a/tests/testthat/test-correctness-regressions.R b/tests/testthat/test-correctness-regressions.R index a84bd42..938956d 100644 --- a/tests/testthat/test-correctness-regressions.R +++ b/tests/testthat/test-correctness-regressions.R @@ -150,6 +150,28 @@ test_that("transition probabilities honor a delayed start time", { ) }) +test_that("identity transition probabilities still validate hazard parameters", { + sys <- MultistateSystem$new(transmat_survival()) + sys$set_knots(c(0, 10)) + log_lambda <- log(0.2) + + expect_error( + solve_trans_prob_matrix(sys, 5, rep(log_lambda, 2), t_start = 5) + ) + expect_error( + solve_trans_prob_matrix( + sys, + 5, + log_lambda, + w = matrix(0, nrow = 2, ncol = sys$num_weights()), + t_start = 5 + ) + ) + expect_error( + solve_trans_prob_matrix(sys, 5, log_lambda, log_m = c(0, 0), t_start = 5) + ) +}) + test_that("B-spline thinning bound dominates the evaluated hazard", { sys <- MultistateSystem$new(transmat_survival()) sys$set_knots(c(0, 0.4999, 0.5001, 1)) @@ -194,6 +216,53 @@ test_that("zero baseline hazard has zero bound and generates no event", { expect_false(any(path$is_event == 1)) }) +test_that("path simulation respects its horizon and censor-row contract", { + sys <- MultistateSystem$new(transmat_survival()) + sys$set_knots(c(0, 1)) + w <- array(0, dim = c(1, 1, sys$num_weights())) + log_w0 <- matrix(log(100), nrow = 1) + log_m <- matrix(0, nrow = 1) + + set.seed(1) + path <- suppressMessages(sys$simulate( + w, + log_w0, + log_m, + t_start = 0, + t_max = 0.5, + min_t_step = 1 + )) + expect_lte(max(path$time), 0.5) + expect_equal(tail(path$time, 1), 0.5) + expect_equal(tail(path$is_event, 1), 0) + expect_equal(tail(path$is_censor, 1), 1) + expect_equal(tail(path$trans_idx, 1), 0) + + zero_horizon <- suppressMessages(sys$simulate( + w, + log_w0, + log_m, + t_start = 0.5, + t_max = 0.5 + )) + expect_equal(nrow(zero_horizon), 1) + expect_equal(zero_horizon$time, 0.5) + expect_equal(zero_horizon$state, 1) + expect_equal(zero_horizon$is_event, 0) + expect_equal(zero_horizon$is_censor, 1) + expect_equal(zero_horizon$trans_idx, 0) + + expect_error( + suppressMessages(sys$simulate( + w, + log_w0, + log_m, + t_start = 0.6, + t_max = 0.5 + )) + ) +}) + test_that("steady-state PK solution is finite when ka equals ke", { pk <- PKModel$new(list(ka = NULL, CL = NULL, V2 = NULL)) times <- c(0, 1, 6, 12) @@ -206,6 +275,29 @@ test_that("steady-state PK solution is finite when ka equals ke", { expect_equal(x, y, tolerance = 1e-7) }) +test_that("near-equal steady-state PK limit retains the first-order ka term", { + pk <- PKModel$new(list(ka = NULL, CL = NULL, V2 = NULL)) + times <- c(0, 1, 6, 12) + dose <- 30 + tau <- 12 + ke <- 0.2 + + for (sign in c(-1, 1)) { + ka <- ke * (1 + sign * 5e-9) + k <- 0.5 * (ka + ke) + tt <- times %% tau + f <- exp(-k * tt) / (-expm1(-k * tau)) + expected <- dose * ka * f * (tt + tau / expm1(k * tau)) + actual <- pk$simulate_ss( + times, + list(ka = ka, CL = ke, V2 = 1), + dose = dose, + tau = tau + ) + expect_equal(actual, expected, tolerance = 1e-13) + } +}) + test_that("small but distinct PK rates are not treated as equal", { pk <- PKModel$new(list(ka = NULL, CL = NULL, V2 = NULL)) times <- c(0, 1, 6, 12) diff --git a/tests/testthat/test-stan.R b/tests/testthat/test-stan.R index 8b5d39a..fcbe9eb 100644 --- a/tests/testthat/test-stan.R +++ b/tests/testthat/test-stan.R @@ -82,3 +82,61 @@ test_that("transient PK solution is stable when elimination is much faster", { expect_true(is.finite(actual)) expect_equal(actual, expected, tolerance = 1e-12) }) + +test_that("Stan PK equal-rate guards match the first-order continuous limits", { + ensure_exposed_stan_functions() + ke <- 0.2 + t <- 6 + tau <- 12 + dose <- 30 + A0 <- 2 + C0 <- 0.4 + + for (sign in c(-1, 1)) { + ka <- ke * (1 + sign * 5e-9) + k <- 0.5 * (ka + ke) + + expected_transient <- + A0 * ka * t * exp(-k * t) + C0 * exp(-ke * t) + expect_equal( + two_cpt_central(t, ka, ke, 1, A0, C0), + expected_transient, + tolerance = 1e-13 + ) + + tt <- t %% tau + f <- exp(-k * tt) / (-expm1(-k * tau)) + expected_ss <- dose * ka * f * (tt + tau / expm1(k * tau)) + expect_equal( + two_cpt_central_ss(t, tau, dose, ka, ke, 1), + expected_ss, + tolerance = 1e-13 + ) + } +}) + +test_that("subject-level Stan hazard multipliers match direct calculation", { + ensure_exposed_stan_functions() + x_haz <- list(c(-1, 0.5, 2), c(0.25, -0.5, 1)) + beta_oth <- list(c(0.4, -0.2), c(-0.1, 0.3)) + x_xpsr <- list(c(1, 2, 4)) + beta_xpsr <- list(c(0.2, -0.4)) + loc <- 2 + scale <- 0.5 + + expected <- + x_haz[[1]] %o% beta_oth[[1]] + + x_haz[[2]] %o% beta_oth[[2]] + + ((x_xpsr[[1]] - loc) / scale) %o% beta_xpsr[[1]] + actual <- compute_log_hazard_multiplier( + 3, + 2, + beta_oth, + beta_xpsr, + x_haz, + x_xpsr, + loc, + scale + ) + expect_equal(actual, expected, tolerance = 1e-14) +}) diff --git a/vignettes/math.Rmd b/vignettes/math.Rmd index 1704e82..ff8f032 100644 --- a/vignettes/math.Rmd +++ b/vignettes/math.Rmd @@ -96,11 +96,19 @@ transitions. Currently the package uses a standard normal prior for the coefficients $\beta$. The prior for the weights of the spline basis functions is -set hierarchically so that transitions -of same type have a shared mean. Also the average log hazard rates -$b_0^{(h)}$ have a hierarchical prior so that transitions of the same type -have a shared mean. This prior mean is estimated from the data using a standard -Cox proportional hazards model fit. +set hierarchically so that transitions of the same type have a shared mean. +The realized spline weights are centered to sum to zero, separating their +time-varying component from the log baseline-hazard intercept. For transition +$h$, the intercept prior is centered at $\log(\widehat r_h)$, where +$\widehat r_h$ is calculated by pooling event counts and at-risk time within +transition type. A half-event is used for an exposed transition type with no +events, while an entirely unexposed type borrows the global pooled rate. +Consequently, $\widehat r_h$ is a natural-scale prior reference rate (a +geometric location or conditional median), not an arithmetic prior mean. +Indeed, under the current $\sigma_0 \sim \operatorname{HalfNormal}(0,3)$ +hyperprior, the marginal arithmetic mean of the baseline hazard is not finite; +prior-predictive checks and sensitivity to this deliberately broad scale are +therefore important. You can view the prior in the Stan code for example like this. @@ -512,5 +520,3 @@ pd$paths$prop_matrix() # References - - From 1c491fd7c9d4d0001697e0b9aaa33daa7976b002 Mon Sep 17 00:00:00 2001 From: jtimonen Date: Tue, 28 Jul 2026 12:32:47 +0300 Subject: [PATCH 4/7] revise --- NEWS.md | 28 +++--- R/MultistateSystem.R | 55 ++++-------- R/stan-data.R | 52 +++++++---- inst/stan/msm.stan | 20 ++--- man/MultistateSystem.Rd | 8 +- tests/testthat/test-correctness-regressions.R | 89 ++++++++++++------- vignettes/math.Rmd | 9 +- 7 files changed, 146 insertions(+), 115 deletions(-) diff --git a/NEWS.md b/NEWS.md index ab7eb4c..c209367 100644 --- a/NEWS.md +++ b/NEWS.md @@ -15,9 +15,11 @@ could leave their empirical prior locations undefined. Reference rates are now computed from event counts divided by exact at-risk time, pooled within transition type. An exposed type with no events uses a half-event - regularization, and an entirely unexposed type borrows the global pooled - rate. Documentation now calls this quantity a natural-scale reference rate, - not an arithmetic prior mean. + regularization for prior calibration. A transition with no positive at-risk + time is rejected explicitly because it contributes neither an event term nor + an integrated-hazard term to the likelihood. 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 @@ -32,14 +34,6 @@ state-by-group grid, joins event counts and denominators separately, and returns explicit zero counts and probabilities. -* **Make event-time thinning valid for flexible spline hazards.** The former - upper bound was the largest hazard on a 100-point grid, which could miss a - narrow spline peak and invalidate thinning. It also discarded transitions - whose bound was at most `1e-9`, introducing rare-event truncation. The new - bound uses non-negativity, partition of unity, and local support of the - B-spline basis to dominate the hazard over the complete requested interval; - all transitions with a positive bound remain eligible. - * **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. @@ -74,11 +68,13 @@ 95th percentiles and returned a 90% interval. It now uses `(1 - ci_alpha) / 2` and its complement directly. -* **Support transitions with zero events or zero at-risk intervals in Stan - data.** Zero counts previously produced zero-width ragged arrays while the - Stan declarations required positive dimensions and at least one risk row. - R now pads storage to width one while retaining a count of zero; Stan accepts - zero counts and skips the corresponding likelihood contribution. +* **Support zero-event transitions without treating unobserved transitions as + estimable.** A transition with no observed events but positive at-risk time + still contributes its integrated-hazard survival term. R now pads only the + zero-width event-index storage while retaining an event count of zero, and + Stan skips only the absent event term. By contrast, a transition with no + positive at-risk time has no hazard-likelihood information and now produces + an explicit error instead of a prior-only fit. * **Separate baseline level from spline shape.** Because the intercept-inclusive B-spline basis partitions unity, adding a constant to diff --git a/R/MultistateSystem.R b/R/MultistateSystem.R index 295111a..9b80bc0 100644 --- a/R/MultistateSystem.R +++ b/R/MultistateSystem.R @@ -97,18 +97,24 @@ MultistateSystem <- R6::R6Class("MultistateSystem", # Generate transition given state and transition intensity functions generate_transition = function(state, t_init, t_max, w, log_w0, log_m) { + tol <- 1.05 + rare_hazard_cutoff <- 1e-9 possible <- private$transmat$possible_transitions_from(state) if (length(possible) == 0) { # Absorbing state, no transitions possible return(list(t = t_max, new_state = 0, idx = 0)) } - UB <- vapply(possible, function(j) { + UB <- tol * vapply(possible, function(j) { self$max_inst_hazard(t_init, t_max, w[j, ], log_w0[j], log_m[j]) }, numeric(1)) if (any(!is.finite(UB))) { stop("non-finite transition-hazard upper bound") } - keep <- UB > 0 + # This is an intentional simulation approximation: transitions below + # this absolute rate threshold are treated as impossible to avoid doing + # proposal work for events that are negligible on the package's intended + # time scale. + keep <- UB > rare_hazard_cutoff possible <- possible[keep] UB <- UB[keep] if (length(possible) == 0) { @@ -363,7 +369,9 @@ MultistateSystem <- R6::R6Class("MultistateSystem", sum(diag(self$tm()$matrix)) > 0 }, - #' @description Max instant hazard on interval (t1, t2) + #' @description Approximate maximum instantaneous hazard on an interval, + #' evaluated on a 100-point grid. Path generation inflates this value by + #' five percent before thinning. #' #' @param t1 Start time point #' @param t2 End time point @@ -381,41 +389,16 @@ MultistateSystem <- R6::R6Class("MultistateSystem", return(0) } - # An intercept-inclusive B-spline basis is nonnegative and partitions - # unity inside its boundary knots. On a given interval, its weighted sum - # is therefore bounded above by the largest coefficient whose compact - # support intersects that interval. - knots <- self$get_knots() - order <- private$spline_k - internal_knots <- if (length(knots) > 2) { - knots[2:(length(knots) - 1)] - } else { - numeric(0) - } - augmented_knots <- c( - rep(knots[1], order), - internal_knots, - rep(knots[length(knots)], order) - ) - weight_idx <- seq_len(self$num_weights()) - support_left <- augmented_knots[weight_idx] - support_right <- augmented_knots[weight_idx + order] - active <- support_right > t1 & support_left < t2 - - # Include basis functions active exactly at either closed endpoint. This - # also handles a zero-length interval without broadening to global support. - endpoint_basis <- self$basisfun_matrix(unique(c(t1, t2))) - active <- active | colSums(endpoint_basis > 0) > 0 - if (!any(active)) { - stop("could not identify an active B-spline basis function") - } - # A negligible upward margin protects thinning from floating-point - # overshoot when the same hazard is subsequently evaluated numerically. - log_margin <- log1p(sqrt(.Machine$double.eps)) - exp(log_w0 + log_m + max(w[active]) + log_margin) + # Fast numerical envelope used by path generation. The caller inflates + # this grid maximum by five percent before thinning. + t_grid <- seq(t1, t2, length.out = 100) + log_haz <- log_m + self$log_baseline_hazard(t_grid, log_w0, w) + exp(max(log_haz)) }, - #' @description Generate paths + #' @description Generate paths. For computational efficiency, a transition + #' whose five-percent-inflated numerical hazard envelope is at most + #' \code{1e-9} is treated as impossible. #' #' @param w An array of shape \code{n_draws} x \code{n_trans} x #' \code{n_weights} diff --git a/R/stan-data.R b/R/stan-data.R index 7794fc0..d26e3d9 100644 --- a/R/stan-data.R +++ b/R/stan-data.R @@ -103,6 +103,13 @@ create_stan_data_timegrid <- function(model) { create_stan_data_transitions <- function(pd) { tm <- pd$transmat dat <- pd$as_transitions() + interval_length <- dat$time - dat$time_prev + if (any(interval_length < 0)) { + stop("transition data contain a negative at-risk interval") + } + if (any(dat$trans_idx > 0 & interval_length <= 0)) { + stop("an observed transition has no positive at-risk time") + } N_int <- nrow(dat) N_trans <- tm$num_trans() transition <- matrix(0, N_trans, N_int) @@ -114,6 +121,16 @@ create_stan_data_transitions <- function(pd) { } at_risk[tm$possible_transitions_from(dat$from[n]), n] <- 1 } + risk_time <- as.vector(at_risk %*% interval_length) + unexposed <- which(risk_time <= 0) + if (length(unexposed) > 0) { + labels <- tm$trans_df()$trans_char[unexposed] + stop( + "cannot fit transition(s) with no positive at-risk time: ", + paste(labels, collapse = ", "), + ". Remove those transitions from the model or provide informative data." + ) + } out <- list( at_risk = at_risk, transition = transition, @@ -306,7 +323,8 @@ create_stan_data_intervalidx <- function(t_start, t_end, t_grid, delta_grid) { # Estimate a constant hazard per transition type by pooling event counts and # at-risk time over all transitions of that type. If an exposed type has no # events, one Jeffreys-style half-event avoids an undefined empirical prior. -# Entirely unexposed types borrow the global pooled transition rate. +# A transition with no positive at-risk time is rejected because its hazard +# parameters have no likelihood contribution. average_haz_per_ttype <- function(pd) { checkmate::assert_class(pd, "PathData") dat <- pd$as_transitions(truncate = TRUE) @@ -323,15 +341,13 @@ average_haz_per_ttype <- function(pd) { } } - total_exposure <- sum(exposure) - total_events <- sum(n_event) - if (total_exposure <= 0) { - stop("cannot calibrate baseline-hazard priors: no transition has positive at-risk time") - } - global_rate <- if (total_events > 0) { - total_events / total_exposure - } else { - 0.5 / total_exposure + unexposed <- which(exposure <= 0) + if (length(unexposed) > 0) { + stop( + "cannot calibrate or fit transition(s) with no positive at-risk time: ", + paste(trans$trans_char[unexposed], collapse = ", "), + ". Remove those transitions from the model or provide informative data." + ) } rate_df <- data.frame( @@ -348,11 +364,8 @@ average_haz_per_ttype <- function(pd) { .groups = "drop" ) |> dplyr::mutate( - rate = ifelse( - .data$exposure > 0, - ifelse(.data$n_event > 0, .data$n_event, 0.5) / .data$exposure, - global_rate - ), + rate = ifelse(.data$n_event > 0, .data$n_event, 0.5) / + .data$exposure, log_h0_avg = log(.data$rate) ) |> dplyr::select("trans_type", "log_h0_avg") @@ -365,7 +378,14 @@ average_haz_per_ttype <- function(pd) { # transition which_format_for_stan <- function(x, name) { H <- nrow(x) - M <- max(1L, max(rowSums(x))) + row_counts <- rowSums(x) + if (identical(name, "risk") && any(row_counts == 0)) { + stop( + "cannot create Stan data: no at-risk intervals for transition(s) ", + paste(which(row_counts == 0), collapse = ", ") + ) + } + M <- max(1L, max(row_counts)) N_sum <- rep(0, H) a <- matrix(0, H, M) for (h in seq_len(H)) { diff --git a/inst/stan/msm.stan b/inst/stan/msm.stan index eb11a1b..dec3406 100644 --- a/inst/stan/msm.stan +++ b/inst/stan/msm.stan @@ -336,7 +336,7 @@ data { int D_trans; // max num of occurrences for a transition int D_risk; // max num of at-risk intervals for a transition array[N_trans] int sum_trans; // total number of occurred - array[N_trans] int sum_risk; // total number at risk + array[N_trans] int sum_risk; // total number at risk array[N_trans, D_trans] int which_trans; array[N_trans, D_risk] int which_risk; @@ -521,16 +521,14 @@ model { weights[1,h], log_w0[1,h] ); } - if(sum_risk[h] > 0){ - array[sum_risk[h]] int idx_atr = which_risk[h, 1:sum_risk[h]]; - vector[N_grid+1] h0_grid = rep_vector(0.0, N_grid+1); - vector[N_grid+1] h0_int; - h0_grid[2:(N_grid+1)] = exp(log_basehaz(SBF_grid, weights[1,h], log_w0[1,h])); - h0_int = cumulative_sum(h0_grid) * delta_grid; - target += - exp(log_C_haz[idx_sub[idx_atr], ttype[h]]) .* - (h0_int[t_end_idx[idx_atr]] - h0_int[t_start_idx_m1[idx_atr]]) .* - correction_multiplier[idx_atr]; - } + array[sum_risk[h]] int idx_atr = which_risk[h, 1:sum_risk[h]]; + vector[N_grid+1] h0_grid = rep_vector(0.0, N_grid+1); + vector[N_grid+1] h0_int; + h0_grid[2:(N_grid+1)] = exp(log_basehaz(SBF_grid, weights[1,h], log_w0[1,h])); + h0_int = cumulative_sum(h0_grid) * delta_grid; + target += - exp(log_C_haz[idx_sub[idx_atr], ttype[h]]) .* + (h0_int[t_end_idx[idx_atr]] - h0_int[t_start_idx_m1[idx_atr]]) .* + correction_multiplier[idx_atr]; } } diff --git a/man/MultistateSystem.Rd b/man/MultistateSystem.Rd index 47a4327..6d9931e 100644 --- a/man/MultistateSystem.Rd +++ b/man/MultistateSystem.Rd @@ -296,7 +296,9 @@ Boolean value \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-MultistateSystem-max_inst_hazard}{}}} \subsection{Method \code{max_inst_hazard()}}{ -Max instant hazard on interval (t1, t2) +Approximate maximum instantaneous hazard on an interval, +evaluated on a 100-point grid. Path generation inflates this value by +five percent before thinning. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{MultistateSystem$max_inst_hazard(t1, t2, w, log_w0, log_m)}\if{html}{\out{
}} } @@ -322,7 +324,9 @@ baseline hazard} \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-MultistateSystem-simulate}{}}} \subsection{Method \code{simulate()}}{ -Generate paths +Generate paths. For computational efficiency, a transition +whose five-percent-inflated numerical hazard envelope is at most +\code{1e-9} is treated as impossible. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{MultistateSystem$simulate( w, diff --git a/tests/testthat/test-correctness-regressions.R b/tests/testthat/test-correctness-regressions.R index 938956d..ed41a3b 100644 --- a/tests/testthat/test-correctness-regressions.R +++ b/tests/testthat/test-correctness-regressions.R @@ -26,7 +26,7 @@ test_that("baseline-hazard prior uses log scale in Stan data", { expect_error(mod$set_prior_mean_h0(0), "strictly positive") }) -test_that("all-censored data get finite prior and padded Stan indices", { +test_that("all-censored data retain survival information and padded event indices", { pd <- make_sparse_survival_paths(event = FALSE) mod <- create_msm(pd$transmat, t_max = 2, n_grid = 10) mod$set_prior_mean_h0_data(pd) @@ -37,9 +37,10 @@ test_that("all-censored data get finite prior and padded Stan indices", { expect_equal(sd$sum_trans, 0) expect_equal(sd$sum_risk, 2) - empty <- which_format_for_stan(matrix(0, nrow = 2, ncol = 3), "risk") - expect_equal(empty$D_risk, 1) - expect_equal(empty$sum_risk, c(0, 0)) + expect_error( + which_format_for_stan(matrix(0, nrow = 2, ncol = 3), "risk"), + "no at-risk intervals" + ) }) test_that("empirical hazard calibration pools within transition type", { @@ -64,7 +65,7 @@ test_that("empirical hazard calibration pools within transition type", { expect_equal(exp(rates$log_h0_avg), c(1 / 3, 0.5 / 5, 0.5 / 5)) }) -test_that("an unexposed transition type borrows the global pooled rate", { +test_that("an unexposed transition is rejected", { tm <- TransitionMatrix$new( matrix(c( 0, 1, 0, 0, @@ -88,17 +89,47 @@ test_that("an unexposed transition type borrows the global pooled rate", { rep_idx = 1L ) pd <- PathData$new(subject_df, path_df, link_df, tm) - rates <- average_haz_per_ttype(pd) |> - dplyr::arrange(.data$trans_idx) - - expect_equal(exp(rates$log_h0_avg), c(1 / 3, 1 / 3)) + expect_error( + average_haz_per_ttype(pd), + "no positive at-risk time: C -> D" + ) mod <- create_msm(tm, t_max = 2, n_grid = 10) - mod$set_prior_mean_h0_data(pd) - sd <- create_stan_data(mod, pd) - expect_equal(sd$sum_trans, c(1, 0)) - expect_equal(sd$sum_risk, c(2, 0)) - expect_equal(sd$D_trans, 1) + mod$set_prior_mean_h0(c(1 / 3, 1 / 3)) + expect_error( + create_stan_data(mod, pd), + "no positive at-risk time: C -> D" + ) +}) + +test_that("zero-duration rows do not count as hazard information", { + tm <- transmat_survival(c("Healthy", "Dead")) + subject_df <- data.frame(subject_id = "a") + link_df <- data.frame( + path_id = 1L, subject_id = "a", draw_idx = 1L, rep_idx = 1L + ) + make_pd <- function(event) { + PathData$new( + subject_df, + data.frame( + path_id = c(1L, 1L), + state = c(1L, if (event) 2L else 1L), + time = c(0, 0), + trans_idx = c(0L, if (event) 1L else 0L) + ), + link_df, + tm + ) + } + + expect_error( + create_stan_data_transitions(make_pd(TRUE)), + "observed transition has no positive at-risk time" + ) + expect_error( + create_stan_data_transitions(make_pd(FALSE)), + "no positive at-risk time" + ) }) test_that("state-visit probabilities retain zero-event subjects", { @@ -172,31 +203,29 @@ test_that("identity transition probabilities still validate hazard parameters", ) }) -test_that("B-spline thinning bound dominates the evaluated hazard", { +test_that("thinning uses the 100-point hazard envelope", { sys <- MultistateSystem$new(transmat_survival()) sys$set_knots(c(0, 0.4999, 0.5001, 1)) w <- c(-3, 3, -3, 3, -3) bound <- sys$max_inst_hazard(0, 1, w, log_w0 = -2, log_m = 0.4) - dense <- exp(sys$log_inst_hazard(seq(0, 1, length.out = 10001), w, -2, 0.4)) + grid <- seq(0, 1, length.out = 100) + evaluated <- exp(sys$log_inst_hazard(grid, w, -2, 0.4)) - expect_gte(bound, max(dense)) - expect_gt(bound, exp(-2 + 0.4 + max(w))) + expect_equal(bound, max(evaluated)) }) -test_that("B-spline thinning bound is local to the requested interval", { +test_that("path generation skips hazards below the deliberate cutoff", { sys <- MultistateSystem$new(transmat_survival()) - sys$set_knots(c(0, 0.4999, 0.5001, 1)) - w <- c(-10, 10, -10, -10, -10) - times <- seq(0.8, 1, length.out = 10001) - - local_bound <- sys$max_inst_hazard( - 0.8, 1, w, log_w0 = 0, log_m = 0 - ) - global_bound <- exp(max(w)) - dense <- exp(sys$log_inst_hazard(times, w, 0, 0)) + sys$set_knots(c(0, 1e9)) + w <- array(0, dim = c(1, 1, sys$num_weights())) + path <- suppressMessages(sys$simulate( + w = w, + log_w0 = matrix(log(0.5e-9), nrow = 1), + log_m = matrix(0, nrow = 1), + t_max = 1e9 + )) - expect_gte(local_bound, max(dense)) - expect_lt(local_bound, global_bound / 1e6) + expect_equal(path$trans_idx, c(0, 0)) }) test_that("zero baseline hazard has zero bound and generates no event", { diff --git a/vignettes/math.Rmd b/vignettes/math.Rmd index ff8f032..788c0c1 100644 --- a/vignettes/math.Rmd +++ b/vignettes/math.Rmd @@ -102,9 +102,11 @@ time-varying component from the log baseline-hazard intercept. For transition $h$, the intercept prior is centered at $\log(\widehat r_h)$, where $\widehat r_h$ is calculated by pooling event counts and at-risk time within transition type. A half-event is used for an exposed transition type with no -events, while an entirely unexposed type borrows the global pooled rate. -Consequently, $\widehat r_h$ is a natural-scale prior reference rate (a -geometric location or conditional median), not an arithmetic prior mean. +events. If any modeled transition has no positive at-risk time, fitting stops +because that transition contributes no hazard-likelihood information; it must +be removed from the model or supported by additional data. Consequently, +$\widehat r_h$ is a natural-scale prior reference rate (a geometric location +or conditional median), not an arithmetic prior mean. Indeed, under the current $\sigma_0 \sim \operatorname{HalfNormal}(0,3)$ hyperprior, the marginal arithmetic mean of the baseline hazard is not finite; prior-predictive checks and sensitivity to this deliberately broad scale are @@ -519,4 +521,3 @@ pd$paths$prop_matrix() # References - From 12f81d704bfc78af87e0fa000a4513e8f7cd0469 Mon Sep 17 00:00:00 2001 From: jtimonen Date: Tue, 28 Jul 2026 12:52:02 +0300 Subject: [PATCH 5/7] revert some things --- NEWS.md | 28 +++------ R/stan-data.R | 36 ++++------- R/utils.R | 6 ++ inst/stan/msm.stan | 14 ++--- tests/testthat/test-correctness-regressions.R | 63 ++++++++++--------- vignettes/math.Rmd | 12 ++-- 6 files changed, 73 insertions(+), 86 deletions(-) diff --git a/NEWS.md b/NEWS.md index c209367..98ed87b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,16 +10,16 @@ programs supplied to `fit_stan(filepath = ...)` must adopt the new name. Fits produced by affected versions must be rerun. -* **Make sparse-transition prior calibration finite and explicit.** The - previous cumulative-hazard summary omitted transitions with no events and - could leave their empirical prior locations undefined. Reference rates are - now computed from event counts divided by exact at-risk time, pooled within - transition type. An exposed type with no events uses a half-event - regularization for prior calibration. A transition with no positive at-risk - time is rejected explicitly because it contributes neither an event term nor - an integrated-hazard term to the likelihood. Documentation now calls the - calibrated quantity a natural-scale reference rate, not an arithmetic prior - mean. +* **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 @@ -68,14 +68,6 @@ 95th percentiles and returned a 90% interval. It now uses `(1 - ci_alpha) / 2` and its complement directly. -* **Support zero-event transitions without treating unobserved transitions as - estimable.** A transition with no observed events but positive at-risk time - still contributes its integrated-hazard survival term. R now pads only the - zero-width event-index storage while retaining an event count of zero, and - Stan skips only the absent event term. By contrast, a transition with no - positive at-risk time has no hazard-likelihood information and now produces - an explicit error instead of a prior-only fit. - * **Separate baseline level from spline shape.** Because the intercept-inclusive B-spline basis partitions unity, adding a constant to every spline coefficient and subtracting it from `log_w0` left the likelihood diff --git a/R/stan-data.R b/R/stan-data.R index d26e3d9..d84cca4 100644 --- a/R/stan-data.R +++ b/R/stan-data.R @@ -121,16 +121,6 @@ create_stan_data_transitions <- function(pd) { } at_risk[tm$possible_transitions_from(dat$from[n]), n] <- 1 } - risk_time <- as.vector(at_risk %*% interval_length) - unexposed <- which(risk_time <= 0) - if (length(unexposed) > 0) { - labels <- tm$trans_df()$trans_char[unexposed] - stop( - "cannot fit transition(s) with no positive at-risk time: ", - paste(labels, collapse = ", "), - ". Remove those transitions from the model or provide informative data." - ) - } out <- list( at_risk = at_risk, transition = transition, @@ -321,10 +311,8 @@ create_stan_data_intervalidx <- function(t_start, t_end, t_grid, delta_grid) { } # Estimate a constant hazard per transition type by pooling event counts and -# at-risk time over all transitions of that type. If an exposed type has no -# events, one Jeffreys-style half-event avoids an undefined empirical prior. -# A transition with no positive at-risk time is rejected because its hazard -# parameters have no likelihood contribution. +# at-risk time over all transitions of that type. Every modeled transition must +# have at least one observed event and positive at-risk time. average_haz_per_ttype <- function(pd) { checkmate::assert_class(pd, "PathData") dat <- pd$as_transitions(truncate = TRUE) @@ -341,6 +329,14 @@ average_haz_per_ttype <- function(pd) { } } + unobserved <- which(n_event == 0) + if (length(unobserved) > 0) { + stop( + "cannot calibrate or fit transition(s) with no observed events: ", + paste(trans$trans_char[unobserved], collapse = ", "), + ". Remove those transitions from the model or provide informative data." + ) + } unexposed <- which(exposure <= 0) if (length(unexposed) > 0) { stop( @@ -364,8 +360,7 @@ average_haz_per_ttype <- function(pd) { .groups = "drop" ) |> dplyr::mutate( - rate = ifelse(.data$n_event > 0, .data$n_event, 0.5) / - .data$exposure, + rate = .data$n_event / .data$exposure, log_h0_avg = log(.data$rate) ) |> dplyr::select("trans_type", "log_h0_avg") @@ -378,14 +373,7 @@ average_haz_per_ttype <- function(pd) { # transition which_format_for_stan <- function(x, name) { H <- nrow(x) - row_counts <- rowSums(x) - if (identical(name, "risk") && any(row_counts == 0)) { - stop( - "cannot create Stan data: no at-risk intervals for transition(s) ", - paste(which(row_counts == 0), collapse = ", ") - ) - } - M <- max(1L, max(row_counts)) + M <- max(rowSums(x)) N_sum <- rep(0, H) a <- matrix(0, H, M) for (h in seq_len(H)) { diff --git a/R/utils.R b/R/utils.R index 469186e..d3ff162 100644 --- a/R/utils.R +++ b/R/utils.R @@ -140,6 +140,12 @@ truncate_after_terminal_events <- function(df, term_state_inds) { # make sense for given data prefit_checks <- function(model, data) { pd <- data$paths + if (!model$pk_only && !model$prior_only) { + # This validates that every transition retained in the fitted hazard model + # has at least one event and positive at-risk time. Prediction data are not + # subject to this requirement. + average_haz_per_ttype(pd) + } lens <- pd$as_transitions() |> dplyr::mutate(time_len = .data$time - .data$time_prev) |> dplyr::pull(.data$time_len) diff --git a/inst/stan/msm.stan b/inst/stan/msm.stan index dec3406..533e440 100644 --- a/inst/stan/msm.stan +++ b/inst/stan/msm.stan @@ -335,7 +335,7 @@ data { int N_int; // total number of intervals int D_trans; // max num of occurrences for a transition int D_risk; // max num of at-risk intervals for a transition - array[N_trans] int sum_trans; // total number of occurred + array[N_trans] int sum_trans; // total number of occurred array[N_trans] int sum_risk; // total number at risk array[N_trans, D_trans] int which_trans; array[N_trans, D_risk] int which_risk; @@ -514,13 +514,11 @@ model { xpsr_loc, xpsr_scale ); for(h in 1:N_trans){ - if(sum_trans[h] > 0){ - array[sum_trans[h]] int idx_occ = which_trans[h, 1:sum_trans[h]]; - target += log_hazard( - log_C_haz[idx_sub[idx_occ], ttype[h]], SBF_end[idx_occ,:], - weights[1,h], log_w0[1,h] - ); - } + array[sum_trans[h]] int idx_occ = which_trans[h, 1:sum_trans[h]]; + target += log_hazard( + log_C_haz[idx_sub[idx_occ], ttype[h]], SBF_end[idx_occ,:], + weights[1,h], log_w0[1,h] + ); array[sum_risk[h]] int idx_atr = which_risk[h, 1:sum_risk[h]]; vector[N_grid+1] h0_grid = rep_vector(0.0, N_grid+1); vector[N_grid+1] h0_int; diff --git a/tests/testthat/test-correctness-regressions.R b/tests/testthat/test-correctness-regressions.R index ed41a3b..d597eeb 100644 --- a/tests/testthat/test-correctness-regressions.R +++ b/tests/testthat/test-correctness-regressions.R @@ -26,35 +26,37 @@ test_that("baseline-hazard prior uses log scale in Stan data", { expect_error(mod$set_prior_mean_h0(0), "strictly positive") }) -test_that("all-censored data retain survival information and padded event indices", { +test_that("a modeled transition with no events is rejected", { pd <- make_sparse_survival_paths(event = FALSE) mod <- create_msm(pd$transmat, t_max = 2, n_grid = 10) - mod$set_prior_mean_h0_data(pd) - sd <- create_stan_data(mod, pd) - - expect_equal(mod$get_prior_mean_h0(), 0.5 / 3) - expect_equal(sd$D_trans, 1) - expect_equal(sd$sum_trans, 0) - expect_equal(sd$sum_risk, 2) expect_error( - which_format_for_stan(matrix(0, nrow = 2, ncol = 3), "risk"), - "no at-risk intervals" + prefit_checks(mod, JointData$new(pd, NULL)), + "no observed events: Healthy -> Dead" + ) + expect_error( + mod$set_prior_mean_h0_data(pd), + "no observed events: Healthy -> Dead" ) + mod$set_prior_mean_h0(1 / 3) + sd <- create_stan_data(mod, pd) + expect_equal(sd$D_trans, 0) + expect_equal(sd$sum_trans, 0) + expect_equal(sd$sum_risk, 2) }) test_that("empirical hazard calibration pools within transition type", { tm <- transmat_illnessdeath(c("Healthy", "Ill", "Dead")) - subject_df <- data.frame(subject_id = c("a", "b")) + subject_df <- data.frame(subject_id = c("a", "b", "c")) path_df <- data.frame( - path_id = c(1L, 1L, 1L, 2L, 2L), - state = c(1L, 2L, 2L, 1L, 1L), - time = c(0, 1, 3, 0, 2), - trans_idx = c(0L, 1L, 0L, 0L, 0L) + path_id = c(1L, 1L, 1L, 2L, 2L, 3L, 3L), + state = c(1L, 2L, 3L, 1L, 3L, 1L, 1L), + time = c(0, 1, 3, 0, 2, 0, 4), + trans_idx = c(0L, 1L, 3L, 0L, 2L, 0L, 0L) ) link_df <- data.frame( - path_id = 1:2, - subject_id = c("a", "b"), + path_id = 1:3, + subject_id = c("a", "b", "c"), draw_idx = 1L, rep_idx = 1L ) @@ -62,10 +64,10 @@ test_that("empirical hazard calibration pools within transition type", { rates <- average_haz_per_ttype(pd) |> dplyr::arrange(.data$trans_idx) - expect_equal(exp(rates$log_h0_avg), c(1 / 3, 0.5 / 5, 0.5 / 5)) + expect_equal(exp(rates$log_h0_avg), c(1 / 7, 2 / 9, 2 / 9)) }) -test_that("an unexposed transition is rejected", { +test_that("fit validation rejects unsupported transitions but prediction data do not", { tm <- TransitionMatrix$new( matrix(c( 0, 1, 0, 0, @@ -89,17 +91,21 @@ test_that("an unexposed transition is rejected", { rep_idx = 1L ) pd <- PathData$new(subject_df, path_df, link_df, tm) + mod <- create_msm(tm, t_max = 2, n_grid = 10) + + expect_error( + prefit_checks(mod, JointData$new(pd, NULL)), + "no observed events: C -> D" + ) expect_error( average_haz_per_ttype(pd), - "no positive at-risk time: C -> D" + "no observed events: C -> D" ) - mod <- create_msm(tm, t_max = 2, n_grid = 10) mod$set_prior_mean_h0(c(1 / 3, 1 / 3)) - expect_error( - create_stan_data(mod, pd), - "no positive at-risk time: C -> D" - ) + sd <- create_stan_data(mod, pd) + expect_equal(sd$sum_trans, c(1, 0)) + expect_equal(sd$sum_risk, c(2, 0)) }) test_that("zero-duration rows do not count as hazard information", { @@ -126,10 +132,9 @@ test_that("zero-duration rows do not count as hazard information", { create_stan_data_transitions(make_pd(TRUE)), "observed transition has no positive at-risk time" ) - expect_error( - create_stan_data_transitions(make_pd(FALSE)), - "no positive at-risk time" - ) + sd <- create_stan_data_transitions(make_pd(FALSE)) + expect_equal(sd$D_trans, 0) + expect_equal(sd$sum_trans, 0) }) test_that("state-visit probabilities retain zero-event subjects", { diff --git a/vignettes/math.Rmd b/vignettes/math.Rmd index 788c0c1..0dc2e6c 100644 --- a/vignettes/math.Rmd +++ b/vignettes/math.Rmd @@ -101,12 +101,11 @@ The realized spline weights are centered to sum to zero, separating their time-varying component from the log baseline-hazard intercept. For transition $h$, the intercept prior is centered at $\log(\widehat r_h)$, where $\widehat r_h$ is calculated by pooling event counts and at-risk time within -transition type. A half-event is used for an exposed transition type with no -events. If any modeled transition has no positive at-risk time, fitting stops -because that transition contributes no hazard-likelihood information; it must -be removed from the model or supported by additional data. Consequently, -$\widehat r_h$ is a natural-scale prior reference rate (a geometric location -or conditional median), not an arithmetic prior mean. +transition type. Every transition retained in the fitted transition matrix +must have at least one observed occurrence and positive at-risk time; fitting +stops otherwise. Consequently, $\widehat r_h$ is a natural-scale prior +reference rate (a geometric location or conditional median), not an arithmetic +prior mean. Indeed, under the current $\sigma_0 \sim \operatorname{HalfNormal}(0,3)$ hyperprior, the marginal arithmetic mean of the baseline hazard is not finite; prior-predictive checks and sensitivity to this deliberately broad scale are @@ -520,4 +519,3 @@ pd$paths$prop_matrix() # References - From b000567d53c95f8cda4141f73d2d5d46189b02ea Mon Sep 17 00:00:00 2001 From: jtimonen Date: Tue, 28 Jul 2026 13:09:03 +0300 Subject: [PATCH 6/7] fix option --- NEWS.md | 3 ++- R/utils.R | 2 +- tests/testthat/test-correctness-regressions.R | 5 +++++ 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 98ed87b..5cb0877 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,7 +7,8 @@ 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 supplied to `fit_stan(filepath = ...)` must adopt the new name. + 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 diff --git a/R/utils.R b/R/utils.R index d3ff162..7942caf 100644 --- a/R/utils.R +++ b/R/utils.R @@ -140,7 +140,7 @@ truncate_after_terminal_events <- function(df, term_state_inds) { # make sense for given data prefit_checks <- function(model, data) { pd <- data$paths - if (!model$pk_only && !model$prior_only) { + if (!model$pk_only) { # This validates that every transition retained in the fitted hazard model # has at least one event and positive at-risk time. Prediction data are not # subject to this requirement. diff --git a/tests/testthat/test-correctness-regressions.R b/tests/testthat/test-correctness-regressions.R index d597eeb..b17aed8 100644 --- a/tests/testthat/test-correctness-regressions.R +++ b/tests/testthat/test-correctness-regressions.R @@ -38,6 +38,11 @@ test_that("a modeled transition with no events is rejected", { mod$set_prior_mean_h0_data(pd), "no observed events: Healthy -> Dead" ) + mod$prior_only <- TRUE + expect_error( + prefit_checks(mod, JointData$new(pd, NULL)), + "no observed events: Healthy -> Dead" + ) mod$set_prior_mean_h0(1 / 3) sd <- create_stan_data(mod, pd) expect_equal(sd$D_trans, 0) From 01aa63f892329416f2f6e9ca24dec1cd89972aac Mon Sep 17 00:00:00 2001 From: jtimonen Date: Tue, 28 Jul 2026 13:29:24 +0300 Subject: [PATCH 7/7] more revert --- NEWS.md | 8 -------- inst/stan/msm.stan | 5 +---- vignettes/math.Rmd | 9 ++++++--- 3 files changed, 7 insertions(+), 15 deletions(-) diff --git a/NEWS.md b/NEWS.md index 5cb0877..634cee7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -69,14 +69,6 @@ 95th percentiles and returned a 90% interval. It now uses `(1 - ci_alpha) / 2` and its complement directly. -* **Separate baseline level from spline shape.** Because the - intercept-inclusive B-spline basis partitions unity, adding a constant to - every spline coefficient and subtracting it from `log_w0` left the likelihood - unchanged. Evaluated spline weights are now centered to sum to zero, so the - baseline intercept alone carries the constant level. The full-dimensional - raw hierarchy still contains prior-only common-shift directions; this change - removes the likelihood confounding, not every latent nuisance direction. - * **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 diff --git a/inst/stan/msm.stan b/inst/stan/msm.stan index 533e440..2cedd2b 100644 --- a/inst/stan/msm.stan +++ b/inst/stan/msm.stan @@ -434,11 +434,8 @@ transformed parameters { // Baseline hazard if(do_haz == 1) { for(j in 1:N_trans){ - vector[N_sbf] weights_raw = mu_weights[1][:,ttype[j]] + + weights[1,j] = mu_weights[1][:,ttype[j]] + sig_weights[1][:,ttype[j]] .* z_weights[1][:,j]; - // The intercept-inclusive B-spline basis partitions unity. Centering - // removes its otherwise exact shift redundancy with log_w0. - weights[1,j] = weights_raw - mean(weights_raw); log_w0[1,j] = mu_log_w0[j] + sig_w0[1] * z_w0[1][j]; } } diff --git a/vignettes/math.Rmd b/vignettes/math.Rmd index 0dc2e6c..c5fe4cb 100644 --- a/vignettes/math.Rmd +++ b/vignettes/math.Rmd @@ -97,9 +97,12 @@ transitions. Currently the package uses a standard normal prior for the coefficients $\beta$. The prior for the weights of the spline basis functions is set hierarchically so that transitions of the same type have a shared mean. -The realized spline weights are centered to sum to zero, separating their -time-varying component from the log baseline-hazard intercept. For transition -$h$, the intercept prior is centered at $\log(\widehat r_h)$, where +Because the B-spline basis partitions unity, the likelihood alone cannot +distinguish a common shift of the spline weights from an opposite shift of the +log baseline-hazard intercept. The proper priors regularize this decomposition, +but inference about baseline hazard should be made from the complete hazard +curve rather than either component in isolation. For transition $h$, the +intercept prior is centered at $\log(\widehat r_h)$, where $\widehat r_h$ is calculated by pooling event counts and at-risk time within transition type. Every transition retained in the fitted transition matrix must have at least one observed occurrence and positive at-risk time; fitting