Enable setting priors for linear models

Author: yoshidk6

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: BayesERtools
 Type: Package
 Title: Bayesian Exposure-Response Analysis Tools
-Version: 0.2.1
+Version: 0.2.1.1001
 Authors@R: 
     c(person("Kenta", "Yoshida", , "yoshida.kenta.6@gmail.com", role = c("aut", "cre"),
              comment = c(ORCID = "0000-0003-4967-3831")),

NAMESPACE

@@ -6,6 +6,7 @@ S3method(as_draws_df,ermod)
 S3method(as_draws_list,ermod)
 S3method(as_draws_matrix,ermod)
 S3method(as_draws_rvars,ermod)
+S3method(bayestestR::p_direction,ermod_bin)
 S3method(coef,ermod)
 S3method(extract_data,ermod)
 S3method(extract_data,ersim)
@@ -26,7 +27,6 @@ S3method(extract_var_selected,ermod_cov_sel)
 S3method(loo,ermod)
 S3method(loo,ermod_bin_emax)
 S3method(loo,ermod_emax)
-S3method(p_direction,ermod_bin)
 S3method(plot,ermod_bin)
 S3method(plot,ermod_cov_sel)
 S3method(plot,ermod_exp_sel)
@@ -41,6 +41,7 @@ S3method(print,ermod)
 S3method(print,ermod_cov_sel)
 S3method(print,ermod_exp_sel)
 S3method(print,kfold_cv_ermod)
+S3method(prior_summary,ermod)
 S3method(summary,ermod)
 export(.dev_ermod_refmodel)
 export(.select_cov_projpred)
@@ -74,7 +75,6 @@ export(extract_var_exposure)
 export(extract_var_resp)
 export(extract_var_selected)
 export(loo)
-export(p_direction)
 export(plot_coveff)
 export(plot_er)
 export(plot_er_exp_sel)
@@ -90,11 +90,11 @@ export(sim_er_curve)
 export(sim_er_curve_marg)
 export(sim_er_new_exp)
 export(sim_er_new_exp_marg)
-importFrom(bayestestR,p_direction)
 importFrom(loo,loo)
 importFrom(posterior,as_draws)
 importFrom(posterior,as_draws_array)
 importFrom(posterior,as_draws_df)
 importFrom(posterior,as_draws_list)
 importFrom(posterior,as_draws_matrix)
 importFrom(posterior,as_draws_rvars)
+importFrom(rstanarm,prior_summary)

NEWS.md

@@ -1,3 +1,10 @@
+# BayesERtools 0.2.2 (in development)
+
+## Minor changes
+
+* Enable setting the prior distribution for linear models
+* Update package dependencies
+
 # BayesERtools 0.2.1
 
 * Update package dependency

R/dev_ermod_lin.R

@@ -37,6 +37,8 @@ dev_ermod_bin <- function(
     var_resp,
     var_exposure,
     var_cov = NULL,
+    prior = rstanarm::default_prior_coef(stats::binomial()),
+    prior_intercept = rstanarm::default_prior_intercept(stats::binomial()),
     verbosity_level = 1,
     chains = 4,
     iter = 2000) {
@@ -68,6 +70,8 @@ dev_ermod_bin <- function(
     formula_final,
     family = stats::binomial(),
     data = data,
+    prior = prior,
+    prior_intercept = prior_intercept,
     QR = dplyr::if_else(length(var_full) > 1, TRUE, FALSE),
     refresh = refresh,
     chains = chains,
@@ -116,9 +120,18 @@ dev_ermod_bin_exp_sel <- function(
     data,
     var_resp,
     var_exp_candidates,
+    prior = rstanarm::default_prior_coef(stats::binomial()),
+    prior_intercept = rstanarm::default_prior_intercept(stats::binomial()),
     verbosity_level = 1,
     chains = 4,
     iter = 2000) {
+  fun_dev_ermod <-
+    purrr::partial(
+      dev_ermod_bin,
+      prior = prior,
+      prior_intercept = prior_intercept
+    )
+
   l_out <-
     .dev_ermod_exp_sel(
       data = data,
@@ -127,7 +140,7 @@ dev_ermod_bin_exp_sel <- function(
       verbosity_level = verbosity_level,
       chains = chains,
       iter = iter,
-      fun_dev_ermod = dev_ermod_bin
+      fun_dev_ermod = fun_dev_ermod
     )
 
   new_ermod_bin_exp_sel(l_out)
@@ -185,9 +198,18 @@ dev_ermod_bin_cov_sel <- function(
     validate_search = FALSE,
     nterms_max = NULL,
     .reduce_obj_size = TRUE,
+    prior = rstanarm::default_prior_coef(stats::binomial()),
+    prior_intercept = rstanarm::default_prior_intercept(stats::binomial()),
     verbosity_level = 1,
     chains = 4,
     iter = 2000) {
+  fun_dev_ermod <-
+    purrr::partial(
+      dev_ermod_bin,
+      prior = prior,
+      prior_intercept = prior_intercept
+    )
+
   ll <- .dev_ermod_cov_sel(
     data = data,
     var_resp = var_resp,
@@ -201,8 +223,10 @@ dev_ermod_bin_cov_sel <- function(
     verbosity_level = verbosity_level,
     chains = chains,
     iter = iter,
-    fun_dev_ermod = dev_ermod_bin,
-    fun_family = quote(stats::binomial())
+    fun_dev_ermod = fun_dev_ermod,
+    fun_family = quote(stats::binomial()),
+    prior = prior,
+    prior_intercept = prior_intercept
   )
 
   with(ll, new_ermod_bin_cov_sel(
@@ -243,6 +267,9 @@ dev_ermod_lin <- function(
     var_resp,
     var_exposure,
     var_cov = NULL,
+    prior = rstanarm::default_prior_coef(stats::binomial()),
+    prior_intercept = rstanarm::default_prior_intercept(stats::binomial()),
+    prior_aux = rstanarm::exponential(autoscale = TRUE),
     verbosity_level = 1,
     chains = 4,
     iter = 2000) {
@@ -272,6 +299,9 @@ dev_ermod_lin <- function(
     formula_final,
     family = stats::gaussian(),
     data = data,
+    prior = prior,
+    prior_intercept = prior_intercept,
+    prior_aux = prior_aux,
     QR = dplyr::if_else(length(var_full) > 1, TRUE, FALSE),
     refresh = refresh,
     chains = chains,
@@ -308,9 +338,20 @@ dev_ermod_lin_exp_sel <- function(
     data,
     var_resp,
     var_exp_candidates,
+    prior = rstanarm::default_prior_coef(stats::binomial()),
+    prior_intercept = rstanarm::default_prior_intercept(stats::binomial()),
+    prior_aux = rstanarm::exponential(autoscale = TRUE),
     verbosity_level = 1,
     chains = 4,
     iter = 2000) {
+  fun_dev_ermod <-
+    purrr::partial(
+      dev_ermod_lin,
+      prior = prior,
+      prior_intercept = prior_intercept,
+      prior_aux = prior_aux
+    )
+
   l_out <-
     .dev_ermod_exp_sel(
       data = data,
@@ -319,7 +360,7 @@ dev_ermod_lin_exp_sel <- function(
       verbosity_level = verbosity_level,
       chains = chains,
       iter = iter,
-      fun_dev_ermod = dev_ermod_lin
+      fun_dev_ermod = fun_dev_ermod
     )
 
   new_ermod_lin_exp_sel(l_out)
@@ -352,9 +393,20 @@ dev_ermod_lin_cov_sel <- function(
     validate_search = FALSE,
     nterms_max = NULL,
     .reduce_obj_size = TRUE,
+    prior = rstanarm::default_prior_coef(stats::binomial()),
+    prior_intercept = rstanarm::default_prior_intercept(stats::binomial()),
+    prior_aux = rstanarm::exponential(autoscale = TRUE),
     verbosity_level = 1,
     chains = 4,
     iter = 2000) {
+  fun_dev_ermod <-
+    purrr::partial(
+      dev_ermod_lin,
+      prior = prior,
+      prior_intercept = prior_intercept,
+      prior_aux = prior_aux
+    )
+
   ll <- .dev_ermod_cov_sel(
     data = data,
     var_resp = var_resp,
@@ -368,8 +420,11 @@ dev_ermod_lin_cov_sel <- function(
     verbosity_level = verbosity_level,
     chains = chains,
     iter = iter,
-    fun_dev_ermod = dev_ermod_lin,
-    fun_family = quote(stats::gaussian())
+    fun_dev_ermod = fun_dev_ermod,
+    fun_family = quote(stats::gaussian()),
+    prior = prior,
+    prior_intercept = prior_intercept,
+    prior_aux = prior_aux
   )
 
   with(ll, new_ermod_lin_cov_sel(
@@ -468,7 +523,10 @@ dev_ermod_lin_cov_sel <- function(
     chains = 4,
     iter = 2000,
     fun_dev_ermod,
-    fun_family) {
+    fun_family,
+    prior = rstanarm::default_prior_coef(stats::binomial()),
+    prior_intercept = rstanarm::default_prior_intercept(stats::binomial()),
+    prior_aux = rstanarm::exponential(autoscale = TRUE)) {
   stopifnot(verbosity_level %in% c(0, 1, 2, 3))
 
   rlang::check_installed("projpred")
@@ -492,7 +550,10 @@ dev_ermod_lin_cov_sel <- function(
     var_cov_candidates = var_cov_candidates,
     verbosity_level = verbosity_level,
     chains = chains, iter = iter,
-    fun_family = fun_family
+    fun_family = fun_family,
+    prior = prior,
+    prior_intercept = prior_intercept,
+    prior_aux = prior_aux
   )
 
   if (verbosity_level >= 1) cli::cli_h2("Step 2: Variable selection")
@@ -567,7 +628,10 @@ NULL
 .dev_ermod_refmodel <- function(
     data, var_resp, var_exposure, var_cov_candidates,
     verbosity_level = 1, chains = 4, iter = 2000,
-    fun_family = quote(stats::binomial())) {
+    fun_family = quote(stats::binomial()),
+    prior = rstanarm::default_prior_coef(stats::binomial()),
+    prior_intercept = rstanarm::default_prior_intercept(stats::binomial()),
+    prior_aux = rstanarm::exponential(autoscale = TRUE)) {
   stopifnot(verbosity_level %in% c(0, 1, 2, 3))
   refresh <- dplyr::if_else(verbosity_level >= 3, iter %/% 4, 0)
 
@@ -592,7 +656,9 @@ NULL
     rlang::call2(rstanarm::stan_glm,
       formula = formula_full,
       family = fun_family, data = quote(data), QR = TRUE,
-      refresh = refresh, chains = chains, iter = iter
+      refresh = refresh, chains = chains, iter = iter,
+      prior = prior, prior_intercept = prior_intercept,
+      prior_aux = prior_aux
     )
   fit_ref <- eval(call_fit_ref)
 

R/ermod-methods.R

@@ -55,8 +55,8 @@ plot.ermod_bin <- function(x, show_orig_data = FALSE, ...) {
 #' @export
 #' @rdname ermod_method
 coef.ermod <- function(object, ...) {
-  if (inherits(object, "ermod_emax")) {
-    stop("coef() not supported for ermod_emax object")
+  if (!inherits(object, c("ermod_bin", "ermod_lin"))) {
+    stop("coef() only supported for linear models")
   }
 
   stats::coef(object$mod, ...)
@@ -65,8 +65,8 @@ coef.ermod <- function(object, ...) {
 #' @export
 #' @rdname ermod_method
 summary.ermod <- function(object, ...) {
-  if (inherits(object, "ermod_emax")) {
-    stop("summary() not supported for ermod_emax object")
+  if (!inherits(object, c("ermod_bin", "ermod_lin"))) {
+    stop("summary() only supported for linear models")
   }
 
   summary(object$mod, ...)
@@ -223,7 +223,7 @@ extract_var_selected.ermod_cov_sel <- function(x) x$var_selected
 #'  credible interval (.lower, .upper)
 #'
 extract_coef_exp_ci <- function(x, ci_width = 0.95) {
-  # Check that input x is ermod object
+  # Check that input x is linear ermod object
   if (!inherits(x, c("ermod_bin", "ermod_lin"))) {
     stop("extract_coef_exp_ci() only supported for linear models")
   }
@@ -367,3 +367,22 @@ as_draws_matrix.ermod <- function(x, ...) {
 as_draws_rvars.ermod <- function(x, ...) {
   posterior::as_draws_rvars(x$mod, ...)
 }
+
+# prior_summary --------------------------------------------------------------
+#' Summarize the priors used for linear or linear logistic regression models
+#'
+#' See [rstanarm::prior_summary()] for details.
+#'
+#' @export
+#' @rdname prior_summary
+#' @importFrom rstanarm prior_summary
+#' @return An object of class `prior_summary.stanreg`
+#'
+prior_summary.ermod <- function(object, ...) {
+  # Check that input x is linear ermod object
+  if (!inherits(object, c("ermod_bin", "ermod_lin"))) {
+    stop("prior_summary.ermod() only supported for linear models")
+  }
+
+  rstanarm::prior_summary(object$mod, ...)
+}

R/p_direction.R

@@ -80,4 +80,3 @@ p_direction.ermod_bin <- function(
 .if_run_ex_p_dir <- function() {
   requireNamespace("bayestestR", quietly = TRUE)
 }
-

R/sim_ermod.R

@@ -181,7 +181,7 @@ sim_er <- function(
 #' )
 #'
 #' ersim_new_exp_med_qi
-#'}
+#' }
 #'
 sim_er_new_exp <- function(
     ermod,