Bayesian Mixture Modelling for Joint Model-Based Clustering/Classification and Batch Correction

Description

Semi-supervised and unsupervised Bayesian mixture models that simultaneously infer the cluster/class structure and a batch correction. Densities available are the multivariate normal and the multivariate t. The model sampler is implemented in C++. This package is aimed at analysis of low-dimensional data generated across several batches. See (Coleman et al. (2022))[https://doi.org/10.1101/2022.01.14.476352] for details of the model.

Author(s)

Stephen Coleman <stcolema@tcd.ie>, Paul D.W. Kirk, Chris Wallace

See Also

Useful links:

• https://github.com/stcolema/batchmix

• Report bugs at https://github.com/stcolema/batchmix/issues

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

# Which labels are observed
fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Sampling parameters
R <- 1000
thin <- 50

# Classification
samples <- runBatchMix(X,
  R,
  thin,
  batch_vec,
  "MVN",
  initial_labels = labels,
  fixed = fixed,
)

# Clustering
samples <- runBatchMix(X, R, thin, batch_vec, "MVT")

Minimum VI lower bound

Description

Local implementation of S. Wade's 'minVI' function from their 'mcclust.ext' package (available from github). Reimplemented here to avoid dependency on a non-CRAN package. Computes the lower bound to the posterior expected Variation of Information. For full details please see the aforementioned package and Wade and Ghahramani, 2018, 'Bayesian Cluster Analysis: Point Estimation and Credible Balls (with Discussion)'.

Usage

VI.lb(cls, psm)

Arguments

Value

A vector of the the lower bound of the Variation of Information for

Examples

## Not run: 
# MCMC samples and BIC vector
mcmc_outputs <- runMCMCChains(
  X,
  n_chains,
  R,
  thin,
  batch_vec,
  type
)

# Note that in this toy example we have not applied a burn in
psm <- createSimilarityMat(mcmc_outputs[[1]]$samples)
VI.lb(mcmc_outputs[[1]]$samples[1, ], psm)

## End(Not run)

Batch semisupervised mixture model

Description

A Bayesian mixture model with batch effects.

Usage

batchSemiSupervisedMixtureModel(
  X,
  R,
  thin,
  initial_labels,
  fixed,
  batch_vec,
  type,
  K_max = length(unique(initial_labels)),
  alpha = NULL,
  concentration = NULL,
  mu_proposal_window = 0.5^2,
  cov_proposal_window = 0.002,
  m_proposal_window = 0.3^2,
  S_proposal_window = 0.01,
  t_df_proposal_window = 0.015,
  m_scale = NULL,
  rho = 3,
  theta = 1,
  initial_class_means = NULL,
  initial_class_covariance = NULL,
  initial_batch_shift = NULL,
  initial_batch_scale = NULL,
  initial_class_df = NULL,
  verbose = TRUE
)

Arguments

Value

A named list containing the sampled partitions, cluster and batch parameters, model fit measures and some details on the model call.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Density choice
type <- "MVN"

# Sampling parameters
R <- 1000
thin <- 50

# MCMC samples and BIC vector
samples <- batchSemiSupervisedMixtureModel(
  X,
  R,
  thin,
  labels,
  fixed,
  batch_vec,
  type
)

# Given an initial value for the parameters
initial_class_means <- matrix(c(1, 1, 3, 4), nrow = 2)
initial_class_covariance <- array(c(1, 0, 0, 1, 1, 0, 0, 1),
  dim = c(2, 2, 2)
)

# We can use values from a previous chain
initial_batch_shift <- samples$batch_shift[, , R / thin]
initial_batch_scale <- matrix(
  c(1.2, 1.3, 1.7, 1.1, 1.4, 1.3, 1.2, 1.2, 1.1, 2.0),
  nrow = 2
)

samples <- batchSemiSupervisedMixtureModel(X,
  R,
  thin,
  labels,
  fixed,
  batch_vec,
  type,
  initial_class_means = initial_class_means,
  initial_class_covariance = initial_class_covariance,
  initial_batch_shift = initial_batch_shift,
  initial_batch_scale = initial_batch_scale
)

Calculate allocation probabilities

Description

Calculate the empirical allocation probability for each class based on the sampled allocation probabilities.

Usage

calcAllocProb(mcmc_samples, burn = 0, method = "median")

Arguments

Value

An N x K matrix of class probabilities.

Examples

# Data in matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Sampling parameters
R <- 1000
thin <- 50

# MCMC samples and BIC vector
samples <- batchSemiSupervisedMixtureModel(X, R, thin, labels, fixed, batch_vec, "MVN")

# Burn in
burn <- 20
eff_burn <- burn / thin

# Probability across classes
probs <- calcAllocProb(samples, burn = burn)

Check data generation inputs

Description

Checks that the inputs for the “generateBatchData“ function are correct. For internal use only.

Usage

checkDataGenerationInputs(
  N,
  P,
  group_means,
  group_std_devs,
  batch_shift,
  batch_scale,
  group_weights,
  batch_weights,
  type,
  group_dfs,
  frac_known,
  permute_variables,
  scale_data
)

Arguments

Value

No return value, called for side effects.

Examples

N <- 500
P <- 2
K <- 2
B <- 5
mean_dist <- 4
batch_dist <- 0.3
group_means <- seq(1, K) * mean_dist
batch_shift <- rnorm(B, mean = batch_dist, sd = batch_dist)
group_std_devs <- rep(2, K)
batch_scale <- rep(1.2, B)
group_weights <- rep(1 / K, K)
batch_weights <- rep(1 / B, B)
type <- "MVT"
group_dfs <- c(4, 7)
frac_known <- 0.3
permute_variables <- TRUE
scale_data <- FALSE

checkDataGenerationInputs(
  N,
  P,
  group_means,
  group_std_devs,
  batch_shift,
  batch_scale,
  group_weights,
  batch_weights,
  type,
  group_dfs,
  frac_known,
  permute_variables,
  scale_data
)

Check proposal windows

Description

Checks the proposal windows are acceptable.

Usage

checkProposalWindows(
  mu_proposal_window,
  cov_proposal_window,
  m_proposal_window,
  S_proposal_window,
  t_df_proposal_window,
  verbose = TRUE
)

Arguments

Value

No return value, called for side effects

Examples

checkProposalWindows(0.1, 0.2, 0.3, 0.1, 0.4, 0.3)

Collect acceptance rate

Description

Collects the acceptance rates for each parameter into a data.frame

Usage

collectAcceptanceRates(samples)

Arguments

Value

A wide data.frame of all the sampled parameters and the iteration.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Sampling parameters
R <- 1000
thin <- 50

# MCMC samples
samples <- runBatchMix(X, R, thin, batch_vec, "MVN",
  initial_labels = labels,
  fixed = fixed
)

# Acceptance rates
collectAcceptanceRates(samples)

Continue chain

Description

Continues sampling from a previous position for a given chain.

Usage

continueChain(mcmc_output, X, fixed, batch_vec, R, keep_old_samples = TRUE)

Arguments

Value

A named list containing the sampled partitions, cluster and batch parameters, model fit measures and some details on the model call.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Density choice
type <- "MVT"

# Sampling parameters
R <- 1000
thin <- 50

# MCMC samples and BIC vector
mcmc_output <- runBatchMix(
  X,
  R,
  thin,
  batch_vec,
  type,
  initial_labels = labels,
  fixed = fixed
)

# Given an initial value for the parameters
mcmc_output <- continueChain(
  mcmc_output,
  X,
  fixed,
  batch_vec,
  R,
)

Continue chains

Description

Continues sampling from a list of previous chains.

Usage

continueChains(mcmc_output, X, fixed, batch_vec, R, keep_old_samples = TRUE)

Arguments

Value

A named list containing the sampled partitions, cluster and batch parameters, model fit measures and some details on the model call.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Density choice
type <- "MVT"

# Sampling parameters
R <- 1000
thin <- 50
n_chains <- 4

# MCMC samples
mcmc_output <- runMCMCChains(
  X,
  n_chains,
  R,
  thin,
  batch_vec,
  type,
  initial_labels = labels,
  fixed = fixed
)

# Given an initial value for the parameters
new_output <- continueChains(
  mcmc_output,
  X,
  fixed,
  batch_vec,
  R,
  keep_old_samples = TRUE
)

Create Similarity Matrix

Description

Constructs a similarity matrix of the pairwise coclustering rate.

Usage

createSimilarityMat(allocations)

Arguments

Value

A symmetric n x n matrix (for n rows in cluster record) describing the fraction of iterations for which each pairwise combination of points are assigned the same label.

Gamma log-likelihood

Description

Used in calculating model probability in Metropolis-Hastings algorithm when proposals are from the Gamma distribution.

Usage

gammaLogLikelihood(x, shape, rate)

Arguments

Value

the unnormalised log-likelihood of x in a Gamma with parameters shape and rate.

Generate batch data

Description

Generate data from K multivaraite normal or multivariate t distributions with additional noise from batches. Assumes independence across columns. In each column the parameters are randomly permuted for both the groups and batches.

Usage

generateBatchData(
  N,
  P,
  group_means,
  group_std_devs,
  batch_shift,
  batch_scale,
  group_weights,
  batch_weights,
  type = "MVN",
  group_dfs = NULL,
  frac_known = 0.2,
  permute_variables = TRUE,
  scale_data = FALSE
)

Arguments

Value

A list of 5 objects; the data generated from the groups with and without batch effects, the label indicating the generating group, the batch label and the vector indicating training versus test.

Examples

N <- 500
P <- 2
K <- 2
B <- 5
mean_dist <- 4
batch_dist <- 0.3
group_means <- seq(1, K) * mean_dist
batch_shift <- rnorm(B, mean = batch_dist, sd = batch_dist)
std_dev <- rep(2, K)
batch_var <- rep(1.2, B)
group_weights <- rep(1 / K, K)
batch_weights <- rep(1 / B, B)
dfs <- c(4, 7)
my_data <- generateBatchData(
  N,
  P,
  group_means,
  std_dev,
  batch_shift,
  batch_var,
  group_weights,
  batch_weights,
  type = "MVT",
  group_dfs = dfs
)

Generate batch data

Description

Generate data from K multivaraite normal or multivariate t distributions with additional noise from batches. Assumes independence across columns. In each column the parameters are randomly permuted for both the groups and batches.

Usage

generateBatchDataLogPoisson(
  N,
  P,
  group_rates,
  batch_rates,
  group_weights,
  batch_weights,
  frac_known = 0.2,
  permute_variables = TRUE,
  scale_data = FALSE
)

Arguments

Value

A list of 5 objects; the data generated from the groups with and without batch effects, the label indicating the generating group, the batch label and the vector indicating training versus test.

Generate batch data from a multivariate t distribution

Description

Generate data from K multivariate t distributions with additional noise from batches. Assumes independence across columns. In each column the parameters are randomly permuted for both the groups and batches.

Usage

generateBatchDataMVT(
  N,
  P,
  group_means,
  group_std_devs,
  batch_shift,
  batch_scale,
  group_weights,
  batch_weights,
  dfs,
  frac_known = 0.2
)

Arguments

Value

A list of 5 objects; the data generated from the groups with and without batch effects, the label indicating the generating group, the batch label and the vector indicating training versus test.

Examples

N <- 500
P <- 2
K <- 2
B <- 5
mean_dist <- 4
batch_dist <- 0.3
group_means <- seq(1, K) * mean_dist
batch_shift <- rnorm(B, mean = batch_dist, sd = batch_dist)
std_dev <- rep(2, K)
batch_var <- rep(1.2, B)
group_weights <- rep(1 / K, K)
batch_weights <- rep(1 / B, B)
dfs <- c(4, 7)
my_data <- generateBatchDataMVT(
  N,
  P,
  group_means,
  std_dev,
  batch_shift,
  batch_var,
  group_weights,
  batch_weights,
  dfs
)

Generate batch data

Description

Generate data from groups across batches. Assumes independence across columns. In each column the parameters are randomly permuted for both the groups and batches.

Usage

generateBatchDataVaryingRepresentation(
  N,
  P,
  group_means,
  group_std_dev,
  batch_shift,
  batch_scale,
  group_weights,
  batch_weights,
  frac_known = 0.2
)

Arguments

Value

A list of 4 objects; the data generated from the groups with and without batch effects, the label indicating the generating group and the batch label.

Examples

N <- 500
P <- 2
K <- 2
B <- 5
mean_dist <- 4
batch_dist <- 0.3
group_means <- seq(1, K) * mean_dist
batch_shift <- rnorm(B, mean = batch_dist, sd = batch_dist)
std_dev <- rep(2, K)
batch_var <- rep(1.2, B)
group_weights <- matrix(
  c(
    0.8, 0.6, 0.4, 0.2, 0.2,
    0.2, 0.4, 0.6, 0.8, 0.8
  ),
  nrow = K, ncol = B, byrow = TRUE
)
batch_weights <- rep(1 / B, B)

my_data <- generateBatchDataVaryingRepresentation(
  N,
  P,
  group_means,
  std_dev,
  batch_shift,
  batch_var,
  group_weights,
  batch_weights
)

Generate group IDs

Description

Generate group IDs within “generateBatchData“.

Usage

generateGroupIDsInSimulator(
  N,
  K,
  B,
  batch_IDs,
  group_weights,
  varying_group_within_batch
)

Arguments

Value

A N-vector of group membership.

Examples

N <- 500
K <- 2
B <- 5
group_weights <- rep(1 / K, K)
batch_weights <- rep(1 / B, B)
batch_IDs <- sample(seq(1, B), N, replace = TRUE, prob = batch_weights)
varying_group_within_batch <- FALSE
group_IDs <- generateGroupIDsInSimulator(
  N,
  K,
  B,
  batch_IDs,
  group_weights,
  varying_group_within_batch
)

Generate initial labels

Description

For simulated data, generates an initial labelling for sampling.

Usage

generateInitialLabels(alpha, K, fixed, labels = NULL)

Arguments

Value

An N vector of labels.

Examples

N <- 500
P <- 2
K <- 2
B <- 5
mean_dist <- 4
batch_dist <- 0.3
cluster_means <- seq(1, K) * mean_dist
batch_shift <- rnorm(B, mean = batch_dist, sd = batch_dist)
std_dev <- rep(2, K)
batch_var <- rep(1.2, B)
cluster_weights <- rep(1 / K, K)
batch_weights <- rep(1 / B, B)

my_data <- generateBatchData(
  N,
  P,
  cluster_means,
  std_dev,
  batch_shift,
  batch_var,
  cluster_weights,
  batch_weights
)

initial_labels <- generateInitialLabels(1, K, my_data$fixed)

Get likelihood

Description

Extracts the model fit score from the mixture model output.

Usage

getLikelihood(mcmc_output, choice = "complete_likelihood")

Arguments

Value

A data.frame containing the model fit score of choice and the iteration.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Sampling parameters
R <- 100
thin <- 5

# MCMC samples and BIC vector
samples <- runBatchMix(X, R, thin, batch_vec, "MVN")

lkl_df <- getLikelihood(samples)

Get sampled batch shift

Description

Given an array of sampled batch scales from the “mixtureModel“ function, acquire a tidy version ready for “ggplot2“ use.

Usage

getSampledBatchScale(
  sampled_batch_scale,
  B = dim(sampled_batch_scale)[2],
  P = dim(sampled_batch_scale)[1],
  R = dim(sampled_batch_scale)[3],
  thin = 1
)

Arguments

Value

A data.frame of three columns; the parameter, the sampled value and the iteration.

Examples

# Data in matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Observed batches represented by integers
batch_vec <- sample(seq(1, 5), size = 100, replace = TRUE)

# MCMC iterations (this is too low for real use)
R <- 100
thin <- 5

# MCMC samples
samples <- runBatchMix(X, R, thin, batch_vec, "MVN")

batch_scale_df <- getSampledBatchShift(samples$batch_scale, R = R, thin = thin)

Get sampled batch shift

Description

Given an array of sampled batch mean shifts from the “mixtureModel“ function, acquire a tidy version ready for “ggplot2“ use.

Usage

getSampledBatchShift(
  sampled_batch_shift,
  B = dim(sampled_batch_shift)[2],
  P = dim(sampled_batch_shift)[1],
  R = dim(sampled_batch_shift)[3],
  thin = 1
)

Arguments

Value

A data.frame of three columns; the parameter, the sampled value and the iteration.

Examples

# Data in matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Observed batches represented by integers
batch_vec <- sample(seq(1, 5), size = 100, replace = TRUE)

# MCMC iterations (this is too low for real use)
R <- 100
thin <- 5

# MCMC samples
samples <- runBatchMix(X, R, thin, batch_vec, "MVN")

batch_shift_df <- getSampledBatchShift(samples$batch_shift,
  R = R,
  thin = thin
)

Get sampled cluster means

Description

Given an array of sampled cluster means from the “mixtureModel“ function, acquire a tidy version ready for “ggplot2“ use.

Usage

getSampledClusterMeans(
  sampled_cluster_means,
  K = dim(sampled_cluster_means)[2],
  P = dim(sampled_cluster_means)[1],
  R = dim(sampled_cluster_means)[3],
  thin = 1
)

Arguments

Value

A data.frame of three columns; the parameter, the sampled value and the iteration.

Examples

# Data in matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Observed batches represented by integers
batch_vec <- sample(seq(1, 5), size = 100, replace = TRUE)

# MCMC iterations (this is too low for real use)
R <- 100
thin <- 5

# MCMC samples
samples <- runBatchMix(X, R, thin, batch_vec, "MVN")

batch_shift_df <- getSampledClusterMeans(samples$means, R = R, thin = thin)

Inverse gamma log-likelihood

Description

Used in calculating model probability in Metropolis-Hastings algorithm when proposals are from the inverse-Gamma distribution.

Usage

invGammaLogLikelihood(x, shape, scale)

Arguments

Value

the unnormalised log-likelihood of x in a inverse-Gamma with parameters shape and scale.

Inverse-Wishart log-likelihood

Description

Used in calculating model probability in Metropolis-Hastings algorithm when proposals are from the Wishart distribution.

Usage

invWishartLogLikelihood(X, Psi, nu, P)

Arguments

Value

the unnormalised log-likelihood of X in a inverse-Wishart with parameters Psi and nu.

Minimium VI

Description

Local implementation of S. Wade's 'minVI' function from their 'mcclust.ext' package (available from github). Reimplemented here to avoid dependency on a non-CRAN package and we have dropped the 'greedy' method. Finds the optimal partition by minimising the lower bound to the Variation of Information obtained from Jensen's inequality where the expectation and log are reversed. For full details please see the aforementioned package and Wade and Ghahramani, 2018, 'Bayesian Cluster Analysis: Point Estimation and Credible Balls (with Discussion)'.

Usage

minVI(psm, cls.draw = NULL, method = "avg", max.k = NULL)

Arguments

Value

If ‘method' is '’all'' returns a matrix of four clusterings, one for each method and a repeat of that which performs best based on minimising the Variation of Information between the clustering and the PSM. Otherwise returns a vector. This is annotated with the attribute '"info"', a named list describing:

* '.$loss': the loss score used (Variation of Information)

* ‘.$maxNClusters': the 'max.k' value used by the '’comp'‘ and '’avg'' methods

* '.$expectedLoss': the estimated minimum Variation of Information for the point clustering(s)

* '.$method': the point method used to infer the clustering(s)

Names are due to legacy reasons - this function is replacing the 'salso::salso' function and name choices are to minimise workflow damage.

Examples

## Not run: 
# MCMC samples and BIC vector
mcmc_outputs <- runMCMCChains(
  X,
  n_chains,
  R,
  thin,
  batch_vec,
  type
)

# Note that in this toy example we have not applied a burn in
psm <- createSimilarityMat(mcmc_outputs[[1]]$samples)
cl_est <- minVI(psm, mcmc_outputs[[1]]$samples)

## End(Not run)

Plot acceptance rates

Description

Plot the acceptance rates for the parameters sampled in a Metropolis-Hastings step. Aiming for acceptance rates in [0.2, 0.5] for the class means, the batch effect on location and scale is a good rule of thumb. The class covariance should be in [0.35, 0.8] based on the authors' experience. The class degree of freedom appears to be prone to high acceptance rates, but aim to keep this above 0.2 at a minimum.

Usage

plotAcceptanceRates(mcmc_lst)

Arguments

Value

A ggplot object of the boxplots of acceptance rates for each parameter across chains.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Sampling parameters
R <- 500
thin <- 10
n_chains <- 4

# MCMC samples and BIC vector
mcmc_lst <- runMCMCChains(X, n_chains, R, thin, batch_vec, "MVN",
  initial_labels = labels,
  fixed = fixed
)

# Plot the acceptance rate of each parameter in the 4 chains
plotAcceptanceRates(mcmc_lst)

Plot likelihoods

Description

Plots the model fit for multiple chains.

Usage

plotLikelihoods(
  mcmc_outputs,
  choice = "complete_likelihood",
  colour_by_chain = TRUE
)

Arguments

Value

A ggplot2 object. Line plot of likelihood across iteration.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Sampling parameters
R <- 1000
thin <- 50
n_chains <- 4

# MCMC samples
samples <- runMCMCChains(X, n_chains, R, thin, batch_vec, "MVN",
  initial_labels = labels,
  fixed = fixed
)

p <- plotLikelihoods(samples)

Plot sampled batch means

Description

Plot the sampled values for the batch mean shifts in each dimension from the output of the mixture model functions. Not recommended for large B or P.

Usage

plotSampledBatchMeans(samples, burn_in = 0)

Arguments

Value

A ggplot object of the values in each sampled batch mean per iteration.

Examples

# Data in matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Observed batches represented by integers
batch_vec <- sample(seq(1, 5), size = 100, replace = TRUE)

# MCMC iterations (this is too low for real use)
R <- 100
thin <- 5

# MCMC samples and BIC vector
samples <- runBatchMix(X, R, thin, batch_vec, "MVN")

# Plot the sampled value of the batch mean shift against MCMC iteration
plotSampledBatchMeans(samples)

Plot sampled batch scales

Description

Plot the sampled values for the batch scale in each dimension from the output of the mixture model functions. Not recommended for large B or P.

Usage

plotSampledBatchScales(samples, burn_in = 0)

Arguments

Value

A ggplot object of the values in each sampled batch mean per iteration.

Examples

# Data in matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Observed batches represented by integers
batch_vec <- sample(seq(1, 5), size = 100, replace = TRUE)

# MCMC iterations (this is too low for real use)
R <- 100
thin <- 5

# MCMC samples and BIC vector
samples <- runBatchMix(X, R, thin, batch_vec, "MVN")

# Plot the sampled value of the batch scales against MCMC iteration
plotSampledBatchScales(samples)

Plot sampled cluster means

Description

Plot the sampled values for the cluster means in each dimension from the output of the mixture model functions. Not recommended for large K or P.

Usage

plotSampledClusterMeans(samples, burn_in = 0)

Arguments

Value

A ggplot object of the values in each sampled cluster mean per iteration.

Examples

# Data in matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Observed batches represented by integers
batch_vec <- sample(seq(1, 5), size = 100, replace = TRUE)

# MCMC iterations (this is too low for real use)
R <- 100
thin <- 5

# MCMC samples and BIC vector
samples <- runBatchMix(X, R, thin, batch_vec, "MVN")

# Plot the sampled value of the cluster means against MCMC iteration
plotSampledClusterMeans(samples)

Plot sampled vector parameter

Description

Plot the sampled values for a sampled vector from the output of the “mixtureModel“ function. Not recommended for large B or P.

Usage

plotSampledParameter(samples, parameter, R = NULL, thin = 1, burn_in = 0)

Arguments

Value

A ggplot object of the values in each sampled batch mean per iteration.

Examples

# Data in matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Observed batches represented by integers
batch_vec <- sample(seq(1, 3), size = 100, replace = TRUE)

# MCMC iterations (this is too low for real use)
R <- 50
thin <- 1

# MCMC samples and BIC vector
samples <- runBatchMix(X, R, thin, batch_vec, "MVN", K_max = 8)

# Plot the sampled value of the cluster means against MCMC iteration
parameter <- "means"
plotSampledParameter(samples, parameter, R, thin)

Predict class

Description

Predicts a final class for each item given a matrix of allocation probabilities.

Usage

predictClass(prob)

Arguments

Value

An N vector of class allocations.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Sampling parameters
R <- 1000
thin <- 50

# MCMC samples and BIC vector
samples <- batchSemiSupervisedMixtureModel(
  X,
  R,
  thin,
  labels,
  fixed,
  batch_vec,
  "MVN"
)

# Burn in
burn <- 200
eff_burn <- burn / thin

# Probability across classes
probs <- calcAllocProb(samples, burn = burn)

# Predict the class
preds <- predictClass(probs)

Predict from multiple MCMC chains

Description

Applies a burn in to and finds a point estimate by combining multiple chains of “callMDI“.

Usage

predictFromMultipleChains(
  mcmc_outputs,
  burn,
  point_estimate_method = "median",
  chains_already_processed = FALSE
)

Arguments

Value

A named list of quantities related to prediction/clustering:

* “allocation_probability“: List with an $(N x K)$ matrix if the model is semi-supervised. The point estimate of the allocation probabilities for each data point to each class.

* “prob“: $N$ vector of the point estimate of the probability of being allocated to the class with the highest probability.

* “pred“: $N$ vector of the predicted class for each sample. If the model is unsupervised then the “salso“ function from Dahl et al. (2021) is used on the sampled partitions using the default settings.

* “samples“: List of sampled allocations for each view. Columns correspond to items being clustered, rows to MCMC samples.

Examples

# Data dimensions
N <- 600
P <- 4
K <- 5
B <- 7

# Generating model parameters
mean_dist <- 2.25
batch_dist <- 0.3
group_means <- seq(1, K) * mean_dist
batch_shift <- rnorm(B, mean = batch_dist, sd = batch_dist)
std_dev <- rep(2, K)
batch_var <- rep(1.2, B)
group_weights <- rep(1 / K, K)
batch_weights <- rep(1 / B, B)
dfs <- c(4, 7, 15, 60, 120)

my_data <- generateBatchData(
  N,
  P,
  group_means,
  std_dev,
  batch_shift,
  batch_var,
  group_weights,
  batch_weights,
  type = "MVT",
  group_dfs = dfs
)


X <- my_data$observed_data

true_labels <- my_data$group_IDs
fixed <- my_data$fixed
batch_vec <- my_data$batch_IDs

alpha <- 1
initial_labels <- generateInitialLabels(alpha, K, fixed, true_labels)

# Sampling parameters
R <- 1000
thin <- 25
burn <- 100
n_chains <- 2

# Density choice
type <- "MVT"

# MCMC samples and BIC vector
mcmc_outputs <- runMCMCChains(
  X,
  n_chains,
  R,
  thin,
  batch_vec,
  type,
  initial_labels = initial_labels,
  fixed = fixed
)
ensemble_mod <- predictFromMultipleChains(mcmc_outputs, burn)

Prepare initial values

Description

Prepares user given values for input into the C++ function.

Usage

prepareInitialParameters(
  initial_class_means,
  initial_class_covariance,
  initial_batch_shift,
  initial_batch_scale,
  initial_class_df,
  P,
  K,
  B,
  type
)

Arguments

Value

A named list containing the different parameters.

Process MCMC chain

Description

Applies a burn in to and finds a point estimate for the output of “batchSemiSupervisedMixtureModel“.

Usage

processMCMCChain(mcmc_output, burn, point_estimate_method = "median")

Arguments

Value

A named list similar to the output of “batchSemiSupervisedMixtureModel“ with some additional entries:

* “mean_est“: $(P x K)$ matrix. The point estimate of the cluster means with columns corresponding to clusters.

* “cov_est“: $(P x P x K)$ array. The point estimate of the cluster covariance matrices with slices corresponding to clusters.

* “shift_est“: $(P x B)$ matrix. The point estimate of the batch shift effect with columns corresponding to batches.

* “scale_est“: $(P x B)$ matrix. The point estimate of the batch scale effects. The $bth$ column contains the diagonal entries of the scaling matrix for the $bth£ batch.

* “mean_sum_est“: $(P x K x B)$ array. The point estimate of the sum of the cluster means and the batch shift effect with columns corresponding to clusters and slices to batches.

* “cov_comb_est“: List of length $B$, with each entry being a $(P x P x K)$ array. The point estimate of the combination of the cluster covariance matrices and the batch scale effect with list entries corresponding to batches and slices of each array corresponding to clusters.

* “inferred_dataset“: $(N x P)$ matrix. The inferred “batch-free” dataset.

* “allocation_probability“: $(N x K)$ matrix. The point estimate of the allocation probabilities for each data point to each class.

* “prob“: $N$ vector. The point estimate of the probability of being allocated to the class with the highest probability.

* “pred“: $N$ vector. The predicted class for each sample.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Sampling parameters
R <- 1000
burn <- 250
thin <- 50

# MCMC samples
samples <- runBatchMix(X, R, thin, batch_vec, "MVN",
  initial_labels = labels,
  fixed = fixed
)

# Process the MCMC samples
processed_samples <- processMCMCChain(samples, burn)

Process MCMC chains

Description

Applies a burn in to and finds a point estimate for each of the chains outputted from “runMCMCChains“.

Usage

processMCMCChains(mcmc_lst, burn, point_estimate_method = "median")

Arguments

Value

A named list similar to the output of “batchSemiSupervisedMixtureModel“ with some additional entries:

* “mean_est“: $(P x K)$ matrix. The point estimate of the cluster means with columns corresponding to clusters.

* “cov_est“: $(P x P x K)$ array. The point estimate of the cluster covariance matrices with slices corresponding to clusters.

* “shift_est“: $(P x B)$ matrix. The point estimate of the batch shift effect with columns corresponding to batches.

* “scale_est“: $(P x B)$ matrix. The point estimate of the batch scale effects. The $bth$ column contains the diagonal entries of the scaling matrix for the $bth£ batch.

* “mean_sum_est“: $(P x K x B)$ array. The point estimate of the sum of the cluster means and the batch shift effect with columns corresponding to clusters and slices to batches.

* “cov_comb_est“: List of length $B$, with each entry being a $(P x P x K)$ array. The point estimate of the combination of the cluster covariance matrices and the batch scale effect with list entries corresponding to batches and slices of each array corresponding to clusters.

* “inferred_dataset“: $(N x P)$ matrix. The inferred “batch-free” dataset.

* “allocation_probability“: $(N x K)$ matrix. The point estimate of the allocation probabilities for each data point to each class.

* “prob“: $N$ vector. The point estimate of the probability of being allocated to the class with the highest probability.

* “pred“: $N$ vector. The predicted class for each sample.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Sampling parameters
R <- 1000
burn <- 250
thin <- 50
n_chains <- 4

# MCMC samples
samples <- runMCMCChains(X, n_chains, R, thin, batch_vec, "MVN",
  initial_labels = labels,
  fixed = fixed
)

# Process the MCMC samples
processed_samples <- processMCMCChains(samples, burn)

Random Draw From Stick Breaking Prior

Description

Draw weights from the stick-breaking prior.

Usage

rStickBreakingPrior(alpha, K)

Arguments

Value

A vector of component weights.

Examples

weights <- rStickBreakingPrior(1, 50)

Run Batch Mixture Model

Description

Runs a MCMC chain for a Bayesian mixture model which models both batch effects and class/cluster structure.

Usage

runBatchMix(
  X,
  R,
  thin,
  batch_vec,
  type,
  K_max = NULL,
  initial_labels = NULL,
  fixed = NULL,
  alpha = 1,
  mu_proposal_window = 0.5^2,
  cov_proposal_window = 0.002,
  m_proposal_window = 0.3^2,
  S_proposal_window = 0.01,
  t_df_proposal_window = 0.015,
  m_scale = NULL,
  rho = 3,
  theta = 1,
  initial_class_means = NULL,
  initial_class_covariance = NULL,
  initial_batch_shift = NULL,
  initial_batch_scale = NULL,
  initial_class_df = NULL,
  verbose = TRUE
)

Arguments

Value

A named list containing the sampled partitions, cluster and batch parameters, model fit measures and some details on the model call.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Density choice
type <- "MVN"

# Sampling parameters
R <- 1000
thin <- 50

# MCMC samples
mcmc_out <- runBatchMix(
  X,
  R,
  thin,
  batch_vec,
  type,
  initial_labels = labels,
  fixed = fixed
)

# Given an initial value for the parameters
initial_class_means <- matrix(c(1, 1, 3, 4), nrow = 2)
initial_class_covariance <- array(c(1, 0, 0, 1, 1, 0, 0, 1),
  dim = c(2, 2, 2)
)

# We can use values from a previous chain
initial_batch_shift <- mcmc_out$batch_shift[, , R / thin]
initial_batch_scale <- matrix(
  c(1.2, 1.3, 1.7, 1.1, 1.4, 1.3, 1.2, 1.2, 1.1, 2.0),
  nrow = 2
)

mcmc_out <- runBatchMix(X,
  R,
  thin,
  batch_vec,
  type,
  initial_labels = labels,
  fixed = fixed,
  initial_class_means = initial_class_means,
  initial_class_covariance = initial_class_covariance,
  initial_batch_shift = initial_batch_shift,
  initial_batch_scale = initial_batch_scale
)

Run MCMC Chains

Description

Run multiple chains of the batch mixture model of the same type.

Usage

runMCMCChains(
  X,
  n_chains,
  R,
  thin,
  batch_vec,
  type,
  K_max = NULL,
  initial_labels = NULL,
  fixed = NULL,
  alpha = 1,
  mu_proposal_window = 0.5^2,
  cov_proposal_window = 0.002,
  m_proposal_window = 0.3^2,
  S_proposal_window = 0.01,
  t_df_proposal_window = 0.015,
  m_scale = 0.01,
  rho = 3,
  theta = 1,
  initial_class_means = NULL,
  initial_class_covariance = NULL,
  initial_batch_shift = NULL,
  initial_batch_scale = NULL,
  initial_class_df = NULL,
  verbose = TRUE
)

Arguments

Value

A list of named lists. Each entry is the output of “runBatchMix“.

Examples

# Data in a matrix format
X <- matrix(c(rnorm(100, 0, 1), rnorm(100, 3, 1)), ncol = 2, byrow = TRUE)

# Initial labelling
labels <- c(
  rep(1, 10),
  sample(c(1, 2), size = 40, replace = TRUE),
  rep(2, 10),
  sample(c(1, 2), size = 40, replace = TRUE)
)

fixed <- c(rep(1, 10), rep(0, 40), rep(1, 10), rep(0, 40))

# Batch
batch_vec <- sample(seq(1, 5), replace = TRUE, size = 100)

# Sampling parameters
R <- 1000
thin <- 50
n_chains <- 4

# MCMC samples
samples <- runMCMCChains(X, n_chains, R, thin, batch_vec, "MVN",
  initial_labels = labels,
  fixed = fixed
)

Sample mixture of multivariate normal distributions with batch effects

Description

Performs MCMC sampling for a mixture model with batch effects.

Usage

sampleMVN(
  X,
  K,
  B,
  labels,
  batch_vec,
  mu_proposal_window,
  cov_proposal_window,
  m_proposal_window,
  S_proposal_window,
  R,
  thin,
  concentration,
  m_scale,
  rho,
  theta,
  initial_mu,
  initial_cov,
  initial_m,
  initial_S,
  mu_initialised,
  cov_initialised,
  m_initialised,
  S_initialised,
  sample_m_scale
)

Arguments

Value

Named list of the different quantities drawn by the sampler.

Sample mixture of multivariate t-distributions with batch effects

Description

Performs MCMC sampling for a MVT mixture model with batch effects.

Usage

sampleMVT(
  X,
  K,
  B,
  labels,
  batch_vec,
  mu_proposal_window,
  cov_proposal_window,
  m_proposal_window,
  S_proposal_window,
  t_df_proposal_window,
  R,
  thin,
  concentration,
  m_scale,
  rho,
  theta,
  initial_mu,
  initial_cov,
  initial_df,
  initial_m,
  initial_S,
  mu_initialised,
  cov_initialised,
  df_initialised,
  m_initialised,
  S_initialised,
  sample_m_scale
)

Arguments

Value

Named list of the different quantities drawn by the sampler.

Sample prior labels

Description

Generate labels from the stick-breaking prior.

Usage

samplePriorLabels(alpha, K, N)

Arguments

Value

A vector of labels.

Examples

initial_labels <- samplePriorLabels(1, 50, 100)

Sample semi-supervised MVN Mixture model

Description

Performs MCMC sampling for a mixture model.

Usage

sampleSemisupervisedMVN(
  X,
  K,
  B,
  labels,
  batch_vec,
  fixed,
  mu_proposal_window,
  cov_proposal_window,
  m_proposal_window,
  S_proposal_window,
  R,
  thin,
  concentration,
  m_scale,
  rho,
  theta,
  initial_mu,
  initial_cov,
  initial_m,
  initial_S,
  mu_initialised,
  cov_initialised,
  m_initialised,
  S_initialised,
  sample_m_scale
)

Arguments

Value

Named list of the different quantities drawn by the sampler.

Sample semi-supervised MVT Mixture model

Description

Performs MCMC sampling for a mixture model.

Usage

sampleSemisupervisedMVT(
  X,
  K,
  B,
  labels,
  batch_vec,
  fixed,
  mu_proposal_window,
  cov_proposal_window,
  m_proposal_window,
  S_proposal_window,
  t_df_proposal_window,
  R,
  thin,
  concentration,
  m_scale,
  rho,
  theta,
  initial_mu,
  initial_cov,
  initial_df,
  initial_m,
  initial_S,
  mu_initialised,
  cov_initialised,
  df_initialised,
  m_initialised,
  S_initialised,
  sample_m_scale
)

Arguments

Value

Named list of the different quantities drawn by the sampler.

Wishart log-likelihood

Description

Used in calculating model probability in Metropolis-Hastings algorithm when proposals are from the Wishart distribution.

Usage

wishartLogLikelihood(X, V, n, P)

Arguments

Value

the unnormalised log-likelihood of X in a Wishart with parameters V and n.