Lasso Method for RCON(V, E) Models

Description

RCON(V, E) models (Hojsgaard, et al., 2008) are a kind of restriction of the Gaussian Graphical Models defined by a set of equality constraints on the entries of the concentration matrix. sglasso package implements the structured graphical lasso (sglasso) estimator proposed in Abbruzzo et al. (2014) for the weighted l1-penalized RCON(V, E) model. Two cyclic coordinate algorithms are implemented to compute the sglasso estimator, i.e. a cyclic coordinate minimization (CCM) and a cyclic coordinate descent (CCD) algorithm.

Details

Author(s)

Luigi Augugliaro

Maintainer:
Luigi Augugliaro <luigi.augugliaro@unipa.it>

References

Abbruzzo, A., Augugliaro, L., Mineo, A. M. and Wit, E. C. (2014) Cyclic coordinate for penalized Gaussian Graphical Models with symmetry restrictions. In Proceeding of COMPSTAT 2014 - 21th International Conference on Computational Statistics, Geneva, August 19-24, 2014.

Hojsgaard, S. and Lauritzen, S. L. (2008) Graphical gaussian models with edge and vertex symmetries. J. Roy. Statist. Soc. Ser. B., Vol. 70(5), 1005–1027.

Extract Sparse Structured Precision Matrices

Description

Function Kh computes the sequence of sparse structured precision matrices estimated by sglasso function.

Usage

Kh(object, rho)

Arguments

Value

Kh returns a named list containing the sequence of estimated sparse structured precision matrices.

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

See Also

sglasso function.

Examples

N <- 100
p <- 5
X <- matrix(rnorm(N * p), N, p)
S <- crossprod(X) / N
mask <- outer(1:p, 1:p, function(i, j) 0.5^abs(i - j))
out.sglasso_path <- sglasso(S, mask, nrho = 5, tol = 1.0e-13)
out.sglasso_path
Kh(out.sglasso_path)
rho <- out.sglasso_path$rho[3]
out.sglasso_single <- sglasso(S, mask, nrho = 1, min_rho = rho, 
   tol = 1.0e-13, algorithm = "ccm")
Kh(out.sglasso_single)

L1-penalized Factorial Graphical Lasso Model

Description

Fit the weight l1-penlized factorial dynamic Gaussian Graphical Model.

Usage

fglasso(S, model, tp, p, ...)

Arguments

Details

The factorial dynamic Gaussian Graphical Model (Abbruzzo et al., 2015) is a special kind of RCON(V, E) model (Hojsgaard, et al.,2008) proposed to study dynamic networks. Let X_t = (X_{it},\ldots,X_{it})' be a p-dimensional random variable at time t. Assuming that X = (X'_1,\ldots,X'_T) follows a multivariate normal distribution, the concentration matrix K has the following block structure

K = \left( \begin{array}{cccc} K_{1,1} & K_{1,2} & \ldots & K_{1,T}\\ K_{2,1} & K_{2,2} & \ldots & K_{2,T}\\ \vdots & \vdots & \ddots & \vdots\\ K_{T,1} & K_{T,2} & \ldots & K_{T,T} \end{array}\right),

where K_{t,t} give information about the conditinal independence structure among the p random variables at time t, and K_{t,t + h} give information about the conditional independence structure between X_t and X_{t + h}. An interpretation of the elements of the submatrices K_{t,t + h} brings to the notion of natural structure, i.e.,

K_{t,t + h} = \left( \begin{array}{cccc} k_{1,1}^{t,t+h} & 0 & \ldots & 0\\ 0 & k_{2,2}^{t,t+h} & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & k_{p,p}^{t,t+h} \end{array}\right) + \left( \begin{array}{cccc} 0 & k_{1,2}^{t,t+h} & \ldots & k_{1,p}^{t,t+h}\\ k_{2,1}^{t,t+h} & 0 & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ k_{p,1}^{t,t+h} & k_{p,2}^{t,t+h} & \ldots & 0 \end{array}\right).

The entries of the first matrix are called self-self conditinal dependences at temporal lag h and represents the (negative) self-similarity of a given random variable across different time points. The entries of the second matrix are the conditional dependence among the p random variables. To make the interpretation of the results more relevant and, at the same time, reduce the number of parameters, the authors propose the following equality constraints:

Argument model is used to specify the restrinctions previously describted. This argument can be a named list or a matrix with dimension nlag\times 2, where nlag\le\code{tp}. To gain more insight, suppose that we want to model only the sub-matrices K_{t,t} and K_{t,t+1}, i.e., the sub-matrices corresponding to the temporal lag zero and one. A possible R code is

