two sets of shape hyperparameters for the Beta hyperprior in the SBM (#68)

* two sets of shape hyperparameters for the Beta hyperprior in the SBM

* Update documentation  for the extra shape parameters

* update documentation

* adress review comments

---------

Co-authored-by: Nikola <[email protected]>

Author: sekulovskin

R/RcppExports.R

@@ -5,8 +5,8 @@ run_bgmCompare_parallel <- function(observations, num_groups, counts_per_categor
     .Call(`_bgms_run_bgmCompare_parallel`, observations, num_groups, counts_per_category, blume_capel_stats, pairwise_stats, num_categories, main_alpha, main_beta, pairwise_scale, difference_scale, difference_selection_alpha, difference_selection_beta, difference_prior, iter, warmup, na_impute, missing_data_indices, is_ordinal_variable, baseline_category, difference_selection, main_effect_indices, pairwise_effect_indices, target_accept, nuts_max_depth, learn_mass_matrix, projection, group_membership, group_indices, interaction_index_matrix, inclusion_probability, num_chains, nThreads, seed, update_method, hmc_num_leapfrogs, progress_type)
 }
 
-run_bgm_parallel <- function(observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type) {
-    .Call(`_bgms_run_bgm_parallel`, observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type)
+run_bgm_parallel <- function(observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type) {
+    .Call(`_bgms_run_bgm_parallel`, observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type)
 }
 
 get_explog_switch <- function() {

R/bgm.R

@@ -207,8 +207,14 @@
 #'
 #' @param beta_bernoulli_alpha,beta_bernoulli_beta Double. Shape parameters
 #'   for the beta distribution in the Beta–Bernoulli and the Stochastic-Block
-#'   priors. Must be positive. Defaults: \code{beta_bernoulli_alpha = 1} and
-#'   \code{beta_bernoulli_beta = 1}.
+#'   priors. Must be positive. For the Stochastic-Block prior these are the shape
+#'   parameters for the within-cluster edge inclusion probabilities.
+#'   Defaults: \code{beta_bernoulli_alpha = 1} and \code{beta_bernoulli_beta = 1}.
+#'
+#' @param beta_bernoulli_alpha_between,beta_bernoulli_beta_between Double.
+#' Shape parameters for the between-cluster edge inclusion probabilities in the
+#' Stochastic-Block prior. Must be positive.
+#' Default: \code{beta_bernoulli_alpha_between = 1} and \code{beta_bernoulli_beta_between = 1}
 #'
 #' @param dirichlet_alpha Double. Concentration parameter of the Dirichlet
 #'   prior on block assignments (used with the Stochastic Block model).
@@ -359,6 +365,8 @@ bgm = function(
     inclusion_probability = 0.5,
     beta_bernoulli_alpha = 1,
     beta_bernoulli_beta = 1,
+    beta_bernoulli_alpha_between = 1,
+    beta_bernoulli_beta_between = 1,
     dirichlet_alpha = 1,
     lambda = 1,
     na_action = c("listwise", "impute"),
@@ -418,7 +426,7 @@ bgm = function(
     } else if(update_method == "hamiltonian-mc") {
       target_accept = 0.65
     } else if(update_method == "nuts") {
-      target_accept = 0.80
+      target_accept = 0.60
     }
   }
 
@@ -444,9 +452,21 @@ bgm = function(
                       inclusion_probability = inclusion_probability,
                       beta_bernoulli_alpha = beta_bernoulli_alpha,
                       beta_bernoulli_beta = beta_bernoulli_beta,
+                      beta_bernoulli_alpha_between = beta_bernoulli_alpha_between,
+                      beta_bernoulli_beta_between = beta_bernoulli_beta_between,
                       dirichlet_alpha = dirichlet_alpha,
                       lambda = lambda)
 
+  # check hyperparameters input
+  # If user left them NULL, pass -1 to C++ (means: ignore between prior)
+  if (is.null(beta_bernoulli_alpha_between) && is.null(beta_bernoulli_beta_between)) {
+    beta_bernoulli_alpha_between <- -1.0
+    beta_bernoulli_beta_between  <- -1.0
+  } else if (is.null(beta_bernoulli_alpha_between) || is.null(beta_bernoulli_beta_between)) {
+    stop("If you wish to specify different between and within cluster probabilites,
+         provide both beta_bernoulli_alpha_between and beta_bernoulli_beta_between,
+         otherwise leave both NULL.")
+  }
   # ----------------------------------------------------------------------------
   # The vector variable_type is now coded as boolean.
   # Ordinal (variable_bool == TRUE) or Blume-Capel (variable_bool == FALSE)
@@ -572,6 +592,8 @@ bgm = function(
     inclusion_probability = inclusion_probability,
     beta_bernoulli_alpha = beta_bernoulli_alpha,
     beta_bernoulli_beta = beta_bernoulli_beta,
+    beta_bernoulli_alpha_between = beta_bernoulli_alpha_between,
+    beta_bernoulli_beta_between = beta_bernoulli_beta_between,
     dirichlet_alpha = dirichlet_alpha, lambda = lambda,
     interaction_index_matrix = interaction_index_matrix, iter = iter,
     warmup = warmup, counts_per_category = counts_per_category,
@@ -603,6 +625,7 @@ bgm = function(
         na_action = na_action, na_impute = na_impute,
         edge_selection = edge_selection, edge_prior = edge_prior, inclusion_probability = inclusion_probability,
         beta_bernoulli_alpha = beta_bernoulli_alpha, beta_bernoulli_beta = beta_bernoulli_beta,
+        beta_bernoulli_alpha_between = beta_bernoulli_alpha_between, beta_bernoulli_beta_between = beta_bernoulli_beta_between,
         dirichlet_alpha = dirichlet_alpha, lambda = lambda,
         variable_type = variable_type,
         update_method = update_method,
@@ -634,6 +657,8 @@ bgm = function(
     edge_selection = edge_selection, edge_prior = edge_prior, inclusion_probability = inclusion_probability,
     beta_bernoulli_alpha = beta_bernoulli_alpha,
     beta_bernoulli_beta = beta_bernoulli_beta,
+    beta_bernoulli_alpha_between = beta_bernoulli_alpha_between,
+    beta_bernoulli_beta_between = beta_bernoulli_beta_between,
     dirichlet_alpha = dirichlet_alpha, lambda = lambda,
     variable_type = variable_type,
     update_method = update_method,

R/function_input_utils.R

@@ -31,6 +31,8 @@ check_model = function(x,
                        inclusion_probability = 0.5,
                        beta_bernoulli_alpha = 1,
                        beta_bernoulli_beta = 1,
+                       beta_bernoulli_alpha_between = 1,
+                       beta_bernoulli_beta_between = 1,
                        dirichlet_alpha = dirichlet_alpha,
                        lambda = lambda) {
 
@@ -204,18 +206,42 @@ check_model = function(x,
          is.null(beta_bernoulli_alpha) || is.null(beta_bernoulli_beta))
         stop("Values for both scale parameters of the beta distribution need to be specified.")
     }
+
     if(edge_prior == "Stochastic-Block") {
       theta = matrix(0.5, nrow = ncol(x), ncol = ncol(x))
-      if(beta_bernoulli_alpha <= 0 || beta_bernoulli_beta <= 0 || dirichlet_alpha <= 0 || lambda <= 0)
-        stop("The scale parameters of the beta and Dirichlet distribution need to be positive.")
-      if(!is.finite(beta_bernoulli_alpha) || !is.finite(beta_bernoulli_beta) || !is.finite(dirichlet_alpha) || !is.finite(lambda))
-        stop("The scale parameters of the beta distribution, the concentration parameter of the Dirichlet distribution, and the rate parameter of the Poisson distribution need to be finite.")
-      if(is.na(beta_bernoulli_alpha) || is.na(beta_bernoulli_beta) ||
-         is.null(beta_bernoulli_alpha) || is.null(beta_bernoulli_beta) ||
-         is.null(dirichlet_alpha) || is.null(dirichlet_alpha) || is.null(lambda) || is.null(lambda))
-        stop("Values for both scale parameters of the beta distribution, the concentration parameter of the Dirichlet distribution, and the rate parameter of the Poisson distribution need to be specified.")
+
+      # Check that all beta parameters are provided
+      if (is.null(beta_bernoulli_alpha) || is.null(beta_bernoulli_beta) ||
+          is.null(beta_bernoulli_alpha_between) || is.null(beta_bernoulli_beta_between)) {
+        stop("The Stochastic-Block prior requires all four beta parameters: ",
+             "beta_bernoulli_alpha, beta_bernoulli_beta, ",
+             "beta_bernoulli_alpha_between, and beta_bernoulli_beta_between.")
+      }
+
+      # Check that all beta parameters are positive
+      if (beta_bernoulli_alpha <= 0 || beta_bernoulli_beta <= 0 ||
+          beta_bernoulli_alpha_between <= 0 || beta_bernoulli_beta_between <= 0 ||
+          dirichlet_alpha <= 0 || lambda <= 0) {
+        stop("The parameters of the beta and Dirichlet distributions need to be positive.")
+      }
+
+      # Check that all beta parameters are finite
+      if (!is.finite(beta_bernoulli_alpha) || !is.finite(beta_bernoulli_beta) ||
+          !is.finite(beta_bernoulli_alpha_between) || !is.finite(beta_bernoulli_beta_between) ||
+          !is.finite(dirichlet_alpha) || !is.finite(lambda)) {
+        stop("The shape parameters of the beta distribution, the concentration parameter of the Dirichlet distribution, ",
+             "and the rate parameter of the Poisson distribution need to be finite.")
+      }
+
+      # Check for NAs
+      if (is.na(beta_bernoulli_alpha) || is.na(beta_bernoulli_beta) ||
+          is.na(beta_bernoulli_alpha_between) || is.na(beta_bernoulli_beta_between) ||
+          is.na(dirichlet_alpha) || is.na(lambda)) {
+        stop("Values for all shape parameters of the beta distribution, the concentration parameter of the Dirichlet distribution, ",
+             "and the rate parameter of the Poisson distribution cannot be NA.")
+      }
     }
-  } else {
+ }else {
     theta = matrix(0.5, nrow = 1, ncol = 1)
     edge_prior = "Not Applicable"
   }

R/output_utils.R

@@ -2,7 +2,8 @@ prepare_output_bgm = function(
     out, x, num_categories, iter, data_columnnames, is_ordinal_variable,
     warmup, pairwise_scale, main_alpha, main_beta,
     na_action, na_impute, edge_selection, edge_prior, inclusion_probability,
-    beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda,
+    beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between,
+    beta_bernoulli_beta_between,dirichlet_alpha, lambda,
     variable_type, update_method, target_accept, hmc_num_leapfrogs,
     nuts_max_depth, learn_mass_matrix, num_chains
 ) {
@@ -22,6 +23,8 @@ prepare_output_bgm = function(
     inclusion_probability = inclusion_probability,
     beta_bernoulli_alpha = beta_bernoulli_alpha,
     beta_bernoulli_beta = beta_bernoulli_beta,
+    beta_bernoulli_alpha_between = beta_bernoulli_alpha_between,
+    beta_bernoulli_beta_between = beta_bernoulli_beta_between,
     dirichlet_alpha = dirichlet_alpha,
     lambda = lambda,
     na_action = na_action,

man/bgm.Rd

@@ -58,8 +58,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // run_bgm_parallel
-Rcpp::List run_bgm_parallel(const arma::imat& observations, const arma::ivec& num_categories, double pairwise_scale, const std::string& edge_prior, const arma::mat& inclusion_probability, double beta_bernoulli_alpha, double beta_bernoulli_beta, double dirichlet_alpha, double lambda, const arma::imat& interaction_index_matrix, int iter, int warmup, const arma::imat& counts_per_category, const arma::imat& blume_capel_stats, double main_alpha, double main_beta, bool na_impute, const arma::imat& missing_index, const arma::uvec& is_ordinal_variable, const arma::ivec& baseline_category, bool edge_selection, const std::string& update_method, const arma::imat& pairwise_effect_indices, double target_accept, const arma::imat& pairwise_stats, int hmc_num_leapfrogs, int nuts_max_depth, bool learn_mass_matrix, int num_chains, int nThreads, int seed, int progress_type);
-RcppExport SEXP _bgms_run_bgm_parallel(SEXP observationsSEXP, SEXP num_categoriesSEXP, SEXP pairwise_scaleSEXP, SEXP edge_priorSEXP, SEXP inclusion_probabilitySEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP interaction_index_matrixSEXP, SEXP iterSEXP, SEXP warmupSEXP, SEXP counts_per_categorySEXP, SEXP blume_capel_statsSEXP, SEXP main_alphaSEXP, SEXP main_betaSEXP, SEXP na_imputeSEXP, SEXP missing_indexSEXP, SEXP is_ordinal_variableSEXP, SEXP baseline_categorySEXP, SEXP edge_selectionSEXP, SEXP update_methodSEXP, SEXP pairwise_effect_indicesSEXP, SEXP target_acceptSEXP, SEXP pairwise_statsSEXP, SEXP hmc_num_leapfrogsSEXP, SEXP nuts_max_depthSEXP, SEXP learn_mass_matrixSEXP, SEXP num_chainsSEXP, SEXP nThreadsSEXP, SEXP seedSEXP, SEXP progress_typeSEXP) {
+Rcpp::List run_bgm_parallel(const arma::imat& observations, const arma::ivec& num_categories, double pairwise_scale, const std::string& edge_prior, const arma::mat& inclusion_probability, double beta_bernoulli_alpha, double beta_bernoulli_beta, double beta_bernoulli_alpha_between, double beta_bernoulli_beta_between, double dirichlet_alpha, double lambda, const arma::imat& interaction_index_matrix, int iter, int warmup, const arma::imat& counts_per_category, const arma::imat& blume_capel_stats, double main_alpha, double main_beta, bool na_impute, const arma::imat& missing_index, const arma::uvec& is_ordinal_variable, const arma::ivec& baseline_category, bool edge_selection, const std::string& update_method, const arma::imat& pairwise_effect_indices, double target_accept, const arma::imat& pairwise_stats, int hmc_num_leapfrogs, int nuts_max_depth, bool learn_mass_matrix, int num_chains, int nThreads, int seed, int progress_type);
+RcppExport SEXP _bgms_run_bgm_parallel(SEXP observationsSEXP, SEXP num_categoriesSEXP, SEXP pairwise_scaleSEXP, SEXP edge_priorSEXP, SEXP inclusion_probabilitySEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP beta_bernoulli_alpha_betweenSEXP, SEXP beta_bernoulli_beta_betweenSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP interaction_index_matrixSEXP, SEXP iterSEXP, SEXP warmupSEXP, SEXP counts_per_categorySEXP, SEXP blume_capel_statsSEXP, SEXP main_alphaSEXP, SEXP main_betaSEXP, SEXP na_imputeSEXP, SEXP missing_indexSEXP, SEXP is_ordinal_variableSEXP, SEXP baseline_categorySEXP, SEXP edge_selectionSEXP, SEXP update_methodSEXP, SEXP pairwise_effect_indicesSEXP, SEXP target_acceptSEXP, SEXP pairwise_statsSEXP, SEXP hmc_num_leapfrogsSEXP, SEXP nuts_max_depthSEXP, SEXP learn_mass_matrixSEXP, SEXP num_chainsSEXP, SEXP nThreadsSEXP, SEXP seedSEXP, SEXP progress_typeSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -70,6 +70,8 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< const arma::mat& >::type inclusion_probability(inclusion_probabilitySEXP);
     Rcpp::traits::input_parameter< double >::type beta_bernoulli_alpha(beta_bernoulli_alphaSEXP);
     Rcpp::traits::input_parameter< double >::type beta_bernoulli_beta(beta_bernoulli_betaSEXP);
+    Rcpp::traits::input_parameter< double >::type beta_bernoulli_alpha_between(beta_bernoulli_alpha_betweenSEXP);
+    Rcpp::traits::input_parameter< double >::type beta_bernoulli_beta_between(beta_bernoulli_beta_betweenSEXP);
     Rcpp::traits::input_parameter< double >::type dirichlet_alpha(dirichlet_alphaSEXP);
     Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
     Rcpp::traits::input_parameter< const arma::imat& >::type interaction_index_matrix(interaction_index_matrixSEXP);
@@ -95,7 +97,7 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< int >::type nThreads(nThreadsSEXP);
     Rcpp::traits::input_parameter< int >::type seed(seedSEXP);
     Rcpp::traits::input_parameter< int >::type progress_type(progress_typeSEXP);
-    rcpp_result_gen = Rcpp::wrap(run_bgm_parallel(observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type));
+    rcpp_result_gen = Rcpp::wrap(run_bgm_parallel(observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -182,7 +184,7 @@ END_RCPP
 
 static const R_CallMethodDef CallEntries[] = {
     {"_bgms_run_bgmCompare_parallel", (DL_FUNC) &_bgms_run_bgmCompare_parallel, 36},
-    {"_bgms_run_bgm_parallel", (DL_FUNC) &_bgms_run_bgm_parallel, 32},
+    {"_bgms_run_bgm_parallel", (DL_FUNC) &_bgms_run_bgm_parallel, 34},
     {"_bgms_get_explog_switch", (DL_FUNC) &_bgms_get_explog_switch, 0},
     {"_bgms_rcpp_ieee754_exp", (DL_FUNC) &_bgms_rcpp_ieee754_exp, 1},
     {"_bgms_rcpp_ieee754_log", (DL_FUNC) &_bgms_rcpp_ieee754_log, 1},