Skip to content

v0.3.1 #49

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jtimonen merged 4 commits into from
Dec 5, 2025
Merged

v0.3.1 #49

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
01ab476
Increment version number to 0.3.0.9000
jtimonen Dec 5, 2025
3e345c9
debug and add warning message
jtimonen Dec 5, 2025
7d43200
better errors
jtimonen Dec 5, 2025
97e3430
Increment version number to 0.3.1
jtimonen Dec 5, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: bmstate
Type: Package
Title: Bayesian multistate modeling
Version: 0.3.0
Version: 0.3.1
Authors@R:
c(person(given = "Juho",
family = "Timonen",
Expand Down
2 changes: 2 additions & 0 deletions NAMESPACE
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,9 @@ export(ensure_exposed_stan_functions)
export(example_sim_setup_illnessdeath)
export(fit_stan)
export(generate_paths)
export(msmfit_check_xpsr_norm)
export(msmfit_exposure)
export(msmfit_exposure_df)
export(msmfit_inst_hazard_param_draws)
export(msmfit_log_hazard_multipliers)
export(msmfit_pk_params)
Expand Down
17 changes: 4 additions & 13 deletions R/DosingData.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,23 +104,14 @@ PSSDosingData <- R6::R6Class(
#'
#' @param t A vector of output times for each subject (a list).
#' @param theta A matrix of parameters.
#' @param skip_assert Skip most assertions and call exposed Stan directly,
#' assuming that it exists?
#' @param MAX_CONC concentration upper bound
#' @return a \code{data.frame}
simulate_pk = function(t, theta, MAX_CONC, skip_assert = FALSE) {
simulate_pk = function(t, theta, MAX_CONC) {
checkmate::assert_list(t, len = self$num_subjects())
checkmate::assert_logical(skip_assert, len = 1)
checkmate::assert_number(MAX_CONC, lower = 0)
if (skip_assert) {
out <- pop_2cpt_partly_ss(
t, self$dose_ss, self$times, self$doses, theta, self$tau_ss, MAX_CONC
)
} else {
out <- pk_2cpt_pss(
t, self$dose_ss, self$times, self$doses, theta, self$tau_ss, MAX_CONC
)
}
out <- pk_2cpt_pss(
t, self$dose_ss, self$times, self$doses, theta, self$tau_ss, MAX_CONC
)

time <- as.numeric(unlist(t))
val <- unlist(out)
Expand Down
21 changes: 0 additions & 21 deletions R/MultistateModel.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -474,24 +474,3 @@ place_internal_knots <- function(t_max, num_knots, t_event) {
knots <- stats::quantile(t_event, probs = seq(0, 1, h))
knots[2:(length(knots) - 1)]
}

# Normalize columns of matrix A
normalize_columns <- function(A) {
K <- ncol(A)
cn <- colnames(A)
for (j in seq_len(K)) {
A[, j] <- normalize(A[, j], cn[j])
}
A
}

# Normalize
normalize <- function(a, name) {
sdd <- stats::sd(a)
if (sdd == 0) {
stop("error in normalization of ", name, ", zero variance")
}
x_norm <- (a - mean(a)) / sdd
check_normalized_covariate(x_norm, name)
x_norm
}
46 changes: 45 additions & 1 deletion R/MultistateModelFit.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ MultistateModelFit <- R6::R6Class("MultistateModelFit",
DF <- list()
for (s in seq_len(S)) {
pb$tick()
df_s <- data$dosing$simulate_pk(ts, pkpar[[s]], MC, skip_assert = TRUE)
df_s <- data$dosing$simulate_pk(ts, pkpar[[s]], MC)
df_s$.draw_idx <- s
DF[[s]] <- df_s
}
Expand Down Expand Up @@ -423,6 +423,50 @@ msmfit_exposure <- function(fit, oos = FALSE, data = NULL) {
out
}

#' Compute exposure and return as rvar in df
#'
#' @export
#' @inheritParams msmfit_exposure
msmfit_exposure_df <- function(fit, oos = FALSE, data = NULL) {
check_oos(oos, data)
xp <- msmfit_exposure(fit, oos, data)
xp <- posterior::rvar(t(sapply(xp, function(x) x)))
if (is.null(data)) {
data <- fit$data
}
df <- data$paths$subject_df
df$xpsr_estimate <- xp
df
}

#' Check exposure normalization
#'
#' @export
#' @inheritParams msmfit_exposure
msmfit_check_xpsr_norm <- function(fit) {
df <- msmfit_exposure_df(fit)
r <- mean(df$xpsr_estimate)
xn <- fit$model$get_xpsr_normalizers()
msg1 <- paste0(
"xpsr normalization loc = ", round(xn$loc, 5),
", mean estimated xpsr = ",
round(mean(r), 5)
)
msg2 <- paste0(
"xpsr normalization scale = ", round(xn$scale, 5),
", estimated xpsr sd = ",
round(stats::sd(r), 5)
)
message(
"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."
)
message(" - ", msg1)
message(" - ", msg2)
}


# Helper
check_oos <- function(oos, data) {
checkmate::assert_logical(oos, len = 1)
Expand Down
3 changes: 2 additions & 1 deletion R/PKModel.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -112,8 +112,9 @@ PKModel <- R6::R6Class("PKModel",
CL <- theta[2]
V2 <- theta[3]
checkmate::assert_number(CL, lower = 0)
checkmate::assert_number(V2, lower = 0)
checkmate::assert_number(dose, lower = 0)
log(dose) - log(CL * V2)
log(dose) - log(CL) - log(V2)
},

#' @description Simulate data with many subjects
Expand Down
10 changes: 2 additions & 8 deletions R/stan-data.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,12 +184,6 @@ create_stan_data_pk <- function(data, model) {
}
t_obs_pk <- pk_obs[, 1:2, drop = FALSE]
conc_pk <- pk_obs[, 3:4, drop = FALSE]
if (max(conc_pk) > 1e5) {
stop(
"Rescale concentration measurements to larger units so that maximum",
"concentration has smaller numeric value"
)
}
pk_lloq <- as.numeric(pk_obs[, 5])

# Return
Expand Down Expand Up @@ -343,11 +337,11 @@ time_since_last_dose <- function(t_obs_pk, last_two_times) {
n_before_post <- sum(last_two_times <= t_obs_pk[2])
if (n_before_pre == 0) {
n_before_pre <- 1
stop("invalid data")
stop("invalid data, , 0 doses before pre-dose time")
}
if (n_before_post == 0) {
n_before_pre <- 1
stop("invalid data")
stop("invalid data, 0 doses before post-dose time")
}
t_before_pre <- last_two_times[n_before_pre]
t_before_post <- last_two_times[n_before_post]
Expand Down
3 changes: 3 additions & 0 deletions R/stan.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,9 @@ fit_stan <- function(model, data,
time = stan_fit$time()
)
fit <- MultistateModelFit$new(data, sd, model, draws, info)
if (fit$model$has_pk()) {
msmfit_check_xpsr_norm(fit)
}
if (!return_stanfit) {
return(fit)
}
Expand Down
22 changes: 22 additions & 0 deletions R/utils.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -175,3 +175,25 @@ prefit_checks <- function(model, data) {
}
TRUE
}


# Normalize columns of matrix A
normalize_columns <- function(A) {
K <- ncol(A)
cn <- colnames(A)
for (j in seq_len(K)) {
A[, j] <- normalize(A[, j], cn[j])
}
A
}

# Normalize
normalize <- function(a, name) {
sdd <- stats::sd(a)
if (sdd == 0) {
stop("error in normalization of ", name, ", zero variance")
}
x_norm <- (a - mean(a)) / sdd
check_normalized_covariate(x_norm, name)
x_norm
}
8 changes: 4 additions & 4 deletions inst/stan/msm.stan
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ functions {

// log SS AUC
vector log_ss_area_under_conc(vector dose_ss, matrix log_theta_pk){
return(log(dose_ss) - log_theta_pk[:,2] - log_theta_pk[:,3]) ; // D/(CL*V2)
return(log(dose_ss) - log_theta_pk[:,2] - log_theta_pk[:,3]) ; // log D/(CL*V2)
}

// log of hazard multiplier
Expand Down Expand Up @@ -304,7 +304,7 @@ data {

// PK options
int<lower=0,upper=1> I_xpsr;
real<lower=0> xpsr_loc;
real xpsr_loc;
real<lower=0> xpsr_scale;
int<lower=0> nc_ka; // num of predictors for ka
int<lower=0> nc_CL; // num of predictors for CL
Expand Down Expand Up @@ -448,7 +448,7 @@ model {
if(do_haz == 1){
if(nc_haz > 0){
for(k in 1:nc_haz){
beta_oth[1, k] ~ normal(0, 1);
beta_oth[1, k] ~ normal(0, 2);
}
}
if(I_xpsr==1){
Expand All @@ -471,7 +471,7 @@ model {
for(n in 1:size(log_z_pk[1])){
log_z_pk[1, n] ~ normal(0, 1);
}
log_mu_pk[1] ~ normal(0, 3);
log_mu_pk[1] ~ normal(0, 2);
log_sig_pk[1] ~ normal(0, 1);
beta_ka[1] ~ normal(0, 1);
beta_CL[1] ~ normal(0, 1);
Expand Down
5 changes: 1 addition & 4 deletions man/PSSDosingData.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

14 changes: 14 additions & 0 deletions man/msmfit_check_xpsr_norm.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

21 changes: 21 additions & 0 deletions man/msmfit_exposure_df.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading