Implement adaptive mala and fisher mala for thresholds and interactions.

Author: MaartenMarsman

R/RcppExports.R

@@ -9,8 +9,8 @@ sample_bcomrf_gibbs <- function(no_states, no_variables, no_categories, interact
     .Call(`_bgms_sample_bcomrf_gibbs`, no_states, no_variables, no_categories, interactions, thresholds, variable_type, reference_category, iter)
 }
 
-run_gibbs_sampler_for_bgm <- function(observations, num_categories, interaction_scale, edge_prior, theta, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, Index, iter, burnin, num_obs_categories, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, is_ordinal_variable, reference_category, save_main = FALSE, save_pairwise = FALSE, save_indicator = FALSE, display_progress = FALSE, edge_selection = TRUE, use_mala = FALSE) {
-    .Call(`_bgms_run_gibbs_sampler_for_bgm`, observations, num_categories, interaction_scale, edge_prior, theta, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, Index, iter, burnin, num_obs_categories, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, is_ordinal_variable, reference_category, save_main, save_pairwise, save_indicator, display_progress, edge_selection, use_mala)
+run_gibbs_sampler_for_bgm <- function(observations, num_categories, interaction_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, burnin, num_obs_categories, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, is_ordinal_variable, reference_category, save_main = FALSE, save_pairwise = FALSE, save_indicator = FALSE, display_progress = FALSE, edge_selection = TRUE, update_method = "adaptive-metropolis") {
+    .Call(`_bgms_run_gibbs_sampler_for_bgm`, observations, num_categories, interaction_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, burnin, num_obs_categories, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, is_ordinal_variable, reference_category, save_main, save_pairwise, save_indicator, display_progress, edge_selection, update_method)
 }
 
 compare_anova_gibbs_sampler <- function(observations, main_effect_indices, pairwise_effect_indices, projection, num_categories, num_groups, group_indices, interaction_scale, pairwise_difference_scale, main_difference_scale, pairwise_difference_prior, main_difference_prior, inclusion_probability_difference, pairwise_beta_bernoulli_alpha, pairwise_beta_bernoulli_beta, main_beta_bernoulli_alpha, main_beta_bernoulli_beta, Index, iter, burnin, num_obs_categories, sufficient_blume_capel, prior_threshold_alpha, prior_threshold_beta, na_impute, missing_data_indices, is_ordinal_variable, baseline_category, independent_thresholds, save_main = FALSE, save_pairwise = FALSE, save_indicator = FALSE, display_progress = FALSE, difference_selection = TRUE) {

R/bgm.R

@@ -152,6 +152,14 @@
 #' @param display_progress Should the function show a progress bar
 #' (\code{display_progress = TRUE})? Or not (\code{display_progress = FALSE})?
 #' The default is \code{TRUE}.
+#' @param update_method Character. Specifies how the MCMC sampler updates the threshold
+#' and interaction parameters:
+#' \describe{
+#'   \item{"adaptive-metropolis"}{Uses adaptive Metropolis-Hastings for both thresholds and interactions.}
+#'   \item{"adaptive-mala"}{Uses Fisher-preconditioned MALA for thresholds and standard MALA for interactions and indicators.}
+#'   \item{"fisher-mala"}{Uses Fisher-preconditioned MALA for both thresholds and interactions.}
+#' }
+#' Defaults to \code{"adaptive-metropolis"}.
 #'
 #' @return If \code{save = FALSE} (the default), the result is a list of class
 #' ``bgms'' containing the following matrices with model-averaged quantities:
@@ -317,7 +325,8 @@ bgm = function(x,
                save_pairwise = FALSE,
                save_indicator = FALSE,
                display_progress = TRUE,
-               mala = FALSE) {
+               update_method = c("adaptive-metropolis", "adaptive-mala", "fisher-mala")
+) {
 
   # Deprecation warning for save parameter
   if(hasArg(save)) {
@@ -333,6 +342,11 @@ bgm = function(x,
   save_indicator = check_logical(save_indicator, "save_indicator")
 
 
+  # Check update method
+  update_method_input = update_method
+  update_method = match.arg(update_method)
+
+
   #Check data input ------------------------------------------------------------
   if(!inherits(x, what = "matrix") && !inherits(x, what = "data.frame"))
     stop("The input x needs to be a matrix or dataframe.")
@@ -368,7 +382,7 @@ bgm = function(x,
   reference_category = model$reference_category
   edge_selection = model$edge_selection
   edge_prior = model$edge_prior
-  theta = model$theta
+  inclusion_probability = model$inclusion_probability
 
   #Check Gibbs input -----------------------------------------------------------
   if(abs(iter - round(iter)) > .Machine$double.eps)
@@ -446,27 +460,28 @@ bgm = function(x,
     }
   }
 
-  # Index vector used to sample interactions in a random order -----------------
-  Index = matrix(0,
+  # Index matrix used in the c++ functions  ------------------------------------
+  interaction_index_matrix = matrix(0,
                  nrow = num_variables * (num_variables - 1) / 2,
                  ncol = 3)
   cntr = 0
   for(variable1 in 1:(num_variables - 1)) {
     for(variable2 in (variable1 + 1):num_variables) {
       cntr =  cntr + 1
-      Index[cntr, 1] = cntr
-      Index[cntr, 2] = variable1 - 1
-      Index[cntr, 3] = variable2 - 1
+      interaction_index_matrix[cntr, 1] = cntr
+      interaction_index_matrix[cntr, 2] = variable1 - 1
+      interaction_index_matrix[cntr, 3] = variable2 - 1
     }
   }
 
   # Call the Rcpp function
   out = run_gibbs_sampler_for_bgm (
     observations = x, num_categories = num_categories,
     interaction_scale = interaction_scale, edge_prior = edge_prior,
-    theta = theta, beta_bernoulli_alpha = beta_bernoulli_alpha,
+    inclusion_probability = inclusion_probability, beta_bernoulli_alpha = beta_bernoulli_alpha,
     beta_bernoulli_beta = beta_bernoulli_beta,
-    dirichlet_alpha = dirichlet_alpha, lambda = lambda, Index = Index,
+    dirichlet_alpha = dirichlet_alpha, lambda = lambda,
+    interaction_index_matrix = interaction_index_matrix,
     iter = iter, burnin = burnin, num_obs_categories = num_obs_categories,
     sufficient_blume_capel = sufficient_blume_capel,
     threshold_alpha = threshold_alpha, threshold_beta = threshold_beta,
@@ -475,7 +490,7 @@ bgm = function(x,
     reference_category = reference_category, save_main = save_main,
     save_pairwise = save_pairwise, save_indicator = save_indicator,
     display_progress = display_progress, edge_selection = edge_selection,
-    use_mala = mala
+    update_method = update_method
   )
 
   # Main output handler in the wrapper function
@@ -488,7 +503,7 @@ bgm = function(x,
     burnin = burnin, interaction_scale = interaction_scale,
     threshold_alpha = threshold_alpha, threshold_beta = threshold_beta,
     na_action = na_action, na_impute = na_impute,
-    edge_selection = edge_selection, edge_prior = edge_prior, theta = theta,
+    edge_selection = edge_selection, edge_prior = edge_prior, inclusion_probability = inclusion_probability,
     beta_bernoulli_alpha = beta_bernoulli_alpha,
     beta_bernoulli_beta = beta_bernoulli_beta,
     dirichlet_alpha = dirichlet_alpha, lambda = lambda,

R/function_input_utils.R

@@ -226,7 +226,7 @@ check_model = function(x,
               reference_category = reference_category,
               edge_selection = edge_selection,
               edge_prior = edge_prior,
-              theta = theta))
+              inclusion_probability = theta))
 }
 
 check_compare_model = function(x,

R/output_utils.R

@@ -3,7 +3,7 @@
 prepare_output_bgm = function (
     out, x, num_categories, iter, data_columnnames, is_ordinal_variable,
     save_options, burnin, interaction_scale, threshold_alpha, threshold_beta,
-    na_action, na_impute, edge_selection, edge_prior, theta,
+    na_action, na_impute, edge_selection, edge_prior, inclusion_probability,
     beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda,
     variable_type) {
 
@@ -15,7 +15,7 @@ prepare_output_bgm = function (
     burnin = burnin, interaction_scale = interaction_scale,
     threshold_alpha = threshold_alpha, threshold_beta = threshold_beta,
     edge_selection = edge_selection, edge_prior = edge_prior,
-    inclusion_probability = theta, beta_bernoulli_alpha = beta_bernoulli_alpha,
+    inclusion_probability = inclusion_probability, beta_bernoulli_alpha = beta_bernoulli_alpha,
     beta_bernoulli_beta =  beta_bernoulli_beta,
     dirichlet_alpha = dirichlet_alpha, lambda = lambda, na_action = na_action,
     save = save, version = packageVersion("bgms")

src/RcppExports.cpp

@@ -46,21 +46,21 @@ BEGIN_RCPP
 END_RCPP
 }
 // run_gibbs_sampler_for_bgm
-List run_gibbs_sampler_for_bgm(arma::imat& observations, const arma::ivec& num_categories, const double interaction_scale, const String& edge_prior, arma::mat& theta, const double beta_bernoulli_alpha, const double beta_bernoulli_beta, const double dirichlet_alpha, const double lambda, const arma::imat& Index, const int iter, const int burnin, arma::imat& num_obs_categories, arma::imat& sufficient_blume_capel, const double threshold_alpha, const double threshold_beta, const bool na_impute, const arma::imat& missing_index, const arma::uvec& is_ordinal_variable, const arma::ivec& reference_category, const bool save_main, const bool save_pairwise, const bool save_indicator, const bool display_progress, bool edge_selection, bool use_mala);
-RcppExport SEXP _bgms_run_gibbs_sampler_for_bgm(SEXP observationsSEXP, SEXP num_categoriesSEXP, SEXP interaction_scaleSEXP, SEXP edge_priorSEXP, SEXP thetaSEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP IndexSEXP, SEXP iterSEXP, SEXP burninSEXP, SEXP num_obs_categoriesSEXP, SEXP sufficient_blume_capelSEXP, SEXP threshold_alphaSEXP, SEXP threshold_betaSEXP, SEXP na_imputeSEXP, SEXP missing_indexSEXP, SEXP is_ordinal_variableSEXP, SEXP reference_categorySEXP, SEXP save_mainSEXP, SEXP save_pairwiseSEXP, SEXP save_indicatorSEXP, SEXP display_progressSEXP, SEXP edge_selectionSEXP, SEXP use_malaSEXP) {
+List run_gibbs_sampler_for_bgm(arma::imat& observations, const arma::ivec& num_categories, const double interaction_scale, const String& edge_prior, arma::mat& inclusion_probability, const double beta_bernoulli_alpha, const double beta_bernoulli_beta, const double dirichlet_alpha, const double lambda, const arma::imat& interaction_index_matrix, const int iter, const int burnin, arma::imat& num_obs_categories, arma::imat& sufficient_blume_capel, const double threshold_alpha, const double threshold_beta, const bool na_impute, const arma::imat& missing_index, const arma::uvec& is_ordinal_variable, const arma::ivec& reference_category, const bool save_main, const bool save_pairwise, const bool save_indicator, const bool display_progress, bool edge_selection, const std::string& update_method);
+RcppExport SEXP _bgms_run_gibbs_sampler_for_bgm(SEXP observationsSEXP, SEXP num_categoriesSEXP, SEXP interaction_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 burninSEXP, SEXP num_obs_categoriesSEXP, SEXP sufficient_blume_capelSEXP, SEXP threshold_alphaSEXP, SEXP threshold_betaSEXP, SEXP na_imputeSEXP, SEXP missing_indexSEXP, SEXP is_ordinal_variableSEXP, SEXP reference_categorySEXP, SEXP save_mainSEXP, SEXP save_pairwiseSEXP, SEXP save_indicatorSEXP, SEXP display_progressSEXP, SEXP edge_selectionSEXP, SEXP update_methodSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< arma::imat& >::type observations(observationsSEXP);
     Rcpp::traits::input_parameter< const arma::ivec& >::type num_categories(num_categoriesSEXP);
     Rcpp::traits::input_parameter< const double >::type interaction_scale(interaction_scaleSEXP);
     Rcpp::traits::input_parameter< const String& >::type edge_prior(edge_priorSEXP);
-    Rcpp::traits::input_parameter< arma::mat& >::type theta(thetaSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type inclusion_probability(inclusion_probabilitySEXP);
     Rcpp::traits::input_parameter< const double >::type beta_bernoulli_alpha(beta_bernoulli_alphaSEXP);
     Rcpp::traits::input_parameter< const double >::type beta_bernoulli_beta(beta_bernoulli_betaSEXP);
     Rcpp::traits::input_parameter< const double >::type dirichlet_alpha(dirichlet_alphaSEXP);
     Rcpp::traits::input_parameter< const double >::type lambda(lambdaSEXP);
-    Rcpp::traits::input_parameter< const arma::imat& >::type Index(IndexSEXP);
+    Rcpp::traits::input_parameter< const arma::imat& >::type interaction_index_matrix(interaction_index_matrixSEXP);
     Rcpp::traits::input_parameter< const int >::type iter(iterSEXP);
     Rcpp::traits::input_parameter< const int >::type burnin(burninSEXP);
     Rcpp::traits::input_parameter< arma::imat& >::type num_obs_categories(num_obs_categoriesSEXP);
@@ -76,8 +76,8 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< const bool >::type save_indicator(save_indicatorSEXP);
     Rcpp::traits::input_parameter< const bool >::type display_progress(display_progressSEXP);
     Rcpp::traits::input_parameter< bool >::type edge_selection(edge_selectionSEXP);
-    Rcpp::traits::input_parameter< bool >::type use_mala(use_malaSEXP);
-    rcpp_result_gen = Rcpp::wrap(run_gibbs_sampler_for_bgm(observations, num_categories, interaction_scale, edge_prior, theta, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, Index, iter, burnin, num_obs_categories, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, is_ordinal_variable, reference_category, save_main, save_pairwise, save_indicator, display_progress, edge_selection, use_mala));
+    Rcpp::traits::input_parameter< const std::string& >::type update_method(update_methodSEXP);
+    rcpp_result_gen = Rcpp::wrap(run_gibbs_sampler_for_bgm(observations, num_categories, interaction_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, burnin, num_obs_categories, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, is_ordinal_variable, reference_category, save_main, save_pairwise, save_indicator, display_progress, edge_selection, update_method));
     return rcpp_result_gen;
 END_RCPP
 }