model.mat <- matrix("", nrow = 2, ncol = 2)
rownames(model.mat) <- c("lag0", "lag1")
colnames(model.mat) <- c("s", "n")
model.mat[1, ] <- c("c", "ut")
model.mat[2, ] <- c("t", ".")

In this example we are modelling the diagonal elements of the sub-matrices K_{t,t} with the constant effect while the off-diagonal elements are modelled by the interaction effect. In the same way, the diagonal elements of the sub-matrices K_{t,t+1} are modelled by the time effect while the remaning elements are equal to zero. The fglasso function passes the matrix model.mat to the internal function fglasso_model2mask, i.e.,

mask <- fglasso_model2mask(model.mat, tp = 3, p = 3)

which returns the mask used in sglasso to fit the specified factorial dynamic Gaussian Graphical model. The same model can be specified by the following named list

model.list <- list(lag0 = c(s = "c", n = "ut"), lag1 = c(s = "t", n = "."))

See the example below for more details.

Value

fglasso returns an obejct with S3 class "sglasso". See the corresponding manual for more details.

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

References

Wit, E. C. and Abbruzzo, A.(2015) Dynamic factorial graphical models for dynamic networks. Network Science, Vol. 3(1), 37– 57

Abbruzzo, A., Augugliaro, L., Mineo, A.M. and Wit, E.C. (2014) Cyclic coordinate for penalized Gaussian Graphical Models with symmetry restrictions. In Proceeding of COMPSTAT 2014 - 21th International Conference on Computational Statistics. Geneva, August 19-24, 2014.

Hojsgaard, S. and Lauritzen, S.L. (2008) Graphical gaussian models with edge and vertex symmetries. J. Roy. Statist. Soc. Ser. B., Vol. 70(5), 1005–1027.

See Also

sglasso function.

Examples

#######################
# fglasso solution path
#######################
N <- 50
tp <- 3
p <- 3
X <- matrix(rnorm(N * p * tp), N, tp * p)
S <- crossprod(X) / N
model <- list(lag0 = c(s = "c", n = "ut"), lag1 = c(s = "t", n = "."))
out.fglasso <- fglasso(S = S, model = model, tp = tp, p = p)
out.fglasso

Plotting Sparse Graph

Description

gplot is a generic function for plotting sparse graphs.

Usage

gplot(object, ...)

Arguments

Details

gplot is a generic function used to plot a graph estimated by sglasso or fglasso. See the method function gplot.sglasso or gplot.fglasso for more details about the specific arguments.

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

See Also

gplot.sglasso and gplot.fglasso method functions.

Plotting Sparse Factorial Dynamic Gaussian Graphical Model

Description

gplot.fglasso shows the sequence of graphs estimated by fglasso.

Usage

## S3 method for class 'fglasso'
gplot(object, rhoid, tp = c(1, 2), sub.tp1, sub.tp2, cex.sub = 1, 
    k = 1.5, layout = layout.circle, ...)

Arguments

Details

For a given value of the tuning parameter, specified by the argument rhoid, gplot.fglasso shows the graphs estimated at the time points tp[1] and tp[2]. By convention, the graph associated to the sub matrix K_{tp[1],tp[2]} is represented by a directed graph where a directed edge is drawn by an arrow from a vertex in the first graph pointing forwards a vertex in the second graph.

Value

gplot.fglasso returns a list with components:

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

See Also

fglasso function.

Examples

N <- 50
tp <- 3
p <- 3
X <- matrix(rnorm(N * p * tp), N, tp * p)
S <- crossprod(X) / N
model <- list(lag0 = c(s = "c", n = "ut"), lag1 = c(s = "t", n = "t"))
out.fglasso <- fglasso(S = S, model = model, tp = tp, p = p)
gplot(out.fglasso, rhoid = 50, sub.tp1 = "First graph", 
   sub.tp2 = "Second graph")

Plotting Sparse Graphs

Description

gplot.sglasso shows the sequence of graphs estimated by sglasso.

Usage

## S3 method for class 'sglasso'
gplot(object, rhoid, layout = layout.circle, ...)

Arguments

Details

gplot.sglasso shows the sequence of topologically different graphs estimated by sglasso. To specify the layout of the graphs, the user can use any layout function available in the R package igraph. The user can also specify the placement of the vertices by a matrix with two columns and the same number of rows as the number of vertices.

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

