Computing prior hyperparameters

---Special considerations: this portion is highly parallelizable---

We are now just about ready to set up our MCMC. First, we need to determine the hyperparameters in the priors of our Gaussian mixture. These are all calculated in an empirical Bayesian manner -- that is, we can recycle information from the pairwise fits to inform our priors in the full-information mixture. This task can be split into 2 sub-tasks:

1. computing the prior hyperparameters for the cluster mixing weights

2. computing every other hyperparameter

The former is the most essential, as it helps us remove more candidate latent classes, ensuring that the number of clusters is fewer than the number of observations. An important note: this is the only step of CLIMB that requires some sort of human intervention, but it does need to happen. A threshold, called $\delta$ in the manuscript, determines how strict one is about including classes in the final model. $\delta\in\{0,1,\ldots,\binom{D}{2}\}$. We will get into selecting $\delta$ shortly.

To get the prior weights on each candidate latent class, use the function get_prior_weights(). This function defaults to the settings used in the CLIMB manuscript. The user can specify:

1. reduced_classes: the matrix of candidate latent classes generated by get_reduced_classes()

2. fits: the list of pairwise fits generated by get_pairwise_fits()

3. parallel: logical specifying if the analysis should be run in parallel (defaults to FALSE)

4. ncores: if in parallel, how many cores to use. Defaults to 20.

5. delta: this is the range of thresholds to try, but it will defaults to a sequence of all possible thresholds.

NB: while parallelization is always available here, it is not always necessary. Speed of this portion depends on sample size, dimension, and the number of candidate latent classes (in reduced_classes).

Now, we are ready to compute the prior weights.

# Read in the candidate latent classes produced in the last step
reduced_classes <- readr::read_tsv("output/red_class.txt", col_names = FALSE)

# load in the pairwise fits from the first step
# (in this example case, I am simply loading the data from the package)
data("fits")

# Compute the prior weights
prior_weights <- get_prior_weights(reduced_classes, fits, parallel = FALSE)

prior_weights is a list of vectors. Each vector corresponds to the computed prior weights for a given value of $\delta$. Here, prior_weights[[j]] corresponds to the prior weights when $\delta = j-1$.

Here comes the human intervention

We can plot how the number of latent classes included in the final model changes as we relax $\delta$.

# this is just grabbing the sample size and dimension
n <- length(fits[[1]]$cluster)
D <- as.numeric(strsplit(tail(names(fits),1), "_")[[1]][2])

# to avoid degenerate distributions, we will only keep clusters such that the prior
# weight times the sample size is greater than the dimension.
plot(
    0:choose(D,2),
    sapply(prior_weights, function(X)
    sum(X * n > D)),
    ylab = "number of retained classes",
    xlab =  expression(delta))

Past versions of unnamed-chunk-2-1.png

This toy example is much cleaner than a real data set, but typically we expect to see that, as we relax $\delta$ away from 0, more classes are included in the final model. We have not identified a uniformly best way to select $\delta$; a decent rule of thumb has simply been to include as many classes as one can while retaining computational feasibility, and selecting the smallest value of $\delta$ that gives this result. In this toy example, we might as well retain all classes, and thus select the prior weights corresponding to $\delta = 1$. Let's store that in the variable p.

# Select out prior weights for delta = 1
p <- prior_weights[[2]]

# Filter out classes which have too small of a prior weight
# (In this toy example, we actually retain everything,
#   but this is not typical for higher-dimensional/empirical analyses)
retained_classes <- reduced_classes[p * n > D, ]

Warning: The `i` argument of ``[`()` can't be a matrix as of tibble 3.0.0.
Convert to a vector.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

p <- p[p * n > D,]

# save the retained classes for downstream analysis
readr::write_tsv(retained_classes, path = "output/retained_classes.txt", col_names = FALSE)

Obtaining the remaining hyperparameters

Now that the human intervention is over, the rest is simple. Just use the function get_hyperparameters() to compute empirical estimates of the remaining prior hyperparameters.

# load the data back in
data("sim")

# obtain the hyperparameters
hyp <- get_hyperparameters(sim$data, fits, retained_classes, p)

# view the output
str(hyp)

List of 4
 $ Psi0  : num [1:3, 1:3, 1:13] 0.816 0 0.488 0 1 ...
 $ mu0   : num [1:13, 1:3] 2.68 2.68 2.68 0 0 ...
 $ alpha : num [1:13] 0.101 0.1403 0.0207 0.0988 0.0502 ...
 $ kappa0: num [1:13] 151 210 31 148 75 178 166 110 46 136 ...

hyp$kappa0 controls the informativity of the priors. To reduce informativity, one can make the elements of kappa0 smaller (but still larger than $D$!). For example, you could use something like hyp$kappa0 <- rep(10, D), instead of the automatic choice (proportional to hyp$alpha) which is returned from the get_hyperparameters function.

We can save these hyperparameters for the next step in the analysis:

save(hyp, file = "output/hyperparameters.Rdata")

After all that work, we have a model to describe our data, and are ready to run the MCMC.

A note on relaxing $\delta$: In particularly challenging/high-dimensional problems, relaxing $\delta$ can start to introduce chunks of classes that look very similar to one another, causing bloat in the candidate latent classes. If reduced_classes contains too many classes--say, in the hundreds of thousands--one can consider the following heuristic before running get_prior_weights(). The function reduce_by_hamming() filters out some classes that are too similar to one another (that is, within a given hamming distance to some other class). This computational advantage comes at the risk of loss of information due to stochasticity/heuristics. However, if certain canonical classes are in reduced_classes (the all 1 class, all 0 class, and all -1 class), you can force them to remain after applying reduce_by_hamming().

The function can be applied as follows:

reduced_by_hamming_classes <- reduce_by_hamming(reduced_classes, hamming_tol = 2)