451: Support for different priors

.github/workflows/check.yaml

@@ -27,7 +27,8 @@ jobs:
     secrets:
       REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
     with:
-      deps-installation-method: |
+      # Use |- so the reusable workflow receives the exact string without a trailing newline.
+      deps-installation-method: |-
         staged-dependencies
       additional-env-vars: |
         CMDSTAN=/root/.cmdstan
@@ -45,7 +46,8 @@ jobs:
     secrets:
       REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
     with:
-      deps-installation-method: |
+      # Use |- so the reusable workflow receives the exact string without a trailing newline.
+      deps-installation-method: |-
         staged-dependencies
       additional-env-vars: |
         CMDSTAN=/root/.cmdstan
@@ -61,7 +63,8 @@ jobs:
     secrets:
       REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
     with:
-      deps-installation-method: |
+      # Use |- so the reusable workflow receives the exact string without a trailing newline.
+      deps-installation-method: |-
         staged-dependencies
       auto-update: true
   gitleaks:

.github/workflows/docs.yaml

@@ -42,12 +42,11 @@ jobs:
     secrets:
       REPO_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
     with:
-      deps-installation-method: |
+      # Use |- so the reusable workflow receives the exact string without a trailing newline.
+      deps-installation-method: |-
         staged-dependencies
       additional-env-vars: |
         CMDSTAN=/root/.cmdstan
         CMDSTAN_PATH=/root/.cmdstan
         CMDSTANR_NO_VER_CHECK=true
-        NOT_CRAN=TRUE
-        JMPOST_GRAPH_SNAPSHOT=TRUE
       default-landing-page: 'main' # Set to 'latest-tag' after CRAN release or when you create a tag

.gitignore

