Sparse structure estimation for multivariate count data with PLN-network

Requirements

The packages required for the analysis are PLNmodels plus some others for data manipulation and representation:

library(PLNmodels)
library(ggplot2)

Data set

We illustrate our point with the trichoptera data set, a full description of which can be found in the corresponding vignette. Data preparation is also detailed in the specific vignette.

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)

The trichoptera data frame stores a matrix of counts (trichoptera$Abundance), a matrix of offsets (trichoptera$Offset) and some vectors of covariates (trichoptera$Wind, trichoptera$Temperature, etc.)

Mathematical background

The network model for multivariate count data that we introduce in Chiquet et al. (2019) is a variant of the Poisson Lognormal model of Aitchison and Ho (1989), see the PLN vignette as a reminder. Compare to the standard PLN model we add a sparsity constraint on the inverse covariance matrix \({\boldsymbol\Sigma}^{-1}\triangleq \boldsymbol\Omega\) by means of the \(\ell_1\)-norm, such that \(\|\boldsymbol\Omega\|_1 < c\). PLN-network is the equivalent of the sparse multivariate Gaussian model (Banerjee et al. 2008) in the PLN framework. It relates some \(p\)-dimensional observation vectors \(\mathbf{Y}_i\) to some \(p\)-dimensional vectors of Gaussian latent variables \(\mathbf{Z}_i\) as follows \[\begin{equation} \begin{array}{rcl} \text{latent space } & \mathbf{Z}_i \sim \mathcal{N}\left({\boldsymbol\mu},\boldsymbol\Omega^{-1}\right) & \|\boldsymbol\Omega\|_1 < c \\ \text{observation space } & Y_{ij} | Z_{ij} \quad \text{indep.} & Y_{ij} | Z_{ij} \sim \mathcal{P}\left(\exp\{Z_{ij}\}\right) \end{array} \end{equation}\]

The parameter \({\boldsymbol\mu}\) corresponds to the main effects and the latent covariance matrix \(\boldsymbol\Sigma\) describes the underlying structure of dependence between the \(p\) variables.

The \(\ell_1\)-penalty on \(\boldsymbol\Omega\) induces sparsity and selection of important direct relationships between entities. Hence, the support of \(\boldsymbol\Omega\) correspond to a network of underlying interactions. The sparsity level (\(c\) in the above mathematical model), which corresponds to the number of edges in the network, is controlled by a penalty parameter in the optimization process sometimes referred to as \(\lambda\). All mathematical details can be found in Chiquet et al. (2019).

Adjusting a collection of network - a.k.a. a regularization path

PLNnetwork finds an hopefully appropriate set of penalties on its own. This set can be controlled by the user, but use it with care and check details in ?PLNnetwork. By default the penalty applies on the scale of the correlations between the latent variables, so that it does not depend on their variances (the correlation-based estimator of Rothman et al. 2008): the penalties lie between 0 and 1, and the path starts from the empty network. The penalty on the covariance scale of versions up to 1.3.2 is obtained with PLNnetwork_param(penalty_scale = "covariance"). The collection of models is fitted as follows:

network_models <- PLNnetwork(Abundance ~ 1 + offset(log(Offset)), data = trichoptera)

## 
##  Initialization...
##  Adjusting 30 PLN with sparse inverse covariance estimation
##  Joint optimization alternating gradient descent and graphical-lasso
##  sparsifying penalty = 0.3985143     sparsifying penalty = 0.368096  sparsifying penalty = 0.3399996     sparsifying penalty = 0.3140477     sparsifying penalty = 0.2900767     sparsifying penalty = 0.2679354     sparsifying penalty = 0.2474841     sparsifying penalty = 0.2285939     sparsifying penalty = 0.2111455     sparsifying penalty = 0.1950289     sparsifying penalty = 0.1801425     sparsifying penalty = 0.1663924     sparsifying penalty = 0.1536918     sparsifying penalty = 0.1419607     sparsifying penalty = 0.1311249     sparsifying penalty = 0.1211163     sparsifying penalty = 0.1118716     sparsifying penalty = 0.1033325     sparsifying penalty = 0.09544523    sparsifying penalty = 0.08815998    sparsifying penalty = 0.0814308     sparsifying penalty = 0.07521526    sparsifying penalty = 0.06947414    sparsifying penalty = 0.06417124    sparsifying penalty = 0.0592731     sparsifying penalty = 0.05474884    sparsifying penalty = 0.05056991    sparsifying penalty = 0.04670995    sparsifying penalty = 0.04314462    sparsifying penalty = 0.03985143 
##  Post-treatments
##  DONE!

Note the use of the formula object to specify the model, similar to the one used in the function PLN.

Structure of PLNnetworkfamily

The network_models variable is an R6 object with class PLNnetworkfamily, which comes with a couple of methods. The most basic is the show/print method, which sends a very basic summary of the estimation process:

network_models

## --------------------------------------------------------
## COLLECTION OF 30 POISSON LOGNORMAL MODELS
## --------------------------------------------------------
##  Task: Network Inference 
## ========================================================
##  - 30 penalties considered: from 0.03985143 to 0.3985143 
##  - Best model (greater BIC): lambda = 0.399 
##  - Best model (greater EBIC): lambda = 0.399

One can also easily access the successive values of the criteria in the collection

network_models$criteria %>% head() %>% knitr::kable()

A diagnostic of the optimization process is available via the convergence field:

network_models$convergence %>% head() %>% knitr::kable()

An nicer view of this output comes with the option “diagnostic” in the plot method:

plot(network_models, "diagnostic")