See Also

sglasso function.

Examples

N <- 100
p <- 5
X <- matrix(rnorm(N * p), N, p)
S <- crossprod(X)/N
mask <- outer(1:p, 1:p, function(i,j) 0.5^abs(i-j))
mask[1,5] <- mask[1,4] <- mask[2,5] <- NA
mask[5,1] <- mask[4,1] <- mask[5,2] <- NA
out.sglasso_path <- sglasso(S, mask, tol = 1.0e-13)
gplot(out.sglasso_path)
gplot(out.sglasso_path, rhoid = 1:5)

Cross-Validated Kullback-Leibler Divergence

Description

Model selection criterion based on the leave-one-out cross-validated Kullback-Leibler divergence.

Usage

klcv(object, X, scale = 1)

Arguments

Details

klcv function implements the leave-one-out cross-validate Kullback-Leibler divergence criterion proposed in Vujacic et al. (2015). For l_1-penalized Gaussian Graphical Models this measure of goodness-of-fit has the following form

klcv(\rho) = -\frac{\ell(\hat K(\rho))}{N} + \frac{\code{scale}}{2N} gdf(\hat K(\rho)),

where \hat K(\rho) is the glasso estimate of the concentration matrix, \ell(\hat K(\rho)) is the corresponding value of the log-likelihood function, scale is a scale factor for the complexity part, i.e. gdf(\hat K(\rho)), which is defined as

gdf(\hat K(\rho)) = \frac{1}{N-1}\sum_{k=1}^N vec\{(\hat K(\rho)^{-1} - S_k)\circ 1_\rho\}'vec[\hat K(\rho)\{(S-S_k)\circ 1_\rho\}\hat K(\rho)].

In the previous expression S is the empirical variance/covariance matrix, S_k = X_k X_k', 1_\rho is a matrix with entries I(\hat k_{ij}(\rho)\ne 0) and \circ is the Hadamard product operator.

Value

klcv returns an S3 object with calls klcv, i.e. a named list with the following components:

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

References

Vujacic, I., Abbruzzo, A. and Wit, E. C. (2015) A computationally fast alternative to cross-validation in penalized Gaussian graphical models. J. Stat. Comput. Simul.

See Also

sglasso, loglik functions and plot.klcv method.

Examples

N <- 100
p <- 5
X <- matrix(rnorm(N * p), N, p)
S <- crossprod(X) / N
mask <- outer(1:p, 1:p, function(i,j) 0.5^abs(i-j))
mask[1,5] <- mask[1,4] <- mask[2,5] <- NA
mask[5,1] <- mask[4,1] <- mask[5,2] <- NA
out.sglasso_path <- sglasso(S, mask, tol = 1.0e-13)
out.klcv <- klcv(out.sglasso_path, X)
out.klcv

Extract Log-Likelihood

Description

This function extracts the log-likelihood for the sequence of weighted l1-penalized RCON(V, E) models estimated by sglasso function.

Usage

loglik(object, N = 2)

Arguments

Details

Denoted with \psi = (\eta',\theta')' the parameter vector of the structured concentration matrix K(\psi), the log-likelihood function of the RCON(V, E) model is equal, up to a constant, to the following expression

\ell(\psi) = \frac{N}{2}[\log det K(\psi) - tr\{S K(\psi)\}],

where S = N^{-1}\sum_{i=1}^NX_iX_i^T, N is the sample size and X_i is the ith observed p-dimensional vector. Denoted with \hat\psi = (\hat\eta',\hat\theta')' the sglasso estimates, straightforward algebra shows that

\ell(\hat\psi) = \frac{N}{2}[\log det K(\hat\psi) - p + \rho\sum_{m=1}^S w_m|\hat\theta_m|],

where \rho is the tuning parameter and w_m are the weights used to define the weighted l1-norm.

Value

loglik returns a vector containing the log-likelihood computed for the sequence of weighted l1-penalized RCON(V, E).

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

See Also

summary.sglasso method and sglasso function.

Examples

N <- 100
p <- 5
X <- matrix(rnorm(N * p), N, p)
S <- crossprod(X) / N
mask <- outer(1:p, 1:p, function(i, j) 0.5^abs(i-j))
out.sglasso_path <- sglasso(S, mask, nrho = 5, tol = 1.0e-13)
out.sglasso_path
loglik(out.sglasso_path, N = N)
rho <- out.sglasso_path$rho[3]
out.sglasso_single <- sglasso(S, mask, nrho = 1, min_rho = rho, 
   tol = 1.0e-13, algorithm = "ccm")
