Merge pull request #48 from generable/develop

v0.3.0

Author: jtimonen

DESCRIPTION

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

R/MultistateModel.R

@@ -47,13 +47,13 @@ MultistateModel <- R6::R6Class("MultistateModel",
     categorical = NULL,
     normalizer_locations = NULL,
     normalizer_scales = NULL,
-    auc_normalizer_loc = 2000,
-    auc_normalizer_scale = 1000,
+    xpsr_normalizer_loc = 7.5,
+    xpsr_normalizer_scale = 0.5,
     n_grid = NULL,
     simulate_log_hazard_multipliers = function(df_subjects, beta) {
       ts <- self$target_states()
       x <- self$covs()
-      auc_norm <- self$get_auc_normalizers()
+      xpsr_norm <- self$get_xpsr_normalizers()
       B <- length(ts)
       K <- length(x)
       checkmate::assert_matrix(beta, nrows = B, ncols = K)
@@ -64,10 +64,10 @@ MultistateModel <- R6::R6Class("MultistateModel",
       X <- df_subjects |> dplyr::select(tidyselect::all_of(x))
       X <- as.matrix(X)
       X_norm <- normalize_columns(X)
-      if ("ss_auc" %in% x) {
-        idx <- which(x == "ss_auc")
-        x_norm <- (X[, idx] - auc_norm$loc) / auc_norm$scale
-        check_normalized_covariate(x_norm, "ss_auc")
+      if ("xpsr" %in% x) {
+        idx <- which(x == "xpsr")
+        x_norm <- (X[, idx] - xpsr_norm$loc) / xpsr_norm$scale
+        check_normalized_covariate(x_norm, "xpsr")
         X_norm[, idx] <- x_norm
       }
       for (s in seq_len(S)) {
@@ -149,28 +149,28 @@ MultistateModel <- R6::R6Class("MultistateModel",
       invisible(NULL)
     },
 
-    #' @description Get normalization constants for AUC (PK)
+    #' @description Get normalization constants for exposure (PK)
     #' @return list
-    get_auc_normalizers = function() {
+    get_xpsr_normalizers = function() {
       list(
-        loc = private$auc_normalizer_loc,
-        scale = private$auc_normalizer_scale
+        loc = private$xpsr_normalizer_loc,
+        scale = private$xpsr_normalizer_scale
       )
     },
 
-    #' @description Set normalization constants for AUC (side effect)
+    #' @description Set normalization constants for exposure (side effect)
     #'
     #' @param loc Location
     #' @param scale Scale
-    set_auc_normalizers = function(loc = 0, scale = 1) {
+    set_xpsr_normalizers = function(loc = 0, scale = 1) {
       checkmate::assert_numeric(loc, lower = 0, len = 1)
       checkmate::assert_numeric(scale, lower = 0, len = 1)
       message(
-        "setting auc normalizers to loc = ",
+        "setting xpsr normalizers to loc = ",
         round(loc, 5), ", scale = ", round(scale, 5)
       )
-      private$auc_normalizer_loc <- loc
-      private$auc_normalizer_scale <- scale
+      private$xpsr_normalizer_loc <- loc
+      private$xpsr_normalizer_scale <- scale
       invisible(NULL)
     },
 
@@ -211,7 +211,7 @@ MultistateModel <- R6::R6Class("MultistateModel",
     #' @param system A \code{\link{MultistateSystem}}
     #' @param covariates The names of the hazard covariates (excluding possible
     #' exposure estimated from PK model). Do not use reserved names
-    #' \code{ss_auc} or \code{dose}.
+    #' \code{xpsr} or \code{dose}.
     #' @param pk_model A \code{\link{PKModel}} or NULL.
     #' @param t_max Max time.
     #' @param num_knots Total number of spline knots.
@@ -232,7 +232,7 @@ MultistateModel <- R6::R6Class("MultistateModel",
       if (!all(categorical %in% covariates)) {
         stop("all categorical covariates should be also in covariates")
       }
-      checkmate::assert_true(!("ss_auc" %in% covariates)) # special name
+      checkmate::assert_true(!("xpsr" %in% covariates)) # special name
       checkmate::assert_true(!("dose" %in% covariates)) # special name
       checkmate::assert_class(system, "MultistateSystem")
       checkmate::assert_integerish(n_grid, lower = 10, len = 1)
@@ -310,7 +310,7 @@ MultistateModel <- R6::R6Class("MultistateModel",
     covs = function() {
       x <- private$hazard_covariates
       if (self$has_pk()) {
-        x <- c(x, "ss_auc")
+        x <- c(x, "xpsr")
       }
       unique(x)
     },

R/MultistateModelFit.R

@@ -292,9 +292,9 @@ msmfit_covariate_effects <- function(fit) {
     df <- rbind(df, df_j)
   }
   if (fit$model$has_pk()) {
-    rv <- as.vector(fit$get_draws("beta_auc")[1, 1, ])
+    rv <- as.vector(fit$get_draws("beta_xpsr")[1, 1, ])
     df_j <- data.frame(
-      covariate = "ss_auc", beta = rv,
+      covariate = "xpsr", beta = rv,
       target_state_idx = fit$model$target_states()
     )
     df <- rbind(df, df_j)
@@ -344,9 +344,11 @@ mat2list <- function(mat) {
 #' @param oos Out-of-sample mode? If \code{FALSE}, the possible subject-specific
 #' fitted parameters are used. If \code{TRUE}, acting
 #' as if the subjects are new.
+#' @param log Get logarithm of params?
 #' @return A list with length equal to number of draws.
-msmfit_pk_params <- function(fit, oos = FALSE, data = NULL) {
+msmfit_pk_params <- function(fit, oos = FALSE, data = NULL, log = FALSE) {
   check_oos(oos, data)
+  checkmate::assert_logical(log, len = 1)
   sd <- msmfit_stan_data(fit, data)
   S <- fit$num_draws()
 
@@ -372,9 +374,9 @@ msmfit_pk_params <- function(fit, oos = FALSE, data = NULL) {
   # Call exposed Stan function for each draw (not optimal)
   out <- list()
   for (s in seq_len(S)) {
-    theta <- NULL
+    log_theta <- NULL
     if (sd$do_pk == 1) {
-      theta <- compute_theta_pk(
+      log_theta <- compute_log_theta_pk(
         mat2list(t(log_z[s, 1, , ])),
         log_mu[s, 1, ],
         log_sig[s, 1, ],
@@ -386,6 +388,11 @@ msmfit_pk_params <- function(fit, oos = FALSE, data = NULL) {
         mat2list(t(sd$x_V2))
       )
     }
+    if (log) {
+      theta <- log_theta
+    } else {
+      theta <- exp(log_theta)
+    }
     out[[s]] <- theta
   }
   out
@@ -400,18 +407,18 @@ msmfit_exposure <- function(fit, oos = FALSE, data = NULL) {
 
   # Get draws
   sd <- msmfit_stan_data(fit, data)
-  pkpar <- msmfit_pk_params(fit, oos, data)
+  log_pkpar <- msmfit_pk_params(fit, oos, data, log = TRUE)
 
   # Call exposed Stan function
   S <- fit$num_draws()
   out <- list()
   for (s in seq_len(S)) {
     if (sd$do_pk == 1) {
-      x_auc <- sd$dose_ss / (pkpar[[s]][, 2] * pkpar[[s]][, 3]) # D/(CL*V2)
+      x_xpsr <- log_ss_area_under_conc(sd$dose_ss, log_pkpar[[s]]) # log D/(CL*V2)
     } else {
-      x_auc <- NULL
+      x_xpsr <- NULL
     }
-    out[[s]] <- x_auc
+    out[[s]] <- x_xpsr
   }
   out
 }
@@ -445,16 +452,16 @@ msmfit_log_hazard_multipliers <- function(fit, oos = FALSE, data = NULL) {
 
   # Get draws
   sd <- msmfit_stan_data(fit, data)
-  auc <- msmfit_exposure(fit, oos, data)
+  xpsr <- msmfit_exposure(fit, oos, data)
   S <- fit$num_draws()
   beta_oth <- fit$get_draws_of("beta_oth")
   if (is.null(beta_oth)) {
     beta_oth <- array(0, dim = c(S, 1, 0, sd$N_trans_types))
   }
   if (sd$do_pk == 1) {
-    beta_auc <- fit$get_draws_of("beta_auc")
+    beta_xpsr <- fit$get_draws_of("beta_xpsr")
   } else {
-    beta_auc <- array(0, dim = c(S, 1, 0, sd$N_trans_types))
+    beta_xpsr <- array(0, dim = c(S, 1, 0, sd$N_trans_types))
   }
 
   # Create x_haz_long (long version of hazard covariates vector)
@@ -465,14 +472,14 @@ msmfit_log_hazard_multipliers <- function(fit, oos = FALSE, data = NULL) {
   } else {
     x_haz_long <- array(0, dim = c(0, sd$N_int))
   }
-  an <- fit$model$get_auc_normalizers()
+  an <- fit$model$get_xpsr_normalizers()
 
   # Call exposed Stan function for each draw (not optimal)
   out <- list()
   for (s in seq_len(S)) {
     if (sd$do_pk == 1) {
-      ba <- list(beta_auc[s, 1, 1, ])
-      aa <- list(auc[[s]][sd$idx_sub])
+      ba <- list(beta_xpsr[s, 1, 1, ])
+      aa <- list(xpsr[[s]][sd$idx_sub])
     } else {
       ba <- NULL
       aa <- NULL

R/PKModel.R

@@ -102,18 +102,18 @@ PKModel <- R6::R6Class("PKModel",
       conc
     },
 
-    #' @description Compute steady-state area under concentration curve over
-    #' one dosing interval
+    #' @description Compute exposure, which is the steady-state
+    #' log 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) {
+    compute_xpsr = function(theta, dose) {
       CL <- theta[2]
       V2 <- theta[3]
       checkmate::assert_number(CL, lower = 0)
       checkmate::assert_number(dose, lower = 0)
-      dose / (CL * V2)
+      log(dose) - log(CL * V2)
     },
 
     #' @description Simulate data with many subjects
@@ -159,17 +159,17 @@ PKModel <- R6::R6Class("PKModel",
       CONC <- dd$simulate_pk(t_obs, THETA, self$get_max_conc())
       df_out <- NULL
       for (n in seq_len(N)) {
-        ss_auc <- self$compute_ss_auc(THETA[n, ], dd$dose_ss[n])
+        xpsr <- self$compute_xpsr(THETA[n, ], dd$dose_ss[n])
         sid <- SUB_ID[n]
         conc <- (CONC |> dplyr::filter(.data$subject_id == sid))[["val"]]
         conc_noisy <- stats::rlnorm(2, meanlog = log(conc), sdlog = sigma)
-        out <- c(t_obs[[n]], conc_noisy, ss_auc)
+        out <- c(t_obs[[n]], conc_noisy, xpsr)
         df_out <- rbind(df_out, out)
       }
       df_out <- data.frame(df_out)
       df_out <- cbind(df_out, THETA)
       colnames(df_out) <- c(
-        "t_pre", "t_post", "conc_pre", "conc_post", "ss_auc",
+        "t_pre", "t_post", "conc_pre", "conc_post", "xpsr",
         "ka", "CL", "V2"
       )
       rownames(df_out) <- NULL

R/stan-data.R

@@ -40,7 +40,7 @@ create_stan_data_model <- function(model) {
   sd_t_grid <- create_stan_data_timegrid(model)
 
   # PK options
-  an <- model$get_auc_normalizers()
+  an <- model$get_xpsr_normalizers()
   MC <- 1e-7
   if (model$has_pk()) {
     MC <- model$pk_model$get_max_conc()
@@ -49,9 +49,9 @@ create_stan_data_model <- function(model) {
     nc_ka = length(model$data_covs("ka")),
     nc_CL = length(model$data_covs("CL")),
     nc_V2 = length(model$data_covs("V2")),
-    I_auc = as.numeric(model$has_pk()),
-    auc_loc = an$loc,
-    auc_scale = an$scale,
+    I_xpsr = as.numeric(model$has_pk()),
+    xpsr_loc = an$loc,
+    xpsr_scale = an$scale,
     MAX_CONC = MC
   )
 

R/stan.R

@@ -101,12 +101,12 @@ fit_stan <- function(model, data,
   if (set_normalizers) {
     model$set_normalizers(data)
     if (!is.null(data$dosing)) {
-      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)
+      log_mu_CL <- -2 # should match msm.stan
+      log_mu_V2 <- -2 # should match msm.stan
+      aaa <- log(data$dosing$dose_ss) - log_mu_CL - log_mu_V2
       loc <- mean(aaa)
       sca <- stats::sd(aaa)
-      model$set_auc_normalizers(loc, sca)
+      model$set_xpsr_normalizers(loc, sca)
     }
   }
 
@@ -133,7 +133,7 @@ fit_stan <- function(model, data,
   # Return
   pars <- c(
     "weights", "log_w0", "beta_ka", "beta_V2", "beta_CL", "beta_oth",
-    "beta_auc", "sigma_pk", "log_z_pk", "log_mu_pk", "log_sig_pk", "lp__"
+    "beta_xpsr", "sigma_pk", "log_z_pk", "log_mu_pk", "log_sig_pk", "lp__"
   )
   draws <- create_rv_list(stan_fit, pars)
   diag <- NULL