--- title: "Analyzing multivariate Gaussian data with the Normal-Block model - First steps" author: "Jeanne Tous, Julien Chiquet" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 4 bibliography: references.bib vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteIndexEntry{Analyzing multivariate Gaussian data with the Normal-Block model - First steps} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # Preliminaries This vignette illustrates the use of the `normal_block()` function and the methods accompanying its R6 classes. From a statistical point of view, the `normal_block()` function fits a multivariate Normal-Block model (a Gaussian graphical model with a latent clustering structure) to a table of observations, possibly after correcting for the effect of covariates. Depending on the arguments given, the function uses a clustering supplied by the user or infers one, includes zero-inflation or not, infers a sparse network or not, and so on. Parameter inference can either be done with an integrated variational Expectation-Maximization approach (recommended) or with a faster heuristic approach. Theoretical explanations for this model can be found in @tous2026. See `inst/normal_block_models.qmd` in the package sources for the full estimation details (criteria, E/M updates). ### Requirements ```{r, set-up} library(normalblockr) library(pheatmap) library(paletteer) ``` # Data simulation We illustrate the analysis with two simulated datasets. Fix the seed to make the results reproducible. ```{r, set-seed} set.seed(1) ``` Simulate data with `generate_normal_block_data()`. ## Fix simulation parameters Several parameters need to be defined to simulate the data: ```{r, fix-simulations-parameters} n = 100 # Number of samples (number of rows in the matrix of observations) p = 40 # Number of entities observed (number of columns in the matrix of observations) d = 2 # Number of covariates q = 3 # Number of clusters kappa = 0 # Mean zero-inflation probability (can also be a vector to define one ZI-probability for each variable). kappa = 0 means that there is no zero-inflation. omega_structure = "erdos-renyi" # Network structure. u_v = c(0.3, 0.1) # Parameters to generate an association matrix from a graph, details given in the bibliography. SNR = 0.75 # Signal to Noise Ratio, defines the relative weight of the covariates and the variance alpha = rep(1/q, q) # Vector giving probabilities of belonging to each cluster range_X = c(0, 10) # Min and max values for the covariates range_D = c(0.5, 1.5) # Min and max values for the individual entities variances ``` ## Simulation 1 (no zero-inflation) `generate_normal_block_var_data()` generates data under the Normal-Block model. It returns a list containing the simulated covariates $X$ and observations $Y$, and the simulation parameters (including the clustering $C$). ```{r, simulation1} my_nb_data <- generate_normal_block_var_data(n, p, d, q, kappa, omega_structure, u_v, SNR, alpha, range_X, range_D) ``` ```{r, simulation1-visualization} pheatmap::pheatmap(my_nb_data$Y, color =paletteer::paletteer_c("ggthemes::Orange-Gold", n = 100), cluster_rows = FALSE, cluster_cols = FALSE, show_rownames = FALSE) ``` ## Simulation 2 (with zero-inflation) To generate zero-inflated data, we set-up kappa as a vector of zero-inflation probabilities. One needs to ensure that these values are between 0 and 1. ```{r, simulation2} kappa_zi <- rnorm(p, mean = 0.7, sd = 0.05) kappa_zi <- unlist(lapply(kappa_zi, f <- function(x) return(max(0, min(x, 0.9))))) my_nb_data_zi <- generate_normal_block_var_data(n = n, p = p, d = d, q = 5, kappa = kappa_zi, omega_structure = omega_structure, u_v = u_v, SNR = SNR, alpha = alpha, range_X = range_X, range_D = range_D) ``` For a better visualization of the zero-inflated data, we set-up the 0 to be shown in white. ```{r, simulation2-visualization} min_val <- min(my_nb_data_zi$Y) ; max_val <- max(my_nb_data_zi$Y) orange_gold_pal <- paletteer::paletteer_c("ggthemes::Orange-Gold", n = 100) n_breaks <- 100 zero_pos <- round((0 - min_val) / (max_val - min_val) * n_breaks) + 1 custom_pal <- c(orange_gold_pal[1:(zero_pos - 1)], "white", orange_gold_pal[zero_pos:n_breaks]) pheatmap::pheatmap(my_nb_data_zi$Y, color = custom_pal, cluster_rows = FALSE, cluster_cols = FALSE, show_rownames = FALSE) ``` # Prepare the data for a Normal-Block analysis A specific data object of class `NormalBlockData` needs to be created to analyse the data with `normalblockr`. One can use $Y$ and $X$ alone to create the `NormalBlockData` object. ```{r, NormalBlockData} my_data <- NormalBlockData$new(my_nb_data$Y, my_nb_data$X) my_data_zi <- NormalBlockData$new(my_nb_data_zi$Y, my_nb_data_zi$X) ``` Alternatively, if one only wants to include some of the covariates contained in $X$, a formula can be used to specify which one: ```{r, NormalBlockData-alt} colnames(my_data$X) <- c("X1", "X2") my_data_alt <- NormalBlockData$new(my_data$Y, my_data$X, formula = ~ 0 + X1) ``` # Run a Normal-Block analysis All Normal-Block analyses are run with the `normal_block()` function, called with different arguments depending on whether the clustering (or the number of clusters) is known, and on the requested level of sparsity, among other factors. See `?normal_block` for full details. By default, the parameters are inferred using a variational Expectation-Maximization approach, with up to 500 iterations. Finer control of the optimization is possible through the `control` argument of `normal_block()`, a list generated by `NB_control()`. ## With non-zero-inflated data ### Fixed clustering When the variables' clustering is known, it can be given directly as an input to `normal_block()`. ```{r, simulation1-NB1} my_NB <- normal_block(data = my_data, blocks = my_nb_data$parameters$C) ``` Convergence of the model can be checked with `plot()`. ```{r, simulation1-NB1-plot} plot(my_NB) ``` A summary of the results is accessible via `print()`. ```{r, simulation1-NB1-print} print(my_NB) ``` The inter-cluster association network can be visualized with the `plot_network()` function. ```{r, simulation1-NB1-plot_network} my_NB$plot_network() ``` ### Fixed number of clusters When the variables' clustering is unknown, the number of clusters can simply be fixed via the `blocks` argument. ```{r, simulation1-NB2} my_NB <- normal_block(data = my_data, blocks = 3) print(my_NB) plot(my_NB) ``` The clustering inferred by `normal_block()` can be compared with the true clustering, when it is known: ```{r, simulation1-ARI} aricode::ARI(my_NB$clustering, apply(my_nb_data$parameters$C, 1, which.max)) ``` ### Unknown number of clusters When the number of clusters is unknown, `normal_block()` can be given a range of candidate values instead, returning a collection of Normal-Block models, one per number of clusters. ```{r, simulation1-NB3} my_NB_unknown <- normal_block(data = my_data, blocks = 2:5) ``` `my_NB_unknown` is a collection of Normal-Block models. A specific one can be selected either by its number of clusters or as the best model for a given criterion (BIC, deviance, EBIC or ICL). The value of each criterion, for every model in the collection, can be visualized with `plot()`. ```{r, simulation1-NB3-plot} plot(my_NB_unknown) ``` Model selection is done with `get_model()` when a specific number of clusters is required, and with `get_best_model()` to select the best model for a given criterion. ```{r, simulation1-NB3-model-selection} myNB_3 <- my_NB_unknown$get_model(3) myNB_BIC <- my_NB_unknown$get_best_model("BIC") ``` ### Playing with sparsity By default, no penalty is applied to the association network. A penalty can be added via the `sparsity` argument of `normal_block()`, with any of the parametrizations seen above (fixed clustering, fixed number of clusters, or unknown number of clusters). The example below uses a fixed clustering. To use a fixed $\ell_1$ penalty [@tous2026] on the network, pass that value to the `sparsity` argument: ```{r, simulation1-NB-fixed-sparsity} my_NB_sparse_low <- normal_block(data = my_data, blocks = my_nb_data$parameters$C, sparsity = 0.1) my_NB_sparse_high <- normal_block(data = my_data, blocks = my_nb_data$parameters$C, sparsity = 10) ``` The larger the penalty, the sparser the network. ```{r, simulation1-NB-low-sparsity-plot} my_NB_sparse_low$plot_network() ``` ```{r, simulation1-NB-high-sparsity-plot} my_NB_sparse_high$plot_network() ``` It is usually hard to know a priori which sparsity penalty is best. `normal_block()` can instead explore a range of sparsity levels by simply setting `sparsity = TRUE`, which returns a collection of Normal-Block models, one per sparsity penalty. ```{r, simulation1-NB-changing-sparsity} my_NB_sparse <- normal_block(data = my_data, blocks = my_nb_data$parameters$C, sparsity = TRUE) ``` The different criteria can then be plotted as a function of the penalty, using `plot()`. ```{r, simulation1-NB-changing-sparsity-plot} plot(my_NB_sparse) ``` Penalty selection can be done similarly to the selection of the number of clusters. ```{r, simulation1-NB-changing-sparsity-model-selection} myNB_sparse_0.1 <- my_NB_sparse$get_model(0.1) myNB_sparse_BIC <- my_NB_sparse$get_best_model("BIC") ``` Both the sparsity penalty and the number of clusters can also be left to vary jointly. ```{r, simulation1-NB-sparse_unknown} my_NB_sparse_unknown <- normal_block(data = my_data, blocks = 2:6, sparsity = TRUE) ``` The result is a collection of Normal-Block models with different numbers of clusters and different penalties. ```{r, simulation1-NB-sparse_unknown-who-am-I} my_NB_sparse_unknown$who_am_I ``` `plot()` allows each criterion to be analysed as a function of both the number of clusters and the penalty. By default, the penalties tested are computed by `normal_block()` separately for each number of clusters, and can differ from one to another -- which is what the blanks in the plot come from. ```{r, simulation1-NB-sparse_unknown-plot} plot(my_NB_sparse_unknown, "BIC") ``` Model selection is also done using `get_model()` and `get_best_model()`. With `get_model()`, one can fix only the number of clusters, getting back a collection of models with different sparsity levels, or fix both the number of clusters and the sparsity level. ```{r, simulation1-NB-sparse_unknown-selection} my_NB_sparse_3 <- my_NB_sparse_unknown$get_model(3) my_NB_sparse_3_0.1 <- my_NB_sparse_unknown$get_model(3, 0.1) ``` ## With zero-inflated data The process is similar with zero-inflated data, but `zero_inflation = TRUE` must be passed to `normal_block()`. The example below uses a fixed number of blocks and no penalty on the network. ```{r, simulation2-NB} my_NB_zi <- normal_block(data = my_data_zi, blocks = 4, zero_inflation = TRUE) ``` When the data is zero-inflated, the inference process may take longer. ```{r, simulation2-NB-plot} plot(my_NB_zi) ``` The clustering inference may also be harder. ```{r, simulation2-ARI} aricode::ARI(my_NB_zi$clustering, apply(my_nb_data_zi$parameters$C, 1, which.max)) ``` # References