Exploring the path of networks

By default, the plot method of PLNnetworkfamily displays evolution of the criteria mentioned above, and is a good starting point for model selection:

plot(network_models)

Note that we use the original definition of the BIC/ICL criterion (\(\texttt{loglik} - \frac{1}{2}\texttt{pen}\)), which is on the same scale as the log-likelihood. A popular alternative consists in using \(-2\texttt{loglik} + \texttt{pen}\) instead. You can do so by specifying reverse = TRUE:

plot(network_models, reverse = TRUE)

In this case, the variational lower bound of the log-likelihood is hopefully strictly increasing (or rather decreasing if using reverse = TRUE) with a lower level of penalty (meaning more edges in the network). The same holds true for the penalized counterpart of the variational surrogate. Generally, smoothness of these criteria is a good sanity check of optimization process. BIC and its extended-version high-dimensional version EBIC are classically used for selecting the correct amount of penalization with sparse estimator like the one used by PLN-network. However, we will consider later a more robust albeit more computationally intensive strategy to chose the appropriate number of edges in the network.

To pursue the analysis, we can represent the coefficient path (i.e., value of the edges in the network according to the penalty level) to see if some edges clearly come off. An alternative and more intuitive view consists in plotting the values of the partial correlations along the path, which can be obtained with the options corr = TRUE. To this end, we provide the S3 function coefficient_path:

coefficient_path(network_models, corr = TRUE) %>%
  ggplot(aes(x = Penalty, y = Coeff, group = Edge, colour = Edge)) +
    geom_line(show.legend = FALSE) +  coord_transform(x="log10") + theme_bw()

Model selection issue: choosing a network

To select a network with a specific level of penalty, one uses the getModel(lambda) S3 method. We can also extract the best model according to the BIC or EBIC with the method getBestModel().

model_pen <- getModel(network_models, network_models$penalties[20]) # give some sparsity
model_BIC <- getBestModel(network_models, "BIC")   # if no criteria is specified, the best BIC is used

An alternative strategy is to use StARS (Liu et al. 2010), which performs resampling to evaluate the robustness of the network along the path of solutions in a similar fashion as the stability selection approach of Meinshausen and Bühlmann (2010), but in a network inference context.

Resampling can be computationally demanding but is easily parallelized: the function stability_selection relies on parallel::mclapply to perform parallel computing. Set the number of workers with the mc.cores option (forking-based, so only effective on Unix-like systems; ignored on Windows):

options(mc.cores = 2)

We first invoke stability_selection explicitly for pedagogical purpose. In this case, we need to build our sub-samples manually:

n <- nrow(trichoptera)
subs <- replicate(10, sample.int(n, size = n/2), simplify = FALSE)
stability_selection(network_models, subsamples = subs)

## 
## Stability Selection for PLNnetwork: 
## subsampling: ++++++++++

Requesting ‘StARS’ in gestBestmodel automatically invokes stability_selection with 20 sub-samples, if it has not yet been run.

model_StARS <- getBestModel(network_models, "StARS")

When “StARS” is requested for the first time, getBestModel automatically calls the method stability_selection with the default parameters. After the first call, the stability path is available from the plot function:

plot(network_models, "stability")

When you are done, do not forget to get back to the default (sequential) behavior.

options(mc.cores = 1)

Structure of a PLNnetworkfit

The variables model_BIC, model_StARS and model_pen are other R6Class objects with class PLNnetworkfit. They all inherits from the class PLNfit and thus own all its methods, with a couple of specific one, mostly for network visualization purposes. Most fields and methods are recalled when such an object is printed:

model_StARS

## Poisson Lognormal with sparse inverse covariance (penalty = 0.29)
## ==================================================================
##  nb_param    loglik       BIC       AIC      ICL n_edges      EBIC pen_loglik
##        37 -1106.539 -1178.538 -1143.539 -2155.23       3 -1187.038   -1109.88
##  density
##    0.022
## ==================================================================
## * Useful fields
##     $model_par, $latent, $latent_pos, $var_par, $optim_par
##     $loglik, $BIC, $ICL, $loglik_vec, $nb_param, $criteria
## * Useful S3 methods
##     print(), coef(), sigma(), vcov(), fitted()
##     predict(), predict_cond(), standard_error()
## * Additional fields for sparse network
##     $EBIC, $density, $penalty 
## * Additional S3 methods for network
##     plot.PLNnetworkfit()

The plot method provides a quick representation of the inferred network, with various options (either as a matrix, a graph, and always send back the plotted object invisibly if users needs to perform additional analyses).

my_graph <- plot(model_StARS, plot = FALSE)
my_graph

## IGRAPH fd54d51 UNW- 17 3 -- 
## + attr: name (v/c), label (v/c), label.cex (v/n), size (v/n),
## | label.color (v/c), weight (e/n), width (e/n), color (e/c)
## + edges from fd54d51 (vertex names):
## [1] Psy--Han Psy--Sta Hfo--Hsp

plot(model_StARS)

plot(model_StARS, type = "support", output = "corrplot")

We can finally check that the fitted value of the counts – even with sparse regularization of the covariance matrix – are close to the observed ones:

data.frame(
  fitted   = as.vector(fitted(model_StARS)),
  observed = as.vector(trichoptera$Abundance)
) %>%
  ggplot(aes(x = observed, y = fitted)) +
    geom_point(size = .5, alpha =.25 ) +
    scale_x_log10(limits = c(1,1000)) +
    scale_y_log10(limits = c(1,1000)) +
    theme_bw() + annotation_logticks()

fitted value vs. observation