BCFM: Bayesian Clustering Factor Models

Description

Implements the Bayesian Clustering Factor Models (BCFM) for simultaneous clustering and latent factor analysis of multivariate longitudinal data. The model accounts for within-cluster dependence through shared latent factors while allowing heterogeneity across clusters, enabling flexible covariance modeling in high-dimensional settings. Inference is performed using Markov chain Monte Carlo (MCMC) methods with computationally intensive steps implemented via 'Rcpp'. Model selection and visualization tools are provided. The methodology is described in Shin, Ferreira, and Tegge (2018) doi:10.1002/sim.70350.

Author(s)

Maintainer: Meriem Touami meriemt@vt.edu

Authors:

• Allison Tegge ategge@vt.edu

• Marco Ferreira marf@vt.edu

• Hwasoo Shin hshin2@hfhs.org

See Also

Useful links:

• https://github.com/ategge/BCFM

• Report bugs at https://github.com/ategge/BCFM/issues

Fit BCFM Model

Description

Fits a single Bayesian Clustering Factor Models (BCFM) using C++ implementation. This function serves as a wrapper for the BCFMcpp function, handling timing and output formatting.

Usage

BCFM.fit(
  data,
  model.attributes,
  hyp.parm,
  n.iter = 50000,
  vague.mu = FALSE,
  covariance = TRUE,
  p.exponent = 2,
  every = 10
)

Arguments

Value

A list containing:

Result

Output from the BCFMcpp C++ function

model.attributes

The model attributes used

hyp.parm

The hyperparameters used

start.time

POSIXct timestamp when model started

end.time

POSIXct timestamp when model completed

run.time

Total elapsed time for model execution

See Also

BCFM.model.selection for fitting multiple models, initialize.model.attributes, initialize.hyp.parm

Examples

# Prepare data using the included simulated dataset
data(sim.data)
data.pre <- init.data(sim.data, paste0("V", 1:5))

# Initialize model components
model.attributes <- initialize.model.attributes(S = nrow(sim.data), times = 1,
                                                 R = 5, L = 2, G = 2)
cluster.hyperparms <- initialize.cluster.hyperparms(data.pre, model.attributes)
hyp.parm <- initialize.hyp.parm(model.attributes, cluster.hyperparms)

# Fit model
result <- BCFM.fit(data.pre, model.attributes, hyp.parm,
                   n.iter = 100, every = 10)
result$run.time

BCFM Model Selection Over Multiple Groups and Factors

Description

Performs Bayesian Clustering Factor Models analysis across a grid of group numbers and factor numbers. For each combination, the function fits the BCFM model, calculates IC, and saves results. This is the primary function for model selection to determine the optimal number of clusters and latent factors.

Usage

BCFM.model.selection(
  data,
  cluster.vars,
  grouplist,
  factorlist,
  n.iter = 50000,
  vague.mu = FALSE,
  covariance = TRUE,
  p.exponent = 2,
  every = 10,
  cluster.size = 0.05,
  burnin = NA,
  output_dir = tempdir(),
  seed = NULL
)

Arguments

Details

The function performs the following steps for each group-factor combination:

1. Preprocesses data using init.data

2. Determines optimal variable ordering using permutation.order

3. Initializes model attributes and hyperparameters

4. Fits BCFM model using BCFM.fit

5. Calculates IC for model comparison

6. Saves individual results and cumulative IC matrix

The IC matrix can be used to identify the optimal model configuration by selecting the combination of groups and factors with the lowest IC value.

Value

Invisibly returns NULL. Results are saved to disk:

results-covarianceF-gX-fY.Rdata

Individual model results for each group-factor combination (where X is the number of groups and Y is the number of factors), containing SDresult and variable order

IC.Rdata

Contains IC.matrix, timing information, and data. Load this file to compare models and identify the optimal configuration.

Note

This function can be computationally intensive as it fits multiple models. Consider running on high-performance computing resources for large datasets or extensive model grids. The function includes error handling to continue execution even if individual models fail to converge.

