Package {normalblockr}

Type: Package
Title: Gaussian Graphical Models with Latent Clustering Structure
Version: 0.2.1
Description: Implements the Normal-Block model, a Gaussian graphical model with a latent clustering structure for the multivariate analysis of continuous data. The model clusters variables and, building on the graphical lasso, infers a network of statistical dependencies between clusters rather than between individual variables, for known or unknown clusterings, with an optional zero-inflation extension for data with an excess of exact zeros. See Tous & Chiquet (2026) <doi:10.1016/j.csda.2026.108347> for the model itself and its variational expectation-maximization estimation procedure.
License: GPL (≥ 3)
URL: https://github.com/jchiquet/normalblockr
BugReports: https://github.com/jchiquet/normalblockr/issues
Encoding: UTF-8
LazyData: true
Depends: R (≥ 4.1.0)
Imports: glassoFast, Matrix, dplyr, tidyr, tibble, purrr, R6, ggplot2, igraph, sbm, corrplot, scales, MASS, stats, Rcpp
LinkingTo: Rcpp, RcppArmadillo
Suggests: testthat (≥ 3.0.0), aricode, covr, Metrics, rmarkdown, pheatmap, paletteer, ggthemes, knitr
Config/testthat/edition: 3
VignetteBuilder: knitr
Config/roxygen2/version: 8.0.0
RoxygenNote: 7.3.3
NeedsCompilation: yes
Packaged: 2026-08-20 11:30:04 UTC; jchiquet
Author: Jeanne Tous [aut], Julien Chiquet ORCID iD [aut, cre]
Maintainer: Julien Chiquet <julien.chiquet@inrae.fr>
Repository: CRAN
Date/Publication: 2026-09-03 12:00:02 UTC

normalblockr: Gaussian Graphical Models with Latent Clustering Structure

Description

Implements the Normal-Block model, a Gaussian graphical model with a latent clustering structure for the multivariate analysis of continuous data. The model clusters variables and, building on the graphical lasso, infers a network of statistical dependencies between clusters rather than between individual variables, for known or unknown clusterings, with an optional zero-inflation extension for data with an excess of exact zeros. See Tous & Chiquet (2026) doi:10.1016/j.csda.2026.108347 for the model itself and its variational Expectation-Maximization estimation procedure.

Author(s)

Maintainer: Julien Chiquet julien.chiquet@inrae.fr (ORCID)

Authors:

See Also

Useful links:

Bayesian Information Criterion for a Normal-Block Model

Description

Extracts the (variational) BIC of a fitted normal-block model, computed as 'deviance + log(n) * nb_param' (lower is better).

Usage

## S3 method for class 'NormalBlockVarBase'
BIC(object, ...)

Arguments

object

An object of class NormalBlockVarBase.

...

not used, only here for S3 compatibility

Value

A scalar: the (variational) BIC.

Examples

ex_data <- generate_normal_block_var_data(n = 50, p = 20, d = 1, q = 3)
data <- NormalBlockData$new(ex_data$Y, ex_data$X)
model <- normal_block(data, blocks = 3, control = NB_control(verbose = FALSE))
BIC(model)

Bayesian Information Criterion for a Collection of Normal-Block Models

Description

Returns the (variational) BIC of every model in a collection of normal-block models.

Usage

## S3 method for class 'NormalBlockVarCollection'
BIC(object, ...)

Arguments

object

An object inheriting from NormalBlockVarCollection.

...

not used, only here for S3 compatibility

Value

A numeric vector of BIC values, one per model in the collection.

Examples

ex_data <- generate_normal_block_var_data(n = 50, p = 20, d = 1, q = 3)
data <- NormalBlockData$new(ex_data$Y, ex_data$X)
models <- normal_block(data, blocks = 2:5, control = NB_control(verbose = FALSE))
BIC(models)

NB_control

Description

Control the model settings and various optimization parameters

Usage

NB_control(
  niter = 500,
  threshold = 1e-04,
  sparsity_weights = NULL,
  sparsity_penalties = NULL,
  n_sparsity_penalties = 30,
  min_ratio = 0.01,
  fixed_tau = FALSE,
  clustering_init = "ward2",
  verbose = TRUE,
  heuristic = FALSE,
  noise_covariance = c("diagonal", "spherical"),
  refine = FALSE
)

Arguments

niter

number of iterations in model optimization

threshold

loglikelihood / elbo threshold under which optimization stops

sparsity_weights

weights with which the penalty should be applied in case sparsity is required, non-0 values on the diagonal mean diagonal shall be penalized too (default is non-penalized diagonal and 1s off-diagonal)

sparsity_penalties

list of penalties the user wants to test, other parameters are only used if penalties is not specified

n_sparsity_penalties

number of penalties to test.

min_ratio

ratio for sparsity between max penalty (0 edge penalty) and min penalty to test

fixed_tau

whether tau should be fixed at clustering_init during optimization useful for calls to fixed_q models in stability_selection

clustering_init

