Gaussian Graphical Models with Latent Clustering Structure • normalblockr

Description

The Normal-Block model¹ is a Gaussian graphical model with a latent clustering structure, designed for the multivariate analysis of continuous data: it clusters variables and, building on the graphical lasso, infers a network of statistical dependencies between clusters rather than between individual variables. This package implements an efficient (variational) EM algorithm to fit it, accompanied by a set of functions for model selection, visualization and diagnostic. See all the dedicated vignettes for a comprehensive introduction.

normalblockr covers the following model variants, all built around the same normal_block()/NormalBlockData interface:

Known clusters: the grouping of variables is given (e.g. from prior knowledge); only the association network between clusters is estimated, by EM.
Unknown clusters: the grouping is inferred jointly with everything else by a variational EM, either for a single number of clusters or over a range explored as a collection, with model selection via BIC/EBIC/ICL.
Sparse (graphical-lasso) network: an $\ell_1$ penalty² on the inter-cluster precision matrix, for a single penalty value or a path explored as a collection.
Zero-inflated extension: an excess-of-zeros layer for data (e.g. abundance/biomass tables) with more exact zeros than a plain Normal model can represent.

Any combination of these is reached through the same normal_block() function – known or unknown clustering, sparse or not, zero-inflated or not are independent choices, not separate model classes to learn.

Installation

normalblockr is not (yet) on CRAN; install the development version from GitHub:

# install.packages("pak")
pak::pak("jchiquet/normalblockr")

Illustration

We illustrate the known-/unknown-clusters and sparse variants on brca_rppa³: reverse-phase protein array measurements of 163 proteins across 346 breast cancer tumor samples from The Cancer Genome Atlas, together with each sample’s PAM50 molecular subtype.

library(normalblockr)
data(brca_rppa)
Y    <- as.matrix(brca_rppa$expr)
X    <- model.matrix(~ 0 + PAM50_SUBTYPE, data = brca_rppa$covariates)
data <- NormalBlockData$new(Y, X)

Known clusters

A simple, fully data-driven grouping (hierarchical clustering on the raw expression profiles, with no reference to the Normal-Block model) handed to normal_block() as fixed – only the network between the 6 blocks is estimated.

hc    <- brca_rppa$expr |> scale() |> t() |> dist() |> hclust("ward.D2")
group <- cutree(hc, 6) |> normalblockr:::as_indicator()
m_known <- normal_block(data, blocks = group)

Fitting a diagonal normal-block model with fixed blocks 

DONE

m_known

A diagonal normal-block model with fixed blocks .
===========================================================================
 nb_param q n_edges sparsity    loglik deviance      BIC      ICL   EBIC niter
      999 6      15        0 -70387.84 140775.7 146616.3 144073.3 146670    14
===========================================================================
* Useful fields
    $model_par, $posterior_par / $var_par, $clustering 
    $loglik, $BIC, $ICL, $objective, $nb_param, $criteria
* Useful S3 methods
    print(), coef(), sigma(), fitted(), predict()

Unknown number of clusters

Leave the clustering for the model to infer, over a range of candidate cluster counts explored as a collection, then select by ICL.

m_unknown <- normal_block(data, blocks = 2:10)

Fitting a diagonal normal-block model with unknown q 
     number of blocks = 2           
     number of blocks = 3           
     number of blocks = 4           
     number of blocks = 5           
     number of blocks = 6           
     number of blocks = 7           
     number of blocks = 8           
     number of blocks = 9           
     number of blocks = 10           
DONE

m_unknown$plot(c("deviance", "BIC", "ICL", "EBIC"))

m_unknown$get_best_model("ICL")

A diagonal normal-block model with 10 unknown blocks .
===========================================================================
 nb_param  q n_edges sparsity    loglik deviance      BIC      ICL     EBIC
     1042 10      45        0 -68179.63 136359.3 142451.2 139402.8 142658.5
 niter
    16
