From 99662ae1b50473e528de2065e73b3636f8a47b96 Mon Sep 17 00:00:00 2001 From: jtimonen Date: Tue, 25 Nov 2025 02:19:51 +0200 Subject: [PATCH 1/6] Increment version number to 0.2.8.9000 --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9775d40..43733ab 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: bmstate Type: Package Title: Bayesian multistate modeling -Version: 0.2.8 +Version: 0.2.8.9000 Authors@R: c(person(given = "Juho", family = "Timonen", From 67e37fbe2b8732d38b526243b915553f565814c6 Mon Sep 17 00:00:00 2001 From: jtimonen Date: Wed, 3 Dec 2025 23:32:57 +0200 Subject: [PATCH 2/6] use concentration auc instead of amount auc --- R/MultistateModel.R | 4 ++-- R/PKModel.R | 6 ++++-- R/stan.R | 5 +++-- inst/stan/msm.stan | 4 ++-- 4 files changed, 11 insertions(+), 8 deletions(-) diff --git a/R/MultistateModel.R b/R/MultistateModel.R index 245b054..632e933 100644 --- a/R/MultistateModel.R +++ b/R/MultistateModel.R @@ -47,8 +47,8 @@ MultistateModel <- R6::R6Class("MultistateModel", categorical = NULL, normalizer_locations = NULL, normalizer_scales = NULL, - auc_normalizer_loc = 300, - auc_normalizer_scale = 100, + auc_normalizer_loc = 2000, + auc_normalizer_scale = 1000, n_grid = NULL, simulate_log_hazard_multipliers = function(df_subjects, beta) { ts <- self$target_states() diff --git a/R/PKModel.R b/R/PKModel.R index 997156d..4136bb4 100644 --- a/R/PKModel.R +++ b/R/PKModel.R @@ -102,16 +102,18 @@ PKModel <- R6::R6Class("PKModel", conc }, - #' @description Compute steady-state area under curve over one dosing interval + #' @description Compute steady-state area under concentration curve over + #' one dosing interval #' #' @param theta parameter values #' @param dose dose #' @return A numeric value compute_ss_auc = function(theta, dose) { CL <- theta[2] + V2 <- theta[3] checkmate::assert_number(CL, lower = 0) checkmate::assert_number(dose, lower = 0) - dose / CL + dose / (CL * V2) }, #' @description Simulate data with many subjects diff --git a/R/stan.R b/R/stan.R index a665556..aff4304 100644 --- a/R/stan.R +++ b/R/stan.R @@ -101,8 +101,9 @@ fit_stan <- function(model, data, if (set_normalizers) { model$set_normalizers(data) if (!is.null(data$dosing)) { - mu_CL <- exp(-2) - aaa <- data$dosing$dose_ss / mu_CL + mu_CL <- exp(-2) # should match msm.stan + mu_V2 <- exp(-2) # should match msm.stan + aaa <- data$dosing$dose_ss / (mu_CL * mu_V2) loc <- mean(aaa) sca <- stats::sd(aaa) model$set_auc_normalizers(loc, sca) diff --git a/inst/stan/msm.stan b/inst/stan/msm.stan index 742c524..357bfa8 100644 --- a/inst/stan/msm.stan +++ b/inst/stan/msm.stan @@ -427,8 +427,8 @@ transformed parameters { MAX_CONC ); - // Auc at steady state - ss_auc[1] = dose_ss ./ theta_pk[1][:,2]; // D/CL + // Concentration AUC at steady state + ss_auc[1] = (dose_ss ./ (theta_pk[1][:,2] .* theta_pk[1][:,3])); // D/(CL*V2) // Set AUC corresponding to each interval for(n in 1:N_int){ From 05c51c7bdab5c22be1df9dfb34c0f0503a3f204c Mon Sep 17 00:00:00 2001 From: jtimonen Date: Wed, 3 Dec 2025 23:42:06 +0200 Subject: [PATCH 3/6] fix vignette --- man/PKModel.Rd | 3 ++- vignettes/math.Rmd | 12 ++++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/man/PKModel.Rd b/man/PKModel.Rd index e035db1..a8a270c 100644 --- a/man/PKModel.Rd +++ b/man/PKModel.Rd @@ -154,7 +154,8 @@ concentration in the central compartment. \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-PKModel-compute_ss_auc}{}}} \subsection{Method \code{compute_ss_auc()}}{ -Compute steady-state area under curve over one dosing interval +Compute steady-state area under concentration curve over +one dosing interval \subsection{Usage}{ \if{html}{\out{
}}\preformatted{PKModel$compute_ss_auc(theta, dose)}\if{html}{\out{
}} } diff --git a/vignettes/math.Rmd b/vignettes/math.Rmd index 3491b56..bf3cdb2 100644 --- a/vignettes/math.Rmd +++ b/vignettes/math.Rmd @@ -300,8 +300,7 @@ print(mod) ``` Now we can simulate death events with known constant hazard rate `lambda`. For this -model, we are simulating from a homogeneous Poisson process with rate $\lambda$, -meaning that the event times should follow an exponential distribution with +model, the event times should follow an exponential distribution with mean $\frac{1}{\lambda}$. ```{r, fig.width=5} @@ -360,10 +359,11 @@ mod <- create_msm(tm, t_max = tmax) ``` This model has 3 possible transitions. -Now we can simulate events with known constant hazard rates `lambda_1`, -`lambda_2`, `lambda_3`. For this -model, we are simulating from independent homogeneous Poisson processes with -rates $\lambda_i$, meaning that the transition probabilities should be +Now we can simulate transition times with known constant hazard rates `lambda_1`, +`lambda_2`, `lambda_3`. Simulating the transition times is mathematically equivalent to +simulating from independent homogeneous Poisson processes with +rates $\lambda_i$, and taking the transition that occurred first. This means that +the transition probabilities should be $$ P(i) = \frac{\lambda_i}{\sum_{j=1}^ 3 \lambda_j}. $$ From bf1bc8a2a5f9431b6d858eff20c3d2707e9abafd Mon Sep 17 00:00:00 2001 From: jtimonen Date: Wed, 3 Dec 2025 23:46:09 +0200 Subject: [PATCH 4/6] allow plotting pk for certain subjects only --- R/DosingData.R | 12 ++++++++++-- R/MultistateModelFit.R | 11 +++++++++-- man/MultistateModelFit.Rd | 5 ++++- man/PSSDosingData.Rd | 9 ++++++++- 4 files changed, 31 insertions(+), 6 deletions(-) diff --git a/R/DosingData.R b/R/DosingData.R index 835f9df..3c2ef1b 100644 --- a/R/DosingData.R +++ b/R/DosingData.R @@ -136,7 +136,9 @@ PSSDosingData <- R6::R6Class( #' and \code{upper}. #' @param subject_df Subject data frame. #' @param max_num_subjects Max number of subjects to plot. - plot = function(df_fit = NULL, subject_df = NULL, max_num_subjects = 12) { + #' @param subject_ids Which subjects to plot? + plot = function(df_fit = NULL, subject_df = NULL, max_num_subjects = 12, + subject_ids = NULL) { dos <- self$as_data_frame() fitcolor <- "steelblue" if (is.null(max_num_subjects)) { @@ -144,7 +146,13 @@ PSSDosingData <- R6::R6Class( } checkmate::assert_integerish(max_num_subjects, len = 1, lower = 1) if (self$num_subjects() > max_num_subjects) { - sid <- sample(unique(dos$subject_id), max_num_subjects) + if (is.null(subject_ids)) { + sid <- sample(unique(dos$subject_id), max_num_subjects) + } else { + checkmate::assert_character(subject_ids) + sid <- subject_ids + } + dos <- dos |> dplyr::filter(.data$subject_id %in% sid) if (!is.null(df_fit)) { df_fit <- df_fit |> dplyr::filter(.data$subject_id %in% sid) diff --git a/R/MultistateModelFit.R b/R/MultistateModelFit.R index 78f1862..82e22ff 100644 --- a/R/MultistateModelFit.R +++ b/R/MultistateModelFit.R @@ -169,8 +169,10 @@ MultistateModelFit <- R6::R6Class("MultistateModelFit", #' @param n_prev number of previous doses to show fit for #' @param oos Out-of-sample subjects? #' @param ci_alpha Width of central credible interval. + #' @param subject_ids Which subjects to plot? plot_pk = function(max_num_subjects = 12, oos = FALSE, data = NULL, L = 100, - timescale = 24, n_prev = 3, ci_alpha = 0.9) { + timescale = 24, n_prev = 3, ci_alpha = 0.9, + subject_ids = NULL) { pksim <- self$simulate_pk(oos, data, L, timescale, n_prev) if (is.null(data)) { data <- self$data @@ -182,13 +184,18 @@ MultistateModelFit <- R6::R6Class("MultistateModelFit", } checkmate::assert_number(ci_alpha, lower = 0, upper = 1) av <- (1 - ci_alpha) / 2 + + # Slow pksim <- pksim |> dplyr::group_by(.data$subject_id, .data$time) |> dplyr::summarise(q = list(quantile(.data$val, probs = c(av / 2, 0.5, 1 - av / 2))), .groups = "drop") |> tidyr::unnest_wider(q, names_sep = "_") colnames(pksim)[3:5] <- c("lower", "val", "upper") - data$plot_dosing(df_fit = pksim, max_num_subjects = max_num_subjects) + + data$plot_dosing( + df_fit = pksim, max_num_subjects = max_num_subjects, + subject_ids = subject_ids + ) + labs(caption = capt) }, diff --git a/man/MultistateModelFit.Rd b/man/MultistateModelFit.Rd index b574e52..0a6e799 100644 --- a/man/MultistateModelFit.Rd +++ b/man/MultistateModelFit.Rd @@ -242,7 +242,8 @@ Plot PK fit. L = 100, timescale = 24, n_prev = 3, - ci_alpha = 0.9 + ci_alpha = 0.9, + subject_ids = NULL )}\if{html}{\out{}} } @@ -263,6 +264,8 @@ training data is used.} \item{\code{n_prev}}{number of previous doses to show fit for} \item{\code{ci_alpha}}{Width of central credible interval.} + +\item{\code{subject_ids}}{Which subjects to plot?} } \if{html}{\out{}} } diff --git a/man/PSSDosingData.Rd b/man/PSSDosingData.Rd index 4811553..6ff85ea 100644 --- a/man/PSSDosingData.Rd +++ b/man/PSSDosingData.Rd @@ -116,7 +116,12 @@ a \code{data.frame} \subsection{Method \code{plot()}}{ Plot dosing (and PK) data \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{PSSDosingData$plot(df_fit = NULL, subject_df = NULL, max_num_subjects = 12)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{PSSDosingData$plot( + df_fit = NULL, + subject_df = NULL, + max_num_subjects = 12, + subject_ids = NULL +)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -128,6 +133,8 @@ and \code{upper}.} \item{\code{subject_df}}{Subject data frame.} \item{\code{max_num_subjects}}{Max number of subjects to plot.} + +\item{\code{subject_ids}}{Which subjects to plot?} } \if{html}{\out{}} } From 47f2e4a37c707f77d0ebea957b36a894ed5b56ae Mon Sep 17 00:00:00 2001 From: jtimonen Date: Wed, 3 Dec 2025 23:55:43 +0200 Subject: [PATCH 5/6] fix --- R/JointData.R | 5 +++-- R/MultistateModelFit.R | 14 ++++---------- R/utils.R | 12 ++++++++++++ man/JointData.Rd | 4 +++- 4 files changed, 22 insertions(+), 13 deletions(-) diff --git a/R/JointData.R b/R/JointData.R index ad84152..bfee466 100644 --- a/R/JointData.R +++ b/R/JointData.R @@ -66,13 +66,14 @@ JointData <- R6::R6Class( #' @description Plot dosing #' @param df_fit Fit data frame. #' @param max_num_subjects Max number of subjects to plot. + #' @param subject_ids Which subjects to plot? #' @return a \code{ggplot} - plot_dosing = function(df_fit = NULL, max_num_subjects = 12) { + plot_dosing = function(df_fit = NULL, max_num_subjects = 12, subject_ids = NULL) { if (is.null(self$dosing)) { stop("No dosing data.") } sdf <- self$paths$subject_df - self$dosing$plot(df_fit, sdf, max_num_subjects) + self$dosing$plot(df_fit, sdf, max_num_subjects, subject_ids) } ) ) diff --git a/R/MultistateModelFit.R b/R/MultistateModelFit.R index 82e22ff..35d3d1f 100644 --- a/R/MultistateModelFit.R +++ b/R/MultistateModelFit.R @@ -182,21 +182,15 @@ MultistateModelFit <- R6::R6Class("MultistateModelFit", } else { capt <- paste0("Point-estimate fit") } - checkmate::assert_number(ci_alpha, lower = 0, upper = 1) - av <- (1 - ci_alpha) / 2 - # Slow - pksim <- pksim |> - dplyr::group_by(.data$subject_id, .data$time) |> - dplyr::summarise(q = list(quantile(.data$val, probs = c(av / 2, 0.5, 1 - av / 2))), .groups = "drop") |> - tidyr::unnest_wider(q, names_sep = "_") - colnames(pksim)[3:5] <- c("lower", "val", "upper") + # Slow for some reason + pksim <- pksim_to_quantiles(pksim, ci_alpha) + # Create plot data$plot_dosing( df_fit = pksim, max_num_subjects = max_num_subjects, subject_ids = subject_ids - ) + - labs(caption = capt) + ) + labs(caption = capt) }, #' @description Plot baseline hazard distribution diff --git a/R/utils.R b/R/utils.R index ec5788a..7210e14 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1,3 +1,15 @@ +# Helper +pksim_to_quantiles <- function(sim, ci_alpha) { + checkmate::assert_number(ci_alpha, lower = 0, upper = 1) + av <- (1 - ci_alpha) / 2 + sim <- sim |> + dplyr::group_by(.data$subject_id, .data$time) |> + dplyr::summarise(q = list(quantile(.data$val, probs = c(av / 2, 0.5, 1 - av / 2))), .groups = "drop") |> + tidyr::unnest_wider(q, names_sep = "_") + colnames(sim)[3:5] <- c("lower", "val", "upper") + sim +} + # Check data frame columns check_columns <- function(df, needed_columns) { checkmate::assert_data_frame(df) diff --git a/man/JointData.Rd b/man/JointData.Rd index f191aaf..af30df3 100644 --- a/man/JointData.Rd +++ b/man/JointData.Rd @@ -78,7 +78,7 @@ Print info \subsection{Method \code{plot_dosing()}}{ Plot dosing \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{JointData$plot_dosing(df_fit = NULL, max_num_subjects = 12)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{JointData$plot_dosing(df_fit = NULL, max_num_subjects = 12, subject_ids = NULL)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -87,6 +87,8 @@ Plot dosing \item{\code{df_fit}}{Fit data frame.} \item{\code{max_num_subjects}}{Max number of subjects to plot.} + +\item{\code{subject_ids}}{Which subjects to plot?} } \if{html}{\out{}} } From 55e69fe8a3605b5c462735251060c76d9c3b620d Mon Sep 17 00:00:00 2001 From: jtimonen Date: Thu, 4 Dec 2025 00:40:32 +0200 Subject: [PATCH 6/6] Increment version number to 0.2.9 --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 43733ab..46e9a4d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: bmstate Type: Package Title: Bayesian multistate modeling -Version: 0.2.8.9000 +Version: 0.2.9 Authors@R: c(person(given = "Juho", family = "Timonen",