how to obtain the initial clustering of the q unknown blocks: a heuristic name ("ward2", the default, "kmeans", "sbm" or "spectral"), an actual clustering (a vector of labels or a p x q indicator matrix, or a list of either per q for a collection), or "best_of_inits" to try several heuristics per model and keep the best-ELBO fit (see [NormalBlockVarBase]'s 'best_of_inits()'; not supported with 'sparsity = TRUE'). See 'inst/methods_initialization_and_refine.md' for the heuristics' rationale, why no single one dominates, and how this interacts with 'refine' (below).

verbose

telling if information should be printed during optimization

heuristic

whether to use the heuristic approach (moment-based, no (V)EM recursion) instead of the full (V)EM. Default is FALSE. In heuristic mode, no likelihood/ELBO is computed, so 'entropy', 'loglik', 'BIC', 'ICL' and 'EBIC' are all 'NA' on the resulting model.

noise_covariance

variance can be variable specific ("diagonal", the default) or common ("spherical")

refine

for [NormalBlockVarCollectionClusters] only: whether 'optimize()' should automatically call 'refine()' afterwards. Default 'FALSE' since it adds real cost; call 'collection$refine()' directly at any point afterwards for the same effect without setting this.

Value

A named list of parameters to pass to [normal_block()]'s 'control' argument.

Data Container for Normal-Block Models

Description

R6 class holding the responses and design matrix used to fit a normal-block model.

Public fields

Y

the matrix of responses (rescaled column-wise if 'scale = TRUE')

Y_scale

the per-column standard deviation Y was divided by (all 1's if 'scale = FALSE')

X

the matrix of covariates

X0

the matrix of zero-inflation covariates, if applicable

formula

describes the relationship between Y and X, and X0 if applicable, useful if not all of X's or X0's covariates should be used, should be formatted ~ X1 + X2... | Z1 + Z2... with the Normal formula before the | and the ZI formula after the |

n

sample size

d

number of covariates

d0

number of zero-inflation covariates, if applicable

p

number of variables

XtXm1

inverse of XtX, useful for inference

XtY

useful for inference

npY

total number of non zeros in Y

nY

total number of non zeros for each column/variable in Y

zeros

where are the zero in Y

zeros_bar

where are the non-zeros in Y

Methods

Public methods

NormalBlockData$new()

Create a new ['NormalBlockData'] object.

Usage
NormalBlockData$new(Y, X, X0 = NULL, formula = NULL, scale = TRUE)
Arguments
Y

the matrix of responses (called Y in the model).

X

design matrix (called X in the model).

X0

zero-inflation design matrix, if applicable.

formula

describes the relationship between Y and X, useful if not all of X's covariates should be used.

scale

whether to rescale each column of Y by its own standard deviation (no centering). Default TRUE – see the class-level documentation for the rationale and its limits.

NormalBlockData$clone()

The objects of this class are cloneable with this method.

Usage
NormalBlockData$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Base Class for Normal-Block Models

Description

R6 abstract class for a generic sparse Normal-Block model.

Public fields

data

object of NormalBlockData class, with responses and design matrix

Active bindings

inference_method

inference procedure used (heuristic or integrated with EM)

n

number of samples

p

number of responses per sample

d

number of variables (dimensions in X)

d0

number of zi variables (dimensions in X0)

q

number of blocks

n_edges

number of edges of the network (non null coefficient of the sparse precision matrix Omegaq)

model_par

a list with the matrices of the model parameters: B (covariates), dm1 (species variance), Omegaq (groups precision matrix)). On the internal fitting scale (‘self$data$Y', possibly column-rescaled by 'NormalBlockData(scale = TRUE)') – use '$B_original'/'$dm1_original' for the same quantities converted back to Y’s original units.

B_original

regression coefficients (d x p), converted back to Y's original units (undoing ‘NormalBlockData(scale = TRUE)'’s column-wise rescaling, if any). Use 'model_par$B' instead for the coefficients on the internal fitting scale.

dm1_original

inverse residual variance per variable (1 / Var(Y_j)), converted back to Y's original units. Use 'model_par$dm1' instead for the internal fitting scale. With 'noise_covariance = "spherical"', 'model_par$dm1' is a single value repeated p times (one shared variance on the fitting scale); once converted back per-variable, the p values returned here generally differ from one another whenever Y's columns were rescaled by different factors – correctly so, since a single shared *scaled* variance does not correspond to a single shared variance in the original, heterogeneous-scale units.

nb_param

number of parameters in the model

objective

evolution of the objective function during (V)EM algorithm

loglik

(or its variational lower bound)

deviance

(or its variational lower bound)

BIC

(or its variational lower bound)

entropy

Entropy of the conditional distribution when applicable

ICL

variational lower bound of the ICL

EBIC

variational lower bound of the EBIC

criteria

a vector with loglik, BIC and number of parameters

sparsity

(overall sparsity parameter)

sparsity_weights

(weights associated to each pair of groups)

sparsity_term

(sparsity_term term in log-likelihood due to sparsity)

get_res_covariance

whether the residual covariance is diagonal or spherical

memberships

cluster memberships

clustering

given as the list of elements contained in each cluster

cluster_sizes

given as a vector of cluster sizes

elements_per_cluster

given as the list of elements contained in each cluster

Methods

Public methods

NormalBlockVarBase$new()

Create a new ['NormalBlockVarBase'] object.

Usage
NormalBlockVarBase$new(
  data,
  q,
  sparsity = 0,
  control = NB_control(),
  zero_inflation = FALSE
)
Arguments
data

object of NormalBlockData class, with responses and design matrix

q

number of block/cluster

sparsity

sparsity penalty on the network density

control

structured list of more specific parameters, to generate with NB_control

zero_inflation

whether the concrete subclass models zero-inflation; set by the ZI subclasses themselves, not meant to be set by the end user. When 'FALSE', the (costly) zero-inflation probability fit ('kappa'/'B0') is skipped entirely, since it would otherwise never be used downstream.

Returns

A new ['NormalBlockVarBase'] object

NormalBlockVarBase$update()

Update a ['NormalBlockVarBase'] object

All possible parameters of the child classes

Usage
NormalBlockVarBase$update(
  B = NA,
  dm1 = NA,
  C = NA,
  Omegaq = NA,
  gamma = NA,
  mu = NA,
  kappa = NA,
  alpha = NA,
  M = NA,
  S = NA,
  ll_list = NA,
  warm_started = NA,
  clustering_init = NA
)
Arguments
B

regression matrix

dm1

diagonal vector of inverse variance matrix (variables level)

C

the matrix of groups memberships (posterior probabilities)

Omegaq

groups inverse variance matrix

gamma

variance of posterior distribution of W

mu

mean for posterior distribution of W

kappa

vector of zero-inflation probabilities

alpha

vector of groups probabilities

M

variational mean for posterior distribution of W

S

variational diagonal of variances for posterior distribution of W

ll_list

list of log-lik (elbo) values

warm_started

whether 'EM_initialize()' should treat the model as already initialized (reuse B/Omegaq/dm1/C/alpha/M/S as they stand) rather than recomputing a fresh heuristic initialization – set by [warm_start_from()] and by [split()]/[merge()].

clustering_init

name of a clustering heuristic to switch to, re-derived at the next 'optimize()' call instead of reusing the current state (see 'NB_control(clustering_init = )'). Used by [best_of_inits()].

Returns

Update the current ['normal'] object

NormalBlockVarBase$best_of_inits()

Try several clustering-initialization heuristics and keep the best-ELBO converged fit (see 'NB_control(clustering_init = )' and 'inst/methods_initialization_and_refine.md' for the rationale). Every candidate is first screened with a short 'trial_niter' run (same idea as 'candidates_split()'/'candidates_merge()'), and only the 'max_training' best-screened ones are fully retrained with 'control'.

Usage
NormalBlockVarBase$best_of_inits(
  inits = c("ward2", "kmeans", "spectral"),
  trial_niter = 10,
  max_training = 2,
  control = list(niter = 500, threshold = 1e-04)
)
Arguments
inits

vector of clustering-heuristic names to try

trial_niter

number of (V)EM iterations used to cheaply screen every candidate in 'inits' before fully retraining the best few

max_training

how many of the screened candidates (best 'loglik' after 'trial_niter' iterations) get fully retrained with 'control'

control

'optimize()' control list ('niter'/'threshold') used for the final full retraining of the 'max_training' best candidates

Returns

a new, already-optimized ['NormalBlockVarBase'] object. Does not mutate the current object; reassign the result ('model <- model$best_of_inits()').

NormalBlockVarBase$optimize()

calls optimization (EM or heuristic) and updates relevant fields

Usage
NormalBlockVarBase$optimize(
  control = list(niter = 500, threshold = 1e-04),
  warn = TRUE
)
Arguments
control

a list for controlling the optimization process

warn

whether to warn when the (V)EM stops at the 'niter' cap without reaching 'threshold' (see 'private$warn_if_not_converged()'). Set to 'FALSE' for deliberately-truncated trial fits (cheap candidate scoring in 'candidates_split()'/'candidates_merge()', the sparsity-path warm-start probe in [NormalBlockVarCollectionSparsity]) where stopping at the cap is expected and not a sign of trouble.

Returns

optimizes the model and updates its parameters

NormalBlockVarBase$warm_start_from()

Seed this model's starting parameters from another, already-optimized model with the same q, instead of a fresh heuristic clustering. Used by [NormalBlockVarCollectionSparsity] to warm-start each penalty in a sparsity path from the previous one's solution. 'B0'/'kappa' (zero-inflation) are left untouched: they depend only on the data, already set correctly and independently on every model.

Usage
NormalBlockVarBase$warm_start_from(other)
Arguments
other

a [NormalBlockVarBase] object, already optimized

Returns

Update the current object in place with ‘other'’s parameters

NormalBlockVarBase$split()

Create a clone of the current ['NormalBlockVarBase'] object after splitting cluster 'cl' We split the cluster according to the species variances

Usage
NormalBlockVarBase$split(index, in_place = FALSE)
Arguments
index

index (integer) of the cluster to split

in_place

should the split applied to the object itself, or should a copy be sent? default FALSE (send a copy)

Returns

A new ['NormalBlockVarBase'] object

NormalBlockVarBase$candidates_split()

generate and select a set of candidate models by splitting the clusters of the current model

Usage
NormalBlockVarBase$candidates_split(trial_niter = 5)
Arguments
trial_niter

number of EM iterations used to cheaply score each candidate before [SelectionNClusters] fully re-optimizes the best few (‘train_best_candidates()'’s 'max_training') – kept short on purpose.

NormalBlockVarBase$candidates_merge()

generate and select a set of candidate models by merging the clusters of the current model

Usage
NormalBlockVarBase$candidates_merge(max_candidates = 30, trial_niter = 2)
Arguments
max_candidates

merge candidates are, unlike split's, quadratic in q ('choose(q, q-2)' pairs) – beyond 'max_candidates' pairs, only the most promising ones (largest '|Omegaq[i, j]|', i.e. the most strongly related cluster pairs in the current fit) are actually built and trial-optimized, since merging two nearly independent blocks is rarely competitive anyway. Set to 'Inf' to always try every pair.

trial_niter

see [candidates_split()]

NormalBlockVarBase$merge()

Create a clone of the current ['NormalBlockVarBase'] object after merging clusters 'cl1' and 'cl2'

Usage
NormalBlockVarBase$merge(indices, in_place = FALSE)
Arguments
indices

indices (couple of integer) of the clusters to merge

in_place

should the split applied to the object itself, or should a copy be sent? default FALSE (send a copy)

Returns

A new ['NormalBlockVarBase'] object

NormalBlockVarBase$predict()

Predicts observations Y for new covariates X.

Usage
NormalBlockVarBase$predict(new_X)
Arguments
new_X

new set of covariates.

Returns

A n*p prediction matrix for new observations

NormalBlockVarBase$latent_network()

Extract interaction network in the latent space, as a matrix rather than a plot – see '$plot_network()' to plot it instead.

Usage
NormalBlockVarBase$latent_network(
  type = c("partial_cor", "support", "precision")
)
Arguments
type

edge value in the network. Can be "support" (binary edges), "precision" (coefficient of the precision matrix) or "partial_cor" (partial correlation between species)

Returns

a square matrix of size 'self$q'

NormalBlockVarBase$plot_loglik()

plots the evolution of the objective (log-likelihood or ELBO) across the (V)EM iterations of the last call to 'optimize()'.

Usage
NormalBlockVarBase$plot_loglik(show_increment = TRUE)
Arguments
show_increment

whether to add a second panel with the (log10) absolute increment between iterations and the convergence 'threshold' – distinguishes true convergence from a flat-looking objective trace.

Returns

a ['ggplot2::ggplot'] graph

NormalBlockVarBase$plot_network()

plot the latent network. To extract the network as a matrix instead of plotting it, use '$latent_network()'.

Usage
NormalBlockVarBase$plot_network(
  type = c("partial_cor", "support"),
  output = c("igraph", "corrplot"),
  edge.color = c("#F8766D", "#00BFC4"),
  remove.isolated = FALSE,
  node.labels = NULL,
  layout = igraph::layout_in_circle,
  plot = TRUE
)
Arguments
type

edge value in the network. Either "precision" (coefficient of the precision matrix) or "partial_cor" (partial correlation between species).

output

Output type. Either 'igraph' (for the network) or 'corrplot' (for the adjacency matrix)

edge.color

Length 2 color vector. Color for positive/negative edges. Default is 'c("#F8766D", "#00BFC4")'. Only relevant for igraph output.

remove.isolated

if 'TRUE', isolated node are remove before plotting. Only relevant for igraph output.

node.labels

vector of character. The labels of the nodes. The default will use the column names ot the response matrix.

layout

an optional igraph layout. Only relevant for igraph output.

plot

logical. Should the final network be displayed or only sent back to the user. Default is 'TRUE'.

NormalBlockVarBase$plot()

plots the evolution of the objective during model optimization (see 'plot_loglik()')

Usage
NormalBlockVarBase$plot()

NormalBlockVarBase$print()

User friendly print method

Usage
NormalBlockVarBase$print(model = paste("A", self$who_am_I, ".\n"))
Arguments
model

First line of the print output

NormalBlockVarBase$clone()

The objects of this class are cloneable with this method.

Usage
NormalBlockVarBase$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Examples

# An internal abstract base class, never instantiated directly -- see
# normal_block() for how concrete models (NormalBlockVarKnownClusters,
# NormalBlockVarUnknownClusters, and their zero-inflated variants) are
# actually created and fitted.

Base Class for a Collection of Normal-Block Models

Description

Shared scaffolding for the collections explored by [get_model()]/ [normal_block()]: a sweep over sparsity penalties (['NormalBlockVarCollectionSparsity']), over the number of clusters (['NormalBlockVarCollectionClusters']), or over both (['NormalBlockVarCollectionClustersSparsity']). Concrete subclasses set 'private$progress_field'/'private$progress_label' in their 'initialize()' and provide their own 'get_best_model()', delegating the (row of 'self$criteria' minimizing a criterion) lookup to 'private$best_id()'.

Public fields

models

list of models (or sub-collections) explored by the collection

control

store the list of user-defined model settings and optimization parameters

Active bindings

criteria

a data frame with the values of some criteria for the collection of models

loglik

not defined for a collection (which of its models?) – accessing it raises an informative error instead of silently returning ‘NULL'. Use 'logLik()' for every model’s log-likelihood, or '$get_best_model()$loglik' for a single one.

Methods

Public methods

NormalBlockVarCollection$optimize()

optimizes every model (or sub-collection) in the collection

Usage
NormalBlockVarCollection$optimize(
  control = list(niter = 500, threshold = 1e-04, verbose = TRUE)
)
Arguments
control

optimization parameters (niter and threshold). When 'control$clustering_init' is '"best_of_inits"', each leaf model is fit via its own 'best_of_inits()' instead of a plain 'optimize()'.

NormalBlockVarCollection$print()

User-friendly print method: model type and the range of q/sparsity explored. See 'summary()' for the full criteria table.

Usage
NormalBlockVarCollection$print()

NormalBlockVarCollection$summary()

Summarize the collection: model type, full criteria table, and the range of q/sparsity explored.

Usage
NormalBlockVarCollection$summary()
Returns

An object of class 'summary.NormalBlockVarCollection', printed with a dedicated [print.summary.NormalBlockVarCollection()] method.

NormalBlockVarCollection$clone()

The objects of this class are cloneable with this method.

Usage
NormalBlockVarCollection$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Examples

# An internal abstract base class, never instantiated directly -- see
# normal_block() for how collections (NormalBlockVarCollectionClusters,
# NormalBlockVarCollectionSparsity, NormalBlockVarCollectionClustersSparsity)
# are actually created and fitted.

Collection of Normal-Block Models over a Range of Cluster Counts

Description

R6 class for a collection of normal-block models with different number of clusters (q) and a fixed sparsity level.

Super class

NormalBlockVarCollection -> NormalBlockVarCollectionClusters

Active bindings

q_list

number of blocks

who_am_I

a method to print what model is being fitted

Methods

Public methods

Inherited methods

NormalBlockVarCollectionClusters$new()

Create a new ['NormalBlockVarCollectionClusters'] object.

Usage
NormalBlockVarCollectionClusters$new(
  mydata,
  q_list,
  zero_inflation = FALSE,
  sparsity = 0,
  control = NB_control()
)
Arguments
mydata

object of NormalBlockData class, with responses and design matrix

q_list

list of q values (number of groups) in the collection

zero_inflation

whether the models in the collection should be zero-inflated or not

sparsity

sparsity penalty on the network density

control

structured list of more specific parameters, to generate with NB_control

Returns

A new ['NormalBlockVarCollectionClusters'] object

NormalBlockVarCollectionClusters$get_model()

returns the NormalBlockVarUnknownClusters model corresponding to given q

Usage
NormalBlockVarCollectionClusters$get_model(q)
Arguments
q

number of blocks asked by user

Returns

A NormalBlockVarUnknownClusters object with given value q

NormalBlockVarCollectionClusters$get_best_model()

Extract best model in the collection

Usage
NormalBlockVarCollectionClusters$get_best_model(
  crit = c("ICL", "BIC", "EBIC", "deviance")
)
Arguments
crit

a character for the criterion used to performed the selection. Either "ICL" or "BIC". "ICL" is the default criterion

Returns

a ['NormalBlockVarUnknownClusters'] object

NormalBlockVarCollectionClusters$plot()

Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of network fits (a ['Networkfamily'])

Usage
NormalBlockVarCollectionClusters$plot(
  criteria = c("deviance", "ICL", "BIC", "EBIC")
)
Arguments
criteria

vector of characters. The criteria to plot in 'c("deviance", "BIC", "ICL")'. Defaults to all of them.

Returns

a ['ggplot2::ggplot'] graph

NormalBlockVarCollectionClusters$optimize()

optimizes every model in the collection, then – if 'control$refine' is 'TRUE' (see [NB_control()], default 'FALSE') – calls [refine()] automatically.

Usage
NormalBlockVarCollectionClusters$optimize(
  control = list(niter = 500, threshold = 1e-04, verbose = TRUE)
)
Arguments
control

optimization parameters (niter, threshold, verbose)

NormalBlockVarCollectionClusters$refine()

Tries to improve every model in the collection with a short split-and-reoptimize trial seeded from its smaller-q neighbor ('"split"'), a short merge-and-reoptimize trial seeded from its larger-q neighbor ('"merge"'), or both (the default); a candidate replaces the original only if it strictly lowers the deviance, so this can only improve (or leave unchanged) each model it touches. Only contiguous q pairs ('q' and 'q -/+ 1', both present in the collection) are refined. See 'inst/methods_initialization_and_refine.md' for the rationale and empirical evidence.

Usage
NormalBlockVarCollectionClusters$refine(
  trial_niter = 2,
  max_candidates = 30,
  directions = c("split", "merge"),
  verbose = self$control$verbose
)
Arguments
trial_niter

number of EM iterations used for the cheap trial candidates (passed to 'candidates_split()'/'candidates_merge()') before fully re-optimizing only the best one.

max_candidates

passed to 'candidates_merge()' (ignored for '"split"', which is never combinatorial in q) – see its documentation.

directions

which neighbor(s) to seed refinement candidates from: '"split"' (smaller-q neighbor), '"merge"' (larger-q neighbor), or both (the default).

verbose

whether to print, for each q attempted, whether the candidate from that neighbor improved on it. Defaults to 'control$verbose' (the value set at construction, see [NB_control()]).

Returns

invisibly returns 'self'; improved models replace the originals in '$models' in place.

NormalBlockVarCollectionClusters$clone()

The objects of this class are cloneable with this method.

Usage
NormalBlockVarCollectionClusters$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Collection of Normal-Block Models over Cluster Counts and Sparsity Levels

Description

R6 class for a collection of normal-block models with different number of clusters (q) and different sparsity levels.

Super class

NormalBlockVarCollection -> NormalBlockVarCollectionClustersSparsity

Active bindings

q_list

number of blocks

sparsity

list of penalties used for each q

who_am_I

a method to print what model is being fitted

Methods

Public methods

Inherited methods

NormalBlockVarCollectionClustersSparsity$new()

Create a new ['NormalBlockVarCollectionClustersSparsity'] object.

Usage
NormalBlockVarCollectionClustersSparsity$new(
  mydata,
  q_list,
  zero_inflation = FALSE,
  control = NB_control()
)
Arguments
mydata

object of NormalBlockData class, with responses and design matrix

q_list

list of q values (number of groups) in the collection

zero_inflation

boolean to specify whether data is zero-inflated

control

structured list of parameters to handle sparsity control

Returns

A new ['NormalBlockVarCollectionClustersSparsity'] object

NormalBlockVarCollectionClustersSparsity$get_model()

returns a collection of models corresponding to given q or one single model if penalty is also given

Usage
NormalBlockVarCollectionClustersSparsity$get_model(q, sparsity = NA)
Arguments
q

number of blocks asked by user.

sparsity

sparsity penalty asked by user

Returns

either a NormalBlockVarCollectionSparsity or a NormalBlockVarUnknownClusters object

NormalBlockVarCollectionClustersSparsity$get_best_model()

Extract best model in the collection

Usage
NormalBlockVarCollectionClustersSparsity$get_best_model(
  crit = c("ICL", "BIC", "EBIC")
)
Arguments
crit

a character for the criterion used to performed the selection. Either "BIC", "EBIC" or "ICL". "ICL" is the default criterion

Returns

a ['NormalBlockVarUnknownClusters'] object

NormalBlockVarCollectionClustersSparsity$plot()

Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of network fits (a ['Networkfamily'])

Usage
NormalBlockVarCollectionClustersSparsity$plot(
  criterion = c("deviance", "ICL", "BIC", "EBIC"),
  n_intervals = NULL
)
Arguments
criterion

The criteria to plot in 'c("deviance", BIC", "EBIC", "ICL")'. Defaults deviance.

n_intervals

number of intervals into which the penalties range should be split

Returns

a ['ggplot'] heatmap

NormalBlockVarCollectionClustersSparsity$clone()

The objects of this class are cloneable with this method.

Usage
NormalBlockVarCollectionClustersSparsity$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Collection of Normal-Block Models over a Sparsity Path

Description

R6 class for a collection of normal-block models with a fixed clustering (blocks) and different sparsity levels.

Super class

NormalBlockVarCollection -> NormalBlockVarCollectionSparsity

Public fields

data

object of NormalBlockData class, with responses and design matrix

Active bindings

q

number of blocks

blocks

group matrix or number of blocks.

sparsity

list of sparsity penalties

sparsity_details

list of information about model's penalties

criteria

a data frame with the values of some criteria ((approximated) log-likelihood, BIC) for the collection of models

stability_path

measure of edges stability based on StARS method

stability

mean edge stability along the sparsity penalties path

who_am_I

a method to print what model is being fitted

Methods

Public methods

Inherited methods

NormalBlockVarCollectionSparsity$new()

Create a new ['NormalBlockVarCollectionSparsity'] object.

Usage
NormalBlockVarCollectionSparsity$new(
  mydata,
  blocks,
  zero_inflation = FALSE,
  control = NB_control()
)
Arguments
mydata

object of NormalBlockData class, with responses and design matrix

blocks

either a clustering matrix (known, fixed clustering) or a single integer (number of blocks to infer)

zero_inflation

boolean to specify whether data is zero-inflated

control

structured list of parameters to handle sparsity control

Returns

A new ['NormalBlockVarCollectionSparsity'] object

NormalBlockVarCollectionSparsity$optimize()

optimizes every model in the sparsity path, warm-starting each one (after the first) from the previous, adjacent penalty's converged parameters (see [NormalBlockVarBase]'s 'warm_start_from()') instead of re-deriving everything from the heuristic clustering, the way the generic [NormalBlockVarCollection] 'optimize()' would. 'blocks' (hence q) is fixed across the whole path, only the sparsity penalty changes, so the warm start is always between models of matching shape.

Usage
NormalBlockVarCollectionSparsity$optimize(
  control = list(niter = 500, threshold = 1e-04, verbose = TRUE)
)
Arguments
control

optimization parameters (niter and threshold)

NormalBlockVarCollectionSparsity$get_model()

returns the NormalBlockVarKnownClusters model corresponding to given penalty

Usage
NormalBlockVarCollectionSparsity$get_model(sparsity)
Arguments
sparsity

sparsity penalty asked by user

Returns

A NormalBlockVarKnownClusters (sparse) object with given value penalty

NormalBlockVarCollectionSparsity$get_best_model()

Extract best model in the collection

Usage
NormalBlockVarCollectionSparsity$get_best_model(
  crit = c("BIC", "EBIC", "ICL", "StARS"),
  stability = 0.9
)
Arguments
crit

a character for the criterion used to performed the selection.

stability

if criterion = "StARS" gives level of stability required. Either "BIC", "EBIC", "ICL" or "StARS". Default is BIC

Returns

a ['NormalBlockVarUnknownClusters'] object

NormalBlockVarCollectionSparsity$plot()

Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of network fits (a ['Networkfamily'])

Usage
NormalBlockVarCollectionSparsity$plot(
  criteria = c("deviance", "BIC", "EBIC", "ICL"),
  log.x = TRUE
)
Arguments
criteria

vector of characters. The criteria to plot in 'c("deviance", BIC", "EBIC", "ICL")'. Defaults to all of them.

log.x

logical: should the x-axis be represented in log-scale? Default is 'TRUE'.

Returns

a ['ggplot'] graph

NormalBlockVarCollectionSparsity$stability_selection()

Compute the stability path by stability selection

Usage
NormalBlockVarCollectionSparsity$stability_selection(
  subsamples = NULL,
  n_subsamples = 10
)
Arguments
subsamples

a list of vectors describing the subsamples. The number of vectors (or list length) determines the number of subsamples used in the stability selection. Automatically set to 20 subsamples with size '10*sqrt(n)' if 'n >= 144' and '0.8*n' otherwise following Liu et al. (2010) recommendations.

n_subsamples

number of subsamples to create if the subsamples are not given

NormalBlockVarCollectionSparsity$clone()

The objects of this class are cloneable with this method.

Usage
NormalBlockVarCollectionSparsity$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Normal-Block Model with Known Clustering

Description

R6 class for a normal-block model with known clustering.

Super class

NormalBlockVarBase -> NormalBlockVarKnownClusters

Active bindings

posterior_par

a list with the parameters of posterior distribution W | Y

entropy

Entropy of the conditional distribution

fitted

Y values predicted by the model, in Y's original units

who_am_I

a method to print what model is being fitted

Methods

Public methods

Inherited methods

NormalBlockVarKnownClusters$new()

Create a new ['NormalBlockVarKnownClusters'] object.

Usage
NormalBlockVarKnownClusters$new(data, C, sparsity = 0, control = NB_control())
Arguments
data

object of NormalBlockVarData class, with responses and design matrix

C

clustering matrix C_jk = 1 if species j belongs to cluster k

sparsity

to apply on variance matrix when calling GLASSO

control

structured list of more specific parameters, to generate with NB_control

Returns

A new ['NormalBlockVarKnownClusters'] object

NormalBlockVarKnownClusters$clone()

The objects of this class are cloneable with this method.

Usage
NormalBlockVarKnownClusters$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Normal-Block Model with Unknown Clustering

Description

R6 class for a normal-block model with a fixed number of clusters (but unknown clustering).

Super class

NormalBlockVarBase -> NormalBlockVarUnknownClusters

Public fields

fixed_tau

whether tau should be fixed at clustering_init during optimization, useful for stability selection

Active bindings

model_par

a list with the matrices of the model parameters: B (covariates), dm1 (species variance), Omegaq (groups precision matrix))

