Tutorials › Multivariate Statistics (Graduate) › Searching Across Cluster Counts and Covariance Families

Model-Based Clustering · Tutorial 666 of 975

Searching Across Cluster Counts and Covariance Families

Build and audit a candidate-model search that compares cluster counts and covariance structures using maximized likelihoods, parameter counts, and selection criteria.

Advanced 11 min read

What You'll Learn

  • Define a candidate model by its cluster count and covariance family
  • Count fitted parameters in Gaussian mixture candidates
  • Organize repeated starts and convergence checks reproducibly
  • Calculate and interpret AIC and BIC comparisons
  • Apply likelihood-improvement thresholds when comparing parameter counts
  • Record failed fits and report a transparent candidate search

Why a Search Needs a Fixed Protocol

Choosing a Gaussian clustering model means making two decisions together: how many components to fit and which covariance constraints to impose. A larger cluster count or a more flexible covariance family can improve the maximized likelihood, but it also adds parameters and can make fitting less stable. A useful search therefore needs a defined candidate set and a consistent fitting and scoring protocol—not just a list of whichever models happened to converge.

The search described here treats each combination of cluster count and covariance family as a candidate. It fits every candidate to the same observations, keeps track of repeated optimization starts, and records the maximized observed-data log-likelihood, parameter count, and selection criteria. Its aim is reproducibility: another analyst with the same data and protocol should be able to reconstruct the candidate table and understand why a candidate was selected.

Definition: For \(n\) observations \(x_1,\ldots,x_n\) in \(p\) dimensions, a candidate Gaussian mixture model is specified by a number of components \(G\) and a covariance family \(c\). It has density $$ f(x_i;\theta)=\sum_{g=1}^{G}\pi_g\,\phi_p(x_i;\mu_g,\Sigma_g), $$ where \(\pi_g\geq 0\), \(\sum_{g=1}^{G}\pi_g=1\), and the means and covariance matrices obey the constraints of family \(c\). Write \(\hat{\ell}_{G,c}\) for the largest observed-data log-likelihood attained by an accepted fit under that candidate, and \(d_{G,c}\) for its number of free scalar parameters.

Here \(\phi_p\) denotes a \(p\)-variate Gaussian density. The component labels are not observed; the likelihood being recorded is the log-likelihood of the mixture density for the observed data, not a likelihood conditional on known cluster memberships. Covariance families such as EII, VVV, and the intermediate parsimonious models impose different constraints on component covariance matrices. The search uses those constraints as candidate specifications rather than deciding the constraints after inspecting the scores.

Counting Parameters Before Comparing Fits

For a mixture with unrestricted component proportions, the mixing weights contribute \(G-1\) free parameters: the last weight is determined by the requirement that the weights sum to one. The component means contribute \(Gp\) parameters. If the covariance family contributes \(q_c(G,p)\) free covariance parameters, the total is

$$ d_{G,c}=(G-1)+Gp+q_c(G,p). $$

The covariance contribution depends on both the family and, for some families, the number of components. For example, EII has one shared spherical variance, so \(q_{\mathrm{EII}}=1\); VII has a separate spherical variance for each component, so \(q_{\mathrm{VII}}=G\). EEE has one unrestricted covariance matrix shared by all components, giving \(q_{\mathrm{EEE}}=p(p+1)/2\). VVV has an unrestricted covariance matrix for each component, giving \(q_{\mathrm{VVV}}=Gp(p+1)/2\).

Parameter-count derivation: A symmetric \(p\times p\) covariance matrix has \(p\) diagonal entries and \(p(p-1)/2\) distinct off-diagonal entries, for \(p(p+1)/2\) free entries in total. A single shared unrestricted matrix therefore contributes \(p(p+1)/2\); \(G\) unrestricted component-specific matrices contribute \(G p(p+1)/2\). Adding the \(G-1\) free weights and \(Gp\) means gives the stated total \(d_{G,c}\). Constrained families must be counted according to their actual constraints, not as if every covariance matrix were unrestricted. \(\square\)

In the standard volume-shape-orientation parameterization, the covariance code indicates whether each feature is equal (E) or variable (V) across components, or whether orientation is spherical (I). For instance, in \(p\) dimensions, the number of covariance parameters for EEV is \(1+(p-1)+G p(p-1)/2\): one common volume, one common shape with \(p-1\) free shape parameters, and \(G\) orientations, each with \(p(p-1)/2\) parameters. Parameter counts should be computed from the specific family definition used by the fitting implementation; they should not be inferred from a short family name alone.

