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