nb_param

number of parameters in the model

var_par

a list with the matrices of the variational parameters: M (means), S (variances), tau (posterior group probabilities)

entropy

Entropy of the conditional distribution

fitted

Y values predicted by the model, in Y's original units

who_am_I

a method to print what model is being fitted

Methods

Public methods

Inherited methods

NormalBlockVarUnknownClusters$new()

Create a new ['NormalBlockVarUnknownClusters'] object.

Usage
NormalBlockVarUnknownClusters$new(
  data,
  q,
  sparsity = 0,
  control = NB_control()
)
Arguments
data

contains the matrix of responses (Y) and the design matrix (X).

q

required number of groups

sparsity

sparsity penalty to add on blocks precision matrix for sparsity

control

structured list for specific parameters

Returns

A new ['NormalBlockVarUnknownClusters'] object

NormalBlockVarUnknownClusters$clone()

The objects of this class are cloneable with this method.

Usage
NormalBlockVarUnknownClusters$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Select the Number of Clusters by Split/Merge Search

Description

R6 class for selecting the number of clusters (q) by forward/backward split-and-merge exploration, comparing models with the ICL.

Kept internal/unexported for reference: empirically (see [NormalBlockVarCollectionClusters]'s 'refine()' documentation) chaining every q from a single starting point this way is no longer recommended over [NormalBlockVarCollectionClusters]'s independent per-q cold starts followed by ‘refine()' – the latter matches or beats this class’s quality at a fraction of the cost on every real dataset tested, and has no equivalent to this class's failure mode (an early bad split/merge propagating to every larger q it explores from there).

Public fields

best_models

list of models explored so far indexed by the cluster sizes

ICL_explored

list of ICL reached so far, indexed by the cluster sizes

n_clusters_range

the range of the cluster sizes to browse

control

store the list of user-defined model settings and optimization parameters

Active bindings

best_model

best model explored so far in term of ICL

Methods

Public methods

SelectionNClusters$new()

Explores and selects the optimal number of classes

Usage
SelectionNClusters$new(
  mydata,
  n_clusters_range,
  zero_inflation = FALSE,
  sparsity = 0,
  control = NB_control()
)
Arguments
mydata

object of NormalBlockData class, with responses and design matrix

n_clusters_range

a size-2 vector with the range of cluster size to browse

zero_inflation

whether the models in the collection should be zero-inflated or not

sparsity

sparsity penalty on the network density

control

structured list of more specific parameters, to generate with NB_control

Returns

The best model in the ICL sense

SelectionNClusters$fit()

perform model selection with forward/backward exploration with split and merge strategy

Usage
SelectionNClusters$fit()

SelectionNClusters$train_best_candidates()

perform model selection with forward/backward exploration with split and merge strategy

Usage
SelectionNClusters$train_best_candidates(model, strategy, max_training = 3)
Arguments
model

a normal-block model

strategy

a character, either "split" or "merge"

max_training

maximal of model fully trained at each step of the exploration

SelectionNClusters$explore_forward()

perform forward exploration with a split strategy

Usage
SelectionNClusters$explore_forward()

SelectionNClusters$explore_backward()

perform backward exploration with a merge strategy

Usage
SelectionNClusters$explore_backward()

SelectionNClusters$plot()

Display the ICL for all the best models explored per number of cluster, and the winner

Usage
SelectionNClusters$plot()
Returns

a ['ggplot2::ggplot'] graph

SelectionNClusters$clone()

The objects of this class are cloneable with this method.

Usage
SelectionNClusters$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Zero-Inflated Normal-Block Model with Known Clustering

Description

R6 class for a zero-inflated normal-block model with a known clustering.

Super class

NormalBlockVarBase -> ZINormalBlockVarKnownClusters

Active bindings

posterior_par

a list with the parameters of posterior distribution W | Y

entropy

Entropy of the conditional distribution

nb_param

number of parameters in the model

model_par

a list with model parameters: B (covariates), dm1 (species variance), Omegaq (groups precision matrix), kappa (zero-inflation probabilities)

fitted

Y values predicted by the model, in Y's original units

who_am_I

a method to print what model is being fitted

Methods

Public methods

Inherited methods

ZINormalBlockVarKnownClusters$new()

Create a new ['ZINormalBlockVarKnownClusters'] object.

Usage
ZINormalBlockVarKnownClusters$new(
  data,
  C,
  sparsity = 0,
  control = NB_control()
)
Arguments
data

object of NormalBlockData class, with responses and design matrix

C

clustering matrix C_jk = 1 if species j belongs to cluster k

sparsity

to apply on variance matrix when calling GLASSO

control

structured list of more specific parameters, to generate with NB_control

Returns

A new ['ZINormalBlockVarKnownClusters'] object

ZINormalBlockVarKnownClusters$clone()

The objects of this class are cloneable with this method.

Usage
ZINormalBlockVarKnownClusters$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Zero-Inflated Normal-Block Model with Unknown Clustering

Description

R6 class for a zero-inflated normal-block model with a fixed number of clusters (but unknown clustering).

Super class

NormalBlockVarBase -> ZINormalBlockVarUnknownClusters

Public fields

fixed_tau

whether tau should be fixed at clustering_init during optimization, useful for stability selection

Active bindings

nb_param

number of parameters in the model

var_par

a list with variational parameters

model_par

a list with model parameters: B (covariates), dm1 (species variance), Omegaq (blocks precision matrix), kappa (zero-inflation probabilities)

entropy

Entropy of the conditional distribution

fitted

Y values predicted by the model, in Y's original units

who_am_I

a method to print what model is being fitted

Methods

Public methods

Inherited methods

ZINormalBlockVarUnknownClusters$new()

Create a new ['ZINormalBlockVarUnknownClusters'] object.

Usage
ZINormalBlockVarUnknownClusters$new(
  data,
  q,
  sparsity = 0,
  control = NB_control()
)
Arguments
data

object of NormalBlockVarData class, with responses and design matrix

q

required number of groups

sparsity

to apply on variance matrix when calling GLASSO

control

structured list of more specific parameters

Returns

A new ['ZINormalBlockVarUnknownClusters'] object

ZINormalBlockVarUnknownClusters$clone()

The objects of this class are cloneable with this method.

Usage
ZINormalBlockVarUnknownClusters$clone(deep = FALSE)
Arguments
deep

Whether to make a deep clone.

Breast cancer proteomics data (TCGA, RPPA)

Description

Reverse-phase protein array (RPPA) measurements of 163 proteins across 346 breast cancer tumor samples from The Cancer Genome Atlas (TCGA), along with the PAM50 molecular subtype of each sample. Used as the running example in Tous & Chiquet (2026).

Usage

brca_rppa

Format

A list with 3 elements:

expr

a 346 x 163 numeric matrix of protein expression levels (normalized log-ratios), samples in rows (named by TCGA sample identifier), proteins in columns.

covariates

a data frame with one row per sample, in the same order as the rows of 'expr': 'sampleId' (matches 'rownames(expr)') and 11 clinical variables, factors unless noted otherwise – 'CN_CLUSTER' (copy-number cluster, 5 levels), 'CONVERTED_STAGE' (tumor stage), 'ER_STATUS' (estrogen receptor status), 'FRACTION_GENOME_ALTERED' (numeric, in \[0, 1\]), 'HER2_STATUS', 'METASTASIS_CODED', 'METHYLATION_CLUSTER' (5 levels), 'MUTATION_COUNT' (numeric), 'PAM50_SUBTYPE' (5 levels: Basal-like, HER2-enriched, Luminal A, Luminal B, Normal-like), 'RPPA_CLUSTER' (7 levels) and 'AGE' (numeric, in years). A few of these have missing values for a handful of samples ('FRACTION_GENOME_ALTERED', 'HER2_STATUS', 'METASTASIS_CODED').

gene_annotation

a data frame with one row per protein, in the same order as the columns of 'expr': 'protein' (matches 'colnames(expr)'), 'entrezGeneId', 'hugoGeneSymbol', 'type' ('"protein-coding"' or '"phosphoprotein"') and 'go_bp_term' (a Gene Ontology Biological Process classification, one term per protein, or '"unknown"' if none could be found – see Details).

Details

'go_bp_term' is derived from 'org.Hs.eg.db', queried by 'entrezGeneId' for '"protein-coding"' rows. Phospho-site antibodies ('type == "phosphoprotein"', e.g. '"EGFR_PY1068"') carry a placeholder, negative 'entrezGeneId' (there is no separate gene for a specific phosphorylation site, only the underlying gene) and are instead queried by gene symbol (the part of 'protein' before the first '"_"'). Each protein typically has dozens of GO Biological Process terms; 'go_bp_term' keeps, for each protein, the one term shared by the most *other* proteins in this dataset – giving a handful of non-trivial groups rather than one near-singleton group per protein, more useful for comparing against an inferred clustering. See 'data-raw/brca_rppa_go_annotation.R' for the full extraction code.

Source

TCGA breast cancer cohort (Cancer Genome Atlas Network, 2012, doi:10.1038/nature11412), downloaded via cBioPortal ('brca_tcga_pub' study, https://www.cbioportal.org). 'go_bp_term' (see Details) was added afterwards from 'org.Hs.eg.db'.

Examples

expr <- brca_rppa$expr
X <- model.matrix(~ 0 + PAM50_SUBTYPE, data = brca_rppa$covariates)
nb_data <- NormalBlockData$new(expr, X)
table(brca_rppa$gene_annotation$go_bp_term)

Extract Model Coefficients

Description

Extract coefficients from a NormalBlockVarBase object.

Usage

## S3 method for class 'NormalBlockVarBase'
coef(object, ...)

Arguments

object

An object of class NormalBlockVarBase.

...

not used, only here for S3 compatibility

Value

A matrix of coefficients extracted from the NormalBlockVarBase model.

Extract Fitted Values

Description

Extract fitted values from a NormalBlockVarBase object.

Usage

## S3 method for class 'NormalBlockVarBase'
fitted(object, ...)

Arguments

object

An object of class NormalBlockVarBase.

...

not used, only here for S3 compatibility

Value

A matrix of fitted values extracted from the object.

Generate Normal Block Data

Description

A function to draw data from the normal block model (see details). The function returns both the generated data and the corresponding model parameters, in a list.

Usage

generate_normal_block_var_data(
  n = 100,
  p = 40,
  d = 1,
  q = 3,
  kappa = 0,
  omega_structure = "erdos-renyi",
  u_v = c(0.3, 0.1),
  SNR = 0.75,
  alpha = rep(1/q, q),
  range_X = c(0, 10),
  range_D = c(0.5, 1.5)
)

Arguments

n

number of individuals. Default to 100.

p

number of variables. Default to 40.

d

number of covariates. Default to 1.

q

number of groups. Default to 3.

kappa

vector (or scalar) of variable-wise probability of zero inflation. Default to 0.

omega_structure

the structure of the graph on which the precision matrix between groups is built. Can be a symmetric matrix with q rows/columns or a character picked in "erdos-renyi", "preferential_attachment", "community" in which case a graph is drawn with sensible generation parameters. See generate_precision_matrix for details.

u_v

two-size vector of positive numbers v and u controlling the generation of the precision matrix Omega: v scales the off-diagonal elements of the precision matrix (magnitude of partial correlations), and u is a positive number added to the diagonal elements to ensure positive-definiteness. The default value is c(0.3, 0.1).

SNR

Signal to noise ratio: magnitude of the regression parameters B will be adjusted so that tr(var(XB)) and tr(Sigma) match the desired SNR.

alpha

the q-size vector of group proportion. Default to rep(1/q, q)

range_X

A 2-size vector defining the range of the uniform distribution used to draw values in X, the regressor matrix. Default is c(0, 10)

range_D

A 2-size vector defining the range of the uniform distribution used to draw values in D, the diagonal matrix of variances of variables. Default is c(0.5, 1.5)

Value

A named list with the following element - Y a matrix of responses - X a regressor/design matrix - a list of model parameters, encompassing - B: matrix of regression coefficients - C: matrix of group membership - D: diagonal matrix of variance of the variables - Omega: precision matrix of the groups - Sigma: covariance matrix of the groups - kappa: vector of ZI inflation probabilities (one per variable)

Create a Normal-Block Model Object

Description

Creates the appropriate normal-block model (or collection of models) depending on the parametrization.

Usage

get_model(
  data,
  blocks,
  sparsity = 0,
  zero_inflation = FALSE,
  control = NB_control()
)

Arguments

data

contains the matrix of responses (Y) and the design matrix (X).

blocks

either an integer (number of blocks), a vector of integer (list of possible number of block) or a p * q matrix (for indicating block membership when its known)

sparsity

boolean to say whether the model should have a changing penalty OR float to run model with a single penalty value

zero_inflation

boolean to indicate if Y is zero-inflated and adjust fitted model as a consequence

control

a list-like structure for detailed control on parameters should be generated with NB_control() for collections of sparse models

Check if an Object is a Normal-Block Model

Description

Checks if a model is of class [NormalBlockVarBase()].

Usage

isNB(object)

Arguments

object

An R object.

Value

A boolean telling whether object inherits from the NormalBlockVarBase class.

Extract Log-Likelihood of a Normal-Block Model

Description

Returns the (variational) log-likelihood of a fitted normal-block model as a '"logLik"' object, compatible with [stats::AIC()] and [stats::BIC()].

Usage

## S3 method for class 'NormalBlockVarBase'
logLik(object, ...)

Arguments

object

An object of class NormalBlockVarBase.

...

not used, only here for S3 compatibility

Value

An object of class '"logLik"'. The numeric value is the log-likelihood or its variational lower bound (ELBO). Attributes 'df' and 'nobs' hold the number of parameters and observations.

Examples

ex_data <- generate_normal_block_var_data(n = 50, p = 20, d = 1, q = 3)
data <- NormalBlockData$new(ex_data$Y, ex_data$X)
model <- normal_block(data, blocks = 3, control = NB_control(verbose = FALSE))
logLik(model)

Extract Log-Likelihood of a Collection of Normal-Block Models

Description

Returns the log-likelihood of every model in a collection of normal-block models (see [NormalBlockVarCollectionClusters], [NormalBlockVarCollectionSparsity], [NormalBlockVarCollectionClustersSparsity]).

Usage

## S3 method for class 'NormalBlockVarCollection'
logLik(object, ...)

Arguments

object

An object inheriting from NormalBlockVarCollection.

...

not used, only here for S3 compatibility

Value

A numeric vector of log-likelihood values, one per model in the collection.

Examples

ex_data <- generate_normal_block_var_data(n = 50, p = 20, d = 1, q = 3)
data <- NormalBlockData$new(ex_data$Y, ex_data$X)
models <- normal_block(data, blocks = 2:5, control = NB_control(verbose = FALSE))
logLik(models)

Normal-block model

Description

Fit a normal-block model with a variational or heuristic algorithm

Usage

normal_block(
  data,
  blocks,
  sparsity = 0,
  zero_inflation = FALSE,
  control = NB_control()
)

Arguments

data

NormalBlockData object, contains the matrix of responses (Y, n x p) and the design matrix (X, n x d), must be created with NormalBlockData$new.

blocks

either a integer (number of blocks), a vector of integer (list of possible number of block) or a p * q matrix (for indicating block membership when its known)

sparsity

either TRUE to run the optimization for different sparsity penalty values OR float to run model with a single sparsity penalty value

zero_inflation

boolean to indicate if Y is zero-inflated and adjust fitted model as a consequence

control

a list-like structure for detailed control on parameters should be generated with NB_control().

Value

an R6 object with one of the model classes (or a collection of model objects).

Examples

## Normal Data
ex_data <- generate_normal_block_var_data(n=50, p=50, d=1, q=3)
data <- NormalBlockData$new(ex_data$Y, ex_data$X)
my_normal_block <- normal_block(data, blocks = 1:6)
my_normal_block$plot(c("deviance", "BIC", "ICL"))
Y_hat <- my_normal_block$get_best_model()$fitted
plot(data$Y, Y_hat, log = "xy"); abline(0,1)
## Normal Data with Zero Inflation
ex_data_zi <- generate_normal_block_var_data(n=50, p=50, d=1, q=3, kappa = rep(0.5,50))
zidata <- NormalBlockData$new(ex_data_zi$Y, ex_data_zi$X)
my_normal_block <- normal_block(zidata, blocks = 1:6, zero_inflation = TRUE)

French stream fish community data (ONEMA / OFB electrofishing surveys)

Description

Fish biomass per species, aggregated by sampling station, together with environmental covariates for each station. Derived from the fish community monitoring of stream sections across metropolitan France (1995-2018) by the French Office for Biodiversity (formerly ONEMA, "Office National de l'Eau et des Milieux Aquatiques"), using electrofishing.

Usage

onema

Format

A list with 2 elements:

biomass

a 399 x 46 numeric matrix of total biomass (grams) per species (3-letter species code, columns) and station (row names).

covariates

a data frame with one row per station, in the same order as the rows of 'biomass': 'station' (matches 'rownames(biomass)') and 14 environmental variables – 'slope', 'alt' (altitude), 'd_source' (distance to source), 'strahler' (Strahler stream order), 'width_river_mean'/'width_river_cv', 'avg_depth_station_mean'/ 'avg_depth_station_cv', 'DBO_med'/'DBO_cv' (biological oxygen demand), 'flow_med'/'flow_cv' and 'temperature_med'/'temperature_cv' ('_med' = median, '_cv' = coefficient of variation, over the monitoring period).

Details

The source data is organized around fishing *operations* ('opcod'): several electrofishing operations can be carried out at the same *station* over time. 'onema' aggregates this to one row per station: every operation is first mapped to its station ('fishing_protocol.csv'), then biomass (in grams) is summed over every operation recorded at a given station, separately for each species (i.e. 'biomass[s, ]' is the total biomass of each species ever caught at station 's', not a single operation's catch). Stations without environmental data are dropped. Species observed at fewer than 10 stations are also dropped: at that level of rarity, a zero-inflated fit's iterative initialization can fit a species' handful of non-zero observations exactly, driving its estimated noise precision to infinity (these species also carry essentially no information for clustering anyway).

Source

Danet, A., Mouchet, M., Bonnaffe, W., Thebault, E., Fontaine, C. (2021) "Species richness and food-web structure jointly drive community biomass and its temporal stability in fish communities", data set, Zenodo, doi:10.5281/zenodo.5095656.

Examples

Y <- log(1 + onema$biomass)
X <- model.matrix(~ 1, data = onema$covariates)
nb_data <- NormalBlockData$new(Y, X)

out <- normal_block(nb_data, 2:15, control = NB_control(clustering_init = "ward2"))

Plot a Normal-Block Model

Description

Plots the evolution of the objective (log-likelihood or ELBO) across the (V)EM iterations of the last call to 'optimize()', see '$plot_loglik()'.

Usage

## S3 method for class 'NormalBlockVarBase'
plot(x, ...)

Arguments

x

An object of class NormalBlockVarBase.

...

not used, only here for S3 compatibility

Value

Invisibly returns the [ggplot2::ggplot] object; called for its side effect of plotting.

Examples

ex_data <- generate_normal_block_var_data(n = 50, p = 20, d = 1, q = 3)
data <- NormalBlockData$new(ex_data$Y, ex_data$X)
model <- normal_block(data, blocks = 3, control = NB_control(verbose = FALSE))
plot(model)

Predict Method for Normal-Block Models

Description

Predicts observations Y for new covariates X.

Usage

## S3 method for class 'NormalBlockVarBase'
predict(object, new_X, ...)

Arguments

object

An object of class NormalBlockVarBase.

new_X

New set of covariates.

...

not used, only here for S3 compatibility

Value

A n*p prediction matrix for new observations

Print a Normal-Block Model

Description

Print a short summary of a fitted normal-block model: model type, goodness-of-fit criteria, and the useful fields/methods to explore it further.

Usage

## S3 method for class 'NormalBlockVarBase'
print(x, ...)

Arguments

x

An object of class NormalBlockVarBase.

...

not used, only here for S3 compatibility

Value

Invisibly returns 'x'.

Examples

ex_data <- generate_normal_block_var_data(n = 50, p = 20, d = 1, q = 3)
data <- NormalBlockData$new(ex_data$Y, ex_data$X)
model <- normal_block(data, blocks = 3, control = NB_control(verbose = FALSE))
print(model)

Print a Collection of Normal-Block Models

Description

Print a short summary of a collection of normal-block models: model type and the range of q/sparsity explored. See [summary.NormalBlockVarCollection()] for the full criteria table.

Usage

## S3 method for class 'NormalBlockVarCollection'
print(x, ...)

Arguments

x

An object inheriting from NormalBlockVarCollection.

...

not used, only here for S3 compatibility

Value

Invisibly returns 'x'.

Examples

ex_data <- generate_normal_block_var_data(n = 50, p = 20, d = 1, q = 3)
data <- NormalBlockData$new(ex_data$Y, ex_data$X)
models <- normal_block(data, blocks = 2:5, control = NB_control(verbose = FALSE))
print(models)

Print a Normal-Block Model Summary

Description

Print method for objects returned by [summary.NormalBlockVarBase()].

Usage

## S3 method for class 'summary.NormalBlockVarBase'
print(x, ...)

Arguments

x

An object of class 'summary.NormalBlockVarBase'.

...

not used, only here for S3 compatibility

Value

Invisibly returns 'x'.

Print a Collection Summary

Description

Print method for objects returned by [summary.NormalBlockVarCollection()].

Usage

## S3 method for class 'summary.NormalBlockVarCollection'
print(x, ...)

Arguments

x

An object of class 'summary.NormalBlockVarCollection'.

...

not used, only here for S3 compatibility

Value

Invisibly returns 'x'.

Extract the Latent-Block Covariance Matrix

Description

Extract the covariance matrix between latent blocks (the inverse of 'Omegaq') from a [NormalBlockVarBase()] object.

Usage

## S3 method for class 'NormalBlockVarBase'
sigma(object, ...)

Arguments

object

An object of class NormalBlockVarBase.

...

not used, only here for S3 compatibility

Value

The q x q covariance matrix between latent blocks.

Summarize a Normal-Block Model

Description

Summarizes a fitted normal-block model: model type, goodness-of-fit criteria, cluster sizes, and the density of the inferred network between blocks.

Usage

## S3 method for class 'NormalBlockVarBase'
summary(object, ...)

Arguments

object

An object of class NormalBlockVarBase.

...

not used, only here for S3 compatibility

Value

An object of class 'summary.NormalBlockVarBase' (a list with the model's 'who_am_I', 'criteria', 'cluster_sizes' and network 'density'), printed with a dedicated [print.summary.NormalBlockVarBase()] method.

Examples

ex_data <- generate_normal_block_var_data(n = 50, p = 20, d = 1, q = 3)
data <- NormalBlockData$new(ex_data$Y, ex_data$X)
model <- normal_block(data, blocks = 3, control = NB_control(verbose = FALSE))
summary(model)

Summarize a Collection of Normal-Block Models

Description

Summarizes a collection of normal-block models: model type, the full criteria table, and the range of q/sparsity explored.

Usage

## S3 method for class 'NormalBlockVarCollection'
summary(object, ...)

Arguments

object

An object inheriting from NormalBlockVarCollection.

...

not used, only here for S3 compatibility

Value

An object of class 'summary.NormalBlockVarCollection' (a list with the collection's 'who_am_I', 'criteria', 'q_range' and 'sparsity_range'), printed with a dedicated [print.summary.NormalBlockVarCollection()] method.

Examples

ex_data <- generate_normal_block_var_data(n = 50, p = 20, d = 1, q = 3)
data <- NormalBlockData$new(ex_data$Y, ex_data$X)
models <- normal_block(data, blocks = 2:5, control = NB_control(verbose = FALSE))
summary(models)

University webpages text data (CMU "4 Universities" / WebKB)

Description

Term-frequency data derived from the student personal webpages of the CMU "World Wide Knowledge Base" (WebKB) project's "4 Universities" dataset (Cornell, Texas, Washington, Wisconsin computer science departments, collected January 1997). Used as one of the running examples in Tous & Chiquet (2026).

Usage

university

Format

A list with 3 elements:

frequencies

a 504 x 1867 numeric matrix of term frequencies (row sums to 1): for each document (row, named by its original file path) and term (column), the fraction of that document's (post-preprocessing) word count made up of that term.

entropies

a named numeric vector of length 1867 (one value per column of 'frequencies', same names/order), the normalized Shannon entropy of each term's distribution across documents.

terms

a character vector of the 100 column names of 'frequencies' with the highest 'entropies', ordered by decreasing entropy.

Details

Each of the 504 pages (rows) is a document; each of the 1867 columns is a term retained after lower-casing, stripping URLs/HTML tags/control characters/punctuation/numbers, removing English stopwords and dropping terms occurring in fewer than 2 documents. 'entropies' ranks every term by how evenly it is spread across documents (Shannon entropy of each term's normalized document distribution, in \[0, 1\], 1 = perfectly uniform); 'terms' is the 100 highest-entropy terms, i.e. the terms that are the most informative for distinguishing documents from one another rather than just reflecting a few documents' idiosyncratic vocabulary – the transformation used by Tan et al. (2015).

Source

CMU Text Learning Group, "World Wide Knowledge Base (Web->KB) project", http://www.cs.cmu.edu/~webkb/; "4 Universities" subset, https://www.cs.cmu.edu/afs/cs/project/theo-20/www/data/. Tan, P.-N., Steinbach, M., Kumar, V. (2015) "Introduction to Data Mining" (transform used to derive 'entropies'/'terms').

Examples

Y <- log(1 + university$frequencies[, university$terms])
nb_data <- NormalBlockData$new(Y, X = matrix(1, nrow(Y), 1))

out <- normal_block(nb_data, 2:15)