===========================================================================
* Useful fields
    $model_par, $posterior_par / $var_par, $clustering 
    $loglik, $BIC, $ICL, $objective, $nb_param, $criteria
* Useful S3 methods
    print(), coef(), sigma(), fitted(), predict()

Sparse network

Treating a clustering as fixed (known, or selected above), explore a path of graphical-lasso penalties on the inter-cluster network and select by BIC. Unlike the dense, unpenalized network, this is the one actually worth visualizing.

m_sparse <- normal_block(data, blocks = group, sparsity = TRUE)

Fitting a Collection of diagonal normal-block models with fixed blocks, with different sparsity penalties. 
     penalty = 0.2565018           
     penalty = 0.2188391           
     penalty = 0.1867065           
     penalty = 0.1592919           
     penalty = 0.1359028           
     penalty = 0.1159479           
     penalty = 0.098923           
     penalty = 0.08439792           
     penalty = 0.07200559           
     penalty = 0.06143286           
     penalty = 0.05241254           
     penalty = 0.04471669           
     penalty = 0.03815085           
     penalty = 0.03254907           
     penalty = 0.02776982           
     penalty = 0.02369232           
     penalty = 0.02021353           
     penalty = 0.01724553           
     penalty = 0.01471333           
     penalty = 0.01255294           
     penalty = 0.01070977           
     penalty = 0.009137229           
     penalty = 0.00779559           
     penalty = 0.006650947           
     penalty = 0.005674374           
     penalty = 0.004841193           
     penalty = 0.004130351           
     penalty = 0.003523882           
     penalty = 0.003006463           
     penalty = 0.002565018           
DONE

sp_best  <- m_sparse$get_best_model("BIC")
sp_best$plot_network()

Zero-inflated data

For data with an excess of exact zeros beyond what a plain Normal model would represent (e.g. abundance/biomass tables), zero_inflation = TRUE adds a per-variable excess-of-zero layer. Illustrated on onema⁴: total biomass of 46 fish species across 399 French stream electrofishing stations, where seven in ten entries are exact zeros.

data(onema)
X_zi    <- model.matrix(~ 1 + temperature_med, data = onema$covariates)
Y_zi    <- log(1 + onema$biomass)
data_zi <- NormalBlockData$new(Y_zi, X_zi)
m_zi    <- normal_block(data_zi, blocks = 2:8, zero_inflation = TRUE)

Fitting a diagonal normal-block model with unknown q 
     number of blocks = 2           
     number of blocks = 3           
     number of blocks = 4           
     number of blocks = 5           
     number of blocks = 6           
     number of blocks = 7           
     number of blocks = 8           
DONE

m_zi$get_best_model("ICL")

A zero-inflated diagonal normal-block model with 6 unknown blocks .
===========================================================================
 nb_param q n_edges sparsity    loglik deviance      BIC      ICL     EBIC
      210 6      15        0 -13519.35 27038.71 28296.39 27090.47 28350.14
 niter
    20
===========================================================================
* Useful fields
    $model_par, $posterior_par / $var_par, $clustering 
    $loglik, $BIC, $ICL, $objective, $nb_param, $criteria
* Useful S3 methods
    print(), coef(), sigma(), fitted(), predict()

Learning more

normal-block (vignette("normal-block")): a general introduction on simulated data – known/unknown clusters, sparsity, zero-inflation.
zero-inflated-normal-block (vignette("zero-inflated-normal-block")): a full worked example on onema, with the math behind the zero-inflation extension.
breast-cancer-proteomics (vignette("breast-cancer-proteomics")): a full worked example on brca_rppa – known vs. inferred clustering, model selection, the post-hoc refine() step, sparsifying the selected network.
inst/normal_block_models.qmd in the package sources: a reference card with every model variant’s exact criterion and E/M (or VE/M) updates, as implemented in src/.

normalblockr: Gaussian graphical models with latent clustering structure for multivariate continuous data