Worked Example: Counting Parameters for Three Candidates

Suppose the data have \(p=3\) variables and a search includes \(G=4\). Count the parameters for EII, EEE, and VVV.

For each candidate the weights contribute \(G-1=3\) parameters and the means contribute \(Gp=4(3)=12\). Thus weights and means contribute \(15\) parameters before counting covariances.

For EII, the covariance is spherical and shared, so \(q_{\mathrm{EII}}=1\). The total is \(d_{4,\mathrm{EII}}=3+12+1=16\).

For EEE, there is one shared unrestricted \(3\times3\) covariance matrix. It has \(3(4)/2=6\) free entries, so \(d_{4,\mathrm{EEE}}=3+12+6=21\).

For VVV, there are four unrestricted covariance matrices, contributing \(4(6)=24\) parameters. Thus \(d_{4,\mathrm{VVV}}=3+12+24=39\). These counts are properties of the model specifications; they do not change with the estimated parameter values.

A Reproducible Fitting and Scoring Protocol

Before fitting, specify the cluster counts and covariance families to be searched, the data and preprocessing, the software and version, and the optimization settings. For each candidate, use multiple starting values because mixture likelihood optimization can reach different local maxima. Retain the largest log-likelihood among accepted converged fits, but preserve the results of all starts so the chosen value can be checked. A failed or nonconverged run is not a valid candidate score.

Also define in advance how inadmissible or numerically degenerate fits will be handled. For example, Gaussian mixture likelihoods can become unbounded when a component covariance collapses around an observation. An implementation may enforce eigenvalue constraints or other safeguards. The search should document such constraints and apply them consistently; otherwise likelihoods from candidates fitted under materially different restrictions are not directly comparable.

For a candidate with \(d\) parameters and maximized log-likelihood \(\hat{\ell}\), the commonly used information criteria are

$$ \mathrm{AIC}=-2\hat{\ell}+2d,\qquad \mathrm{BIC}=-2\hat{\ell}+d\log n. $$

Smaller values are preferred within a criterion. These formulae make it essential to record both the likelihood and parameter count: a score alone conceals whether a candidate was rewarded for fit or penalized for complexity. If classification-oriented selection is part of the prespecified analysis, the ICL-BIC version used in the earlier material can also be recorded. With posterior membership probabilities \(\tau_{ig}\), its common form is \(\mathrm{ICL}_{\mathrm{BIC}}=\mathrm{BIC}-2\sum_i\sum_g\tau_{ig}\log\tau_{ig}\). The exact convention should be named in the protocol.

Result: Likelihood-improvement threshold. Suppose candidate B has \(d_B=d_A+r\) parameters, where \(r>0\), and its maximized log-likelihood exceeds candidate A's by \(\Delta=\hat{\ell}_B-\hat{\ell}_A\). Then B has lower AIC exactly when \(2\Delta>2r\), and lower BIC exactly when \(2\Delta>r\log n\).

To verify the result, subtract A's AIC from B's AIC:

$$ \mathrm{AIC}_B-\mathrm{AIC}_A =-2(\hat{\ell}_B-\hat{\ell}_A)+2(d_B-d_A) =-2\Delta+2r. $$

This difference is negative precisely when \(2\Delta>2r\). Similarly,

$$ \mathrm{BIC}_B-\mathrm{BIC}_A =-2\Delta+r\log n, $$

which is negative precisely when \(2\Delta>r\log n\). These are algebraic comparisons of the stated scores; they do not assert that either criterion is an exact finite-sample model-selection guarantee for mixtures. \(\square\)

Worked Example: Comparing Two Cluster Counts

A data set has \(n=180\) observations in \(p=2\) dimensions. Two EEE candidates have log-likelihoods \(-510\) for \(G=2\) and \(-500\) for \(G=3\). Compute their parameter counts, AIC, and BIC.

For EEE in two dimensions, a shared unrestricted covariance has \(p(p+1)/2=3\) parameters. For \(G=2\), \(d=(2-1)+2(2)+3=8\). For \(G=3\), \(d=(3-1)+3(2)+3=11\).

For \(G=2\), \(\mathrm{AIC}=-2(-510)+2(8)=1020+16=1036\). Since \(\log(180)\approx5.193\), \(\mathrm{BIC}=1020+8(5.193)\approx1061.54\).

For \(G=3\), \(\mathrm{AIC}=-2(-500)+2(11)=1000+22=1022\). Its BIC is \(1000+11(5.193)\approx1057.12\). Both criteria are lower for \(G=3\).

