Exposure-hazard multistate modeling with bmstate
Juho Timonen
20th Nov 2025
exposure-hazard.RmdThis vignette is work in progress.
Data simulation
Setup
library(bmstate)
#> Attached bmstate 0.4.0. Type ?bmstate to get started.
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(tidyr)
library(tibble)
library(ggplot2)True data-generating model
# True beta
create_true_covariate_effects <- function(mod) {
NTT <- mod$system$num_states() - 1
C <- length(mod$covs())
bh_true <- matrix(0, NTT, C)
bh_true[1, 2] <- 1
bh_true[2, 2] <- -1
bh_true[3, 1] <- 0.3
sn <- event_state_names(mod)
rownames(bh_true) <- paste0("Effect on ", sn)
colnames(bh_true) <- mod$covs()
df <- data.frame(bh_true) |>
rownames_to_column("event") |>
pivot_longer(cols = -event, names_to = "covariate", values_to = "beta")
beta_pk <- list(CL = c(0.3, -0.3), V2 = c(0.3))
list(df = df, matrix = bh_true, pk = beta_pk)
}
# True baseline hazard parameters
create_true_baseline_hazard <- function(mod) {
tm <- mod$system$tm()
# Spline weights
w_true <- matrix(0, 7, 9)
for (j in 1:7) {
ww <- rep(0, 9)
if (tm$trans_df()$trans_type[j] == 2) {
ww <- rep(-1, 9)
ww[3:5] <- 1
}
if (tm$trans_df()$trans_type[j] == 3) {
ww <- rep(-1, 9)
ww[4:8] <- 1
}
w_true[j, ] <- ww
}
# Intercept
w0_true <- 0.5 * 1e-3
w0_true_vec <- rep(w0_true, 7)
w0_true_vec[3] <- 0.1 * w0_true
w0_true_vec[5] <- 5 * w0_true
w0_true_vec[7] <- 20 * w0_true
# Return
list(w0 = w0_true_vec, w = w_true)
}
# True data-generating model
create_true_model <- function() {
# Create models
sn <- c("Healthy", "Bleed", "Stroke", "Dead")
tm <- transmat_diamond(state_names = sn)
t3yr <- 3 * 365.25
haz_covs <- c("age")
pk_covs <- list(
CL = c("CrCL", "age"), V2 = "weight"
)
create_msm(
tm,
hazard_covs = haz_covs, pk_covs = pk_covs, num_knots = 8, t_max = t3yr
)
}
# Create oracle fit draws
create_oracle_fit <- function(similar_fit, beta_true, h0_true) {
checkmate::assert_true(similar_fit$is_point_estimate())
weights <- similar_fit$get_draws("weights")
weights[1, , ] <- h0_true$w
log_w0 <- similar_fit$get_draws("log_w0")
log_w0[1, ] <- log(h0_true$w0)
beta_oth <- similar_fit$get_draws("beta_oth")
beta_oth[1, , ] <- t(beta_true$matrix[, 1])
beta_oth <- posterior::rvar(beta_oth)
beta_xpsr <- similar_fit$get_draws("beta_xpsr")
beta_xpsr[1, , ] <- t(beta_true$matrix[, 2])
lp <- similar_fit$get_draws("lp__")
lp[] <- NA
log_z_pk <- similar_fit$get_draws("log_z_pk")
log_z_pk[, , ] <- 0
log_sig_pk <- similar_fit$get_draws("log_sig_pk")
log_sig_pk[, ] <- 0
log_mu_pk <- similar_fit$get_draws("log_mu_pk")
log_mu_pk[, ] <- 0
sigma_pk <- similar_fit$get_draws("sigma_pk")
sigma_pk[] <- 0.3
beta_CL <- similar_fit$get_draws("beta_CL")
beta_CL[1, ] <- beta_true$pk$CL
beta_V2 <- similar_fit$get_draws("beta_V2")
beta_V2[1, ] <- beta_true$pk$V2
# Create fit
draws <- list(
beta_oth = beta_oth,
beta_xpsr = beta_xpsr,
beta_CL = beta_CL,
beta_V2 = beta_V2,
weights = weights,
log_z_pk = log_z_pk,
log_sig_pk = log_sig_pk,
log_mu_pk = log_mu_pk,
sigma_pk = sigma_pk,
log_w0 = log_w0,
lp__ = lp
)
sd <- similar_fit$get_data()
mod <- similar_fit$model
MultistateModelFit$new(similar_fit$data, sd, mod, draws,
info = "Oracle fit"
)
}Data simulation
mod_true <- create_true_model()
beta_true <- create_true_covariate_effects(mod_true)
h0_true <- create_true_baseline_hazard(mod_true)
simdat <- mod_true$simulate_data(
params$N_subject,
beta_haz = beta_true$matrix,
beta_pk = beta_true$pk,
w0 = h0_true$w0,
w = h0_true$w
)
#> Recompiling Stan model
#> Using stan file at /home/runner/work/_temp/Library/bmstate/stan/msm.stan
#> Generating 600 paths
covs_dh <- unique(c(mod_true$data_covs(), "dose_amt"))
simdat_dh <- simdat$paths$subset_covariates(covs_dh, renamed_old = "dose", renamed_new = "dose_amt")
simdat_death <- as_single_event(simdat_dh, "Dead", null_state = "Healthy")
sa <- simdat$paths$subject_df$xpsr
mod_true$set_xpsr_normalizers(loc = mean(sa), scale = stats::sd(sa))
#> setting xpsr normalizers to loc = 5.34993, scale = 0.69161
simdat <- mod_true$simulate_data(
params$N_subject,
beta_haz = beta_true$matrix,
beta_pk = beta_true$pk,
w0 = h0_true$w0,
w = h0_true$w
)
#> Generating 600 paths
print(simdat)
#> A JointData object:
#> PathData object with 600 paths
#> * States = {Healthy, Bleed, Stroke, Dead}
#> * Covariates = {age, CrCL, weight, dose, t_pre, t_post, conc_pre, conc_post, xpsr, ka, CL, V2, pk_lloq}
#>
#> A DosingData object with 600 subjectsModeling
Defining models
NK <- 4
tm <- mod_true$system$tm()
pk_covs <- list(
ka = mod_true$pk_model$ka_covs(),
CL = mod_true$pk_model$CL_covs(),
V2 = mod_true$pk_model$V2_covs()
)
# Exposure-hazard multistate model
mod_ms_eh <- create_msm(
tm,
hazard_covs = setdiff(mod_true$covs(), "xpsr"),
pk_covs = pk_covs,
num_knots = NK, t_max = mod_true$get_tmax()
)
# Dose-hazard multistate model
mod_ms_dh <- create_msm(
tm,
hazard_covs = covs_dh, num_knots = NK, t_max = mod_true$get_tmax()
)
# Survival model
tm0 <- transmat_survival(tm$states[c(1, 4)])
mod_death <- create_msm(
tm0,
hazard_covs = covs_dh, num_knots = NK, t_max = mod_true$get_tmax()
)
# Inference model knots
t3yr <- mod_true$get_tmax()
tt1 <- simdat$paths$transition_times()
tt2 <- simdat_dh$transition_times()
tt3 <- simdat_death$transition_times()
mod_ms_eh$set_knots(t3yr, tt1, NK)
mod_ms_dh$set_knots(t3yr, tt2, NK)
mod_death$set_knots(t3yr, tt3, NK)
# Oracle fit
fit_prelim <- fit_stan(mod_true, simdat, method = "optimize", init = 0, iter = 100)
#> Shortest time interval (0.0359547288841782) is smaller than delta_grid (1.09575). Consider increasing n_grid or decreasing t_max of the model.
#> Using stan file at /home/runner/work/_temp/Library/bmstate/stan/msm.stan
#> setting xpsr normalizers to loc = 5.42315, scale = 0.56671
#> setting max conc = 9565.56038
#> Initial log joint probability = -11306.4
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 99 -7684.29 0.0366536 1007.03 1 1 104
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 100 -7670.48 0.111742 2081.25 1 1 105
#> Optimization terminated normally:
#> Maximum number of iterations hit, may not be at an optima
#> Finished in 0.2 seconds.
#> If these don't roughly match, consider refitting after setting xpsr normalizer loc and scale closer to estimated mean and sd. Otherwise interpret baseline hazards and xpsr effect size accodingly.
#> - xpsr normalization loc = 5.42315, mean estimated xpsr = 5.4245
#> - xpsr normalization scale = 0.56671, estimated xpsr sd = 0.73647
fit_true <- create_oracle_fit(fit_prelim, beta_true, h0_true)
fit_true$covariate_effects()
#> covariate beta target_state_idx target_state
#> 1 age 0.0 ± NA 2 Bleed
#> 2 age 0.0 ± NA 3 Stroke
#> 3 age 0.3 ± NA 4 Dead
#> 4 xpsr 1.0 ± NA 2 Bleed
#> 5 xpsr -1.0 ± NA 3 Stroke
#> 6 xpsr 0.0 ± NA 4 Dead
fit_true$plot_pk()
#> PK simulation
fit_true$plot_h0()
Fitting various models
fit_ms_eh <- fit_stan(mod_ms_eh, simdat, method = "optimize", init = 0)
#> Shortest time interval (0.0359547288841782) is smaller than delta_grid (1.09575). Consider increasing n_grid or decreasing t_max of the model.
#> Using stan file at /home/runner/work/_temp/Library/bmstate/stan/msm.stan
#> setting xpsr normalizers to loc = 5.42315, scale = 0.56671
#> setting max conc = 9565.56038
#> Initial log joint probability = -11156.4
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 99 -7717.64 0.0729145 1171.84 1 1 109
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 199 -7132.53 0.0404134 870.515 1 1 229
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 299 -6913.43 0.00734514 1507.19 1 1 340
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 399 -6735.44 0.00503622 2088.9 1 1 453
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 499 -6660.31 0.0355828 859.492 1 1 565
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 599 -6584.9 0.0163218 852.574 1 1 675
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 699 -6529.04 0.00542709 465.525 1 1 785
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 799 -6503.43 0.0055988 951.853 0.3624 1 898
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 899 -6483.59 0.0265756 1063.93 1 1 1006
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 999 -6460.4 0.0035246 552.015 1 1 1116
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 1099 -6448.94 0.00170797 289.855 1 1 1227
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 1199 -6441.94 0.000532439 327.47 0.3737 0.3737 1339
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 1299 -6435.49 0.00150055 269.552 1 1 1446
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 1399 -6431.76 0.00254161 192.315 1 1 1557
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 1499 -6424.11 0.00085538 490.992 0.8571 0.8571 1662
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 1599 -6412.74 0.0116704 646.135 1 1 1772
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 1699 -6404.36 0.00407611 350.42 1 1 1883
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 1799 -6398.57 0.0016371 234.811 1 1 1990
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 1899 -6397.5 0.00183939 142.567 1 1 2096
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 1999 -6396.44 0.000560921 95.6849 1 1 2209
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 2000 -6396.43 0.00056563 231.704 1 1 2210
#> Optimization terminated normally:
#> Maximum number of iterations hit, may not be at an optima
#> Finished in 4.6 seconds.
#> If these don't roughly match, consider refitting after setting xpsr normalizer loc and scale closer to estimated mean and sd. Otherwise interpret baseline hazards and xpsr effect size accodingly.
#> - xpsr normalization loc = 5.42315, mean estimated xpsr = 5.50741
#> - xpsr normalization scale = 0.56671, estimated xpsr sd = 0.96188
fit_ms_dh <- fit_stan(mod_ms_dh, simdat_dh, method = "optimize", init = 0)
#> Shortest time interval (0.0336870958247317) is smaller than delta_grid (1.09575). Consider increasing n_grid or decreasing t_max of the model.
#> Using stan file at /home/runner/work/_temp/Library/bmstate/stan/msm.stan
#> Initial log joint probability = -7578.43
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 99 -6096.63 0.00509309 3.96255 1 1 106
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 199 -6095.58 0.0438785 2.63666 1 1 216
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 299 -6095.38 0.00517853 0.72308 1 1 326
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 399 -6095.36 0.0116761 0.575551 1 1 431
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 499 -6095.32 0.0270586 1.26346 1 1 539
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 599 -6095.3 0.00505974 0.578409 1 1 641
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 626 -6095.29 0.00145307 0.0857739 1 1 671
#> Optimization terminated normally:
#> Convergence detected: relative gradient magnitude is below tolerance
#> Finished in 0.4 seconds.
fit_death <- fit_stan(mod_death, simdat_death, method = "optimize", init = 0)
#> Using stan file at /home/runner/work/_temp/Library/bmstate/stan/msm.stan
#> Initial log joint probability = -3594.95
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 99 -3379.61 0.0192462 0.637589 1 1 107
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 199 -3379.46 0.0480659 0.416013 0.8681 0.8681 218
#> Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
#> 260 -3379.45 0.000684573 0.0415297 0.8453 0.8453 284
#> Optimization terminated normally:
#> Convergence detected: relative gradient magnitude is below tolerance
#> Finished in 0.1 seconds.
fit_ms_eh$plot_h0()
fit_ms_eh$covariate_effects()
#> covariate beta target_state_idx target_state
#> 1 age 0.21 ± NA 2 Bleed
#> 2 age -0.14 ± NA 3 Stroke
#> 3 age 0.23 ± NA 4 Dead
#> 4 xpsr 0.14 ± NA 2 Bleed
#> 5 xpsr -0.45 ± NA 3 Stroke
#> 6 xpsr 0.02 ± NA 4 Dead
fit_ms_dh$covariate_effects()
#> covariate beta target_state_idx target_state
#> 1 age 0.705 ± NA 2 Bleed
#> 2 age -0.619 ± NA 3 Stroke
#> 3 age 0.260 ± NA 4 Dead
#> 4 CrCL -0.602 ± NA 2 Bleed
#> 5 CrCL 0.514 ± NA 3 Stroke
#> 6 CrCL 0.073 ± NA 4 Dead
#> 7 weight 0.086 ± NA 2 Bleed
#> 8 weight 0.022 ± NA 3 Stroke
#> 9 weight 0.033 ± NA 4 Dead
#> 10 dose_amt 1.020 ± NA 2 Bleed
#> 11 dose_amt -1.168 ± NA 3 Stroke
#> 12 dose_amt 0.028 ± NA 4 Dead
fit_death$covariate_effects()
#> covariate beta target_state_idx target_state
#> 1 age -0.126 ± NA 2 Dead
#> 2 CrCL 0.181 ± NA 2 Dead
#> 3 weight 0.043 ± NA 2 Dead
#> 4 dose_amt -0.222 ± NA 2 Dead
fit_ms_eh$plot_pk()
#> PK simulation