From 4c892b6b182212671f773051b314a950c692d4f2 Mon Sep 17 00:00:00 2001 From: gowerc Date: Fri, 2 Feb 2024 16:57:30 +0000 Subject: [PATCH 01/21] Different initial values for each mcmc chain Fixes #4 --- NAMESPACE | 3 + R/JointModel.R | 72 ++++++-- R/Link.R | 5 +- R/LinkGSF.R | 6 +- R/LinkRandomSlope.R | 2 +- R/LongitudinalGSF.R | 40 ++--- R/LongitudinalRandomSlope.R | 10 +- R/Parameter.R | 2 +- R/ParameterList.R | 16 +- R/Prior.R | 157 +++++++++++++----- R/StanModel.R | 6 + R/SurvivalLoglogistic.R | 4 +- R/generics.R | 15 +- design/initial_values.md | 42 ----- man/LinkRandomSlope-class.Rd | 2 +- man/Local_Sample.Rd | 56 +++++++ man/LongitudinalGSF-class.Rd | 34 +--- man/LongitudinalRandomSlope-class.Rd | 10 +- man/Parameter-Getter-Methods.Rd | 2 +- man/ParameterList-Getter-Methods.Rd | 4 +- man/Prior-Getter-Methods.Rd | 4 +- man/Prior-class.Rd | 6 +- man/SurvivalLogLogistic-class.Rd | 4 +- man/initialValues.Rd | 26 ++- man/link_gsf_dsld-class.Rd | 2 +- man/link_gsf_identity-class.Rd | 2 +- man/link_gsf_ttg-class.Rd | 2 +- man/prior_beta.Rd | 4 +- man/prior_cauchy.Rd | 4 +- man/prior_gamma.Rd | 4 +- man/prior_invgamma.Rd | 4 +- man/prior_logistic.Rd | 4 +- man/prior_loglogistic.Rd | 4 +- man/prior_lognormal.Rd | 4 +- man/prior_none.Rd | 5 +- man/prior_normal.Rd | 4 +- man/prior_std_normal.Rd | 5 +- man/prior_student_t.Rd | 4 +- man/prior_uniform.Rd | 4 +- tests/testthat/_snaps/JointModel.md | 8 +- .../_snaps/LongitudinalRandomSlope.md | 6 +- tests/testthat/_snaps/Prior.md | 4 +- tests/testthat/_snaps/SurvivalExponential.md | 2 +- tests/testthat/_snaps/SurvivalLoglogistic.md | 3 +- tests/testthat/_snaps/SurvivalWeibullPH.md | 3 +- tests/testthat/test-LinkRandomSlope.R | 13 +- tests/testthat/test-Parameter.R | 19 ++- tests/testthat/test-ParameterList.R | 20 ++- tests/testthat/test-Prior.R | 40 ++--- tests/testthat/test-StanModel.R | 2 +- tests/testthat/test-SurvivalExponential.R | 2 +- tests/testthat/test-SurvivalLoglogistic.R | 2 +- tests/testthat/test-SurvivalWeibullPH.R | 2 +- tests/testthat/test-initialValues.R | 83 +++++++++ vignettes/model_fitting.Rmd | 42 +++++ 55 files changed, 551 insertions(+), 284 deletions(-) delete mode 100644 design/initial_values.md create mode 100644 man/Local_Sample.Rd create mode 100644 tests/testthat/test-initialValues.R diff --git a/NAMESPACE b/NAMESPACE index 2dcd3bfe..799c4438 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -52,9 +52,11 @@ S3method(extractVariableNames,DataSurvival) S3method(generateQuantities,JointModelSamples) S3method(getParameters,default) S3method(initialValues,JointModel) +S3method(initialValues,Link) S3method(initialValues,Parameter) S3method(initialValues,ParameterList) S3method(initialValues,Prior) +S3method(initialValues,StanModel) S3method(names,Parameter) S3method(names,ParameterList) S3method(sampleStanModel,JointModel) @@ -100,6 +102,7 @@ export(generateQuantities) export(gsf_dsld) export(gsf_sld) export(gsf_ttg) +export(initialValues) export(link_gsf_abstract) export(link_gsf_dsld) export(link_gsf_identity) diff --git a/R/JointModel.R b/R/JointModel.R index 20d634de..c09816e6 100755 --- a/R/JointModel.R +++ b/R/JointModel.R @@ -136,6 +136,8 @@ compileStanModel.JointModel <- function(object) { #' @export sampleStanModel.JointModel <- function(object, data, ...) { + assert_class(data, "DataJoint") + if (!is.null(object@survival)) { assert_that( !is.null(data@survival), @@ -150,21 +152,36 @@ sampleStanModel.JointModel <- function(object, data, ...) { } args <- list(...) + args[["data"]] <- append( as_stan_list(data), as_stan_list(object@parameters) ) - if (!"init" %in% names(args)) { - values_initial <- initialValues(object) - values_sizes <- size(object@parameters) - values_sizes_complete <- replace_with_lookup(values_sizes, args[["data"]]) - values_initial_expanded <- expand_initial_values(values_initial, values_sizes_complete) - args[["init"]] <- function() values_initial_expanded + args[["chains"]] <- if ("chains" %in% names(args)) { + args[["chains"]] + } else { + 4 } + initial_values <- if ("init" %in% names(args)) { + args[["init"]] + } else { + initialValues(object, nchains = args[["chains"]]) + } + + args[["init"]] <- ensure_initial_values( + initial_values, + args[["data"]], + object@parameters + ) + model <- compileStanModel(object) - results <- do.call(model$sample, args) + + results <- do.call( + model$sample, + args + ) .JointModelSamples( model = object, @@ -174,12 +191,47 @@ sampleStanModel.JointModel <- function(object, data, ...) { } -# initialValues-JointModel ---- +#' Ensure that initial values are correctly specified +#' +#' @param initial_values (`list`)\cr A list of lists containing the initial values +#' must be 1 list per desired chain. All elements should have identical names +#' @param data (`list`)\cr specifies the size to expand each of our initial values to be. +#' That is elements of size 1 in `initial_values` will be expanded to be the same +#' size as the corresponding element in `data` by broadcasting the value. +#' @param parameters ([`ParameterList`])\cr the parameters object +#' +#' @details +#' This function is mostly a thin wrapper around `expand_initial_values` to +#' enable easier unit testing. +#' +#' @keywords internal +ensure_initial_values <- function(initial_values, data, parameters) { + if (is.function(initial_values)) { + return(initial_values) + } + + assert_class(data, "list") + assert_class(parameters, "ParameterList") + assert_class(initial_values, "list") + + values_sizes <- size(parameters) + values_sizes_complete <- replace_with_lookup( + values_sizes, + data + ) + lapply( + initial_values, + expand_initial_values, + sizes = values_sizes_complete + ) +} + + #' @rdname initialValues #' @export -initialValues.JointModel <- function(object) { - initialValues(object@parameters) +initialValues.JointModel <- function(object, nchains, ...) { + initialValues(object@parameters, nchains) } diff --git a/R/Link.R b/R/Link.R index 663bb950..ea845168 100755 --- a/R/Link.R +++ b/R/Link.R @@ -59,8 +59,9 @@ setMethod( # initialValues-Link ---- #' @rdname initialValues -initialValues.Link <- function(object) { - initialValues(object@parameters) +#' @export +initialValues.Link <- function(object, nchains, ...) { + initialValues(object@parameters, nchains) } diff --git a/R/LinkGSF.R b/R/LinkGSF.R index 07fe2126..5a154eb7 100644 --- a/R/LinkGSF.R +++ b/R/LinkGSF.R @@ -149,7 +149,7 @@ link_gsf_abstract <- function( #' #' @export link_gsf_ttg <- function( - gamma = prior_normal(0, 5, init = 0) + gamma = prior_normal(0, 5) ) { .link_gsf_ttg( name = "TTG", @@ -182,7 +182,7 @@ link_gsf_ttg <- function( #' #' @export link_gsf_dsld <- function( - beta = prior_normal(0, 5, init = 0) + beta = prior_normal(0, 5) ) { .link_gsf_dsld( name = "dSLD", @@ -215,7 +215,7 @@ link_gsf_dsld <- function( #' @param tau (`Prior`)\cr prior for the link coefficient `tau`. #' #' @export -link_gsf_identity <- function(tau = prior_normal(0, 5, init = 0)) { +link_gsf_identity <- function(tau = prior_normal(0, 5)) { .link_gsf_identity( name = "Identity", stan = StanModule("lm-gsf/link_identity.stan"), diff --git a/R/LinkRandomSlope.R b/R/LinkRandomSlope.R index c5ebb482..79c5f24f 100644 --- a/R/LinkRandomSlope.R +++ b/R/LinkRandomSlope.R @@ -23,7 +23,7 @@ NULL #' #' @export LinkRandomSlope <- function( - link_lm_phi = prior_normal(0.2, 0.5, init = 0.02) + link_lm_phi = prior_normal(0.2, 0.5) ) { .LinkRandomSlope( Link( diff --git a/R/LongitudinalGSF.R b/R/LongitudinalGSF.R index f0dd0e97..6a85fa55 100755 --- a/R/LongitudinalGSF.R +++ b/R/LongitudinalGSF.R @@ -36,37 +36,23 @@ NULL #' @param a_phi (`Prior`)\cr for the alpha parameter for the fraction of cells that respond to treatment. #' @param b_phi (`Prior`)\cr for the beta parameter for the fraction of cells that respond to treatment. #' -#' @param psi_bsld (`Prior`)\cr for the baseline value random effect `psi_bsld`. Only used in the -#' centered parameterization to set the initial value. -#' @param psi_ks (`Prior`)\cr for the shrinkage rate random effect `psi_ks`. Only used in the -#' centered parameterization to set the initial value. -#' @param psi_kg (`Prior`)\cr for the growth rate random effect `psi_kg`. Only used in the -#' centered parameterization to set the initial value. -#' @param psi_phi (`Prior`)\cr for the shrinkage proportion random effect `psi_phi`. Only used in the -#' centered parameterization to set the initial value. -#' #' @param centered (`logical`)\cr whether to use the centered parameterization. #' #' @export LongitudinalGSF <- function( - mu_bsld = prior_normal(log(60), 1, init = 60), - mu_ks = prior_normal(log(0.5), 1, init = 0.5), - mu_kg = prior_normal(log(0.3), 1, init = 0.3), - - omega_bsld = prior_lognormal(log(0.2), 1, init = 0.2), - omega_ks = prior_lognormal(log(0.2), 1, init = 0.2), - omega_kg = prior_lognormal(log(0.2), 1, init = 0.2), + mu_bsld = prior_normal(log(60), 1), + mu_ks = prior_normal(log(0.5), 1), + mu_kg = prior_normal(log(0.3), 1), - a_phi = prior_lognormal(log(5), 1, init = 5), - b_phi = prior_lognormal(log(5), 1, init = 5), + omega_bsld = prior_lognormal(log(0.2), 1), + omega_ks = prior_lognormal(log(0.2), 1), + omega_kg = prior_lognormal(log(0.2), 1), - sigma = prior_lognormal(log(0.1), 1, init = 0.1), + a_phi = prior_lognormal(log(5), 1), + b_phi = prior_lognormal(log(5), 1), - psi_bsld = prior_none(init = 60), - psi_ks = prior_none(init = 0.5), - psi_kg = prior_none(init = 0.5), - psi_phi = prior_none(init = 0.5), + sigma = prior_lognormal(log(0.1), 1), centered = FALSE ) { @@ -87,7 +73,7 @@ LongitudinalGSF <- function( Parameter(name = "lm_gsf_a_phi", prior = a_phi, size = "n_arms"), Parameter(name = "lm_gsf_b_phi", prior = b_phi, size = "n_arms"), - Parameter(name = "lm_gsf_psi_phi", prior = psi_phi, size = "Nind"), + Parameter(name = "lm_gsf_psi_phi", prior = prior_none(), size = "Nind"), Parameter(name = "lm_gsf_sigma", prior = sigma, size = 1) ) @@ -95,9 +81,9 @@ LongitudinalGSF <- function( assert_flag(centered) parameters_extra <- if (centered) { list( - Parameter(name = "lm_gsf_psi_bsld", prior = psi_bsld, size = "Nind"), - Parameter(name = "lm_gsf_psi_ks", prior = psi_ks, size = "Nind"), - Parameter(name = "lm_gsf_psi_kg", prior = psi_kg, size = "Nind") + Parameter(name = "lm_gsf_psi_bsld", prior = prior_none(), size = "Nind"), + Parameter(name = "lm_gsf_psi_ks", prior = prior_none(), size = "Nind"), + Parameter(name = "lm_gsf_psi_kg", prior = prior_none(), size = "Nind") ) } else { list( diff --git a/R/LongitudinalRandomSlope.R b/R/LongitudinalRandomSlope.R index 9002f9fe..a2fafc0d 100755 --- a/R/LongitudinalRandomSlope.R +++ b/R/LongitudinalRandomSlope.R @@ -28,11 +28,11 @@ NULL #' #' @export LongitudinalRandomSlope <- function( - intercept = prior_normal(30, 10, init = 30), - slope_mu = prior_normal(0, 15, init = 0.001), - slope_sigma = prior_lognormal(1, 5, init = 1), - sigma = prior_lognormal(1, 5, init = 1), - random_slope = prior_none(init = 0) + intercept = prior_normal(30, 10), + slope_mu = prior_normal(0, 15), + slope_sigma = prior_lognormal(0, 1.5), + sigma = prior_lognormal(0, 1.5), + random_slope = prior_none() ) { stan <- StanModule( diff --git a/R/Parameter.R b/R/Parameter.R index 9061ad4c..8bdf5bac 100755 --- a/R/Parameter.R +++ b/R/Parameter.R @@ -115,7 +115,7 @@ names.Parameter <- function(x) x@name #' @describeIn Parameter-Getter-Methods The parameter's initial values #' @export -initialValues.Parameter <- function(object) initialValues(object@prior) +initialValues.Parameter <- function(object, ...) initialValues(object@prior) #' @describeIn Parameter-Getter-Methods The parameter's dimensionality #' @export diff --git a/R/ParameterList.R b/R/ParameterList.R index d820129c..16fb4961 100644 --- a/R/ParameterList.R +++ b/R/ParameterList.R @@ -145,11 +145,17 @@ names.ParameterList <- function(x) { #' @describeIn ParameterList-Getter-Methods The parameter-list's parameter initial values #' @export -initialValues.ParameterList <- function(object) { - vals <- lapply(object@parameters, initialValues) - name <- vapply(object@parameters, names, character(1)) - names(vals) <- name - return(vals) +initialValues.ParameterList <- function(object, nchains, ...) { + x <- lapply( + seq_len(nchains), + \(i) { + vals <- lapply(object@parameters, initialValues) + name <- vapply(object@parameters, names, character(1)) + names(vals) <- name + vals + } + ) + return(x) } diff --git a/R/Prior.R b/R/Prior.R index 4135ba9f..7abac649 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -29,6 +29,7 @@ NULL #' @slot init (`numeric`)\cr See arguments. #' @slot validation (`list`)\cr See arguments. #' @slot display (`string`)\cr See arguments. +#' @slot sample (`function`)\cr See arguments. #' #' @family Prior-internal #' @export Prior @@ -41,7 +42,8 @@ NULL "repr_model" = "character", "repr_data" = "character", "init" = "numeric", - "validation" = "list" + "validation" = "list", + "sample" = "function" ) ) @@ -53,15 +55,25 @@ NULL #' @param init (`numeric`)\cr the initial value. #' @param validation (`list`)\cr the prior distribution parameter validation functions. Must have #' the same names as the `paramaters` slot. +#' @param sample (`function`)\cr a function to sample from the prior distribution. #' @rdname Prior-class -Prior <- function(parameters, display, repr_model, repr_data, init, validation) { +Prior <- function( + parameters, + display, + repr_model, + repr_data, + init, + validation, + sample +) { .Prior( parameters = parameters, repr_model = repr_model, repr_data = repr_data, init = init, display = display, - validation = validation + validation = validation, + sample = sample ) } @@ -175,7 +187,9 @@ NULL #' @describeIn Prior-Getter-Methods The prior's initial value #' @export -initialValues.Prior <- function(object) object@init +initialValues.Prior <- function(object, ...) { + 0.5 * object@init + 0.5 * object@sample(1) +} # Prior-constructors ---- @@ -184,10 +198,9 @@ initialValues.Prior <- function(object) object@init #' #' @param mu (`number`)\cr mean. #' @param sigma (`number`)\cr standard deviation. -#' @inheritParams Prior-Shared #' @family Prior #' @export -prior_normal <- function(mu, sigma, init = mu) { +prior_normal <- function(mu, sigma) { Prior( parameters = list(mu = mu, sigma = sigma), display = "normal(mu = {mu}, sigma = {sigma})", @@ -196,7 +209,8 @@ prior_normal <- function(mu, sigma, init = mu) { "real prior_mu_{name};", "real prior_sigma_{name};" ), - init = init, + init = mu, + sample = \(n) local_rnorm(n, mu, sigma), validation = list( mu = is.numeric, sigma = \(x) x > 0 @@ -207,17 +221,17 @@ prior_normal <- function(mu, sigma, init = mu) { #' Standard Normal Prior Distribution #' -#' @inheritParams Prior-Shared #' #' @family Prior #' @export -prior_std_normal <- function(init = 0) { +prior_std_normal <- function() { Prior( parameters = list(), display = "std_normal()", repr_model = "{name} ~ std_normal();", repr_data = "", - init = init, + init = 0, + sample = \(n) local_rnorm(n), validation = list() ) } @@ -226,11 +240,10 @@ prior_std_normal <- function(init = 0) { #' #' @param mu (`number`)\cr mean. #' @param sigma (`number`)\cr scale. -#' @inheritParams Prior-Shared #' @family Prior #' #' @export -prior_cauchy <- function(mu, sigma, init = mu) { +prior_cauchy <- function(mu, sigma) { Prior( parameters = list(mu = mu, sigma = sigma), display = "cauchy(mu = {mu}, sigma = {sigma})", @@ -239,7 +252,8 @@ prior_cauchy <- function(mu, sigma, init = mu) { "real prior_mu_{name};", "real prior_sigma_{name};" ), - init = init, + init = mu, + sample = \(n) local_rcauchy(n, mu, sigma), validation = list( mu = is.numeric, sigma = \(x) x > 0 @@ -251,11 +265,10 @@ prior_cauchy <- function(mu, sigma, init = mu) { #' #' @param alpha (`number`)\cr shape. #' @param beta (`number`)\cr inverse scale. -#' @inheritParams Prior-Shared #' @family Prior #' #' @export -prior_gamma <- function(alpha, beta, init = alpha / beta) { +prior_gamma <- function(alpha, beta) { Prior( parameters = list(alpha = alpha, beta = beta), repr_model = "{name} ~ gamma(prior_alpha_{name}, prior_beta_{name});", @@ -264,7 +277,8 @@ prior_gamma <- function(alpha, beta, init = alpha / beta) { "real prior_alpha_{name};", "real prior_beta_{name};" ), - init = init, + init = alpha / beta, + sample = \(n) local_rgamma(n, shape = alpha, rate = beta), validation = list( alpha = \(x) x > 0, beta = \(x) x > 0 @@ -276,11 +290,10 @@ prior_gamma <- function(alpha, beta, init = alpha / beta) { #' #' @param mu (`number`)\cr mean of the logarithm. #' @param sigma (`number`)\cr standard deviation of the logarithm. -#' @inheritParams Prior-Shared #' @family Prior #' #' @export -prior_lognormal <- function(mu, sigma, init = exp(mu + (sigma^2) / 2)) { +prior_lognormal <- function(mu, sigma) { Prior( parameters = list(mu = mu, sigma = sigma), display = "lognormal(mu = {mu}, sigma = {sigma})", @@ -289,7 +302,8 @@ prior_lognormal <- function(mu, sigma, init = exp(mu + (sigma^2) / 2)) { "real prior_mu_{name};", "real prior_sigma_{name};" ), - init = init, + init = exp(mu + (sigma^2) / 2), + sample = \(n) local_rlnorm(n, mu, sigma), validation = list( mu = is.numeric, sigma = \(x) x > 0 @@ -301,11 +315,10 @@ prior_lognormal <- function(mu, sigma, init = exp(mu + (sigma^2) / 2)) { #' #' @param a (`number`)\cr first parameter. #' @param b (`number`)\cr second parameter -#' @inheritParams Prior-Shared #' @family Prior #' #' @export -prior_beta <- function(a, b, init = a / (a + b)) { +prior_beta <- function(a, b) { Prior( parameters = list(a = a, b = b), display = "beta(a = {a}, b = {b})", @@ -314,7 +327,8 @@ prior_beta <- function(a, b, init = a / (a + b)) { "real prior_a_{name};", "real prior_b_{name};" ), - init = init, + init = a / (a + b), + sample = \(n) local_rbeta(n, a, b), validation = list( a = \(x) x > 0, b = \(x) x > 0 @@ -324,17 +338,17 @@ prior_beta <- function(a, b, init = a / (a + b)) { #' Only Initial Values Specification #' -#' @inheritParams Prior-Shared #' @family Prior #' #' @export -prior_none <- function(init = 0.00001) { +prior_none <- function() { Prior( parameters = list(), display = "", repr_model = "", repr_data = "", - init = init, + sample = \(n) local_runif(n, -4, 4), + init = 0, validation = list() ) } @@ -346,11 +360,10 @@ prior_none <- function(init = 0.00001) { #' #' @param alpha (`number`)\cr minimum value parameter. #' @param beta (`number`)\cr maximum value parameter. -#' @inheritParams Prior-Shared #' @family Prior #' #' @export -prior_uniform <- function(alpha, beta, init = 0.5 * (alpha + beta)) { +prior_uniform <- function(alpha, beta) { Prior( parameters = list(alpha = alpha, beta = beta), display = "uniform(alpha = {alpha}, beta = {beta})", @@ -359,7 +372,8 @@ prior_uniform <- function(alpha, beta, init = 0.5 * (alpha + beta)) { "real prior_alpha_{name};", "real prior_beta_{name};" ), - init = init, + init = 0.5 * (alpha + beta), + sample = \(n) local_runif(n, alpha, beta), validation = list( alpha = is.numeric, beta = is.numeric @@ -373,11 +387,10 @@ prior_uniform <- function(alpha, beta, init = 0.5 * (alpha + beta)) { #' @param nu (`number`)\cr Degrees of freedom parameter. #' @param mu (`number`)\cr Location parameter. #' @param sigma (`number`)\cr Scale parameter. -#' @inheritParams Prior-Shared #' @family Prior #' #' @export -prior_student_t <- function(nu, mu, sigma, init = mu) { +prior_student_t <- function(nu, mu, sigma) { Prior( parameters = list( nu = nu, @@ -391,7 +404,8 @@ prior_student_t <- function(nu, mu, sigma, init = mu) { "real prior_mu_{name};", "real prior_sigma_{name};" ), - init = init, + init = mu, + sample = \(n) local_rt(n, nu, mu, sigma), validation = list( nu = \(x) x > 0, mu = is.numeric, @@ -401,15 +415,15 @@ prior_student_t <- function(nu, mu, sigma, init = mu) { } + #' Logistic Prior Distribution #' #' @param mu (`number`)\cr Location parameter. #' @param sigma (`number`)\cr Scale parameter. -#' @inheritParams Prior-Shared #' @family Prior #' #' @export -prior_logistic <- function(mu, sigma, init = mu) { +prior_logistic <- function(mu, sigma) { Prior( parameters = list( mu = mu, @@ -421,7 +435,8 @@ prior_logistic <- function(mu, sigma, init = mu) { "real prior_mu_{name};", "real prior_sigma_{name};" ), - init = init, + init = mu, + sample = \(n) local_rlogis(n, mu, sigma), validation = list( mu = is.numeric, sigma = \(x) x > 0 @@ -434,11 +449,10 @@ prior_logistic <- function(mu, sigma, init = mu) { #' #' @param alpha (`number`)\cr Scale parameter. #' @param beta (`number`)\cr Shape parameter. -#' @inheritParams Prior-Shared #' @family Prior #' #' @export -prior_loglogistic <- function(alpha, beta, init = alpha * pi / (beta * sin(pi / beta))) { +prior_loglogistic <- function(alpha, beta) { Prior( parameters = list( alpha = alpha, @@ -450,7 +464,10 @@ prior_loglogistic <- function(alpha, beta, init = alpha * pi / (beta * sin(pi / "real prior_alpha_{name};", "real prior_beta_{name};" ), - init = init, + init = alpha * pi / (beta * sin(pi / beta)), + sample = \(n) { + local_rloglogis(n, alpha, beta) + }, validation = list( alpha = \(x) x > 0, beta = \(x) x > 0 @@ -463,11 +480,10 @@ prior_loglogistic <- function(alpha, beta, init = alpha * pi / (beta * sin(pi / #' #' @param alpha (`number`)\cr Shape parameter. #' @param beta (`number`)\cr Scale parameter. -#' @inheritParams Prior-Shared #' @family Prior #' #' @export -prior_invgamma <- function(alpha, beta, init = beta / (alpha - 1)) { +prior_invgamma <- function(alpha, beta) { Prior( parameters = list( alpha = alpha, @@ -479,10 +495,71 @@ prior_invgamma <- function(alpha, beta, init = beta / (alpha - 1)) { "real prior_alpha_{name};", "real prior_beta_{name};" ), - init = init, + init = beta / (alpha - 1), + sample = \(n) local_rinvgamma(n, alpha, beta), validation = list( alpha = \(x) x > 0, beta = \(x) x > 0 ) ) } + + + + + +#' Stub functions for sampling from distributions +#' +#' @description +#' These functions only exist so that they can be mocked during unit +#' tests in order to provide deterministic values. In most cases +#' these are just straight forward pass throughs for the underlying +#' distributions. +#' +#' @param alpha (`number`)\cr Parameter for underlying distribution. +#' @param beta (`number`)\cr Parameter for underlying distribution. +#' @param mu (`number`)\cr Parameter for underlying distribution. +#' @param sigma (`number`)\cr Parameter for underlying distribution. +#' @param nu (`number`)\cr Parameter for underlying distribution. +#' @param ... Pass any additional arguments to the underlying distribution. +#' +#' @name Local_Sample +#' @keywords internal +NULL + +#' @rdname Local_Sample +local_rnorm <- \(...) rnorm(...) + +#' @rdname Local_Sample +local_rcauchy <- \(...) rcauchy(...) + +#' @rdname Local_Sample +local_rgamma <- \(...) rgamma(...) + +#' @rdname Local_Sample +local_rlnorm <- \(...) rlnorm(...) + +#' @rdname Local_Sample +local_rbeta <- \(...) rbeta(...) + +#' @rdname Local_Sample +local_runif <- \(...) runif(...) + +#' @rdname Local_Sample +local_rt <- \(n, nu, mu, sigma) { + rt(n, nu) * sigma + mu +} + +#' @rdname Local_Sample +local_rlogis <- \(...) rlogis(...) + +#' @rdname Local_Sample +local_rloglogis <- \(n, alpha, beta) { + r <- runif(n) + ((r / (1 - r))^(1 / beta)) * alpha +} + +#' @rdname Local_Sample +local_rinvgamma <- \(n, alpha, beta) { + 1 / rgamma(n, alpha, rate = beta) +} diff --git a/R/StanModel.R b/R/StanModel.R index c3c274b7..703d59c3 100644 --- a/R/StanModel.R +++ b/R/StanModel.R @@ -90,3 +90,9 @@ setMethod( cat(as_print_string(object)) } ) + +#' @rdname initialValues +#' @export +initialValues.StanModel <- function(object, nchains, ...) { + initialValues(object@parameters, nchains) +} diff --git a/R/SurvivalLoglogistic.R b/R/SurvivalLoglogistic.R index 50230c4b..f872b1a7 100644 --- a/R/SurvivalLoglogistic.R +++ b/R/SurvivalLoglogistic.R @@ -24,8 +24,8 @@ NULL #' #' @export SurvivalLogLogistic <- function( - lambda = prior_lognormal(log(0.1), 5, init = 0.1), - p = prior_gamma(2, 5, init = 0.5), + lambda = prior_lognormal(log(0.1), 5), + p = prior_gamma(2, 5), beta = prior_normal(0, 5) ) { .SurvivalLogLogistic( diff --git a/R/generics.R b/R/generics.R index e273f953..3bfd281e 100755 --- a/R/generics.R +++ b/R/generics.R @@ -151,9 +151,20 @@ extractVariableNames <- function(object) { #' Obtain the `list` of initial values to be passed to the Stan sampler. #' #' @param object where to get the initial values from. +#' @param nchains the number of initial values to generate. See details. +#' @param ... Not currently used. #' -#' @keywords internal -initialValues <- function(object) { +#' @details +#' There are multiple ways of specifying initial values to Stan, see the `init` argument +#' in [cmdstanr::model-method-sample] for full details. Within this package we supply +#' initial values via a list of lists where each inner list contains the initial values +#' for a single chain. As such the `nchains` argument specifies the number of inner lists +#' to generate. +#' +#' See the Vignette for further details of how to specify initial values. +#' +#' @export +initialValues <- function(object, ...) { UseMethod("initialValues") } diff --git a/design/initial_values.md b/design/initial_values.md deleted file mode 100644 index e08e91e8..00000000 --- a/design/initial_values.md +++ /dev/null @@ -1,42 +0,0 @@ - -# Initial Values - -You can specify initial values in one of three ways: - -1. Use the default. By default all initial values are set to the mean of the prior distribution. That is if you use a `prior_gamma(1, 3)` the initial value will be set to `1/3`. - -2. Manually specify in the prior function. That is you can set initial values via `prior_normal(0, 1, init = 20)`. Note that if the parameter represents multiple parameters (e.g. 1 per arm) then single values are replicated for each of the parameters. You can specify different initival values by instead specifying a vector e.g. `prior_normal(0, 1, init = c(10, 20))`. Care is needed to ensure the correct dimensionality. - -3. Provide a list of initial values to the sampling function. Instead of relying on `jmpost` to generate the initial values for you, you can instead specify a list of values directly to the sampling function e.g. - -```R -initial_values <- list("mu" = c(10, 15), "sigma" = 0.4) -samples <- mp <- sampleStanModel( - JointModel(...), - data = DataJoint(...), - init = initial_values -) -``` - -However ensuring you have set initial values for all parameters and ensuring you have the parameter names specified correctly can be tricky. To assist with this it is recommend to use the `jmpost` generated initial values as a starting point e.g. - -```R -jm <- JointModel(...) -initial_values <- initialValues(jm) -initial_values$mu <- c(10, 15) -samples <- mp <- sampleStanModel( - jm, - data = DataJoint(...), - init = initial_values -) -``` - - - - - - - - - - diff --git a/man/LinkRandomSlope-class.Rd b/man/LinkRandomSlope-class.Rd index 709c5f2a..56d39f74 100644 --- a/man/LinkRandomSlope-class.Rd +++ b/man/LinkRandomSlope-class.Rd @@ -7,7 +7,7 @@ \alias{LinkRandomSlope} \title{\code{LinkRandomSlope}} \usage{ -LinkRandomSlope(link_lm_phi = prior_normal(0.2, 0.5, init = 0.02)) +LinkRandomSlope(link_lm_phi = prior_normal(0.2, 0.5)) } \arguments{ \item{link_lm_phi}{(\code{Prior})\cr prior for the link coefficient for the diff --git a/man/Local_Sample.Rd b/man/Local_Sample.Rd new file mode 100644 index 00000000..ddf7f12b --- /dev/null +++ b/man/Local_Sample.Rd @@ -0,0 +1,56 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Prior.R +\name{Local_Sample} +\alias{Local_Sample} +\alias{local_rnorm} +\alias{local_rcauchy} +\alias{local_rgamma} +\alias{local_rlnorm} +\alias{local_rbeta} +\alias{local_runif} +\alias{local_rt} +\alias{local_rlogis} +\alias{local_rloglogis} +\alias{local_rinvgamma} +\title{Stub functions for sampling from distributions} +\usage{ +local_rnorm(...) + +local_rcauchy(...) + +local_rgamma(...) + +local_rlnorm(...) + +local_rbeta(...) + +local_runif(...) + +local_rt(n, nu, mu, sigma) + +local_rlogis(...) + +local_rloglogis(n, alpha, beta) + +local_rinvgamma(n, alpha, beta) +} +\arguments{ +\item{...}{Pass any additional arguments to the underlying distribution.} + +\item{nu}{(\code{number})\cr Parameter for underlying distribution.} + +\item{mu}{(\code{number})\cr Parameter for underlying distribution.} + +\item{sigma}{(\code{number})\cr Parameter for underlying distribution.} + +\item{alpha}{(\code{number})\cr Parameter for underlying distribution.} + +\item{beta}{(\code{number})\cr Parameter for underlying distribution.} +} +\description{ +These functions only exist so that they can be mocked during unit +tests in order to provide deterministic values. In most cases +these are just straight forward pass throughs for the underlying +distributions. +} +\keyword{internal} diff --git a/man/LongitudinalGSF-class.Rd b/man/LongitudinalGSF-class.Rd index bb384114..dd026fd2 100644 --- a/man/LongitudinalGSF-class.Rd +++ b/man/LongitudinalGSF-class.Rd @@ -8,19 +8,15 @@ \title{\code{LongitudinalGSF}} \usage{ LongitudinalGSF( - mu_bsld = prior_normal(log(60), 1, init = 60), - mu_ks = prior_normal(log(0.5), 1, init = 0.5), - mu_kg = prior_normal(log(0.3), 1, init = 0.3), - omega_bsld = prior_lognormal(log(0.2), 1, init = 0.2), - omega_ks = prior_lognormal(log(0.2), 1, init = 0.2), - omega_kg = prior_lognormal(log(0.2), 1, init = 0.2), - a_phi = prior_lognormal(log(5), 1, init = 5), - b_phi = prior_lognormal(log(5), 1, init = 5), - sigma = prior_lognormal(log(0.1), 1, init = 0.1), - psi_bsld = prior_none(init = 60), - psi_ks = prior_none(init = 0.5), - psi_kg = prior_none(init = 0.5), - psi_phi = prior_none(init = 0.5), + mu_bsld = prior_normal(log(60), 1), + mu_ks = prior_normal(log(0.5), 1), + mu_kg = prior_normal(log(0.3), 1), + omega_bsld = prior_lognormal(log(0.2), 1), + omega_ks = prior_lognormal(log(0.2), 1), + omega_kg = prior_lognormal(log(0.2), 1), + a_phi = prior_lognormal(log(5), 1), + b_phi = prior_lognormal(log(5), 1), + sigma = prior_lognormal(log(0.1), 1), centered = FALSE ) } @@ -43,18 +39,6 @@ LongitudinalGSF( \item{sigma}{(\code{Prior})\cr for the variance of the longitudinal values \code{sigma}.} -\item{psi_bsld}{(\code{Prior})\cr for the baseline value random effect \code{psi_bsld}. Only used in the -centered parameterization to set the initial value.} - -\item{psi_ks}{(\code{Prior})\cr for the shrinkage rate random effect \code{psi_ks}. Only used in the -centered parameterization to set the initial value.} - -\item{psi_kg}{(\code{Prior})\cr for the growth rate random effect \code{psi_kg}. Only used in the -centered parameterization to set the initial value.} - -\item{psi_phi}{(\code{Prior})\cr for the shrinkage proportion random effect \code{psi_phi}. Only used in the -centered parameterization to set the initial value.} - \item{centered}{(\code{logical})\cr whether to use the centered parameterization.} } \description{ diff --git a/man/LongitudinalRandomSlope-class.Rd b/man/LongitudinalRandomSlope-class.Rd index ab46fd39..cb99a5c9 100644 --- a/man/LongitudinalRandomSlope-class.Rd +++ b/man/LongitudinalRandomSlope-class.Rd @@ -8,11 +8,11 @@ \title{\code{LongitudinalRandomSlope}} \usage{ LongitudinalRandomSlope( - intercept = prior_normal(30, 10, init = 30), - slope_mu = prior_normal(0, 15, init = 0.001), - slope_sigma = prior_lognormal(1, 5, init = 1), - sigma = prior_lognormal(1, 5, init = 1), - random_slope = prior_none(init = 0) + intercept = prior_normal(30, 10), + slope_mu = prior_normal(0, 15), + slope_sigma = prior_lognormal(0, 1.5), + sigma = prior_lognormal(0, 1.5), + random_slope = prior_none() ) } \arguments{ diff --git a/man/Parameter-Getter-Methods.Rd b/man/Parameter-Getter-Methods.Rd index 4cc1e88b..f9b65e2a 100644 --- a/man/Parameter-Getter-Methods.Rd +++ b/man/Parameter-Getter-Methods.Rd @@ -9,7 +9,7 @@ \usage{ \method{names}{Parameter}(x) -\method{initialValues}{Parameter}(object) +\method{initialValues}{Parameter}(object, ...) \method{size}{Parameter}(object) } diff --git a/man/ParameterList-Getter-Methods.Rd b/man/ParameterList-Getter-Methods.Rd index be530325..4e85c466 100644 --- a/man/ParameterList-Getter-Methods.Rd +++ b/man/ParameterList-Getter-Methods.Rd @@ -9,7 +9,7 @@ \usage{ \method{names}{ParameterList}(x) -\method{initialValues}{ParameterList}(object) +\method{initialValues}{ParameterList}(object, nchains, ...) \method{size}{ParameterList}(object) } @@ -17,6 +17,8 @@ \item{x}{(\code{ParameterList}) \cr A List of \code{\link{Parameter}} Objects.} \item{object}{(\code{ParameterList}) \cr A List of \code{\link{Parameter}} Objects.} + +\item{...}{Not Used.} } \description{ Getter functions for the slots of a \code{\link{ParameterList}} object diff --git a/man/Prior-Getter-Methods.Rd b/man/Prior-Getter-Methods.Rd index e86d2ac1..b9774f85 100644 --- a/man/Prior-Getter-Methods.Rd +++ b/man/Prior-Getter-Methods.Rd @@ -5,10 +5,12 @@ \alias{initialValues.Prior} \title{Prior Getter Functions} \usage{ -\method{initialValues}{Prior}(object) +\method{initialValues}{Prior}(object, ...) } \arguments{ \item{object}{(\code{\link{Prior}})\cr a prior Distribution} + +\item{...}{Not Used.} } \description{ Getter functions for the slots of a \code{\link{Prior}} object diff --git a/man/Prior-class.Rd b/man/Prior-class.Rd index 92be6732..d3e8a522 100644 --- a/man/Prior-class.Rd +++ b/man/Prior-class.Rd @@ -7,7 +7,7 @@ \alias{Prior} \title{Prior Object and Constructor Function} \usage{ -Prior(parameters, display, repr_model, repr_data, init, validation) +Prior(parameters, display, repr_model, repr_data, init, validation, sample) } \arguments{ \item{parameters}{(\code{list})\cr the prior distribution parameters.} @@ -22,6 +22,8 @@ Prior(parameters, display, repr_model, repr_data, init, validation) \item{validation}{(\code{list})\cr the prior distribution parameter validation functions. Must have the same names as the \code{paramaters} slot.} + +\item{sample}{(\code{function})\cr a function to sample from the prior distribution.} } \description{ Specifies the prior distribution in a Stan Model @@ -40,6 +42,8 @@ Specifies the prior distribution in a Stan Model \item{\code{validation}}{(\code{list})\cr See arguments.} \item{\code{display}}{(\code{string})\cr See arguments.} + +\item{\code{sample}}{(\code{function})\cr See arguments.} }} \seealso{ diff --git a/man/SurvivalLogLogistic-class.Rd b/man/SurvivalLogLogistic-class.Rd index 9af00d0d..0ebd7886 100644 --- a/man/SurvivalLogLogistic-class.Rd +++ b/man/SurvivalLogLogistic-class.Rd @@ -8,8 +8,8 @@ \title{\code{SurvivalLogLogistic}} \usage{ SurvivalLogLogistic( - lambda = prior_lognormal(log(0.1), 5, init = 0.1), - p = prior_gamma(2, 5, init = 0.5), + lambda = prior_lognormal(log(0.1), 5), + p = prior_gamma(2, 5), beta = prior_normal(0, 5) ) } diff --git a/man/initialValues.Rd b/man/initialValues.Rd index 65197d40..617a9e7a 100644 --- a/man/initialValues.Rd +++ b/man/initialValues.Rd @@ -1,21 +1,37 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/generics.R, R/Link.R, R/JointModel.R +% Please edit documentation in R/generics.R, R/StanModel.R, R/Link.R, +% R/JointModel.R \name{initialValues} \alias{initialValues} +\alias{initialValues.StanModel} \alias{initialValues.Link} \alias{initialValues.JointModel} \title{\code{initialValues}} \usage{ -initialValues(object) +initialValues(object, ...) -\method{initialValues}{Link}(object) +\method{initialValues}{StanModel}(object, nchains, ...) -\method{initialValues}{JointModel}(object) +\method{initialValues}{Link}(object, nchains, ...) + +\method{initialValues}{JointModel}(object, nchains, ...) } \arguments{ \item{object}{where to get the initial values from.} + +\item{...}{Not currently used.} + +\item{nchains}{the number of initial values to generate. See details.} } \description{ Obtain the \code{list} of initial values to be passed to the Stan sampler. } -\keyword{internal} +\details{ +There are multiple ways of specifying initial values to Stan, see the \code{init} argument +in \link[cmdstanr:model-method-sample]{cmdstanr::model-method-sample} for full details. Within this package we supply +initial values via a list of lists where each inner list contains the initial values +for a single chain. As such the \code{nchains} argument specifies the number of inner lists +to generate. + +See the Vignette for further details of how to specify initial values. +} diff --git a/man/link_gsf_dsld-class.Rd b/man/link_gsf_dsld-class.Rd index aac6c75d..a58df51d 100644 --- a/man/link_gsf_dsld-class.Rd +++ b/man/link_gsf_dsld-class.Rd @@ -7,7 +7,7 @@ \alias{link_gsf_dsld} \title{\code{link_gsf_dsld}} \usage{ -link_gsf_dsld(beta = prior_normal(0, 5, init = 0)) +link_gsf_dsld(beta = prior_normal(0, 5)) } \arguments{ \item{beta}{(\code{Prior})\cr prior for the link coefficient \code{beta}.} diff --git a/man/link_gsf_identity-class.Rd b/man/link_gsf_identity-class.Rd index 28d76eb0..9f586e1d 100644 --- a/man/link_gsf_identity-class.Rd +++ b/man/link_gsf_identity-class.Rd @@ -7,7 +7,7 @@ \alias{link_gsf_identity} \title{\code{link_gsf_identity}} \usage{ -link_gsf_identity(tau = prior_normal(0, 5, init = 0)) +link_gsf_identity(tau = prior_normal(0, 5)) } \arguments{ \item{tau}{(\code{Prior})\cr prior for the link coefficient \code{tau}.} diff --git a/man/link_gsf_ttg-class.Rd b/man/link_gsf_ttg-class.Rd index 96946c7f..1aab6bc6 100644 --- a/man/link_gsf_ttg-class.Rd +++ b/man/link_gsf_ttg-class.Rd @@ -7,7 +7,7 @@ \alias{link_gsf_ttg} \title{\code{link_gsf_ttg}} \usage{ -link_gsf_ttg(gamma = prior_normal(0, 5, init = 0)) +link_gsf_ttg(gamma = prior_normal(0, 5)) } \arguments{ \item{gamma}{(\code{Prior})\cr prior for the link coefficient \code{gamma}.} diff --git a/man/prior_beta.Rd b/man/prior_beta.Rd index 64da7843..f04aaf01 100644 --- a/man/prior_beta.Rd +++ b/man/prior_beta.Rd @@ -4,14 +4,12 @@ \alias{prior_beta} \title{Beta Prior Distribution} \usage{ -prior_beta(a, b, init = a/(a + b)) +prior_beta(a, b) } \arguments{ \item{a}{(\code{number})\cr first parameter.} \item{b}{(\code{number})\cr second parameter} - -\item{init}{(\code{number})\cr initial value.} } \description{ Beta Prior Distribution diff --git a/man/prior_cauchy.Rd b/man/prior_cauchy.Rd index 85614625..f57b71ab 100644 --- a/man/prior_cauchy.Rd +++ b/man/prior_cauchy.Rd @@ -4,14 +4,12 @@ \alias{prior_cauchy} \title{Cauchy Prior Distribution} \usage{ -prior_cauchy(mu, sigma, init = mu) +prior_cauchy(mu, sigma) } \arguments{ \item{mu}{(\code{number})\cr mean.} \item{sigma}{(\code{number})\cr scale.} - -\item{init}{(\code{number})\cr initial value.} } \description{ Cauchy Prior Distribution diff --git a/man/prior_gamma.Rd b/man/prior_gamma.Rd index a1e7667c..a1ddaced 100644 --- a/man/prior_gamma.Rd +++ b/man/prior_gamma.Rd @@ -4,14 +4,12 @@ \alias{prior_gamma} \title{Gamma Prior Distribution} \usage{ -prior_gamma(alpha, beta, init = alpha/beta) +prior_gamma(alpha, beta) } \arguments{ \item{alpha}{(\code{number})\cr shape.} \item{beta}{(\code{number})\cr inverse scale.} - -\item{init}{(\code{number})\cr initial value.} } \description{ Gamma Prior Distribution diff --git a/man/prior_invgamma.Rd b/man/prior_invgamma.Rd index 82fdaaa7..dbc94b9b 100644 --- a/man/prior_invgamma.Rd +++ b/man/prior_invgamma.Rd @@ -4,14 +4,12 @@ \alias{prior_invgamma} \title{Inverse-Gamma Prior Distribution} \usage{ -prior_invgamma(alpha, beta, init = beta/(alpha - 1)) +prior_invgamma(alpha, beta) } \arguments{ \item{alpha}{(\code{number})\cr Shape parameter.} \item{beta}{(\code{number})\cr Scale parameter.} - -\item{init}{(\code{number})\cr initial value.} } \description{ Inverse-Gamma Prior Distribution diff --git a/man/prior_logistic.Rd b/man/prior_logistic.Rd index cdaadae0..74f420d5 100644 --- a/man/prior_logistic.Rd +++ b/man/prior_logistic.Rd @@ -4,14 +4,12 @@ \alias{prior_logistic} \title{Logistic Prior Distribution} \usage{ -prior_logistic(mu, sigma, init = mu) +prior_logistic(mu, sigma) } \arguments{ \item{mu}{(\code{number})\cr Location parameter.} \item{sigma}{(\code{number})\cr Scale parameter.} - -\item{init}{(\code{number})\cr initial value.} } \description{ Logistic Prior Distribution diff --git a/man/prior_loglogistic.Rd b/man/prior_loglogistic.Rd index 89581e94..3f84ff18 100644 --- a/man/prior_loglogistic.Rd +++ b/man/prior_loglogistic.Rd @@ -4,14 +4,12 @@ \alias{prior_loglogistic} \title{Log-Logistic Prior Distribution} \usage{ -prior_loglogistic(alpha, beta, init = alpha * pi/(beta * sin(pi/beta))) +prior_loglogistic(alpha, beta) } \arguments{ \item{alpha}{(\code{number})\cr Scale parameter.} \item{beta}{(\code{number})\cr Shape parameter.} - -\item{init}{(\code{number})\cr initial value.} } \description{ Log-Logistic Prior Distribution diff --git a/man/prior_lognormal.Rd b/man/prior_lognormal.Rd index 9aebbcae..d4140540 100644 --- a/man/prior_lognormal.Rd +++ b/man/prior_lognormal.Rd @@ -4,14 +4,12 @@ \alias{prior_lognormal} \title{Log-Normal Prior Distribution} \usage{ -prior_lognormal(mu, sigma, init = exp(mu + (sigma^2)/2)) +prior_lognormal(mu, sigma) } \arguments{ \item{mu}{(\code{number})\cr mean of the logarithm.} \item{sigma}{(\code{number})\cr standard deviation of the logarithm.} - -\item{init}{(\code{number})\cr initial value.} } \description{ Log-Normal Prior Distribution diff --git a/man/prior_none.Rd b/man/prior_none.Rd index b5e272f0..48787a80 100644 --- a/man/prior_none.Rd +++ b/man/prior_none.Rd @@ -4,10 +4,7 @@ \alias{prior_none} \title{Only Initial Values Specification} \usage{ -prior_none(init = 1e-05) -} -\arguments{ -\item{init}{(\code{number})\cr initial value.} +prior_none() } \description{ Only Initial Values Specification diff --git a/man/prior_normal.Rd b/man/prior_normal.Rd index 6f311e90..1869a014 100644 --- a/man/prior_normal.Rd +++ b/man/prior_normal.Rd @@ -4,14 +4,12 @@ \alias{prior_normal} \title{Normal Prior Distribution} \usage{ -prior_normal(mu, sigma, init = mu) +prior_normal(mu, sigma) } \arguments{ \item{mu}{(\code{number})\cr mean.} \item{sigma}{(\code{number})\cr standard deviation.} - -\item{init}{(\code{number})\cr initial value.} } \description{ Normal Prior Distribution diff --git a/man/prior_std_normal.Rd b/man/prior_std_normal.Rd index 2dcee02d..614e8375 100644 --- a/man/prior_std_normal.Rd +++ b/man/prior_std_normal.Rd @@ -4,10 +4,7 @@ \alias{prior_std_normal} \title{Standard Normal Prior Distribution} \usage{ -prior_std_normal(init = 0) -} -\arguments{ -\item{init}{(\code{number})\cr initial value.} +prior_std_normal() } \description{ Standard Normal Prior Distribution diff --git a/man/prior_student_t.Rd b/man/prior_student_t.Rd index 939d1d94..bfcb6162 100644 --- a/man/prior_student_t.Rd +++ b/man/prior_student_t.Rd @@ -4,7 +4,7 @@ \alias{prior_student_t} \title{Student-t Prior Distribution} \usage{ -prior_student_t(nu, mu, sigma, init = mu) +prior_student_t(nu, mu, sigma) } \arguments{ \item{nu}{(\code{number})\cr Degrees of freedom parameter.} @@ -12,8 +12,6 @@ prior_student_t(nu, mu, sigma, init = mu) \item{mu}{(\code{number})\cr Location parameter.} \item{sigma}{(\code{number})\cr Scale parameter.} - -\item{init}{(\code{number})\cr initial value.} } \description{ Student-t Prior Distribution diff --git a/man/prior_uniform.Rd b/man/prior_uniform.Rd index c9b926c2..af57bfe8 100644 --- a/man/prior_uniform.Rd +++ b/man/prior_uniform.Rd @@ -4,14 +4,12 @@ \alias{prior_uniform} \title{Uniform Prior Distribution} \usage{ -prior_uniform(alpha, beta, init = 0.5 * (alpha + beta)) +prior_uniform(alpha, beta) } \arguments{ \item{alpha}{(\code{number})\cr minimum value parameter.} \item{beta}{(\code{number})\cr maximum value parameter.} - -\item{init}{(\code{number})\cr initial value.} } \description{ Uniform Prior Distribution diff --git a/tests/testthat/_snaps/JointModel.md b/tests/testthat/_snaps/JointModel.md index f2fd02ab..6dbac6bb 100644 --- a/tests/testthat/_snaps/JointModel.md +++ b/tests/testthat/_snaps/JointModel.md @@ -18,8 +18,8 @@ Random Slope Longitudinal Model with parameters: lm_rs_intercept ~ normal(mu = 30, sigma = 10) lm_rs_slope_mu ~ normal(mu = 0, sigma = 15) - lm_rs_slope_sigma ~ lognormal(mu = 1, sigma = 5) - lm_rs_sigma ~ lognormal(mu = 1, sigma = 5) + lm_rs_slope_sigma ~ lognormal(mu = 0, sigma = 1.5) + lm_rs_sigma ~ lognormal(mu = 0, sigma = 1.5) lm_rs_ind_rnd_slope ~ Link: @@ -101,8 +101,8 @@ Random Slope Longitudinal Model with parameters: lm_rs_intercept ~ normal(mu = 30, sigma = 10) lm_rs_slope_mu ~ normal(mu = 0, sigma = 15) - lm_rs_slope_sigma ~ lognormal(mu = 1, sigma = 5) - lm_rs_sigma ~ lognormal(mu = 1, sigma = 5) + 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/LongitudinalRandomSlope.md b/tests/testthat/_snaps/LongitudinalRandomSlope.md index 2b6e34e1..49563b6e 100644 --- a/tests/testthat/_snaps/LongitudinalRandomSlope.md +++ b/tests/testthat/_snaps/LongitudinalRandomSlope.md @@ -8,8 +8,8 @@ Random Slope Longitudinal Model with parameters: lm_rs_intercept ~ normal(mu = 30, sigma = 10) lm_rs_slope_mu ~ normal(mu = 0, sigma = 15) - lm_rs_slope_sigma ~ lognormal(mu = 1, sigma = 5) - lm_rs_sigma ~ lognormal(mu = 1, sigma = 5) + 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,7 +24,7 @@ Random Slope Longitudinal Model with parameters: lm_rs_intercept ~ normal(mu = 0, sigma = 1) lm_rs_slope_mu ~ normal(mu = 0, sigma = 15) - lm_rs_slope_sigma ~ lognormal(mu = 1, sigma = 5) + 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/Prior.md b/tests/testthat/_snaps/Prior.md index b9c26610..656ecbda 100644 --- a/tests/testthat/_snaps/Prior.md +++ b/tests/testthat/_snaps/Prior.md @@ -1,7 +1,7 @@ # show() works for Prior objects Code - print(prior_cauchy(0, 0.8, init = 4)) + print(prior_cauchy(0, 0.8)) Output Prior Object: @@ -81,7 +81,7 @@ --- Code - print(prior_logistic(sigma = 2, init = 1, 10)) + print(prior_logistic(sigma = 2, 10)) Output Prior Object: diff --git a/tests/testthat/_snaps/SurvivalExponential.md b/tests/testthat/_snaps/SurvivalExponential.md index 9ea06394..ba5502f5 100644 --- a/tests/testthat/_snaps/SurvivalExponential.md +++ b/tests/testthat/_snaps/SurvivalExponential.md @@ -13,7 +13,7 @@ --- Code - x <- SurvivalExponential(beta = prior_gamma(3, 4, init = 10)) + x <- SurvivalExponential(beta = prior_gamma(3, 4)) print(x) Output diff --git a/tests/testthat/_snaps/SurvivalLoglogistic.md b/tests/testthat/_snaps/SurvivalLoglogistic.md index d0cd51f6..3f999a12 100644 --- a/tests/testthat/_snaps/SurvivalLoglogistic.md +++ b/tests/testthat/_snaps/SurvivalLoglogistic.md @@ -14,8 +14,7 @@ --- Code - x <- SurvivalLogLogistic(beta = prior_gamma(3, 4, init = 10), p = prior_cauchy( - 0, 1)) + x <- SurvivalLogLogistic(beta = prior_gamma(3, 4), p = prior_cauchy(0, 1)) print(x) Output diff --git a/tests/testthat/_snaps/SurvivalWeibullPH.md b/tests/testthat/_snaps/SurvivalWeibullPH.md index d0920963..a13a8c5d 100644 --- a/tests/testthat/_snaps/SurvivalWeibullPH.md +++ b/tests/testthat/_snaps/SurvivalWeibullPH.md @@ -14,8 +14,7 @@ --- Code - x <- SurvivalWeibullPH(beta = prior_gamma(3, 4, init = 10), gamma = prior_cauchy( - 0, 1)) + x <- SurvivalWeibullPH(beta = prior_gamma(3, 4), gamma = prior_cauchy(0, 1)) print(x) Output diff --git a/tests/testthat/test-LinkRandomSlope.R b/tests/testthat/test-LinkRandomSlope.R index 61e92703..dc578b89 100644 --- a/tests/testthat/test-LinkRandomSlope.R +++ b/tests/testthat/test-LinkRandomSlope.R @@ -1,10 +1,17 @@ test_that("LinkRandomSlope smoke tests", { linkobject <- LinkRandomSlope( - link_lm_phi = prior_normal(0, 5, init = 0.01) + link_lm_phi = prior_normal(0, 5) + ) + iv <- with_mocked_bindings( + initialValues(linkobject, nchains = 2), + local_rnorm = \(...) 5 ) expect_equal( - initialValues(linkobject), - list("link_lm_phi" = 0.01) + iv, + list( + list("link_lm_phi" = 5 / 2), + list("link_lm_phi" = 5 / 2) + ) ) expect_true( is(as.StanModule(linkobject), "StanModule") diff --git a/tests/testthat/test-Parameter.R b/tests/testthat/test-Parameter.R index ad7b0ceb..1637b469 100644 --- a/tests/testthat/test-Parameter.R +++ b/tests/testthat/test-Parameter.R @@ -1,19 +1,24 @@ test_that("Parameters smoke tests", { p <- Parameter(name = "intercept", prior = prior_beta(5, 4)) - expect_equal(initialValues(p), 5 / (5 + 4)) + + expected_mu <- 5 / (5 + 4) + with_mocked_bindings( + expect_equal( + initialValues(p), + expected_mu * 0.5 + ), + local_rbeta = \(...) 0 + ) expect_equal(names(p), "intercept") - p <- Parameter(name = "myp", prior = prior_beta(5, 4, init = 9)) - expect_equal(initialValues(p), 9) - expect_equal(names(p), "myp") }) test_that("show() works for Paramneter objects", { - x <- Parameter(prior_normal(1, 3, 4), "bob", "ben") + x <- Parameter(prior_normal(1, 3), "bob", "size1") expect_snapshot(print(x)) - x <- Parameter(prior_beta(0.5, 0.2), "var1", "var2") + x <- Parameter(prior_beta(0.5, 0.2), "var1", "size1") expect_snapshot(print(x)) - x <- Parameter(prior_none(4), "x", "y") + x <- Parameter(prior_none(), "x", "size1") expect_snapshot(print(x)) }) diff --git a/tests/testthat/test-ParameterList.R b/tests/testthat/test-ParameterList.R index ea2fcf8f..a8ba2c79 100644 --- a/tests/testthat/test-ParameterList.R +++ b/tests/testthat/test-ParameterList.R @@ -5,19 +5,25 @@ test_that("ParameterList smoke tests", { pl <- ParameterList( Parameter(name = "inter", prior = prior_gamma(1, 2)), - Parameter(name = "myp", prior = prior_normal(1, 4, init = 8)) + Parameter(name = "myp", prior = prior_normal(1, 4)) ) # Can extract parameter names expect_equal(names(pl), c("inter", "myp")) # Can extract initial values - actual <- initialValues(pl) - expected <- list( - "inter" = 1 / 2, - "myp" = 8 + with_mocked_bindings( + { + actual <- initialValues(pl, nchains = 1) + expected <- list(list( + "inter" = (1 / 2) * 0.5, + "myp" = 1 * 0.5 + )) + expect_equal(actual, expected) + }, + local_rgamma = \(...) 0, + local_rnorm = \(...) 0 ) - expect_equal(actual, expected) # Can render to character expect_equal( @@ -40,7 +46,7 @@ test_that("show() works for ParameterList objects", { x <- ParameterList( Parameter(name = "bob", prior = prior_normal(1, 4)), Parameter(name = "sam", prior = prior_beta(3, 1)), - Parameter(name = "dave", prior = prior_lognormal(3, 2, init = 3), size = 4), + Parameter(name = "dave", prior = prior_lognormal(3, 2), size = 4), Parameter(name = "steve", prior = prior_none()) ) diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index e2ee6032..4266d441 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -1,18 +1,13 @@ test_that("Priors work as expected", { x <- prior_normal(4, 10) - expect_equal(initialValues(x), 4) - expect_equal( - as.StanModule(x, name = "bob"), - StanModule(test_path("models", "Prior_1.stan")) + with_mocked_bindings( + expect_equal( + initialValues(x), + 4 * 0.5 + ), + local_rnorm = \(...) 0 ) - expect_equal( - as_stan_list(x, name = "bob"), - list(prior_mu_bob = 4, prior_sigma_bob = 10) - ) - - x <- prior_normal(4, 10, 20) - expect_equal(initialValues(x), 20) expect_equal( as.StanModule(x, name = "bob"), StanModule(test_path("models", "Prior_1.stan")) @@ -22,20 +17,14 @@ test_that("Priors work as expected", { list(prior_mu_bob = 4, prior_sigma_bob = 10) ) - x <- prior_lognormal(log(4), 2) - expect_equal(initialValues(x), exp(log(4) + 2)) - expect_equal( - as.StanModule(x, name = "tim"), - StanModule(test_path("models", "Prior_2.stan")) + with_mocked_bindings( + expect_equal( + initialValues(x), + exp(log(4) + 2) * 0.5 + ), + local_rlnorm = \(...) 0 ) - expect_equal( - as_stan_list(x, name = "tim"), - list(prior_mu_tim = log(4), prior_sigma_tim = 2) - ) - - x <- prior_lognormal(log(4), 2, 20) - expect_equal(initialValues(x), 20) expect_equal( as.StanModule(x, name = "tim"), StanModule(test_path("models", "Prior_2.stan")) @@ -46,7 +35,6 @@ test_that("Priors work as expected", { ) - tom <- prior_logistic(1, 2) dave <- prior_loglogistic(3, 4) jim <- prior_invgamma(5, 6) @@ -151,7 +139,7 @@ test_that("Invalid prior parameters are rejected", { test_that("show() works for Prior objects", { - expect_snapshot(print(prior_cauchy(0, 0.8, init = 4))) + expect_snapshot(print(prior_cauchy(0, 0.8))) expect_snapshot(print(prior_normal(0, 0.8))) expect_snapshot(print(prior_std_normal())) expect_snapshot(print(prior_beta(5, 1))) @@ -159,7 +147,7 @@ test_that("show() works for Prior objects", { expect_snapshot(print(prior_none())) expect_snapshot(print(prior_uniform(8, 10))) expect_snapshot(print(prior_student_t(3, 10, 4))) - expect_snapshot(print(prior_logistic(sigma = 2, init = 1, 10))) + expect_snapshot(print(prior_logistic(sigma = 2, 10))) expect_snapshot(print(prior_loglogistic(1, 2))) expect_snapshot(print(prior_invgamma(alpha = 1, beta = 2))) }) diff --git a/tests/testthat/test-StanModel.R b/tests/testthat/test-StanModel.R index ec9b275c..f40cfc42 100644 --- a/tests/testthat/test-StanModel.R +++ b/tests/testthat/test-StanModel.R @@ -4,7 +4,7 @@ test_that("StanModel print function works as expected", { stan = StanModule(), parameters = ParameterList( Parameter(name = "x", prior = prior_normal(3, 1)), - Parameter(name = "z", prior = prior_gamma(3, 1, init = 5)) + Parameter(name = "z", prior = prior_gamma(3, 1)) ), name = "MyModel" ) diff --git a/tests/testthat/test-SurvivalExponential.R b/tests/testthat/test-SurvivalExponential.R index 1721c24d..a89a7751 100644 --- a/tests/testthat/test-SurvivalExponential.R +++ b/tests/testthat/test-SurvivalExponential.R @@ -119,7 +119,7 @@ test_that("Print method for SurvivalExponential works as expected", { }) expect_snapshot({ - x <- SurvivalExponential(beta = prior_gamma(3, 4, init = 10)) + x <- SurvivalExponential(beta = prior_gamma(3, 4)) print(x) }) }) diff --git a/tests/testthat/test-SurvivalLoglogistic.R b/tests/testthat/test-SurvivalLoglogistic.R index 25ff06c5..41ba8159 100644 --- a/tests/testthat/test-SurvivalLoglogistic.R +++ b/tests/testthat/test-SurvivalLoglogistic.R @@ -8,7 +8,7 @@ test_that("Print method for SurvivalLogLogistic works as expected", { expect_snapshot({ x <- SurvivalLogLogistic( - beta = prior_gamma(3, 4, init = 10), + beta = prior_gamma(3, 4), p = prior_cauchy(0, 1) ) print(x) diff --git a/tests/testthat/test-SurvivalWeibullPH.R b/tests/testthat/test-SurvivalWeibullPH.R index 0c50ab6f..da03fc0e 100644 --- a/tests/testthat/test-SurvivalWeibullPH.R +++ b/tests/testthat/test-SurvivalWeibullPH.R @@ -9,7 +9,7 @@ test_that("Print method for SurvivalWeibullPH works as expected", { expect_snapshot({ x <- SurvivalWeibullPH( - beta = prior_gamma(3, 4, init = 10), + beta = prior_gamma(3, 4), gamma = prior_cauchy(0, 1) ) print(x) diff --git a/tests/testthat/test-initialValues.R b/tests/testthat/test-initialValues.R new file mode 100644 index 00000000..c98ffd34 --- /dev/null +++ b/tests/testthat/test-initialValues.R @@ -0,0 +1,83 @@ + + +test_that("initialValues() works as expected", { + jm <- JointModel( + longitudinal = LongitudinalRandomSlope(), + survival = SurvivalWeibullPH(), + link = LinkRandomSlope() + ) + + set.seed(341) + initial_values <- initialValues(jm, nchains = 2) + + # Ensure that we actually got 2 chains worth of initial values + expect_length(initial_values, 2) + + # Ensure that we get different initial values per chain + expect_true( + initial_values[[1]]$lm_rs_intercept != initial_values[[2]]$lm_rs_intercept + ) + + # Ensure each inner list has the same parameters + expect_equal( + names(initial_values[[1]]), + names(initial_values[[2]]) + ) + + # Explicit test to ensure we got all the expected parameters + expect_equal( + c( + "lm_rs_intercept", "lm_rs_slope_mu", "lm_rs_slope_sigma", "lm_rs_sigma", + "lm_rs_ind_rnd_slope", "link_lm_phi", "sm_weibull_ph_lambda", + "sm_weibull_ph_gamma", "beta_os_cov" + ), + names(initial_values[[1]]) + ) + + + # show that if we mock the random number generator, we get the same initial values + ivs <- testthat::with_mocked_bindings( + initialValues(jm, nchains = 2), + local_rnorm = \(...) 0, + local_rbeta = \(...) 0, + local_rlnorm = \(...) 0, + local_rgamma = \(...) 0, + local_runif = \(...) 0, + local_rlogis = \(...) 0, + ) + expect_equal(ivs[[1]], ivs[[2]]) +}) + + + +test_that("ensure_initial_values() works as expected", { + pars <- ParameterList( + Parameter(name = "p1", prior = prior_beta(4, 2), size = 1), + Parameter(name = "p2", prior = prior_normal(5, 1), size = "n_arms"), + Parameter(name = "p3", prior = prior_lognormal(0, 1), size = 3) + ) + + ivs <- initialValues(pars, nchains = 2) + ivs[[1]]$p1 <- c(1) + ivs[[1]]$p2 <- c(2, 3) + ivs[[1]]$p3 <- c(4, 5, 6) + + dta <- list(n_arms = 2) + + res <- ensure_initial_values(ivs, dta, pars) + + # First check that the lengths have been expanded out to be the expected + # sizes + expect_length(res[[1]]$p1, 1) + expect_length(res[[1]]$p2, 2) + expect_length(res[[1]]$p3, 3) + expect_length(res[[2]]$p1, 1) + expect_length(res[[2]]$p2, 2) + expect_length(res[[2]]$p3, 3) + + # Now double check that our pre-set values in the first element + # have not been modified + expect_equal(res[[1]]$p1, 1) + expect_equal(res[[1]]$p2, array(c(2, 3), dim = c(2))) + expect_equal(res[[1]]$p3, array(c(4, 5, 6), dim = c(3))) +}) diff --git a/vignettes/model_fitting.Rmd b/vignettes/model_fitting.Rmd index 62e6716f..ec5b43d8 100644 --- a/vignettes/model_fitting.Rmd +++ b/vignettes/model_fitting.Rmd @@ -373,3 +373,45 @@ sq <- SurvivalQuantities( ) brierScore(sq) ``` + + +## Initial Values + +By default `jmpost` will set the initial values for all parameters to be a random value +drawn from the prior distribution that has been shrunk towards the +mean of said distribution e.g. for a `prior_normal(4, 2)` +the initial value for each chain will be: `(4 + rnorm(1, 4, 2)) * 0.5`. + +If you wish to manually specify the initial values you can do so via the `initialValues()` function. +For example: + +```{r, eval=FALSE} +joint_model <- JointModel( + LongitudinalRandomSlope(), + SurvivalExponential(), + LinkNone() +) +initial_values <- initialValues(joint_model, nchains = 2) + +initial_values[[1]]$lm_rs_intercept <- 0.2 +initial_values[[2]]$lm_rs_intercept <- 0.3 + +mcmc_results <- sampleStanModel( + joint_model, + data = DataJoint(...), + init = initial_values, + chains = 2 +) +``` + +Note the following: +- `initialValues()` will return a list of lists where each sublist contains the initial +values for the corresponding chain index. +- `initialValues()` by default just returns 1 initial value for each parameter; the `sampleStanModel()` +function will then broadcast this number as many times as required e.g. if you have 3 covariates +in your survival model then the initial value for the Beta coeficient will be repeated 3 times. If +however you want to specify individual initial values for each covariate you can do so by passing +in a vector of the same length as the number of covariates. + + + From f29172ad1ab09d47434d5a724f32798bd6e44e61 Mon Sep 17 00:00:00 2001 From: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Date: Fri, 2 Feb 2024 17:06:07 +0000 Subject: [PATCH 02/21] [skip actions] Roxygen Man Pages Auto Update --- man/ensure_initial_values.Rd | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 man/ensure_initial_values.Rd diff --git a/man/ensure_initial_values.Rd b/man/ensure_initial_values.Rd new file mode 100644 index 00000000..d9bfc623 --- /dev/null +++ b/man/ensure_initial_values.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/JointModel.R +\name{ensure_initial_values} +\alias{ensure_initial_values} +\title{Ensure that initial values are correctly specified} +\usage{ +ensure_initial_values(initial_values, data, parameters) +} +\arguments{ +\item{initial_values}{(\code{list})\cr A list of lists containing the initial values +must be 1 list per desired chain. All elements should have identical names} + +\item{data}{(\code{list})\cr specifies the size to expand each of our initial values to be. +That is elements of size 1 in \code{initial_values} will be expanded to be the same +size as the corresponding element in \code{data} by broadcasting the value.} + +\item{parameters}{(\code{\link{ParameterList}})\cr the parameters object} +} +\description{ +Ensure that initial values are correctly specified +} +\details{ +This function is mostly a thin wrapper around \code{expand_initial_values} to +enable easier unit testing. +} +\keyword{internal} From 18c9720d433c50145017393ab54f3633694e9613 Mon Sep 17 00:00:00 2001 From: gowerc Date: Fri, 2 Feb 2024 17:13:26 +0000 Subject: [PATCH 03/21] Update to fix broken CI tests --- R/LongitudinalRandomSlope.R | 13 ++----------- _pkgdown.yml | 1 + man/LongitudinalRandomSlope-class.Rd | 5 +---- man/ensure_initial_values.Rd | 26 ++++++++++++++++++++++++++ vignettes/model_fitting.Rmd | 17 +++++++---------- 5 files changed, 37 insertions(+), 25 deletions(-) create mode 100644 man/ensure_initial_values.Rd diff --git a/R/LongitudinalRandomSlope.R b/R/LongitudinalRandomSlope.R index a2fafc0d..eb5a866f 100755 --- a/R/LongitudinalRandomSlope.R +++ b/R/LongitudinalRandomSlope.R @@ -24,28 +24,19 @@ NULL #' @param slope_mu (`Prior`)\cr for the population slope `slope_mu`. #' @param slope_sigma (`Prior`)\cr for the random slope standard deviation `slope_sigma`. #' @param sigma (`Prior`)\cr for the variance of the longitudinal values `sigma`. -#' @param random_slope (`Prior`)\cr must be `prior_none()`, just used to specify initial values. #' #' @export LongitudinalRandomSlope <- function( intercept = prior_normal(30, 10), slope_mu = prior_normal(0, 15), slope_sigma = prior_lognormal(0, 1.5), - sigma = prior_lognormal(0, 1.5), - random_slope = prior_none() + sigma = prior_lognormal(0, 1.5) ) { stan <- StanModule( x = "lm-random-slope/model.stan" ) - assert_that( - inherits(random_slope, "Prior"), - random_slope@repr_data == "", - random_slope@repr_model == "", - msg = "`random_slope` must be a `prior_none()`" - ) - .LongitudinalRandomSlope( LongitudinalModel( name = "Random Slope", @@ -55,7 +46,7 @@ LongitudinalRandomSlope <- function( Parameter(name = "lm_rs_slope_mu", prior = slope_mu, size = "n_arms"), Parameter(name = "lm_rs_slope_sigma", prior = slope_sigma, size = 1), Parameter(name = "lm_rs_sigma", prior = sigma, size = 1), - Parameter(name = "lm_rs_ind_rnd_slope", prior = random_slope, size = "Nind") + Parameter(name = "lm_rs_ind_rnd_slope", prior = prior_none(), size = "Nind") ) ) ) diff --git a/_pkgdown.yml b/_pkgdown.yml index e626c09c..9ab3ae93 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -122,6 +122,7 @@ reference: - title: Convenience Functions contents: + - initialValues - merge - show - as_stan_list diff --git a/man/LongitudinalRandomSlope-class.Rd b/man/LongitudinalRandomSlope-class.Rd index cb99a5c9..38d093d4 100644 --- a/man/LongitudinalRandomSlope-class.Rd +++ b/man/LongitudinalRandomSlope-class.Rd @@ -11,8 +11,7 @@ LongitudinalRandomSlope( intercept = prior_normal(30, 10), slope_mu = prior_normal(0, 15), slope_sigma = prior_lognormal(0, 1.5), - sigma = prior_lognormal(0, 1.5), - random_slope = prior_none() + sigma = prior_lognormal(0, 1.5) ) } \arguments{ @@ -23,8 +22,6 @@ LongitudinalRandomSlope( \item{slope_sigma}{(\code{Prior})\cr for the random slope standard deviation \code{slope_sigma}.} \item{sigma}{(\code{Prior})\cr for the variance of the longitudinal values \code{sigma}.} - -\item{random_slope}{(\code{Prior})\cr must be \code{prior_none()}, just used to specify initial values.} } \description{ This class extends the general \code{\link{LongitudinalModel}} class for using the diff --git a/man/ensure_initial_values.Rd b/man/ensure_initial_values.Rd new file mode 100644 index 00000000..d9bfc623 --- /dev/null +++ b/man/ensure_initial_values.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/JointModel.R +\name{ensure_initial_values} +\alias{ensure_initial_values} +\title{Ensure that initial values are correctly specified} +\usage{ +ensure_initial_values(initial_values, data, parameters) +} +\arguments{ +\item{initial_values}{(\code{list})\cr A list of lists containing the initial values +must be 1 list per desired chain. All elements should have identical names} + +\item{data}{(\code{list})\cr specifies the size to expand each of our initial values to be. +That is elements of size 1 in \code{initial_values} will be expanded to be the same +size as the corresponding element in \code{data} by broadcasting the value.} + +\item{parameters}{(\code{\link{ParameterList}})\cr the parameters object} +} +\description{ +Ensure that initial values are correctly specified +} +\details{ +This function is mostly a thin wrapper around \code{expand_initial_values} to +enable easier unit testing. +} +\keyword{internal} diff --git a/vignettes/model_fitting.Rmd b/vignettes/model_fitting.Rmd index ec5b43d8..4f639e19 100644 --- a/vignettes/model_fitting.Rmd +++ b/vignettes/model_fitting.Rmd @@ -56,12 +56,12 @@ from the arguments of the constructors (or from the help page): args(LongitudinalRandomSlope) ``` -So here we see that there are 4 model parameters with a proper prior and an -additional `random_slope` one where just the initial value can be specified. +So here we see that the Longitudinal Random Slope model has 4 parameters that we can +define a prior for. -## Specifying Prior and Initial Values +## Specifying Priors -We can alternatively also specify the prior distributions and the initial values +We can alternatively also specify the prior distributions for the parameters manually. This is important in practice to obtain a meaningful model specification and hence converging MCMC chains that allow to estimate the posterior distributions. @@ -71,15 +71,12 @@ For the random slope model for the longitudinal outcome, we can e.g. say: ```{r random_slope_prior_inits} random_slope_model <- LongitudinalRandomSlope( intercept = prior_normal(40, 5), - slope_mu = prior_normal(10, 2, init = 30) + slope_mu = prior_normal(10, 2) ) ``` -This sets the prior for the `intercept`, without stating an explicit initial value. -In this case the `init` value is just set to the mean of the normal distribution, -so 40. The `slope_mu` argument then also sets the prior and a custom initial value -for the `slope_mu` parameter. The other parameters are left with their default priors -and initial values. +This sets the prior for the `intercept` to be a $N(40, 5)$ distribution and +the prior for the `slope_mu` parameter to be a $N(10, 2)$ distribution. ## Separate Models From 114aeeadf8f384a54f0c7a3bb03ed1ed529c1120 Mon Sep 17 00:00:00 2001 From: gowerc Date: Fri, 2 Feb 2024 17:18:36 +0000 Subject: [PATCH 04/21] fix spellings --- inst/WORDLIST | 1 + 1 file changed, 1 insertion(+) diff --git a/inst/WORDLIST b/inst/WORDLIST index ccdcf4df..c9a5df06 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -115,3 +115,4 @@ gk xk LogNormal AST +throughs From 78cf8074b2dfd21214c520a3de3415419a28a5a7 Mon Sep 17 00:00:00 2001 From: gowerc Date: Fri, 2 Feb 2024 17:19:29 +0000 Subject: [PATCH 05/21] fix lintr error --- R/generics.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/generics.R b/R/generics.R index 3bfd281e..ce0bc2ba 100755 --- a/R/generics.R +++ b/R/generics.R @@ -156,7 +156,7 @@ extractVariableNames <- function(object) { #' #' @details #' There are multiple ways of specifying initial values to Stan, see the `init` argument -#' in [cmdstanr::model-method-sample] for full details. Within this package we supply +#' in [cmdstanr::model-method-sample] for full details. Within this package we supply #' initial values via a list of lists where each inner list contains the initial values #' for a single chain. As such the `nchains` argument specifies the number of inner lists #' to generate. From d13a0905348c3f1b7ed425a57f55af63beee750f Mon Sep 17 00:00:00 2001 From: gowerc Date: Fri, 2 Feb 2024 17:39:27 +0000 Subject: [PATCH 06/21] Fix documentation --- NAMESPACE | 8 ++++++++ R/Parameter.R | 1 + R/ParameterList.R | 1 + R/Prior.R | 2 ++ man/Parameter-Getter-Methods.Rd | 2 ++ man/ParameterList-Getter-Methods.Rd | 2 ++ 6 files changed, 16 insertions(+) diff --git a/NAMESPACE b/NAMESPACE index 799c4438..bb6e1912 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -165,5 +165,13 @@ importFrom(ggplot2,autoplot) importFrom(ggplot2.utils,geom_km) importFrom(glue,as_glue) importFrom(stats,acf) +importFrom(stats,rbeta) +importFrom(stats,rcauchy) +importFrom(stats,rgamma) +importFrom(stats,rlnorm) +importFrom(stats,rlogis) +importFrom(stats,rnorm) +importFrom(stats,rt) +importFrom(stats,runif) importFrom(survival,Surv) importFrom(tibble,add_case) diff --git a/R/Parameter.R b/R/Parameter.R index 8bdf5bac..cec2251a 100755 --- a/R/Parameter.R +++ b/R/Parameter.R @@ -102,6 +102,7 @@ as_stan_list.Parameter <- function(object, ...) { #' #' @param x (`Paramater`) \cr A model parameter #' @param object (`Paramater`) \cr A model parameter +#' @param ... Not Used. #' #' @description #' Getter functions for the slots of a [`Parameter`] object diff --git a/R/ParameterList.R b/R/ParameterList.R index 16fb4961..a3e4e471 100644 --- a/R/ParameterList.R +++ b/R/ParameterList.R @@ -132,6 +132,7 @@ as.list.ParameterList <- function(x, ...) { #' Getter functions for the slots of a [`ParameterList`] object #' @inheritParams ParameterList-Shared #' @family ParameterList +#' @param nchains (`integer`) \cr the number of chains. #' @name ParameterList-Getter-Methods NULL diff --git a/R/Prior.R b/R/Prior.R index 7abac649..1a7d4300 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -523,6 +523,8 @@ prior_invgamma <- function(alpha, beta) { #' @param nu (`number`)\cr Parameter for underlying distribution. #' @param ... Pass any additional arguments to the underlying distribution. #' +#' @importFrom stats rbeta rcauchy rgamma rlnorm rlogis rnorm rt runif +#' #' @name Local_Sample #' @keywords internal NULL diff --git a/man/Parameter-Getter-Methods.Rd b/man/Parameter-Getter-Methods.Rd index f9b65e2a..01a1db30 100644 --- a/man/Parameter-Getter-Methods.Rd +++ b/man/Parameter-Getter-Methods.Rd @@ -17,6 +17,8 @@ \item{x}{(\code{Paramater}) \cr A model parameter} \item{object}{(\code{Paramater}) \cr A model parameter} + +\item{...}{Not Used.} } \description{ Getter functions for the slots of a \code{\link{Parameter}} object diff --git a/man/ParameterList-Getter-Methods.Rd b/man/ParameterList-Getter-Methods.Rd index 4e85c466..3d2771f9 100644 --- a/man/ParameterList-Getter-Methods.Rd +++ b/man/ParameterList-Getter-Methods.Rd @@ -18,6 +18,8 @@ \item{object}{(\code{ParameterList}) \cr A List of \code{\link{Parameter}} Objects.} +\item{nchains}{(\code{integer}) \cr the number of chains.} + \item{...}{Not Used.} } \description{ From 3381560c35a5ff12099089d6121fe99f3b196028 Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 16:07:34 +0000 Subject: [PATCH 07/21] Switch to nchains to n_chains --- R/JointModel.R | 6 +++--- R/Link.R | 4 ++-- R/ParameterList.R | 6 +++--- R/StanModel.R | 4 ++-- R/generics.R | 4 ++-- man/ParameterList-Getter-Methods.Rd | 4 ++-- man/initialValues.Rd | 10 +++++----- tests/testthat/test-LinkRandomSlope.R | 2 +- tests/testthat/test-ParameterList.R | 2 +- tests/testthat/test-initialValues.R | 6 +++--- 10 files changed, 24 insertions(+), 24 deletions(-) diff --git a/R/JointModel.R b/R/JointModel.R index c09816e6..8a061e1c 100755 --- a/R/JointModel.R +++ b/R/JointModel.R @@ -167,7 +167,7 @@ sampleStanModel.JointModel <- function(object, data, ...) { initial_values <- if ("init" %in% names(args)) { args[["init"]] } else { - initialValues(object, nchains = args[["chains"]]) + initialValues(object, n_chains = args[["chains"]]) } args[["init"]] <- ensure_initial_values( @@ -230,8 +230,8 @@ ensure_initial_values <- function(initial_values, data, parameters) { #' @rdname initialValues #' @export -initialValues.JointModel <- function(object, nchains, ...) { - initialValues(object@parameters, nchains) +initialValues.JointModel <- function(object, n_chains, ...) { + initialValues(object@parameters, n_chains) } diff --git a/R/Link.R b/R/Link.R index ea845168..e85ecd4c 100755 --- a/R/Link.R +++ b/R/Link.R @@ -60,8 +60,8 @@ setMethod( #' @rdname initialValues #' @export -initialValues.Link <- function(object, nchains, ...) { - initialValues(object@parameters, nchains) +initialValues.Link <- function(object, n_chains, ...) { + initialValues(object@parameters, n_chains) } diff --git a/R/ParameterList.R b/R/ParameterList.R index a3e4e471..14b9e4ca 100644 --- a/R/ParameterList.R +++ b/R/ParameterList.R @@ -132,7 +132,7 @@ as.list.ParameterList <- function(x, ...) { #' Getter functions for the slots of a [`ParameterList`] object #' @inheritParams ParameterList-Shared #' @family ParameterList -#' @param nchains (`integer`) \cr the number of chains. +#' @param n_chains (`integer`) \cr the number of chains. #' @name ParameterList-Getter-Methods NULL @@ -146,9 +146,9 @@ names.ParameterList <- function(x) { #' @describeIn ParameterList-Getter-Methods The parameter-list's parameter initial values #' @export -initialValues.ParameterList <- function(object, nchains, ...) { +initialValues.ParameterList <- function(object, n_chains, ...) { x <- lapply( - seq_len(nchains), + seq_len(n_chains), \(i) { vals <- lapply(object@parameters, initialValues) name <- vapply(object@parameters, names, character(1)) diff --git a/R/StanModel.R b/R/StanModel.R index 703d59c3..f05e07af 100644 --- a/R/StanModel.R +++ b/R/StanModel.R @@ -93,6 +93,6 @@ setMethod( #' @rdname initialValues #' @export -initialValues.StanModel <- function(object, nchains, ...) { - initialValues(object@parameters, nchains) +initialValues.StanModel <- function(object, n_chains, ...) { + initialValues(object@parameters, n_chains) } diff --git a/R/generics.R b/R/generics.R index ce0bc2ba..dcd640a2 100755 --- a/R/generics.R +++ b/R/generics.R @@ -151,14 +151,14 @@ extractVariableNames <- function(object) { #' Obtain the `list` of initial values to be passed to the Stan sampler. #' #' @param object where to get the initial values from. -#' @param nchains the number of initial values to generate. See details. +#' @param n_chains the number of initial values to generate. See details. #' @param ... Not currently used. #' #' @details #' There are multiple ways of specifying initial values to Stan, see the `init` argument #' in [cmdstanr::model-method-sample] for full details. Within this package we supply #' initial values via a list of lists where each inner list contains the initial values -#' for a single chain. As such the `nchains` argument specifies the number of inner lists +#' for a single chain. As such the `n_chains` argument specifies the number of inner lists #' to generate. #' #' See the Vignette for further details of how to specify initial values. diff --git a/man/ParameterList-Getter-Methods.Rd b/man/ParameterList-Getter-Methods.Rd index 3d2771f9..a922f0d2 100644 --- a/man/ParameterList-Getter-Methods.Rd +++ b/man/ParameterList-Getter-Methods.Rd @@ -9,7 +9,7 @@ \usage{ \method{names}{ParameterList}(x) -\method{initialValues}{ParameterList}(object, nchains, ...) +\method{initialValues}{ParameterList}(object, n_chains, ...) \method{size}{ParameterList}(object) } @@ -18,7 +18,7 @@ \item{object}{(\code{ParameterList}) \cr A List of \code{\link{Parameter}} Objects.} -\item{nchains}{(\code{integer}) \cr the number of chains.} +\item{n_chains}{(\code{integer}) \cr the number of chains.} \item{...}{Not Used.} } diff --git a/man/initialValues.Rd b/man/initialValues.Rd index 617a9e7a..25536f89 100644 --- a/man/initialValues.Rd +++ b/man/initialValues.Rd @@ -10,18 +10,18 @@ \usage{ initialValues(object, ...) -\method{initialValues}{StanModel}(object, nchains, ...) +\method{initialValues}{StanModel}(object, n_chains, ...) -\method{initialValues}{Link}(object, nchains, ...) +\method{initialValues}{Link}(object, n_chains, ...) -\method{initialValues}{JointModel}(object, nchains, ...) +\method{initialValues}{JointModel}(object, n_chains, ...) } \arguments{ \item{object}{where to get the initial values from.} \item{...}{Not currently used.} -\item{nchains}{the number of initial values to generate. See details.} +\item{n_chains}{the number of initial values to generate. See details.} } \description{ Obtain the \code{list} of initial values to be passed to the Stan sampler. @@ -30,7 +30,7 @@ Obtain the \code{list} of initial values to be passed to the Stan sampler. There are multiple ways of specifying initial values to Stan, see the \code{init} argument in \link[cmdstanr:model-method-sample]{cmdstanr::model-method-sample} for full details. Within this package we supply initial values via a list of lists where each inner list contains the initial values -for a single chain. As such the \code{nchains} argument specifies the number of inner lists +for a single chain. As such the \code{n_chains} argument specifies the number of inner lists to generate. See the Vignette for further details of how to specify initial values. diff --git a/tests/testthat/test-LinkRandomSlope.R b/tests/testthat/test-LinkRandomSlope.R index dc578b89..55fda0f3 100644 --- a/tests/testthat/test-LinkRandomSlope.R +++ b/tests/testthat/test-LinkRandomSlope.R @@ -3,7 +3,7 @@ test_that("LinkRandomSlope smoke tests", { link_lm_phi = prior_normal(0, 5) ) iv <- with_mocked_bindings( - initialValues(linkobject, nchains = 2), + initialValues(linkobject, n_chains = 2), local_rnorm = \(...) 5 ) expect_equal( diff --git a/tests/testthat/test-ParameterList.R b/tests/testthat/test-ParameterList.R index a8ba2c79..eec6abab 100644 --- a/tests/testthat/test-ParameterList.R +++ b/tests/testthat/test-ParameterList.R @@ -14,7 +14,7 @@ test_that("ParameterList smoke tests", { # Can extract initial values with_mocked_bindings( { - actual <- initialValues(pl, nchains = 1) + actual <- initialValues(pl, n_chains = 1) expected <- list(list( "inter" = (1 / 2) * 0.5, "myp" = 1 * 0.5 diff --git a/tests/testthat/test-initialValues.R b/tests/testthat/test-initialValues.R index c98ffd34..2a52918f 100644 --- a/tests/testthat/test-initialValues.R +++ b/tests/testthat/test-initialValues.R @@ -8,7 +8,7 @@ test_that("initialValues() works as expected", { ) set.seed(341) - initial_values <- initialValues(jm, nchains = 2) + initial_values <- initialValues(jm, n_chains = 2) # Ensure that we actually got 2 chains worth of initial values expect_length(initial_values, 2) @@ -37,7 +37,7 @@ test_that("initialValues() works as expected", { # show that if we mock the random number generator, we get the same initial values ivs <- testthat::with_mocked_bindings( - initialValues(jm, nchains = 2), + initialValues(jm, n_chains = 2), local_rnorm = \(...) 0, local_rbeta = \(...) 0, local_rlnorm = \(...) 0, @@ -57,7 +57,7 @@ test_that("ensure_initial_values() works as expected", { Parameter(name = "p3", prior = prior_lognormal(0, 1), size = 3) ) - ivs <- initialValues(pars, nchains = 2) + ivs <- initialValues(pars, n_chains = 2) ivs[[1]]$p1 <- c(1) ivs[[1]]$p2 <- c(2, 3) ivs[[1]]$p3 <- c(4, 5, 6) From 031ef752876a29615280565ea0925a3e8112f621 Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 16:19:02 +0000 Subject: [PATCH 08/21] Resolve Typo --- R/Parameter.R | 2 +- man/Parameter-Getter-Methods.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/Parameter.R b/R/Parameter.R index cec2251a..f9fa0ca7 100755 --- a/R/Parameter.R +++ b/R/Parameter.R @@ -102,7 +102,7 @@ as_stan_list.Parameter <- function(object, ...) { #' #' @param x (`Paramater`) \cr A model parameter #' @param object (`Paramater`) \cr A model parameter -#' @param ... Not Used. +#' @param ... Not used. #' #' @description #' Getter functions for the slots of a [`Parameter`] object diff --git a/man/Parameter-Getter-Methods.Rd b/man/Parameter-Getter-Methods.Rd index 01a1db30..06d2b76e 100644 --- a/man/Parameter-Getter-Methods.Rd +++ b/man/Parameter-Getter-Methods.Rd @@ -18,7 +18,7 @@ \item{object}{(\code{Paramater}) \cr A model parameter} -\item{...}{Not Used.} +\item{...}{Not used.} } \description{ Getter functions for the slots of a \code{\link{Parameter}} object From c587a69f837f4efdd74eb1f1f4066ec0ef712b50 Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 16:21:14 +0000 Subject: [PATCH 09/21] convert lapply to vapply --- R/ParameterList.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/ParameterList.R b/R/ParameterList.R index 14b9e4ca..f2a5aca9 100644 --- a/R/ParameterList.R +++ b/R/ParameterList.R @@ -147,16 +147,16 @@ names.ParameterList <- function(x) { #' @describeIn ParameterList-Getter-Methods The parameter-list's parameter initial values #' @export initialValues.ParameterList <- function(object, n_chains, ...) { - x <- lapply( + vapply( seq_len(n_chains), \(i) { vals <- lapply(object@parameters, initialValues) name <- vapply(object@parameters, names, character(1)) names(vals) <- name vals - } + }, + FUN.VALUE = list() ) - return(x) } From 1b6ffe5029b0cb9ee65fc5685e62bf11b8ad7f6c Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 16:26:06 +0000 Subject: [PATCH 10/21] Ensure alpha < beta for uniform prior --- R/Prior.R | 4 ++++ tests/testthat/test-Prior.R | 5 +++++ 2 files changed, 9 insertions(+) diff --git a/R/Prior.R b/R/Prior.R index 1a7d4300..8ef9f09a 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -364,6 +364,10 @@ prior_none <- function() { #' #' @export prior_uniform <- function(alpha, beta) { + assert_that( + alpha < beta, + msg = "`alpha`` must be less than `beta`" + ) Prior( parameters = list(alpha = alpha, beta = beta), display = "uniform(alpha = {alpha}, beta = {beta})", diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 4266d441..e5536f34 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -130,6 +130,11 @@ test_that("Invalid prior parameters are rejected", { regexp = "Invalid.*`beta`" ) + expect_error( + prior_uniform(10, 9), + regexp = "`alpha`` must be less than `beta`" + ) + # Ensure that validation doesn't wrongly reject priors with no user specified parameters expect_s4_class(prior_none(), "Prior") From fba92f82e3a66a735d9f8e2fa5236bc095e9fc1d Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 16:58:52 +0000 Subject: [PATCH 11/21] reverted lappy code --- R/ParameterList.R | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/R/ParameterList.R b/R/ParameterList.R index f2a5aca9..1c435d7c 100644 --- a/R/ParameterList.R +++ b/R/ParameterList.R @@ -147,15 +147,17 @@ names.ParameterList <- function(x) { #' @describeIn ParameterList-Getter-Methods The parameter-list's parameter initial values #' @export initialValues.ParameterList <- function(object, n_chains, ...) { - vapply( + # 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 - }, - FUN.VALUE = list() + } ) } From af2518ab9d4d27beefa58fc6ca76b8d90cbf8e1f Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 17:31:24 +0000 Subject: [PATCH 12/21] Added option to change shrinkage factor --- DESCRIPTION | 1 + R/Prior.R | 3 ++- R/settings.R | 35 +++++++++++++++++++++++++++++++++++ man/jmpost-settings.Rd | 31 +++++++++++++++++++++++++++++++ tests/testthat/test-Prior.R | 30 ++++++++++++++++++++++++++++++ 5 files changed, 99 insertions(+), 1 deletion(-) create mode 100644 R/settings.R create mode 100644 man/jmpost-settings.Rd diff --git a/DESCRIPTION b/DESCRIPTION index cff59d90..37e6ee08 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -88,6 +88,7 @@ Collate: 'defaults.R' 'external-exports.R' 'jmpost-package.R' + 'settings.R' 'simulations.R' 'simulations_gsf.R' 'simulations_os.R' diff --git a/R/Prior.R b/R/Prior.R index 8ef9f09a..0d00d186 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -188,7 +188,8 @@ NULL #' @describeIn Prior-Getter-Methods The prior's initial value #' @export initialValues.Prior <- function(object, ...) { - 0.5 * object@init + 0.5 * object@sample(1) + getOption("jmpost.prior_shrinkage") * object@init + + (1 - getOption("jmpost.prior_shrinkage")) * object@sample(1) } diff --git a/R/settings.R b/R/settings.R new file mode 100644 index 00000000..099acb9f --- /dev/null +++ b/R/settings.R @@ -0,0 +1,35 @@ + + + +#' jmpost settings +#' +#' @description +#' Define settings that modify the behaviour of the `jmpost` package +#' +#' Each of the following are the name of options that can be set via: +#' ``` +#' options( = ) +#' ``` +#' +#' ## `jmpost.prior_shrinkage` +#' +#' Default = `0.5` +#' +#' By default all initial values are drawn as random sample from the respective prior +#' distribution with a shrinkage factor towards the mean. That is: +#' ``` +#' initial_value = prior_sample * prior_shrinkage + (1 - prior_shrinkage) * prior_mean +#' ``` +#' This setting controls the shrinkage factor. A value of 0 means no shrinkage (i.e. +#' pure random draw) whilst a value of 1 means the intial value is just the mean. +#' +#' @examples +#' \dontrun{ +#' options(jmpost.prior_shrinkage = 0.5) +#' } +#' @name jmpost-settings +NULL + +if (is.null(getOption("jmpost.prior_shrinkage", default = NULL))) { + options("jmpost.prior_shrinkage" = 0.5) +} diff --git a/man/jmpost-settings.Rd b/man/jmpost-settings.Rd new file mode 100644 index 00000000..19c31663 --- /dev/null +++ b/man/jmpost-settings.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/settings.R +\name{jmpost-settings} +\alias{jmpost-settings} +\title{jmpost settings} +\description{ +Define settings that modify the behaviour of the \code{jmpost} package + +Each of the following are the name of options that can be set via: + +\if{html}{\out{
}}\preformatted{options( = ) +}\if{html}{\out{
}} +\subsection{\code{jmpost.prior_shrinkage}}{ + +Default = \code{0.5} + +By default all initial values are drawn as random sample from the respective prior +distribution with a shrinkage factor towards the mean. That is: + +\if{html}{\out{
}}\preformatted{initial_value = prior_sample * prior_shrinkage + (1 - prior_shrinkage) * prior_mean +}\if{html}{\out{
}} + +This setting controls the shrinkage factor. A value of 0 means no shrinkage (i.e. +pure random draw) whilst a value of 1 means the intial value is just the mean. +} +} +\examples{ +\dontrun{ +options(jmpost.prior_shrinkage = 0.5) +} +} diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index e5536f34..859f9223 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -156,3 +156,33 @@ test_that("show() works for Prior objects", { expect_snapshot(print(prior_loglogistic(1, 2))) expect_snapshot(print(prior_invgamma(alpha = 1, beta = 2))) }) + + +test_that("jmpost.prior_shrinkage works as expected", { + x <- prior_normal(1, 2) + with_mocked_bindings( + { + options("jmpost.prior_shrinkage" = 0.5) + expect_equal( + initialValues(x), + 1 * 0.5 + 4 * 0.5 + ) + + options("jmpost.prior_shrinkage" = 0.9) + expect_equal( + initialValues(x), + 1 * 0.9 + 4 * 0.1 + ) + + options("jmpost.prior_shrinkage" = 0.1) + expect_equal( + initialValues(x), + 1 * 0.1 + 4 * 0.9 + ) + + ## Reset Shrinkage factor + options("jmpost.prior_shrinkage" = 0.1) + }, + local_rnorm = \(...) 4 + ) +}) From 38ee15f57ee461c9c5a9f76d9da801707b5a882a Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 17:35:18 +0000 Subject: [PATCH 13/21] add shrinkage information to vignette --- tests/testthat/test-Prior.R | 2 +- vignettes/model_fitting.Rmd | 12 +++++++++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 859f9223..1defd442 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -181,7 +181,7 @@ test_that("jmpost.prior_shrinkage works as expected", { ) ## Reset Shrinkage factor - options("jmpost.prior_shrinkage" = 0.1) + options("jmpost.prior_shrinkage" = 0.5) }, local_rnorm = \(...) 4 ) diff --git a/vignettes/model_fitting.Rmd b/vignettes/model_fitting.Rmd index 4f639e19..5ea6db0c 100644 --- a/vignettes/model_fitting.Rmd +++ b/vignettes/model_fitting.Rmd @@ -377,7 +377,17 @@ brierScore(sq) By default `jmpost` will set the initial values for all parameters to be a random value drawn from the prior distribution that has been shrunk towards the mean of said distribution e.g. for a `prior_normal(4, 2)` -the initial value for each chain will be: `(4 + rnorm(1, 4, 2)) * 0.5`. +the initial value for each chain will be: +``` +4 * shrinkage_factor + rnorm(1, 4, 2) * (1 - shrinkage_factor) +``` + +Note that the shrinkage factor is set to 0.5 by default and can be changed via the +`jmpost.prior_shrinkage` option e.g. + +``` +options("jmpost.prior_shrinkage" = 0.7) +``` If you wish to manually specify the initial values you can do so via the `initialValues()` function. For example: From 730ed0adee48828f4c53910c4490f0d9570f4142 Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 17:49:38 +0000 Subject: [PATCH 14/21] formally documented magic number --- DESCRIPTION | 1 + R/JointModel.R | 4 +++- R/constants.R | 14 ++++++++++++++ 3 files changed, 18 insertions(+), 1 deletion(-) create mode 100644 R/constants.R diff --git a/DESCRIPTION b/DESCRIPTION index 37e6ee08..25496c02 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -63,6 +63,7 @@ Collate: 'DataLongitudinal.R' 'DataSurvival.R' 'DataJoint.R' + 'constants.R' 'StanModule.R' 'Prior.R' 'Parameter.R' diff --git a/R/JointModel.R b/R/JointModel.R index 8a061e1c..61386db6 100755 --- a/R/JointModel.R +++ b/R/JointModel.R @@ -8,6 +8,7 @@ #' @include LongitudinalModel.R #' @include SurvivalModel.R #' @include Link.R +#' @include constants.R NULL @@ -161,7 +162,8 @@ sampleStanModel.JointModel <- function(object, data, ...) { args[["chains"]] <- if ("chains" %in% names(args)) { args[["chains"]] } else { - 4 + # Magic constant from R/constants.R + CMDSTAN_DEFAULT_CHAINS } initial_values <- if ("init" %in% names(args)) { diff --git a/R/constants.R b/R/constants.R new file mode 100644 index 00000000..82592b63 --- /dev/null +++ b/R/constants.R @@ -0,0 +1,14 @@ + + + +# This file defines any magic numbers or constants used in the package. + +# The default number of chains to use when sampling from a Stan model. +# This is the default set by the `cmdstanr` package however we need to fix it +# as if they were to change it it would break our code as we create the number +# of initial values based on this number. +# +# Unfortunately the default value is in a method of an unexported object (only the constructor +# is exported) so there is no way for us to access it without digging through the +# internals of `cmdstanr`. +CMDSTAN_DEFAULT_CHAINS <- 4 From 048a242771e77e29f2234710924a8a222da69bae Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 18:00:01 +0000 Subject: [PATCH 15/21] Fix spelling mistakes --- R/settings.R | 2 +- inst/WORDLIST | 1 + man/jmpost-settings.Rd | 2 +- 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/R/settings.R b/R/settings.R index 099acb9f..a5fa2368 100644 --- a/R/settings.R +++ b/R/settings.R @@ -21,7 +21,7 @@ #' initial_value = prior_sample * prior_shrinkage + (1 - prior_shrinkage) * prior_mean #' ``` #' This setting controls the shrinkage factor. A value of 0 means no shrinkage (i.e. -#' pure random draw) whilst a value of 1 means the intial value is just the mean. +#' pure random draw) whilst a value of 1 means the initial value is just the mean. #' #' @examples #' \dontrun{ diff --git a/inst/WORDLIST b/inst/WORDLIST index c9a5df06..8f0966fc 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -116,3 +116,4 @@ xk LogNormal AST throughs +behaviour diff --git a/man/jmpost-settings.Rd b/man/jmpost-settings.Rd index 19c31663..2ca801f6 100644 --- a/man/jmpost-settings.Rd +++ b/man/jmpost-settings.Rd @@ -21,7 +21,7 @@ distribution with a shrinkage factor towards the mean. That is: }\if{html}{\out{}} This setting controls the shrinkage factor. A value of 0 means no shrinkage (i.e. -pure random draw) whilst a value of 1 means the intial value is just the mean. +pure random draw) whilst a value of 1 means the initial value is just the mean. } } \examples{ From 519d049c796cf51fb957b2c255dc3ce1b625baad Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 18:06:33 +0000 Subject: [PATCH 16/21] add missing settings page to pkgdown site --- _pkgdown.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index 9ab3ae93..969880b7 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -27,6 +27,7 @@ reference: - title: Package contents: - jmpost-package + - jmpost-settings - title: Data Specification contents: From daafa567d7ee2c86927fa50a446a73837f65b601 Mon Sep 17 00:00:00 2001 From: gowerc Date: Mon, 5 Feb 2024 18:31:10 +0000 Subject: [PATCH 17/21] Fix issue of options not being set in installed package; document cache directory --- R/StanModule.R | 2 +- R/settings.R | 27 ++++++++++++++++++++++++--- R/zzz.R | 19 +------------------ design/tests/generated_quantities.R | 2 +- design/tests/gsf-no-link.R | 2 +- design/tests/gsf-with-link-identity.R | 2 +- design/tests/gsf-with-link.R | 2 +- design/tests/keri-no-link.R | 2 +- design/tests/kerioui-replication.R | 2 +- design/tests/os-brier-score.R | 2 +- design/tests/os-exponential.R | 2 +- design/tests/os-loglogistic.R | 2 +- design/tests/os-weibull.R | 2 +- design/tests/os-zero-covariates.R | 2 +- design/tests/rs-no-link.R | 2 +- design/tests/rs-with-link.R | 2 +- man/jmpost-settings.Rd | 12 ++++++++++++ tests/testthat/test-compile.R | 4 ++-- vignettes/model_fitting.Rmd | 6 +++--- 19 files changed, 56 insertions(+), 40 deletions(-) diff --git a/R/StanModule.R b/R/StanModule.R index 571c271a..4d291fee 100755 --- a/R/StanModule.R +++ b/R/StanModule.R @@ -167,7 +167,7 @@ setMethod( #' @rdname compileStanModel compileStanModel.StanModule <- function(object) { - exe_dir <- getOption("jmpost.cache.dir") + exe_dir <- getOption("jmpost.cache_dir") if (!dir.exists(exe_dir)) { dir.create(exe_dir, recursive = TRUE) } diff --git a/R/settings.R b/R/settings.R index a5fa2368..9004a0d7 100644 --- a/R/settings.R +++ b/R/settings.R @@ -23,13 +23,34 @@ #' This setting controls the shrinkage factor. A value of 0 means no shrinkage (i.e. #' pure random draw) whilst a value of 1 means the initial value is just the mean. #' +#' ## `jmpost.cache_dir` +#' +#' Default = `tempfile()` +#' +#' Directory to store compiled stan models in. If not set, a temporary directory is used for +#' the given R session. Can also be set via the environment variable `JMPOST_CACHE_DIR`. +#' #' @examples #' \dontrun{ #' options(jmpost.prior_shrinkage = 0.5) #' } #' @name jmpost-settings -NULL +set_options <- function() { + + cache_dir <- Sys.getenv("JMPOST_CACHE_DIR") + + if (cache_dir == "" || is.null(cache_dir)) { + cache_dir <- tempfile() + } -if (is.null(getOption("jmpost.prior_shrinkage", default = NULL))) { - options("jmpost.prior_shrinkage" = 0.5) + current_opts <- names(options()) + jmpost_opts <- list( + jmpost.cache_dir = cache_dir, + jmpost.prior_shrinkage = 0.5 + ) + for (opt in names(jmpost_opts)) { + if (!opt %in% current_opts) { + options(jmpost_opts[opt]) + } + } } diff --git a/R/zzz.R b/R/zzz.R index 013e0cf6..dd419b1a 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -1,23 +1,6 @@ .onLoad <- function(libname, pkgname) { - - cache_dir <- Sys.getenv("JMPOST_CACHE_DIR") - - if (cache_dir == "" || is.null(cache_dir)) { - cache_dir <- tempfile() - } - - current_opts <- names(options()) - jmpost_opts <- list( - jmpost.cache.dir = cache_dir - ) - for (opt in names(jmpost_opts)) { - if (!opt %in% current_opts) { - options(jmpost_opts[opt]) - } - } - - + set_options() } diff --git a/design/tests/generated_quantities.R b/design/tests/generated_quantities.R index 93d0f941..e3b7ba43 100644 --- a/design/tests/generated_quantities.R +++ b/design/tests/generated_quantities.R @@ -7,7 +7,7 @@ library(ggplot2) library(stringr) library(tidyr) -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) #### Example 1 - Fully specified model diff --git a/design/tests/gsf-no-link.R b/design/tests/gsf-no-link.R index d33393d0..c2d94b59 100644 --- a/design/tests/gsf-no-link.R +++ b/design/tests/gsf-no-link.R @@ -9,7 +9,7 @@ library(cmdstanr) # devtools::document() devtools::load_all(export_all = FALSE) -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) #### Example 1 - Fully specified model - using the defaults for everything diff --git a/design/tests/gsf-with-link-identity.R b/design/tests/gsf-with-link-identity.R index d18b2aa1..549147d3 100644 --- a/design/tests/gsf-with-link-identity.R +++ b/design/tests/gsf-with-link-identity.R @@ -9,7 +9,7 @@ library(cmdstanr) devtools::document() devtools::load_all(export_all = FALSE) -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) #### Example 1 - Fully specified model - using the defaults for everything diff --git a/design/tests/gsf-with-link.R b/design/tests/gsf-with-link.R index 859da202..70c62560 100644 --- a/design/tests/gsf-with-link.R +++ b/design/tests/gsf-with-link.R @@ -9,7 +9,7 @@ library(cmdstanr) devtools::document() devtools::load_all(export_all = FALSE) -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) #### Example 1 - Fully specified model - using the defaults for everything diff --git a/design/tests/keri-no-link.R b/design/tests/keri-no-link.R index 93f7ade2..b70b7022 100644 --- a/design/tests/keri-no-link.R +++ b/design/tests/keri-no-link.R @@ -10,7 +10,7 @@ library(cmdstanr) devtools::document() devtools::load_all(export_all = FALSE) -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) ### Try to re-create analysis of Kerioui et. al 2020 diff --git a/design/tests/kerioui-replication.R b/design/tests/kerioui-replication.R index a4fac552..5a19e1fb 100644 --- a/design/tests/kerioui-replication.R +++ b/design/tests/kerioui-replication.R @@ -10,7 +10,7 @@ library(cmdstanr) devtools::document() devtools::load_all(export_all = FALSE) -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) ### Try to re-create analysis of Kerioui et. al 2020 diff --git a/design/tests/os-brier-score.R b/design/tests/os-brier-score.R index c1112bd1..2c4166f4 100644 --- a/design/tests/os-brier-score.R +++ b/design/tests/os-brier-score.R @@ -2,7 +2,7 @@ devtools::document() devtools::load_all() -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) library(bayesplot) diff --git a/design/tests/os-exponential.R b/design/tests/os-exponential.R index 70ece677..77ef5980 100644 --- a/design/tests/os-exponential.R +++ b/design/tests/os-exponential.R @@ -2,7 +2,7 @@ devtools::document() devtools::load_all() -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) library(bayesplot) diff --git a/design/tests/os-loglogistic.R b/design/tests/os-loglogistic.R index 1be7535e..272c67f5 100644 --- a/design/tests/os-loglogistic.R +++ b/design/tests/os-loglogistic.R @@ -3,7 +3,7 @@ devtools::document() devtools::load_all() library(bayesplot) -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) jlist <- simulate_joint_data( n_arm = c(1000, 1000), diff --git a/design/tests/os-weibull.R b/design/tests/os-weibull.R index 2ad062fa..3f90aae6 100644 --- a/design/tests/os-weibull.R +++ b/design/tests/os-weibull.R @@ -3,7 +3,7 @@ devtools::document() devtools::load_all() library(bayesplot) -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) set.seed(6042) diff --git a/design/tests/os-zero-covariates.R b/design/tests/os-zero-covariates.R index 431072c1..cb28d9ac 100644 --- a/design/tests/os-zero-covariates.R +++ b/design/tests/os-zero-covariates.R @@ -1,7 +1,7 @@ devtools::document() devtools::load_all() -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) true_lambda <- 1/100 diff --git a/design/tests/rs-no-link.R b/design/tests/rs-no-link.R index f8a662c0..26c64bba 100644 --- a/design/tests/rs-no-link.R +++ b/design/tests/rs-no-link.R @@ -6,7 +6,7 @@ library(tidyr) devtools::document() devtools::load_all(export_all = FALSE) -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) #### Example 1 - Fully specified model diff --git a/design/tests/rs-with-link.R b/design/tests/rs-with-link.R index a6166e3b..da97c748 100644 --- a/design/tests/rs-with-link.R +++ b/design/tests/rs-with-link.R @@ -6,7 +6,7 @@ library(tidyr) devtools::document() devtools::load_all(export_all = FALSE) -options("jmpost.cache.dir" = file.path("local", "models")) +options("jmpost.cache_dir" = file.path("local", "models")) #### Example 1 - Fully specified model - using the defaults for everything jm <- JointModel( diff --git a/man/jmpost-settings.Rd b/man/jmpost-settings.Rd index 2ca801f6..b9d11ee3 100644 --- a/man/jmpost-settings.Rd +++ b/man/jmpost-settings.Rd @@ -2,7 +2,11 @@ % Please edit documentation in R/settings.R \name{jmpost-settings} \alias{jmpost-settings} +\alias{set_options} \title{jmpost settings} +\usage{ +set_options() +} \description{ Define settings that modify the behaviour of the \code{jmpost} package @@ -23,6 +27,14 @@ distribution with a shrinkage factor towards the mean. That is: This setting controls the shrinkage factor. A value of 0 means no shrinkage (i.e. pure random draw) whilst a value of 1 means the initial value is just the mean. } + +\subsection{\code{jmpost.cache_dir}}{ + +Default = \code{tempfile()} + +Directory to store compiled stan models in. If not set, a temporary directory is used for +the given R session. Can also be set via the environment variable \code{JMPOST_CACHE_DIR}. +} } \examples{ \dontrun{ diff --git a/tests/testthat/test-compile.R b/tests/testthat/test-compile.R index 10e33180..1608a9f4 100644 --- a/tests/testthat/test-compile.R +++ b/tests/testthat/test-compile.R @@ -18,8 +18,8 @@ model { test_that("compileStanModel doesn't error if the directory doesn't exist", { smod <- StanModule(model) fpath <- file.path(tempdir(), "abcd", "efg", "model") - old_path <- options("jmpost.cache.dir") - options("jmpost.cache.dir" = fpath) + old_path <- options("jmpost.cache_dir") + options("jmpost.cache_dir" = fpath) z <- compileStanModel(smod) options(old_path) expect_true(file.exists(fpath)) diff --git a/vignettes/model_fitting.Rmd b/vignettes/model_fitting.Rmd index 5ea6db0c..e82fcfac 100644 --- a/vignettes/model_fitting.Rmd +++ b/vignettes/model_fitting.Rmd @@ -322,7 +322,7 @@ long_quantities <- LongitudinalQuantities( mcmc_results, groups = selected_patients ) -as.data.frame(long_quantities) +as.data.frame(long_quantities) |> head() summary(long_quantities) autoplot(long_quantities) ``` @@ -392,13 +392,13 @@ options("jmpost.prior_shrinkage" = 0.7) If you wish to manually specify the initial values you can do so via the `initialValues()` function. For example: -```{r, eval=FALSE} +```r joint_model <- JointModel( LongitudinalRandomSlope(), SurvivalExponential(), LinkNone() ) -initial_values <- initialValues(joint_model, nchains = 2) +initial_values <- initialValues(joint_model, n_chains = 2) initial_values[[1]]$lm_rs_intercept <- 0.2 initial_values[[2]]$lm_rs_intercept <- 0.3 From 310ae6a7352d30503d91a5ef403dbaad3f6ec76c Mon Sep 17 00:00:00 2001 From: gowerc Date: Tue, 6 Feb 2024 10:16:26 +0000 Subject: [PATCH 18/21] Implemented prior_init_only --- NAMESPACE | 1 - R/LongitudinalGSF.R | 24 +++++++++++++++++++---- R/LongitudinalRandomSlope.R | 6 +++++- R/Prior.R | 14 +++++++++---- _pkgdown.yml | 2 +- man/prior_beta.Rd | 2 +- man/prior_cauchy.Rd | 2 +- man/prior_gamma.Rd | 2 +- man/{prior_none.Rd => prior_init_only.Rd} | 13 ++++++++---- man/prior_invgamma.Rd | 2 +- man/prior_logistic.Rd | 2 +- man/prior_loglogistic.Rd | 2 +- man/prior_lognormal.Rd | 2 +- man/prior_normal.Rd | 2 +- man/prior_std_normal.Rd | 2 +- man/prior_student_t.Rd | 2 +- man/prior_uniform.Rd | 2 +- tests/testthat/_snaps/Prior.md | 2 +- tests/testthat/test-Parameter.R | 2 +- tests/testthat/test-ParameterList.R | 2 +- tests/testthat/test-Prior.R | 4 ++-- tests/testthat/test-initialValues.R | 15 ++++++++++++++ 22 files changed, 76 insertions(+), 31 deletions(-) rename man/{prior_none.Rd => prior_init_only.Rd} (61%) diff --git a/NAMESPACE b/NAMESPACE index bb6e1912..d0166622 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -115,7 +115,6 @@ export(prior_invgamma) export(prior_logistic) export(prior_loglogistic) export(prior_lognormal) -export(prior_none) export(prior_normal) export(prior_std_normal) export(prior_student_t) diff --git a/R/LongitudinalGSF.R b/R/LongitudinalGSF.R index 6a85fa55..f7eb59b1 100755 --- a/R/LongitudinalGSF.R +++ b/R/LongitudinalGSF.R @@ -73,7 +73,11 @@ LongitudinalGSF <- function( Parameter(name = "lm_gsf_a_phi", prior = a_phi, size = "n_arms"), Parameter(name = "lm_gsf_b_phi", prior = b_phi, size = "n_arms"), - Parameter(name = "lm_gsf_psi_phi", prior = prior_none(), size = "Nind"), + Parameter( + name = "lm_gsf_psi_phi", + prior = prior_init_only(prior_beta(a_phi@init, b_phi@init)), + size = "Nind" + ), Parameter(name = "lm_gsf_sigma", prior = sigma, size = 1) ) @@ -81,9 +85,21 @@ LongitudinalGSF <- function( assert_flag(centered) parameters_extra <- if (centered) { list( - Parameter(name = "lm_gsf_psi_bsld", prior = prior_none(), size = "Nind"), - Parameter(name = "lm_gsf_psi_ks", prior = prior_none(), size = "Nind"), - Parameter(name = "lm_gsf_psi_kg", prior = prior_none(), size = "Nind") + Parameter( + name = "lm_gsf_psi_bsld", + prior = prior_init_only(prior_lognormal(mu_bsld@init, omega_bsld@init)), + size = "Nind" + ), + Parameter( + name = "lm_gsf_psi_ks", + prior = prior_init_only(prior_lognormal(mu_ks@init, omega_ks@init)), + size = "Nind" + ), + Parameter( + name = "lm_gsf_psi_kg", + prior = prior_init_only(prior_lognormal(mu_kg@init, omega_kg@init)), + size = "Nind" + ) ) } else { list( diff --git a/R/LongitudinalRandomSlope.R b/R/LongitudinalRandomSlope.R index eb5a866f..70e24340 100755 --- a/R/LongitudinalRandomSlope.R +++ b/R/LongitudinalRandomSlope.R @@ -46,7 +46,11 @@ LongitudinalRandomSlope <- function( Parameter(name = "lm_rs_slope_mu", prior = slope_mu, size = "n_arms"), Parameter(name = "lm_rs_slope_sigma", prior = slope_sigma, size = 1), Parameter(name = "lm_rs_sigma", prior = sigma, size = 1), - Parameter(name = "lm_rs_ind_rnd_slope", prior = prior_none(), size = "Nind") + Parameter( + name = "lm_rs_ind_rnd_slope", + prior = prior_init_only(prior_normal(slope_mu@init, slope_sigma@init)), + size = "Nind" + ) ) ) ) diff --git a/R/Prior.R b/R/Prior.R index 0d00d186..41b397be 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -339,17 +339,23 @@ prior_beta <- function(a, b) { #' Only Initial Values Specification #' +#' @param dist (`Prior`)\cr a prior Distribution #' @family Prior +#' @description +#' This function is used to specify only the initial values for a parameter. +#' This is primarily used for heiracrhical parameters whose distributions +#' are fixed within the model and cannot be altered by the user. #' -#' @export -prior_none <- function() { +prior_init_only <- function(dist) { Prior( parameters = list(), display = "", repr_model = "", repr_data = "", - sample = \(n) local_runif(n, -4, 4), - init = 0, + sample = \(n) { + dist@sample(n) + }, + init = dist@init, validation = list() ) } diff --git a/_pkgdown.yml b/_pkgdown.yml index 969880b7..65edab30 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -55,11 +55,11 @@ reference: - prior_logistic - prior_loglogistic - prior_lognormal - - prior_none - prior_normal - prior_uniform - prior_std_normal - prior_student_t + - prior_init_only - title: Longitudinal Model Specification contents: diff --git a/man/prior_beta.Rd b/man/prior_beta.Rd index f04aaf01..b51e5142 100644 --- a/man/prior_beta.Rd +++ b/man/prior_beta.Rd @@ -18,11 +18,11 @@ Beta Prior Distribution Other Prior: \code{\link{prior_cauchy}()}, \code{\link{prior_gamma}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_invgamma}()}, \code{\link{prior_logistic}()}, \code{\link{prior_loglogistic}()}, \code{\link{prior_lognormal}()}, -\code{\link{prior_none}()}, \code{\link{prior_normal}()}, \code{\link{prior_std_normal}()}, \code{\link{prior_student_t}()}, diff --git a/man/prior_cauchy.Rd b/man/prior_cauchy.Rd index f57b71ab..ac5272b6 100644 --- a/man/prior_cauchy.Rd +++ b/man/prior_cauchy.Rd @@ -18,11 +18,11 @@ Cauchy Prior Distribution Other Prior: \code{\link{prior_beta}()}, \code{\link{prior_gamma}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_invgamma}()}, \code{\link{prior_logistic}()}, \code{\link{prior_loglogistic}()}, \code{\link{prior_lognormal}()}, -\code{\link{prior_none}()}, \code{\link{prior_normal}()}, \code{\link{prior_std_normal}()}, \code{\link{prior_student_t}()}, diff --git a/man/prior_gamma.Rd b/man/prior_gamma.Rd index a1ddaced..f573e736 100644 --- a/man/prior_gamma.Rd +++ b/man/prior_gamma.Rd @@ -18,11 +18,11 @@ Gamma Prior Distribution Other Prior: \code{\link{prior_beta}()}, \code{\link{prior_cauchy}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_invgamma}()}, \code{\link{prior_logistic}()}, \code{\link{prior_loglogistic}()}, \code{\link{prior_lognormal}()}, -\code{\link{prior_none}()}, \code{\link{prior_normal}()}, \code{\link{prior_std_normal}()}, \code{\link{prior_student_t}()}, diff --git a/man/prior_none.Rd b/man/prior_init_only.Rd similarity index 61% rename from man/prior_none.Rd rename to man/prior_init_only.Rd index 48787a80..e04bb9c0 100644 --- a/man/prior_none.Rd +++ b/man/prior_init_only.Rd @@ -1,13 +1,18 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/Prior.R -\name{prior_none} -\alias{prior_none} +\name{prior_init_only} +\alias{prior_init_only} \title{Only Initial Values Specification} \usage{ -prior_none() +prior_init_only(dist) +} +\arguments{ +\item{dist}{(\code{Prior})\cr a prior Distribution} } \description{ -Only Initial Values Specification +This function is used to specify only the initial values for a parameter. +This is primarily used for heiracrhical parameters whose distributions +are fixed within the model and cannot be altered by the user. } \seealso{ Other Prior: diff --git a/man/prior_invgamma.Rd b/man/prior_invgamma.Rd index dbc94b9b..f59e610d 100644 --- a/man/prior_invgamma.Rd +++ b/man/prior_invgamma.Rd @@ -19,10 +19,10 @@ Other Prior: \code{\link{prior_beta}()}, \code{\link{prior_cauchy}()}, \code{\link{prior_gamma}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_logistic}()}, \code{\link{prior_loglogistic}()}, \code{\link{prior_lognormal}()}, -\code{\link{prior_none}()}, \code{\link{prior_normal}()}, \code{\link{prior_std_normal}()}, \code{\link{prior_student_t}()}, diff --git a/man/prior_logistic.Rd b/man/prior_logistic.Rd index 74f420d5..0e0aac7f 100644 --- a/man/prior_logistic.Rd +++ b/man/prior_logistic.Rd @@ -19,10 +19,10 @@ Other Prior: \code{\link{prior_beta}()}, \code{\link{prior_cauchy}()}, \code{\link{prior_gamma}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_invgamma}()}, \code{\link{prior_loglogistic}()}, \code{\link{prior_lognormal}()}, -\code{\link{prior_none}()}, \code{\link{prior_normal}()}, \code{\link{prior_std_normal}()}, \code{\link{prior_student_t}()}, diff --git a/man/prior_loglogistic.Rd b/man/prior_loglogistic.Rd index 3f84ff18..97a961da 100644 --- a/man/prior_loglogistic.Rd +++ b/man/prior_loglogistic.Rd @@ -19,10 +19,10 @@ Other Prior: \code{\link{prior_beta}()}, \code{\link{prior_cauchy}()}, \code{\link{prior_gamma}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_invgamma}()}, \code{\link{prior_logistic}()}, \code{\link{prior_lognormal}()}, -\code{\link{prior_none}()}, \code{\link{prior_normal}()}, \code{\link{prior_std_normal}()}, \code{\link{prior_student_t}()}, diff --git a/man/prior_lognormal.Rd b/man/prior_lognormal.Rd index d4140540..525c55a2 100644 --- a/man/prior_lognormal.Rd +++ b/man/prior_lognormal.Rd @@ -19,10 +19,10 @@ Other Prior: \code{\link{prior_beta}()}, \code{\link{prior_cauchy}()}, \code{\link{prior_gamma}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_invgamma}()}, \code{\link{prior_logistic}()}, \code{\link{prior_loglogistic}()}, -\code{\link{prior_none}()}, \code{\link{prior_normal}()}, \code{\link{prior_std_normal}()}, \code{\link{prior_student_t}()}, diff --git a/man/prior_normal.Rd b/man/prior_normal.Rd index 1869a014..a597e573 100644 --- a/man/prior_normal.Rd +++ b/man/prior_normal.Rd @@ -19,11 +19,11 @@ Other Prior: \code{\link{prior_beta}()}, \code{\link{prior_cauchy}()}, \code{\link{prior_gamma}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_invgamma}()}, \code{\link{prior_logistic}()}, \code{\link{prior_loglogistic}()}, \code{\link{prior_lognormal}()}, -\code{\link{prior_none}()}, \code{\link{prior_std_normal}()}, \code{\link{prior_student_t}()}, \code{\link{prior_uniform}()} diff --git a/man/prior_std_normal.Rd b/man/prior_std_normal.Rd index 614e8375..c35525c7 100644 --- a/man/prior_std_normal.Rd +++ b/man/prior_std_normal.Rd @@ -14,11 +14,11 @@ Other Prior: \code{\link{prior_beta}()}, \code{\link{prior_cauchy}()}, \code{\link{prior_gamma}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_invgamma}()}, \code{\link{prior_logistic}()}, \code{\link{prior_loglogistic}()}, \code{\link{prior_lognormal}()}, -\code{\link{prior_none}()}, \code{\link{prior_normal}()}, \code{\link{prior_student_t}()}, \code{\link{prior_uniform}()} diff --git a/man/prior_student_t.Rd b/man/prior_student_t.Rd index bfcb6162..b22f3bec 100644 --- a/man/prior_student_t.Rd +++ b/man/prior_student_t.Rd @@ -21,11 +21,11 @@ Other Prior: \code{\link{prior_beta}()}, \code{\link{prior_cauchy}()}, \code{\link{prior_gamma}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_invgamma}()}, \code{\link{prior_logistic}()}, \code{\link{prior_loglogistic}()}, \code{\link{prior_lognormal}()}, -\code{\link{prior_none}()}, \code{\link{prior_normal}()}, \code{\link{prior_std_normal}()}, \code{\link{prior_uniform}()} diff --git a/man/prior_uniform.Rd b/man/prior_uniform.Rd index af57bfe8..d0a29028 100644 --- a/man/prior_uniform.Rd +++ b/man/prior_uniform.Rd @@ -19,11 +19,11 @@ Other Prior: \code{\link{prior_beta}()}, \code{\link{prior_cauchy}()}, \code{\link{prior_gamma}()}, +\code{\link{prior_init_only}()}, \code{\link{prior_invgamma}()}, \code{\link{prior_logistic}()}, \code{\link{prior_loglogistic}()}, \code{\link{prior_lognormal}()}, -\code{\link{prior_none}()}, \code{\link{prior_normal}()}, \code{\link{prior_std_normal}()}, \code{\link{prior_student_t}()} diff --git a/tests/testthat/_snaps/Prior.md b/tests/testthat/_snaps/Prior.md index 656ecbda..908f17bb 100644 --- a/tests/testthat/_snaps/Prior.md +++ b/tests/testthat/_snaps/Prior.md @@ -51,7 +51,7 @@ --- Code - print(prior_none()) + print(prior_init_only(prior_normal(1, 4))) Output Prior Object: diff --git a/tests/testthat/test-Parameter.R b/tests/testthat/test-Parameter.R index 1637b469..b81bb926 100644 --- a/tests/testthat/test-Parameter.R +++ b/tests/testthat/test-Parameter.R @@ -19,6 +19,6 @@ test_that("show() works for Paramneter objects", { x <- Parameter(prior_beta(0.5, 0.2), "var1", "size1") expect_snapshot(print(x)) - x <- Parameter(prior_none(), "x", "size1") + x <- Parameter(prior_init_only(prior_normal(0, 1)), "x", "size1") expect_snapshot(print(x)) }) diff --git a/tests/testthat/test-ParameterList.R b/tests/testthat/test-ParameterList.R index eec6abab..968fea45 100644 --- a/tests/testthat/test-ParameterList.R +++ b/tests/testthat/test-ParameterList.R @@ -47,7 +47,7 @@ test_that("show() works for ParameterList objects", { Parameter(name = "bob", prior = prior_normal(1, 4)), Parameter(name = "sam", prior = prior_beta(3, 1)), Parameter(name = "dave", prior = prior_lognormal(3, 2), size = 4), - Parameter(name = "steve", prior = prior_none()) + Parameter(name = "steve", prior = prior_init_only(prior_normal(4, 2))) ) expect_snapshot(print(x)) diff --git a/tests/testthat/test-Prior.R b/tests/testthat/test-Prior.R index 1defd442..3e434193 100644 --- a/tests/testthat/test-Prior.R +++ b/tests/testthat/test-Prior.R @@ -137,7 +137,7 @@ test_that("Invalid prior parameters are rejected", { # Ensure that validation doesn't wrongly reject priors with no user specified parameters - expect_s4_class(prior_none(), "Prior") + expect_s4_class(prior_init_only(prior_normal(3, 1)), "Prior") expect_s4_class(prior_std_normal(), "Prior") }) @@ -149,7 +149,7 @@ test_that("show() works for Prior objects", { expect_snapshot(print(prior_std_normal())) expect_snapshot(print(prior_beta(5, 1))) expect_snapshot(print(prior_gamma(2.56, 12))) - expect_snapshot(print(prior_none())) + expect_snapshot(print(prior_init_only(prior_normal(1, 4)))) expect_snapshot(print(prior_uniform(8, 10))) expect_snapshot(print(prior_student_t(3, 10, 4))) expect_snapshot(print(prior_logistic(sigma = 2, 10))) diff --git a/tests/testthat/test-initialValues.R b/tests/testthat/test-initialValues.R index 2a52918f..8a6332cb 100644 --- a/tests/testthat/test-initialValues.R +++ b/tests/testthat/test-initialValues.R @@ -81,3 +81,18 @@ test_that("ensure_initial_values() works as expected", { expect_equal(res[[1]]$p2, array(c(2, 3), dim = c(2))) expect_equal(res[[1]]$p3, array(c(4, 5, 6), dim = c(3))) }) + + +test_that("intial values for fixed distributions gives valid values", { + + set.seed(3150) + gsfmodel <- LongitudinalGSF(centered = TRUE) + ivs <- initialValues(gsfmodel, n_chains = 100) + + for (values in ivs) { + expect_true(values$lm_gsf_psi_phi > 0 & values$lm_gsf_psi_phi < 1) + expect_true(values$lm_gsf_psi_bsld > 0) + expect_true(values$lm_gsf_psi_ks > 0) + expect_true(values$lm_gsf_psi_kg > 0) + } +}) From 44df406d89e0008fd0aa1326ca1573851e492aa5 Mon Sep 17 00:00:00 2001 From: gowerc Date: Tue, 6 Feb 2024 10:22:05 +0000 Subject: [PATCH 19/21] Fix spelling mistake --- R/Prior.R | 2 +- man/prior_init_only.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/Prior.R b/R/Prior.R index 41b397be..abb3149d 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -343,7 +343,7 @@ prior_beta <- function(a, b) { #' @family Prior #' @description #' This function is used to specify only the initial values for a parameter. -#' This is primarily used for heiracrhical parameters whose distributions +#' This is primarily used for hierarchical parameters whose distributions #' are fixed within the model and cannot be altered by the user. #' prior_init_only <- function(dist) { diff --git a/man/prior_init_only.Rd b/man/prior_init_only.Rd index e04bb9c0..a7a27ffe 100644 --- a/man/prior_init_only.Rd +++ b/man/prior_init_only.Rd @@ -11,7 +11,7 @@ prior_init_only(dist) } \description{ This function is used to specify only the initial values for a parameter. -This is primarily used for heiracrhical parameters whose distributions +This is primarily used for hierarchical parameters whose distributions are fixed within the model and cannot be altered by the user. } \seealso{ From 08e5b21a88113746f7fec7a8097c711068d7331f Mon Sep 17 00:00:00 2001 From: gowerc Date: Tue, 6 Feb 2024 11:09:19 +0000 Subject: [PATCH 20/21] Add documentation to priors --- R/Prior.R | 28 ++++++++++++++++++++++++++++ man/Local_Sample.Rd | 30 ++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+) diff --git a/R/Prior.R b/R/Prior.R index abb3149d..e15375aa 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -536,6 +536,34 @@ prior_invgamma <- function(alpha, beta) { #' #' @importFrom stats rbeta rcauchy rgamma rlnorm rlogis rnorm rt runif #' +#' @details +#' +#' ## Log-Logistic +#' +#' There is no log-logistic sampling function within base R so it was implemented +#' in terms of sampling from the CDF distribution. Using the Stan parameterisation +#' the CDF is defined as: +#' \deqn{ +#' u = F(x) = \frac{1}{1 + (x/ \alpha)^{-\beta}} +#' } +#' The inverse of this function is: +#' \deqn{ +#' x = ((u / (1 - u))^(1 / beta)) * alpha +#' } +#' +#' Thus we can sample u from a \eqn{Uni(0, 1)} distribution and then derive x from this. +#' +#' ## Inverse-Gamma +#' +#' The inverse Gamma distribution is defined as 1/Gamma thus we calculate this simply +#' by sampling sampling from the Gamma distribution and then taking the reciprocal. +#' +#' ## Student-t +#' +#' R's sampling functions only produce the standard Student-t distribution so in order +#' to match Stan's implementation we multiply by the scale parameter and add the location +#' parameter. See https://stats.stackexchange.com/a/623611 for details +#' #' @name Local_Sample #' @keywords internal NULL diff --git a/man/Local_Sample.Rd b/man/Local_Sample.Rd index ddf7f12b..9927c9f1 100644 --- a/man/Local_Sample.Rd +++ b/man/Local_Sample.Rd @@ -53,4 +53,34 @@ tests in order to provide deterministic values. In most cases these are just straight forward pass throughs for the underlying distributions. } +\details{ +\subsection{Log-Logistic}{ + +There is no log-logistic sampling function within base R so it was implemented +in terms of sampling from the CDF distribution. Using the Stan parameterisation +the CDF is defined as: +\deqn{ +u = F(x) = \frac{1}{1 + (x/ \alpha)^{-\beta}} +} +The inverse of this function is: +\deqn{ +x = ((u / (1 - u))^(1 / beta)) * alpha +} + +Thus we can sample u from a \eqn{Uni(0, 1)} distribution and then derive x from this. +} + +\subsection{Inverse-Gamma}{ + +The inverse Gamma distribution is defined as 1/Gamma thus we calculate this simply +by sampling sampling from the Gamma distribution and then taking the reciprocal. +} + +\subsection{Student-t}{ + +R's sampling functions only produce the standard Student-t distribution so in order +to match Stan's implementation we multiply by the scale parameter and add the location +parameter. See https://stats.stackexchange.com/a/623611 for details +} +} \keyword{internal} From 18aaae5dc239e61c73f4bff1dd461caf2969c443 Mon Sep 17 00:00:00 2001 From: gowerc Date: Tue, 6 Feb 2024 11:13:50 +0000 Subject: [PATCH 21/21] add link syntax to docs --- R/Prior.R | 3 ++- man/Local_Sample.Rd | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/R/Prior.R b/R/Prior.R index e15375aa..9708844a 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -562,7 +562,8 @@ prior_invgamma <- function(alpha, beta) { #' #' R's sampling functions only produce the standard Student-t distribution so in order #' to match Stan's implementation we multiply by the scale parameter and add the location -#' parameter. See https://stats.stackexchange.com/a/623611 for details +#' parameter. See this \href{https://stats.stackexchange.com/a/623611}{Stack Overflow} post +#' for details #' #' @name Local_Sample #' @keywords internal diff --git a/man/Local_Sample.Rd b/man/Local_Sample.Rd index 9927c9f1..8624d95e 100644 --- a/man/Local_Sample.Rd +++ b/man/Local_Sample.Rd @@ -80,7 +80,8 @@ by sampling sampling from the Gamma distribution and then taking the reciprocal. R's sampling functions only produce the standard Student-t distribution so in order to match Stan's implementation we multiply by the scale parameter and add the location -parameter. See https://stats.stackexchange.com/a/623611 for details +parameter. See this \href{https://stats.stackexchange.com/a/623611}{Stack Overflow} post +for details } } \keyword{internal}