See Also

BCFM.fit for fitting a single model, init.data, initialize.model.attributes, initialize.hyp.parm

Examples

# Run model selection using the included simulated dataset
data(sim.data)
BCFM.model.selection(
  data = sim.data,
  cluster.vars = paste0("V", 1:5),
  grouplist = 2:3,
  factorlist = 1:2,
  n.iter = 100,
  every = 10,
  burnin = 10
)

# Load and examine IC results
load(file.path(tempdir(), "IC.Rdata"))
print(IC.matrix)

Gibbs sampler of BCFM

Description

It runs a Gibbs sampler for common factors X, factor loadings B, group mean mu, group covariate Omega, idiosyncratic variances sigma^2, group assignment Z and group probabilities probs. This function uses Rcpp.

Usage

BCFMcpp(data, model.attributes, hyp.parm, n.iter, every = 1, verbose = FALSE)

Arguments

Value

List of Gibbs sampler of parameters described in description

Information Criterion. Very close to the original BIC method, but this uses the integrated likelihood instead.

Description

It finds a Laplace-Metropolis marginal density of likelihood using posterior mean. It also uses Woodbury lemma for fast calculation

Usage

IC(data.row, Gibbs, model.attributes, cluster.size = 0.05, burnin = NA)

Arguments

Value

The value of Laplace-Metropolis marginal density

Get the mode of a vector

Description

The function returns the mode of a vector.

Usage

getmode(v)

Arguments

Value

The mode of v.

Build factor loadings plot

Description

The function builds a column-wise plots of factor loadings. The parameters fixed at 1 are displayed with red dashed vertical lines.

Usage

ggplot_B.CI(
  Gibbs,
  true.val = NA,
  burnin = NA,
  permutation = NA,
  main.bool = TRUE,
  var_labels = NULL
)

Arguments

Value

A ggplot object

Trace plot for posterior of factor loadings

Description

It returns a trace plot of factor loadings, B, showing MCMC convergence

Usage

ggplot_B.trace(
  Gibbs,
  burnin = NA,
  permutation = NA,
  true.val = NA,
  factor.num = 1,
  var_labels = NULL
)

Arguments

Value

A ggplot object showing trace plots

Plot IC Matrix from Model Selection

Description

Creates two plots showing IC values across different numbers of groups and factors to help identify the optimal model configuration.

Usage

ggplot_IC(matrix, factor_list = 2:4, group_list = 2:4, combine = TRUE)

Arguments

Details

The function creates two complementary visualizations:

• Plot 1: IC vs. number of groups, with lines for each number of factors

• Plot 2: IC vs. number of factors, with lines for each number of groups

Lower IC values indicate better model fit.

Value

If combine = TRUE, a combined ggplot object. If combine = FALSE, a list with two elements:

by_groups

ggplot showing IC vs. number of groups

by_factors

ggplot showing IC vs. number of factors

Examples

# Create a toy IC matrix for demonstration
IC.matrix <- matrix(c(100, 95, 90, 92, 88, 85), nrow = 2, ncol = 3)
rownames(IC.matrix) <- paste0("G", 2:3)
colnames(IC.matrix) <- paste0("F", 1:3)

# Combined plot (default)
ggplot_IC(IC.matrix, factor_list = 1:3, group_list = 2:3)

# Separate ggplot objects
plots <- ggplot_IC(IC.matrix, factor_list = 1:3, group_list = 2:3,
                    combine = FALSE)
plots$by_groups
plots$by_factors

A heatmap of group assignments, Z using ggplot2

Description

A heatmap of group assignments, Z using ggplot2. It first sorts by the largest group with the most assigned.

Usage

ggplot_Zit.heatmap(
  Gibbs,
  true.val = NA,
  burnin = NA,
  main.title = "Cluster Assignment Heatmap",
  x.label = "Cluster",
  y.label = "Subject-Time"
)

Arguments

Value

A ggplot object

Plot Latent Factor Profiles by Cluster

Description

