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


<!-- badges: start -->

[![R-CMD-check](https://github.com/jchiquet/normalblockr/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/jchiquet/normalblockr/actions/workflows/R-CMD-check.yaml)
[![Coverage
status](https://codecov.io/gh/jchiquet/normalblockr/branch/master/graph/badge.svg)](https://app.codecov.io/github/jchiquet/normalblockr?branch=master)
[![Lifecycle:
experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html)
[![License: GPL
v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
[![](https://img.shields.io/github/last-commit/jchiquet/normalblockr.svg)](https://github.com/jchiquet/normalblockr/commits/master)
<!-- badges: end -->

## Description

> The Normal-Block model[^1] 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](https://jchiquet.github.io/normalblockr/articles/) 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[^2] 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](https://github.com/jchiquet/normalblockr):

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

## Illustration

We illustrate the known-/unknown-clusters and sparse variants on
`brca_rppa`[^3]: 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.

``` r
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.

``` r
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

``` r
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.

``` r
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

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

![](man/figures/README-unknown-1.png)

``` r
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.

``` r
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

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

![](man/figures/README-sparse-1.png)

### 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`[^4]:
total biomass of 46 fish species across 399 French stream electrofishing
stations, where seven in ten entries are exact zeros.

``` r
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

``` r
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`](https://jchiquet.github.io/normalblockr/articles/normal-block.html)
  (`vignette("normal-block")`): a general introduction on simulated data
  – known/unknown clusters, sparsity, zero-inflation.
- [`zero-inflated-normal-block`](https://jchiquet.github.io/normalblockr/articles/zero-inflated-normal-block.html)
  (`vignette("zero-inflated-normal-block")`): a full worked example on
  `onema`, with the math behind the zero-inflation extension.
- [`breast-cancer-proteomics`](https://jchiquet.github.io/normalblockr/articles/breast-cancer-proteomics.html)
  (`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/`.

## References

[^1]: Tous, J., & Chiquet, J. (2026). An integrated method for
    clustering and association network inference. Computational
    Statistics & Data Analysis, 219, 108347.
    [doi:10.1016/j.csda.2026.108347](https://doi.org/10.1016/j.csda.2026.108347)

[^2]: Friedman, J., Hastie, T. and Tibshirani, R. Sparse inverse
    covariance estimation with the graphical lasso. Biostatistics, 9(3),
    2008.

[^3]: Cancer Genome Atlas Network. Comprehensive molecular portraits of
    human breast tumours. Nature, 490, 2012, 61-70.
    [doi:10.1038/nature11412](https://doi.org/10.1038/nature11412)

[^4]: Danet, A., Mouchet, M., Bonnaffé, W., Thébault, E., Fontaine, C.
    Species richness and food-web structure jointly drive community
    biomass and its temporal stability in fish communities, data set,
    2021, Zenodo.
    [doi:10.5281/zenodo.5095656](https://doi.org/10.5281/zenodo.5095656)
