v0.2.9

DESCRIPTION

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

R/DosingData.R

@@ -136,15 +136,23 @@ 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)) {
         max_num_subjects <- self$num_subjects()
       }
       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)

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)
     }
   )
 )

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()

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
@@ -180,16 +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
-      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) +
-        labs(caption = capt)
+
+      # 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)
     },
 
     #' @description Plot baseline hazard distribution

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

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)

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)

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){

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