Visualizes the latent factor profile means (mu) for each cluster, similar to Latent Profile Analysis (LPA) plots

Usage

ggplot_latent.profiles(
  Gibbs,
  burnin = NA,
  factor_labels = NULL,
  cluster_names = NULL,
  colors = NULL,
  title = "Latent Factor Profiles by Cluster",
  x_label = "Factor",
  y_label = "Posterior Mean"
)

Arguments

Value

A ggplot object

Examples

# Fit a model first using the included simulated dataset
data(sim.data)
data.pre <- init.data(sim.data, paste0("V", 1:5))
model.attributes <- initialize.model.attributes(S = nrow(sim.data), times = 1,
                                                 R = 5, L = 2, G = 2)
cluster.hyperparms <- initialize.cluster.hyperparms(data.pre, model.attributes)
hyp.parm <- initialize.hyp.parm(model.attributes, cluster.hyperparms)
result <- BCFM.fit(data.pre, model.attributes, hyp.parm,
                   n.iter = 100, every = 10)

# Plot latent profiles
ggplot_latent.profiles(Gibbs = result$Result)

Density of group means mu using ggplot2

Description

The function returns multiple density plots of group mean parameter mu.

Usage

ggplot_mu.density(
  Gibbs,
  true.val = NA,
  add.legend = FALSE,
  burnin = NA,
  layout.dim = NA,
  main.title = "Posterior Densities of Group Means (mu)",
  x.label = "mu",
  y.label = "Density"
)

Arguments

Value

A ggplot object (grob from grid.arrange)

The density plot of the diagonal of group covariance, Omega, with ggplot2

Description

It returns multiple plots of the diagonal of group covariance, Omega using ggplot2 and gridExtra. It returns the result by each factor, different colors representing different factors

Usage

ggplot_omega.density(
  Gibbs,
  group.select = 1,
  true.val = NA,
  burnin = NA,
  main.title = "Posterior Densities of Omega",
  x.label = "Value",
  y.label = "Density",
  factor_labels = NULL,
  show.offdiag = TRUE
)

Arguments

Value

A ggplot object showing densities of Omega elements

Density plot for posterior of probabilities

Description

The function returns a density plot of the cluster assignment probability, p. Different color represent different Clusters.

Usage

ggplot_probs.density(
  Gibbs,
  burnin = NA,
  truep = NA,
  main.title = "Posterior Densities of Cluster Probabilities",
  x.label = "Probability",
  y.label = "Density",
  cluster_names = NULL
)

Arguments

Value

A ggplot object (grob from grid.arrange) with plot and legend

Trace plot of probabilities parameter

Description

It returns a trace plot of probabilities parameter after burn-in. Different colors represent different groups.

Usage

ggplot_probs.trace(
  Gibbs,
  burnin = NA,
  main.title = "Trace Plot: Cluster Probabilities",
  x.label = "BCFM Iteration (post burn-in)",
  y.label = "Probability",
  cluster_names = NULL
)

Arguments

Value

A ggplot object showing trace plots

A credible interval plot of posterior of sigma squared

Description

It returns a credible interval plot of idiosyncratic variance, sigma squared. The lines are 95% intervals, while the circles are posterior mean.

Usage

ggplot_sigma2.CI(
  Gibbs,
  burnin = NA,
  permutation = NA,
  main.bool = TRUE,
  main.title = NULL,
  x.label = "Variables",
  y.label = "Variance",
  var_labels = NULL
)

Arguments

Value

A ggplot object showing the 95\ posterior means (points) of the idiosyncratic variances sigma squared for each variable.

A credible interval plot of posterior of factor loadings covariance, tau

Description

It returns a credible interval plot of factor loadings covariance, tau. The lines are 95% intervals, while the circles are posterior mean.

Usage

ggplot_tau.CI(
  Gibbs,
  burnin = NA,
  true.val = NA,
  main.bool = TRUE,
  main.title = NULL,
  x.label = "Factor",
  y.label = "Tau",
  factor_labels = NULL
)