loglik(out.sglasso_single, N = N)

Neisseria Data Set

Description

This data set contains the gene expression data from a high-resolution time-course experiment besed on the sequenced Neisseria meningitidis serogroup strain B strain MC58. Specifically, the expression level of ten genes is measured at ten different time points. Each column is standardized to have zero mean and standard deviation equal to one.

Usage

data("neisseria")

Plot Method for Leave-One-Out Cross-Validated Kullback-Leibler Divergence

Description

plot.klcv produces a plot to study the sequence of leave-one-out cross-validated Kullback-Leibler divergences computed by klcv.

Usage

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

Arguments

Details

This method function produces a plot showing the sequence of leave-one-out cross-validated Kullback-Leibler as function of the tuning parameter rho. The optimal value of the tuning parameter is identified by a vertical dashed line.

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

See Also

klcv function.

Examples

N <- 100
p <- 5
X <- matrix(rnorm(N * p), N, p)
S <- crossprod(X) / N
mask <- outer(1:p, 1:p, function(i,j) 0.5^abs(i-j))
mask[1,5] <- mask[1,4] <- mask[2,5] <- NA
mask[5,1] <- mask[4,1] <- mask[5,2] <- NA
out.sglasso_path <- sglasso(S, mask, tol = 1.0e-13)
out.klcv <- klcv(out.sglasso_path, X)
plot(out.klcv)

Plot Method for the Weighted l1-Penalized RCON(V, E) Model

Description

plot.sglasso produces two plots to study the sequence of models estimates by sglasso or fglasso.

Usage

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

Arguments

Details

This function produces two different plots. The first one shows the path of the estimated parameters as function of the tuning parameter \rho. In the same way, the second plot shows the path of the weighted scores as function of \rho.

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

See Also

sglasso function and summary.sglasso method.

Examples

N <- 100
p <- 5
X <- matrix(rnorm(N * p), N, p)
S <- crossprod(X)/N
mask <- outer(1:p, 1:p, function(i,j) 0.5^abs(i-j))
mask[1,5] <- mask[1,4] <- mask[2,5] <- NA
mask[5,1] <- mask[4,1] <- mask[5,2] <- NA
out.sglasso_path <- sglasso(S, mask, tol = 1.0e-13)
plot(out.sglasso_path)

Lasso Method for the RCON(V, E) Models

Description

Fit the weighted l1-penalized RCON(V, E) models using a cyclic coordinate algorithm.

Usage

sglasso(S, mask, w = NULL, flg = NULL, min_rho = 1.0e-02, nrho = 50,  
        nstep = 1.0e+05, algorithm = c("ccd","ccm"), truncate = 1e-05, 
        tol = 1.0e-03)

Arguments

Details

The RCON(V, E) model (Hojsgaard et al., 2008) is a kind of restriction of the Gaussian Graphical Model defined using a coloured graph to specify a set of equality constraints on the entries of the concentration matrix. Roughly speaking, a coloured graph implies a partition of the vertex set into R disjoint subsets, called vertex colour classes, and a partition of the edge set into S disjoint subsets, called edge colour classes. At each vertex/edge colour class is associated a specific colour. If we denote by K = (k_{ij}) the concentration matrix, i.e. the inverse of the variance/covariance matrix \Sigma, the coloured graph implies the following equality constraints:

1. k_{ii} = \eta_n for any index i belonging to the nth vertex colour class;

2. k_{ij} = \theta_m for any pair (i,j) belonging to the mth edge colour class.

Denoted with \psi = (\eta',\theta')' the (R+S)-dimensional parameter vector, the concentration matrix can be defined as

K(\psi) = \sum_{n=1}^R\eta_nD_n + \sum_{m=1}^S\theta_mT_m,

where D_n is a diagonal matrix with entries D^n_{ii} = 1 if the index i belongs to the nth vertex colour class and zero otherwise. In the same way, T_m is a symmetrix matrix with entries T^m_{ij} = 1 if the pair (i,j) belongs to the mth edge colour class. Using the previous specification of the concentration matrix, the structured graphical lasso (sglasso) estimator (Abbruzzo et al., 2014) is defined as

\hat\psi = \arg\max_{\psi} \log det K(\psi) - tr\{Sk(\psi)\} - \rho\sum_{m=1}^Sw_m|\theta_m|,

