From 2d740a75cf41842a08bf274780ed070358de8166 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 19 Jun 2026 18:31:47 +0800 Subject: [PATCH 01/19] first stab, next stan module conversion --- R/Prior.R | 136 +++++++++++++++++++--- tests/testthat/test-Prior.R | 25 ++++ tests/testthat/test-SurvivalExponential.R | 16 +++ 3 files changed, 160 insertions(+), 17 deletions(-) diff --git a/R/Prior.R b/R/Prior.R index d7e5c997a..b66c11b60 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -50,7 +50,8 @@ NULL "centre" = "numeric", "validation" = "list", "sample" = "function", - "limits" = "numeric" + "limits" = "numeric", + ".allow_vectors" = "logical" ) ) @@ -72,6 +73,8 @@ NULL #' a function to sample from the prior distribution. #' @typed limits: numeric #' the lower and upper limits for a truncated distribution +#' @typed .allow_vectors: flag +#' whether to allow vector parameters. #' @rdname Prior-class Prior <- function( parameters, @@ -81,7 +84,8 @@ Prior <- function( centre, validation, sample, - limits = c(-Inf, Inf) + limits = c(-Inf, Inf), + .allow_vectors = FALSE ) { .Prior( parameters = parameters, @@ -91,7 +95,8 @@ Prior <- function( display = display, validation = validation, sample = sample, - limits = limits + limits = limits, + .allow_vectors = .allow_vectors ) } @@ -106,7 +111,10 @@ setValidity( param )) } - if (length(object@parameters[[param]]) != 1) { + if ( + !object@.allow_vectors && + length(object@parameters[[param]]) != 1 + ) { return(sprintf("Parameter `%s` must be a single value", param)) } if (!object@validation[[param]](object@parameters[[param]])) { @@ -118,6 +126,18 @@ setValidity( return(return_message) } } + if (object@.allow_vectors) { + # Check that all parameters have either same length or are + # length 1. + all_lengths <- sapply( + object@parameters, + length + ) + none_one_lengths <- setdiff(unique(all_lengths), 1) + if (length(none_one_lengths) > 1) { + return("All parameters must be the same length or length 1") + } + } if (length(object@limits) != 2) { return("Limits must be a vector of length 2") } @@ -266,21 +286,63 @@ NULL #' @describeIn Prior-Getter-Methods The prior's initial value #' @export initialValues.Prior <- function(object, ...) { - samples <- getOption("jmpost.prior_shrinkage") * - object@centre + - (1 - getOption("jmpost.prior_shrinkage")) * object@sample(100) + n_samples <- 100 + centre_value <- object@centre + sample_values <- object@sample(n_samples) + + is_scalar <- length(centre_value) == 1 + if (is_scalar) { + assert_that( + length(sample_values) == n_samples, + msg = "Sample function must return a vector of length n_samples" + ) - valid_samples <- samples[ - samples >= min(object@limits) & samples <= max(object@limits) - ] - assert_that( - length(valid_samples) >= 1, - msg = "Unable to generate an initial value that meets the required constraints" - ) - if (length(valid_samples) == 1) { - return(valid_samples) + samples <- getOption("jmpost.prior_shrinkage") * + object@centre + + (1 - getOption("jmpost.prior_shrinkage")) * object@sample(n_samples) + + valid_samples <- samples[ + samples >= min(object@limits) & samples <= max(object@limits) + ] + assert_that( + length(valid_samples) >= 1, + msg = "Unable to generate an initial value that meets the required constraints" + ) + if (length(valid_samples) == 1) { + return(valid_samples) + } + return(sample(valid_samples, 1)) + } else { + n_centre_vals <- length(centre_value) + assert_that( + ncol(sample_values) == n_centre_vals, + msg = "Sample function must return a matrix with n_samples rows and the same number of columns as the length of the centre value" + ) + samples <- getOption("jmpost.prior_shrinkage") * + matrix( + centre_value, + nrow = n_samples, + ncol = n_centre_vals, + byrow = TRUE + ) + + (1 - getOption("jmpost.prior_shrinkage")) * sample_values + valid_samples <- samples[ + apply(samples, 1, function(row) { + all(row >= min(object@limits) & row <= max(object@limits)) + }), + , + drop = FALSE + ] + assert_that( + nrow(valid_samples) >= 1, + msg = "Unable to generate an initial value that meets the required constraints" + ) + if (nrow(valid_samples) == 1) { + return(valid_samples[1, ]) + } else { + return(valid_samples[sample(nrow(valid_samples), 1), ]) + } } - return(sample(valid_samples, 1)) } @@ -313,6 +375,41 @@ prior_normal <- function(mu, sigma) { } +#' Normal Prior for a Vector Distribution +#' +#' @typed mus: numeric +#' means. +#' @typed sigmas: numeric +#' standard deviations. +#' @family Prior +#' @export +prior_normal_vector <- function(mus, sigmas) { + Prior( + parameters = list( + mus = mus, + sigmas = sigmas, + dim_mus = length(mus), + dim_sigmas = length(sigmas) + ), + display = "normal(mus = [{toString(mus)}], sigmas = [{toString(sigmas)}])", + repr_model = "{name} ~ normal(prior_mus_{name}, prior_sigmas_{name})", + repr_data = c( + "vector[{dim_mus}] prior_mus_{name};", + "vector[{dim_sigmas}] prior_sigmas_{name};" + ), + centre = mus, + sample = \(n) local_rnorm_vector(n, mus, sigmas), + validation = list( + mus = \(x) all(is.numeric(x)), + sigmas = \(x) all(x > 0), + dim_mus = is.count, + dim_sigmas = is.count + ), + .allow_vectors = TRUE + ) +} + + #' Standard Normal Prior Distribution #' #' @@ -719,6 +816,11 @@ NULL #' @rdname Local_Sample local_rnorm <- \(...) rnorm(...) +#' @rdname Local_Sample +local_rnorm_vector <- \(n, mus, sigmas) { + mapply(local_rnorm, n = n, mean = mus, sd = sigmas) +} + #' @rdname Local_Sample local_rcauchy <- \(...) rcauchy(...) diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index bff18597c..1edd8bf21 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -270,3 +270,28 @@ test_that("Parameters in priors must be length 1 #422", { "Parameter `alpha`" ) }) + +test_that("prior_normal_vector works as expected", { + x <- prior_normal_vector(c(1, 2, 3), c(4, 5, 6)) + x_char <- as.character(x) + + expect_equal( + x_char, + "normal(mus = [1, 2, 3], sigmas = [4, 5, 6])" + ) + + x_inits <- initialValues(x) + expect_numeric(x_inits, len = 3) + + expect_equal( + as.StanModule(x, name = "bob"), + StanModule(test_path("models", "Prior_4.stan")) + ) + expect_equal( + as_stan_list(x, name = "bob"), + list( + prior_mu_bob = c(1, 2, 3), + prior_sigma_bob = c(4, 5, 6) + ) + ) +}) diff --git a/tests/testthat/test-SurvivalExponential.R b/tests/testthat/test-SurvivalExponential.R index b4ed93084..c08d439a5 100644 --- a/tests/testthat/test-SurvivalExponential.R +++ b/tests/testthat/test-SurvivalExponential.R @@ -94,3 +94,19 @@ test_that("Print method for SurvivalExponential works as expected", { print(x) }) }) + +test_that("Different priors for the beta components are possible", { + # Same iid prior for all beta components: + expect_snapshot({ + x <- SurvivalExponential(beta = prior_normal(0, 1)) + print(x) + }) + + # Different priors for each beta component: + expect_snapshot({ + x <- SurvivalExponential( + beta = prior_normal_vector(c(0, 1, 2), c(1, 2, 3)) + ) + print(x) + }) +}) From aa3abbb561fd611f2c73f4babb8dcd7b0d8f75b4 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 19 Jun 2026 20:34:41 +0800 Subject: [PATCH 02/19] it works already now! --- R/Prior.R | 6 ++- tests/testthat/test-Prior.R | 9 ++-- tests/testthat/test-SurvivalExponential.R | 57 +++++++++++++++++++++++ 3 files changed, 67 insertions(+), 5 deletions(-) diff --git a/R/Prior.R b/R/Prior.R index b66c11b60..c4a305117 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -394,8 +394,10 @@ prior_normal_vector <- function(mus, sigmas) { display = "normal(mus = [{toString(mus)}], sigmas = [{toString(sigmas)}])", repr_model = "{name} ~ normal(prior_mus_{name}, prior_sigmas_{name})", repr_data = c( - "vector[{dim_mus}] prior_mus_{name};", - "vector[{dim_sigmas}] prior_sigmas_{name};" + "int prior_dim_mus_{name};", + "int prior_dim_sigmas_{name};", + "vector[prior_dim_mus_{name}] prior_mus_{name};", + "vector[prior_dim_sigmas_{name}] prior_sigmas_{name};" ), centre = mus, sample = \(n) local_rnorm_vector(n, mus, sigmas), diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 1edd8bf21..382ea8d46 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -283,15 +283,18 @@ test_that("prior_normal_vector works as expected", { x_inits <- initialValues(x) expect_numeric(x_inits, len = 3) + x_stan_module <- as.StanModule(x, name = "bob") expect_equal( - as.StanModule(x, name = "bob"), + x_stan_module, StanModule(test_path("models", "Prior_4.stan")) ) expect_equal( as_stan_list(x, name = "bob"), list( - prior_mu_bob = c(1, 2, 3), - prior_sigma_bob = c(4, 5, 6) + prior_mus_bob = c(1, 2, 3), + prior_sigmas_bob = c(4, 5, 6), + prior_dim_mus_bob = 3, + prior_dim_sigmas_bob = 3 ) ) }) diff --git a/tests/testthat/test-SurvivalExponential.R b/tests/testthat/test-SurvivalExponential.R index c08d439a5..5592cb45b 100644 --- a/tests/testthat/test-SurvivalExponential.R +++ b/tests/testthat/test-SurvivalExponential.R @@ -109,4 +109,61 @@ test_that("Different priors for the beta components are possible", { ) print(x) }) + + skip_if_not(is_full_test()) + + true_lambda <- 1 / 100 + true_beta <- c(0.5, -0.2, 0.1) + set.seed(2034) + jdat <- SimJointData( + design = list(SimGroup(700, "Arm-A", "Study-X")), + survival = SimSurvivalExponential( + lambda = true_lambda, + lambda_censor = 1 / 9000, + beta_cat = c("A" = 0, "B" = true_beta[1], "C" = true_beta[2]), + beta_cont = true_beta[3], + ), + longitudinal = SimLongitudinalRandomSlope( + slope_mu = 0, + slope_sigma = 0.5 + ) + ) + + dat_os <- jdat@survival + + jm <- JointModel( + survival = SurvivalExponential( + beta = prior_normal_vector(c(0, 1, 2), c(1, 2, 3)) + ) + ) + + jdat <- DataJoint( + subject = DataSubject( + data = dat_os, + subject = "subject", + arm = "arm", + study = "study" + ), + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont + ) + ) + + mp <- run_quietly({ + sampleStanModel( + jm, + data = jdat, + iter_sampling = 600, + iter_warmup = 500, + chains = 1, + refresh = 0, + parallel_chains = 1, + seed = 123 + ) + }) + + # Variables to extract (order important) + vars <- c("sm_exp_lambda", "beta_os_cov") + results_summary <- cmdstanr::as.CmdStanMCMC(mp)$summary(vars) }) From 08e7c93bc8c1ddeb973ca072f6de563347971818 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 19 Jun 2026 20:35:25 +0800 Subject: [PATCH 03/19] run docs --- NAMESPACE | 1 + man/Local_Sample.Rd | 3 +++ man/Prior-class.Rd | 5 ++++- man/prior_beta.Rd | 1 + man/prior_cauchy.Rd | 1 + man/prior_gamma.Rd | 1 + man/prior_init_only.Rd | 1 + man/prior_invgamma.Rd | 1 + man/prior_logistic.Rd | 1 + man/prior_loglogistic.Rd | 1 + man/prior_lognormal.Rd | 1 + man/prior_normal.Rd | 1 + man/prior_normal_vector.Rd | 32 ++++++++++++++++++++++++++++++++ man/prior_std_normal.Rd | 1 + man/prior_student_t.Rd | 1 + man/prior_uniform.Rd | 1 + 16 files changed, 52 insertions(+), 1 deletion(-) create mode 100644 man/prior_normal_vector.Rd diff --git a/NAMESPACE b/NAMESPACE index e79038600..0651c9853 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -272,6 +272,7 @@ export(prior_logistic) export(prior_loglogistic) export(prior_lognormal) export(prior_normal) +export(prior_normal_vector) export(prior_std_normal) export(prior_student_t) export(prior_uniform) diff --git a/man/Local_Sample.Rd b/man/Local_Sample.Rd index 000958413..6bbe1459c 100644 --- a/man/Local_Sample.Rd +++ b/man/Local_Sample.Rd @@ -3,6 +3,7 @@ \name{Local_Sample} \alias{Local_Sample} \alias{local_rnorm} +\alias{local_rnorm_vector} \alias{local_rcauchy} \alias{local_rgamma} \alias{local_rlnorm} @@ -16,6 +17,8 @@ \usage{ local_rnorm(...) +local_rnorm_vector(n, mus, sigmas) + local_rcauchy(...) local_rgamma(...) diff --git a/man/Prior-class.Rd b/man/Prior-class.Rd index 985c35f4f..6604e4e2f 100644 --- a/man/Prior-class.Rd +++ b/man/Prior-class.Rd @@ -15,7 +15,8 @@ Prior( centre, validation, sample, - limits = c(-Inf, Inf) + limits = c(-Inf, Inf), + .allow_vectors = FALSE ) } \arguments{ @@ -35,6 +36,8 @@ the same names as the \code{paramaters} slot.} \item{sample}{(\code{function}) a function to sample from the prior distribution.} \item{limits}{(\code{numeric}) the lower and upper limits for a truncated distribution} + +\item{.allow_vectors}{(\code{flag}) whether to allow vector parameters.} } \description{ Specifies the prior distribution in a Stan Model diff --git a/man/prior_beta.Rd b/man/prior_beta.Rd index b4d068592..1479a2e41 100644 --- a/man/prior_beta.Rd +++ b/man/prior_beta.Rd @@ -24,6 +24,7 @@ Other Prior: \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_student_t]{prior_student_t()}}, \code{\link[=prior_uniform]{prior_uniform()}} diff --git a/man/prior_cauchy.Rd b/man/prior_cauchy.Rd index 271f6e5d3..c559b096c 100644 --- a/man/prior_cauchy.Rd +++ b/man/prior_cauchy.Rd @@ -24,6 +24,7 @@ Other Prior: \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_student_t]{prior_student_t()}}, \code{\link[=prior_uniform]{prior_uniform()}} diff --git a/man/prior_gamma.Rd b/man/prior_gamma.Rd index cc01229a4..cc386c882 100644 --- a/man/prior_gamma.Rd +++ b/man/prior_gamma.Rd @@ -24,6 +24,7 @@ Other Prior: \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_student_t]{prior_student_t()}}, \code{\link[=prior_uniform]{prior_uniform()}} diff --git a/man/prior_init_only.Rd b/man/prior_init_only.Rd index 2bc1df2f7..b36e9153d 100644 --- a/man/prior_init_only.Rd +++ b/man/prior_init_only.Rd @@ -24,6 +24,7 @@ Other Prior: \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_student_t]{prior_student_t()}}, \code{\link[=prior_uniform]{prior_uniform()}} diff --git a/man/prior_invgamma.Rd b/man/prior_invgamma.Rd index a6d41b86a..c3a7f3d04 100644 --- a/man/prior_invgamma.Rd +++ b/man/prior_invgamma.Rd @@ -24,6 +24,7 @@ Other Prior: \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_student_t]{prior_student_t()}}, \code{\link[=prior_uniform]{prior_uniform()}} diff --git a/man/prior_logistic.Rd b/man/prior_logistic.Rd index 4148facda..fc0447c14 100644 --- a/man/prior_logistic.Rd +++ b/man/prior_logistic.Rd @@ -24,6 +24,7 @@ Other Prior: \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_student_t]{prior_student_t()}}, \code{\link[=prior_uniform]{prior_uniform()}} diff --git a/man/prior_loglogistic.Rd b/man/prior_loglogistic.Rd index 5bff32f68..62a46e5d5 100644 --- a/man/prior_loglogistic.Rd +++ b/man/prior_loglogistic.Rd @@ -24,6 +24,7 @@ Other Prior: \code{\link[=prior_logistic]{prior_logistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_student_t]{prior_student_t()}}, \code{\link[=prior_uniform]{prior_uniform()}} diff --git a/man/prior_lognormal.Rd b/man/prior_lognormal.Rd index 372ada3ad..863357980 100644 --- a/man/prior_lognormal.Rd +++ b/man/prior_lognormal.Rd @@ -24,6 +24,7 @@ Other Prior: \code{\link[=prior_logistic]{prior_logistic()}}, \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_student_t]{prior_student_t()}}, \code{\link[=prior_uniform]{prior_uniform()}} diff --git a/man/prior_normal.Rd b/man/prior_normal.Rd index 319632732..2f532e8c2 100644 --- a/man/prior_normal.Rd +++ b/man/prior_normal.Rd @@ -24,6 +24,7 @@ Other Prior: \code{\link[=prior_logistic]{prior_logistic()}}, \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_student_t]{prior_student_t()}}, \code{\link[=prior_uniform]{prior_uniform()}} diff --git a/man/prior_normal_vector.Rd b/man/prior_normal_vector.Rd new file mode 100644 index 000000000..5d7fbb7fa --- /dev/null +++ b/man/prior_normal_vector.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Prior.R +\name{prior_normal_vector} +\alias{prior_normal_vector} +\title{Normal Prior for a Vector Distribution} +\usage{ +prior_normal_vector(mus, sigmas) +} +\arguments{ +\item{mus}{(\code{numeric}) means.} + +\item{sigmas}{(\code{numeric}) standard deviations.} +} +\description{ +Normal Prior for a Vector Distribution +} +\seealso{ +Other Prior: +\code{\link[=prior_beta]{prior_beta()}}, +\code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_init_only]{prior_init_only()}}, +\code{\link[=prior_invgamma]{prior_invgamma()}}, +\code{\link[=prior_logistic]{prior_logistic()}}, +\code{\link[=prior_loglogistic]{prior_loglogistic()}}, +\code{\link[=prior_lognormal]{prior_lognormal()}}, +\code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_std_normal]{prior_std_normal()}}, +\code{\link[=prior_student_t]{prior_student_t()}}, +\code{\link[=prior_uniform]{prior_uniform()}} +} +\concept{Prior} diff --git a/man/prior_std_normal.Rd b/man/prior_std_normal.Rd index f0c4de840..46f67d4f2 100644 --- a/man/prior_std_normal.Rd +++ b/man/prior_std_normal.Rd @@ -20,6 +20,7 @@ Other Prior: \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_student_t]{prior_student_t()}}, \code{\link[=prior_uniform]{prior_uniform()}} } diff --git a/man/prior_student_t.Rd b/man/prior_student_t.Rd index 119f96b62..1101043e8 100644 --- a/man/prior_student_t.Rd +++ b/man/prior_student_t.Rd @@ -27,6 +27,7 @@ Other Prior: \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_uniform]{prior_uniform()}} } diff --git a/man/prior_uniform.Rd b/man/prior_uniform.Rd index 7c715e48c..15efa1719 100644 --- a/man/prior_uniform.Rd +++ b/man/prior_uniform.Rd @@ -25,6 +25,7 @@ Other Prior: \code{\link[=prior_loglogistic]{prior_loglogistic()}}, \code{\link[=prior_lognormal]{prior_lognormal()}}, \code{\link[=prior_normal]{prior_normal()}}, +\code{\link[=prior_normal_vector]{prior_normal_vector()}}, \code{\link[=prior_std_normal]{prior_std_normal()}}, \code{\link[=prior_student_t]{prior_student_t()}} } From b39f3ef5e943d3e12e24b2ddd5b57d23bbfcf282 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 19 Jun 2026 20:39:44 +0800 Subject: [PATCH 04/19] explain in vignette --- vignettes/quickstart.Rmd | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/vignettes/quickstart.Rmd b/vignettes/quickstart.Rmd index a9fa408ea..c1e88b8e9 100644 --- a/vignettes/quickstart.Rmd +++ b/vignettes/quickstart.Rmd @@ -98,6 +98,16 @@ random_slope_model <- LongitudinalRandomSlope( This sets the prior for the `intercept` to be a $N(40, 5)$ distribution and the prior for the `slope_mu` parameter to be a $N(10, 2)$ distribution. +For the survival models, we can set different independent normal priors for the coefficients of the covariates in the linear predictor using the `prior_normal_vector()` specification. Here we can give a vector of the same length as the number of design matrix columns for both the means or standard deviations of the normal distributions, or we can give a single value that will be repeated for all covariates. +For example, say we have 3 design matrix columns, i.e. coefficients, in the linear predictor of the survival model, then we can set the priors for these coefficients as follows: + +```{r survival_prior_inits} +survival_model <- SurvivalWeibullPH( + beta = prior_normal_vector(mus = c(-1, 0, 5), sigmas = c(1, 1, 10)) +) +survival_model +``` + ## Separate Models It is also possible to not link the longitudinal and the survival models, by using From 1dcbe881619116f757ed54ac10022179293db5e0 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 19 Jun 2026 20:49:02 +0800 Subject: [PATCH 05/19] add NEWS entry --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 9cb74c727..c3436de1d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,7 @@ - Included new `LongitudinalRandomEffects()` function which can be used to extract the patient-level random effects parameter samples from a `JointModelSample` object (#423). - Introduced the `saveObject()` method for `JointModelSample` objects in order to serialise them to disk (#431). - Added support for truncated prior distributions e.g. you can now apply a normal prior to a strictly positive parameter and jmpost will take care of adjusting the density accordingly (#429). +- Added `prior_normal_vector()` which can be used to assign different normal prior distributions to the coefficients of the survival model covariates. - Included new Gamma distribution survival model (#411). - Reworked LOO calculations to apply to each individual submodel and disabled LOO calculations for the overall joint model (#402). - Added support for additive variance (#403). From 3c59accf984003e98400d843cfa02b01d8571172 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 19 Jun 2026 21:12:02 +0800 Subject: [PATCH 06/19] comply with new assertion for sampler --- tests/testthat/test-Parameter.R | 2 +- tests/testthat/test-ParameterList.R | 4 +-- tests/testthat/test-Prior.R | 6 ++-- tests/testthat/test-initialValues.R | 54 +++++++++++++++++++++++++---- 4 files changed, 53 insertions(+), 13 deletions(-) diff --git a/tests/testthat/test-Parameter.R b/tests/testthat/test-Parameter.R index b81bb9262..58d43112b 100644 --- a/tests/testthat/test-Parameter.R +++ b/tests/testthat/test-Parameter.R @@ -7,7 +7,7 @@ test_that("Parameters smoke tests", { initialValues(p), expected_mu * 0.5 ), - local_rbeta = \(...) 0 + local_rbeta = \(n, ...) rep(0, n) ) expect_equal(names(p), "intercept") }) diff --git a/tests/testthat/test-ParameterList.R b/tests/testthat/test-ParameterList.R index 192e207ce..de39ae9cf 100644 --- a/tests/testthat/test-ParameterList.R +++ b/tests/testthat/test-ParameterList.R @@ -17,8 +17,8 @@ test_that("ParameterList smoke tests", { )) expect_equal(actual, expected) }, - local_rgamma = \(...) 0, - local_rnorm = \(...) 0 + local_rgamma = \(n, ...) rep(0, n), + local_rnorm = \(n, ...) rep(0, n) ) # Can render to character diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 382ea8d46..28a4d6e84 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -5,7 +5,7 @@ test_that("Priors work as expected", { initialValues(x), 4 * 0.5 ), - local_rnorm = \(...) 0 + local_rnorm = \(n, ...) rep(0, n) ) expect_equal( as.StanModule(x, name = "bob"), @@ -22,7 +22,7 @@ test_that("Priors work as expected", { initialValues(x), exp(log(4) + 2) * 0.5 ), - local_rlnorm = \(...) 0 + local_rlnorm = \(n, ...) rep(0, n) ) expect_equal( as.StanModule(x, name = "tim"), @@ -181,7 +181,7 @@ test_that("jmpost.prior_shrinkage works as expected", { ## Reset Shrinkage factor options("jmpost.prior_shrinkage" = 0.5) }, - local_rnorm = \(...) 4 + local_rnorm = \(n, ...) rep(4, n) ) }) diff --git a/tests/testthat/test-initialValues.R b/tests/testthat/test-initialValues.R index b2f0a88bc..333f0a102 100644 --- a/tests/testthat/test-initialValues.R +++ b/tests/testthat/test-initialValues.R @@ -42,16 +42,56 @@ test_that("initialValues() works as expected", { # show that if we mock the random number generator, we get the same initial values ivs <- testthat::with_mocked_bindings( initialValues(jm, n_chains = 2), - local_rnorm = \(...) 0, - local_rbeta = \(...) 0, - local_rlnorm = \(...) 0, - local_rgamma = \(...) 0, - local_runif = \(...) 0, - local_rlogis = \(...) 0, + local_rnorm = \(n, ...) rep(0, n), + local_rbeta = \(n, ...) rep(0, n), + local_rlnorm = \(n, ...) rep(0, n), + local_rgamma = \(n, ...) rep(0, n), + local_runif = \(n, ...) rep(0, n), + local_rlogis = \(n, ...) rep(0, n) ) expect_equal(ivs[[1]], ivs[[2]]) }) +test_that("initialValues() works also for vectorized parameters", { + jm <- JointModel( + longitudinal = LongitudinalRandomSlope(), + survival = SurvivalWeibullPH( + beta = prior_normal_vector(mus = c(-1, 0, 5), sigmas = c(1, 1, 10)) + ), + link = linkDSLD() + ) + + set.seed(341) + initial_values <- initialValues(jm, n_chains = 2) + + # Ensure that we actually got 2 chains worth of initial values + expect_length(initial_values, 2) + + # Ensure that we get different initial values per chain + expect_numeric(initial_values[[1]]$beta_os_cov, len = 3) + expect_numeric(initial_values[[2]]$beta_os_cov, len = 3) + expect_true( + all(initial_values[[1]]$beta_os_cov != initial_values[[2]]$beta_os_cov) + ) + + # Ensure each inner list has the same parameters + expect_equal( + names(initial_values[[1]]), + names(initial_values[[2]]) + ) + + # Show that if we mock the random number generator, we get the same initial values: + ivs <- testthat::with_mocked_bindings( + initialValues(jm, n_chains = 2), + local_rnorm = \(n, ...) rep(0, n), + local_rlnorm = \(n, ...) rep(0, n), + local_rgamma = \(n, ...) rep(0, n), + local_rnorm_vector = \(n, mus, sigmas) { + matrix(0, nrow = n, ncol = length(mus)) + } + ) + expect_equal(ivs[[1]], ivs[[2]]) +}) test_that("ensure_initial_values() works as expected", { pars <- ParameterList( @@ -86,7 +126,7 @@ test_that("ensure_initial_values() works as expected", { }) -test_that("intial values for fixed distributions gives valid values", { +test_that("initial values for fixed distributions gives valid values", { set.seed(3150) gsfmodel <- LongitudinalGSF(centred = TRUE) ivs <- initialValues(gsfmodel, n_chains = 100) From 5a9d94f045a99e34490895579f57cec5e6ee86dc Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 19 Jun 2026 21:12:47 +0800 Subject: [PATCH 07/19] add prior_normal_vector to reference section --- _pkgdown.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index be237c1f3..0c46ed614 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -64,6 +64,7 @@ reference: - prior_loglogistic - prior_lognormal - prior_normal + - prior_normal_vector - prior_uniform - prior_std_normal - prior_student_t From dab1dabcaf297698b0ce0f70fe24db3c4a3e714e Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 19 Jun 2026 21:13:41 +0800 Subject: [PATCH 08/19] fix linter --- R/Prior.R | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/R/Prior.R b/R/Prior.R index c4a305117..07e0df3e9 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -316,7 +316,10 @@ initialValues.Prior <- function(object, ...) { n_centre_vals <- length(centre_value) assert_that( ncol(sample_values) == n_centre_vals, - msg = "Sample function must return a matrix with n_samples rows and the same number of columns as the length of the centre value" + msg = paste( + "Sample function must return a matrix with n_samples rows and the same", + "number of columns as the length of the centre value" + ) ) samples <- getOption("jmpost.prior_shrinkage") * matrix( From 35ca5cc650e5d2e4fb0bf7a652c7827aa568a04b Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Sat, 20 Jun 2026 20:55:47 +0800 Subject: [PATCH 09/19] fix some tests --- tests/testthat/_snaps/SurvivalExponential.md | 24 ++++++++++++++++++++ tests/testthat/test-misc_models.R | 3 ++- tests/testthat/test-populationHR.R | 2 +- 3 files changed, 27 insertions(+), 2 deletions(-) diff --git a/tests/testthat/_snaps/SurvivalExponential.md b/tests/testthat/_snaps/SurvivalExponential.md index abf6a39d7..844e993b4 100644 --- a/tests/testthat/_snaps/SurvivalExponential.md +++ b/tests/testthat/_snaps/SurvivalExponential.md @@ -22,3 +22,27 @@ beta_os_cov ~ gamma(alpha = 3, beta = 4) +# Different priors for the beta components are possible + + Code + x <- SurvivalExponential(beta = prior_normal(0, 1)) + print(x) + Output + + Exponential Survival Model with parameters: + sm_exp_lambda ~ gamma(alpha = 2, beta = 5) + beta_os_cov ~ normal(mu = 0, sigma = 1) + + +--- + + Code + x <- SurvivalExponential(beta = prior_normal_vector(c(0, 1, 2), c(1, 2, 3))) + print(x) + Output + + Exponential Survival Model with parameters: + sm_exp_lambda ~ gamma(alpha = 2, beta = 5) + beta_os_cov ~ normal(mus = [0, 1, 2], sigmas = [1, 2, 3]) + + diff --git a/tests/testthat/test-misc_models.R b/tests/testthat/test-misc_models.R index 813593309..95a8ca15f 100644 --- a/tests/testthat/test-misc_models.R +++ b/tests/testthat/test-misc_models.R @@ -14,7 +14,8 @@ test_that("Longitudinal Model doesn't print sampler rejection messages", { iter_warmup = 3, chains = 1, refresh = 0, - parallel_chains = 1 + parallel_chains = 1, + seed = 1234 ) }) }) diff --git a/tests/testthat/test-populationHR.R b/tests/testthat/test-populationHR.R index 3f907c533..ab7225b64 100644 --- a/tests/testthat/test-populationHR.R +++ b/tests/testthat/test-populationHR.R @@ -16,7 +16,7 @@ test_that("populationHR works as expected for default parameters", { expect_equal( 0.2 - -0.4, result[[1]][["cov_catC", "mean"]] - result[[1]][["cov_catB", "mean"]], - tolerance = 0.2 + tolerance = 0.3 ) expect_equal( From 1f5e411e98460c2eb037a71409e5ef6b2a5940cd Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Sun, 21 Jun 2026 21:41:17 +0800 Subject: [PATCH 10/19] add .stan model --- .gitignore | 1 + tests/testthat/models/Prior_4.stan | 10 ++++++++++ 2 files changed, 11 insertions(+) create mode 100644 tests/testthat/models/Prior_4.stan diff --git a/.gitignore b/.gitignore index ef1095472..bb48685d8 100755 --- a/.gitignore +++ b/.gitignore @@ -76,6 +76,7 @@ local/ # Compiled stan model inst/stanmodels/stan_model **/models/* +!**/models/*.stan diff --git a/tests/testthat/models/Prior_4.stan b/tests/testthat/models/Prior_4.stan new file mode 100644 index 000000000..6490d40e9 --- /dev/null +++ b/tests/testthat/models/Prior_4.stan @@ -0,0 +1,10 @@ +data { + int prior_dim_mus_bob; + int prior_dim_sigmas_bob; + vector[prior_dim_mus_bob] prior_mus_bob; + vector[prior_dim_sigmas_bob] prior_sigmas_bob; +} + +model { + bob ~ normal(prior_mus_bob, prior_sigmas_bob); +} From 3afed61561853dc749b9a5eb010285f0be1bb6b4 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Mon, 22 Jun 2026 08:14:55 +0800 Subject: [PATCH 11/19] introduce lower bound settings --- NEWS.md | 2 +- R/settings.R | 15 ++++++++++++++- R/utilities.R | 4 ++-- man/jmpost-settings.Rd | 12 ++++++++++++ tests/testthat/test-jinjar_render.R | 2 +- tests/testthat/test-misc_models.R | 4 ++-- 6 files changed, 32 insertions(+), 7 deletions(-) diff --git a/NEWS.md b/NEWS.md index c3436de1d..74bd69161 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,3 @@ - # jmpost (development version) - `LongitudinalModel` objects now have a new slot `scaled_variance` which stores the used variance option, and this is part of the `print` output now thereby transparently communicating to the user the choice of the multiplicative or additive error model. @@ -13,6 +12,7 @@ - Added support for additive variance (#403). - Added support for independent variances per study/arm (#389). - Miscellaneous bug fixes. +- Introduce new package options for bounds close to zero which are used to avoid MCMC sampler starting warnings. # jmpost 0.0.1 diff --git a/R/settings.R b/R/settings.R index c0cc820c1..555c21a89 100644 --- a/R/settings.R +++ b/R/settings.R @@ -42,6 +42,17 @@ #' higher values of `n` lead to better accuracy of the approximation but at the cost of #' increased computational time. #' +#' ## `jmpost.double_eps` and `jmpost.double_neg_eps` +#' +#' Default = `sqrt(.Machine$double.eps)` and `sqrt(.Machine$double.neg.eps)` +#' +#' These are the smallest positive floating-point numbers `x` which are used +#' in the parameter definitions as lower boundaries to prevent numerical issues +#' with the MCMC sampler in the beginning of the sampling process. +#' The default values are set to the square root of the machine epsilon for +#' double precision numbers, which is a common choice to balance numerical stability +#' and the ability to explore the parameter space effectively. +#' #' @examples #' \dontrun{ #' options(jmpost.prior_shrinkage = 0.5) @@ -58,7 +69,9 @@ set_options <- function() { jmpost_opts <- list( jmpost.cache_dir = cache_dir, jmpost.prior_shrinkage = 0.5, - jmpost.gauss_quad_n = 15 + jmpost.gauss_quad_n = 15, + jmpost.double_eps = sqrt(.Machine$double.eps), + jmpost.double_neg_eps = sqrt(.Machine$double.neg.eps) ) for (opt in names(jmpost_opts)) { if (!opt %in% current_opts) { diff --git a/R/utilities.R b/R/utilities.R index ae92df3cb..b5a238c56 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -207,8 +207,8 @@ samples_median_ci <- function(samples, level = 0.95) { decorated_render <- function(...) { jinjar::render( ..., - machine_double_eps = 0, - machine_double_neg_eps = 0 + machine_double_eps = getOption("jmpost.double_eps"), + machine_double_neg_eps = getOption("jmpost.double_neg_eps"), ) } diff --git a/man/jmpost-settings.Rd b/man/jmpost-settings.Rd index c847e7f4a..5557ab6d3 100644 --- a/man/jmpost-settings.Rd +++ b/man/jmpost-settings.Rd @@ -49,6 +49,18 @@ This option specifies the \code{n} argument in the call to \code{\link[statmod:g higher values of \code{n} lead to better accuracy of the approximation but at the cost of increased computational time. } + +\subsection{\code{jmpost.double_eps} and \code{jmpost.double_neg_eps}}{ + +Default = \code{sqrt(.Machine$double.eps)} and \code{sqrt(.Machine$double.neg.eps)} + +These are the smallest positive floating-point numbers \code{x} which are used +in the parameter definitions as lower boundaries to prevent numerical issues +with the MCMC sampler in the beginning of the sampling process. +The default values are set to the square root of the machine epsilon for +double precision numbers, which is a common choice to balance numerical stability +and the ability to explore the parameter space effectively. +} } \examples{ \dontrun{ diff --git a/tests/testthat/test-jinjar_render.R b/tests/testthat/test-jinjar_render.R index 5157b7463..919877370 100644 --- a/tests/testthat/test-jinjar_render.R +++ b/tests/testthat/test-jinjar_render.R @@ -21,7 +21,7 @@ test_that("Global variables defined by `decorated_render() are usable", { "%d hi there %s and %s", 5, "bob", - 0 + scales::scientific(getOption("jmpost.double_eps"), digits = 5) ) expect_equal(observed, expected) }) diff --git a/tests/testthat/test-misc_models.R b/tests/testthat/test-misc_models.R index 95a8ca15f..8eeddc2f7 100644 --- a/tests/testthat/test-misc_models.R +++ b/tests/testthat/test-misc_models.R @@ -3,7 +3,7 @@ test_data_1 <- ensure_test_data_1() test_that("Longitudinal Model doesn't print sampler rejection messages", { # These rejections typically happen when the sampler samples a # 0 value for the variance parameter. Sensible initial values + - # setting near 0 limits (as opposed to 0) should avoid this + # setting near 0 limits (as opposed to 0) should avoid this. mp <- capture_messages({ devnull_out <- capture.output({ @@ -15,7 +15,7 @@ test_that("Longitudinal Model doesn't print sampler rejection messages", { chains = 1, refresh = 0, parallel_chains = 1, - seed = 1234 + seed = 324 ) }) }) From 45d68eaacd0042ffb476fb033f9fdb2de71ea66d Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Mon, 22 Jun 2026 08:19:01 +0800 Subject: [PATCH 12/19] be more tolerant --- tests/testthat/test-populationHR.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-populationHR.R b/tests/testthat/test-populationHR.R index ab7225b64..c226f6952 100644 --- a/tests/testthat/test-populationHR.R +++ b/tests/testthat/test-populationHR.R @@ -22,7 +22,7 @@ test_that("populationHR works as expected for default parameters", { expect_equal( 0.2, result[[1]][["cov_cont", "mean"]], - tolerance = 0.1 + tolerance = 0.3 ) # Summary calculations are match expectations From 617b2fba3c56c6de3a986f29b794fc8f9a4a146b Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Mon, 22 Jun 2026 08:26:58 +0800 Subject: [PATCH 13/19] remove NOT_CRAN environment variable from docs workflow just to see if it makes a difference --- .github/workflows/docs.yaml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml index 3c4093b49..53887c7c7 100644 --- a/.github/workflows/docs.yaml +++ b/.github/workflows/docs.yaml @@ -48,6 +48,5 @@ jobs: CMDSTAN=/root/.cmdstan CMDSTAN_PATH=/root/.cmdstan CMDSTANR_NO_VER_CHECK=true - NOT_CRAN=TRUE JMPOST_GRAPH_SNAPSHOT=TRUE default-landing-page: 'main' # Set to 'latest-tag' after CRAN release or when you create a tag From 25a7b4f61f2b415f172c286f763b882dd1bc0af5 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Mon, 22 Jun 2026 10:44:27 +0800 Subject: [PATCH 14/19] fix class doc --- R/Prior.R | 1 + man/Prior-class.Rd | 2 ++ 2 files changed, 3 insertions(+) diff --git a/R/Prior.R b/R/Prior.R index 07e0df3e9..b96339846 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -36,6 +36,7 @@ NULL #' @slot display (`string`)\cr See arguments. #' @slot sample (`function`)\cr See arguments. #' @slot limits (`numeric`)\cr See arguments. +#' @slot .allow_vectors (`logical`)\cr See arguments. #' #' @family Prior-internal #' @export Prior diff --git a/man/Prior-class.Rd b/man/Prior-class.Rd index 6604e4e2f..490da7d57 100644 --- a/man/Prior-class.Rd +++ b/man/Prior-class.Rd @@ -60,6 +60,8 @@ Specifies the prior distribution in a Stan Model \item{\code{sample}}{(\code{function})\cr See arguments.} \item{\code{limits}}{(\code{numeric})\cr See arguments.} + +\item{\code{.allow_vectors}}{(\code{logical})\cr See arguments.} }} \seealso{ From 0025efcbcdf31e5eb5d58b500147937a9b62cdaf Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Mon, 22 Jun 2026 10:45:14 +0800 Subject: [PATCH 15/19] remove one more flag from docs.yaml to see if we can trigger staged dep workflow --- .github/workflows/docs.yaml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml index 53887c7c7..3a5807f78 100644 --- a/.github/workflows/docs.yaml +++ b/.github/workflows/docs.yaml @@ -48,5 +48,4 @@ jobs: CMDSTAN=/root/.cmdstan CMDSTAN_PATH=/root/.cmdstan CMDSTANR_NO_VER_CHECK=true - JMPOST_GRAPH_SNAPSHOT=TRUE default-landing-page: 'main' # Set to 'latest-tag' after CRAN release or when you create a tag From 7af8049d217676a09bf70943ec596b58ca2095c1 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Mon, 22 Jun 2026 10:54:46 +0800 Subject: [PATCH 16/19] refactor: update deps-installation-method to use |- for trailing newline consistency --- .github/workflows/check.yaml | 9 ++++++--- .github/workflows/docs.yaml | 3 ++- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/.github/workflows/check.yaml b/.github/workflows/check.yaml index 86e006c22..da35e26a0 100644 --- a/.github/workflows/check.yaml +++ b/.github/workflows/check.yaml @@ -27,7 +27,8 @@ jobs: secrets: REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} with: - deps-installation-method: | + # Use |- so the reusable workflow receives the exact string without a trailing newline. + deps-installation-method: |- staged-dependencies additional-env-vars: | CMDSTAN=/root/.cmdstan @@ -45,7 +46,8 @@ jobs: secrets: REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} with: - deps-installation-method: | + # Use |- so the reusable workflow receives the exact string without a trailing newline. + deps-installation-method: |- staged-dependencies additional-env-vars: | CMDSTAN=/root/.cmdstan @@ -61,7 +63,8 @@ jobs: secrets: REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} with: - deps-installation-method: | + # Use |- so the reusable workflow receives the exact string without a trailing newline. + deps-installation-method: |- staged-dependencies auto-update: true gitleaks: diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml index 3a5807f78..fa4e83165 100644 --- a/.github/workflows/docs.yaml +++ b/.github/workflows/docs.yaml @@ -42,7 +42,8 @@ jobs: secrets: REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} with: - deps-installation-method: | + # Use |- so the reusable workflow receives the exact string without a trailing newline. + deps-installation-method: |- staged-dependencies additional-env-vars: | CMDSTAN=/root/.cmdstan From 026172a12f7b01afe5ec1fe8b2c3095ce431990b Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 24 Jun 2026 18:53:48 +0800 Subject: [PATCH 17/19] add model.matrix() method for DataSurvival objects and explain it --- NAMESPACE | 1 + R/DataSurvival.R | 21 +++++++++++++++++---- man/as_stan_list.DataObject.Rd | 8 ++++++++ tests/testthat/_snaps/DataSurvival.md | 7 +++++++ tests/testthat/test-DataSurvival.R | 18 ++++++++++++++++++ vignettes/.gitignore | 1 + vignettes/quickstart.Rmd | 12 ++++++++++-- 7 files changed, 62 insertions(+), 6 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 0651c9853..34924da3c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -161,6 +161,7 @@ S3method(linkTTG,LongitudinalSteinFojo) S3method(linkTTG,PromiseLongitudinalModel) S3method(linkTTG,default) S3method(median,Prior) +S3method(model.matrix,DataSurvival) S3method(names,LinkComponent) S3method(names,Parameter) S3method(names,ParameterList) diff --git a/R/DataSurvival.R b/R/DataSurvival.R index 7f416b69f..4175e49be 100644 --- a/R/DataSurvival.R +++ b/R/DataSurvival.R @@ -118,14 +118,27 @@ as.data.frame.DataSurvival <- function(x, ...) { #' @rdname as_stan_list.DataObject #' @family DataSurvival #' @export -as_stan_list.DataSurvival <- function(object, ...) { - df <- as.data.frame(object) - vars <- extractVariableNames(object) - +model.matrix.DataSurvival <- function( + object, + df = as.data.frame(object), + vars = extractVariableNames(object), + ... +) { design_mat <- stats::model.matrix(vars$frm, data = df) remove_index <- grep("(Intercept)", colnames(design_mat), fixed = TRUE) design_mat <- design_mat[, -remove_index, drop = FALSE] rownames(design_mat) <- NULL + design_mat +} + +#' @rdname as_stan_list.DataObject +#' @family DataSurvival +#' @export +as_stan_list.DataSurvival <- function(object, ...) { + df <- as.data.frame(object) + vars <- extractVariableNames(object) + + design_mat <- model.matrix(object, df = df, vars = vars) # Parameters for efficient integration of hazard function -> survival function gh_parameters <- statmod::gauss.quad( diff --git a/man/as_stan_list.DataObject.Rd b/man/as_stan_list.DataObject.Rd index 32c9a791c..f88e9987e 100644 --- a/man/as_stan_list.DataObject.Rd +++ b/man/as_stan_list.DataObject.Rd @@ -6,6 +6,7 @@ \alias{as.list.DataSubject} \alias{as_stan_list.DataLongitudinal} \alias{as.list.DataLongitudinal} +\alias{model.matrix.DataSurvival} \alias{as_stan_list.DataSurvival} \alias{as.list.DataSurvival} \alias{as_stan_list.DataObject} @@ -21,6 +22,13 @@ \method{as.list}{DataLongitudinal}(x, ...) +\method{model.matrix}{DataSurvival}( + object, + df = as.data.frame(object), + vars = extractVariableNames(object), + ... +) + \method{as_stan_list}{DataSurvival}(object, ...) \method{as.list}{DataSurvival}(x, ...) diff --git a/tests/testthat/_snaps/DataSurvival.md b/tests/testthat/_snaps/DataSurvival.md index 479e1ec24..7b1f3caa7 100644 --- a/tests/testthat/_snaps/DataSurvival.md +++ b/tests/testthat/_snaps/DataSurvival.md @@ -16,3 +16,10 @@ Formula = Surv(vtime, vevent) ~ vcov1 * vcov2 + vcov1:vcov2 + vcov1^2 + vcov2^2 +# model.matrix for DataSurvival works as expected + + structure(c(0, 0, 1, 1, 0, -0.560475646552213, -0.23017748948328, + 1.55870831414912, 0.070508391424576, 0.129287735160946, 0, 0, + 1.55870831414912, 0.070508391424576, 0), dim = c(5L, 3L), dimnames = list( + NULL, c("vcov1B", "vcov2", "vcov1B:vcov2"))) + diff --git a/tests/testthat/test-DataSurvival.R b/tests/testthat/test-DataSurvival.R index 85bab6d00..807400cb2 100644 --- a/tests/testthat/test-DataSurvival.R +++ b/tests/testthat/test-DataSurvival.R @@ -105,6 +105,24 @@ test_that("DataSurvival print method works as expected", { }) }) +test_that("model.matrix for DataSurvival works as expected", { + set.seed(123) + x <- data.frame( + vpt = c("b", "a", "c", "d", "e"), + vtime = c(10, 20, 30, 25, 15), + vevent = c(1, 1, 0, 1, 0), + vcov1 = c("A", "A", "B", "B", "A"), + vcov2 = rnorm(5) + ) + + obj <- DataSurvival( + data = x, + formula = Surv(vtime, vevent) ~ vcov1 * vcov2 + ) + + res <- model.matrix(obj) + expect_snapshot_value(res, tolerance = 1e-2, style = "deparse") +}) test_that("mirror_design_matrix() works as expected", { set.seed(3102) diff --git a/vignettes/.gitignore b/vignettes/.gitignore index 097b24163..0a1389f7a 100644 --- a/vignettes/.gitignore +++ b/vignettes/.gitignore @@ -1,2 +1,3 @@ *.html *.R +*_files/ diff --git a/vignettes/quickstart.Rmd b/vignettes/quickstart.Rmd index c1e88b8e9..b342e3a20 100644 --- a/vignettes/quickstart.Rmd +++ b/vignettes/quickstart.Rmd @@ -108,6 +108,8 @@ survival_model <- SurvivalWeibullPH( survival_model ``` +In order to know which position corresponds to which design matrix column and thus coefficient, we can use the `model.matrix()` function on the `DataSurvival` object. This is explained in the section "Formatting Data" below in detail. + ## Separate Models It is also possible to not link the longitudinal and the survival models, by using @@ -256,6 +258,14 @@ joint_data <- DataJoint( ) ``` +Note that we can also actively look into the design matrix of the `DataSurvival` object. This can be helpful if we want to know which coefficients are at which position, e.g. for setting the right prior means and standard deviations with `prior_normal_vector()`. It works like this: + +```{r design_matrix} +head(model.matrix(joint_data@survival)) +``` + +So here we see that the first coefficient corresponds to category `B` and the second coefficient corresponds to category `C` of the categorical covariate `cov_cat`, while the third coefficient corresponds to the continuous covariate `cov_cont`. + # Model Fitting Now let's have a look how we can fit the (joint) models. @@ -289,8 +299,6 @@ mcmc_results <- sampleStanModel( ) ``` - - ## Convergence checks After the sampling finishes, we can inspect the parameter distributions. From cbf05528597c2f9d2c8efabc6a9766e5e22527b0 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 24 Jun 2026 21:38:00 +0800 Subject: [PATCH 18/19] fix docs --- R/DataSurvival.R | 2 ++ man/as_stan_list.DataObject.Rd | 4 ++++ 2 files changed, 6 insertions(+) diff --git a/R/DataSurvival.R b/R/DataSurvival.R index 4175e49be..3d9c3d8bb 100644 --- a/R/DataSurvival.R +++ b/R/DataSurvival.R @@ -117,6 +117,8 @@ as.data.frame.DataSurvival <- function(x, ...) { #' @rdname as_stan_list.DataObject #' @family DataSurvival +#' @param df (`data.frame`)\cr The data frame to be used for generating the design matrix. +#' @param vars (`list`)\cr The list of variable names to be used for generating the design matrix. #' @export model.matrix.DataSurvival <- function( object, diff --git a/man/as_stan_list.DataObject.Rd b/man/as_stan_list.DataObject.Rd index f88e9987e..ead323471 100644 --- a/man/as_stan_list.DataObject.Rd +++ b/man/as_stan_list.DataObject.Rd @@ -46,6 +46,10 @@ \item{subject_var}{(\code{character}) the name of the variable containing the subject identifier.} + +\item{df}{(\code{data.frame})\cr The data frame to be used for generating the design matrix.} + +\item{vars}{(\code{list})\cr The list of variable names to be used for generating the design matrix.} } \description{ Coerces a data object into a \code{list} of data components required From 7b0d59d642f20929165a8e28cb4ea855517216c2 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 24 Jun 2026 21:48:13 +0800 Subject: [PATCH 19/19] fix typo --- tests/testthat/test-populationHR.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-populationHR.R b/tests/testthat/test-populationHR.R index c226f6952..7860c7725 100644 --- a/tests/testthat/test-populationHR.R +++ b/tests/testthat/test-populationHR.R @@ -50,7 +50,7 @@ test_that("populationHR fails for bad input", { }) -test_that("populationHR works as expected for alternative specfications", { +test_that("populationHR works as expected for alternative specifications", { mp <- test_data_1$jsamples set.seed(1231) # Arm + continuous covariate