450: Add horseshoe prior, as well as constant prior

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

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,13 +81,16 @@ 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)
 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)

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).

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

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

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
+}

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"
             )
         )

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(

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(

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"
             )
         )