where S is the empirical variance/covariance matrix, \rho is the tuning parameter used to control the ammount of shrinkage and w_m are weights used to define the weighted \ell_1-norm. By default, the sglasso function sets the weights equal to the cardinality of the edge colour classes.

Value

sglasso returns an obejct with S3 class "sglasso", i.e. a named list containing the following components:

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

References

Abbruzzo, A., Augugliaro, L., Mineo, A. M. and Wit, E. C. (2014) Cyclic coordinate for penalized Gaussian Graphical Models with symmetry restrictions. In Proceeding of COMPSTAT 2014 - 21th International Conference on Computational Statistics, Geneva, August 19-24, 2014.

Hojsgaard, S. and Lauritzen, S. L. (2008) Graphical gaussian models with edge and vertex symmetries. J. Roy. Statist. Soc. Ser. B., Vol. 70(5), 1005–1027.

See Also

summary.sglasso, plot.sglasso gplot.sglasso and methods.

The function Kh extracts the estimated sparse structured concentration matrices.

Examples

########################################################
# sglasso solution path
#
## structural zeros:
## there are two ways to specify structural zeros which are 
## related to the kind of mask. If mask is a numeric matrix
## NA is used to identify the structural zero. If mask is a
## character matrix then the structural zeros are specified
## using NA or ".".
N <- 100
p <- 5
X <- matrix(rnorm(N * p), N, p)
S <- crossprod(X) / N
mask <- outer(1:p, 1:p, function(i,j) 0.5^abs(i-j))
mask[1,5] <- mask[1,4] <- mask[2,5] <- NA
mask[5,1] <- mask[4,1] <- mask[5,2] <- NA
mask

out.sglasso_path <- sglasso(S, mask, tol = 1.0e-13)
out.sglasso_path

rho <- out.sglasso_path$rho[20]
out.sglasso <- sglasso(S, mask, nrho = 1, min_rho = rho, tol = 1.0e-13, algorithm = "ccm")
out.sglasso

out.sglasso_path$theta[, 20]
out.sglasso$theta[, 1]

Internal sglasso functions

Description

Internal sglasso functions

Details

These are not intended for use by users

Author(s)

Luigi Augugliaro

Summarizing sglasso Fits

Description

summary method for class "sglasso".

Usage

## S3 method for class 'sglasso'
summary(object, N, k = c("bic","aic"), 
        digits = max(3, getOption("digits") - 3), ...)

Arguments

Details

summary.sglasso gives us information about the sequence of models estimated by the sglasso estimator. To select the best model, summary method uses a measure of Goodness-of-Fit (GoF) defined as follows:

-2\ell(\hat\psi) + k \times df,

where \ell(\hat\psi) is the log-likelihood of the estimated weighted l1-penalized RCON(V, E) model, df is the number of nonzero estimated parameters and k is a non-negative value used to weight the complexity part in the measure of goodness-of-fit. By default the summary method computes the BIC criterion to select the best model (k = "bic"). The AIC criterion can be easily computed setting k = "aic". The user can also define other measures of goodness-of-fit specifying k as any non-negative value.

The output of the summary method is divided in two sections. First section shows the call producing the argument object followed by a data.frame. The column named rho shows the sequence of the \rho values used to compute the solution curve, while the column log-lik shows the corresponding values of the log-likelihood function. The remaining columns show the number of estimated non-zero parameters, the values of the GoF and the asscoated ranking of the estimated models. Finally, the second section shows the estimated parameters of the best model identified by the used GoF criterion. Informations about the algorithm and the corresponding convergence are also provided.

Value

A list with components table and theta_gof is silently returned. The table component is the data.frame previously described while the component theta_gof is the vector of the estimated parameters corresponding to the best models identified by the GoF criterion.

Author(s)

Luigi Augugliaro
Maintainer: Luigi Augugliaro luigi.augugliaro@unipa.it

See Also

sglasso and loglik functions.

Examples

N <- 100
p <- 5
X <- matrix(rnorm(N * p), N, p)
S <- crossprod(X) / N
mask <- outer(1:p, 1:p, function(i,j) 0.5^abs(i-j))
mask[1,5] <- mask[1,4] <- mask[2,5] <- NA
mask[5,1] <- mask[4,1] <- mask[5,2] <- NA
out.sglasso_path <- sglasso(S, mask, tol = 1.0e-13)
summary(out.sglasso_path, N)
rho <- out.sglasso_path$rho[20]
out.sglasso <- sglasso(S, mask, nrho = 1, min_rho = rho, tol = 1.0e-13)
summary(out.sglasso, N)