@@ -76,6 +76,7 @@ local/
 # Compiled stan model
 inst/stanmodels/stan_model
 **/models/*
+!**/models/*.stan
 
 
 

NAMESPACE

@@ -161,6 +161,7 @@ S3method(linkTTG,LongitudinalSteinFojo)
 S3method(linkTTG,PromiseLongitudinalModel)
 S3method(linkTTG,default)
 S3method(median,Prior)
+S3method(model.matrix,DataSurvival)
 S3method(names,LinkComponent)
 S3method(names,Parameter)
 S3method(names,ParameterList)
@@ -272,6 +273,7 @@ export(prior_logistic)
 export(prior_loglogistic)
 export(prior_lognormal)
 export(prior_normal)
+export(prior_normal_vector)
 export(prior_std_normal)
 export(prior_student_t)
 export(prior_uniform)

NEWS.md

@@ -1,4 +1,3 @@
-
 # jmpost (development version)
 
 - `LongitudinalModel` objects now have a new slot `scaled_variance` which stores the used variance option, and this is part of the `print` output now thereby transparently communicating to the user the choice of the multiplicative or additive error model.
@@ -7,11 +6,13 @@
 - Included new `LongitudinalRandomEffects()` function which can be used to extract the patient-level random effects parameter samples from a `JointModelSample` object (#423).
 - Introduced the `saveObject()` method for `JointModelSample` objects in order to serialise them to disk (#431).
 - Added support for truncated prior distributions e.g. you can now apply a normal prior to a strictly positive parameter and jmpost will take care of adjusting the density accordingly (#429).
+- Added `prior_normal_vector()` which can be used to assign different normal prior distributions to the coefficients of the survival model covariates.
 - Included new Gamma distribution survival model (#411).
 - Reworked LOO calculations to apply to each individual submodel and disabled LOO calculations for the overall joint model (#402).
 - Added support for additive variance (#403).
 - Added support for independent variances per study/arm (#389).
 - Miscellaneous bug fixes.
+- Introduce new package options for bounds close to zero which are used to avoid MCMC sampler starting warnings.
 
 # jmpost 0.0.1
 

R/DataSurvival.R

@@ -117,15 +117,30 @@ as.data.frame.DataSurvival <- function(x, ...) {
 
 #' @rdname as_stan_list.DataObject
 #' @family DataSurvival
+#' @param df (`data.frame`)\cr The data frame to be used for generating the design matrix.
+#' @param vars (`list`)\cr The list of variable names to be used for generating the design matrix.
 #' @export
-as_stan_list.DataSurvival <- function(object, ...) {
-    df <- as.data.frame(object)
-    vars <- extractVariableNames(object)
-
+model.matrix.DataSurvival <- function(
+    object,
+    df = as.data.frame(object),
+    vars = extractVariableNames(object),
+    ...
+) {
     design_mat <- stats::model.matrix(vars$frm, data = df)
     remove_index <- grep("(Intercept)", colnames(design_mat), fixed = TRUE)
     design_mat <- design_mat[, -remove_index, drop = FALSE]
     rownames(design_mat) <- NULL
+    design_mat
+}
+
+#' @rdname as_stan_list.DataObject
+#' @family DataSurvival
+#' @export
+as_stan_list.DataSurvival <- function(object, ...) {
+    df <- as.data.frame(object)
+    vars <- extractVariableNames(object)
+
+    design_mat <- model.matrix(object, df = df, vars = vars)
 
     # Parameters for efficient integration of hazard function -> survival function
     gh_parameters <- statmod::gauss.quad(

R/Prior.R

@@ -36,6 +36,7 @@ NULL
 #' @slot display (`string`)\cr See arguments.
 #' @slot sample (`function`)\cr See arguments.
 #' @slot limits (`numeric`)\cr See arguments.
+#' @slot .allow_vectors (`logical`)\cr See arguments.
 #'
 #' @family Prior-internal
 #' @export Prior
@@ -50,7 +51,8 @@ NULL
         "centre" = "numeric",
         "validation" = "list",
         "sample" = "function",
-        "limits" = "numeric"
+        "limits" = "numeric",
+        ".allow_vectors" = "logical"
     )
 )
 
@@ -72,6 +74,8 @@ NULL
 #'   a function to sample from the prior distribution.
 #' @typed limits: numeric
 #'   the lower and upper limits for a truncated distribution
+#' @typed .allow_vectors: flag
+#'   whether to allow vector parameters.
 #' @rdname Prior-class
 Prior <- function(
     parameters,
@@ -81,7 +85,8 @@ Prior <- function(
     centre,
     validation,
     sample,
-    limits = c(-Inf, Inf)
+    limits = c(-Inf, Inf),
+    .allow_vectors = FALSE
 ) {
     .Prior(
         parameters = parameters,
@@ -91,7 +96,8 @@ Prior <- function(
         display = display,
         validation = validation,
         sample = sample,
-        limits = limits
+        limits = limits,
+        .allow_vectors = .allow_vectors
     )
 }
 
@@ -106,7 +112,10 @@ setValidity(
                     param
                 ))
             }
-            if (length(object@parameters[[param]]) != 1) {
+            if (
+                !object@.allow_vectors &&
+                    length(object@parameters[[param]]) != 1
+            ) {
                 return(sprintf("Parameter `%s` must be a single value", param))
             }
             if (!object@validation[[param]](object@parameters[[param]])) {
@@ -118,6 +127,18 @@ setValidity(
                 return(return_message)
             }
         }
+        if (object@.allow_vectors) {
+            # Check that all parameters have either same length or are
+            # length 1.
+            all_lengths <- sapply(
+                object@parameters,
+                length
+            )
+            none_one_lengths <- setdiff(unique(all_lengths), 1)
+            if (length(none_one_lengths) > 1) {
+                return("All parameters must be the same length or length 1")
+            }
+        }
         if (length(object@limits) != 2) {
             return("Limits must be a vector of length 2")
         }
@@ -266,21 +287,66 @@ NULL
 #' @describeIn Prior-Getter-Methods The prior's initial value
 #' @export
 initialValues.Prior <- function(object, ...) {
-    samples <- getOption("jmpost.prior_shrinkage") *
-        object@centre +
-        (1 - getOption("jmpost.prior_shrinkage")) * object@sample(100)
+    n_samples <- 100
+    centre_value <- object@centre
+    sample_values <- object@sample(n_samples)
+
+    is_scalar <- length(centre_value) == 1
+    if (is_scalar) {
+        assert_that(
+            length(sample_values) == n_samples,
+            msg = "Sample function must return a vector of length n_samples"
+        )
 
-    valid_samples <- samples[
-        samples >= min(object@limits) & samples <= max(object@limits)
-    ]
-    assert_that(
-        length(valid_samples) >= 1,
-        msg = "Unable to generate an initial value that meets the required constraints"
-    )
-    if (length(valid_samples) == 1) {
-        return(valid_samples)
+        samples <- getOption("jmpost.prior_shrinkage") *
+            object@centre +
+            (1 - getOption("jmpost.prior_shrinkage")) * object@sample(n_samples)
+
+        valid_samples <- samples[
+            samples >= min(object@limits) & samples <= max(object@limits)
+        ]
+        assert_that(
+            length(valid_samples) >= 1,
+            msg = "Unable to generate an initial value that meets the required constraints"
+        )
+        if (length(valid_samples) == 1) {
+            return(valid_samples)
+        }
+        return(sample(valid_samples, 1))
+    } else {
+        n_centre_vals <- length(centre_value)
+        assert_that(
+            ncol(sample_values) == n_centre_vals,
+            msg = paste(
+                "Sample function must return a matrix with n_samples rows and the same",
+                "number of columns as the length of the centre value"
+            )
+        )
+        samples <- getOption("jmpost.prior_shrinkage") *
+            matrix(
+                centre_value,
+                nrow = n_samples,
+                ncol = n_centre_vals,
+                byrow = TRUE
+            ) +
+            (1 - getOption("jmpost.prior_shrinkage")) * sample_values
+        valid_samples <- samples[
+            apply(samples, 1, function(row) {
+                all(row >= min(object@limits) & row <= max(object@limits))
+            }),
+            ,
+            drop = FALSE
+        ]
+        assert_that(
+            nrow(valid_samples) >= 1,
+            msg = "Unable to generate an initial value that meets the required constraints"
+        )
+        if (nrow(valid_samples) == 1) {
+            return(valid_samples[1, ])
+        } else {
+            return(valid_samples[sample(nrow(valid_samples), 1), ])
+        }
     }
-    return(sample(valid_samples, 1))
 }
 
 
@@ -313,6 +379,43 @@ prior_normal <- function(mu, sigma) {
 }
 
 
+#' Normal Prior for a Vector Distribution
+#'
+#' @typed mus: numeric
+#'   means.
+#' @typed sigmas: numeric
+#'   standard deviations.
+#' @family Prior
+#' @export
+prior_normal_vector <- function(mus, sigmas) {
+    Prior(
+        parameters = list(
+            mus = mus,
+            sigmas = sigmas,
+            dim_mus = length(mus),
+            dim_sigmas = length(sigmas)
+        ),
+        display = "normal(mus = [{toString(mus)}], sigmas = [{toString(sigmas)}])",
+        repr_model = "{name} ~ normal(prior_mus_{name}, prior_sigmas_{name})",
+        repr_data = c(
+            "int<lower=1> prior_dim_mus_{name};",
+            "int<lower=1> prior_dim_sigmas_{name};",
+            "vector[prior_dim_mus_{name}] prior_mus_{name};",
+            "vector<lower=0>[prior_dim_sigmas_{name}] prior_sigmas_{name};"
+        ),
+        centre = mus,
+        sample = \(n) local_rnorm_vector(n, mus, sigmas),
+        validation = list(
+            mus = \(x) all(is.numeric(x)),
+            sigmas = \(x) all(x > 0),
+            dim_mus = is.count,
+            dim_sigmas = is.count
+        ),
+        .allow_vectors = TRUE
+    )
+}
+
+
 #' Standard Normal Prior Distribution
 #'
 #'
@@ -719,6 +822,11 @@ NULL
 #' @rdname Local_Sample
 local_rnorm <- \(...) rnorm(...)
 
+#' @rdname Local_Sample
+local_rnorm_vector <- \(n, mus, sigmas) {
+    mapply(local_rnorm, n = n, mean = mus, sd = sigmas)
+}
+
 #' @rdname Local_Sample
 local_rcauchy <- \(...) rcauchy(...)
 

R/settings.R

@@ -42,6 +42,17 @@
 #' higher values of `n` lead to better accuracy of the approximation but at the cost of
 #' increased computational time.
 #'
+#' ## `jmpost.double_eps` and `jmpost.double_neg_eps`
+#'
+#' Default = `sqrt(.Machine$double.eps)` and `sqrt(.Machine$double.neg.eps)`
+#'
+#' These are the smallest positive floating-point numbers `x` which are used
+#' in the parameter definitions as lower boundaries to prevent numerical issues
+#' with the MCMC sampler in the beginning of the sampling process.
+#' The default values are set to the square root of the machine epsilon for
+#' double precision numbers, which is a common choice to balance numerical stability
+#' and the ability to explore the parameter space effectively.
+#'
 #' @examples
 #' \dontrun{
 #' options(jmpost.prior_shrinkage = 0.5)
@@ -58,7 +69,9 @@ set_options <- function() {
     jmpost_opts <- list(
         jmpost.cache_dir = cache_dir,
         jmpost.prior_shrinkage = 0.5,
-        jmpost.gauss_quad_n = 15
+        jmpost.gauss_quad_n = 15,
+        jmpost.double_eps = sqrt(.Machine$double.eps),
+        jmpost.double_neg_eps = sqrt(.Machine$double.neg.eps)
     )
     for (opt in names(jmpost_opts)) {
         if (!opt %in% current_opts) {

R/utilities.R

@@ -207,8 +207,8 @@ samples_median_ci <- function(samples, level = 0.95) {
 decorated_render <- function(...) {
     jinjar::render(
         ...,
-        machine_double_eps = 0,
-        machine_double_neg_eps = 0
+        machine_double_eps = getOption("jmpost.double_eps"),
+        machine_double_neg_eps = getOption("jmpost.double_neg_eps"),
     )
 }
 

_pkgdown.yml

@@ -64,6 +64,7 @@ reference:
       - prior_loglogistic
       - prior_lognormal
       - prior_normal
+      - prior_normal_vector
       - prior_uniform
       - prior_std_normal
       - prior_student_t