From 32e9cb510bc3393cd53bca9da71d875596051a81 Mon Sep 17 00:00:00 2001 From: Ven Popov Date: Fri, 22 Mar 2024 19:58:47 +0100 Subject: [PATCH 1/3] add failing test for no fixed parameters --- tests/testthat/test-helpers-prior.R | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/testthat/test-helpers-prior.R b/tests/testthat/test-helpers-prior.R index 4fdaf6fa..c9c78798 100644 --- a/tests/testthat/test-helpers-prior.R +++ b/tests/testthat/test-helpers-prior.R @@ -95,3 +95,18 @@ test_that("no check for sort_data with default_priors function", { sdmSimple('dev_rad'))) expect_false(any(grepl("sort", res))) }) + + +test_that("default priors work when there are no fixed parameters", { + formula <- bmf(mu ~ 1, + c ~ 1, + kappa ~ 1) + if (utils::packageVersion("brms") >= "2.20.14") { + prior_fn <- default_prior + } else { + prior_fn <- get_model_prior + } + + pr <- prior_fn(formula, OberauerLin_2017, sdmSimple('dev_rad')) + expect_s3_class(pr, 'brmsprior') +}) From 45af221c58e21fdde7a9f2419a87a5d802242a9d Mon Sep 17 00:00:00 2001 From: Ven Popov Date: Sat, 23 Mar 2024 01:17:40 +0100 Subject: [PATCH 2/3] improve stability of k2sd --- R/helpers-parameters.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/helpers-parameters.R b/R/helpers-parameters.R index af6f34ec..f363f581 100644 --- a/R/helpers-parameters.R +++ b/R/helpers-parameters.R @@ -27,7 +27,7 @@ k2sd <- function(K) { if (K[j] == 0) S[j] <- Inf if (is.infinite(K[j])) S[j] <- 0 if (K[j] >= 0 & !is.infinite(K[j])) { - S[j] <- sqrt(-2 * log(besselI(K[j], 1) / besselI(K[j], 0))) + S[j] <- sqrt(-2 * log(besselI(K[j], 1, expon.scaled = T) / besselI(K[j], 0, expon.scaled = T))) } } as.numeric(S) From 4ca431c892b004d436a418c3bf162a5b144e8fd4 Mon Sep 17 00:00:00 2001 From: Ven Popov Date: Sat, 23 Mar 2024 17:23:28 +0100 Subject: [PATCH 3/3] fix need for at least one constant parameter - add default priors for mu - tan_half link for mu in SDM --- R/bmm_model_IMM.R | 3 +++ R/bmm_model_mixture2p.R | 1 + R/bmm_model_mixture3p.R | 1 + R/bmm_model_sdmSimple.R | 9 +++++---- R/helpers-prior.R | 13 +++++++++---- inst/stan_chunks/sdmSimple_funs.stan | 2 ++ man/SDM.Rd | 2 +- tests/testthat/test-summary.R | 2 +- 8 files changed, 23 insertions(+), 10 deletions(-) diff --git a/R/bmm_model_IMM.R b/R/bmm_model_IMM.R index 8d8acec1..57847673 100644 --- a/R/bmm_model_IMM.R +++ b/R/bmm_model_IMM.R @@ -39,6 +39,7 @@ ), fixed_parameters = list(mu1 = 0, mu2 = 0, kappa2 = -100), default_priors = list( + mu1 = list(main = "student_t(1, 0, 1)"), kappa = list(main = "normal(2,1)", effects = "normal(0,1)"), a = list(main = "normal(0,1)", effects = "normal(0,1)"), c = list(main = "normal(0,1)", effects = "normal(0,1)") @@ -90,6 +91,7 @@ s = "log"), fixed_parameters = list(mu1 = 0, mu2 = 0, kappa2 = -100), default_priors = list( + mu1 = list(main = "student_t(1, 0, 1)"), kappa = list(main = "normal(2,1)", effects = "normal(0,1)"), c = list(main = "normal(0,1)", effects = "normal(0,1)"), s = list(main = "normal(0,1)", effects = "normal(0,1)") @@ -144,6 +146,7 @@ ), fixed_parameters = list(mu1 = 0, mu2 = 0, kappa2 = -100), default_priors = list( + mu1 = list(main = "student_t(1, 0, 1)"), kappa = list(main = "normal(2,1)", effects = "normal(0,1)"), a = list(main = "normal(0,1)", effects = "normal(0,1)"), c = list(main = "normal(0,1)", effects = "normal(0,1)"), diff --git a/R/bmm_model_mixture2p.R b/R/bmm_model_mixture2p.R index 7684b58d..bb358cec 100644 --- a/R/bmm_model_mixture2p.R +++ b/R/bmm_model_mixture2p.R @@ -36,6 +36,7 @@ ), fixed_parameters = list(mu1 = 0, mu2 = 0, kappa2 = -100), default_priors = list( + mu1 = list(main = "student_t(1, 0, 1)"), kappa = list(main = "normal(2, 1)", effects = "normal(0, 1)"), thetat = list(main = "logistic(0, 1)") ), diff --git a/R/bmm_model_mixture3p.R b/R/bmm_model_mixture3p.R index 6ea96906..a1d7d857 100644 --- a/R/bmm_model_mixture3p.R +++ b/R/bmm_model_mixture3p.R @@ -41,6 +41,7 @@ ), fixed_parameters = list(mu1 = 0, mu2 = 0, kappa2 = -100), default_priors = list( + mu1 = list(main = "student_t(1, 0, 1)"), kappa = list(main = "normal(2,1)", effects = "normal(0,1)"), thetat = list(main = "logistic(0, 1)"), thetant = list(main = "logistic(0, 1)") diff --git a/R/bmm_model_sdmSimple.R b/R/bmm_model_sdmSimple.R index 6efbac65..666a70cd 100644 --- a/R/bmm_model_sdmSimple.R +++ b/R/bmm_model_sdmSimple.R @@ -25,12 +25,13 @@ kappa = 'Precision parameter of the SDM distribution' ), links = list( - mu = 'identity', + mu = 'tan_half', c = 'log', kappa = 'log' ), fixed_parameters = list(mu = 0), default_priors = list( + mu = list(main = "student_t(1, 0, 1)"), kappa = list(main = "student_t(5,1.75,0.75)", effects = "normal(0,1)"), c = list(main = "student_t(5,2,0.75)", effects = "normal(0,1)") ), @@ -123,9 +124,9 @@ configure_model.sdmSimple <- function(model, data, formula) { sdm_simple <- brms::custom_family( "sdm_simple", dpars = c("mu", "c", "kappa"), - links = c("identity", "identity", "log"), - lb = c(-pi, NA, NA), - ub = c(pi, NA, NA), + links = c("tan_half", "identity", "log"), + lb = c(NA, NA, NA), + ub = c(NA, NA, NA), type = "real", loop = FALSE, log_lik = log_lik_sdm_simple, posterior_predict = posterior_predict_sdm_simple diff --git a/R/helpers-prior.R b/R/helpers-prior.R index 41f1e4ab..21763d70 100644 --- a/R/helpers-prior.R +++ b/R/helpers-prior.R @@ -95,10 +95,10 @@ get_model_prior <- function(object, data, model, formula = object, ...) { #' model #' @noRd fixed_pars_priors <- function(model, formula, additional_pars = list()) { - # determine type of parameters - bterms <- brms::brmsterms(formula) - dpars <- names(bterms$dpars) - nlpars <- names(bterms$nlpars) + fix_pars <- model$fixed_parameters + if (length(fix_pars) == 0) { + return(brms::empty_prior()) + } # construct parameter names and prior values par_list <- c(model$fixed_parameters, additional_pars) @@ -106,6 +106,11 @@ fixed_pars_priors <- function(model, formula, additional_pars = list()) { values <- unlist(par_list) priors <- glue::glue("constant({values})") + # determine type of parameters + bterms <- brms::brmsterms(formula) + dpars <- names(bterms$dpars) + nlpars <- names(bterms$nlpars) + # flexibly set the variables for set_prior classes <- ifelse(pars %in% dpars, "Intercept", "b") coefs <- ifelse(pars %in% dpars, "", "Intercept") diff --git a/inst/stan_chunks/sdmSimple_funs.stan b/inst/stan_chunks/sdmSimple_funs.stan index a260fac2..d3068b7e 100644 --- a/inst/stan_chunks/sdmSimple_funs.stan +++ b/inst/stan_chunks/sdmSimple_funs.stan @@ -1,3 +1,5 @@ + #include 'fun_tan_half.stan' + // utility function trick for converting real to integer type int bin_search(real x, int min_val, int max_val) { int mid_p; diff --git a/man/SDM.Rd b/man/SDM.Rd index 9a65ce4b..7f393a6a 100644 --- a/man/SDM.Rd +++ b/man/SDM.Rd @@ -52,7 +52,7 @@ and how to use it. } \item \strong{Default parameter links:} -mu = identity; c = log; kappa = log +mu = tan_half; c = log; kappa = log } } \examples{ diff --git a/tests/testthat/test-summary.R b/tests/testthat/test-summary.R index 8b08de54..7d1368ed 100644 --- a/tests/testthat/test-summary.R +++ b/tests/testthat/test-summary.R @@ -9,6 +9,6 @@ test_that("summary has reasonable outputs", { "u-95% CI", "Rhat", "Bulk_ESS", "Tail_ESS")) expect_output(print(summary1), "Constant Parameters:") expect_output(print(summary1), "Model: sdmSimple") - expect_output(print(summary1), "Links: mu = identity; c = log; kappa = log") + expect_output(print(summary1), "Links: mu = tan_half; c = log; kappa = log") expect_output(print(summary1), "Formula: mu = 0") })