The likelihood gain is \(\Delta=10\), and the increase in parameter count is \(r=3\). The AIC threshold is \(2r=6\), and the BIC threshold is \(r\log(180)\approx15.58\), each compared with \(2\Delta=20\). The gain exceeds both thresholds, confirming the computed rankings.

Putting the Candidate Search on Record

A candidate table should have one row per combination of \(G\) and covariance family. It should include the parameter count, best accepted log-likelihood, criteria, convergence status, and enough information to identify how that score was obtained. Where repeated starts matter, either report the best score and number of successful starts or provide a linked run-level record. A concise summary might look like this:

Worked Example: Auditing a Small Candidate Table

Consider \(n=120\), \(p=2\), and three candidates. The one-component EEE model has \(d=0+2+3=5\) and best log-likelihood \(-340\). The two-component EII model has \(d=1+4+1=6\) and log-likelihood \(-320\). The two-component VVV model has \(d=1+4+2(3)=11\) and log-likelihood \(-315\).

For one-component EEE, AIC is \(680+10=690\), and BIC is \(680+5\log(120)\approx703.94\). For two-component EII, AIC is \(640+12=652\), and BIC is \(640+6\log(120)\approx668.72\). For two-component VVV, AIC is \(630+22=652\), and BIC is \(630+11\log(120)\approx682.66\).

Thus AIC ties between the two \(G=2\) candidates, while BIC is smallest for \(G=2\), EII. The VVV candidate has the largest likelihood, but its extra covariance parameters offset that improvement under BIC. A responsible record would add the number of starts, the best and other converged log-likelihoods, the convergence tolerance, and any covariance constraints. The arithmetic alone does not establish that the fitted clusters are substantively meaningful.

1
Freeze the candidate grid.
List the planned values of \(G\), covariance families, data transformations, and any constraints before reviewing the winning criterion.
2
Fit each candidate repeatedly.
Use documented initialization and convergence settings. Retain the best accepted observed-data log-likelihood, and keep run-level information for auditing.
3
Compute and verify counts and scores.
Calculate \(d_{G,c}\), AIC, and BIC from the recorded values. Check that all candidates use the same observations and the same stated likelihood convention.
4
Report ranking and limitations.
Name the selected criterion, identify ties or failed candidates, and report the best candidate alongside enough of the search table to show the comparison.

Interpretation and Common Pitfalls

A search across covariance families helps distinguish a need for more components from a need for more flexible covariance structure. If a variable-volume family improves fit substantially at fixed \(G\), that is a different model comparison from increasing \(G\) while keeping the covariance family fixed. Reporting both dimensions of the grid makes that distinction visible. Selection criteria rank the fitted candidates under their penalties; they do not prove that a mixture's components are real, well separated, or scientifically useful.

  • Comparing likelihoods from different data. A missing-data or filtering difference changes the observations contributing to the likelihood. Fit and score candidates on the same analyzed data if the goal is a direct comparison.
  • Using a parameter count that ignores constraints. Counting every covariance as unrestricted over-penalizes constrained models; counting shared parameters repeatedly over-penalizes shared-covariance models. Derive the count from the fitted specification.
  • Keeping the first converged run. A converged run may not be the best available maximum. Use multiple starts and keep the highest accepted log-likelihood, while recording the run history.
  • Changing safeguards between candidates. Different covariance floors or convergence rules can alter fitted likelihoods. Specify the safeguards and apply them consistently.
  • Treating a criterion winner as a final scientific conclusion. AIC or BIC selects among the candidates and assumptions supplied. It does not replace checking whether the model is adequate for the intended use.
Key takeaway: A defensible candidate search is a documented comparison of a prespecified grid. For each candidate, report the same-data maximized likelihood, a correct parameter count, consistently computed criteria, and the fitting details needed to reproduce the result.

Check Your Understanding

Use these questions to check whether you can design and audit a candidate-model search.

  1. For a \(G\)-component mixture in \(p\) dimensions, explain why the unrestricted mixing weights contribute \(G-1\) parameters rather than \(G\).
  2. For \(p=4\), how many covariance parameters are contributed by EEE and by VVV when \(G=3\)?
  3. Two candidates differ by five parameters and are fitted to \(n=100\) observations. What likelihood gain \(\Delta\) makes the larger candidate's BIC lower?
  4. Why should the search retain results from multiple starts even when it reports only one log-likelihood per candidate?
  5. Name two fitting details that should be documented to make a candidate search reproducible.
Clinical Biostats