From 6311e73f36f1bacc27a099c613aa16869e8b6587 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 7 Jul 2026 14:29:13 +0800 Subject: [PATCH 01/30] fix math printing in design doc --- design/covariate-selection.qmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/design/covariate-selection.qmd b/design/covariate-selection.qmd index eacad4c7..93dec3ce 100644 --- a/design/covariate-selection.qmd +++ b/design/covariate-selection.qmd @@ -192,7 +192,7 @@ Then the model assumes: - `b` are the actual regression coefficients, which are obtained by multiplying `zb` with `scales` - The intercept has a Student-t prior - `global` has a half-Student-t prior with `df_global` degrees of freedom and scale `scale_global` -- `slab` has an inverse gamma prior with `df_slab/2` parameters. Note that after multiplying with `scale_slab^2` and obtaining `c2`, this corresponds to $c^2 \sim \text{Inv-Gamma}(\alpha, \beta)$ with $\alpha = \textrm{df_slab}/2$ and $\beta = \textrm{df_slab} * \textrm{scale_slab}^2 / 2$ in the above equations +- `slab` has an inverse gamma prior with `df_slab/2` parameters. Note that after multiplying with `scale_slab^2` and obtaining `c2`, this corresponds to $c^2 \sim \text{Inv-Gamma}(\alpha, \beta)$ with $\alpha = \text{df\_slab}/2$ and $\beta = \text{df\_slab} * \text{scale\_slab}^2 / 2$ in the above equations - `local` parameters have half-Student-t priors with `df` degrees of freedom and scale 1 So mapping this back to the math equations above, this is more general than what was defined above. This model has: From 7e4286635a6689ec8f4a9f8c141050222dfc8ab9 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 7 Jul 2026 15:08:13 +0800 Subject: [PATCH 02/30] add stats spec --- design/covariate-selection.qmd | 5 +- inst/REFERENCES.bib | 51 +++++++++++++++++-- vignettes/statistical-specification.Rmd | 68 ++++++++++++++++++++++++- 3 files changed, 117 insertions(+), 7 deletions(-) diff --git a/design/covariate-selection.qmd b/design/covariate-selection.qmd index 93dec3ce..ea6c4587 100644 --- a/design/covariate-selection.qmd +++ b/design/covariate-selection.qmd @@ -205,10 +205,10 @@ z_{j} &\sim \text{Normal}(0, 1), \\ \xi &\sim \text{Inv-Gamma}\!\left(\frac{\nu_{c}}{2}, \frac{\nu_{c}}{2}\right), \\ c^{2} &= s_{c}^{2} \xi, \\ \tilde{\lambda}_{j} &= -\sqrt{ +\left( \frac{c^{2}\lambda_{j}^{2}} {c^{2} + \tau^{2}\lambda_{j}^{2}} -}, \\ +\right)^{1/2}, \\ \beta_{j} &= z_{j}\tau\tilde{\lambda}_{j}, \quad j = 1, \dotsc, p \end{align*} $$ @@ -255,4 +255,3 @@ In order to make sure that everything is correctly implemented, in addition to s We can also provide the shrinkage factors $\kappa_{j} = 1/(1 + \tau^2\lambda_j^2/c^2)$ for each coefficient, which are a measure of how much shrinkage is applied to each coefficient: $\kappa_{j} \approx 1$ means that the coefficient is shrunk to zero, while $\kappa_{j} \approx 0$ means that the coefficient is not shrunk. Density plots of the shrinkage factors, or even just the median shrinkage factor for each coefficient, can be a useful diagnostic tool to assess which covariates are important and which are not. - diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index 1c300f55..f8c6a48d 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -1,6 +1,3 @@ - - - @article{blanche2015, author = {Blanche, Paul and Proust-Lima, Cécile and Loubère, Lucie and Berr, Claudine and Dartigues, Jean-François and Jacqmin-Gadda, Hélène}, title = {Quantifying and comparing dynamic predictive accuracy of joint models for longitudinal marker and time-to-event in presence of censoring and competing risks}, @@ -51,3 +48,51 @@ @article{oudenhoven2020marginal eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/sim.8713}, year = {2020} } + +@article{CarvalhoEtAl2010, + author = {Carvalho, Carlos M. and Polson, Nicholas G. and Scott, James G.}, + title = {The Horseshoe Estimator for Sparse Signals}, + journal = {Biometrika}, + year = {2010}, + volume = {97}, + number = {2}, + pages = {465--480}, + doi = {10.1093/biomet/asq017}, + url = {https://doi.org/10.1093/biomet/asq017} +} + +@article{PiironenVehtari2017, + author = {Piironen, Juho and Vehtari, Aki}, + title = {Sparsity Information and Regularization in the Horseshoe and Other Shrinkage Priors}, + journal = {Electronic Journal of Statistics}, + year = {2017}, + volume = {11}, + number = {2}, + pages = {5018--5051}, + doi = {10.1214/17-EJS1337SI}, + url = {https://doi.org/10.1214/17-EJS1337SI} +} + +@article{PourzanjaniDavis2026, + author = {Pourzanjani, Arya and Davis, Casey}, + title = {The Regularized Horseshoe for Covariate Selection Improves Convenience and Predictive Performance in Population PK/PD Models}, + journal = {CPT: Pharmacometrics \& Systems Pharmacology}, + year = {2026}, + volume = {15}, + number = {2}, + pages = {e70198}, + doi = {10.1002/psp4.70198}, + url = {https://doi.org/10.1002/psp4.70198} +} + +@article{Buerkner2017, + title = {{brms}: An {R} Package for {Bayesian} Multilevel Models Using {Stan}}, + author = {Paul-Christian Bürkner}, + journal = {Journal of Statistical Software}, + year = {2017}, + volume = {80}, + number = {1}, + pages = {1--28}, + doi = {10.18637/jss.v080.i01}, + encoding = {UTF-8} +} \ No newline at end of file diff --git a/vignettes/statistical-specification.Rmd b/vignettes/statistical-specification.Rmd index 7e7ced79..4ef87995 100644 --- a/vignettes/statistical-specification.Rmd +++ b/vignettes/statistical-specification.Rmd @@ -109,6 +109,73 @@ Where: - $\Gamma(.)$ is the complete gamma function - $\gamma(., .)$ is the lower incomplete gamma function +## Covariate Selection with Regularized Horseshoe Prior + +In the survival model the covariate coefficients $\beta$ typically just have a normal prior distribution, using `prior_normal()` or `prior_normal_vector()`. However, if the user wishes to perform covariate selection then they can instead use the regularized horseshoe prior via `prior_horseshoe()`. Here we explain how the prior is defined and how the hyperparameters can be set. Note that we use the same parametrization and notation as in the `brms` package [@Buerkner2017]. + +The horseshoe prior is a global-local shrinkage prior [@CarvalhoEtAl2010]. The global parameter shrinks all coefficients towards zero, while coefficient-specific local parameters allow some coefficients to escape the shrinkage. The regularized horseshoe prior [@PiironenVehtari2017] additionally introduces a slab component that regularizes coefficients that are far from zero. This is useful for covariate selection because all candidate covariates can be included in one model fit, with irrelevant coefficients strongly shrunk towards zero. + +For $p$ covariate coefficients $\beta_{j}, j = 1, \dotsc, p$, the regularized horseshoe prior applied with `prior_horseshoe()`, and with positive hyperparameters $\nu_{\lambda}, \nu_{\tau}, \nu_c, s_{\tau}, s_c$, can be written as + +$$ +\begin{align*} +z_j &\sim \operatorname{Normal}(0, 1), \\ +\lambda_j &\sim t^+_{\nu_{\lambda}}(0, 1), \\ +\tau &\sim t^+_{\nu_{\tau}}(0, s_{\tau}), \\ +\xi &\sim \text{Inv-Gamma}\!\left(\frac{\nu_c}{2}, \frac{\nu_c}{2}\right), \\ +c^2 &= s_c^2 \xi, \\ +\tilde{\lambda}_j &= +\left( +\frac{c^2\lambda_{j}^2} +{c^2 + \tau^2\lambda_{j}^2} +\right)^{1/2}, \\ +\beta_j &= z_j\tau\tilde{\lambda}_j, \qquad j = 1, \dotsc, p . +\end{align*} +$$ + +Here $t^+_{\nu}(0, s)$ denotes a half-Student-$t$ distribution with degrees of freedom $\nu$ and scale $s$. The $p$ local shrinkage parameters $\lambda_j$ control the amount of shrinkage for each coefficient, while the global shrinkage parameter $\tau$ controls the overall amount of shrinkage across all coefficients. The slab scale $s_c$ and slab degrees of freedom $\nu_c$ determine how strongly large coefficients are regularized. + +Equivalently, conditional on $\lambda_j$, $\tau$, and $c$, this implies + +$$ +\beta_j \mid \lambda_j, \tau, c +\sim +\operatorname{Normal}\!\left( +0,\, +\tau^2 +\frac{c^2\lambda_j^2}{c^2 + \tau^2\lambda_j^2} +\right). +$$ + +The arguments of `prior_horseshoe()` correspond to these hyperparameters: + +- `df`: local degrees of freedom, $\nu_{\lambda}$ +- `df_global`: global degrees of freedom, $\nu_{\tau}$ +- `df_slab`: slab degrees of freedom, $\nu_c$ +- `scale_global`: global scale, $s_{\tau}$ +- `scale_slab`: slab scale, $s_c$ + +As a special case, setting $\nu_{\lambda} = 1$ and $\nu_{\tau} = 1$ gives half-Cauchy priors for the local and global shrinkage parameters. The slab part has the interpretation + +$$ +c^2 \sim \text{Inv-Gamma}\!\left( +\frac{\nu_c}{2}, +\frac{\nu_c s_c^2}{2} +\right), +$$ + +which corresponds to a Student-$t$ slab with $\nu_c$ degrees of freedom and scale $s_c$ for coefficients that escape the main shrinkage. Smaller values of $s_{\tau}$ imply stronger prior shrinkage towards zero, while larger values allow more coefficients to remain non-zero. The slab scale should be chosen large enough to permit clinically or scientifically plausible covariate effects, but not so large that unrealistically large effects are left essentially unregularized. + +A useful summary of the amount of shrinkage applied to coefficient $j$ is the shrinkage factor + +$$ +\kappa_j = +\frac{1}{1 + \tau^2\lambda_j^2 / c^2}. +$$ + +Values of $\kappa_j$ close to 1 indicate strong shrinkage towards zero, while values close to 0 indicate little shrinkage. Regularized horseshoe priors have been applied in PK/PD covariate selection by @PourzanjaniDavis2026, where the main practical advantage is the ability to assess all candidate covariates in a single Bayesian model fit. + + # Longitudinal Model Specification @@ -568,4 +635,3 @@ the censoring distribution). Both of these default options can be changed if req # References - From a6bd3eca45665f3b2b2413c3706bdc66fb2723a4 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 7 Jul 2026 15:42:05 +0800 Subject: [PATCH 03/30] add prior_horseshoe - first stab --- NAMESPACE | 1 + R/Parameter.R | 2 +- R/Prior.R | 137 ++++++++++++++++++++++++++++++++++-- inst/stan/base/base.stan | 3 +- man/prior_horseshoe.Rd | 45 ++++++++++++ tests/testthat/test-Prior.R | 89 +++++++++++++++++++++++ 6 files changed, 271 insertions(+), 6 deletions(-) create mode 100644 man/prior_horseshoe.Rd diff --git a/NAMESPACE b/NAMESPACE index 34924da3..87a98bfe 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -267,6 +267,7 @@ export(populationHR) export(prior_beta) export(prior_cauchy) export(prior_gamma) +export(prior_horseshoe) export(prior_init_only) export(prior_invgamma) export(prior_logistic) diff --git a/R/Parameter.R b/R/Parameter.R index 1db50bba..73e12b9f 100755 --- a/R/Parameter.R +++ b/R/Parameter.R @@ -86,7 +86,7 @@ setValidity( #' @family as.StanModule #' @export as.StanModule.Parameter <- function(object, ...) { - as.StanModule(object@prior, name = object@name) + as.StanModule(object@prior, name = object@name, size = object@size) } diff --git a/R/Prior.R b/R/Prior.R index b9633984..f1b4ee06 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -31,6 +31,8 @@ NULL #' @slot parameters (`list`)\cr See arguments. #' @slot repr_model (`string`)\cr See arguments. #' @slot repr_data (`string`)\cr See arguments. +#' @slot repr_parameters (`string`)\cr See arguments. +#' @slot repr_transformed_parameters (`string`)\cr See arguments. #' @slot centre (`numeric`)\cr See arguments. #' @slot validation (`list`)\cr See arguments. #' @slot display (`string`)\cr See arguments. @@ -48,6 +50,8 @@ NULL "display" = "character", "repr_model" = "character", "repr_data" = "character", + "repr_parameters" = "character", + "repr_transformed_parameters" = "character", "centre" = "numeric", "validation" = "list", "sample" = "function", @@ -63,6 +67,10 @@ NULL #' the Stan code representation for the model block. #' @typed repr_data: string #' the Stan code representation for the data block. +#' @typed repr_parameters: string +#' the Stan code representation for the parameters block. +#' @typed repr_transformed_parameters: string +#' the Stan code representation for the transformed parameters block. #' @typed display: string #' the string to display when object is printed. #' @typed centre: numeric @@ -85,6 +93,8 @@ Prior <- function( centre, validation, sample, + repr_parameters = "", + repr_transformed_parameters = "", limits = c(-Inf, Inf), .allow_vectors = FALSE ) { @@ -92,6 +102,8 @@ Prior <- function( parameters = parameters, repr_model = repr_model, repr_data = repr_data, + repr_parameters = repr_parameters, + repr_transformed_parameters = repr_transformed_parameters, centre = centre, display = display, validation = validation, @@ -150,6 +162,20 @@ setValidity( ) { return("Model representation must be length 1 string") } + if ( + length(object@repr_parameters) < 1 || + !is.character(object@repr_parameters) + ) { + return("Parameters representation must be a character vector") + } + if ( + length(object@repr_transformed_parameters) < 1 || + !is.character(object@repr_transformed_parameters) + ) { + return( + "Transformed parameters representation must be a character vector" + ) + } return(TRUE) } ) @@ -234,7 +260,14 @@ setMethod( #' @family Prior-internal #' @family as.StanModule #' @export -as.StanModule.Prior <- function(object, name, ...) { +as.StanModule.Prior <- function(object, name, size = 1, ...) { + indent_stan <- function(x) { + x <- x[nchar(x) >= 1] + if (length(x) == 0) { + return("") + } + paste0(" ", x, collapse = "\n") + } trunctation <- if (object@repr_model != "") { paste0(render_stan_limits(object@limits), ";") } else { @@ -242,14 +275,20 @@ as.StanModule.Prior <- function(object, name, ...) { } string <- paste( "data {{", - paste0(" ", object@repr_data, collapse = "\n"), + indent_stan(object@repr_data), + "}}", + "parameters {{", + indent_stan(object@repr_parameters), + "}}", + "transformed parameters {{", + indent_stan(object@repr_transformed_parameters), "}}", "model {{", - paste0(" ", object@repr_model, trunctation), + indent_stan(paste0(object@repr_model, trunctation)), "}}", sep = "\n" ) - StanModule(glue::glue(string, name = name)) + StanModule(glue::glue(string, name = name, size = size)) } @@ -415,6 +454,86 @@ prior_normal_vector <- function(mus, sigmas) { ) } +#' Regularized Horseshoe Prior for a Vector Distribution +#' +#' @typed df: number +#' degrees of freedom of the half-Student-t prior for local shrinkage +#' parameters. +#' @typed df_global: number +#' degrees of freedom of the half-Student-t prior for the global shrinkage +#' parameter. +#' @typed df_slab: number +#' degrees of freedom of the Student-t slab. +#' @typed scale_global: number +#' scale of the half-Student-t prior for the global shrinkage parameter. +#' @typed scale_slab: number +#' scale of the Student-t slab. +#' @family Prior +#' @export +prior_horseshoe <- function( + df = 1, + df_global = 1, + df_slab = 4, + scale_global = 1, + scale_slab = 2 +) { + Prior( + parameters = list( + df = df, + df_global = df_global, + df_slab = df_slab, + scale_global = scale_global, + scale_slab = scale_slab + ), + display = paste0( + "horseshoe(df = {df}, df_global = {df_global}, ", + "df_slab = {df_slab}, scale_global = {scale_global}, ", + "scale_slab = {scale_slab})" + ), + repr_model = paste( + "prior_local_{name} ~ student_t(prior_df_{name}, 0, 1);", + "prior_global_{name} ~ student_t(prior_df_global_{name}, 0, prior_scale_global_{name});", + "prior_slab_{name} ~ inv_gamma(prior_df_slab_{name} / 2, prior_df_slab_{name} / 2);", + "{name} ~ normal(rep_vector(0, {size}), prior_scales_{name})", + sep = "\n " + ), + repr_data = c( + "real prior_df_{name};", + "real prior_df_global_{name};", + "real prior_df_slab_{name};", + "real prior_scale_global_{name};", + "real prior_scale_slab_{name};" + ), + repr_parameters = c( + "vector[{size}] prior_local_{name};", + "real prior_global_{name};", + "real prior_slab_{name};" + ), + repr_transformed_parameters = c( + "real prior_c2_{name} = square(prior_scale_slab_{name}) * prior_slab_{name};", + "vector[{size}] prior_scales_{name} = prior_global_{name} * sqrt((prior_c2_{name} * square(prior_local_{name})) ./ (prior_c2_{name} + square(prior_global_{name}) * square(prior_local_{name})));" + ), + centre = 0, + sample = \(n) { + local_rhorseshoe( + n = n, + df = df, + df_global = df_global, + df_slab = df_slab, + scale_global = scale_global, + scale_slab = scale_slab + ) + }, + validation = list( + df = \(x) x > 0, + df_global = \(x) x > 0, + df_slab = \(x) x > 0, + scale_global = \(x) x > 0, + scale_slab = \(x) x > 0 + ) + ) +} + #' Standard Normal Prior Distribution #' @@ -827,6 +946,16 @@ local_rnorm_vector <- \(n, mus, sigmas) { mapply(local_rnorm, n = n, mean = mus, sd = sigmas) } +#' @rdname Local_Sample +local_rhorseshoe <- \(n, df, df_global, df_slab, scale_global, scale_slab) { + local <- abs(local_rt(n, df, 0, 1)) + global <- abs(local_rt(n, df_global, 0, scale_global)) + slab <- local_rinvgamma(n, df_slab / 2, df_slab / 2) + c2 <- scale_slab^2 * slab + scales <- global * sqrt((c2 * local^2) / (c2 + global^2 * local^2)) + local_rnorm(n, mean = 0, sd = scales) +} + #' @rdname Local_Sample local_rcauchy <- \(...) rcauchy(...) diff --git a/inst/stan/base/base.stan b/inst/stan/base/base.stan index a446cfdf..d611b966 100755 --- a/inst/stan/base/base.stan +++ b/inst/stan/base/base.stan @@ -38,6 +38,7 @@ parameters{ {{ longitudinal.parameters }} {{ link.parameters }} {{ survival.parameters }} +{{ priors.parameters }} } @@ -46,6 +47,7 @@ transformed parameters{ {{ longitudinal.transformed_parameters }} {{ link.transformed_parameters }} {{ survival.transformed_parameters }} +{{ priors.transformed_parameters }} // // Source - base/base.stan @@ -71,4 +73,3 @@ generated quantities{ {{ survival.generated_quantities }} } - diff --git a/man/prior_horseshoe.Rd b/man/prior_horseshoe.Rd new file mode 100644 index 00000000..ef27d1a7 --- /dev/null +++ b/man/prior_horseshoe.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Prior.R +\name{prior_horseshoe} +\alias{prior_horseshoe} +\title{Regularized Horseshoe Prior for a Vector Distribution} +\usage{ +prior_horseshoe( + df = 1, + df_global = 1, + df_slab = 4, + scale_global = 1, + scale_slab = 2 +) +} +\arguments{ +\item{df}{(\code{number}) degrees of freedom of the half-Student-t prior for local shrinkage parameters.} + +\item{df_global}{(\code{number}) degrees of freedom of the half-Student-t prior for the global shrinkage parameter.} + +\item{df_slab}{(\code{number}) degrees of freedom of the Student-t slab.} + +\item{scale_global}{(\code{number}) scale of the half-Student-t prior for the global shrinkage parameter.} + +\item{scale_slab}{(\code{number}) scale of the Student-t slab.} +} +\description{ +Regularized Horseshoe 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_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()}} +} +\concept{Prior} diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 28a4d6e8..7d340845 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -269,6 +269,16 @@ test_that("Parameters in priors must be length 1 #422", { prior_gamma(c(1, 2), 2), "Parameter `alpha`" ) + + expect_error( + prior_horseshoe(df = -1), + regexp = "Invalid.*`df`" + ) + + expect_error( + prior_horseshoe(scale_global = -1), + regexp = "Invalid.*`scale_global`" + ) }) test_that("prior_normal_vector works as expected", { @@ -298,3 +308,82 @@ test_that("prior_normal_vector works as expected", { ) ) }) + +test_that("prior_horseshoe works as expected", { + x <- prior_horseshoe( + df = 1, + df_global = 2, + df_slab = 3, + scale_global = 0.4, + scale_slab = 2 + ) + + expect_equal( + as.character(x), + paste0( + "horseshoe(df = 1, df_global = 2, df_slab = 3, ", + "scale_global = 0.4, scale_slab = 2)" + ) + ) + + x_inits <- initialValues(x) + expect_numeric(x_inits, len = 1) + + x_stan_module <- as.StanModule(x, name = "beta", size = "p") + expect_equal( + x_stan_module@data, + c( + " real prior_df_beta;", + " real prior_df_global_beta;", + " real prior_df_slab_beta;", + " real prior_scale_global_beta;", + " real prior_scale_slab_beta;" + ) + ) + expect_equal( + x_stan_module@parameters, + c( + " vector[p] prior_local_beta;", + " real prior_global_beta;", + " real prior_slab_beta;" + ) + ) + expect_true( + any(grepl( + "vector\\[p\\] prior_scales_beta", + x_stan_module@transformed_parameters + )) + ) + expect_true( + any(grepl( + "beta ~ normal(rep_vector(0, p), prior_scales_beta);", + x_stan_module@model, + fixed = TRUE + )) + ) + expect_equal( + as_stan_list(x, name = "beta"), + list( + prior_df_beta = 1, + prior_df_global_beta = 2, + prior_df_slab_beta = 3, + prior_scale_global_beta = 0.4, + prior_scale_slab_beta = 2 + ) + ) + + header <- StanModule( + "data { + int p; +} +parameters { + vector[p] beta; +}" + ) + stan_file <- cmdstanr::write_stan_file( + as.character(merge(header, x_stan_module)), + dir = tempdir() + ) + model_obj <- cmdstanr::cmdstan_model(stan_file, compile = FALSE) + expect_true(model_obj$check_syntax(quiet = TRUE)) +}) From 874f7c662d8f2d51e6f0727e990f31acfda94fdb Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 7 Jul 2026 15:49:11 +0800 Subject: [PATCH 04/30] explain syntax test part --- tests/testthat/test-Prior.R | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 7d340845..c2b3f153 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -372,6 +372,7 @@ test_that("prior_horseshoe works as expected", { ) ) + # Check the model syntax. header <- StanModule( "data { int p; From a5f31d2c1ccf1d003592e3ca262eb0dfb834f127 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 7 Jul 2026 15:57:15 +0800 Subject: [PATCH 05/30] add a test for it --- tests/testthat/test-misc_models.R | 38 +++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/tests/testthat/test-misc_models.R b/tests/testthat/test-misc_models.R index 8eeddc2f..818c755e 100644 --- a/tests/testthat/test-misc_models.R +++ b/tests/testthat/test-misc_models.R @@ -25,3 +25,41 @@ test_that("Longitudinal Model doesn't print sampler rejection messages", { mp ))) }) + +test_that("Survival model can be fit with a horseshoe prior", { + jm <- JointModel( + survival = SurvivalWeibullPH( + beta = prior_horseshoe( + df = 1, + df_global = 1, + df_slab = 4, + scale_global = 0.3, + scale_slab = 2 + ) + ) + ) + + mp <- run_quietly({ + sampleStanModel( + jm, + data = test_data_1$jdata, + iter_sampling = 3, + iter_warmup = 3, + chains = 1, + refresh = 0, + parallel_chains = 1, + seed = 325 + ) + }) + + draws <- cmdstanr::as.CmdStanMCMC(mp)$draws() + variables <- posterior::variables(draws) + + expect_true(all(c( + "beta_os_cov[1]", + "prior_local_beta_os_cov[1]", + "prior_global_beta_os_cov", + "prior_slab_beta_os_cov", + "prior_scales_beta_os_cov[1]" + ) %in% variables)) +}) From f0fe6147a96301f1dcc37fa84e35e621290d5013 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 7 Jul 2026 16:21:08 +0800 Subject: [PATCH 06/30] add shrinkage factors as transformed parameters --- R/Prior.R | 1 + tests/testthat/test-Prior.R | 13 +++++++++++++ tests/testthat/test-misc_models.R | 1 + 3 files changed, 15 insertions(+) diff --git a/R/Prior.R b/R/Prior.R index f1b4ee06..f9ce1315 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -511,6 +511,7 @@ prior_horseshoe <- function( ), repr_transformed_parameters = c( "real prior_c2_{name} = square(prior_scale_slab_{name}) * prior_slab_{name};", + "vector[{size}] prior_shrinkage_factors_{name} = rep_vector(1, {size}) ./ (rep_vector(1, {size}) + square(prior_global_{name}) * square(prior_local_{name}) / prior_c2_{name});", "vector[{size}] prior_scales_{name} = prior_global_{name} * sqrt((prior_c2_{name} * square(prior_local_{name})) ./ (prior_c2_{name} + square(prior_global_{name}) * square(prior_local_{name})));" ), centre = 0, diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index c2b3f153..1caffc41 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -354,6 +354,19 @@ test_that("prior_horseshoe works as expected", { x_stan_module@transformed_parameters )) ) + expect_true( + any(grepl( + "vector\\[p\\] prior_shrinkage_factors_beta", + x_stan_module@transformed_parameters + )) + ) + expect_true( + any(grepl( + "rep_vector(1, p) ./ (rep_vector(1, p) + square(prior_global_beta) * square(prior_local_beta) / prior_c2_beta)", + x_stan_module@transformed_parameters, + fixed = TRUE + )) + ) expect_true( any(grepl( "beta ~ normal(rep_vector(0, p), prior_scales_beta);", diff --git a/tests/testthat/test-misc_models.R b/tests/testthat/test-misc_models.R index 818c755e..e0ad9ac1 100644 --- a/tests/testthat/test-misc_models.R +++ b/tests/testthat/test-misc_models.R @@ -60,6 +60,7 @@ test_that("Survival model can be fit with a horseshoe prior", { "prior_local_beta_os_cov[1]", "prior_global_beta_os_cov", "prior_slab_beta_os_cov", + "prior_shrinkage_factors_beta_os_cov[1]", "prior_scales_beta_os_cov[1]" ) %in% variables)) }) From daabd5390d3de1f5eb3c039b0d54dc50b82dd908 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 7 Jul 2026 16:39:56 +0800 Subject: [PATCH 07/30] make sure all shrinkage factors in samples --- tests/testthat/test-misc_models.R | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/tests/testthat/test-misc_models.R b/tests/testthat/test-misc_models.R index e0ad9ac1..8cfe5284 100644 --- a/tests/testthat/test-misc_models.R +++ b/tests/testthat/test-misc_models.R @@ -55,12 +55,17 @@ test_that("Survival model can be fit with a horseshoe prior", { draws <- cmdstanr::as.CmdStanMCMC(mp)$draws() variables <- posterior::variables(draws) - expect_true(all(c( - "beta_os_cov[1]", - "prior_local_beta_os_cov[1]", - "prior_global_beta_os_cov", - "prior_slab_beta_os_cov", - "prior_shrinkage_factors_beta_os_cov[1]", - "prior_scales_beta_os_cov[1]" - ) %in% variables)) + expect_true(all( + c( + "beta_os_cov[1]", + "prior_local_beta_os_cov[1]", + "prior_global_beta_os_cov", + "prior_slab_beta_os_cov", + "prior_shrinkage_factors_beta_os_cov[1]", + "prior_shrinkage_factors_beta_os_cov[2]", + "prior_shrinkage_factors_beta_os_cov[3]", + "prior_scales_beta_os_cov[1]" + ) %in% + variables + )) }) From 064c551dd789a5aeaf0ab6ff24d9d6ab692f2781 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 7 Jul 2026 21:29:05 +0800 Subject: [PATCH 08/30] add covariates and shrinkage gen. functions and methods, run docs --- NAMESPACE | 4 ++++ R/DataSurvival.R | 12 ++++++++++++ R/JointModelSamples.R | 26 ++++++++++++++++++++++++++ R/generics.R | 22 ++++++++++++++++++++++ man/Local_Sample.Rd | 3 +++ man/Prior-class.Rd | 10 ++++++++++ man/as.StanModule.Prior.Rd | 2 +- man/covariates.DataSurvival.Rd | 19 +++++++++++++++++++ man/covariates.Rd | 19 +++++++++++++++++++ man/prior_beta.Rd | 1 + man/prior_cauchy.Rd | 1 + man/prior_gamma.Rd | 1 + man/prior_horseshoe.Rd | 6 ++++-- 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 | 1 + man/prior_std_normal.Rd | 1 + man/prior_student_t.Rd | 1 + man/prior_uniform.Rd | 1 + man/shrinkage.JointModelSamples.Rd | 22 ++++++++++++++++++++++ man/shrinkage.Rd | 19 +++++++++++++++++++ tests/testthat/test-misc_models.R | 25 +++++++++++++++++++++---- 26 files changed, 195 insertions(+), 7 deletions(-) create mode 100644 man/covariates.DataSurvival.Rd create mode 100644 man/covariates.Rd create mode 100644 man/shrinkage.JointModelSamples.Rd create mode 100644 man/shrinkage.Rd diff --git a/NAMESPACE b/NAMESPACE index 87a98bfe..78e44b94 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -87,6 +87,7 @@ S3method(coalesceGridTime,GridPopulation) S3method(coalesceGridTime,GridPrediction) S3method(coalesceGridTime,default) S3method(compileStanModel,JointModel) +S3method(covariates,DataSurvival) S3method(createLongitudinalSimObject,LongitudinalClaretBruno) S3method(createLongitudinalSimObject,LongitudinalGSF) S3method(createLongitudinalSimObject,LongitudinalRandomSlope) @@ -181,6 +182,7 @@ S3method(sampleSubjects,SimLongitudinalSteinFojo) S3method(sampleSubjects,SimSurvival) S3method(saveObject,JointModelSamples) S3method(set_limits,Prior) +S3method(shrinkage,JointModelSamples) S3method(simulate,JointModelSamples) S3method(size,Parameter) S3method(size,ParameterList) @@ -248,6 +250,7 @@ export(as_stan_list) export(autoplot) export(brierScore) export(compileStanModel) +export(covariates) export(cut_data) export(enableGQ) export(enableLink) @@ -285,6 +288,7 @@ export(sampleSubjects) export(saveObject) export(set_limits) export(show) +export(shrinkage) export(write_stan) exportClasses(DataJoint) exportClasses(DataLongitudinal) diff --git a/R/DataSurvival.R b/R/DataSurvival.R index 3d9c3d8b..32c0916c 100644 --- a/R/DataSurvival.R +++ b/R/DataSurvival.R @@ -133,6 +133,18 @@ model.matrix.DataSurvival <- function( design_mat } +#' Extract Covariate Names for Survival Data +#' +#' @typed object: DataSurvival +#' Survival data object to extract covariate names from. +#' +#' @family covariates +#' @export +covariates.DataSurvival <- function(object, ...) { + design_mat <- model.matrix(object, ...) + colnames(design_mat) +} + #' @rdname as_stan_list.DataObject #' @family DataSurvival #' @export diff --git a/R/JointModelSamples.R b/R/JointModelSamples.R index 697ff21e..131990b2 100644 --- a/R/JointModelSamples.R +++ b/R/JointModelSamples.R @@ -183,3 +183,29 @@ saveObject.JointModelSamples <- function(object, file, ...) { try(object@results$profiles(), silent = TRUE) saveRDS(object, file, ...) } + +#' Extract shrinkage factors from a `JointModelSamples` object +#' +#' @typed object: JointModelSamples +#' the object to extract shrinkage factors from. +#' @return the shrinkage factors correctly named after the covariate names. +#' +#' @family shrinkage +#' @export +shrinkage.JointModelSamples <- function(object, ...) { + draws <- cmdstanr::as.CmdStanMCMC(object)$draws() + assert_class(object@data@survival, "DataSurvival") + covariate_names <- covariates(object@data@survival) + assert_subset( + "prior_shrinkage_factors_beta_os_cov[1]", + posterior::variables(draws) + ) + shrinkage_draws <- subset( + draws, + variable = "prior_shrinkage_factors_beta_os_cov" + ) + shrinkage_vars <- posterior::variables(shrinkage_draws) + assert_true(identical(length(covariate_names), length(shrinkage_vars))) + posterior::variables(shrinkage_draws) <- covariate_names + shrinkage_draws +} diff --git a/R/generics.R b/R/generics.R index 07f63451..f11fad8b 100755 --- a/R/generics.R +++ b/R/generics.R @@ -526,3 +526,25 @@ set_limits <- function(object, lower = -Inf, upper = Inf) { saveObject <- function(object, file, ...) { UseMethod("saveObject") } + + +#' Extract Covariate Names +#' +#' @typed object: ANY +#' +#' @family covariates +#' @export +covariates <- function(object, ...) { + UseMethod("covariates") +} + +#' Extract Covariate Shrinkage Factors from the Posterior Samples +#' +#' @typed object: ANY +#' the object to extract shrinkage factors from. +#' +#' @family shrinkage +#' @export +shrinkage <- function(object, ...) { + UseMethod("shrinkage") +} diff --git a/man/Local_Sample.Rd b/man/Local_Sample.Rd index 6bbe1459..c3371fc0 100644 --- a/man/Local_Sample.Rd +++ b/man/Local_Sample.Rd @@ -4,6 +4,7 @@ \alias{Local_Sample} \alias{local_rnorm} \alias{local_rnorm_vector} +\alias{local_rhorseshoe} \alias{local_rcauchy} \alias{local_rgamma} \alias{local_rlnorm} @@ -19,6 +20,8 @@ local_rnorm(...) local_rnorm_vector(n, mus, sigmas) +local_rhorseshoe(n, df, df_global, df_slab, scale_global, scale_slab) + local_rcauchy(...) local_rgamma(...) diff --git a/man/Prior-class.Rd b/man/Prior-class.Rd index 490da7d5..9b964890 100644 --- a/man/Prior-class.Rd +++ b/man/Prior-class.Rd @@ -15,6 +15,8 @@ Prior( centre, validation, sample, + repr_parameters = "", + repr_transformed_parameters = "", limits = c(-Inf, Inf), .allow_vectors = FALSE ) @@ -35,6 +37,10 @@ the same names as the \code{paramaters} slot.} \item{sample}{(\code{function}) a function to sample from the prior distribution.} +\item{repr_parameters}{(\code{string}) the Stan code representation for the parameters block.} + +\item{repr_transformed_parameters}{(\code{string}) the Stan code representation for the transformed parameters block.} + \item{limits}{(\code{numeric}) the lower and upper limits for a truncated distribution} \item{.allow_vectors}{(\code{flag}) whether to allow vector parameters.} @@ -51,6 +57,10 @@ Specifies the prior distribution in a Stan Model \item{\code{repr_data}}{(\code{string})\cr See arguments.} +\item{\code{repr_parameters}}{(\code{string})\cr See arguments.} + +\item{\code{repr_transformed_parameters}}{(\code{string})\cr See arguments.} + \item{\code{centre}}{(\code{numeric})\cr See arguments.} \item{\code{validation}}{(\code{list})\cr See arguments.} diff --git a/man/as.StanModule.Prior.Rd b/man/as.StanModule.Prior.Rd index e018daf1..ad90e95f 100644 --- a/man/as.StanModule.Prior.Rd +++ b/man/as.StanModule.Prior.Rd @@ -4,7 +4,7 @@ \alias{as.StanModule.Prior} \title{\code{Prior} -> \code{StanModule}} \usage{ -\method{as.StanModule}{Prior}(object, name, ...) +\method{as.StanModule}{Prior}(object, name, size = 1, ...) } \arguments{ \item{object}{(\code{Prior}) a prior Distribution} diff --git a/man/covariates.DataSurvival.Rd b/man/covariates.DataSurvival.Rd new file mode 100644 index 00000000..0cea798d --- /dev/null +++ b/man/covariates.DataSurvival.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DataSurvival.R +\name{covariates.DataSurvival} +\alias{covariates.DataSurvival} +\title{Extract Covariate Names for Survival Data} +\usage{ +\method{covariates}{DataSurvival}(object, ...) +} +\arguments{ +\item{object}{(\code{DataSurvival}) Survival data object to extract covariate names from.} +} +\description{ +Extract Covariate Names for Survival Data +} +\seealso{ +Other covariates: +\code{\link[=covariates]{covariates()}} +} +\concept{covariates} diff --git a/man/covariates.Rd b/man/covariates.Rd new file mode 100644 index 00000000..6d183a02 --- /dev/null +++ b/man/covariates.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/generics.R +\name{covariates} +\alias{covariates} +\title{Extract Covariate Names} +\usage{ +covariates(object, ...) +} +\arguments{ +\item{}{(\if{html}{\out{}}) \if{html}{\out{}}} +} +\description{ +Extract Covariate Names +} +\seealso{ +Other covariates: +\code{\link[=covariates.DataSurvival]{covariates.DataSurvival()}} +} +\concept{covariates} diff --git a/man/prior_beta.Rd b/man/prior_beta.Rd index 1479a2e4..e2bf46ce 100644 --- a/man/prior_beta.Rd +++ b/man/prior_beta.Rd @@ -18,6 +18,7 @@ Beta Prior Distribution Other Prior: \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_cauchy.Rd b/man/prior_cauchy.Rd index c559b096..22864916 100644 --- a/man/prior_cauchy.Rd +++ b/man/prior_cauchy.Rd @@ -18,6 +18,7 @@ Cauchy Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_gamma.Rd b/man/prior_gamma.Rd index cc386c88..daa80b47 100644 --- a/man/prior_gamma.Rd +++ b/man/prior_gamma.Rd @@ -18,6 +18,7 @@ Gamma Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_horseshoe.Rd b/man/prior_horseshoe.Rd index ef27d1a7..d5dfceb5 100644 --- a/man/prior_horseshoe.Rd +++ b/man/prior_horseshoe.Rd @@ -13,9 +13,11 @@ prior_horseshoe( ) } \arguments{ -\item{df}{(\code{number}) degrees of freedom of the half-Student-t prior for local shrinkage parameters.} +\item{df}{(\code{number}) degrees of freedom of the half-Student-t prior for local shrinkage +parameters.} -\item{df_global}{(\code{number}) degrees of freedom of the half-Student-t prior for the global shrinkage parameter.} +\item{df_global}{(\code{number}) degrees of freedom of the half-Student-t prior for the global shrinkage +parameter.} \item{df_slab}{(\code{number}) degrees of freedom of the Student-t slab.} diff --git a/man/prior_init_only.Rd b/man/prior_init_only.Rd index b36e9153..f3c0019d 100644 --- a/man/prior_init_only.Rd +++ b/man/prior_init_only.Rd @@ -19,6 +19,7 @@ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, \code{\link[=prior_loglogistic]{prior_loglogistic()}}, diff --git a/man/prior_invgamma.Rd b/man/prior_invgamma.Rd index c3a7f3d0..7ad4f91a 100644 --- a/man/prior_invgamma.Rd +++ b/man/prior_invgamma.Rd @@ -19,6 +19,7 @@ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_logistic]{prior_logistic()}}, \code{\link[=prior_loglogistic]{prior_loglogistic()}}, diff --git a/man/prior_logistic.Rd b/man/prior_logistic.Rd index fc0447c1..478873ba 100644 --- a/man/prior_logistic.Rd +++ b/man/prior_logistic.Rd @@ -19,6 +19,7 @@ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_loglogistic]{prior_loglogistic()}}, diff --git a/man/prior_loglogistic.Rd b/man/prior_loglogistic.Rd index 62a46e5d..986907dc 100644 --- a/man/prior_loglogistic.Rd +++ b/man/prior_loglogistic.Rd @@ -19,6 +19,7 @@ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_lognormal.Rd b/man/prior_lognormal.Rd index 86335798..3354db94 100644 --- a/man/prior_lognormal.Rd +++ b/man/prior_lognormal.Rd @@ -19,6 +19,7 @@ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_normal.Rd b/man/prior_normal.Rd index 2f532e8c..d9b1c2ba 100644 --- a/man/prior_normal.Rd +++ b/man/prior_normal.Rd @@ -19,6 +19,7 @@ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_normal_vector.Rd b/man/prior_normal_vector.Rd index 5d7fbb7f..8c10cd38 100644 --- a/man/prior_normal_vector.Rd +++ b/man/prior_normal_vector.Rd @@ -19,6 +19,7 @@ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_std_normal.Rd b/man/prior_std_normal.Rd index 46f67d4f..3fe3f752 100644 --- a/man/prior_std_normal.Rd +++ b/man/prior_std_normal.Rd @@ -14,6 +14,7 @@ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_student_t.Rd b/man/prior_student_t.Rd index 1101043e..b8607b21 100644 --- a/man/prior_student_t.Rd +++ b/man/prior_student_t.Rd @@ -21,6 +21,7 @@ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_uniform.Rd b/man/prior_uniform.Rd index 15efa171..6f3ae289 100644 --- a/man/prior_uniform.Rd +++ b/man/prior_uniform.Rd @@ -19,6 +19,7 @@ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, \code{\link[=prior_gamma]{prior_gamma()}}, +\code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/shrinkage.JointModelSamples.Rd b/man/shrinkage.JointModelSamples.Rd new file mode 100644 index 00000000..72c92543 --- /dev/null +++ b/man/shrinkage.JointModelSamples.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/JointModelSamples.R +\name{shrinkage.JointModelSamples} +\alias{shrinkage.JointModelSamples} +\title{Extract shrinkage factors from a \code{JointModelSamples} object} +\usage{ +\method{shrinkage}{JointModelSamples}(object, ...) +} +\arguments{ +\item{object}{(\code{JointModelSamples}) the object to extract shrinkage factors from.} +} +\value{ +the shrinkage factors correctly named after the covariate names. +} +\description{ +Extract shrinkage factors from a \code{JointModelSamples} object +} +\seealso{ +Other shrinkage: +\code{\link[=shrinkage]{shrinkage()}} +} +\concept{shrinkage} diff --git a/man/shrinkage.Rd b/man/shrinkage.Rd new file mode 100644 index 00000000..9a6d9b07 --- /dev/null +++ b/man/shrinkage.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/generics.R +\name{shrinkage} +\alias{shrinkage} +\title{Extract Covariate Shrinkage Factors from the Posterior Samples} +\usage{ +shrinkage(object, ...) +} +\arguments{ +\item{object}{(\code{ANY}) the object to extract shrinkage factors from.} +} +\description{ +Extract Covariate Shrinkage Factors from the Posterior Samples +} +\seealso{ +Other shrinkage: +\code{\link[=shrinkage.JointModelSamples]{shrinkage.JointModelSamples()}} +} +\concept{shrinkage} diff --git a/tests/testthat/test-misc_models.R b/tests/testthat/test-misc_models.R index 8cfe5284..58fce01c 100644 --- a/tests/testthat/test-misc_models.R +++ b/tests/testthat/test-misc_models.R @@ -43,11 +43,11 @@ test_that("Survival model can be fit with a horseshoe prior", { sampleStanModel( jm, data = test_data_1$jdata, - iter_sampling = 3, - iter_warmup = 3, - chains = 1, + iter_sampling = 1000, + iter_warmup = 100, + chains = 4, refresh = 0, - parallel_chains = 1, + parallel_chains = 4, seed = 325 ) }) @@ -68,4 +68,21 @@ test_that("Survival model can be fit with a horseshoe prior", { ) %in% variables )) + + # Look at the shrinkage factors in more detail now. + shrinkage_draws <- draws |> + subset( + variable = "prior_shrinkage_factors_beta_os_cov" + ) + shrinkage_summary <- summary(shrinkage_draws) + expect_true(all(shrinkage_summary$rhat < 1.05)) + shrinkage_with_method <- shrinkage(mp) + expect_equal( + posterior::variables(shrinkage_with_method), + c("cov_catB", "cov_catC", "cov_cont") + ) + shrinkage_summary2 <- summary(shrinkage_with_method) + expect_equal(shrinkage_summary$median, shrinkage_summary2$median) + # Now one could also do density plots etc. to visualize the shrinkage factors + # with the bayesplot package. }) From 5fa4f38efb7da086c98af0b0a2eb085ac3039cb6 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 7 Jul 2026 21:58:57 +0800 Subject: [PATCH 09/30] first stab at comparison with brms but does not run yet --- tests/testthat/test-brms_compare.R | 128 +++++++++++++++++++++++++++++ 1 file changed, 128 insertions(+) diff --git a/tests/testthat/test-brms_compare.R b/tests/testthat/test-brms_compare.R index 8d1ec9c8..67b8dbc8 100644 --- a/tests/testthat/test-brms_compare.R +++ b/tests/testthat/test-brms_compare.R @@ -266,3 +266,131 @@ test_that("jmpost and brms get similar loo for survival models", { cor(b_est$pointwise[, "looic"], j_est$pointwise[, "looic"]) > 0.95 ) }) + + +test_that("jmpost and brms get similar horseshoe estimates for survival models", { + skip_if_not(is_full_test()) + + set.seed(9826) + n <- 500 + dat_surv <- dplyr::tibble( + pt = sprintf("pt_%06i", seq_len(n)), + lambda_0 = 365 / 150, + ones = 1, + cov1 = rnorm(n), + cov2 = rnorm(n), + cov3 = rnorm(n), + lambda = lambda_0 * exp(cov1 * 0.5 + cov2 * -0.3), + time = rexp(n, 1 / lambda), + event = 1 + ) + + mp_brms <- brms::brm( + time ~ 0 + ones + cov1 + cov2 + cov3, # No intercept and no centering, please. + family = brms::exponential(), + data = dat_surv, + prior = brms::prior( + "horseshoe(df = 1, df_global = 1, df_slab = 4, scale_global = 0.3, scale_slab = 2)", + class = "b" + ), + warmup = 5000, + iter = 10000, + chains = 4, + cores = 4, + seed = 9826, + backend = "cmdstanr" + ) + + dat_surv2 <- dat_surv |> + dplyr::mutate(arm = "A") |> + dplyr::mutate(study = "A") + + dat_surv_bl <- dat_surv2 |> + dplyr::select(pt, arm, study) |> + dplyr::group_by(pt) |> + dplyr::slice(1) |> + dplyr::ungroup() + + jdat <- DataJoint( + subject = DataSubject( + data = dat_surv_bl, + subject = "pt", + arm = "arm", + study = "study" + ), + survival = DataSurvival( + data = dat_surv2, + formula = Surv(time, event) ~ 0 + ones + cov1 + cov2 + cov3 + ) + ) + + jm <- JointModel( + survival = SurvivalExponential( + beta = prior_horseshoe( + df = 1, + df_global = 1, + df_slab = 4, + scale_global = 0.3, + scale_slab = 2 + ) + ) + ) + + mp <- sampleStanModel( + jm, + data = jdat, + iter_warmup = 5000, + iter_sampling = 5000, + chains = 4, + refresh = 200, + parallel_chains = 4, + seed = 9826 + ) + + b_draws <- posterior::as_draws(mp_brms) + stanmod <- cmdstanr::as.CmdStanMCMC(mp) + j_draws <- stanmod$draws() + + summary_quantiles <- function(draws) { + apply( + as.matrix(draws), + 2, + quantile, + probs = c(0.1, 0.5, 0.9), + names = FALSE + ) + } + + b_beta <- b_draws |> + posterior::subset_draws( + variable = c("b_ones", "b_cov1", "b_cov2", "b_cov3") + ) |> + posterior::as_draws_df() + + j_beta <- j_draws |> + posterior::subset_draws(variable = "beta_os_cov") |> + posterior::as_draws_df() + + expect_equal( + summary_quantiles(j_beta), + summary_quantiles(b_beta), + tolerance = 0.1, + ignore_attr = TRUE + ) + + b_c2 <- 2^2 * b_draws$hs_slab + b_shrinkage <- dplyr::tibble( + cov1 = 1 / (1 + b_draws$hs_global^2 * b_draws$`hs_local[1]`^2 / b_c2), + cov2 = 1 / (1 + b_draws$hs_global^2 * b_draws$`hs_local[2]`^2 / b_c2), + cov3 = 1 / (1 + b_draws$hs_global^2 * b_draws$`hs_local[3]`^2 / b_c2) + ) + + j_shrinkage <- shrinkage(mp) + + expect_equal( + summary_quantiles(j_shrinkage), + summary_quantiles(b_shrinkage), + tolerance = 0.1, + ignore_attr = TRUE + ) +}) From bd6099a11b8035fb68e3ee304bd792d18ff1782e Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 7 Jul 2026 22:04:03 +0800 Subject: [PATCH 10/30] fix model.matrix for models without intercept --- R/DataSurvival.R | 4 +++- tests/testthat/test-DataSurvival.R | 23 +++++++++++++++++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/R/DataSurvival.R b/R/DataSurvival.R index 32c0916c..0640821c 100644 --- a/R/DataSurvival.R +++ b/R/DataSurvival.R @@ -128,7 +128,9 @@ model.matrix.DataSurvival <- function( ) { 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] + if (length(remove_index) > 0) { + design_mat <- design_mat[, -remove_index, drop = FALSE] + } rownames(design_mat) <- NULL design_mat } diff --git a/tests/testthat/test-DataSurvival.R b/tests/testthat/test-DataSurvival.R index 807400cb..cbfbd3e9 100644 --- a/tests/testthat/test-DataSurvival.R +++ b/tests/testthat/test-DataSurvival.R @@ -124,6 +124,29 @@ test_that("model.matrix for DataSurvival works as expected", { expect_snapshot_value(res, tolerance = 1e-2, style = "deparse") }) +test_that("model.matrix for DataSurvival preserves no-intercept formulas", { + x <- data.frame( + vtime = c(10, 20, 30, 25, 15), + vevent = c(1, 1, 0, 1, 0), + ones = 1, + vcov1 = c(-1, -0.5, 0, 0.5, 1), + vcov2 = c(1, 2, 3, 4, 5) + ) + + obj <- DataSurvival( + data = x, + formula = Surv(vtime, vevent) ~ 0 + ones + vcov1 + vcov2 + ) + + expected <- model.matrix(~ 0 + ones + vcov1 + vcov2, data = x) + rownames(expected) <- NULL + + res <- as_stan_list(obj) + expect_equal(res$p_os_cov_design, 3) + expect_equal(res$os_cov_design, expected) + expect_equal(covariates(obj), c("ones", "vcov1", "vcov2")) +}) + test_that("mirror_design_matrix() works as expected", { set.seed(3102) N <- 50 From 0ad5437312af26f3a331ed24be2a8b67735caeb2 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 8 Jul 2026 11:37:17 +0800 Subject: [PATCH 11/30] set lambda effectively to 1 in jmpost comparison with brms --- tests/testthat/test-brms_compare.R | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/tests/testthat/test-brms_compare.R b/tests/testthat/test-brms_compare.R index 67b8dbc8..68722776 100644 --- a/tests/testthat/test-brms_compare.R +++ b/tests/testthat/test-brms_compare.R @@ -293,8 +293,8 @@ test_that("jmpost and brms get similar horseshoe estimates for survival models", "horseshoe(df = 1, df_global = 1, df_slab = 4, scale_global = 0.3, scale_slab = 2)", class = "b" ), - warmup = 5000, - iter = 10000, + warmup = 50, + iter = 100, chains = 4, cores = 4, seed = 9826, @@ -326,6 +326,9 @@ test_that("jmpost and brms get similar horseshoe estimates for survival models", jm <- JointModel( survival = SurvivalExponential( + # Let's use an approximate point mass + # at lambda = 1 to get the same model as brms. + lambda = prior_normal(1, 0.00001), beta = prior_horseshoe( df = 1, df_global = 1, @@ -339,8 +342,8 @@ test_that("jmpost and brms get similar horseshoe estimates for survival models", mp <- sampleStanModel( jm, data = jdat, - iter_warmup = 5000, - iter_sampling = 5000, + iter_warmup = 50, + iter_sampling = 50, chains = 4, refresh = 200, parallel_chains = 4, From 4654f95cab7cf72149f978ea3624bfabb34f0f5f Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 8 Jul 2026 12:09:01 +0800 Subject: [PATCH 12/30] truncation, functions --- R/Prior.R | 17 +++++++++++++---- inst/stan/base/functions.stan | 16 ++++++++++++++++ 2 files changed, 29 insertions(+), 4 deletions(-) diff --git a/R/Prior.R b/R/Prior.R index f9ce1315..64cd2dd7 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -491,8 +491,9 @@ prior_horseshoe <- function( "scale_slab = {scale_slab})" ), repr_model = paste( - "prior_local_{name} ~ student_t(prior_df_{name}, 0, 1);", - "prior_global_{name} ~ student_t(prior_df_global_{name}, 0, prior_scale_global_{name});", + # Important: T[0, ] for truncating to positive values! + "prior_local_{name} ~ student_t(prior_df_{name}, 0, 1) T[0, ];", + "prior_global_{name} ~ student_t(prior_df_global_{name}, 0, prior_scale_global_{name}) T[0, ];", "prior_slab_{name} ~ inv_gamma(prior_df_slab_{name} / 2, prior_df_slab_{name} / 2);", "{name} ~ normal(rep_vector(0, {size}), prior_scales_{name})", sep = "\n " @@ -511,8 +512,14 @@ prior_horseshoe <- function( ), repr_transformed_parameters = c( "real prior_c2_{name} = square(prior_scale_slab_{name}) * prior_slab_{name};", - "vector[{size}] prior_shrinkage_factors_{name} = rep_vector(1, {size}) ./ (rep_vector(1, {size}) + square(prior_global_{name}) * square(prior_local_{name}) / prior_c2_{name});", - "vector[{size}] prior_scales_{name} = prior_global_{name} * sqrt((prior_c2_{name} * square(prior_local_{name})) ./ (prior_c2_{name} + square(prior_global_{name}) * square(prior_local_{name})));" + paste( + "vector[{size}] prior_shrinkage_factors_{name} =", + "shrinkage_horseshoe(prior_local_{name}, prior_global_{name}, prior_c2_{name});" + ), + paste( + "vector[{size}] prior_scales_{name} =", + "scales_horseshoe(prior_local_{name}, prior_global_{name}, prior_c2_{name});" + ) ), centre = 0, sample = \(n) { @@ -949,6 +956,8 @@ local_rnorm_vector <- \(n, mus, sigmas) { #' @rdname Local_Sample local_rhorseshoe <- \(n, df, df_global, df_slab, scale_global, scale_slab) { + # Note: Half-t distribution is here the same as truncated t + # distribution because the t distribution is symmetric around 0. local <- abs(local_rt(n, df, 0, 1)) global <- abs(local_rt(n, df_global, 0, scale_global)) slab <- local_rinvgamma(n, df_slab / 2, df_slab / 2) diff --git a/inst/stan/base/functions.stan b/inst/stan/base/functions.stan index 39528bc0..d3f91411 100644 --- a/inst/stan/base/functions.stan +++ b/inst/stan/base/functions.stan @@ -130,5 +130,21 @@ functions { } return result; } + + // Efficient computation of the horseshoe scale parameters. + vector scales_horseshoe(vector lambda, real tau, real c2) { + int K = rows(lambda); + vector[K] lambda2 = square(lambda); + vector[K] lambda_tilde = sqrt(c2 * lambda2 ./ (c2 + tau^2 * lambda2)); + return lambda_tilde * tau; + } + + // Similar for shrinkage factors calculation. + vector shrinkage_horseshoe(vector lambda, real tau, real c2) { + int K = rows(lambda); + vector[K] lambda2 = square(lambda); + vector[K] shrinkage = 1 ./ (1 + (tau^2 * lambda2) ./ c2); + return shrinkage; + } } From 59d5a4deef9717c5c8fb7f0e51b0f5adab77b81d Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 8 Jul 2026 12:11:52 +0800 Subject: [PATCH 13/30] move shrinkage factor calc to generated quantities to save sampling time --- R/Prior.R | 19 +++++++++++++++---- inst/stan/base/base.stan | 2 +- man/Prior-class.Rd | 5 +++++ tests/testthat/test-Prior.R | 11 +++++++---- 4 files changed, 28 insertions(+), 9 deletions(-) diff --git a/R/Prior.R b/R/Prior.R index 64cd2dd7..b645ce0c 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -33,6 +33,7 @@ NULL #' @slot repr_data (`string`)\cr See arguments. #' @slot repr_parameters (`string`)\cr See arguments. #' @slot repr_transformed_parameters (`string`)\cr See arguments. +#' @slot repr_generated_quantities (`string`)\cr See arguments. #' @slot centre (`numeric`)\cr See arguments. #' @slot validation (`list`)\cr See arguments. #' @slot display (`string`)\cr See arguments. @@ -52,6 +53,7 @@ NULL "repr_data" = "character", "repr_parameters" = "character", "repr_transformed_parameters" = "character", + "repr_generated_quantities" = "character", "centre" = "numeric", "validation" = "list", "sample" = "function", @@ -71,6 +73,8 @@ NULL #' the Stan code representation for the parameters block. #' @typed repr_transformed_parameters: string #' the Stan code representation for the transformed parameters block. +#' @typed repr_generated_quantities: string +#' the Stan code representation for the generated quantities block. #' @typed display: string #' the string to display when object is printed. #' @typed centre: numeric @@ -95,6 +99,7 @@ Prior <- function( sample, repr_parameters = "", repr_transformed_parameters = "", + repr_generated_quantities = "", limits = c(-Inf, Inf), .allow_vectors = FALSE ) { @@ -104,6 +109,7 @@ Prior <- function( repr_data = repr_data, repr_parameters = repr_parameters, repr_transformed_parameters = repr_transformed_parameters, + repr_generated_quantities = repr_generated_quantities, centre = centre, display = display, validation = validation, @@ -286,6 +292,9 @@ as.StanModule.Prior <- function(object, name, size = 1, ...) { "model {{", indent_stan(paste0(object@repr_model, trunctation)), "}}", + "generated quantities {{", + indent_stan(object@repr_generated_quantities), + "}}", sep = "\n" ) StanModule(glue::glue(string, name = name, size = size)) @@ -512,15 +521,17 @@ prior_horseshoe <- function( ), repr_transformed_parameters = c( "real prior_c2_{name} = square(prior_scale_slab_{name}) * prior_slab_{name};", - paste( - "vector[{size}] prior_shrinkage_factors_{name} =", - "shrinkage_horseshoe(prior_local_{name}, prior_global_{name}, prior_c2_{name});" - ), paste( "vector[{size}] prior_scales_{name} =", "scales_horseshoe(prior_local_{name}, prior_global_{name}, prior_c2_{name});" ) ), + repr_generated_quantities = c( + paste( + "vector[{size}] prior_shrinkage_factors_{name} =", + "shrinkage_horseshoe(prior_local_{name}, prior_global_{name}, prior_c2_{name});" + ) + ), centre = 0, sample = \(n) { local_rhorseshoe( diff --git a/inst/stan/base/base.stan b/inst/stan/base/base.stan index d611b966..f949f0bc 100755 --- a/inst/stan/base/base.stan +++ b/inst/stan/base/base.stan @@ -71,5 +71,5 @@ generated quantities{ {{ longitudinal.generated_quantities }} {{ link.generated_quantities }} {{ survival.generated_quantities }} +{{ priors.generated_quantities }} } - diff --git a/man/Prior-class.Rd b/man/Prior-class.Rd index 9b964890..6ee77f5b 100644 --- a/man/Prior-class.Rd +++ b/man/Prior-class.Rd @@ -17,6 +17,7 @@ Prior( sample, repr_parameters = "", repr_transformed_parameters = "", + repr_generated_quantities = "", limits = c(-Inf, Inf), .allow_vectors = FALSE ) @@ -41,6 +42,8 @@ the same names as the \code{paramaters} slot.} \item{repr_transformed_parameters}{(\code{string}) the Stan code representation for the transformed parameters block.} +\item{repr_generated_quantities}{(\code{string}) the Stan code representation for the generated quantities block.} + \item{limits}{(\code{numeric}) the lower and upper limits for a truncated distribution} \item{.allow_vectors}{(\code{flag}) whether to allow vector parameters.} @@ -61,6 +64,8 @@ Specifies the prior distribution in a Stan Model \item{\code{repr_transformed_parameters}}{(\code{string})\cr See arguments.} +\item{\code{repr_generated_quantities}}{(\code{string})\cr See arguments.} + \item{\code{centre}}{(\code{numeric})\cr See arguments.} \item{\code{validation}}{(\code{list})\cr See arguments.} diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 1caffc41..94908c48 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -357,13 +357,13 @@ test_that("prior_horseshoe works as expected", { expect_true( any(grepl( "vector\\[p\\] prior_shrinkage_factors_beta", - x_stan_module@transformed_parameters + x_stan_module@generated_quantities )) ) expect_true( any(grepl( - "rep_vector(1, p) ./ (rep_vector(1, p) + square(prior_global_beta) * square(prior_local_beta) / prior_c2_beta)", - x_stan_module@transformed_parameters, + "shrinkage_horseshoe(prior_local_beta, prior_global_beta, prior_c2_beta)", + x_stan_module@generated_quantities, fixed = TRUE )) ) @@ -386,13 +386,16 @@ test_that("prior_horseshoe works as expected", { ) # Check the model syntax. - header <- StanModule( + header <- merge( + StanModule("base/functions.stan"), + StanModule( "data { int p; } parameters { vector[p] beta; }" + ) ) stan_file <- cmdstanr::write_stan_file( as.character(merge(header, x_stan_module)), From 6fe50966b3f877eecc89b8c02d95cbebd78e8773 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 8 Jul 2026 12:15:55 +0800 Subject: [PATCH 14/30] ensure to return empty stan module for empty prior to avoid warning --- R/Prior.R | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/R/Prior.R b/R/Prior.R index b645ce0c..cddd1701 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -267,6 +267,17 @@ setMethod( #' @family as.StanModule #' @export as.StanModule.Prior <- function(object, name, size = 1, ...) { + stan_repr <- c( + object@repr_data, + object@repr_parameters, + object@repr_transformed_parameters, + object@repr_model, + object@repr_generated_quantities + ) + if (all(nchar(stan_repr) == 0)) { + return(StanModule()) + } + indent_stan <- function(x) { x <- x[nchar(x) >= 1] if (length(x) == 0) { From 8180d0571d3f6247cc31b1fa17ddb5049f5a5984 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 09:31:42 +0800 Subject: [PATCH 15/30] Add prior const (#470) * slight tangent: add prior_const for surv exp first as a try * first matching comparison between jmpost and brms for horseshoe prior usage * first try to reduce to one model.stan file * let Prior handle parameter declaration in the Stan code together with limits * adapt documentation * add snapshot tests for model stan code construction * add vignette * few fixes --- DESCRIPTION | 4 +- NAMESPACE | 2 + NEWS.md | 2 + R/JointModel.R | 2 +- R/LongitudinalClaretBruno.R | 26 +- R/LongitudinalGSF.R | 22 +- R/LongitudinalRandomSlope.R | 4 +- R/LongitudinalSteinFojo.R | 20 +- R/Parameter.R | 116 ++++++- R/ParameterList.R | 13 +- R/Prior.R | 40 ++- R/StanModel.R | 6 +- R/SurvivalExponential.R | 1 + R/SurvivalGamma.R | 4 +- R/SurvivalLoglogistic.R | 3 + R/SurvivalWeibullPH.R | 4 +- R/generics.R | 1 + _pkgdown.yml | 9 +- inst/WORDLIST | 7 + inst/stan/base/base.stan | 2 +- inst/stan/base/link.stan | 11 - inst/stan/base/survival.stan | 5 - inst/stan/lm-claret-bruno/model.stan | 36 --- inst/stan/lm-gsf/model.stan | 37 --- inst/stan/lm-random-slope/model.stan | 14 - inst/stan/lm-stein-fojo/model.stan | 33 -- inst/stan/sm-exponential/model.stan | 7 - inst/stan/sm-gamma/model.stan | 10 - inst/stan/sm-loglogistic/model.stan | 8 - inst/stan/sm-weibull-ph/model.stan | 8 - man/JointModel-class.Rd | 2 +- man/LongitudinalModel-class.Rd | 4 +- man/ParameterList-class.Rd | 3 +- man/Prior-class.Rd | 7 +- man/StanModel-class.Rd | 4 +- man/SurvivalModel-class.Rd | 4 +- man/as.StanModule.ParameterDeclaration.Rd | 19 ++ man/covariates.Rd | 2 +- man/prior_beta.Rd | 1 + man/prior_cauchy.Rd | 1 + man/prior_const.Rd | 32 ++ man/prior_gamma.Rd | 1 + man/prior_horseshoe.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 | 1 + man/prior_std_normal.Rd | 1 + man/prior_student_t.Rd | 1 + man/prior_uniform.Rd | 1 + man/render_stan_const_declaration.Rd | 23 ++ man/render_stan_parameter_declaration.Rd | 22 ++ man/render_stan_parameter_limits.Rd | 19 ++ man/stanmodel_arguments.Rd | 6 +- tests/testthat/_snaps/JointModel.md | 63 ++++ tests/testthat/_snaps/Prior.md | 10 + tests/testthat/_snaps/SurvivalExponential.md | 12 + tests/testthat/models/ParameterList.stan | 6 +- tests/testthat/test-JointModel.R | 63 ++++ tests/testthat/test-Link.R | 7 +- tests/testthat/test-Parameter.R | 62 ++++ tests/testthat/test-Prior.R | 35 ++ tests/testthat/test-SurvivalExponential.R | 49 +++ tests/testthat/test-brms_compare.R | 61 +++- vignettes/covariate-selection-horseshoe.Rmd | 318 +++++++++++++++++++ vignettes/custom-model.Rmd | 64 +++- vignettes/custom-model.stan | 19 -- vignettes/extending-jmpost.Rmd | 40 ++- 71 files changed, 1118 insertions(+), 309 deletions(-) create mode 100644 man/as.StanModule.ParameterDeclaration.Rd create mode 100644 man/prior_const.Rd create mode 100644 man/render_stan_const_declaration.Rd create mode 100644 man/render_stan_parameter_declaration.Rd create mode 100644 man/render_stan_parameter_limits.Rd create mode 100644 vignettes/covariate-selection-horseshoe.Rmd diff --git a/DESCRIPTION b/DESCRIPTION index 22baeab7..617d7867 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -26,7 +26,7 @@ Encoding: UTF-8 Language: en-GB URL: https://genentech.github.io/jmpost/ BugReports: https://github.com/Genentech/jmpost/issues -Roxygen: list(markdown = TRUE, packages = "roxytypes", roclets = c("namespace", "rd", "roxylint::roxylint")) +Roxygen: list(markdown = TRUE, packages = "roxytypes", roclets = c("namespace", "rd")) Depends: R (>= 4.1.0) Imports: @@ -134,4 +134,4 @@ RdMacros: Rdpack LazyData: true Config/roxygen2/version: 8.0.0 Config/Needs/documentation: - roxytypes, roxylint + roxytypes diff --git a/NAMESPACE b/NAMESPACE index 78e44b94..2b3da074 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -21,6 +21,7 @@ S3method(as.StanModule,JointModelSamples) S3method(as.StanModule,Link) S3method(as.StanModule,LinkComponent) S3method(as.StanModule,Parameter) +S3method(as.StanModule,ParameterDeclaration) S3method(as.StanModule,ParameterList) S3method(as.StanModule,Prior) S3method(as.StanModule,PromiseLinkComponent) @@ -269,6 +270,7 @@ export(merge) export(populationHR) export(prior_beta) export(prior_cauchy) +export(prior_const) export(prior_gamma) export(prior_horseshoe) export(prior_init_only) diff --git a/NEWS.md b/NEWS.md index 74bd6916..206641b3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,8 @@ - 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. +- Added `prior_const()` which fixes a parameter at a constant value, i.e. uses a point-mass prior distribution. +- The Stan files associated with models no longer hardcode the parameter declarations. This task is now handled with the `Prior` objects. If a constant value is used for a parameter instead of a prior distribution, then this parameter is declared in the `data` block and reassigned in the `transformed_parameters` block. - 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). diff --git a/R/JointModel.R b/R/JointModel.R index a8ea2e15..74ccfb6d 100755 --- a/R/JointModel.R +++ b/R/JointModel.R @@ -33,7 +33,7 @@ setClassUnion("SurvivalModel_OR_NULL", c("SurvivalModel", "NULL")) #' @slot longitudinal ([`LongitudinalModel`] or `NULL`)\cr the longitudinal model. #' @slot survival ([`SurvivalModel`] or `NULL`)\cr the survival model. #' @slot link (`Link`)\cr the link. -#' @slot parameters (`ParameterList`)\cr the parameter specification. +#' @slot parameters (`ParameterList`)\cr the parameter declaration and prior specification. #' #' @family JointModel #' @export JointModel diff --git a/R/LongitudinalClaretBruno.R b/R/LongitudinalClaretBruno.R index 99e487b5..15e49919 100755 --- a/R/LongitudinalClaretBruno.R +++ b/R/LongitudinalClaretBruno.R @@ -78,11 +78,11 @@ LongitudinalClaretBruno <- function( )) # Apply constraints - omega_b <- set_limits(omega_b, lower = 0) - omega_g <- set_limits(omega_g, lower = 0) - omega_c <- set_limits(omega_c, lower = 0) - omega_p <- set_limits(omega_p, lower = 0) - sigma <- set_limits(sigma, lower = 0) + omega_b <- set_limits(omega_b, lower = getOption("jmpost.double_eps")) + omega_g <- set_limits(omega_g, lower = getOption("jmpost.double_eps")) + omega_c <- set_limits(omega_c, lower = getOption("jmpost.double_eps")) + omega_p <- set_limits(omega_p, lower = getOption("jmpost.double_eps")) + sigma <- set_limits(sigma, lower = getOption("jmpost.double_eps")) parameters <- list( Parameter(name = "lm_clbr_mu_b", prior = mu_b, size = "n_studies"), @@ -107,34 +107,34 @@ LongitudinalClaretBruno <- function( list( Parameter( name = "lm_clbr_ind_b", - prior = prior_init_only(prior_lognormal( + prior = set_limits(prior_init_only(prior_lognormal( median(mu_b), median(omega_b) - )), + )), lower = getOption("jmpost.double_eps")), size = "n_subjects" ), Parameter( name = "lm_clbr_ind_g", - prior = prior_init_only(prior_lognormal( + prior = set_limits(prior_init_only(prior_lognormal( median(mu_g), median(omega_g) - )), + )), lower = getOption("jmpost.double_eps")), size = "n_subjects" ), Parameter( name = "lm_clbr_ind_c", - prior = prior_init_only(prior_lognormal( + prior = set_limits(prior_init_only(prior_lognormal( median(mu_c), median(omega_c) - )), + )), lower = getOption("jmpost.double_eps")), size = "n_subjects" ), Parameter( name = "lm_clbr_ind_p", - prior = prior_init_only(prior_lognormal( + prior = set_limits(prior_init_only(prior_lognormal( median(mu_p), median(omega_p) - )), + )), lower = getOption("jmpost.double_eps")), size = "n_subjects" ) ) diff --git a/R/LongitudinalGSF.R b/R/LongitudinalGSF.R index 83175ae6..46462501 100755 --- a/R/LongitudinalGSF.R +++ b/R/LongitudinalGSF.R @@ -82,11 +82,11 @@ LongitudinalGSF <- function( )) # Apply constraints - omega_bsld <- set_limits(omega_bsld, lower = 0) - omega_ks <- set_limits(omega_ks, lower = 0) - omega_kg <- set_limits(omega_kg, lower = 0) - omega_phi <- set_limits(omega_phi, lower = 0) - sigma <- set_limits(sigma, lower = 0) + omega_bsld <- set_limits(omega_bsld, lower = getOption("jmpost.double_eps")) + omega_ks <- set_limits(omega_ks, lower = getOption("jmpost.double_eps")) + omega_kg <- set_limits(omega_kg, lower = getOption("jmpost.double_eps")) + omega_phi <- set_limits(omega_phi, lower = getOption("jmpost.double_eps")) + sigma <- set_limits(sigma, lower = getOption("jmpost.double_eps")) parameters <- list( Parameter(name = "lm_gsf_mu_bsld", prior = mu_bsld, size = "n_studies"), @@ -115,26 +115,26 @@ LongitudinalGSF <- function( list( Parameter( name = "lm_gsf_psi_bsld", - prior = prior_init_only(prior_lognormal( + prior = set_limits(prior_init_only(prior_lognormal( median(mu_bsld), median(omega_bsld) - )), + )), lower = getOption("jmpost.double_eps")), size = "n_subjects" ), Parameter( name = "lm_gsf_psi_ks", - prior = prior_init_only(prior_lognormal( + prior = set_limits(prior_init_only(prior_lognormal( median(mu_ks), median(omega_ks) - )), + )), lower = getOption("jmpost.double_eps")), size = "n_subjects" ), Parameter( name = "lm_gsf_psi_kg", - prior = prior_init_only(prior_lognormal( + prior = set_limits(prior_init_only(prior_lognormal( median(mu_kg), median(omega_kg) - )), + )), lower = getOption("jmpost.double_eps")), size = "n_subjects" ), Parameter( diff --git a/R/LongitudinalRandomSlope.R b/R/LongitudinalRandomSlope.R index 67f74d21..2964d91a 100755 --- a/R/LongitudinalRandomSlope.R +++ b/R/LongitudinalRandomSlope.R @@ -50,8 +50,8 @@ LongitudinalRandomSlope <- function( )) # Apply constriants - sigma <- set_limits(sigma, lower = 0) - slope_sigma <- set_limits(slope_sigma, lower = 0) + sigma <- set_limits(sigma, lower = getOption("jmpost.double_eps")) + slope_sigma <- set_limits(slope_sigma, lower = getOption("jmpost.double_eps")) .LongitudinalRandomSlope( LongitudinalModel( diff --git a/R/LongitudinalSteinFojo.R b/R/LongitudinalSteinFojo.R index 8ae2b39a..361c830c 100755 --- a/R/LongitudinalSteinFojo.R +++ b/R/LongitudinalSteinFojo.R @@ -72,10 +72,10 @@ LongitudinalSteinFojo <- function( )) # Apply constriants - omega_bsld <- set_limits(omega_bsld, lower = 0) - omega_ks <- set_limits(omega_ks, lower = 0) - omega_kg <- set_limits(omega_kg, lower = 0) - sigma <- set_limits(sigma, lower = 0) + omega_bsld <- set_limits(omega_bsld, lower = getOption("jmpost.double_eps")) + omega_ks <- set_limits(omega_ks, lower = getOption("jmpost.double_eps")) + omega_kg <- set_limits(omega_kg, lower = getOption("jmpost.double_eps")) + sigma <- set_limits(sigma, lower = getOption("jmpost.double_eps")) parameters <- list( Parameter(name = "lm_sf_mu_bsld", prior = mu_bsld, size = "n_studies"), @@ -98,26 +98,26 @@ LongitudinalSteinFojo <- function( list( Parameter( name = "lm_sf_psi_bsld", - prior = prior_init_only(prior_lognormal( + prior = set_limits(prior_init_only(prior_lognormal( median(mu_bsld), median(omega_bsld) - )), + )), lower = getOption("jmpost.double_eps")), size = "n_subjects" ), Parameter( name = "lm_sf_psi_ks", - prior = prior_init_only(prior_lognormal( + prior = set_limits(prior_init_only(prior_lognormal( median(mu_ks), median(omega_ks) - )), + )), lower = getOption("jmpost.double_eps")), size = "n_subjects" ), Parameter( name = "lm_sf_psi_kg", - prior = prior_init_only(prior_lognormal( + prior = set_limits(prior_init_only(prior_lognormal( median(mu_kg), median(omega_kg) - )), + )), lower = getOption("jmpost.double_eps")), size = "n_subjects" ) ) diff --git a/R/Parameter.R b/R/Parameter.R index 73e12b9f..9618191d 100755 --- a/R/Parameter.R +++ b/R/Parameter.R @@ -86,7 +86,118 @@ setValidity( #' @family as.StanModule #' @export as.StanModule.Parameter <- function(object, ...) { - as.StanModule(object@prior, name = object@name, size = object@size) + merge( + as.StanModule.ParameterDeclaration(object), + as.StanModule(object@prior, name = object@name, size = object@size) + ) +} + +#' Render Stan Parameter Constraints +#' +#' Converts a lower/upper limit vector from a [`Prior`] into a Stan declaration +#' constraint, e.g. ``. +#' +#' @typed limits: numeric +#' length-two vector containing lower and upper parameter limits. +#' +#' @return A length-one character vector. +#' +#' @keywords internal +render_stan_parameter_limits <- function(limits) { + limits_names <- c("lower", "upper") + limits <- stats::setNames(limits, limits_names) + limits <- limits[is.finite(limits)] + if (length(limits) == 0) { + return("") + } + constraints <- paste(names(limits), limits, sep = "=", collapse = ", ") + paste0("<", constraints, ">") +} + +#' Render Stan Parameter Declaration +#' +#' Creates a Stan declaration for a sampled scalar or vector parameter. +#' +#' @typed name: string +#' parameter name. +#' @typed size: numeric_OR_character +#' parameter size. +#' @typed limits: numeric +#' lower and upper parameter limits. +#' +#' @return A length-one character vector. +#' +#' @keywords internal +render_stan_parameter_declaration <- function(name, size, limits) { + constraints <- render_stan_parameter_limits(limits) + if (length(size) == 1 && is.numeric(size) && size == 1) { + return(glue::glue("real{constraints} {name};")) + } + glue::glue("vector{constraints}[{size}] {name};") +} + +#' Render Stan Constant Parameter Declaration +#' +#' Creates a Stan transformed-parameter declaration that fixes a scalar or vector +#' parameter at a data-supplied constant value. +#' +#' @typed name: string +#' parameter name. +#' @typed size: numeric_OR_character +#' parameter size. +#' @typed limits: numeric +#' lower and upper parameter limits. +#' +#' @return A length-one character vector. +#' +#' @keywords internal +render_stan_const_declaration <- function(name, size, limits) { + constraints <- render_stan_parameter_limits(limits) + value <- glue::glue("prior_const_{name}") + if (length(size) == 1 && is.numeric(size) && size == 1) { + return(glue::glue("real{constraints} {name} = {value};")) + } + glue::glue( + "vector{constraints}[{size}] {name} = rep_vector({value}, {size});" + ) +} + +#' `Parameter` Declaration -> `StanModule` +#' +#' Creates only the Stan declaration for a [`Parameter`]. Sampled parameters are +#' declared in the `parameters` block; constant parameters are declared in the +#' `transformed parameters` block. +#' +#' @inheritParams Parameter-Shared +#' +#' @return A [`StanModule`] object. +#' +#' @export +as.StanModule.ParameterDeclaration <- function(object) { + declaration <- if (object@prior@.is_const) { + render_stan_const_declaration( + name = object@name, + size = object@size, + limits = object@prior@limits + ) + } else { + render_stan_parameter_declaration( + name = object@name, + size = object@size, + limits = object@prior@limits + ) + } + + block <- if (object@prior@.is_const) { + "transformed parameters" + } else { + "parameters" + } + StanModule(glue::glue( + "{block} {{ + {declaration} +}}" + )) } @@ -139,6 +250,9 @@ size.Parameter <- function(object) object@size #' @family Parameter #' @export as.character.Parameter <- function(x, ...) { + if (x@prior@.is_const) { + return(paste0(x@name, " = ", as.character(x@prior))) + } paste0(x@name, " ~ ", as.character(x@prior)) } diff --git a/R/ParameterList.R b/R/ParameterList.R index 2e6e0b54..8f30fd9b 100644 --- a/R/ParameterList.R +++ b/R/ParameterList.R @@ -20,7 +20,8 @@ NULL #' `ParameterList` #' #' This class extends the general [`list`] type for containing [`Parameter`] -#' specifications. +#' specifications. When converted to a [`StanModule`], the list supplies the +#' Stan declarations and prior statements for those parameters. #' #' @@ -151,14 +152,15 @@ names.ParameterList <- function(x) { #' @describeIn ParameterList-Getter-Methods The parameter-list's parameter initial values #' @export initialValues.ParameterList <- function(object, n_chains, ...) { + parameters <- Filter(\(x) !x@prior@.is_const, object@parameters) # Generate initial values as a list of lists. This is to ensure it is in the required # format as specified by cmdstanr see the `init` argument of # `help("model-method-sample", "cmdstanr")` for more details lapply( seq_len(n_chains), \(i) { - vals <- lapply(object@parameters, initialValues) - name <- vapply(object@parameters, names, character(1)) + vals <- lapply(parameters, initialValues) + name <- vapply(parameters, names, character(1)) names(vals) <- name vals } @@ -169,8 +171,9 @@ initialValues.ParameterList <- function(object, n_chains, ...) { #' @describeIn ParameterList-Getter-Methods The parameter-list's parameter dimensionality #' @export size.ParameterList <- function(object) { - x <- lapply(object@parameters, size) - names(x) <- names(object) + parameters <- Filter(\(x) !x@prior@.is_const, object@parameters) + x <- lapply(parameters, size) + names(x) <- vapply(parameters, names, character(1)) return(x) } diff --git a/R/Prior.R b/R/Prior.R index cddd1701..26bf86b5 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -40,6 +40,7 @@ NULL #' @slot sample (`function`)\cr See arguments. #' @slot limits (`numeric`)\cr See arguments. #' @slot .allow_vectors (`logical`)\cr See arguments. +#' @slot .is_const (`logical`)\cr See arguments. #' #' @family Prior-internal #' @export Prior @@ -58,7 +59,8 @@ NULL "validation" = "list", "sample" = "function", "limits" = "numeric", - ".allow_vectors" = "logical" + ".allow_vectors" = "logical", + ".is_const" = "logical" ) ) @@ -88,6 +90,8 @@ NULL #' the lower and upper limits for a truncated distribution #' @typed .allow_vectors: flag #' whether to allow vector parameters. +#' @typed .is_const: flag +#' whether this prior fixes the parameter at a constant value. #' @rdname Prior-class Prior <- function( parameters, @@ -101,7 +105,8 @@ Prior <- function( repr_transformed_parameters = "", repr_generated_quantities = "", limits = c(-Inf, Inf), - .allow_vectors = FALSE + .allow_vectors = FALSE, + .is_const = FALSE ) { .Prior( parameters = parameters, @@ -115,7 +120,8 @@ Prior <- function( validation = validation, sample = sample, limits = limits, - .allow_vectors = .allow_vectors + .allow_vectors = .allow_vectors, + .is_const = .is_const ) } @@ -324,6 +330,10 @@ as.StanModule.Prior <- function(object, name, size = 1, ...) { #' @export as_stan_list.Prior <- function(object, name, ...) { vals <- object@parameters + if (object@.is_const) { + names(vals) <- paste0("prior_const_", name) + return(vals) + } vals_names <- names(vals) if (length(vals_names) >= 1) { names(vals) <- paste0("prior_", vals_names, "_", name) @@ -582,6 +592,27 @@ prior_std_normal <- function() { ) } +#' Constant Prior Distribution +#' +#' @typed value: number +#' the fixed parameter value. +#' @family Prior +#' @export +prior_const <- function(value) { + Prior( + parameters = list(value = value), + display = "const(value = {value})", + repr_model = "", + repr_data = "real prior_const_{name};", + centre = value, + sample = \(n) rep(value, n), + validation = list( + value = is.numeric + ), + .is_const = TRUE + ) +} + #' Cauchy Prior Distribution #' #' @typed mu: number @@ -712,7 +743,8 @@ prior_init_only <- function(dist) { dist@sample(n) }, centre = dist@centre, - validation = list() + validation = list(), + limits = dist@limits ) } diff --git a/R/StanModel.R b/R/StanModel.R index 957d7d2d..b17696d4 100644 --- a/R/StanModel.R +++ b/R/StanModel.R @@ -9,11 +9,11 @@ NULL #' [StanModel()]. #' #' @typed stan: StanModule -#' code containing the Stan code specification. +#' code containing the model-specific Stan code specification. #' @typed parameters: ParameterList -#' the parameter specification. +#' the parameter declaration and prior specification. #' @typed parameter: ParameterList -#' the (single) parameter specification. +#' the (single) parameter declaration and prior specification. #' @typed name: character #' display name for the model object. #' @param ... additional arguments for [StanModel()]. diff --git a/R/SurvivalExponential.R b/R/SurvivalExponential.R index 396468a2..e63e877c 100755 --- a/R/SurvivalExponential.R +++ b/R/SurvivalExponential.R @@ -28,6 +28,7 @@ SurvivalExponential <- function( lambda = prior_gamma(2, 5), beta = prior_normal(0, 2) ) { + lambda <- set_limits(lambda, lower = getOption("jmpost.double_eps")) .SurvivalExponential( SurvivalModel( name = "Exponential", diff --git a/R/SurvivalGamma.R b/R/SurvivalGamma.R index d2cc2014..efcfef4d 100755 --- a/R/SurvivalGamma.R +++ b/R/SurvivalGamma.R @@ -30,8 +30,8 @@ SurvivalGamma <- function( theta = prior_gamma(2, 0.5), beta = prior_normal(0, 2) ) { - k <- set_limits(k, lower = 0) - theta <- set_limits(theta, lower = 0) + k <- set_limits(k, lower = getOption("jmpost.double_eps")) + theta <- set_limits(theta, lower = getOption("jmpost.double_eps")) .SurvivalGamma( SurvivalModel( diff --git a/R/SurvivalLoglogistic.R b/R/SurvivalLoglogistic.R index 721008a9..d7f1c7ae 100644 --- a/R/SurvivalLoglogistic.R +++ b/R/SurvivalLoglogistic.R @@ -31,6 +31,9 @@ SurvivalLogLogistic <- function( b = prior_gamma(2, 5), beta = prior_normal(0, 2) ) { + a <- set_limits(a, lower = getOption("jmpost.double_eps")) + b <- set_limits(b, lower = getOption("jmpost.double_eps")) + .SurvivalLogLogistic( SurvivalModel( name = "Log-Logistic", diff --git a/R/SurvivalWeibullPH.R b/R/SurvivalWeibullPH.R index b47562a7..b2bfe5cc 100755 --- a/R/SurvivalWeibullPH.R +++ b/R/SurvivalWeibullPH.R @@ -31,8 +31,8 @@ SurvivalWeibullPH <- function( gamma = prior_gamma(2, 0.5), beta = prior_normal(0, 2) ) { - lambda <- set_limits(lambda, lower = 0) - gamma <- set_limits(gamma, lower = 0) + lambda <- set_limits(lambda, lower = getOption("jmpost.double_eps")) + gamma <- set_limits(gamma, lower = getOption("jmpost.double_eps")) .SurvivalWeibullPH( SurvivalModel( diff --git a/R/generics.R b/R/generics.R index f11fad8b..3a7f26eb 100755 --- a/R/generics.R +++ b/R/generics.R @@ -531,6 +531,7 @@ saveObject <- function(object, file, ...) { #' Extract Covariate Names #' #' @typed object: ANY +#' the object to extract covariate names from. #' #' @family covariates #' @export diff --git a/_pkgdown.yml b/_pkgdown.yml index 0c46ed61..2b409885 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -58,6 +58,7 @@ reference: - Prior - prior_beta - prior_cauchy + - prior_const - prior_gamma - prior_invgamma - prior_logistic @@ -65,6 +66,7 @@ reference: - prior_lognormal - prior_normal - prior_normal_vector + - prior_horseshoe - prior_uniform - prior_std_normal - prior_student_t @@ -124,6 +126,8 @@ reference: - brierScore - brierScore.SurvivalQuantities - populationHR + - shrinkage + - shrinkage.JointModelSamples - title: Stan Code contents: @@ -149,6 +153,7 @@ reference: - as.CmdStanMCMC.JointModelSamples - as.character.JointModel - as.StanModule.Parameter + - as.StanModule.ParameterDeclaration - as.StanModule.ParameterList - as.StanModule.Prior - as.StanModule.JointModel @@ -185,6 +190,8 @@ reference: - hazardWindows - sampleObservations - sampleSubjects + - covariates + - covariates.DataSurvival - enableLink - enableGQ - as_formula @@ -218,7 +225,7 @@ articles: contents: - quickstart - '`statistical-specification`' + - '`covariate-selection-horseshoe`' - '`custom-model`' - '`extending-jmpost`' - '`simulating-data`' - diff --git a/inst/WORDLIST b/inst/WORDLIST index 599eeebd..9c20732b 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -21,6 +21,7 @@ Grajek HS Hoffmann IPCW +Inv JointModel JointModelSamples Kerioui @@ -63,6 +64,7 @@ coeficient csv dSLD docetaxel +dotsc du etc extractVariableNames @@ -71,6 +73,7 @@ geq gk groupwise gt +hardcode hardcoded ie ij @@ -80,16 +83,20 @@ ln lt mathbb nonlinear +operatorname parameterise parameterised pk prior's +qquad saveObject sim sk sld stan +standalone throughs tumor u +unregularized urothelial diff --git a/inst/stan/base/base.stan b/inst/stan/base/base.stan index f949f0bc..34da7415 100755 --- a/inst/stan/base/base.stan +++ b/inst/stan/base/base.stan @@ -44,10 +44,10 @@ parameters{ transformed parameters{ +{{ priors.transformed_parameters }} {{ longitudinal.transformed_parameters }} {{ link.transformed_parameters }} {{ survival.transformed_parameters }} -{{ priors.transformed_parameters }} // // Source - base/base.stan diff --git a/inst/stan/base/link.stan b/inst/stan/base/link.stan index 967bf417..73bde668 100644 --- a/inst/stan/base/link.stan +++ b/inst/stan/base/link.stan @@ -24,20 +24,9 @@ functions { } } -parameters{ - // - // Source - base/link.stan - // - {% for item in items -%} - real {{ item }}; - {% endfor -%} - -} - transformed parameters { vector[{{ length(items) }}] link_coefficients; {% for item in items -%} link_coefficients[{{loop.index1}}] = {{ item }}; {% endfor -%} } - diff --git a/inst/stan/base/survival.stan b/inst/stan/base/survival.stan index f74f5f4a..dd4eeb4f 100755 --- a/inst/stan/base/survival.stan +++ b/inst/stan/base/survival.stan @@ -120,11 +120,6 @@ transformed data { parameters { - // - // Source - base/survival.stan - // - // Covariate coefficients. - vector[p_os_cov_design] beta_os_cov; {{ stan.parameters }} } diff --git a/inst/stan/lm-claret-bruno/model.stan b/inst/stan/lm-claret-bruno/model.stan index ce371a16..01e2f150 100755 --- a/inst/stan/lm-claret-bruno/model.stan +++ b/inst/stan/lm-claret-bruno/model.stan @@ -2,42 +2,6 @@ -parameters{ - // - // Source - lm-claret-bruno/model.stan - // - - vector[n_studies] lm_clbr_mu_b; - vector[n_arms] lm_clbr_mu_g; - vector[n_arms] lm_clbr_mu_c; - vector[n_arms] lm_clbr_mu_p; - - vector[n_studies] lm_clbr_omega_b; - vector[n_arms] lm_clbr_omega_g; - vector[n_arms] lm_clbr_omega_c; - vector[n_arms] lm_clbr_omega_p; - -{% if centred -%} - vector[n_subjects] lm_clbr_ind_b; - vector[n_subjects] lm_clbr_ind_g; - vector[n_subjects] lm_clbr_ind_c; - vector[n_subjects] lm_clbr_ind_p; -{% else -%} - vector[n_subjects] lm_clbr_eta_b; - vector[n_subjects] lm_clbr_eta_g; - vector[n_subjects] lm_clbr_eta_c; - vector[n_subjects] lm_clbr_eta_p; -{%- endif -%} - - // Standard deviation of the error term - real lm_clbr_sigma; - -} - - - - - transformed parameters{ // // Source - lm-claret-bruno/model.stan diff --git a/inst/stan/lm-gsf/model.stan b/inst/stan/lm-gsf/model.stan index f05e7045..a7179947 100755 --- a/inst/stan/lm-gsf/model.stan +++ b/inst/stan/lm-gsf/model.stan @@ -2,42 +2,6 @@ -parameters{ - // - // Source - lm-gsf/model.stan - // - - vector[n_studies] lm_gsf_mu_bsld; - vector[n_arms] lm_gsf_mu_ks; - vector[n_arms] lm_gsf_mu_kg; - vector[n_arms] lm_gsf_mu_phi; - - vector[n_studies] lm_gsf_omega_bsld; - vector[n_arms] lm_gsf_omega_ks; - vector[n_arms] lm_gsf_omega_kg; - vector[n_arms] lm_gsf_omega_phi; - -{% if centred -%} - vector[n_subjects] lm_gsf_psi_bsld; - vector[n_subjects] lm_gsf_psi_ks; - vector[n_subjects] lm_gsf_psi_kg; - vector[n_subjects] lm_gsf_psi_phi_logit; -{% else -%} - vector[n_subjects] lm_gsf_eta_tilde_bsld; - vector[n_subjects] lm_gsf_eta_tilde_ks; - vector[n_subjects] lm_gsf_eta_tilde_kg; - vector[n_subjects] lm_gsf_eta_tilde_phi; -{%- endif -%} - - // Standard deviation of the error term - real lm_gsf_sigma; - -} - - - - - transformed parameters{ // // Source - lm-gsf/model.stan @@ -107,4 +71,3 @@ model { lm_gsf_psi_phi_logit ~ normal(lm_gsf_mu_phi[subject_arm_index], lm_gsf_omega_phi[subject_arm_index]); {%- endif -%} } - diff --git a/inst/stan/lm-random-slope/model.stan b/inst/stan/lm-random-slope/model.stan index 11b79b6b..1da9cf3b 100755 --- a/inst/stan/lm-random-slope/model.stan +++ b/inst/stan/lm-random-slope/model.stan @@ -1,19 +1,6 @@ -parameters { - // - // Source - lm-random-slope/model.stan - // - - array [n_studies] real lm_rs_intercept; - array [n_arms] real lm_rs_slope_mu; - vector[n_arms] lm_rs_slope_sigma; - real lm_rs_sigma; - vector[n_subjects] lm_rs_ind_rnd_slope; -} - - transformed parameters { // // Source - lm-random-slope/model.stan @@ -57,4 +44,3 @@ model { ); } - diff --git a/inst/stan/lm-stein-fojo/model.stan b/inst/stan/lm-stein-fojo/model.stan index 889c8760..a0bf5d82 100755 --- a/inst/stan/lm-stein-fojo/model.stan +++ b/inst/stan/lm-stein-fojo/model.stan @@ -2,38 +2,6 @@ -parameters{ - // - // Source - lm-stein-fojo/model.stan - // - - vector[n_studies] lm_sf_mu_bsld; - vector[n_arms] lm_sf_mu_ks; - vector[n_arms] lm_sf_mu_kg; - - vector[n_studies] lm_sf_omega_bsld; - vector[n_arms] lm_sf_omega_ks; - vector[n_arms] lm_sf_omega_kg; - -{% if centred -%} - vector[n_subjects] lm_sf_psi_bsld; - vector[n_subjects] lm_sf_psi_ks; - vector[n_subjects] lm_sf_psi_kg; -{% else -%} - vector[n_subjects] lm_sf_eta_tilde_bsld; - vector[n_subjects] lm_sf_eta_tilde_ks; - vector[n_subjects] lm_sf_eta_tilde_kg; -{%- endif -%} - - // Standard deviation of the error term - real lm_sf_sigma; - -} - - - - - transformed parameters{ // // Source - lm-stein-fojo/model.stan @@ -93,4 +61,3 @@ model { lm_sf_psi_kg ~ lognormal(lm_sf_mu_kg[subject_arm_index], lm_sf_omega_kg[subject_arm_index]); {%- endif -%} } - diff --git a/inst/stan/sm-exponential/model.stan b/inst/stan/sm-exponential/model.stan index e551674d..d5e0e650 100755 --- a/inst/stan/sm-exponential/model.stan +++ b/inst/stan/sm-exponential/model.stan @@ -15,14 +15,7 @@ functions { } -parameters { - // SurvivalExponential - real sm_exp_lambda; -} - - transformed parameters { // SurvivalExponential vector[1] pars_os = [sm_exp_lambda]'; } - diff --git a/inst/stan/sm-gamma/model.stan b/inst/stan/sm-gamma/model.stan index 5d9e6eb1..f28b43e1 100755 --- a/inst/stan/sm-gamma/model.stan +++ b/inst/stan/sm-gamma/model.stan @@ -15,19 +15,9 @@ functions { } -parameters { - // - // Source - sm-gamma/model.stan - // - real sm_gamma_k; - real sm_gamma_theta; -} - - transformed parameters { // // Source - sm-gamma/model.stan // vector[2] pars_os = [sm_gamma_k, sm_gamma_theta]'; } - diff --git a/inst/stan/sm-loglogistic/model.stan b/inst/stan/sm-loglogistic/model.stan index e95f9a0e..5b89683b 100755 --- a/inst/stan/sm-loglogistic/model.stan +++ b/inst/stan/sm-loglogistic/model.stan @@ -11,15 +11,7 @@ functions { } } -parameters { - // SurvivalLogLogistic - real sm_loglogis_a; - real sm_loglogis_b; -} - - transformed parameters { // SurvivalLogLogistic vector[2] pars_os = [sm_loglogis_a, sm_loglogis_b]'; } - diff --git a/inst/stan/sm-weibull-ph/model.stan b/inst/stan/sm-weibull-ph/model.stan index 629a730c..96bc2e56 100755 --- a/inst/stan/sm-weibull-ph/model.stan +++ b/inst/stan/sm-weibull-ph/model.stan @@ -12,15 +12,7 @@ functions { } -parameters { - // SurvivalWeibullPH - real sm_weibull_ph_lambda; - real sm_weibull_ph_gamma; -} - - transformed parameters { // SurvivalWeibullPH vector[2] pars_os = [sm_weibull_ph_lambda, sm_weibull_ph_gamma]'; } - diff --git a/man/JointModel-class.Rd b/man/JointModel-class.Rd index b2b845b2..d42ecb3a 100644 --- a/man/JointModel-class.Rd +++ b/man/JointModel-class.Rd @@ -28,7 +28,7 @@ Joint Model Object and Constructor Function \item{\code{link}}{(\code{Link})\cr the link.} -\item{\code{parameters}}{(\code{ParameterList})\cr the parameter specification.} +\item{\code{parameters}}{(\code{ParameterList})\cr the parameter declaration and prior specification.} }} \seealso{ diff --git a/man/LongitudinalModel-class.Rd b/man/LongitudinalModel-class.Rd index 35090db8..8b8bab4b 100644 --- a/man/LongitudinalModel-class.Rd +++ b/man/LongitudinalModel-class.Rd @@ -16,9 +16,9 @@ LongitudinalModel( ) } \arguments{ -\item{stan}{(\code{StanModule}) code containing the Stan code specification.} +\item{stan}{(\code{StanModule}) code containing the model-specific Stan code specification.} -\item{parameters}{(\code{ParameterList}) the parameter specification.} +\item{parameters}{(\code{ParameterList}) the parameter declaration and prior specification.} \item{name}{(\code{character}) display name for the model object.} diff --git a/man/ParameterList-class.Rd b/man/ParameterList-class.Rd index 2f74bdcd..8ab68009 100644 --- a/man/ParameterList-class.Rd +++ b/man/ParameterList-class.Rd @@ -14,7 +14,8 @@ ParameterList(...) } \description{ This class extends the general \code{\link{list}} type for containing \code{\link{Parameter}} -specifications. +specifications. When converted to a \code{\link{StanModule}}, the list supplies the +Stan declarations and prior statements for those parameters. } \section{Slots}{ diff --git a/man/Prior-class.Rd b/man/Prior-class.Rd index 6ee77f5b..b04bc554 100644 --- a/man/Prior-class.Rd +++ b/man/Prior-class.Rd @@ -19,7 +19,8 @@ Prior( repr_transformed_parameters = "", repr_generated_quantities = "", limits = c(-Inf, Inf), - .allow_vectors = FALSE + .allow_vectors = FALSE, + .is_const = FALSE ) } \arguments{ @@ -47,6 +48,8 @@ the same names as the \code{paramaters} slot.} \item{limits}{(\code{numeric}) the lower and upper limits for a truncated distribution} \item{.allow_vectors}{(\code{flag}) whether to allow vector parameters.} + +\item{.is_const}{(\code{flag}) whether this prior fixes the parameter at a constant value.} } \description{ Specifies the prior distribution in a Stan Model @@ -77,6 +80,8 @@ Specifies the prior distribution in a Stan Model \item{\code{limits}}{(\code{numeric})\cr See arguments.} \item{\code{.allow_vectors}}{(\code{logical})\cr See arguments.} + +\item{\code{.is_const}}{(\code{logical})\cr See arguments.} }} \seealso{ diff --git a/man/StanModel-class.Rd b/man/StanModel-class.Rd index 6c29171d..91cf9f73 100644 --- a/man/StanModel-class.Rd +++ b/man/StanModel-class.Rd @@ -10,9 +10,9 @@ StanModel(stan, parameters, name = "") } \arguments{ -\item{stan}{(\code{StanModule}) code containing the Stan code specification.} +\item{stan}{(\code{StanModule}) code containing the model-specific Stan code specification.} -\item{parameters}{(\code{ParameterList}) the parameter specification.} +\item{parameters}{(\code{ParameterList}) the parameter declaration and prior specification.} \item{name}{(\code{character}) display name for the model object.} } diff --git a/man/SurvivalModel-class.Rd b/man/SurvivalModel-class.Rd index 745d5c0b..2d217c0d 100644 --- a/man/SurvivalModel-class.Rd +++ b/man/SurvivalModel-class.Rd @@ -15,9 +15,9 @@ SurvivalModel( ) } \arguments{ -\item{stan}{(\code{StanModule}) code containing the Stan code specification.} +\item{stan}{(\code{StanModule}) code containing the model-specific Stan code specification.} -\item{parameters}{(\code{ParameterList}) the parameter specification.} +\item{parameters}{(\code{ParameterList}) the parameter declaration and prior specification.} \item{name}{(\code{character}) display name for the model object.} diff --git a/man/as.StanModule.ParameterDeclaration.Rd b/man/as.StanModule.ParameterDeclaration.Rd new file mode 100644 index 00000000..f603830c --- /dev/null +++ b/man/as.StanModule.ParameterDeclaration.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Parameter.R +\name{as.StanModule.ParameterDeclaration} +\alias{as.StanModule.ParameterDeclaration} +\title{\code{Parameter} Declaration -> \code{StanModule}} +\usage{ +\method{as.StanModule}{ParameterDeclaration}(object) +} +\arguments{ +\item{object}{(\code{Parameter}) a prior Distribution} +} +\value{ +A \code{\link{StanModule}} object. +} +\description{ +Creates only the Stan declaration for a \code{\link{Parameter}}. Sampled parameters are +declared in the \code{parameters} block; constant parameters are declared in the +\verb{transformed parameters} block. +} diff --git a/man/covariates.Rd b/man/covariates.Rd index 6d183a02..5837075a 100644 --- a/man/covariates.Rd +++ b/man/covariates.Rd @@ -7,7 +7,7 @@ covariates(object, ...) } \arguments{ -\item{}{(\if{html}{\out{}}) \if{html}{\out{}}} +\item{object}{(\code{ANY}) the object to extract covariate names from.} } \description{ Extract Covariate Names diff --git a/man/prior_beta.Rd b/man/prior_beta.Rd index e2bf46ce..d53b62a5 100644 --- a/man/prior_beta.Rd +++ b/man/prior_beta.Rd @@ -17,6 +17,7 @@ Beta Prior Distribution \seealso{ Other Prior: \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/prior_cauchy.Rd b/man/prior_cauchy.Rd index 22864916..4204820e 100644 --- a/man/prior_cauchy.Rd +++ b/man/prior_cauchy.Rd @@ -17,6 +17,7 @@ Cauchy Prior Distribution \seealso{ Other Prior: \code{\link[=prior_beta]{prior_beta()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/prior_const.Rd b/man/prior_const.Rd new file mode 100644 index 00000000..7012a30b --- /dev/null +++ b/man/prior_const.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Prior.R +\name{prior_const} +\alias{prior_const} +\title{Constant Prior Distribution} +\usage{ +prior_const(value) +} +\arguments{ +\item{value}{(\code{number}) the fixed parameter value.} +} +\description{ +Constant Prior 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_horseshoe]{prior_horseshoe()}}, +\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_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()}} +} +\concept{Prior} diff --git a/man/prior_gamma.Rd b/man/prior_gamma.Rd index daa80b47..78bee628 100644 --- a/man/prior_gamma.Rd +++ b/man/prior_gamma.Rd @@ -18,6 +18,7 @@ Gamma Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, diff --git a/man/prior_horseshoe.Rd b/man/prior_horseshoe.Rd index d5dfceb5..046c3f21 100644 --- a/man/prior_horseshoe.Rd +++ b/man/prior_horseshoe.Rd @@ -32,6 +32,7 @@ Regularized Horseshoe Prior for a Vector Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_init_only]{prior_init_only()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, diff --git a/man/prior_init_only.Rd b/man/prior_init_only.Rd index f3c0019d..68e3c8c7 100644 --- a/man/prior_init_only.Rd +++ b/man/prior_init_only.Rd @@ -18,6 +18,7 @@ are fixed within the model and cannot be altered by the user. Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, diff --git a/man/prior_invgamma.Rd b/man/prior_invgamma.Rd index 7ad4f91a..89547e7a 100644 --- a/man/prior_invgamma.Rd +++ b/man/prior_invgamma.Rd @@ -18,6 +18,7 @@ Inverse-Gamma Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/prior_logistic.Rd b/man/prior_logistic.Rd index 478873ba..0b05f70e 100644 --- a/man/prior_logistic.Rd +++ b/man/prior_logistic.Rd @@ -18,6 +18,7 @@ Logistic Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/prior_loglogistic.Rd b/man/prior_loglogistic.Rd index 986907dc..48739510 100644 --- a/man/prior_loglogistic.Rd +++ b/man/prior_loglogistic.Rd @@ -18,6 +18,7 @@ Log-Logistic Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/prior_lognormal.Rd b/man/prior_lognormal.Rd index 3354db94..db21e500 100644 --- a/man/prior_lognormal.Rd +++ b/man/prior_lognormal.Rd @@ -18,6 +18,7 @@ Log-Normal Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/prior_normal.Rd b/man/prior_normal.Rd index d9b1c2ba..ec1abaac 100644 --- a/man/prior_normal.Rd +++ b/man/prior_normal.Rd @@ -18,6 +18,7 @@ Normal Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/prior_normal_vector.Rd b/man/prior_normal_vector.Rd index 8c10cd38..e522cade 100644 --- a/man/prior_normal_vector.Rd +++ b/man/prior_normal_vector.Rd @@ -18,6 +18,7 @@ Normal Prior for a Vector Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/prior_std_normal.Rd b/man/prior_std_normal.Rd index 3fe3f752..b41985dc 100644 --- a/man/prior_std_normal.Rd +++ b/man/prior_std_normal.Rd @@ -13,6 +13,7 @@ Standard Normal Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/prior_student_t.Rd b/man/prior_student_t.Rd index b8607b21..f1bbbc09 100644 --- a/man/prior_student_t.Rd +++ b/man/prior_student_t.Rd @@ -20,6 +20,7 @@ Student-t Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/prior_uniform.Rd b/man/prior_uniform.Rd index 6f3ae289..2b74bb1b 100644 --- a/man/prior_uniform.Rd +++ b/man/prior_uniform.Rd @@ -18,6 +18,7 @@ Uniform Prior Distribution Other Prior: \code{\link[=prior_beta]{prior_beta()}}, \code{\link[=prior_cauchy]{prior_cauchy()}}, +\code{\link[=prior_const]{prior_const()}}, \code{\link[=prior_gamma]{prior_gamma()}}, \code{\link[=prior_horseshoe]{prior_horseshoe()}}, \code{\link[=prior_init_only]{prior_init_only()}}, diff --git a/man/render_stan_const_declaration.Rd b/man/render_stan_const_declaration.Rd new file mode 100644 index 00000000..a8e38dbd --- /dev/null +++ b/man/render_stan_const_declaration.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Parameter.R +\name{render_stan_const_declaration} +\alias{render_stan_const_declaration} +\title{Render Stan Constant Parameter Declaration} +\usage{ +render_stan_const_declaration(name, size, limits) +} +\arguments{ +\item{name}{(\code{string}) parameter name.} + +\item{size}{(\code{numeric_OR_character}) parameter size.} + +\item{limits}{(\code{numeric}) lower and upper parameter limits.} +} +\value{ +A length-one character vector. +} +\description{ +Creates a Stan transformed-parameter declaration that fixes a scalar or vector +parameter at a data-supplied constant value. +} +\keyword{internal} diff --git a/man/render_stan_parameter_declaration.Rd b/man/render_stan_parameter_declaration.Rd new file mode 100644 index 00000000..1b38fc6e --- /dev/null +++ b/man/render_stan_parameter_declaration.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Parameter.R +\name{render_stan_parameter_declaration} +\alias{render_stan_parameter_declaration} +\title{Render Stan Parameter Declaration} +\usage{ +render_stan_parameter_declaration(name, size, limits) +} +\arguments{ +\item{name}{(\code{string}) parameter name.} + +\item{size}{(\code{numeric_OR_character}) parameter size.} + +\item{limits}{(\code{numeric}) lower and upper parameter limits.} +} +\value{ +A length-one character vector. +} +\description{ +Creates a Stan declaration for a sampled scalar or vector parameter. +} +\keyword{internal} diff --git a/man/render_stan_parameter_limits.Rd b/man/render_stan_parameter_limits.Rd new file mode 100644 index 00000000..7d7122f0 --- /dev/null +++ b/man/render_stan_parameter_limits.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Parameter.R +\name{render_stan_parameter_limits} +\alias{render_stan_parameter_limits} +\title{Render Stan Parameter Constraints} +\usage{ +render_stan_parameter_limits(limits) +} +\arguments{ +\item{limits}{(\code{numeric}) length-two vector containing lower and upper parameter limits.} +} +\value{ +A length-one character vector. +} +\description{ +Converts a lower/upper limit vector from a \code{\link{Prior}} into a Stan declaration +constraint, e.g. \verb{}. +} +\keyword{internal} diff --git a/man/stanmodel_arguments.Rd b/man/stanmodel_arguments.Rd index 80157cf4..8cd038b3 100644 --- a/man/stanmodel_arguments.Rd +++ b/man/stanmodel_arguments.Rd @@ -4,11 +4,11 @@ \alias{stanmodel_arguments} \title{\code{StanModel} Function Arguments} \arguments{ -\item{stan}{(\code{StanModule}) code containing the Stan code specification.} +\item{stan}{(\code{StanModule}) code containing the model-specific Stan code specification.} -\item{parameters}{(\code{ParameterList}) the parameter specification.} +\item{parameters}{(\code{ParameterList}) the parameter declaration and prior specification.} -\item{parameter}{(\code{ParameterList}) the (single) parameter specification.} +\item{parameter}{(\code{ParameterList}) the (single) parameter declaration and prior specification.} \item{name}{(\code{character}) display name for the model object.} diff --git a/tests/testthat/_snaps/JointModel.md b/tests/testthat/_snaps/JointModel.md index 9a9555fe..cf991fda 100644 --- a/tests/testthat/_snaps/JointModel.md +++ b/tests/testthat/_snaps/JointModel.md @@ -1,3 +1,66 @@ +# JointModel snapshots assembled parameter declarations + + Code + snapshot_joint_model_parameter_declarations(JointModel(longitudinal = LongitudinalRandomSlope(), + survival = SurvivalWeibullPH(), link = linkDSLD())) + Output + parameters + vector[n_studies] lm_rs_intercept; + vector[n_arms] lm_rs_slope_mu; + vector[n_arms] lm_rs_slope_sigma; + real lm_rs_sigma; + vector[n_subjects] lm_rs_ind_rnd_slope; + real sm_weibull_ph_lambda; + real sm_weibull_ph_gamma; + vector[p_os_cov_design] beta_os_cov; + real link_dsld; + + transformed parameter constants + + +--- + + Code + snapshot_joint_model_parameter_declarations(JointModel(longitudinal = LongitudinalGSF( + centred = FALSE), survival = SurvivalWeibullPH(), link = Link(linkTTG(), + linkDSLD(), linkGrowth()))) + Output + parameters + vector[n_studies] lm_gsf_mu_bsld; + vector[n_arms] lm_gsf_mu_ks; + vector[n_arms] lm_gsf_mu_kg; + vector[n_arms] lm_gsf_mu_phi; + vector[n_studies] lm_gsf_omega_bsld; + vector[n_arms] lm_gsf_omega_ks; + vector[n_arms] lm_gsf_omega_kg; + vector[n_arms] lm_gsf_omega_phi; + real lm_gsf_sigma; + vector[n_subjects] lm_gsf_eta_tilde_bsld; + vector[n_subjects] lm_gsf_eta_tilde_ks; + vector[n_subjects] lm_gsf_eta_tilde_kg; + vector[n_subjects] lm_gsf_eta_tilde_phi; + real sm_weibull_ph_lambda; + real sm_weibull_ph_gamma; + vector[p_os_cov_design] beta_os_cov; + real link_ttg; + real link_dsld; + real link_growth; + + transformed parameter constants + + +--- + + Code + snapshot_joint_model_parameter_declarations(JointModel(survival = SurvivalExponential( + lambda = prior_const(0.5)))) + Output + parameters + vector[p_os_cov_design] beta_os_cov; + + transformed parameter constants + real sm_exp_lambda = prior_const_sm_exp_lambda; + # JointModel print method works as expected Code diff --git a/tests/testthat/_snaps/Prior.md b/tests/testthat/_snaps/Prior.md index 908f17bb..c453b0e5 100644 --- a/tests/testthat/_snaps/Prior.md +++ b/tests/testthat/_snaps/Prior.md @@ -8,6 +8,16 @@ cauchy(mu = 0, sigma = 0.8) +--- + + Code + print(prior_const(1)) + Output + + Prior Object: + const(value = 1) + + --- Code diff --git a/tests/testthat/_snaps/SurvivalExponential.md b/tests/testthat/_snaps/SurvivalExponential.md index 844e993b..24805437 100644 --- a/tests/testthat/_snaps/SurvivalExponential.md +++ b/tests/testthat/_snaps/SurvivalExponential.md @@ -22,6 +22,18 @@ beta_os_cov ~ gamma(alpha = 3, beta = 4) +--- + + Code + x <- SurvivalExponential(lambda = prior_const(1)) + print(x) + Output + + Exponential Survival Model with parameters: + sm_exp_lambda = const(value = 1) + beta_os_cov ~ normal(mu = 0, sigma = 2) + + # Different priors for the beta components are possible Code diff --git a/tests/testthat/models/ParameterList.stan b/tests/testthat/models/ParameterList.stan index d017b72c..68316a28 100644 --- a/tests/testthat/models/ParameterList.stan +++ b/tests/testthat/models/ParameterList.stan @@ -6,10 +6,14 @@ data { real prior_sigma_myp; } +parameters { + real inter; + real myp; +} + model { inter ~ gamma(prior_alpha_inter, prior_beta_inter); myp ~ normal(prior_mu_myp, prior_sigma_myp); } - diff --git a/tests/testthat/test-JointModel.R b/tests/testthat/test-JointModel.R index 6c447a66..0e6f9583 100644 --- a/tests/testthat/test-JointModel.R +++ b/tests/testthat/test-JointModel.R @@ -11,7 +11,70 @@ test_that("JointModel smoke tests", { }) +# Keep declaration snapshots focused on semantic Stan declarations, rather than +# churn from blank lines or source comments introduced by StanModule merging. +normalise_stan_declarations <- function(x) { + x <- trimws(x) + x <- x[nzchar(x)] + x[!grepl("^//", x)] +} + +# Snapshot only the assembled parameter declarations and constant-parameter +# declarations. This catches regressions in dynamic parameter insertion without +# snapshotting the full generated Stan program. +snapshot_joint_model_parameter_declarations <- function(object) { + stan <- as.StanModule(object) + cat("parameters\n") + cat(normalise_stan_declarations(stan@parameters), sep = "\n") + + const_declarations <- normalise_stan_declarations( + stan@transformed_parameters + ) + const_declarations <- const_declarations[grepl( + "prior_const_", + const_declarations, + fixed = TRUE + )] + + cat("\ntransformed parameter constants\n") + if (length(const_declarations) == 0) { + cat("\n") + } else { + cat(const_declarations, sep = "\n") + } +} + + +test_that("JointModel snapshots assembled parameter declarations", { + withr::local_options(list(jmpost.double_eps = 0)) + + expect_snapshot( + snapshot_joint_model_parameter_declarations(JointModel( + longitudinal = LongitudinalRandomSlope(), + survival = SurvivalWeibullPH(), + link = linkDSLD() + )) + ) + + expect_snapshot( + snapshot_joint_model_parameter_declarations(JointModel( + longitudinal = LongitudinalGSF(centred = FALSE), + survival = SurvivalWeibullPH(), + link = Link(linkTTG(), linkDSLD(), linkGrowth()) + )) + ) + + expect_snapshot( + snapshot_joint_model_parameter_declarations(JointModel( + survival = SurvivalExponential(lambda = prior_const(0.5)) + )) + ) +}) + + test_that("JointModel print method works as expected", { + withr::local_options(list(jmpost.double_eps = 0)) + expect_snapshot({ x <- JointModel( longitudinal = LongitudinalRandomSlope(), diff --git a/tests/testthat/test-Link.R b/tests/testthat/test-Link.R index 2f8f31f1..3d8dcbb2 100644 --- a/tests/testthat/test-Link.R +++ b/tests/testthat/test-Link.R @@ -57,8 +57,11 @@ test_that("Link works as expected", { ## Fully resolved stan code pass's the syntax checker expect_stan_syntax( merge( - load_with_base_stan("lm-gsf/functions.stan"), - as.StanModule(x) + merge( + load_with_base_stan("lm-gsf/functions.stan"), + as.StanModule(x) + ), + as.StanModule(getParameters(x)) ) ) diff --git a/tests/testthat/test-Parameter.R b/tests/testthat/test-Parameter.R index 58d43112..1e382807 100644 --- a/tests/testthat/test-Parameter.R +++ b/tests/testthat/test-Parameter.R @@ -12,6 +12,68 @@ test_that("Parameters smoke tests", { expect_equal(names(p), "intercept") }) +test_that("Parameter declaration helpers render Stan code", { + eps <- getOption("jmpost.double_eps") + + expect_equal( + render_stan_parameter_limits(c(-Inf, Inf)), + "" + ) + expect_equal( + render_stan_parameter_limits(c(eps, 1)), + paste0("") + ) + + expect_equal( + render_stan_parameter_declaration("theta", 1, c(eps, Inf)), + paste0("real theta;") + ) + expect_equal( + render_stan_parameter_declaration("theta", "n_subjects", c(-Inf, Inf)), + "vector[n_subjects] theta;" + ) + + expect_equal( + render_stan_const_declaration("theta", 1, c(eps, Inf)), + paste0("real theta = prior_const_theta;") + ) + expect_equal( + render_stan_const_declaration("theta", "n_subjects", c(eps, Inf)), + paste0( + "vector[n_subjects] theta = ", + "rep_vector(prior_const_theta, n_subjects);" + ) + ) +}) + +test_that("Parameter declaration blocks depend on prior_const", { + eps <- getOption("jmpost.double_eps") + + sampled <- Parameter( + name = "theta", + prior = set_limits(prior_normal(0, 1), lower = eps), + size = 1 + ) + sampled_module <- as.StanModule.ParameterDeclaration(sampled) + expect_equal( + sampled_module@parameters, + paste0(" real theta;") + ) + expect_equal(sampled_module@transformed_parameters, "") + + fixed <- Parameter( + name = "theta", + prior = set_limits(prior_const(1), lower = eps), + size = 1 + ) + fixed_module <- as.StanModule.ParameterDeclaration(fixed) + expect_equal(fixed_module@parameters, "") + expect_equal( + fixed_module@transformed_parameters, + paste0(" real theta = prior_const_theta;") + ) +}) + test_that("show() works for Paramneter objects", { x <- Parameter(prior_normal(1, 3), "bob", "size1") expect_snapshot(print(x)) diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 94908c48..2ff65062 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -138,11 +138,46 @@ test_that("Invalid prior parameters are rejected", { # Ensure that validation doesn't wrongly reject priors with no user specified parameters expect_s4_class(prior_init_only(prior_normal(3, 1)), "Prior") expect_s4_class(prior_std_normal(), "Prior") + expect_s4_class(prior_const(1), "Prior") +}) + +test_that("prior_const() works as expected", { + x <- prior_const(1) + + expect_equal(initialValues(x), 1) + expect_equal( + as_stan_list(x, name = "sm_exp_lambda"), + list(prior_const_sm_exp_lambda = 1) + ) + expect_equal( + as.StanModule(x, name = "sm_exp_lambda")@data, + " real prior_const_sm_exp_lambda;" + ) + expect_equal(as.StanModule(x, name = "sm_exp_lambda")@model, "") + expect_equal( + as.StanModule(x, name = "bob")@data, + " real prior_const_bob;" + ) + + pars <- ParameterList( + Parameter(name = "fixed", prior = x, size = 1), + Parameter(name = "sampled", prior = prior_normal(0, 1), size = 1) + ) + expect_equal(names(initialValues(pars, n_chains = 1)[[1]]), "sampled") + expect_equal( + as_stan_list(pars), + list( + prior_const_fixed = 1, + prior_mu_sampled = 0, + prior_sigma_sampled = 1 + ) + ) }) test_that("show() works for Prior objects", { expect_snapshot(print(prior_cauchy(0, 0.8))) + expect_snapshot(print(prior_const(1))) expect_snapshot(print(prior_normal(0, 0.8))) expect_snapshot(print(prior_std_normal())) expect_snapshot(print(prior_beta(5, 1))) diff --git a/tests/testthat/test-SurvivalExponential.R b/tests/testthat/test-SurvivalExponential.R index 5592cb45..805cc897 100644 --- a/tests/testthat/test-SurvivalExponential.R +++ b/tests/testthat/test-SurvivalExponential.R @@ -16,6 +16,50 @@ test_that("Can load and compile SurvivalExponential() model", { expect_stan_syntax(x) }) +test_that("SurvivalExponential() can fix lambda with prior_const()", { + default_model <- JointModel(survival = SurvivalExponential()) + default_stan_code <- paste( + as.character(as.StanModule(default_model)), + collapse = "\n" + ) + expect_match( + default_stan_code, + "real]+> sm_exp_lambda;" + ) + + jm <- JointModel( + survival = SurvivalExponential(lambda = prior_const(1)) + ) + + x <- as.StanModule(jm) + expect_stan_syntax(x) + + stan_code <- paste(as.character(x), collapse = "\n") + expect_match(stan_code, "real prior_const_sm_exp_lambda;", fixed = TRUE) + expect_match( + stan_code, + "real]+> sm_exp_lambda = prior_const_sm_exp_lambda;" + ) + expect_false(grepl( + "real]+> sm_exp_lambda;", + stan_code, + perl = TRUE + )) + + expect_equal( + as_stan_list(jm@parameters)$prior_const_sm_exp_lambda, + 1 + ) + expect_false( + "sm_exp_lambda" %in% names(initialValues(jm, n_chains = 1)[[1]]) + ) + + beta_const <- JointModel( + survival = SurvivalExponential(beta = prior_const(0)) + ) + expect_stan_syntax(as.StanModule(beta_const)) +}) + test_that("SurvivalExponential can recover true parameter (including covariates)", { skip_if_not(is_full_test()) @@ -93,6 +137,11 @@ test_that("Print method for SurvivalExponential works as expected", { x <- SurvivalExponential(beta = prior_gamma(3, 4)) print(x) }) + + expect_snapshot({ + x <- SurvivalExponential(lambda = prior_const(1)) + print(x) + }) }) test_that("Different priors for the beta components are possible", { diff --git a/tests/testthat/test-brms_compare.R b/tests/testthat/test-brms_compare.R index 68722776..13318f7c 100644 --- a/tests/testthat/test-brms_compare.R +++ b/tests/testthat/test-brms_compare.R @@ -293,12 +293,13 @@ test_that("jmpost and brms get similar horseshoe estimates for survival models", "horseshoe(df = 1, df_global = 1, df_slab = 4, scale_global = 0.3, scale_slab = 2)", class = "b" ), - warmup = 50, - iter = 100, + warmup = 5000, + iter = 10000, chains = 4, cores = 4, seed = 9826, - backend = "cmdstanr" + backend = "cmdstanr", + save_pars = brms::save_pars(manual = "hs_local") ) dat_surv2 <- dat_surv |> @@ -328,7 +329,7 @@ test_that("jmpost and brms get similar horseshoe estimates for survival models", survival = SurvivalExponential( # Let's use an approximate point mass # at lambda = 1 to get the same model as brms. - lambda = prior_normal(1, 0.00001), + lambda = prior_const(1), beta = prior_horseshoe( df = 1, df_global = 1, @@ -342,8 +343,8 @@ test_that("jmpost and brms get similar horseshoe estimates for survival models", mp <- sampleStanModel( jm, data = jdat, - iter_warmup = 50, - iter_sampling = 50, + iter_warmup = 5000, + iter_sampling = 5000, chains = 4, refresh = 200, parallel_chains = 4, @@ -368,7 +369,13 @@ test_that("jmpost and brms get similar horseshoe estimates for survival models", posterior::subset_draws( variable = c("b_ones", "b_cov1", "b_cov2", "b_cov3") ) |> - posterior::as_draws_df() + posterior::as_draws_df() |> + dplyr::mutate( + b_ones = -b_ones, + b_cov1 = -b_cov1, + b_cov2 = -b_cov2, + b_cov3 = -b_cov3 + ) j_beta <- j_draws |> posterior::subset_draws(variable = "beta_os_cov") |> @@ -377,23 +384,49 @@ test_that("jmpost and brms get similar horseshoe estimates for survival models", expect_equal( summary_quantiles(j_beta), summary_quantiles(b_beta), - tolerance = 0.1, + tolerance = 0.01, ignore_attr = TRUE ) - b_c2 <- 2^2 * b_draws$hs_slab + # Also compare shrinkage factors. + b_shrink_pars <- b_draws |> + posterior::subset_draws( + variable = c("hs_slab", "hs_global", "hs_local") + ) |> + posterior::as_draws_df() + + b_c2 <- 2^2 * b_shrink_pars$hs_slab b_shrinkage <- dplyr::tibble( - cov1 = 1 / (1 + b_draws$hs_global^2 * b_draws$`hs_local[1]`^2 / b_c2), - cov2 = 1 / (1 + b_draws$hs_global^2 * b_draws$`hs_local[2]`^2 / b_c2), - cov3 = 1 / (1 + b_draws$hs_global^2 * b_draws$`hs_local[3]`^2 / b_c2) + cov1 = 1 / + (1 + + b_shrink_pars$hs_global^2 * + b_shrink_pars$`hs_local[1]`^2 / + b_c2), + cov2 = 1 / + (1 + + b_shrink_pars$hs_global^2 * + b_shrink_pars$`hs_local[2]`^2 / + b_c2), + cov3 = 1 / + (1 + + b_shrink_pars$hs_global^2 * + b_shrink_pars$`hs_local[3]`^2 / + b_c2), + cov4 = 1 / + (1 + + b_shrink_pars$hs_global^2 * + b_shrink_pars$`hs_local[4]`^2 / + b_c2) ) - j_shrinkage <- shrinkage(mp) + j_shrinkage <- shrinkage(mp) |> + posterior::as_draws_df() |> + dplyr::select(ones, cov1, cov2, cov3) expect_equal( summary_quantiles(j_shrinkage), summary_quantiles(b_shrinkage), - tolerance = 0.1, + tolerance = 0.01, ignore_attr = TRUE ) }) diff --git a/vignettes/covariate-selection-horseshoe.Rmd b/vignettes/covariate-selection-horseshoe.Rmd new file mode 100644 index 00000000..7a4eb3c5 --- /dev/null +++ b/vignettes/covariate-selection-horseshoe.Rmd @@ -0,0 +1,318 @@ +--- +title: "Covariate Selection with the Horseshoe Prior" +package: jmpost +bibliography: '`r system.file("REFERENCES.bib", package = "jmpost")`' +csl: '`r system.file("jss.csl", package = "jmpost")`' +output: + rmarkdown::html_vignette: + toc: true +link-citations: true +vignette: > + %\VignetteIndexEntry{Covariate Selection with the Horseshoe Prior} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +editor_options: + chunk_output_type: console +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +library(jmpost) +``` + +## Overview + +Covariate selection in joint models is often difficult because each candidate +covariate can interact with the survival model, the longitudinal model, and the +association between them. The regularized horseshoe prior provides a Bayesian +shrinkage approach: include the candidate covariates in one model, shrink most +survival coefficients strongly towards zero, and allow a smaller number of +coefficients to escape that shrinkage when supported by the data. + +In `jmpost`, this workflow is available for the survival covariate coefficients +through `prior_horseshoe()`. This vignette starts from the prior specification +and then shows a complete joint model example in which the survival covariates +use the horseshoe prior. The fitted model can be inspected through the posterior +coefficient draws and through the shrinkage factors extracted with `shrinkage()`. + +## Theory + +The horseshoe prior is a global-local shrinkage prior [@CarvalhoEtAl2010]. A +global parameter shrinks all coefficients towards zero, while coefficient-level +local parameters allow selected coefficients to remain large. The regularized +horseshoe [@PiironenVehtari2017] adds a slab component that regularizes +coefficients that escape the main shrinkage. This is the parameterization used +by `jmpost`, following the same notation as `brms` [@Buerkner2017]. +Details are provided in the statistical specifications vignette [here](statistical-specification.html). + +The hyperparameters map to `prior_horseshoe()` +arguments as follows: + +- `df`: local degrees of freedom, $\nu_{\lambda}$ +- `df_global`: global degrees of freedom, $\nu_{\tau}$ +- `df_slab`: slab degrees of freedom, $\nu_c$ +- `scale_global`: global scale, $s_{\tau}$ +- `scale_slab`: slab scale, $s_c$ + +The local shrinkage parameters $\lambda_j$ let individual coefficients escape +the global shrinkage $\tau$. The slab scale $s_c$ and slab degrees of freedom +$\nu_c$ control how much very large coefficients are regularized. Smaller +`scale_global` values encode stronger prior belief that only a small number of +candidate covariates are relevant. + +A useful posterior diagnostic is the shrinkage factor + +$$ +\kappa_j = +\frac{1}{1 + \tau^2\lambda_j^2 / c^2}. +$$ + +Values of $\kappa_j$ close to 1 indicate strong shrinkage towards zero. Values +close to 0 indicate little shrinkage. The shrinkage factor should be read +together with the coefficient posterior: it is evidence about how strongly the +prior-likelihood combination pulled a coefficient towards zero, not a +standalone posterior inclusion probability. + +## Simulate Example Data + +We first simulate data from a joint model with a random-slope longitudinal +submodel and a Weibull proportional hazards survival submodel. The survival +model contains one categorical and one continuous candidate covariate. The +categorical covariate contributes two columns to the design matrix, so the +survival model has three covariate coefficients in total. + +```{r simulate-data} +set.seed(129) +sim_data <- SimJointData( + design = list( + SimGroup(50, "Arm-A", "Study-X"), + SimGroup(50, "Arm-B", "Study-X") + ), + longitudinal = SimLongitudinalRandomSlope( + times = c(1, 20, 50, 100, 150, 200, 250, 300), + intercept = 30, + slope_mu = c(1, 2), + slope_sigma = 0.2, + sigma = 20, + link_dsld = 0.1 + ), + survival = SimSurvivalWeibullPH( + lambda = 1 / 300, + gamma = 0.97, + time_max = 2000, + time_step = 1, + lambda_cen = 1 / 9000, + beta_cat = c( + "A" = 0, + "B" = -0.1, + "C" = 0.5 + ), + beta_cont = 0.3 + ) +) +``` + +```{r prepare-data} +os_data <- sim_data@survival +long_data <- sim_data@longitudinal + +joint_data <- DataJoint( + subject = DataSubject( + data = os_data, + subject = "subject", + arm = "arm", + study = "study" + ), + survival = DataSurvival( + data = os_data, + formula = Surv(time, event) ~ cov_cat + cov_cont + ), + longitudinal = DataLongitudinal( + data = long_data, + formula = sld ~ time, + threshold = 5 + ) +) +``` + +Before setting priors or interpreting the output, check the survival design +matrix. With `covariates()` we can access the column names, which map to the +corresponding survival coefficients that receive the horseshoe prior. + +```{r design-matrix} +head(model.matrix(joint_data@survival)) +survival_covariates <- covariates(joint_data@survival) +survival_covariates +``` + +## Fit a Joint Model + +Now we define a joint model. The longitudinal model and the association link use +standard priors. The survival model uses `prior_horseshoe()` for the vector of +survival covariate coefficients. + +```{r model-specification} +joint_model <- JointModel( + longitudinal = LongitudinalRandomSlope(), + survival = SurvivalWeibullPH( + beta = prior_horseshoe( + df = 1, + df_global = 1, + df_slab = 4, + scale_global = 0.3, + scale_slab = 2 + ) + ), + link = linkDSLD() +) +``` + +The choice above uses half-Cauchy priors for the local and global shrinkage +parameters, because `df = 1` and `df_global = 1`. The relatively small +`scale_global = 0.3` favours sparse survival effects in this small example, +while `scale_slab = 2` still allows meaningfully large log-hazard coefficients +when the data support them. + +The following code fits the model. In a serious analysis, increase the number +of warmup and sampling iterations and check convergence carefully. + +```{r fit-model} +fit <- sampleStanModel( + joint_model, + data = joint_data, + iter_warmup = 500, + iter_sampling = 500, + chains = 4, + parallel_chains = 4, + seed = 325, + refresh = 200 +) +``` + +The iteration counts above are chosen to keep this vignette reasonably fast to +render. Refit with longer chains before making inferential claims. + +## Inspect Coefficients + +After fitting, inspect the coefficient posterior and standard MCMC diagnostics. +The survival covariate coefficients are stored in `beta_os_cov`. + +```{r inspect-coefficients} +stan_fit <- cmdstanr::as.CmdStanMCMC(fit) + +stan_fit$summary( + variables = c( + "beta_os_cov", + "prior_global_beta_os_cov", + "prior_slab_beta_os_cov" + ) +) +``` + +The coefficient names come from the design matrix: + +```{r name-coefficients} +beta_summary <- stan_fit$summary("beta_os_cov") +beta_summary$covariate <- survival_covariates +beta_summary[, c("covariate", "median", "q5", "q95", "rhat", "ess_bulk")] +``` + +Coefficients whose posterior remains close to zero and whose shrinkage factors +are close to 1 are natural candidates to treat as weakly supported covariates. +Coefficients whose posterior is away from zero and whose shrinkage factors are +closer to 0 have escaped shrinkage and are more strongly supported by the model. + +## Extract and Plot Shrinkage Factors + +Use `shrinkage()` to extract the posterior draws of $\kappa_j$. The returned +draws are named with the survival covariate names, so they can be plotted or +summarised directly. + +```{r shrinkage-factors} +shrinkage_draws <- shrinkage(fit) +posterior::variables(shrinkage_draws) + +summary(shrinkage_draws) +``` + +One compact visual summary is a density plot of the shrinkage factors. + +```{r shrinkage-density, fig.width=7, fig.height=4} +library(bayesplot) + +mcmc_dens_overlay(shrinkage_draws) + + ggplot2::labs( + x = "Shrinkage factor", + y = "Posterior density" + ) +``` + +A second useful display is the median shrinkage factor with an interval for each +covariate. + +```{r shrinkage-intervals, fig.width=7, fig.height=4} +library(dplyr) +library(ggplot2) + +shrinkage_summary <- posterior::summarise_draws( + shrinkage_draws, + median, + ~ quantile(.x, 0.1), + ~ quantile(.x, 0.9) +) |> + dplyr::rename(q10 = `10%`, q90 = `90%`) + +ggplot(shrinkage_summary, aes(x = median, y = variable)) + + geom_errorbar(aes(xmin = q10, xmax = q90), width = 0) + + geom_point() + + scale_x_continuous(limits = c(0, 1)) + + labs( + x = "Shrinkage factor", + y = NULL + ) +``` + +## Interpret the Output + +In this simulated example, the design matrix contains `cov_catB`, `cov_catC`, +and `cov_cont`. The data were generated with a small effect for `cov_catB`, a +larger positive effect for `cov_catC`, and a moderate positive effect for +`cov_cont`. + +The horseshoe analysis should therefore be interpreted along these lines: + +- A shrinkage factor near 1 for `cov_catB`, combined with a coefficient posterior + close to zero, would indicate that the model sees little evidence for that + weak covariate. +- Smaller shrinkage factors for `cov_catC` and `cov_cont`, combined with + coefficient posteriors away from zero, would indicate that these covariates + escaped the global shrinkage. +- Intermediate shrinkage factors should be treated as uncertainty, not as an + automatic include/exclude decision. + +The practical decision is still scientific. The regularized horseshoe helps rank +and regularize many candidate covariates in a single joint model fit, but final +selection should also consider prior clinical plausibility, multiplicity of +candidate transformations, model diagnostics, and posterior predictive +performance. + +## Practical Guidance + +The global scale is the main sparsity control. Smaller `scale_global` values +expect fewer relevant covariates; larger values allow more coefficients to +remain away from zero. In applied work, this choice should reflect the number of +candidate covariates and the expected number of non-negligible effects. + +The slab scale should be large enough for plausible survival effects but not so +large that implausible log-hazard ratios are effectively unregularized. For +standardized continuous covariates, a `scale_slab` around 2 is often already +wide on the log-hazard scale. For unstandardized covariates, first consider +whether the covariate scale itself should be transformed or standardized. + +Finally, the horseshoe prior is a shrinkage prior, not a replacement for model +checking. Always inspect convergence diagnostics, posterior predictive fit, and +sensitivity to reasonable prior choices before using the selected covariates for +scientific conclusions. diff --git a/vignettes/custom-model.Rmd b/vignettes/custom-model.Rmd index 5046d713..341b4792 100644 --- a/vignettes/custom-model.Rmd +++ b/vignettes/custom-model.Rmd @@ -201,28 +201,65 @@ longmodel <- WangModel( stan = StanModule("custom-model.stan"), scaled_variance = FALSE, parameters = ParameterList( - Parameter(name = "mu_baseline", prior = prior_lognormal(log(60), 1), size = 1), - Parameter(name = "mu_shrinkage", prior = prior_lognormal(log(2), 1), size = 1), - Parameter(name = "mu_growth", prior = prior_lognormal(log(10), 1), size = 1), - Parameter(name = "sigma_baseline", prior = prior_lognormal(0.3, 1), size = 1), - Parameter(name = "sigma_shrinkage", prior = prior_lognormal(0.3, 1), size = 1), - Parameter(name = "sigma_growth", prior = prior_lognormal(0.3, 1), size = 1), - Parameter(name = "sigma", prior = prior_lognormal(1.5, 1), size = 1), + Parameter( + name = "mu_baseline", + prior = set_limits(prior_lognormal(log(60), 1), lower = getOption("jmpost.double_eps")), + size = 1 + ), + Parameter( + name = "mu_shrinkage", + prior = set_limits(prior_lognormal(log(2), 1), lower = getOption("jmpost.double_eps")), + size = 1 + ), + Parameter( + name = "mu_growth", + prior = set_limits(prior_lognormal(log(10), 1), lower = getOption("jmpost.double_eps")), + size = 1 + ), + Parameter( + name = "sigma_baseline", + prior = set_limits(prior_lognormal(0.3, 1), lower = getOption("jmpost.double_eps")), + size = 1 + ), + Parameter( + name = "sigma_shrinkage", + prior = set_limits(prior_lognormal(0.3, 1), lower = getOption("jmpost.double_eps")), + size = 1 + ), + Parameter( + name = "sigma_growth", + prior = set_limits(prior_lognormal(0.3, 1), lower = getOption("jmpost.double_eps")), + size = 1 + ), + Parameter( + name = "sigma", + prior = set_limits(prior_lognormal(1.5, 1), lower = getOption("jmpost.double_eps")), + size = 1 + ), # The following is only required if we want jmpost to generate # initial values automatically for us Parameter( name = "baseline_idv", - prior = prior_init_only(prior_lognormal(log(60), 1)), + prior = set_limits( + prior_init_only(prior_lognormal(log(60), 1)), + lower = getOption("jmpost.double_eps") + ), size = "n_subjects" ), Parameter( name = "shrinkage_idv", - prior = prior_init_only(prior_lognormal(log(2), 1)), + prior = set_limits( + prior_init_only(prior_lognormal(log(2), 1)), + lower = getOption("jmpost.double_eps") + ), size = "n_subjects" ), Parameter( name = "growth_idv", - prior = prior_init_only(prior_lognormal(log(10), 1)), + prior = set_limits( + prior_init_only(prior_lognormal(log(10), 1)), + lower = getOption("jmpost.double_eps") + ), size = "n_subjects" ) ) @@ -233,9 +270,12 @@ longmodel <- WangModel( Please note that: - The `scaled_variance` flag needs to be specified in the `LongitudinalModel()` construction, in order to indicate whether a multiplicative error model is used. Here we have an additive error model, therefore `scaled_variance` is set to `FALSE`. -- The `parameters` argument is used to specify the priors for the model and the `name` argument for the `Parameter`'s objects must match the name of the parameter used within the corresponding Stan code. +- The `parameters` argument declares the model parameters and specifies their priors. The `name` argument for each `Parameter()` must match the name used within the corresponding Stan code. This includes hierarchical subject-level parameters such as `baseline_idv`, where `prior_init_only()` is used so that `jmpost` declares the parameter and can generate initial values, but does not add a prior statement. +- Since this model requires positive parameters, the priors are wrapped in `set_limits()` so that +the generated Stan declarations include lower bounds. -The `StanModule` object contains all of the stan code used to implement the model. For this +The `StanModule` object contains the model-specific Stan code. Parameters listed in +`ParameterList()` should not be declared again in this Stan code. For this particular model the Stan code specified in the `custom-model.stan` file is as follows: ```{R, results='asis', echo=FALSE} x <- readLines("./custom-model.stan") diff --git a/vignettes/custom-model.stan b/vignettes/custom-model.stan index 7f7e66db..6a2ea07c 100644 --- a/vignettes/custom-model.stan +++ b/vignettes/custom-model.stan @@ -9,24 +9,6 @@ functions { } } -parameters{ - // Declare individual subject parameters - vector[n_subjects] baseline_idv; - vector[n_subjects] shrinkage_idv; - vector[n_subjects] growth_idv; - - // Declare population level parameters - real mu_baseline; - real mu_shrinkage; - real mu_growth; - real sigma_baseline; - real sigma_shrinkage; - real sigma_growth; - - // Declare standard deviation for the overall model error - real sigma; -} - transformed parameters{ // Calculated the fitted Tumour values @@ -53,4 +35,3 @@ model { shrinkage_idv ~ lognormal(log(mu_shrinkage), sigma_shrinkage); growth_idv ~ lognormal(log(mu_growth), sigma_growth); } - diff --git a/vignettes/extending-jmpost.Rmd b/vignettes/extending-jmpost.Rmd index 9003ce2c..5b5ebf4d 100644 --- a/vignettes/extending-jmpost.Rmd +++ b/vignettes/extending-jmpost.Rmd @@ -40,20 +40,18 @@ distribution only. Survival distributions are implemented as S4 classes that inherit from `SurvivalModel`. The two main components of the `SurvivalModel` object are the `stan` and `parameters` slots. -The `parameters` slot is a `ParametersList()` object which is used to specify the prior distributions -of each parameter used by the survival distribution as well as for the design matrix coefficients +The `parameters` slot is a `ParameterList()` object which is used to declare each parameter +used by the survival distribution and specify its prior distribution, as well as to declare +the design matrix coefficients (please see the "Prior Specification" section below). The `stan` slot is a `StanModule()` object which needs to define the following: -1) In the `parameters` block each of the distributions parameters must be formally declared. - - -2) In the `transformed parameters` block there must define a vector called `pars_os` which contains one +1) In the `transformed parameters` block you must define a vector called `pars_os` which contains one element for each parameter used by the survival distribution. -3) In the `functions` block there must define a function called `log_h0` which returns the +2) In the `functions` block you must define a function called `log_h0` which returns the log-hazard contribution. This function must have the following signature: ```stan @@ -64,7 +62,7 @@ Where: - The return matrix must be of the same dimensionality as the `time` argument matrix where each element is the log-hazard contribution for the corresponding timepoint in the `time` matrix. - The `time` matrix contains one row per unique subject and one column per timepoint for that subject -- The `pars_os` vector is as defined from (2) above and contains the parameters for that +- The `pars_os` vector is as defined from (1) above and contains the parameters for that particular distribution. For example the following is roughly equivalent to how the Weibull distribution is implemented: @@ -74,6 +72,9 @@ SurvivalWeibullPH <- function( gamma = prior_gamma(2, 0.5), beta = prior_normal(0, 2) ) { + lambda <- set_limits(lambda, lower = getOption("jmpost.double_eps")) + gamma <- set_limits(gamma, lower = getOption("jmpost.double_eps")) + stan <- StanModule(" functions { matrix log_h0(matrix time, vector pars_os) { @@ -83,11 +84,6 @@ SurvivalWeibullPH <- function( } } - parameters { - real sm_weibull_ph_lambda; - real sm_weibull_ph_gamma; - } - transformed parameters { vector[2] pars_os = [sm_weibull_ph_lambda, sm_weibull_ph_gamma]'; } @@ -109,8 +105,8 @@ SurvivalWeibullPH <- function( Similar to the survival model the longitudinal models are implemented as S4 classes that inherit from the `LongitudinalModel` class. The two main components of the `LongitudinalModel` object are -the `stan` and `parameters` slots which specify the underlying Stan code and prior distributions -for the models parameters respectively (please see the "Prior Specification" section below). +the `stan` and `parameters` slots which specify the underlying Stan code and the model parameter +declarations and prior distributions respectively (please see the "Prior Specification" section below). Unlike the survival distributions, the longitudinal models are a lot more flexible and have less constraints on how they are implemented. That is there aren't any specific variables or functions @@ -226,12 +222,14 @@ custom generated quantity code. ## Prior Specification When writing your own custom longitudinal or survival model it is important to understand -how the prior definitions are specified. By default `jmpost` will insert the Stan statements -for any prior distributions based on the `parameters` slot of the model object. -The importance of this is that it means you should not define the prior distributions in the -Stan code itself. Note that this does not apply to hierarchical parameters who must have their -distributions specified in the Stan code. For example in the `LongitudinalRandomSlope` model -their is a different random slope for each treatment arm which is specified in the Stan code as: +how model parameters and prior definitions are specified. By default `jmpost` will insert +the Stan declarations and prior statements for entries in the `parameters` slot of the model +object. This means that parameters listed in `ParameterList()` should not be declared in the +custom Stan code itself, and their prior distributions should not be defined there either. +Note that hierarchical relationships between parameters are still model-specific and must be +specified in the Stan code. For example in the `LongitudinalRandomSlope` model there is a +different random slope for each treatment arm, and this relationship is specified in the Stan +code as: ```stan model { From caa20e6f832cd21f430fa6c46f9161a4a63bdd68 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 12:49:31 +0800 Subject: [PATCH 16/30] some test fixes --- tests/testthat/test-GridPrediction.R | 6 +++--- tests/testthat/test-LongitudinalClaretBruno.R | 2 ++ tests/testthat/test-LongitudinalGSF.R | 4 ++++ tests/testthat/test-LongitudinalRandomSlope.R | 2 ++ tests/testthat/test-LongitudinalSteinFojo.R | 1 + tests/testthat/test-SurvivalExponential.R | 2 ++ tests/testthat/test-SurvivalGamma.R | 2 ++ tests/testthat/test-SurvivalLoglogistic.R | 2 ++ tests/testthat/test-SurvivalQuantities.R | 2 ++ tests/testthat/test-SurvivalWeibullPH.R | 2 ++ tests/testthat/test-brms_compare.R | 1 + 11 files changed, 23 insertions(+), 3 deletions(-) diff --git a/tests/testthat/test-GridPrediction.R b/tests/testthat/test-GridPrediction.R index f955d7b1..47637abe 100644 --- a/tests/testthat/test-GridPrediction.R +++ b/tests/testthat/test-GridPrediction.R @@ -522,9 +522,9 @@ test_that("GridPrediction() works for survival only models", { ) |> dplyr::arrange(group, time) - expect_gt(cor(actual$median, expected$median), 0.999999) - expect_gt(cor(actual$lower, expected$lower), 0.999999) - expect_gt(cor(actual$upper, expected$upper), 0.999999) + expect_gt(cor(actual$median, expected$median), 0.99999) + expect_gt(cor(actual$lower, expected$lower), 0.99999) + expect_gt(cor(actual$upper, expected$upper), 0.99999) expect_equal(actual$time, expected$time) expect_equal(actual$group, expected$group) }) diff --git a/tests/testthat/test-LongitudinalClaretBruno.R b/tests/testthat/test-LongitudinalClaretBruno.R index 67f6c2d7..23e3518c 100644 --- a/tests/testthat/test-LongitudinalClaretBruno.R +++ b/tests/testthat/test-LongitudinalClaretBruno.R @@ -5,6 +5,8 @@ test_that("LongitudinalClaretBruno works as expected with default arguments", { test_that("Print method for LongitudinalClaretBruno works as expected", { + withr::local_options(list(jmpost.double_eps = 0)) + expect_snapshot({ x <- LongitudinalClaretBruno() print(x) diff --git a/tests/testthat/test-LongitudinalGSF.R b/tests/testthat/test-LongitudinalGSF.R index c7850d96..783df469 100644 --- a/tests/testthat/test-LongitudinalGSF.R +++ b/tests/testthat/test-LongitudinalGSF.R @@ -7,6 +7,8 @@ test_that("LongitudinalGSF works as expected with default arguments", { test_that("Print method for LongitudinalGSF works as expected", { + withr::local_options(list(jmpost.double_eps = 0)) + expect_snapshot({ x <- LongitudinalGSF() print(x) @@ -310,6 +312,8 @@ test_that("Quantity models pass the parser", { test_that("Can generate valid initial values", { + withr::local_options(list(jmpost.double_eps = 0)) + pars <- c( "lm_gsf_omega_bsld", "lm_gsf_omega_ks", diff --git a/tests/testthat/test-LongitudinalRandomSlope.R b/tests/testthat/test-LongitudinalRandomSlope.R index 92f769c7..bfae33d3 100644 --- a/tests/testthat/test-LongitudinalRandomSlope.R +++ b/tests/testthat/test-LongitudinalRandomSlope.R @@ -1,4 +1,6 @@ test_that("Print method for LongitudinalRandomSlope works as expected", { + withr::local_options(list(jmpost.double_eps = 0)) + expect_snapshot({ x <- LongitudinalRandomSlope() print(x) diff --git a/tests/testthat/test-LongitudinalSteinFojo.R b/tests/testthat/test-LongitudinalSteinFojo.R index d9302080..2d513942 100644 --- a/tests/testthat/test-LongitudinalSteinFojo.R +++ b/tests/testthat/test-LongitudinalSteinFojo.R @@ -5,6 +5,7 @@ test_that("LongitudinalSteinFojo works as expected with default arguments", { test_that("Print method for LongitudinalSteinFojo works as expected", { + withr::local_options(list(jmpost.double_eps = 0)) expect_snapshot({ x <- LongitudinalSteinFojo() print(x) diff --git a/tests/testthat/test-SurvivalExponential.R b/tests/testthat/test-SurvivalExponential.R index 805cc897..38fe3ebb 100644 --- a/tests/testthat/test-SurvivalExponential.R +++ b/tests/testthat/test-SurvivalExponential.R @@ -128,6 +128,8 @@ test_that("SurvivalExponential can recover true parameter (including covariates) test_that("Print method for SurvivalExponential works as expected", { + withr::local_options(list(jmpost.double_eps = 0)) + expect_snapshot({ x <- SurvivalExponential() print(x) diff --git a/tests/testthat/test-SurvivalGamma.R b/tests/testthat/test-SurvivalGamma.R index ee86a157..febec990 100644 --- a/tests/testthat/test-SurvivalGamma.R +++ b/tests/testthat/test-SurvivalGamma.R @@ -1,4 +1,6 @@ test_that("Print method for SurvivalGamma works as expected", { + withr::local_options(list(jmpost.double_eps = 0)) + expect_snapshot({ x <- SurvivalGamma() print(x) diff --git a/tests/testthat/test-SurvivalLoglogistic.R b/tests/testthat/test-SurvivalLoglogistic.R index 1e9c26fb..5e10bb1f 100644 --- a/tests/testthat/test-SurvivalLoglogistic.R +++ b/tests/testthat/test-SurvivalLoglogistic.R @@ -102,6 +102,8 @@ test_that("SurvivalLogLogistic can recover known values", { test_that("Print method for SurvivalLogLogistic works as expected", { + withr::local_options(list(jmpost.double_eps = 0)) + expect_snapshot({ x <- SurvivalLogLogistic() print(x) diff --git a/tests/testthat/test-SurvivalQuantities.R b/tests/testthat/test-SurvivalQuantities.R index f58d2b55..197d92a0 100644 --- a/tests/testthat/test-SurvivalQuantities.R +++ b/tests/testthat/test-SurvivalQuantities.R @@ -142,6 +142,8 @@ test_that("autoplot.SurvivalSamples works as expected", { test_that("SurvivalQuantities print method works as expected", { + withr::local_options(list(jmpost.double_eps = 0)) + expect_snapshot({ set.seed(3219) subjectgroups <- list( diff --git a/tests/testthat/test-SurvivalWeibullPH.R b/tests/testthat/test-SurvivalWeibullPH.R index d3d3de56..a15f170e 100644 --- a/tests/testthat/test-SurvivalWeibullPH.R +++ b/tests/testthat/test-SurvivalWeibullPH.R @@ -1,4 +1,6 @@ test_that("Print method for SurvivalWeibullPH works as expected", { + withr::local_options(list(jmpost.double_eps = 0)) + expect_snapshot({ x <- SurvivalWeibullPH() print(x) diff --git a/tests/testthat/test-brms_compare.R b/tests/testthat/test-brms_compare.R index 13318f7c..ed235af0 100644 --- a/tests/testthat/test-brms_compare.R +++ b/tests/testthat/test-brms_compare.R @@ -421,6 +421,7 @@ test_that("jmpost and brms get similar horseshoe estimates for survival models", j_shrinkage <- shrinkage(mp) |> posterior::as_draws_df() |> + dplyr::as_tibble() |> dplyr::select(ones, cov1, cov2, cov3) expect_equal( From cef0a6e476990a7a0f0b7f5bda07bd0864ba9021 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 13:04:58 +0800 Subject: [PATCH 17/30] some fixes --- .../testthat/_snaps/LongitudinalSteinFojo.md | 2 -- tests/testthat/test-LongitudinalClaretBruno.R | 2 ++ tests/testthat/test-LongitudinalGSF.R | 18 ++++------ tests/testthat/test-LongitudinalSteinFojo.R | 19 +++++----- tests/testthat/test-SurvivalExponential.R | 36 ++++++++----------- 5 files changed, 31 insertions(+), 46 deletions(-) diff --git a/tests/testthat/_snaps/LongitudinalSteinFojo.md b/tests/testthat/_snaps/LongitudinalSteinFojo.md index 715db37a..7ebfa542 100644 --- a/tests/testthat/_snaps/LongitudinalSteinFojo.md +++ b/tests/testthat/_snaps/LongitudinalSteinFojo.md @@ -1,7 +1,6 @@ # Print method for LongitudinalSteinFojo works as expected Code - x <- LongitudinalSteinFojo() print(x) Output @@ -21,7 +20,6 @@ --- Code - x <- LongitudinalSteinFojo(sigma = prior_normal(0, 1), mu_kg = prior_gamma(2, 1)) print(x) Output diff --git a/tests/testthat/test-LongitudinalClaretBruno.R b/tests/testthat/test-LongitudinalClaretBruno.R index 23e3518c..c52d2009 100644 --- a/tests/testthat/test-LongitudinalClaretBruno.R +++ b/tests/testthat/test-LongitudinalClaretBruno.R @@ -305,6 +305,8 @@ test_that("Quantity models pass the parser", { test_that("Can generate valid initial values", { + withr::local_options(list(jmpost.double_eps = 0)) + pars <- c( "lm_clbr_omega_b", "lm_clbr_omega_g", diff --git a/tests/testthat/test-LongitudinalGSF.R b/tests/testthat/test-LongitudinalGSF.R index 783df469..69f79a71 100644 --- a/tests/testthat/test-LongitudinalGSF.R +++ b/tests/testthat/test-LongitudinalGSF.R @@ -9,18 +9,14 @@ test_that("LongitudinalGSF works as expected with default arguments", { test_that("Print method for LongitudinalGSF works as expected", { withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot({ - x <- LongitudinalGSF() - print(x) - }) + x <- LongitudinalGSF() + expect_snapshot(print(x)) - expect_snapshot({ - x <- LongitudinalGSF( - sigma = prior_normal(0, 1), - mu_kg = prior_gamma(2, 1) - ) - print(x) - }) + x <- LongitudinalGSF( + sigma = prior_normal(0, 1), + mu_kg = prior_gamma(2, 1) + ) + expect_snapshot(print(x)) }) diff --git a/tests/testthat/test-LongitudinalSteinFojo.R b/tests/testthat/test-LongitudinalSteinFojo.R index 2d513942..e1463a85 100644 --- a/tests/testthat/test-LongitudinalSteinFojo.R +++ b/tests/testthat/test-LongitudinalSteinFojo.R @@ -6,18 +6,15 @@ test_that("LongitudinalSteinFojo works as expected with default arguments", { test_that("Print method for LongitudinalSteinFojo works as expected", { withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot({ - x <- LongitudinalSteinFojo() - print(x) - }) - expect_snapshot({ - x <- LongitudinalSteinFojo( - sigma = prior_normal(0, 1), - mu_kg = prior_gamma(2, 1) - ) - print(x) - }) + x <- LongitudinalSteinFojo() + expect_snapshot(print(x)) + + x <- LongitudinalSteinFojo( + sigma = prior_normal(0, 1), + mu_kg = prior_gamma(2, 1) + ) + expect_snapshot(print(x)) }) diff --git a/tests/testthat/test-SurvivalExponential.R b/tests/testthat/test-SurvivalExponential.R index 38fe3ebb..2c9a5852 100644 --- a/tests/testthat/test-SurvivalExponential.R +++ b/tests/testthat/test-SurvivalExponential.R @@ -130,36 +130,28 @@ test_that("SurvivalExponential can recover true parameter (including covariates) test_that("Print method for SurvivalExponential works as expected", { withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot({ - x <- SurvivalExponential() - print(x) - }) + x <- SurvivalExponential() + expect_snapshot(print(x)) - expect_snapshot({ - x <- SurvivalExponential(beta = prior_gamma(3, 4)) - print(x) - }) + x <- SurvivalExponential(beta = prior_gamma(3, 4)) + expect_snapshot(print(x)) - expect_snapshot({ - x <- SurvivalExponential(lambda = prior_const(1)) - print(x) - }) + x <- SurvivalExponential(lambda = prior_const(1)) + expect_snapshot(print(x)) }) test_that("Different priors for the beta components are possible", { + withr::local_options(list(jmpost.double_eps = 0)) + # Same iid prior for all beta components: - expect_snapshot({ - x <- SurvivalExponential(beta = prior_normal(0, 1)) - print(x) - }) + x <- SurvivalExponential(beta = prior_normal(0, 1)) + expect_snapshot(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) - }) + x <- SurvivalExponential( + beta = prior_normal_vector(c(0, 1, 2), c(1, 2, 3)) + ) + expect_snapshot(print(x)) skip_if_not(is_full_test()) From 355db42ad0b909ef5ebe0bc4e49c267f8e4c3dec Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 13:08:19 +0800 Subject: [PATCH 18/30] fix prior_const printing --- R/Prior.R | 2 +- tests/testthat/test-Prior.R | 1 + tests/testthat/test-SurvivalExponential.R | 8 ++++++++ 3 files changed, 10 insertions(+), 1 deletion(-) diff --git a/R/Prior.R b/R/Prior.R index 26bf86b5..a487af83 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -215,7 +215,7 @@ as.character.Prior <- function(x, ...) { glue::glue, append(x@display, parameters_rounded) ) - display_limits <- render_stan_limits(x@limits) + display_limits <- if (x@.is_const) "" else render_stan_limits(x@limits) if ( display_limits != "" && display_string != "" && diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 2ff65062..cc844eb0 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -145,6 +145,7 @@ test_that("prior_const() works as expected", { x <- prior_const(1) expect_equal(initialValues(x), 1) + expect_equal(as.character(set_limits(x, lower = 0)), "const(value = 1)") expect_equal( as_stan_list(x, name = "sm_exp_lambda"), list(prior_const_sm_exp_lambda = 1) diff --git a/tests/testthat/test-SurvivalExponential.R b/tests/testthat/test-SurvivalExponential.R index 2c9a5852..8dcaeb15 100644 --- a/tests/testthat/test-SurvivalExponential.R +++ b/tests/testthat/test-SurvivalExponential.R @@ -60,6 +60,14 @@ test_that("SurvivalExponential() can fix lambda with prior_const()", { expect_stan_syntax(as.StanModule(beta_const)) }) +test_that("SurvivalExponential() does not print truncation for prior_const()", { + x <- SurvivalExponential(lambda = prior_const(1)) + expect_equal( + as.character(x@parameters@parameters[[1]]), + "sm_exp_lambda = const(value = 1)" + ) +}) + test_that("SurvivalExponential can recover true parameter (including covariates)", { skip_if_not(is_full_test()) From ea9007bab5150ba4fd271c01f66358747c367fa3 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 13:38:33 +0800 Subject: [PATCH 19/30] prior-aware dispatch for render_stan_limits(), so T[0, ] is omitted for: prior_cauchy(), prior_gamma(), prior_invgamma(), prior_loglogistic(), prior_lognormal() --- R/Prior.R | 51 ++++++++++++++----- R/settings.R | 13 ++--- _pkgdown.yml | 3 -- man/Prior-class.Rd | 5 ++ man/jmpost-settings.Rd | 13 ++--- man/render_stan_limits.Rd | 5 +- tests/testthat/_snaps/JointModel.md | 38 +++++++------- .../_snaps/LongitudinalClaretBruno.md | 18 +++---- tests/testthat/_snaps/LongitudinalGSF.md | 20 ++++---- .../_snaps/LongitudinalRandomSlope.md | 8 +-- .../testthat/_snaps/LongitudinalSteinFojo.md | 14 ++--- tests/testthat/_snaps/SurvivalGamma.md | 8 +-- tests/testthat/_snaps/SurvivalWeibullPH.md | 8 +-- tests/testthat/test-Prior.R | 40 ++++++++++++++- 14 files changed, 152 insertions(+), 92 deletions(-) diff --git a/R/Prior.R b/R/Prior.R index a487af83..be222ea7 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -39,6 +39,7 @@ NULL #' @slot display (`string`)\cr See arguments. #' @slot sample (`function`)\cr See arguments. #' @slot limits (`numeric`)\cr See arguments. +#' @slot .omit_zero_lower_truncation (`logical`)\cr See arguments. #' @slot .allow_vectors (`logical`)\cr See arguments. #' @slot .is_const (`logical`)\cr See arguments. #' @@ -59,6 +60,7 @@ NULL "validation" = "list", "sample" = "function", "limits" = "numeric", + ".omit_zero_lower_truncation" = "logical", ".allow_vectors" = "logical", ".is_const" = "logical" ) @@ -87,7 +89,9 @@ NULL #' @typed sample: function #' a function to sample from the prior distribution. #' @typed limits: numeric -#' the lower and upper limits for a truncated distribution +#' the lower and upper limits for a truncated distribution. +#' @typed .omit_zero_lower_truncation: flag +#' whether to omit a lower-zero truncation adjustment. #' @typed .allow_vectors: flag #' whether to allow vector parameters. #' @typed .is_const: flag @@ -105,6 +109,7 @@ Prior <- function( repr_transformed_parameters = "", repr_generated_quantities = "", limits = c(-Inf, Inf), + .omit_zero_lower_truncation = FALSE, .allow_vectors = FALSE, .is_const = FALSE ) { @@ -120,6 +125,7 @@ Prior <- function( validation = validation, sample = sample, limits = limits, + .omit_zero_lower_truncation = .omit_zero_lower_truncation, .allow_vectors = .allow_vectors, .is_const = .is_const ) @@ -215,7 +221,7 @@ as.character.Prior <- function(x, ...) { glue::glue, append(x@display, parameters_rounded) ) - display_limits <- if (x@.is_const) "" else render_stan_limits(x@limits) + display_limits <- if (x@.is_const) "" else render_stan_limits(x) if ( display_limits != "" && display_string != "" && @@ -230,12 +236,27 @@ as.character.Prior <- function(x, ...) { #' Creates Stan Syntax for Truncated distributions #' @description #' This function creates the Stan syntax for truncated distributions -#' @typed limits: numeric -#' the lower and upper limits for a truncated distribution +#' @typed object: Prior | numeric +#' prior or lower and upper limits for a truncated distribution. #' @keywords internal #' @typedreturn character #' the Stan syntax for truncated distributions -render_stan_limits <- function(limits) { +render_stan_limits <- function(object) { + UseMethod("render_stan_limits") +} + +render_stan_limits.Prior <- function(object) { + if ( + object@.omit_zero_lower_truncation && + identical(object@limits, c(0, Inf)) + ) { + return("") + } + render_stan_limits(object@limits) +} + +render_stan_limits.numeric <- function(object) { + limits <- object l_bound <- if (limits[[1]] > -Inf) limits[[1]] else "" u_bound <- if (limits[[2]] < Inf) limits[[2]] else "" string <- "" @@ -292,7 +313,7 @@ as.StanModule.Prior <- function(object, name, size = 1, ...) { paste0(" ", x, collapse = "\n") } trunctation <- if (object@repr_model != "") { - paste0(render_stan_limits(object@limits), ";") + paste0(render_stan_limits(object), ";") } else { "" } @@ -636,7 +657,8 @@ prior_cauchy <- function(mu, sigma) { validation = list( mu = is.numeric, sigma = \(x) x > 0 - ) + ), + .omit_zero_lower_truncation = TRUE ) } @@ -664,7 +686,8 @@ prior_gamma <- function(alpha, beta) { validation = list( alpha = \(x) x > 0, beta = \(x) x > 0 - ) + ), + .omit_zero_lower_truncation = TRUE ) } @@ -691,7 +714,8 @@ prior_lognormal <- function(mu, sigma) { validation = list( mu = is.numeric, sigma = \(x) x > 0 - ) + ), + .omit_zero_lower_truncation = TRUE ) } @@ -744,7 +768,8 @@ prior_init_only <- function(dist) { }, centre = dist@centre, validation = list(), - limits = dist@limits + limits = dist@limits, + .omit_zero_lower_truncation = dist@.omit_zero_lower_truncation ) } @@ -876,7 +901,8 @@ prior_loglogistic <- function(alpha, beta) { validation = list( alpha = \(x) x > 0, beta = \(x) x > 0 - ) + ), + .omit_zero_lower_truncation = TRUE ) } @@ -907,7 +933,8 @@ prior_invgamma <- function(alpha, beta) { validation = list( alpha = \(x) x > 0, beta = \(x) x > 0 - ) + ), + .omit_zero_lower_truncation = TRUE ) } diff --git a/R/settings.R b/R/settings.R index 555c21a8..373f4adb 100644 --- a/R/settings.R +++ b/R/settings.R @@ -44,14 +44,11 @@ #' #' ## `jmpost.double_eps` and `jmpost.double_neg_eps` #' -#' Default = `sqrt(.Machine$double.eps)` and `sqrt(.Machine$double.neg.eps)` +#' Default = `0` 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. +#' These are small floating-point numbers `x` which are used in parameter +#' definitions as lower boundaries to prevent numerical issues with the MCMC +#' sampler in the beginning of the sampling process. #' #' @examples #' \dontrun{ @@ -70,7 +67,7 @@ set_options <- function() { jmpost.cache_dir = cache_dir, jmpost.prior_shrinkage = 0.5, jmpost.gauss_quad_n = 15, - jmpost.double_eps = sqrt(.Machine$double.eps), + jmpost.double_eps = 0, jmpost.double_neg_eps = sqrt(.Machine$double.neg.eps) ) for (opt in names(jmpost_opts)) { diff --git a/_pkgdown.yml b/_pkgdown.yml index 2b409885..6bb901bc 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -211,14 +211,11 @@ reference: contents: - STAN_BLOCKS - - title: Data Sets contents: - os_data - tumor_data - - articles: - title: Articles navbar: ~ diff --git a/man/Prior-class.Rd b/man/Prior-class.Rd index b04bc554..95693e17 100644 --- a/man/Prior-class.Rd +++ b/man/Prior-class.Rd @@ -19,6 +19,7 @@ Prior( repr_transformed_parameters = "", repr_generated_quantities = "", limits = c(-Inf, Inf), + .omit_zero_lower_truncation = FALSE, .allow_vectors = FALSE, .is_const = FALSE ) @@ -47,6 +48,8 @@ the same names as the \code{paramaters} slot.} \item{limits}{(\code{numeric}) the lower and upper limits for a truncated distribution} +\item{.omit_zero_lower_truncation}{(\code{flag}) whether to omit a lower-zero truncation adjustment.} + \item{.allow_vectors}{(\code{flag}) whether to allow vector parameters.} \item{.is_const}{(\code{flag}) whether this prior fixes the parameter at a constant value.} @@ -79,6 +82,8 @@ Specifies the prior distribution in a Stan Model \item{\code{limits}}{(\code{numeric})\cr See arguments.} +\item{\code{.omit_zero_lower_truncation}}{(\code{logical})\cr See arguments.} + \item{\code{.allow_vectors}}{(\code{logical})\cr See arguments.} \item{\code{.is_const}}{(\code{logical})\cr See arguments.} diff --git a/man/jmpost-settings.Rd b/man/jmpost-settings.Rd index 5557ab6d..0dbb5527 100644 --- a/man/jmpost-settings.Rd +++ b/man/jmpost-settings.Rd @@ -52,14 +52,11 @@ 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. +Default = \code{0} and \code{sqrt(.Machine$double.neg.eps)} + +These are small floating-point numbers \code{x} which are used in parameter +definitions as lower boundaries to prevent numerical issues with the MCMC +sampler in the beginning of the sampling process. } } \examples{ diff --git a/man/render_stan_limits.Rd b/man/render_stan_limits.Rd index 3428e7c4..0ca4d492 100644 --- a/man/render_stan_limits.Rd +++ b/man/render_stan_limits.Rd @@ -4,10 +4,11 @@ \alias{render_stan_limits} \title{Creates Stan Syntax for Truncated distributions} \usage{ -render_stan_limits(limits) +render_stan_limits(object) } \arguments{ -\item{limits}{(\code{numeric}) the lower and upper limits for a truncated distribution} +\item{object}{a \code{\link{Prior}} or numeric lower/upper limits for a +truncated distribution} } \value{ (\code{character}) the Stan syntax for truncated distributions diff --git a/tests/testthat/_snaps/JointModel.md b/tests/testthat/_snaps/JointModel.md index cf991fda..df22db35 100644 --- a/tests/testthat/_snaps/JointModel.md +++ b/tests/testthat/_snaps/JointModel.md @@ -73,16 +73,16 @@ Survival: Weibull-PH Survival Model with parameters: - sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) T[0, ] - sm_weibull_ph_gamma ~ gamma(alpha = 2, beta = 0.5) T[0, ] + sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) + sm_weibull_ph_gamma ~ gamma(alpha = 2, beta = 0.5) beta_os_cov ~ normal(mu = 0, sigma = 2) Longitudinal: Random Slope Longitudinal Model (additive error) with parameters: lm_rs_intercept ~ normal(mu = 30, sigma = 10) lm_rs_slope_mu ~ normal(mu = 1, sigma = 3) - lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) T[0, ] - lm_rs_sigma ~ lognormal(mu = 0, sigma = 1.5) T[0, ] + lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) + lm_rs_sigma ~ lognormal(mu = 0, sigma = 1.5) lm_rs_ind_rnd_slope ~ Link: @@ -101,16 +101,16 @@ Survival: Weibull-PH Survival Model with parameters: - sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) T[0, ] - sm_weibull_ph_gamma ~ gamma(alpha = 2, beta = 0.5) T[0, ] + sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) + sm_weibull_ph_gamma ~ gamma(alpha = 2, beta = 0.5) beta_os_cov ~ normal(mu = 0, sigma = 2) Longitudinal: Random Slope Longitudinal Model (additive error) with parameters: lm_rs_intercept ~ normal(mu = 30, sigma = 10) lm_rs_slope_mu ~ normal(mu = 1, sigma = 3) - lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) T[0, ] - lm_rs_sigma ~ lognormal(mu = 0, sigma = 1.5) T[0, ] + lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) + lm_rs_sigma ~ lognormal(mu = 0, sigma = 1.5) lm_rs_ind_rnd_slope ~ Link: @@ -129,8 +129,8 @@ Survival: Weibull-PH Survival Model with parameters: - sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) T[0, ] - sm_weibull_ph_gamma ~ gamma(alpha = 2, beta = 0.5) T[0, ] + sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) + sm_weibull_ph_gamma ~ gamma(alpha = 2, beta = 0.5) beta_os_cov ~ normal(mu = 0, sigma = 2) Longitudinal: @@ -157,11 +157,11 @@ lm_gsf_mu_ks ~ normal(mu = -0.69315, sigma = 1) lm_gsf_mu_kg ~ normal(mu = -1.20397, sigma = 1) lm_gsf_mu_phi ~ normal(mu = 0, sigma = 1) - lm_gsf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_omega_phi ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_sigma ~ lognormal(mu = -2.30259, sigma = 1) T[0, ] + lm_gsf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_omega_phi ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_sigma ~ lognormal(mu = -2.30259, sigma = 1) lm_gsf_eta_tilde_bsld ~ std_normal() lm_gsf_eta_tilde_ks ~ std_normal() lm_gsf_eta_tilde_kg ~ std_normal() @@ -182,16 +182,16 @@ Survival: Weibull-PH Survival Model with parameters: - sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) T[0, ] - sm_weibull_ph_gamma ~ gamma(alpha = 2, beta = 0.5) T[0, ] + sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) + sm_weibull_ph_gamma ~ gamma(alpha = 2, beta = 0.5) beta_os_cov ~ normal(mu = 0, sigma = 2) Longitudinal: Random Slope Longitudinal Model (additive error) with parameters: lm_rs_intercept ~ normal(mu = 30, sigma = 10) lm_rs_slope_mu ~ normal(mu = 1, sigma = 3) - lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) T[0, ] - lm_rs_sigma ~ lognormal(mu = 0, sigma = 1.5) T[0, ] + lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) + lm_rs_sigma ~ lognormal(mu = 0, sigma = 1.5) lm_rs_ind_rnd_slope ~ Link: diff --git a/tests/testthat/_snaps/LongitudinalClaretBruno.md b/tests/testthat/_snaps/LongitudinalClaretBruno.md index 1942da06..29601911 100644 --- a/tests/testthat/_snaps/LongitudinalClaretBruno.md +++ b/tests/testthat/_snaps/LongitudinalClaretBruno.md @@ -10,11 +10,11 @@ lm_clbr_mu_g ~ normal(mu = 0, sigma = 0.5) lm_clbr_mu_c ~ normal(mu = -0.91629, sigma = 0.5) lm_clbr_mu_p ~ normal(mu = 0.69315, sigma = 0.5) - lm_clbr_omega_b ~ lognormal(mu = -1.60944, sigma = 0.5) T[0, ] - lm_clbr_omega_g ~ lognormal(mu = -1.60944, sigma = 0.5) T[0, ] - lm_clbr_omega_c ~ lognormal(mu = -1.60944, sigma = 0.5) T[0, ] - lm_clbr_omega_p ~ lognormal(mu = -1.60944, sigma = 0.5) T[0, ] - lm_clbr_sigma ~ lognormal(mu = -2.30259, sigma = 0.5) T[0, ] + lm_clbr_omega_b ~ lognormal(mu = -1.60944, sigma = 0.5) + lm_clbr_omega_g ~ lognormal(mu = -1.60944, sigma = 0.5) + lm_clbr_omega_c ~ lognormal(mu = -1.60944, sigma = 0.5) + lm_clbr_omega_p ~ lognormal(mu = -1.60944, sigma = 0.5) + lm_clbr_sigma ~ lognormal(mu = -2.30259, sigma = 0.5) lm_clbr_eta_b ~ std_normal() lm_clbr_eta_g ~ std_normal() lm_clbr_eta_c ~ std_normal() @@ -34,10 +34,10 @@ lm_clbr_mu_g ~ gamma(alpha = 2, beta = 1) lm_clbr_mu_c ~ normal(mu = -0.91629, sigma = 0.5) lm_clbr_mu_p ~ normal(mu = 0.69315, sigma = 0.5) - lm_clbr_omega_b ~ lognormal(mu = -1.60944, sigma = 0.5) T[0, ] - lm_clbr_omega_g ~ lognormal(mu = -1.60944, sigma = 0.5) T[0, ] - lm_clbr_omega_c ~ lognormal(mu = -1.60944, sigma = 0.5) T[0, ] - lm_clbr_omega_p ~ lognormal(mu = -1.60944, sigma = 0.5) T[0, ] + lm_clbr_omega_b ~ lognormal(mu = -1.60944, sigma = 0.5) + lm_clbr_omega_g ~ lognormal(mu = -1.60944, sigma = 0.5) + lm_clbr_omega_c ~ lognormal(mu = -1.60944, sigma = 0.5) + lm_clbr_omega_p ~ lognormal(mu = -1.60944, sigma = 0.5) lm_clbr_sigma ~ normal(mu = 0, sigma = 1) T[0, ] lm_clbr_eta_b ~ std_normal() lm_clbr_eta_g ~ std_normal() diff --git a/tests/testthat/_snaps/LongitudinalGSF.md b/tests/testthat/_snaps/LongitudinalGSF.md index a531a6b2..51cb48cd 100644 --- a/tests/testthat/_snaps/LongitudinalGSF.md +++ b/tests/testthat/_snaps/LongitudinalGSF.md @@ -1,7 +1,6 @@ # Print method for LongitudinalGSF works as expected Code - x <- LongitudinalGSF() print(x) Output @@ -10,11 +9,11 @@ lm_gsf_mu_ks ~ normal(mu = -0.69315, sigma = 1) lm_gsf_mu_kg ~ normal(mu = -1.20397, sigma = 1) lm_gsf_mu_phi ~ normal(mu = 0, sigma = 1) - lm_gsf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_omega_phi ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_sigma ~ lognormal(mu = -2.30259, sigma = 1) T[0, ] + lm_gsf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_omega_phi ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_sigma ~ lognormal(mu = -2.30259, sigma = 1) lm_gsf_eta_tilde_bsld ~ std_normal() lm_gsf_eta_tilde_ks ~ std_normal() lm_gsf_eta_tilde_kg ~ std_normal() @@ -24,7 +23,6 @@ --- Code - x <- LongitudinalGSF(sigma = prior_normal(0, 1), mu_kg = prior_gamma(2, 1)) print(x) Output @@ -33,10 +31,10 @@ lm_gsf_mu_ks ~ normal(mu = -0.69315, sigma = 1) lm_gsf_mu_kg ~ gamma(alpha = 2, beta = 1) lm_gsf_mu_phi ~ normal(mu = 0, sigma = 1) - lm_gsf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_gsf_omega_phi ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] + lm_gsf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1) + lm_gsf_omega_phi ~ lognormal(mu = -1.60944, sigma = 1) lm_gsf_sigma ~ normal(mu = 0, sigma = 1) T[0, ] lm_gsf_eta_tilde_bsld ~ std_normal() lm_gsf_eta_tilde_ks ~ std_normal() diff --git a/tests/testthat/_snaps/LongitudinalRandomSlope.md b/tests/testthat/_snaps/LongitudinalRandomSlope.md index cdbb06b8..93900a76 100644 --- a/tests/testthat/_snaps/LongitudinalRandomSlope.md +++ b/tests/testthat/_snaps/LongitudinalRandomSlope.md @@ -8,8 +8,8 @@ Random Slope Longitudinal Model (additive error) with parameters: lm_rs_intercept ~ normal(mu = 30, sigma = 10) lm_rs_slope_mu ~ normal(mu = 1, sigma = 3) - lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) T[0, ] - lm_rs_sigma ~ lognormal(mu = 0, sigma = 1.5) T[0, ] + lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) + lm_rs_sigma ~ lognormal(mu = 0, sigma = 1.5) lm_rs_ind_rnd_slope ~ @@ -24,8 +24,8 @@ Random Slope Longitudinal Model (additive error) with parameters: lm_rs_intercept ~ normal(mu = 0, sigma = 1) lm_rs_slope_mu ~ normal(mu = 1, sigma = 3) - lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) T[0, ] - lm_rs_sigma ~ gamma(alpha = 2, beta = 1) T[0, ] + lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) + lm_rs_sigma ~ gamma(alpha = 2, beta = 1) lm_rs_ind_rnd_slope ~ diff --git a/tests/testthat/_snaps/LongitudinalSteinFojo.md b/tests/testthat/_snaps/LongitudinalSteinFojo.md index 7ebfa542..85eccd07 100644 --- a/tests/testthat/_snaps/LongitudinalSteinFojo.md +++ b/tests/testthat/_snaps/LongitudinalSteinFojo.md @@ -8,10 +8,10 @@ lm_sf_mu_bsld ~ normal(mu = 4.09434, sigma = 1) lm_sf_mu_ks ~ normal(mu = -0.69315, sigma = 1) lm_sf_mu_kg ~ normal(mu = -1.20397, sigma = 1) - lm_sf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_sf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_sf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_sf_sigma ~ lognormal(mu = -2.30259, sigma = 1) T[0, ] + lm_sf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1) + lm_sf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1) + lm_sf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1) + lm_sf_sigma ~ lognormal(mu = -2.30259, sigma = 1) lm_sf_eta_tilde_bsld ~ std_normal() lm_sf_eta_tilde_ks ~ std_normal() lm_sf_eta_tilde_kg ~ std_normal() @@ -27,9 +27,9 @@ lm_sf_mu_bsld ~ normal(mu = 4.09434, sigma = 1) lm_sf_mu_ks ~ normal(mu = -0.69315, sigma = 1) lm_sf_mu_kg ~ gamma(alpha = 2, beta = 1) - lm_sf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_sf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] - lm_sf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1) T[0, ] + lm_sf_omega_bsld ~ lognormal(mu = -1.60944, sigma = 1) + lm_sf_omega_ks ~ lognormal(mu = -1.60944, sigma = 1) + lm_sf_omega_kg ~ lognormal(mu = -1.60944, sigma = 1) lm_sf_sigma ~ normal(mu = 0, sigma = 1) T[0, ] lm_sf_eta_tilde_bsld ~ std_normal() lm_sf_eta_tilde_ks ~ std_normal() diff --git a/tests/testthat/_snaps/SurvivalGamma.md b/tests/testthat/_snaps/SurvivalGamma.md index ca86ce0b..2cb71967 100644 --- a/tests/testthat/_snaps/SurvivalGamma.md +++ b/tests/testthat/_snaps/SurvivalGamma.md @@ -6,8 +6,8 @@ Output Gamma Survival Model with parameters: - sm_gamma_k ~ gamma(alpha = 2, beta = 0.5) T[0, ] - sm_gamma_theta ~ gamma(alpha = 2, beta = 0.5) T[0, ] + sm_gamma_k ~ gamma(alpha = 2, beta = 0.5) + sm_gamma_theta ~ gamma(alpha = 2, beta = 0.5) beta_os_cov ~ normal(mu = 0, sigma = 2) @@ -19,8 +19,8 @@ Output Gamma Survival Model with parameters: - sm_gamma_k ~ gamma(alpha = 3, beta = 4) T[0, ] - sm_gamma_theta ~ cauchy(mu = 0, sigma = 1) T[0, ] + sm_gamma_k ~ gamma(alpha = 3, beta = 4) + sm_gamma_theta ~ cauchy(mu = 0, sigma = 1) beta_os_cov ~ normal(mu = 0, sigma = 2) diff --git a/tests/testthat/_snaps/SurvivalWeibullPH.md b/tests/testthat/_snaps/SurvivalWeibullPH.md index 97e832c3..326b229e 100644 --- a/tests/testthat/_snaps/SurvivalWeibullPH.md +++ b/tests/testthat/_snaps/SurvivalWeibullPH.md @@ -6,8 +6,8 @@ Output Weibull-PH Survival Model with parameters: - sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) T[0, ] - sm_weibull_ph_gamma ~ gamma(alpha = 2, beta = 0.5) T[0, ] + sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) + sm_weibull_ph_gamma ~ gamma(alpha = 2, beta = 0.5) beta_os_cov ~ normal(mu = 0, sigma = 2) @@ -19,8 +19,8 @@ Output Weibull-PH Survival Model with parameters: - sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) T[0, ] - sm_weibull_ph_gamma ~ cauchy(mu = 0, sigma = 1) T[0, ] + sm_weibull_ph_lambda ~ gamma(alpha = 2, beta = 0.5) + sm_weibull_ph_gamma ~ cauchy(mu = 0, sigma = 1) beta_os_cov ~ gamma(alpha = 3, beta = 4) diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index cc844eb0..3cba84fb 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -246,7 +246,7 @@ test_that("Limits work as expected", { expect_true(all(ivs > 0)) expect_equal( as.StanModule(x, name = "tim")@model, - " tim ~ cauchy(prior_mu_tim, prior_sigma_tim) T[0, ];" + " tim ~ cauchy(prior_mu_tim, prior_sigma_tim);" ) ## Put an impossible constraint on the distribution @@ -260,6 +260,44 @@ test_that("Limits work as expected", { }) +test_that("redundant positive-support prior truncation is omitted", { + positive_priors <- list( + cauchy = prior_cauchy(0, 1), + gamma = prior_gamma(2, 1), + invgamma = prior_invgamma(2, 1), + loglogistic = prior_loglogistic(2, 3), + lognormal = prior_lognormal(0, 1) + ) + + for (prior in positive_priors) { + x <- set_limits(prior, lower = 0) + expect_false(grepl("T\\[0, \\]", as.character(x))) + expect_false(grepl( + "T\\[0, \\]", + as.StanModule(x, name = "theta")@model + )) + + x <- set_limits(prior, lower = 1e-8) + expect_true(grepl("T\\[1e-08, \\]", as.character(x))) + expect_true(grepl( + "T\\[1e-08, \\]", + as.StanModule(x, name = "theta")@model + )) + + x <- set_limits(prior, lower = 0, upper = 1) + expect_true(grepl("T\\[0, 1\\]", as.character(x))) + expect_true(grepl( + "T\\[0, 1\\]", + as.StanModule(x, name = "theta")@model + )) + } + + x <- set_limits(prior_normal(0, 1), lower = 0) + expect_true(grepl("T\\[0, \\]", as.character(x))) + expect_true(grepl("T\\[0, \\]", as.StanModule(x, name = "theta")@model)) +}) + + test_that("median(Prior) works as expected", { set.seed(2410) From 1dde042aeaba9fae4904f5f84fa8cf0aa89289b6 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 13:49:21 +0800 Subject: [PATCH 20/30] updates to tests --- tests/testthat/_snaps/SurvivalExponential.md | 5 ----- tests/testthat/test-JointModel.R | 4 ---- tests/testthat/test-LongitudinalClaretBruno.R | 4 ---- tests/testthat/test-LongitudinalGSF.R | 4 ---- tests/testthat/test-LongitudinalRandomSlope.R | 2 -- tests/testthat/test-LongitudinalSteinFojo.R | 2 -- tests/testthat/test-SurvivalExponential.R | 4 ---- tests/testthat/test-SurvivalGamma.R | 2 -- tests/testthat/test-SurvivalLoglogistic.R | 2 -- tests/testthat/test-SurvivalQuantities.R | 2 -- tests/testthat/test-SurvivalWeibullPH.R | 2 -- 11 files changed, 33 deletions(-) diff --git a/tests/testthat/_snaps/SurvivalExponential.md b/tests/testthat/_snaps/SurvivalExponential.md index 24805437..4d5b7d0d 100644 --- a/tests/testthat/_snaps/SurvivalExponential.md +++ b/tests/testthat/_snaps/SurvivalExponential.md @@ -1,7 +1,6 @@ # Print method for SurvivalExponential works as expected Code - x <- SurvivalExponential() print(x) Output @@ -13,7 +12,6 @@ --- Code - x <- SurvivalExponential(beta = prior_gamma(3, 4)) print(x) Output @@ -25,7 +23,6 @@ --- Code - x <- SurvivalExponential(lambda = prior_const(1)) print(x) Output @@ -37,7 +34,6 @@ # Different priors for the beta components are possible Code - x <- SurvivalExponential(beta = prior_normal(0, 1)) print(x) Output @@ -49,7 +45,6 @@ --- Code - x <- SurvivalExponential(beta = prior_normal_vector(c(0, 1, 2), c(1, 2, 3))) print(x) Output diff --git a/tests/testthat/test-JointModel.R b/tests/testthat/test-JointModel.R index 0e6f9583..e1a77dfa 100644 --- a/tests/testthat/test-JointModel.R +++ b/tests/testthat/test-JointModel.R @@ -46,8 +46,6 @@ snapshot_joint_model_parameter_declarations <- function(object) { test_that("JointModel snapshots assembled parameter declarations", { - withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot( snapshot_joint_model_parameter_declarations(JointModel( longitudinal = LongitudinalRandomSlope(), @@ -73,8 +71,6 @@ test_that("JointModel snapshots assembled parameter declarations", { test_that("JointModel print method works as expected", { - withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot({ x <- JointModel( longitudinal = LongitudinalRandomSlope(), diff --git a/tests/testthat/test-LongitudinalClaretBruno.R b/tests/testthat/test-LongitudinalClaretBruno.R index c52d2009..67f6c2d7 100644 --- a/tests/testthat/test-LongitudinalClaretBruno.R +++ b/tests/testthat/test-LongitudinalClaretBruno.R @@ -5,8 +5,6 @@ test_that("LongitudinalClaretBruno works as expected with default arguments", { test_that("Print method for LongitudinalClaretBruno works as expected", { - withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot({ x <- LongitudinalClaretBruno() print(x) @@ -305,8 +303,6 @@ test_that("Quantity models pass the parser", { test_that("Can generate valid initial values", { - withr::local_options(list(jmpost.double_eps = 0)) - pars <- c( "lm_clbr_omega_b", "lm_clbr_omega_g", diff --git a/tests/testthat/test-LongitudinalGSF.R b/tests/testthat/test-LongitudinalGSF.R index 69f79a71..7377a89d 100644 --- a/tests/testthat/test-LongitudinalGSF.R +++ b/tests/testthat/test-LongitudinalGSF.R @@ -7,8 +7,6 @@ test_that("LongitudinalGSF works as expected with default arguments", { test_that("Print method for LongitudinalGSF works as expected", { - withr::local_options(list(jmpost.double_eps = 0)) - x <- LongitudinalGSF() expect_snapshot(print(x)) @@ -308,8 +306,6 @@ test_that("Quantity models pass the parser", { test_that("Can generate valid initial values", { - withr::local_options(list(jmpost.double_eps = 0)) - pars <- c( "lm_gsf_omega_bsld", "lm_gsf_omega_ks", diff --git a/tests/testthat/test-LongitudinalRandomSlope.R b/tests/testthat/test-LongitudinalRandomSlope.R index bfae33d3..92f769c7 100644 --- a/tests/testthat/test-LongitudinalRandomSlope.R +++ b/tests/testthat/test-LongitudinalRandomSlope.R @@ -1,6 +1,4 @@ test_that("Print method for LongitudinalRandomSlope works as expected", { - withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot({ x <- LongitudinalRandomSlope() print(x) diff --git a/tests/testthat/test-LongitudinalSteinFojo.R b/tests/testthat/test-LongitudinalSteinFojo.R index e1463a85..05f1f26e 100644 --- a/tests/testthat/test-LongitudinalSteinFojo.R +++ b/tests/testthat/test-LongitudinalSteinFojo.R @@ -5,8 +5,6 @@ test_that("LongitudinalSteinFojo works as expected with default arguments", { test_that("Print method for LongitudinalSteinFojo works as expected", { - withr::local_options(list(jmpost.double_eps = 0)) - x <- LongitudinalSteinFojo() expect_snapshot(print(x)) diff --git a/tests/testthat/test-SurvivalExponential.R b/tests/testthat/test-SurvivalExponential.R index 8dcaeb15..a5435754 100644 --- a/tests/testthat/test-SurvivalExponential.R +++ b/tests/testthat/test-SurvivalExponential.R @@ -136,8 +136,6 @@ test_that("SurvivalExponential can recover true parameter (including covariates) test_that("Print method for SurvivalExponential works as expected", { - withr::local_options(list(jmpost.double_eps = 0)) - x <- SurvivalExponential() expect_snapshot(print(x)) @@ -149,8 +147,6 @@ test_that("Print method for SurvivalExponential works as expected", { }) test_that("Different priors for the beta components are possible", { - withr::local_options(list(jmpost.double_eps = 0)) - # Same iid prior for all beta components: x <- SurvivalExponential(beta = prior_normal(0, 1)) expect_snapshot(print(x)) diff --git a/tests/testthat/test-SurvivalGamma.R b/tests/testthat/test-SurvivalGamma.R index febec990..ee86a157 100644 --- a/tests/testthat/test-SurvivalGamma.R +++ b/tests/testthat/test-SurvivalGamma.R @@ -1,6 +1,4 @@ test_that("Print method for SurvivalGamma works as expected", { - withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot({ x <- SurvivalGamma() print(x) diff --git a/tests/testthat/test-SurvivalLoglogistic.R b/tests/testthat/test-SurvivalLoglogistic.R index 5e10bb1f..1e9c26fb 100644 --- a/tests/testthat/test-SurvivalLoglogistic.R +++ b/tests/testthat/test-SurvivalLoglogistic.R @@ -102,8 +102,6 @@ test_that("SurvivalLogLogistic can recover known values", { test_that("Print method for SurvivalLogLogistic works as expected", { - withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot({ x <- SurvivalLogLogistic() print(x) diff --git a/tests/testthat/test-SurvivalQuantities.R b/tests/testthat/test-SurvivalQuantities.R index 197d92a0..f58d2b55 100644 --- a/tests/testthat/test-SurvivalQuantities.R +++ b/tests/testthat/test-SurvivalQuantities.R @@ -142,8 +142,6 @@ test_that("autoplot.SurvivalSamples works as expected", { test_that("SurvivalQuantities print method works as expected", { - withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot({ set.seed(3219) subjectgroups <- list( diff --git a/tests/testthat/test-SurvivalWeibullPH.R b/tests/testthat/test-SurvivalWeibullPH.R index a15f170e..d3d3de56 100644 --- a/tests/testthat/test-SurvivalWeibullPH.R +++ b/tests/testthat/test-SurvivalWeibullPH.R @@ -1,6 +1,4 @@ test_that("Print method for SurvivalWeibullPH works as expected", { - withr::local_options(list(jmpost.double_eps = 0)) - expect_snapshot({ x <- SurvivalWeibullPH() print(x) From 77f7733f1bf6bb634aa1e7987b9ab36849895d73 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 13:55:06 +0800 Subject: [PATCH 21/30] fix one test --- tests/testthat/test-jinjar_render.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-jinjar_render.R b/tests/testthat/test-jinjar_render.R index 91987737..2f33481a 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", - scales::scientific(getOption("jmpost.double_eps"), digits = 5) + getOption("jmpost.double_eps") ) expect_equal(observed, expected) }) From ff585b5e9801d6159cd502ced2c839adccbb90f6 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 13:58:56 +0800 Subject: [PATCH 22/30] use positive epsilon for the one test --- tests/testthat/helper-example_data.R | 7 ++++++- tests/testthat/test-misc_models.R | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/tests/testthat/helper-example_data.R b/tests/testthat/helper-example_data.R index 5f9d826c..7f372d53 100644 --- a/tests/testthat/helper-example_data.R +++ b/tests/testthat/helper-example_data.R @@ -1,4 +1,5 @@ -ensure_test_data_1 <- function() { +ensure_test_data_1 <- function(positive_epsilon = FALSE) { + assert_flag(positive_epsilon) set.seed(739) simjdat <- SimJointData( design = list( @@ -20,6 +21,10 @@ ensure_test_data_1 <- function() { .silent = TRUE ) + withr::local_options(list( + jmpost.double_eps = ifelse(positive_epsilon, 1e-10, 0) + )) + dat_os <- simjdat@survival dat_lm <- simjdat@longitudinal diff --git a/tests/testthat/test-misc_models.R b/tests/testthat/test-misc_models.R index 58fce01c..3128acd7 100644 --- a/tests/testthat/test-misc_models.R +++ b/tests/testthat/test-misc_models.R @@ -1,4 +1,4 @@ -test_data_1 <- ensure_test_data_1() +test_data_1 <- ensure_test_data_1(positive_epsilon = TRUE) test_that("Longitudinal Model doesn't print sampler rejection messages", { # These rejections typically happen when the sampler samples a From d9fe6d440fa20bdc449e25dccd01de7ba031d072 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 14:38:10 +0800 Subject: [PATCH 23/30] make test less rigid --- tests/testthat/test-simulate.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-simulate.R b/tests/testthat/test-simulate.R index bb87b448..138950fb 100644 --- a/tests/testthat/test-simulate.R +++ b/tests/testthat/test-simulate.R @@ -139,7 +139,8 @@ test_that("simulate works with lambda_censor", { joint_results, times = (1:10), lambda_censor = 1 / 6, - time_step = 0.5 + time_step = 0.5, + seed = 9382 ) expect_data_frame( results@survival, @@ -147,7 +148,7 @@ test_that("simulate works with lambda_censor", { ncols = 9 ) - expect_equal(mean(results@survival$event), 0.6601, tolerance = 0.001) + expect_equal(mean(results@survival$event), 0.6601, tolerance = 0.01) expect_number(mean(results@survival$time), lower = 1, upper = 3) expect_number( mean(results@survival$time[results@survival$event == 0]), From e7511690de1420285163a21df5d0ce697bd17b90 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 15:14:22 +0800 Subject: [PATCH 24/30] polish vignette --- vignettes/covariate-selection-horseshoe.Rmd | 30 ++++++++++++--------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/vignettes/covariate-selection-horseshoe.Rmd b/vignettes/covariate-selection-horseshoe.Rmd index 7a4eb3c5..e93d8274 100644 --- a/vignettes/covariate-selection-horseshoe.Rmd +++ b/vignettes/covariate-selection-horseshoe.Rmd @@ -60,7 +60,7 @@ arguments as follows: The local shrinkage parameters $\lambda_j$ let individual coefficients escape the global shrinkage $\tau$. The slab scale $s_c$ and slab degrees of freedom $\nu_c$ control how much very large coefficients are regularized. Smaller -`scale_global` values encode stronger prior belief that only a small number of +$s_{\tau}$ values encode stronger prior belief that only a small number of candidate covariates are relevant. A useful posterior diagnostic is the shrinkage factor @@ -107,14 +107,16 @@ sim_data <- SimJointData( lambda_cen = 1 / 9000, beta_cat = c( "A" = 0, - "B" = -0.1, - "C" = 0.5 + "B" = 0.3, + "C" = 0.7 ), beta_cont = 0.3 ) ) ``` +Then we prepare the corresponding objects: + ```{r prepare-data} os_data <- sim_data@survival long_data <- sim_data@longitudinal @@ -138,7 +140,7 @@ joint_data <- DataJoint( ) ``` -Before setting priors or interpreting the output, check the survival design +Before setting priors or interpreting the output, it is a good habit to check the survival design matrix. With `covariates()` we can access the column names, which map to the corresponding survival coefficients that receive the horseshoe prior. @@ -152,7 +154,7 @@ survival_covariates Now we define a joint model. The longitudinal model and the association link use standard priors. The survival model uses `prior_horseshoe()` for the vector of -survival covariate coefficients. +survival covariate coefficients: ```{r model-specification} joint_model <- JointModel( @@ -188,16 +190,13 @@ fit <- sampleStanModel( chains = 4, parallel_chains = 4, seed = 325, - refresh = 200 + refresh = 0 ) ``` -The iteration counts above are chosen to keep this vignette reasonably fast to -render. Refit with longer chains before making inferential claims. - ## Inspect Coefficients -After fitting, inspect the coefficient posterior and standard MCMC diagnostics. +After fitting, let's inspect the coefficient posterior and standard MCMC diagnostics. The survival covariate coefficients are stored in `beta_os_cov`. ```{r inspect-coefficients} @@ -227,7 +226,7 @@ closer to 0 have escaped shrinkage and are more strongly supported by the model. ## Extract and Plot Shrinkage Factors -Use `shrinkage()` to extract the posterior draws of $\kappa_j$. The returned +Let's use the `shrinkage()` function to extract the posterior draws of $\kappa_j$. The returned draws are named with the survival covariate names, so they can be plotted or summarised directly. @@ -250,10 +249,13 @@ mcmc_dens_overlay(shrinkage_draws) + ) ``` +We see that both coefficients for the categorical covariate have shrinkage factors near 1, while the continuous +covariate's shrinkage factor is a bit more distributed towards smaller values. + A second useful display is the median shrinkage factor with an interval for each -covariate. +covariate: -```{r shrinkage-intervals, fig.width=7, fig.height=4} +```{r shrinkage-intervals, fig.width=7, fig.height=4, message=FALSE} library(dplyr) library(ggplot2) @@ -316,3 +318,5 @@ Finally, the horseshoe prior is a shrinkage prior, not a replacement for model checking. Always inspect convergence diagnostics, posterior predictive fit, and sensitivity to reasonable prior choices before using the selected covariates for scientific conclusions. + +## References From bf1771c7ad943a163216fe97daa52291dd24b9f7 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 15:34:21 +0800 Subject: [PATCH 25/30] add initial values also for auxiliary prior parameters (needed for horseshoe prior) to avoid warnings about uninitialized parameters --- NAMESPACE | 2 ++ R/ParameterList.R | 32 +++++++++++++++----- R/Prior.R | 47 +++++++++++++++++++++++++++++ R/generics.R | 40 ++++++++++++++++++++++++ _pkgdown.yml | 2 ++ man/Prior-Getter-Methods.Rd | 14 +++++++++ man/Prior-class.Rd | 12 ++++++++ man/auxiliaryInitialValues.Rd | 22 ++++++++++++++ man/auxiliarySize.Rd | 22 ++++++++++++++ tests/testthat/test-Prior.R | 31 +++++++++++++++++++ tests/testthat/test-initialValues.R | 39 ++++++++++++++++++++++++ 11 files changed, 256 insertions(+), 7 deletions(-) create mode 100644 man/auxiliaryInitialValues.Rd create mode 100644 man/auxiliarySize.Rd diff --git a/NAMESPACE b/NAMESPACE index 2b3da074..7b08324b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -79,6 +79,8 @@ S3method(as_stan_list,QuantityGeneratorPopulation) S3method(as_stan_list,QuantityGeneratorPrediction) S3method(as_stan_list,QuantityGeneratorSubject) S3method(as_stan_list,default) +S3method(auxiliaryInitialValues,Prior) +S3method(auxiliarySize,Prior) S3method(autoplot,LongitudinalQuantities) S3method(autoplot,SurvivalQuantities) S3method(brierScore,SurvivalQuantities) diff --git a/R/ParameterList.R b/R/ParameterList.R index 8f30fd9b..f46b6dc0 100644 --- a/R/ParameterList.R +++ b/R/ParameterList.R @@ -159,10 +159,20 @@ initialValues.ParameterList <- function(object, n_chains, ...) { lapply( seq_len(n_chains), \(i) { - vals <- lapply(parameters, initialValues) - name <- vapply(parameters, names, character(1)) - names(vals) <- name - vals + vals <- lapply(parameters, function(parameter) { + c( + stats::setNames( + list(initialValues(parameter)), + names(parameter) + ), + auxiliaryInitialValues( + parameter@prior, + name = names(parameter), + size = size(parameter) + ) + ) + }) + unlist(vals, recursive = FALSE) } ) } @@ -172,9 +182,17 @@ initialValues.ParameterList <- function(object, n_chains, ...) { #' @export size.ParameterList <- function(object) { parameters <- Filter(\(x) !x@prior@.is_const, object@parameters) - x <- lapply(parameters, size) - names(x) <- vapply(parameters, names, character(1)) - return(x) + sizes <- lapply(parameters, function(parameter) { + c( + stats::setNames(list(size(parameter)), names(parameter)), + auxiliarySize( + parameter@prior, + name = names(parameter), + size = size(parameter) + ) + ) + }) + unlist(sizes, recursive = FALSE) } diff --git a/R/Prior.R b/R/Prior.R index be222ea7..eb718bfa 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -16,6 +16,8 @@ NULL #' a prior Distribution #' @typed name: character #' the name of the parameter the prior distribution is for +#' @typed size: "`numeric` or `character`" +#' the parameter size. #' @param ... Not Used. #' #' @name Prior-Shared @@ -34,6 +36,8 @@ NULL #' @slot repr_parameters (`string`)\cr See arguments. #' @slot repr_transformed_parameters (`string`)\cr See arguments. #' @slot repr_generated_quantities (`string`)\cr See arguments. +#' @slot auxiliary_initial_values (`function`)\cr See arguments. +#' @slot auxiliary_size (`function`)\cr See arguments. #' @slot centre (`numeric`)\cr See arguments. #' @slot validation (`list`)\cr See arguments. #' @slot display (`string`)\cr See arguments. @@ -56,6 +60,8 @@ NULL "repr_parameters" = "character", "repr_transformed_parameters" = "character", "repr_generated_quantities" = "character", + "auxiliary_initial_values" = "function", + "auxiliary_size" = "function", "centre" = "numeric", "validation" = "list", "sample" = "function", @@ -79,6 +85,12 @@ NULL #' the Stan code representation for the transformed parameters block. #' @typed repr_generated_quantities: string #' the Stan code representation for the generated quantities block. +#' @typed auxiliary_initial_values: function +#' a function that returns initial values for extra Stan parameters +#' introduced by the prior. +#' @typed auxiliary_size: function +#' a function that returns sizes for extra Stan parameters introduced +#' by the prior. #' @typed display: string #' the string to display when object is printed. #' @typed centre: numeric @@ -108,6 +120,8 @@ Prior <- function( repr_parameters = "", repr_transformed_parameters = "", repr_generated_quantities = "", + auxiliary_initial_values = \(name, size) list(), + auxiliary_size = \(name, size) list(), limits = c(-Inf, Inf), .omit_zero_lower_truncation = FALSE, .allow_vectors = FALSE, @@ -120,6 +134,8 @@ Prior <- function( repr_parameters = repr_parameters, repr_transformed_parameters = repr_transformed_parameters, repr_generated_quantities = repr_generated_quantities, + auxiliary_initial_values = auxiliary_initial_values, + auxiliary_size = auxiliary_size, centre = centre, display = display, validation = validation, @@ -440,6 +456,20 @@ initialValues.Prior <- function(object, ...) { } +#' @describeIn Prior-Getter-Methods The prior's auxiliary initial values +#' @export +auxiliaryInitialValues.Prior <- function(object, name, size, ...) { + object@auxiliary_initial_values(name = name, size = size) +} + + +#' @describeIn Prior-Getter-Methods The prior's auxiliary parameter sizes +#' @export +auxiliarySize.Prior <- function(object, name, size, ...) { + object@auxiliary_size(name = name, size = size) +} + + # Prior-constructors ---- #' Normal Prior Distribution @@ -574,6 +604,23 @@ prior_horseshoe <- function( "shrinkage_horseshoe(prior_local_{name}, prior_global_{name}, prior_c2_{name});" ) ), + auxiliary_initial_values = \(name, size) { + local_size <- if (is.numeric(size)) size else 1 + stats::setNames( + list( + abs(local_rt(local_size, df, 0, 1)), + abs(local_rt(1, df_global, 0, scale_global)), + local_rinvgamma(1, df_slab / 2, df_slab / 2) + ), + paste0(c("prior_local_", "prior_global_", "prior_slab_"), name) + ) + }, + auxiliary_size = \(name, size) { + stats::setNames( + list(size, 1, 1), + paste0(c("prior_local_", "prior_global_", "prior_slab_"), name) + ) + }, centre = 0, sample = \(n) { local_rhorseshoe( diff --git a/R/generics.R b/R/generics.R index 3a7f26eb..4c361df0 100755 --- a/R/generics.R +++ b/R/generics.R @@ -166,6 +166,46 @@ size <- function(object) { } +# auxiliaryInitialValues ---- + +#' `auxiliaryInitialValues` +#' +#' Obtain initial values for any additional Stan parameters introduced by an object. +#' +#' @param object where to get the auxiliary initial values from. +#' @param ... additional options. +#' +#' @details +#' Some objects, such as shrinkage priors, declare additional Stan parameters +#' beyond the main model parameter. This helper returns named initial values for +#' those additional parameters. +#' +#' @keywords internal +auxiliaryInitialValues <- function(object, ...) { + UseMethod("auxiliaryInitialValues") +} + + +# auxiliarySize ---- + +#' `auxiliarySize` +#' +#' Obtain sizes for any additional Stan parameters introduced by an object. +#' +#' @param object where to get the auxiliary parameter sizes from. +#' @param ... additional options. +#' +#' @details +#' Some objects, such as shrinkage priors, declare additional Stan parameters +#' beyond the main model parameter. This helper returns named sizes for those +#' additional parameters so their initial values can be expanded consistently. +#' +#' @keywords internal +auxiliarySize <- function(object, ...) { + UseMethod("auxiliarySize") +} + + # generateQuantities ---- #' `generateQuantities` diff --git a/_pkgdown.yml b/_pkgdown.yml index 6bb901bc..eca5c685 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -136,6 +136,8 @@ reference: - ParameterList - ParameterList-Getter-Methods - Prior-Getter-Methods + - auxiliaryInitialValues + - auxiliarySize - read_stan - StanModel - StanModule diff --git a/man/Prior-Getter-Methods.Rd b/man/Prior-Getter-Methods.Rd index e2bfc0dd..f7627217 100644 --- a/man/Prior-Getter-Methods.Rd +++ b/man/Prior-Getter-Methods.Rd @@ -3,14 +3,24 @@ \name{Prior-Getter-Methods} \alias{Prior-Getter-Methods} \alias{initialValues.Prior} +\alias{auxiliaryInitialValues.Prior} +\alias{auxiliarySize.Prior} \title{Prior Getter Functions} \usage{ \method{initialValues}{Prior}(object, ...) + +\method{auxiliaryInitialValues}{Prior}(object, name, size, ...) + +\method{auxiliarySize}{Prior}(object, name, size, ...) } \arguments{ \item{object}{(\code{Prior}) a prior Distribution} \item{...}{Not Used.} + +\item{name}{(\code{character}) the name of the parameter the prior distribution is for} + +\item{size}{the parameter size.} } \description{ Getter functions for the slots of a \code{\link{Prior}} object @@ -19,6 +29,10 @@ Getter functions for the slots of a \code{\link{Prior}} object \itemize{ \item \code{initialValues(Prior)}: The prior's initial value +\item \code{auxiliaryInitialValues(Prior)}: The prior's auxiliary initial values + +\item \code{auxiliarySize(Prior)}: The prior's auxiliary parameter sizes + }} \seealso{ Other Prior-internal: diff --git a/man/Prior-class.Rd b/man/Prior-class.Rd index 95693e17..d61ac529 100644 --- a/man/Prior-class.Rd +++ b/man/Prior-class.Rd @@ -18,6 +18,8 @@ Prior( repr_parameters = "", repr_transformed_parameters = "", repr_generated_quantities = "", + auxiliary_initial_values = function(name, size) list(), + auxiliary_size = function(name, size) list(), limits = c(-Inf, Inf), .omit_zero_lower_truncation = FALSE, .allow_vectors = FALSE, @@ -46,6 +48,12 @@ the same names as the \code{paramaters} slot.} \item{repr_generated_quantities}{(\code{string}) the Stan code representation for the generated quantities block.} +\item{auxiliary_initial_values}{(\code{function}) a function that returns initial values for extra Stan parameters +introduced by the prior.} + +\item{auxiliary_size}{(\code{function}) a function that returns sizes for extra Stan parameters introduced +by the prior.} + \item{limits}{(\code{numeric}) the lower and upper limits for a truncated distribution} \item{.omit_zero_lower_truncation}{(\code{flag}) whether to omit a lower-zero truncation adjustment.} @@ -72,6 +80,10 @@ Specifies the prior distribution in a Stan Model \item{\code{repr_generated_quantities}}{(\code{string})\cr See arguments.} +\item{\code{auxiliary_initial_values}}{(\code{function})\cr See arguments.} + +\item{\code{auxiliary_size}}{(\code{function})\cr See arguments.} + \item{\code{centre}}{(\code{numeric})\cr See arguments.} \item{\code{validation}}{(\code{list})\cr See arguments.} diff --git a/man/auxiliaryInitialValues.Rd b/man/auxiliaryInitialValues.Rd new file mode 100644 index 00000000..111f8292 --- /dev/null +++ b/man/auxiliaryInitialValues.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/generics.R +\name{auxiliaryInitialValues} +\alias{auxiliaryInitialValues} +\title{\code{auxiliaryInitialValues}} +\usage{ +auxiliaryInitialValues(object, ...) +} +\arguments{ +\item{object}{where to get the auxiliary initial values from.} + +\item{...}{additional options.} +} +\description{ +Obtain initial values for any additional Stan parameters introduced by an object. +} +\details{ +Some objects, such as shrinkage priors, declare additional Stan parameters +beyond the main model parameter. This helper returns named initial values for +those additional parameters. +} +\keyword{internal} diff --git a/man/auxiliarySize.Rd b/man/auxiliarySize.Rd new file mode 100644 index 00000000..8c71e54b --- /dev/null +++ b/man/auxiliarySize.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/generics.R +\name{auxiliarySize} +\alias{auxiliarySize} +\title{\code{auxiliarySize}} +\usage{ +auxiliarySize(object, ...) +} +\arguments{ +\item{object}{where to get the auxiliary parameter sizes from.} + +\item{...}{additional options.} +} +\description{ +Obtain sizes for any additional Stan parameters introduced by an object. +} +\details{ +Some objects, such as shrinkage priors, declare additional Stan parameters +beyond the main model parameter. This helper returns named sizes for those +additional parameters so their initial values can be expanded consistently. +} +\keyword{internal} diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 3cba84fb..f3bb1f18 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -403,6 +403,24 @@ test_that("prior_horseshoe works as expected", { x_inits <- initialValues(x) expect_numeric(x_inits, len = 1) + auxiliary_inits <- auxiliaryInitialValues(x, name = "beta", size = 3) + expect_equal( + names(auxiliary_inits), + c("prior_local_beta", "prior_global_beta", "prior_slab_beta") + ) + expect_numeric(auxiliary_inits$prior_local_beta, len = 3, lower = 0) + expect_numeric(auxiliary_inits$prior_global_beta, len = 1, lower = 0) + expect_numeric(auxiliary_inits$prior_slab_beta, len = 1, lower = 0) + + expect_equal( + auxiliarySize(x, name = "beta", size = "p"), + list( + prior_local_beta = "p", + prior_global_beta = 1, + prior_slab_beta = 1 + ) + ) + x_stan_module <- as.StanModule(x, name = "beta", size = "p") expect_equal( x_stan_module@data, @@ -478,3 +496,16 @@ parameters { model_obj <- cmdstanr::cmdstan_model(stan_file, compile = FALSE) expect_true(model_obj$check_syntax(quiet = TRUE)) }) + +test_that("default auxiliary prior methods return empty lists", { + x <- prior_normal(0, 1) + + expect_equal( + auxiliaryInitialValues(x, name = "beta", size = 3), + list() + ) + expect_equal( + auxiliarySize(x, name = "beta", size = 3), + list() + ) +}) diff --git a/tests/testthat/test-initialValues.R b/tests/testthat/test-initialValues.R index 333f0a10..a8a7c596 100644 --- a/tests/testthat/test-initialValues.R +++ b/tests/testthat/test-initialValues.R @@ -93,6 +93,45 @@ test_that("initialValues() works also for vectorized parameters", { expect_equal(ivs[[1]], ivs[[2]]) }) +test_that("initialValues() includes auxiliary prior parameters", { + pars <- ParameterList( + Parameter( + name = "beta", + prior = prior_horseshoe( + df = 1, + df_global = 1, + df_slab = 4, + scale_global = 0.3, + scale_slab = 2 + ), + size = "p" + ) + ) + + set.seed(342) + initial_values <- initialValues(pars, n_chains = 2) + + expect_equal( + names(initial_values[[1]]), + c( + "beta", + "prior_local_beta", + "prior_global_beta", + "prior_slab_beta" + ) + ) + expect_equal(names(initial_values[[1]]), names(initial_values[[2]])) + expect_true(initial_values[[1]]$prior_global_beta > 0) + expect_true(initial_values[[1]]$prior_slab_beta > 0) + expect_false(identical(initial_values[[1]], initial_values[[2]])) + + expanded <- ensure_initial_values(initial_values, list(p = 3), pars) + expect_length(expanded[[1]]$beta, 3) + expect_length(expanded[[1]]$prior_local_beta, 3) + expect_length(expanded[[1]]$prior_global_beta, 1) + expect_length(expanded[[1]]$prior_slab_beta, 1) +}) + test_that("ensure_initial_values() works as expected", { pars <- ParameterList( Parameter(name = "p1", prior = prior_beta(4, 2), size = 1), From 792b3860350013e01e7f66a2f446a091fc7108c1 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 15:41:24 +0800 Subject: [PATCH 26/30] polish vignette --- NAMESPACE | 4 ++-- R/Prior.R | 6 ++++-- man/Prior-Getter-Methods.Rd | 2 +- man/Prior-Shared.Rd | 2 ++ man/Prior-class.Rd | 2 +- man/as.StanModule.Prior.Rd | 2 ++ man/render_stan_limits.Prior.Rd | 11 +++++++++++ man/render_stan_limits.Rd | 3 +-- man/render_stan_limits.numeric.Rd | 11 +++++++++++ vignettes/covariate-selection-horseshoe.Rmd | 4 +++- 10 files changed, 38 insertions(+), 9 deletions(-) create mode 100644 man/render_stan_limits.Prior.Rd create mode 100644 man/render_stan_limits.numeric.Rd diff --git a/NAMESPACE b/NAMESPACE index 7b08324b..f8dd44c7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -79,10 +79,10 @@ S3method(as_stan_list,QuantityGeneratorPopulation) S3method(as_stan_list,QuantityGeneratorPrediction) S3method(as_stan_list,QuantityGeneratorSubject) S3method(as_stan_list,default) -S3method(auxiliaryInitialValues,Prior) -S3method(auxiliarySize,Prior) S3method(autoplot,LongitudinalQuantities) S3method(autoplot,SurvivalQuantities) +S3method(auxiliaryInitialValues,Prior) +S3method(auxiliarySize,Prior) S3method(brierScore,SurvivalQuantities) S3method(coalesceGridTime,GridFixed) S3method(coalesceGridTime,GridGrouped) diff --git a/R/Prior.R b/R/Prior.R index eb718bfa..b98823c7 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -261,6 +261,8 @@ render_stan_limits <- function(object) { UseMethod("render_stan_limits") } +#' describeIn render_stan_limits method for [Prior] objects. +#' @exportS3Method NULL render_stan_limits.Prior <- function(object) { if ( object@.omit_zero_lower_truncation && @@ -271,6 +273,8 @@ render_stan_limits.Prior <- function(object) { render_stan_limits(object@limits) } +#' describeIn render_stan_limits method for numeric vectors. +#' @exportS3Method NULL render_stan_limits.numeric <- function(object) { limits <- object l_bound <- if (limits[[1]] > -Inf) limits[[1]] else "" @@ -286,7 +290,6 @@ render_stan_limits.numeric <- function(object) { return(string) } - #' @rdname show-object #' @export setMethod( @@ -299,7 +302,6 @@ setMethod( } ) - #' `Prior` -> `StanModule` #' #' Converts a [`Prior`] object to a [`StanModule`] object diff --git a/man/Prior-Getter-Methods.Rd b/man/Prior-Getter-Methods.Rd index f7627217..f74032d7 100644 --- a/man/Prior-Getter-Methods.Rd +++ b/man/Prior-Getter-Methods.Rd @@ -20,7 +20,7 @@ \item{name}{(\code{character}) the name of the parameter the prior distribution is for} -\item{size}{the parameter size.} +\item{size}{(\code{numeric} or \code{character}) the parameter size.} } \description{ Getter functions for the slots of a \code{\link{Prior}} object diff --git a/man/Prior-Shared.Rd b/man/Prior-Shared.Rd index d578d7b7..8b1cdcfa 100644 --- a/man/Prior-Shared.Rd +++ b/man/Prior-Shared.Rd @@ -13,6 +13,8 @@ \item{name}{(\code{character}) the name of the parameter the prior distribution is for} +\item{size}{(\code{numeric} or \code{character}) the parameter size.} + \item{...}{Not Used.} } \description{ diff --git a/man/Prior-class.Rd b/man/Prior-class.Rd index d61ac529..432900f4 100644 --- a/man/Prior-class.Rd +++ b/man/Prior-class.Rd @@ -54,7 +54,7 @@ introduced by the prior.} \item{auxiliary_size}{(\code{function}) a function that returns sizes for extra Stan parameters introduced by the prior.} -\item{limits}{(\code{numeric}) the lower and upper limits for a truncated distribution} +\item{limits}{(\code{numeric}) the lower and upper limits for a truncated distribution.} \item{.omit_zero_lower_truncation}{(\code{flag}) whether to omit a lower-zero truncation adjustment.} diff --git a/man/as.StanModule.Prior.Rd b/man/as.StanModule.Prior.Rd index ad90e95f..0f2b8410 100644 --- a/man/as.StanModule.Prior.Rd +++ b/man/as.StanModule.Prior.Rd @@ -11,6 +11,8 @@ \item{name}{(\code{character}) the name of the parameter the prior distribution is for} +\item{size}{(\code{numeric} or \code{character}) the parameter size.} + \item{...}{Not Used.} } \description{ diff --git a/man/render_stan_limits.Prior.Rd b/man/render_stan_limits.Prior.Rd new file mode 100644 index 00000000..57e14e2b --- /dev/null +++ b/man/render_stan_limits.Prior.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Prior.R +\name{render_stan_limits.Prior} +\alias{render_stan_limits.Prior} +\title{describeIn render_stan_limits method for \link{Prior} objects.} +\usage{ +\method{render_stan_limits}{Prior}(object) +} +\description{ +describeIn render_stan_limits method for \link{Prior} objects. +} diff --git a/man/render_stan_limits.Rd b/man/render_stan_limits.Rd index 0ca4d492..c4e87e25 100644 --- a/man/render_stan_limits.Rd +++ b/man/render_stan_limits.Rd @@ -7,8 +7,7 @@ render_stan_limits(object) } \arguments{ -\item{object}{a \code{\link{Prior}} or numeric lower/upper limits for a -truncated distribution} +\item{object}{(\code{Prior | numeric}) prior or lower and upper limits for a truncated distribution.} } \value{ (\code{character}) the Stan syntax for truncated distributions diff --git a/man/render_stan_limits.numeric.Rd b/man/render_stan_limits.numeric.Rd new file mode 100644 index 00000000..8807f3f2 --- /dev/null +++ b/man/render_stan_limits.numeric.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Prior.R +\name{render_stan_limits.numeric} +\alias{render_stan_limits.numeric} +\title{describeIn render_stan_limits method for numeric vectors.} +\usage{ +\method{render_stan_limits}{numeric}(object) +} +\description{ +describeIn render_stan_limits method for numeric vectors. +} diff --git a/vignettes/covariate-selection-horseshoe.Rmd b/vignettes/covariate-selection-horseshoe.Rmd index e93d8274..bd4db096 100644 --- a/vignettes/covariate-selection-horseshoe.Rmd +++ b/vignettes/covariate-selection-horseshoe.Rmd @@ -190,7 +190,9 @@ fit <- sampleStanModel( chains = 4, parallel_chains = 4, seed = 325, - refresh = 0 + refresh = 0, + show_exceptions = FALSE, + show_messages = FALSE ) ``` From b3aa17d5eda527d18a33b00ff5204655e6aceeea Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 20:21:27 +0800 Subject: [PATCH 27/30] fix linters --- tests/testthat/test-JointModel.R | 8 ++++---- tests/testthat/test-Prior.R | 2 +- tests/testthat/test-brms_compare.R | 28 ++++++++-------------------- 3 files changed, 13 insertions(+), 25 deletions(-) diff --git a/tests/testthat/test-JointModel.R b/tests/testthat/test-JointModel.R index e1a77dfa..637ddb1b 100644 --- a/tests/testthat/test-JointModel.R +++ b/tests/testthat/test-JointModel.R @@ -22,7 +22,7 @@ normalise_stan_declarations <- function(x) { # Snapshot only the assembled parameter declarations and constant-parameter # declarations. This catches regressions in dynamic parameter insertion without # snapshotting the full generated Stan program. -snapshot_joint_model_parameter_declarations <- function(object) { +snapshot_joint_model_parameter_decl <- function(object) { stan <- as.StanModule(object) cat("parameters\n") cat(normalise_stan_declarations(stan@parameters), sep = "\n") @@ -47,7 +47,7 @@ snapshot_joint_model_parameter_declarations <- function(object) { test_that("JointModel snapshots assembled parameter declarations", { expect_snapshot( - snapshot_joint_model_parameter_declarations(JointModel( + snapshot_joint_model_parameter_decl(JointModel( longitudinal = LongitudinalRandomSlope(), survival = SurvivalWeibullPH(), link = linkDSLD() @@ -55,7 +55,7 @@ test_that("JointModel snapshots assembled parameter declarations", { ) expect_snapshot( - snapshot_joint_model_parameter_declarations(JointModel( + snapshot_joint_model_parameter_decl(JointModel( longitudinal = LongitudinalGSF(centred = FALSE), survival = SurvivalWeibullPH(), link = Link(linkTTG(), linkDSLD(), linkGrowth()) @@ -63,7 +63,7 @@ test_that("JointModel snapshots assembled parameter declarations", { ) expect_snapshot( - snapshot_joint_model_parameter_declarations(JointModel( + snapshot_joint_model_parameter_decl(JointModel( survival = SurvivalExponential(lambda = prior_const(0.5)) )) ) diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index f3bb1f18..1de8755a 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -481,7 +481,7 @@ test_that("prior_horseshoe works as expected", { header <- merge( StanModule("base/functions.stan"), StanModule( - "data { + "data { int p; } parameters { diff --git a/tests/testthat/test-brms_compare.R b/tests/testthat/test-brms_compare.R index ed235af0..86376287 100644 --- a/tests/testthat/test-brms_compare.R +++ b/tests/testthat/test-brms_compare.R @@ -396,27 +396,15 @@ test_that("jmpost and brms get similar horseshoe estimates for survival models", posterior::as_draws_df() b_c2 <- 2^2 * b_shrink_pars$hs_slab + b_g2 <- b_shrink_pars$hs_global^2 + shrink_fun <- function(hs_local) { + 1 / (1 + b_g2 * hs_local^2 / b_c2) + } b_shrinkage <- dplyr::tibble( - cov1 = 1 / - (1 + - b_shrink_pars$hs_global^2 * - b_shrink_pars$`hs_local[1]`^2 / - b_c2), - cov2 = 1 / - (1 + - b_shrink_pars$hs_global^2 * - b_shrink_pars$`hs_local[2]`^2 / - b_c2), - cov3 = 1 / - (1 + - b_shrink_pars$hs_global^2 * - b_shrink_pars$`hs_local[3]`^2 / - b_c2), - cov4 = 1 / - (1 + - b_shrink_pars$hs_global^2 * - b_shrink_pars$`hs_local[4]`^2 / - b_c2) + cov1 = shrink_fun(b_shrink_pars$`hs_local[1]`), + cov2 = shrink_fun(b_shrink_pars$`hs_local[2]`), + cov3 = shrink_fun(b_shrink_pars$`hs_local[3]`), + cov4 = shrink_fun(b_shrink_pars$`hs_local[4]`) ) j_shrinkage <- shrinkage(mp) |> From 1a3d9a7153a26cede514b029330375bc31236214 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 20:23:11 +0800 Subject: [PATCH 28/30] fix docs and spelling --- R/Prior.R | 4 ++-- inst/WORDLIST | 1 + man/render_stan_limits.Prior.Rd | 11 ----------- man/render_stan_limits.Rd | 13 +++++++++++++ man/render_stan_limits.numeric.Rd | 11 ----------- 5 files changed, 16 insertions(+), 24 deletions(-) delete mode 100644 man/render_stan_limits.Prior.Rd delete mode 100644 man/render_stan_limits.numeric.Rd diff --git a/R/Prior.R b/R/Prior.R index b98823c7..c0346fc2 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -261,7 +261,7 @@ render_stan_limits <- function(object) { UseMethod("render_stan_limits") } -#' describeIn render_stan_limits method for [Prior] objects. +#' @describeIn render_stan_limits method for [Prior] objects. #' @exportS3Method NULL render_stan_limits.Prior <- function(object) { if ( @@ -273,7 +273,7 @@ render_stan_limits.Prior <- function(object) { render_stan_limits(object@limits) } -#' describeIn render_stan_limits method for numeric vectors. +#' @describeIn render_stan_limits method for numeric vectors. #' @exportS3Method NULL render_stan_limits.numeric <- function(object) { limits <- object diff --git a/inst/WORDLIST b/inst/WORDLIST index 9c20732b..4011823f 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -61,6 +61,7 @@ bl brierScore cdot coeficient +covariate's csv dSLD docetaxel diff --git a/man/render_stan_limits.Prior.Rd b/man/render_stan_limits.Prior.Rd deleted file mode 100644 index 57e14e2b..00000000 --- a/man/render_stan_limits.Prior.Rd +++ /dev/null @@ -1,11 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Prior.R -\name{render_stan_limits.Prior} -\alias{render_stan_limits.Prior} -\title{describeIn render_stan_limits method for \link{Prior} objects.} -\usage{ -\method{render_stan_limits}{Prior}(object) -} -\description{ -describeIn render_stan_limits method for \link{Prior} objects. -} diff --git a/man/render_stan_limits.Rd b/man/render_stan_limits.Rd index c4e87e25..8b58f62d 100644 --- a/man/render_stan_limits.Rd +++ b/man/render_stan_limits.Rd @@ -2,9 +2,15 @@ % Please edit documentation in R/Prior.R \name{render_stan_limits} \alias{render_stan_limits} +\alias{render_stan_limits.Prior} +\alias{render_stan_limits.numeric} \title{Creates Stan Syntax for Truncated distributions} \usage{ render_stan_limits(object) + +\method{render_stan_limits}{Prior}(object) + +\method{render_stan_limits}{numeric}(object) } \arguments{ \item{object}{(\code{Prior | numeric}) prior or lower and upper limits for a truncated distribution.} @@ -15,4 +21,11 @@ render_stan_limits(object) \description{ This function creates the Stan syntax for truncated distributions } +\section{Methods (by class)}{ +\itemize{ +\item \code{render_stan_limits(Prior)}: method for \link{Prior} objects. + +\item \code{render_stan_limits(numeric)}: method for numeric vectors. + +}} \keyword{internal} diff --git a/man/render_stan_limits.numeric.Rd b/man/render_stan_limits.numeric.Rd deleted file mode 100644 index 8807f3f2..00000000 --- a/man/render_stan_limits.numeric.Rd +++ /dev/null @@ -1,11 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Prior.R -\name{render_stan_limits.numeric} -\alias{render_stan_limits.numeric} -\title{describeIn render_stan_limits method for numeric vectors.} -\usage{ -\method{render_stan_limits}{numeric}(object) -} -\description{ -describeIn render_stan_limits method for numeric vectors. -} From 6bcbbc100310c53a905af6724ab0a8281478f936 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 9 Jul 2026 20:26:37 +0800 Subject: [PATCH 29/30] fix docs --- R/DataSurvival.R | 1 + R/JointModelSamples.R | 1 + R/Parameter.R | 3 ++- R/generics.R | 2 ++ man/as.StanModule.ParameterDeclaration.Rd | 4 +++- man/covariates.DataSurvival.Rd | 2 ++ man/covariates.Rd | 2 ++ man/shrinkage.JointModelSamples.Rd | 2 ++ man/shrinkage.Rd | 2 ++ 9 files changed, 17 insertions(+), 2 deletions(-) diff --git a/R/DataSurvival.R b/R/DataSurvival.R index 0640821c..51a8ffda 100644 --- a/R/DataSurvival.R +++ b/R/DataSurvival.R @@ -139,6 +139,7 @@ model.matrix.DataSurvival <- function( #' #' @typed object: DataSurvival #' Survival data object to extract covariate names from. +#' @param ... not used. #' #' @family covariates #' @export diff --git a/R/JointModelSamples.R b/R/JointModelSamples.R index 131990b2..f67b9b2c 100644 --- a/R/JointModelSamples.R +++ b/R/JointModelSamples.R @@ -188,6 +188,7 @@ saveObject.JointModelSamples <- function(object, file, ...) { #' #' @typed object: JointModelSamples #' the object to extract shrinkage factors from. +#' @param ... not used. #' @return the shrinkage factors correctly named after the covariate names. #' #' @family shrinkage diff --git a/R/Parameter.R b/R/Parameter.R index 9618191d..6d1ee22f 100755 --- a/R/Parameter.R +++ b/R/Parameter.R @@ -169,11 +169,12 @@ render_stan_const_declaration <- function(name, size, limits) { #' `transformed parameters` block. #' #' @inheritParams Parameter-Shared +#' @param ... not used. #' #' @return A [`StanModule`] object. #' #' @export -as.StanModule.ParameterDeclaration <- function(object) { +as.StanModule.ParameterDeclaration <- function(object, ...) { declaration <- if (object@prior@.is_const) { render_stan_const_declaration( name = object@name, diff --git a/R/generics.R b/R/generics.R index 4c361df0..b4548527 100755 --- a/R/generics.R +++ b/R/generics.R @@ -572,6 +572,7 @@ saveObject <- function(object, file, ...) { #' #' @typed object: ANY #' the object to extract covariate names from. +#' @param ... additional arguments added by methods. #' #' @family covariates #' @export @@ -583,6 +584,7 @@ covariates <- function(object, ...) { #' #' @typed object: ANY #' the object to extract shrinkage factors from. +#' @param ... additional arguments added by methods. #' #' @family shrinkage #' @export diff --git a/man/as.StanModule.ParameterDeclaration.Rd b/man/as.StanModule.ParameterDeclaration.Rd index f603830c..7ec036dd 100644 --- a/man/as.StanModule.ParameterDeclaration.Rd +++ b/man/as.StanModule.ParameterDeclaration.Rd @@ -4,10 +4,12 @@ \alias{as.StanModule.ParameterDeclaration} \title{\code{Parameter} Declaration -> \code{StanModule}} \usage{ -\method{as.StanModule}{ParameterDeclaration}(object) +\method{as.StanModule}{ParameterDeclaration}(object, ...) } \arguments{ \item{object}{(\code{Parameter}) a prior Distribution} + +\item{...}{not used.} } \value{ A \code{\link{StanModule}} object. diff --git a/man/covariates.DataSurvival.Rd b/man/covariates.DataSurvival.Rd index 0cea798d..935099b3 100644 --- a/man/covariates.DataSurvival.Rd +++ b/man/covariates.DataSurvival.Rd @@ -8,6 +8,8 @@ } \arguments{ \item{object}{(\code{DataSurvival}) Survival data object to extract covariate names from.} + +\item{...}{not used.} } \description{ Extract Covariate Names for Survival Data diff --git a/man/covariates.Rd b/man/covariates.Rd index 5837075a..ef3e5cc8 100644 --- a/man/covariates.Rd +++ b/man/covariates.Rd @@ -8,6 +8,8 @@ covariates(object, ...) } \arguments{ \item{object}{(\code{ANY}) the object to extract covariate names from.} + +\item{...}{additional arguments added by methods.} } \description{ Extract Covariate Names diff --git a/man/shrinkage.JointModelSamples.Rd b/man/shrinkage.JointModelSamples.Rd index 72c92543..ab9d3d48 100644 --- a/man/shrinkage.JointModelSamples.Rd +++ b/man/shrinkage.JointModelSamples.Rd @@ -8,6 +8,8 @@ } \arguments{ \item{object}{(\code{JointModelSamples}) the object to extract shrinkage factors from.} + +\item{...}{not used.} } \value{ the shrinkage factors correctly named after the covariate names. diff --git a/man/shrinkage.Rd b/man/shrinkage.Rd index 9a6d9b07..cbb65177 100644 --- a/man/shrinkage.Rd +++ b/man/shrinkage.Rd @@ -8,6 +8,8 @@ shrinkage(object, ...) } \arguments{ \item{object}{(\code{ANY}) the object to extract shrinkage factors from.} + +\item{...}{additional arguments added by methods.} } \description{ Extract Covariate Shrinkage Factors from the Posterior Samples From 9576e1d8aa072f9f9bfbb540b3b652ee74c3bcaf Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 10 Jul 2026 14:08:22 +0800 Subject: [PATCH 30/30] fix snaps --- tests/testthat/_snaps/JointModel.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/testthat/_snaps/JointModel.md b/tests/testthat/_snaps/JointModel.md index df22db35..48de8683 100644 --- a/tests/testthat/_snaps/JointModel.md +++ b/tests/testthat/_snaps/JointModel.md @@ -1,7 +1,7 @@ # JointModel snapshots assembled parameter declarations Code - snapshot_joint_model_parameter_declarations(JointModel(longitudinal = LongitudinalRandomSlope(), + snapshot_joint_model_parameter_decl(JointModel(longitudinal = LongitudinalRandomSlope(), survival = SurvivalWeibullPH(), link = linkDSLD())) Output parameters @@ -21,7 +21,7 @@ --- Code - snapshot_joint_model_parameter_declarations(JointModel(longitudinal = LongitudinalGSF( + snapshot_joint_model_parameter_decl(JointModel(longitudinal = LongitudinalGSF( centred = FALSE), survival = SurvivalWeibullPH(), link = Link(linkTTG(), linkDSLD(), linkGrowth()))) Output @@ -52,7 +52,7 @@ --- Code - snapshot_joint_model_parameter_declarations(JointModel(survival = SurvivalExponential( + snapshot_joint_model_parameter_decl(JointModel(survival = SurvivalExponential( lambda = prior_const(0.5)))) Output parameters