diff --git a/DESCRIPTION b/DESCRIPTION index cff59d90..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' @@ -88,6 +89,7 @@ Collate: 'defaults.R' 'external-exports.R' 'jmpost-package.R' + 'settings.R' 'simulations.R' 'simulations_gsf.R' 'simulations_os.R' diff --git a/NAMESPACE b/NAMESPACE index 2dcd3bfe..d0166622 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) @@ -112,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) @@ -162,5 +164,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/JointModel.R b/R/JointModel.R index 20d634de..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 @@ -136,6 +137,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 +153,37 @@ 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 { + # Magic constant from R/constants.R + CMDSTAN_DEFAULT_CHAINS } + initial_values <- if ("init" %in% names(args)) { + args[["init"]] + } else { + initialValues(object, n_chains = 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 +193,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, n_chains, ...) { + initialValues(object@parameters, n_chains) } diff --git a/R/Link.R b/R/Link.R index 663bb950..e85ecd4c 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, n_chains, ...) { + initialValues(object@parameters, n_chains) } 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..f7eb59b1 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,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 = psi_phi, 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) ) @@ -95,9 +85,21 @@ 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_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 9002f9fe..70e24340 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, 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) ) { 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,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 = random_slope, 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/Parameter.R b/R/Parameter.R index 9061ad4c..f9fa0ca7 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 @@ -115,7 +116,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..1c435d7c 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 n_chains (`integer`) \cr the number of chains. #' @name ParameterList-Getter-Methods NULL @@ -145,11 +146,19 @@ 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, n_chains, ...) { + # 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 + } + ) } diff --git a/R/Prior.R b/R/Prior.R index 4135ba9f..9708844a 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,10 @@ NULL #' @describeIn Prior-Getter-Methods The prior's initial value #' @export -initialValues.Prior <- function(object) object@init +initialValues.Prior <- function(object, ...) { + getOption("jmpost.prior_shrinkage") * object@init + + (1 - getOption("jmpost.prior_shrinkage")) * object@sample(1) +} # Prior-constructors ---- @@ -184,10 +199,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 +210,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 +222,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 +241,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 +253,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 +266,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 +278,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 +291,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 +303,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 +316,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 +328,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 +339,23 @@ prior_beta <- function(a, b, init = a / (a + b)) { #' Only Initial Values Specification #' -#' @inheritParams Prior-Shared +#' @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 hierarchical parameters whose distributions +#' are fixed within the model and cannot be altered by the user. #' -#' @export -prior_none <- function(init = 0.00001) { +prior_init_only <- function(dist) { Prior( parameters = list(), display = "", repr_model = "", repr_data = "", - init = init, + sample = \(n) { + dist@sample(n) + }, + init = dist@init, validation = list() ) } @@ -346,11 +367,14 @@ 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) { + assert_that( + alpha < beta, + msg = "`alpha`` must be less than `beta`" + ) Prior( parameters = list(alpha = alpha, beta = beta), display = "uniform(alpha = {alpha}, beta = {beta})", @@ -359,7 +383,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 +398,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 +415,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 +426,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 +446,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 +460,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 +475,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 +491,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 +506,102 @@ 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. +#' +#' @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 this \href{https://stats.stackexchange.com/a/623611}{Stack Overflow} post +#' for details +#' +#' @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..f05e07af 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, n_chains, ...) { + initialValues(object@parameters, n_chains) +} 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/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/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 diff --git a/R/generics.R b/R/generics.R index e273f953..dcd640a2 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 n_chains 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 `n_chains` 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/R/settings.R b/R/settings.R new file mode 100644 index 00000000..9004a0d7 --- /dev/null +++ b/R/settings.R @@ -0,0 +1,56 @@ + + + +#' 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 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 +set_options <- function() { + + 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, + 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/_pkgdown.yml b/_pkgdown.yml index e626c09c..65edab30 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -27,6 +27,7 @@ reference: - title: Package contents: - jmpost-package + - jmpost-settings - title: Data Specification contents: @@ -54,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: @@ -122,6 +123,7 @@ reference: - title: Convenience Functions contents: + - initialValues - merge - show - as_stan_list 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/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/inst/WORDLIST b/inst/WORDLIST index ccdcf4df..8f0966fc 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -115,3 +115,5 @@ gk xk LogNormal AST +throughs +behaviour 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..8624d95e --- /dev/null +++ b/man/Local_Sample.Rd @@ -0,0 +1,87 @@ +% 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. +} +\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 this \href{https://stats.stackexchange.com/a/623611}{Stack Overflow} post +for details +} +} +\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..38d093d4 100644 --- a/man/LongitudinalRandomSlope-class.Rd +++ b/man/LongitudinalRandomSlope-class.Rd @@ -8,11 +8,10 @@ \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) ) } \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/Parameter-Getter-Methods.Rd b/man/Parameter-Getter-Methods.Rd index 4cc1e88b..06d2b76e 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) } @@ -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 be530325..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) +\method{initialValues}{ParameterList}(object, n_chains, ...) \method{size}{ParameterList}(object) } @@ -17,6 +17,10 @@ \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{n_chains}{(\code{integer}) \cr the number of chains.} + +\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/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/man/initialValues.Rd b/man/initialValues.Rd index 65197d40..25536f89 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, n_chains, ...) -\method{initialValues}{JointModel}(object) +\method{initialValues}{Link}(object, n_chains, ...) + +\method{initialValues}{JointModel}(object, n_chains, ...) } \arguments{ \item{object}{where to get the initial values from.} + +\item{...}{Not currently used.} + +\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. } -\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{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/jmpost-settings.Rd b/man/jmpost-settings.Rd new file mode 100644 index 00000000..b9d11ee3 --- /dev/null +++ b/man/jmpost-settings.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% 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 + +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 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{ +options(jmpost.prior_shrinkage = 0.5) +} +} 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..b51e5142 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 @@ -20,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 85614625..ac5272b6 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 @@ -20,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 a1e7667c..f573e736 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 @@ -20,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 62% rename from man/prior_none.Rd rename to man/prior_init_only.Rd index b5e272f0..a7a27ffe 100644 --- a/man/prior_none.Rd +++ b/man/prior_init_only.Rd @@ -1,16 +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(init = 1e-05) +prior_init_only(dist) } \arguments{ -\item{init}{(\code{number})\cr initial value.} +\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 hierarchical 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 82fdaaa7..f59e610d 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 @@ -21,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 cdaadae0..0e0aac7f 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 @@ -21,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 89581e94..97a961da 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 @@ -21,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 9aebbcae..525c55a2 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 @@ -21,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 6f311e90..a597e573 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 @@ -21,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 2dcee02d..c35525c7 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 @@ -17,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 939d1d94..b22f3bec 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 @@ -23,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 c9b926c2..d0a29028 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 @@ -21,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/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..908f17bb 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: @@ -51,7 +51,7 @@ --- Code - print(prior_none()) + print(prior_init_only(prior_normal(1, 4))) 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..55fda0f3 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, n_chains = 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..b81bb926 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_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 ea2fcf8f..968fea45 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, n_chains = 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,8 +46,8 @@ 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 = "steve", prior = prior_none()) + Parameter(name = "dave", prior = prior_lognormal(3, 2), size = 4), + 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 e2ee6032..3e434193 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) @@ -142,24 +130,59 @@ 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") + expect_s4_class(prior_init_only(prior_normal(3, 1)), "Prior") expect_s4_class(prior_std_normal(), "Prior") }) 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))) 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, 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))) }) + + +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.5) + }, + local_rnorm = \(...) 4 + ) +}) 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-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/tests/testthat/test-initialValues.R b/tests/testthat/test-initialValues.R new file mode 100644 index 00000000..8a6332cb --- /dev/null +++ b/tests/testthat/test-initialValues.R @@ -0,0 +1,98 @@ + + +test_that("initialValues() works as expected", { + jm <- JointModel( + longitudinal = LongitudinalRandomSlope(), + survival = SurvivalWeibullPH(), + link = LinkRandomSlope() + ) + + set.seed(341) + initial_values <- initialValues(jm, n_chains = 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, n_chains = 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, n_chains = 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))) +}) + + +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) + } +}) diff --git a/vignettes/model_fitting.Rmd b/vignettes/model_fitting.Rmd index 62e6716f..e82fcfac 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 @@ -325,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) ``` @@ -373,3 +370,55 @@ 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 * 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: + +```r +joint_model <- JointModel( + LongitudinalRandomSlope(), + SurvivalExponential(), + LinkNone() +) +initial_values <- initialValues(joint_model, n_chains = 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. + + +