Arguments

Value

A ggplot object

Variability explained by factors

Description

Plots the proportion of variability explained by each factor based on eigenvalues of the correlation matrix. Useful for determining the number of factors.

Usage

ggplot_variability(
  data,
  nfactors = 5,
  main.title = "Proportion of Variability Explained by Factors",
  x.label = "Number of Factors",
  y.label = "Proportion of Variability"
)

Arguments

Value

A ggplot object

Initialize Data Array for BCFM Model

Description

Prepares the input data by converting it into a 3D array format required by the BCFM model. If data is already in the correct 3D array format, it returns the data as-is. Takes selected clustering variables and creates an array with dimensions (observations, variables, time points).

Usage

init.data(data, cluster.vars = NULL)

Arguments

Value

A 3D array with dimensions (n, p, t) where n is the number of observations, p is the number of clustering variables, and t is the number of time points (defaults to 1 for cross-sectional data).

Examples

# Example 1: Data frame with variable selection
data <- data.frame(x1 = rnorm(100), x2 = rnorm(100), x3 = rnorm(100))
cluster.vars <- c("x1", "x2", "x3")
data.pre <- init.data(data, cluster.vars)

# Example 2: Matrix (uses all columns)
data_matrix <- matrix(rnorm(300), nrow = 100, ncol = 3)
data.pre <- init.data(data_matrix)

# Example 3: 3D array (returns as-is)
data_3d <- array(rnorm(1500), dim = c(100, 3, 5))
data.pre <- init.data(data_3d)  # Returns unchanged

Initialize cluster hyperparameters

Description

The function returns list of hyperparmeters for Omegas and mus.

Usage

initialize.cluster.hyperparms(
  data,
  model.attributes,
  covariance = FALSE,
  diag.Psi = FALSE,
  vague.mu = FALSE,
  zero.mu = FALSE,
  seed = NULL
)

Arguments

Value

A list of mean and variance hyperparameter of mu, and scale hyperparameter of Omega

Initialize hyperparmeters for BCFM model

Description

The function returns a list of hyperparameters of Omega, sigma^2 and mu. It also calls the results from cluster.hyperparms and information from model.attributes.

Usage

initialize.hyp.parm(
  model.attributes,
  cluster.hyperparms,
  n.sigma = 2.2,
  n.s2.sigma = 0.1,
  n.tau = 1,
  n.s2.tau = 1,
  omega.diag.nu = 2,
  p.exponent = NA
)

Arguments

Value

A list of fixed hyperparmeters of mu, Omega and sigma squared.

Build model attributes from the dataset

Description

Basic setting and information of the dataset: number of subjects, timepoints, variables, factors and groups.

Usage

initialize.model.attributes(S = 216, times = 5, R = 9, L = 3, G = 4)

Arguments

Value

A list of model attributes

Order of permutation by the largest absolute value in each eigenvector

Description

It finds the vector of permutation to permute data by its largest absolute value in each eigenvector. It sets the order by specified number of factors, and the rest is ordered as they were.

Usage

permutation.order(data, covariance = FALSE, L = 1, fa = TRUE)

Arguments

Value

The vector of permutation

Permute the dataset by the largest absolute value in each eigenvector, and scale

Description

It finds the vector of permutation to permute data by its largest absolute value in each eigenvector. It sets the order by specified number of factors, and the rest is ordered as they were. The data is permuted, and if needed, scaled.

Usage

permutation.scale(
  data,
  permutation = NA,
  covariance = FALSE,
  return.array = TRUE,
  num.layers = 1
)

Arguments

Value

The dataset that is permuted, either in matrix or array

Simulated data for BCFM model

Description

A simulated dataset for demonstrating the Bayesian Consensus Factor Model.

Usage

sim.data

Format

A data frame with 200 rows and 20 variables (V1 through V20): All variables are simulated numeric values.

Source

Simulated data generated for package examples

Examples

data(sim.data)
dim(sim.data)
head(sim.data)