diff --git a/DESCRIPTION b/DESCRIPTION index 22baeab7f..617d78679 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 34924da3c..f8dd44c7b 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) @@ -80,6 +81,8 @@ S3method(as_stan_list,QuantityGeneratorSubject) S3method(as_stan_list,default) S3method(autoplot,LongitudinalQuantities) S3method(autoplot,SurvivalQuantities) +S3method(auxiliaryInitialValues,Prior) +S3method(auxiliarySize,Prior) S3method(brierScore,SurvivalQuantities) S3method(coalesceGridTime,GridFixed) S3method(coalesceGridTime,GridGrouped) @@ -87,6 +90,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 +185,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 +253,7 @@ export(as_stan_list) export(autoplot) export(brierScore) export(compileStanModel) +export(covariates) export(cut_data) export(enableGQ) export(enableLink) @@ -266,7 +272,9 @@ export(merge) export(populationHR) export(prior_beta) export(prior_cauchy) +export(prior_const) export(prior_gamma) +export(prior_horseshoe) export(prior_init_only) export(prior_invgamma) export(prior_logistic) @@ -284,6 +292,7 @@ export(sampleSubjects) export(saveObject) export(set_limits) export(show) +export(shrinkage) export(write_stan) exportClasses(DataJoint) exportClasses(DataLongitudinal) diff --git a/NEWS.md b/NEWS.md index 74bd69161..206641b3a 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/DataSurvival.R b/R/DataSurvival.R index 3d9c3d8bb..51a8ffda8 100644 --- a/R/DataSurvival.R +++ b/R/DataSurvival.R @@ -128,11 +128,26 @@ 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 } +#' Extract Covariate Names for Survival Data +#' +#' @typed object: DataSurvival +#' Survival data object to extract covariate names from. +#' @param ... not used. +#' +#' @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/JointModel.R b/R/JointModel.R index a8ea2e153..74ccfb6d8 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/JointModelSamples.R b/R/JointModelSamples.R index 697ff21e1..f67b9b2c9 100644 --- a/R/JointModelSamples.R +++ b/R/JointModelSamples.R @@ -183,3 +183,30 @@ 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. +#' @param ... not used. +#' @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/LongitudinalClaretBruno.R b/R/LongitudinalClaretBruno.R index 99e487b5a..15e499194 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 83175ae68..46462501a 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 67f74d219..2964d91ac 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 8ae2b39a9..361c830cc 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 1db50bbaa..6d1ee22ff 100755 --- a/R/Parameter.R +++ b/R/Parameter.R @@ -86,7 +86,119 @@ setValidity( #' @family as.StanModule #' @export as.StanModule.Parameter <- function(object, ...) { - as.StanModule(object@prior, name = object@name) + 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 +#' @param ... not used. +#' +#' @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 +251,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 2e6e0b54e..f46b6dc0e 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,16 +152,27 @@ 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)) - 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) } ) } @@ -169,9 +181,18 @@ 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) - return(x) + parameters <- Filter(\(x) !x@prior@.is_const, object@parameters) + 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 b96339846..c0346fc2e 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 @@ -31,12 +33,19 @@ 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 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. #' @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. #' #' @family Prior-internal #' @export Prior @@ -48,11 +57,18 @@ NULL "display" = "character", "repr_model" = "character", "repr_data" = "character", + "repr_parameters" = "character", + "repr_transformed_parameters" = "character", + "repr_generated_quantities" = "character", + "auxiliary_initial_values" = "function", + "auxiliary_size" = "function", "centre" = "numeric", "validation" = "list", "sample" = "function", "limits" = "numeric", - ".allow_vectors" = "logical" + ".omit_zero_lower_truncation" = "logical", + ".allow_vectors" = "logical", + ".is_const" = "logical" ) ) @@ -63,6 +79,18 @@ 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 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 @@ -73,9 +101,13 @@ 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 +#' whether this prior fixes the parameter at a constant value. #' @rdname Prior-class Prior <- function( parameters, @@ -85,19 +117,33 @@ Prior <- function( centre, validation, sample, + repr_parameters = "", + repr_transformed_parameters = "", + repr_generated_quantities = "", + auxiliary_initial_values = \(name, size) list(), + auxiliary_size = \(name, size) list(), limits = c(-Inf, Inf), - .allow_vectors = FALSE + .omit_zero_lower_truncation = FALSE, + .allow_vectors = FALSE, + .is_const = FALSE ) { .Prior( parameters = parameters, repr_model = repr_model, repr_data = repr_data, + 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, sample = sample, limits = limits, - .allow_vectors = .allow_vectors + .omit_zero_lower_truncation = .omit_zero_lower_truncation, + .allow_vectors = .allow_vectors, + .is_const = .is_const ) } @@ -150,6 +196,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) } ) @@ -177,7 +237,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) if ( display_limits != "" && display_string != "" && @@ -192,12 +252,31 @@ 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") +} + +#' @describeIn render_stan_limits method for [Prior] objects. +#' @exportS3Method NULL +render_stan_limits.Prior <- function(object) { + if ( + object@.omit_zero_lower_truncation && + identical(object@limits, c(0, Inf)) + ) { + return("") + } + 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 "" u_bound <- if (limits[[2]] < Inf) limits[[2]] else "" string <- "" @@ -211,7 +290,6 @@ render_stan_limits <- function(limits) { return(string) } - #' @rdname show-object #' @export setMethod( @@ -224,7 +302,6 @@ setMethod( } ) - #' `Prior` -> `StanModule` #' #' Converts a [`Prior`] object to a [`StanModule`] object @@ -234,22 +311,49 @@ setMethod( #' @family Prior-internal #' @family as.StanModule #' @export -as.StanModule.Prior <- function(object, name, ...) { +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) { + return("") + } + paste0(" ", x, collapse = "\n") + } trunctation <- if (object@repr_model != "") { - paste0(render_stan_limits(object@limits), ";") + paste0(render_stan_limits(object), ";") } else { "" } 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)), + "}}", + "generated quantities {{", + indent_stan(object@repr_generated_quantities), "}}", sep = "\n" ) - StanModule(glue::glue(string, name = name)) + StanModule(glue::glue(string, name = name, size = size)) } @@ -265,6 +369,10 @@ as.StanModule.Prior <- function(object, name, ...) { #' @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) @@ -350,6 +458,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 @@ -415,6 +537,113 @@ 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( + # 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 " + ), + 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};", + 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});" + ) + ), + 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( + 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 #' @@ -433,6 +662,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 @@ -456,7 +706,8 @@ prior_cauchy <- function(mu, sigma) { validation = list( mu = is.numeric, sigma = \(x) x > 0 - ) + ), + .omit_zero_lower_truncation = TRUE ) } @@ -484,7 +735,8 @@ prior_gamma <- function(alpha, beta) { validation = list( alpha = \(x) x > 0, beta = \(x) x > 0 - ) + ), + .omit_zero_lower_truncation = TRUE ) } @@ -511,7 +763,8 @@ prior_lognormal <- function(mu, sigma) { validation = list( mu = is.numeric, sigma = \(x) x > 0 - ) + ), + .omit_zero_lower_truncation = TRUE ) } @@ -563,7 +816,9 @@ prior_init_only <- function(dist) { dist@sample(n) }, centre = dist@centre, - validation = list() + validation = list(), + limits = dist@limits, + .omit_zero_lower_truncation = dist@.omit_zero_lower_truncation ) } @@ -695,7 +950,8 @@ prior_loglogistic <- function(alpha, beta) { validation = list( alpha = \(x) x > 0, beta = \(x) x > 0 - ) + ), + .omit_zero_lower_truncation = TRUE ) } @@ -726,7 +982,8 @@ prior_invgamma <- function(alpha, beta) { validation = list( alpha = \(x) x > 0, beta = \(x) x > 0 - ) + ), + .omit_zero_lower_truncation = TRUE ) } @@ -827,6 +1084,18 @@ 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) { + # 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) + 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/R/StanModel.R b/R/StanModel.R index 957d7d2d8..b17696d42 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 396468a2c..e63e877c5 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 d2cc2014b..efcfef4de 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 721008a92..d7f1c7ae1 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 b47562a72..b2bfe5ccb 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 07f63451f..b4548527f 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` @@ -526,3 +566,28 @@ set_limits <- function(object, lower = -Inf, upper = Inf) { saveObject <- function(object, file, ...) { UseMethod("saveObject") } + + +#' Extract Covariate Names +#' +#' @typed object: ANY +#' the object to extract covariate names from. +#' @param ... additional arguments added by methods. +#' +#' @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. +#' @param ... additional arguments added by methods. +#' +#' @family shrinkage +#' @export +shrinkage <- function(object, ...) { + UseMethod("shrinkage") +} diff --git a/R/settings.R b/R/settings.R index 555c21a89..373f4adba 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 0c46ed614..eca5c6856 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: @@ -132,6 +136,8 @@ reference: - ParameterList - ParameterList-Getter-Methods - Prior-Getter-Methods + - auxiliaryInitialValues + - auxiliarySize - read_stan - StanModel - StanModule @@ -149,6 +155,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 +192,8 @@ reference: - hazardWindows - sampleObservations - sampleSubjects + - covariates + - covariates.DataSurvival - enableLink - enableGQ - as_formula @@ -204,21 +213,18 @@ reference: contents: - STAN_BLOCKS - - title: Data Sets contents: - os_data - tumor_data - - articles: - title: Articles navbar: ~ contents: - quickstart - '`statistical-specification`' + - '`covariate-selection-horseshoe`' - '`custom-model`' - '`extending-jmpost`' - '`simulating-data`' - diff --git a/design/covariate-selection.qmd b/design/covariate-selection.qmd index eacad4c7c..ea6c45876 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: @@ -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 1c300f551..f8c6a48da 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/inst/WORDLIST b/inst/WORDLIST index 599eeebdd..4011823f1 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -21,6 +21,7 @@ Grajek HS Hoffmann IPCW +Inv JointModel JointModelSamples Kerioui @@ -60,9 +61,11 @@ bl brierScore cdot coeficient +covariate's csv dSLD docetaxel +dotsc du etc extractVariableNames @@ -71,6 +74,7 @@ geq gk groupwise gt +hardcode hardcoded ie ij @@ -80,16 +84,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 a446cfdfd..34da7415e 100755 --- a/inst/stan/base/base.stan +++ b/inst/stan/base/base.stan @@ -38,11 +38,13 @@ parameters{ {{ longitudinal.parameters }} {{ link.parameters }} {{ survival.parameters }} +{{ priors.parameters }} } transformed parameters{ +{{ priors.transformed_parameters }} {{ longitudinal.transformed_parameters }} {{ link.transformed_parameters }} {{ survival.transformed_parameters }} @@ -69,6 +71,5 @@ generated quantities{ {{ longitudinal.generated_quantities }} {{ link.generated_quantities }} {{ survival.generated_quantities }} +{{ priors.generated_quantities }} } - - diff --git a/inst/stan/base/functions.stan b/inst/stan/base/functions.stan index 39528bc01..d3f914110 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; + } } diff --git a/inst/stan/base/link.stan b/inst/stan/base/link.stan index 967bf4179..73bde668d 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 f74f5f4a5..dd4eeb4fe 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 ce371a160..01e2f150a 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 f05e7045c..a7179947c 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 11b79b6b0..1da9cf3bc 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 889c8760a..a0bf5d82a 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 e551674d0..d5e0e650a 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 5d9e6eb12..f28b43e17 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 e95f9a0e6..5b89683b8 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 629a730c8..96bc2e56d 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 b2b845b2a..d42ecb3a5 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/Local_Sample.Rd b/man/Local_Sample.Rd index 6bbe1459c..c3371fc09 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/LongitudinalModel-class.Rd b/man/LongitudinalModel-class.Rd index 35090db81..8b8bab4b2 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 2f74bdcd3..8ab68009d 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-Getter-Methods.Rd b/man/Prior-Getter-Methods.Rd index e2bfc0ddd..f74032d7e 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}{(\code{numeric} or \code{character}) 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-Shared.Rd b/man/Prior-Shared.Rd index d578d7b7e..8b1cdcfab 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 490da7d57..432900f4e 100644 --- a/man/Prior-class.Rd +++ b/man/Prior-class.Rd @@ -15,8 +15,15 @@ Prior( centre, validation, sample, + 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), - .allow_vectors = FALSE + .omit_zero_lower_truncation = FALSE, + .allow_vectors = FALSE, + .is_const = FALSE ) } \arguments{ @@ -35,9 +42,25 @@ the same names as the \code{paramaters} slot.} \item{sample}{(\code{function}) a function to sample from the prior distribution.} -\item{limits}{(\code{numeric}) the lower and upper limits for a truncated distribution} +\item{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{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.} \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 @@ -51,6 +74,16 @@ 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{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.} @@ -61,7 +94,11 @@ 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.} }} \seealso{ diff --git a/man/StanModel-class.Rd b/man/StanModel-class.Rd index 6c29171d0..91cf9f73c 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 745d5c0bf..2d217c0d6 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 000000000..7ec036dd5 --- /dev/null +++ b/man/as.StanModule.ParameterDeclaration.Rd @@ -0,0 +1,21 @@ +% 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} + +\item{...}{not used.} +} +\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/as.StanModule.Prior.Rd b/man/as.StanModule.Prior.Rd index e018daf18..0f2b84108 100644 --- a/man/as.StanModule.Prior.Rd +++ b/man/as.StanModule.Prior.Rd @@ -4,13 +4,15 @@ \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} \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/auxiliaryInitialValues.Rd b/man/auxiliaryInitialValues.Rd new file mode 100644 index 000000000..111f82920 --- /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 000000000..8c71e54ba --- /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/man/covariates.DataSurvival.Rd b/man/covariates.DataSurvival.Rd new file mode 100644 index 000000000..935099b39 --- /dev/null +++ b/man/covariates.DataSurvival.Rd @@ -0,0 +1,21 @@ +% 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.} + +\item{...}{not used.} +} +\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 000000000..ef3e5cc8d --- /dev/null +++ b/man/covariates.Rd @@ -0,0 +1,21 @@ +% 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{object}{(\code{ANY}) the object to extract covariate names from.} + +\item{...}{additional arguments added by methods.} +} +\description{ +Extract Covariate Names +} +\seealso{ +Other covariates: +\code{\link[=covariates.DataSurvival]{covariates.DataSurvival()}} +} +\concept{covariates} diff --git a/man/jmpost-settings.Rd b/man/jmpost-settings.Rd index 5557ab6d3..0dbb55278 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/prior_beta.Rd b/man/prior_beta.Rd index 1479a2e41..d53b62a55 100644 --- a/man/prior_beta.Rd +++ b/man/prior_beta.Rd @@ -17,7 +17,9 @@ 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()}}, \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 c559b096c..4204820ed 100644 --- a/man/prior_cauchy.Rd +++ b/man/prior_cauchy.Rd @@ -17,7 +17,9 @@ 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()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_const.Rd b/man/prior_const.Rd new file mode 100644 index 000000000..7012a30bc --- /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 cc386c882..78bee628c 100644 --- a/man/prior_gamma.Rd +++ b/man/prior_gamma.Rd @@ -18,6 +18,8 @@ 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()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/prior_horseshoe.Rd b/man/prior_horseshoe.Rd new file mode 100644 index 000000000..046c3f210 --- /dev/null +++ b/man/prior_horseshoe.Rd @@ -0,0 +1,48 @@ +% 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_const]{prior_const()}}, +\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/man/prior_init_only.Rd b/man/prior_init_only.Rd index b36e9153d..68e3c8c76 100644 --- a/man/prior_init_only.Rd +++ b/man/prior_init_only.Rd @@ -18,7 +18,9 @@ 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()}}, \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 c3a7f3d04..89547e7ad 100644 --- a/man/prior_invgamma.Rd +++ b/man/prior_invgamma.Rd @@ -18,7 +18,9 @@ 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()}}, \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 fc0447c14..0b05f70e9 100644 --- a/man/prior_logistic.Rd +++ b/man/prior_logistic.Rd @@ -18,7 +18,9 @@ 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()}}, \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 62a46e5d5..487395101 100644 --- a/man/prior_loglogistic.Rd +++ b/man/prior_loglogistic.Rd @@ -18,7 +18,9 @@ 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()}}, \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 863357980..db21e5009 100644 --- a/man/prior_lognormal.Rd +++ b/man/prior_lognormal.Rd @@ -18,7 +18,9 @@ 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()}}, \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 2f532e8c2..ec1abaac5 100644 --- a/man/prior_normal.Rd +++ b/man/prior_normal.Rd @@ -18,7 +18,9 @@ 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()}}, \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 5d7fbb7fa..e522cade4 100644 --- a/man/prior_normal_vector.Rd +++ b/man/prior_normal_vector.Rd @@ -18,7 +18,9 @@ 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()}}, \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 46f67d4f2..b41985dc7 100644 --- a/man/prior_std_normal.Rd +++ b/man/prior_std_normal.Rd @@ -13,7 +13,9 @@ 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()}}, \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 1101043e8..f1bbbc099 100644 --- a/man/prior_student_t.Rd +++ b/man/prior_student_t.Rd @@ -20,7 +20,9 @@ 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()}}, \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 15efa1719..2b74bb1b7 100644 --- a/man/prior_uniform.Rd +++ b/man/prior_uniform.Rd @@ -18,7 +18,9 @@ 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()}}, \code{\link[=prior_invgamma]{prior_invgamma()}}, \code{\link[=prior_logistic]{prior_logistic()}}, diff --git a/man/render_stan_const_declaration.Rd b/man/render_stan_const_declaration.Rd new file mode 100644 index 000000000..a8e38dbd8 --- /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_limits.Rd b/man/render_stan_limits.Rd index 3428e7c47..8b58f62de 100644 --- a/man/render_stan_limits.Rd +++ b/man/render_stan_limits.Rd @@ -2,12 +2,18 @@ % 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(limits) +render_stan_limits(object) + +\method{render_stan_limits}{Prior}(object) + +\method{render_stan_limits}{numeric}(object) } \arguments{ -\item{limits}{(\code{numeric}) the lower and 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 @@ -15,4 +21,11 @@ render_stan_limits(limits) \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_parameter_declaration.Rd b/man/render_stan_parameter_declaration.Rd new file mode 100644 index 000000000..1b38fc6eb --- /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 000000000..7d7122f07 --- /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/shrinkage.JointModelSamples.Rd b/man/shrinkage.JointModelSamples.Rd new file mode 100644 index 000000000..ab9d3d487 --- /dev/null +++ b/man/shrinkage.JointModelSamples.Rd @@ -0,0 +1,24 @@ +% 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.} + +\item{...}{not used.} +} +\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 000000000..cbb65177d --- /dev/null +++ b/man/shrinkage.Rd @@ -0,0 +1,21 @@ +% 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.} + +\item{...}{additional arguments added by methods.} +} +\description{ +Extract Covariate Shrinkage Factors from the Posterior Samples +} +\seealso{ +Other shrinkage: +\code{\link[=shrinkage.JointModelSamples]{shrinkage.JointModelSamples()}} +} +\concept{shrinkage} diff --git a/man/stanmodel_arguments.Rd b/man/stanmodel_arguments.Rd index 80157cf4a..8cd038b3f 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 9a9555fea..48de8683c 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_decl(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_decl(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_decl(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 @@ -10,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: @@ -38,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: @@ -66,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: @@ -94,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() @@ -119,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 1942da062..296019116 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 a531a6b29..51cb48cdf 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 cdbb06b8c..93900a768 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 715db37ad..85eccd07c 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 @@ -9,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() @@ -21,7 +20,6 @@ --- Code - x <- LongitudinalSteinFojo(sigma = prior_normal(0, 1), mu_kg = prior_gamma(2, 1)) print(x) Output @@ -29,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/Prior.md b/tests/testthat/_snaps/Prior.md index 908f17bb5..c453b0e59 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 844e993b4..4d5b7d0d7 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 @@ -22,10 +20,20 @@ beta_os_cov ~ gamma(alpha = 3, beta = 4) +--- + + Code + 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 - x <- SurvivalExponential(beta = prior_normal(0, 1)) print(x) Output @@ -37,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/_snaps/SurvivalGamma.md b/tests/testthat/_snaps/SurvivalGamma.md index ca86ce0ba..2cb719671 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 97e832c3e..326b229e0 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/helper-example_data.R b/tests/testthat/helper-example_data.R index 5f9d826cb..7f372d532 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/models/ParameterList.stan b/tests/testthat/models/ParameterList.stan index d017b72ce..68316a28f 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-DataSurvival.R b/tests/testthat/test-DataSurvival.R index 807400cb2..cbfbd3e9a 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 diff --git a/tests/testthat/test-GridPrediction.R b/tests/testthat/test-GridPrediction.R index f955d7b15..47637abe9 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-JointModel.R b/tests/testthat/test-JointModel.R index 6c447a660..637ddb1bd 100644 --- a/tests/testthat/test-JointModel.R +++ b/tests/testthat/test-JointModel.R @@ -11,6 +11,65 @@ 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_decl <- 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", { + expect_snapshot( + snapshot_joint_model_parameter_decl(JointModel( + longitudinal = LongitudinalRandomSlope(), + survival = SurvivalWeibullPH(), + link = linkDSLD() + )) + ) + + expect_snapshot( + snapshot_joint_model_parameter_decl(JointModel( + longitudinal = LongitudinalGSF(centred = FALSE), + survival = SurvivalWeibullPH(), + link = Link(linkTTG(), linkDSLD(), linkGrowth()) + )) + ) + + expect_snapshot( + snapshot_joint_model_parameter_decl(JointModel( + survival = SurvivalExponential(lambda = prior_const(0.5)) + )) + ) +}) + + test_that("JointModel print method works as expected", { expect_snapshot({ x <- JointModel( diff --git a/tests/testthat/test-Link.R b/tests/testthat/test-Link.R index 2f8f31f12..3d8dcbb2a 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-LongitudinalGSF.R b/tests/testthat/test-LongitudinalGSF.R index c7850d96e..7377a89d8 100644 --- a/tests/testthat/test-LongitudinalGSF.R +++ b/tests/testthat/test-LongitudinalGSF.R @@ -7,18 +7,14 @@ test_that("LongitudinalGSF works as expected with default arguments", { test_that("Print method for LongitudinalGSF works as expected", { - 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 d93020802..05f1f26eb 100644 --- a/tests/testthat/test-LongitudinalSteinFojo.R +++ b/tests/testthat/test-LongitudinalSteinFojo.R @@ -5,18 +5,14 @@ test_that("LongitudinalSteinFojo works as expected with default arguments", { test_that("Print method for LongitudinalSteinFojo works as expected", { - expect_snapshot({ - x <- LongitudinalSteinFojo() - print(x) - }) + x <- LongitudinalSteinFojo() + expect_snapshot(print(x)) - expect_snapshot({ - x <- LongitudinalSteinFojo( - sigma = prior_normal(0, 1), - mu_kg = prior_gamma(2, 1) - ) - 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-Parameter.R b/tests/testthat/test-Parameter.R index 58d43112b..1e3828076 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 28a4d6e84..1de8755aa 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -138,11 +138,47 @@ 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.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) + ) + 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))) @@ -210,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 @@ -224,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) @@ -269,6 +343,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 +382,130 @@ 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) + + 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, + 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( + "vector\\[p\\] prior_shrinkage_factors_beta", + x_stan_module@generated_quantities + )) + ) + expect_true( + any(grepl( + "shrinkage_horseshoe(prior_local_beta, prior_global_beta, prior_c2_beta)", + x_stan_module@generated_quantities, + fixed = TRUE + )) + ) + 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 + ) + ) + + # Check the model syntax. + 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)), + dir = tempdir() + ) + 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-SurvivalExponential.R b/tests/testthat/test-SurvivalExponential.R index 5592cb45b..a5435754a 100644 --- a/tests/testthat/test-SurvivalExponential.R +++ b/tests/testthat/test-SurvivalExponential.R @@ -16,6 +16,58 @@ 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() 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()) @@ -84,31 +136,26 @@ test_that("SurvivalExponential can recover true parameter (including covariates) test_that("Print method for SurvivalExponential works as expected", { - 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)) + + x <- SurvivalExponential(lambda = prior_const(1)) + expect_snapshot(print(x)) }) test_that("Different priors for the beta components are possible", { # Same iid prior for all beta components: - expect_snapshot({ - x <- SurvivalExponential(beta = prior_normal(0, 1)) - print(x) - }) + 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()) diff --git a/tests/testthat/test-brms_compare.R b/tests/testthat/test-brms_compare.R index 8d1ec9c80..863762877 100644 --- a/tests/testthat/test-brms_compare.R +++ b/tests/testthat/test-brms_compare.R @@ -266,3 +266,156 @@ 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", + save_pars = brms::save_pars(manual = "hs_local") + ) + + 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( + # Let's use an approximate point mass + # at lambda = 1 to get the same model as brms. + lambda = prior_const(1), + 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() |> + 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") |> + posterior::as_draws_df() + + expect_equal( + summary_quantiles(j_beta), + summary_quantiles(b_beta), + tolerance = 0.01, + ignore_attr = TRUE + ) + + # 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_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 = 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) |> + posterior::as_draws_df() |> + dplyr::as_tibble() |> + dplyr::select(ones, cov1, cov2, cov3) + + expect_equal( + summary_quantiles(j_shrinkage), + summary_quantiles(b_shrinkage), + tolerance = 0.01, + ignore_attr = TRUE + ) +}) diff --git a/tests/testthat/test-initialValues.R b/tests/testthat/test-initialValues.R index 333f0a102..a8a7c5964 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), diff --git a/tests/testthat/test-jinjar_render.R b/tests/testthat/test-jinjar_render.R index 919877370..2f33481a9 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) }) diff --git a/tests/testthat/test-misc_models.R b/tests/testthat/test-misc_models.R index 8eeddc2f7..3128acd7e 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 @@ -25,3 +25,64 @@ 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 = 1000, + iter_warmup = 100, + chains = 4, + refresh = 0, + parallel_chains = 4, + 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_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 + )) + + # 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. +}) diff --git a/tests/testthat/test-simulate.R b/tests/testthat/test-simulate.R index bb87b4481..138950fb6 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]), diff --git a/vignettes/covariate-selection-horseshoe.Rmd b/vignettes/covariate-selection-horseshoe.Rmd new file mode 100644 index 000000000..bd4db0964 --- /dev/null +++ b/vignettes/covariate-selection-horseshoe.Rmd @@ -0,0 +1,324 @@ +--- +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 +$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 + +$$ +\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.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 + +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, 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. + +```{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 = 0, + show_exceptions = FALSE, + show_messages = FALSE +) +``` + +## Inspect Coefficients + +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} +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 + +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. + +```{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" + ) +``` + +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: + +```{r shrinkage-intervals, fig.width=7, fig.height=4, message=FALSE} +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. + +## References diff --git a/vignettes/custom-model.Rmd b/vignettes/custom-model.Rmd index 5046d7135..341b47921 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 7f7e66dba..6a2ea07c6 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 9003ce2cf..5b5ebf4da 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 { diff --git a/vignettes/statistical-specification.Rmd b/vignettes/statistical-specification.Rmd index 7e7ced79b..4ef879955 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 -