Package 'rags2ridges'

Description

Package contains proper L2-penalized ML estimators for the precision matrix as well as supporting functions to employ these estimators in a (integrative or meta-analytic) graphical modeling setting.

Details

The main function of the package is ridgeP which enables archetypal and proper alternative ML ridge estimation of the precision matrix. The alternative ridge estimators can be found in van Wieringen and Peeters (2015) and encapsulate both target and non-target shrinkage for the multivariate normal precision matrix. The estimators are analytic and enable estimation in large $p$ small $n$ settings. Supporting functions to employ these estimators in a graphical modeling setting are also given. These supporting functions enable, a.o., the determination of the optimal value of the penalty parameter, the determination of the support of a shrunken precision estimate, as well as various visualization options.

The package has a modular setup. The core module (rags2ridges.R) contains the functionality stated above. The fused module (rags2ridgesFused.R) extends the functionality of the core module to the joint estimation of multiple precision matrices from (aggregated) high-dimensional data consisting of distinct classes. The result is a targeted fused ridge estimator that is of use when the precision matrices of the constituent classes are believed to chiefly share the same structure while potentially differing in a number of locations of interest. The fused module also contains supporting functions for integrative or meta-analytic Gaussian graphical modeling. The third module is the miscellaneous module (rags2RidgesMisc.R) which contains assorted hidden functions.

Function overview core module:

• Function for (proper) ridge estimation of the precision matrix

  • ridgeP

• Functions for penalty parameter selection

  • CNplot

  • optPenalty.aLOOCV

  • optPenalty.kCV

  • optPenalty.kCVauto

• Functions for loss/entropy/fit evaluation

  • evaluateSfit

  • KLdiv

  • loss

• Functions for block-independence testing

  • GGMblockNullPenalty

  • GGMblockTest

• Function for support determination

  • sparsify

• Functions for (network) visualization

  • edgeHeat

  • ridgePathS

  • Ugraph

• Functions for topology statistics

  • GGMmutualInfo

  • GGMnetworkStats

  • GGMpathStats

• Wrapper function

  • fullMontyS

• Support functions

  • adjacentMat

  • covML

  • covMLknown

  • default.target

  • evaluateS

  • pcor

  • symm

Function overview fused module:

• Function for targeted fused ridge estimation of multiple precision matrices

  • ridgeP.fused

• Function for fused penalty parameter selection

  • optPenalty.fused

• Functions for loss/entropy/fit evaluation

  • KLdiv.fused

  • NLL

• Function for testing the necessity of fusion

  • fused.test

• Function for support determination

  • sparsify.fused

• Functions for topology statistics

  • GGMnetworkStats.fused

  • GGMpathStats.fused

• Support functions

  • createS

  • default.penalty

  • default.target.fused

  • getKEGGPathway

  • isSymmetricPD

  • is.Xlist

  • kegg.target

  • plot.ptest

  • pooledS

  • print.optPenaltyFusedGrid

  • print.ptest

  • rmvnormal

Calls of interest to miscellaneous module:

• rags2ridges:::.TwoCents() ~~(Unsolicited advice)

• rags2ridges:::.Brooke() ~~(Endorsement)

• rags2ridges:::.JayZScore() ~~(The truth)

• rags2ridges:::.theHoff() ~~(Wish)

• rags2ridges:::.rags2logo() ~~(Warm welcome)

Author(s)

Carel F.W. Peeters, Anders Ellern Bilgrau, Wessel, N. van Wieringen
Maintainer: Carel F.W. Peeters <[email protected]>

References

Peeters, C.F.W., Bilgrau, A.E., and van Wieringen, W.N. (2022). rags2ridges: A One-Stop-l2-Shop for Graphical Modeling of High-Dimensional Precision Matrices. Journal of Statistical Software, vol. 102(4): 1-32.

Bilgrau, A.E., Peeters, C.F.W., Eriksen, P.S., Boegsted, M., and van Wieringen, W.N. (2020). Targeted Fused Ridge Estimation of Inverse Covariance Matrices from Multiple High-Dimensional Data Classes. Journal of Machine Learning Research, 21(26): 1-52. Also available as arXiv:1509.07982v2 [stat.ME].

Peeters, C.F.W., van de Wiel, M.A., & van Wieringen, W.N. (2020). The Spectral Condition Number Plot for Regularization Parameter Evaluation. Computational Statistics, 35: 629-646. Also available as arXiv:1608.04123 [stat.CO].

van Wieringen, W.N. & Peeters, C.F.W. (2016). Ridge Estimation of Inverse Covariance Matrices from High-Dimensional Data. Computational Statistics & Data Analysis, vol. 103: 284-303. Also available as arXiv:1403.0904v3 [stat.ME].

van Wieringen, W.N. & Peeters, C.F.W. (2015). Application of a New Ridge Estimator of the Inverse Covariance Matrix to the Reconstruction of Gene-Gene Interaction Networks. In: di Serio, C., Lio, P., Nonis, A., and Tagliaferri, R. (Eds.) 'Computational Intelligence Methods for Bioinformatics and Biostatistics'. Lecture Notes in Computer Science, vol. 8623. Springer, pp. 170-179.

Description

ADdata contains 3 objects related to metabolomics data on patients with Alzheimer's Disease (AD).

Details

ADmetabolites is a matrix containing metabolic expressions of 230 metabolites (rows) on 127 samples (columns).

sampleInfo is a data.frame containing information on the samples. Information pertains to diagnosis: AD class 1 or AD class 2.

variableInfo is a data.frame containing information on the metabolic features. Information pertains to compound families: Amines, organic acids, lipids, and oxidative stress compounds.

See description.

Author(s)

Carel F.W. Peeters <[email protected]>

Source

de Leeuw, F., Peeters, C.F.W., Kester, M.I., Harms, A.C., Struys, E., Hankemeijer, T., van Vlijmen, H.W.T., van Duijn, C.M., Scheltens, P., Demirkan, A., van de Wiel, M.A., van der Flier, W.M., and Teunissen, C.E. (2017). Blood-based metabolic signatures in Alzheimer's Disease. Alzheimer's & Dementia: Diagnosis, Assessment & Disease Monitoring, 8: 196-207.

Examples

data(ADdata)

## Look at sample information
sampleInfo

## Look at feature information
variableInfo

Description

Function that transforms a real matrix into an adjacency matrix. Intended use: Turn sparsified precision matrix into an adjacency matrix for undirected graphical representation.

Usage

adjacentMat(M, diag = FALSE)

Arguments

Value

Function returns an adjacency matrix.

Author(s)

Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

See Also

ridgeP, covML, sparsify, edgeHeat, Ugraph

Examples

## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]
Cx <- covML(X)

## Obtain regularized precision matrix
P <- ridgeP(Cx, lambda = 10, type = "Alt")

## Obtain sparsified partial correlation matrix
PC0 <- sparsify(P, threshold = "localFDR", FDRcut = .8)

## Obtain adjacency matrix
adjacentMat(PC0$sparsePrecision)

Description

Function that visualizes the spectral condition number of the regularized precision matrix against the domain of the regularization parameter. The function can be used to heuristically determine an acceptable (minimal) value for the penalty parameter.

Usage

CNplot(
  S,
  lambdaMin,
  lambdaMax,
  step,
  type = "Alt",
  target = default.target(S, type = "DUPV"),
  norm = "2",
  Iaids = FALSE,
  vertical = FALSE,
  value = 1e-100,
  main = "",
  nOutput = FALSE,
  verbose = TRUE,
  suppressChecks = FALSE
)

Arguments

Details

Under certain target choices the proposed ridge estimators (see ridgeP) are rotation equivariant, i.e., the eigenvectors of $\mathbf{S}$ are left intact. Such rotation equivariant situations help to understand the effect of the ridge penalty on the precision estimate: The effect can be understood in terms of shrinkage of the eigenvalues of the unpenalized precision estimate $\mathbf{S}^{-1}$. Maximum shrinkage implies that all eigenvalues are forced to be equal (in the rotation equivariant situation). The spectral condition number w.r.t. inversion (ratio of maximum to minimum eigenvalue) of the regularized precision matrix may function as a heuristic in determining the (minimal) value of the penalty parameter. A matrix with a high condition number is near-singular (the relative distance to the set of singular matrices equals the reciprocal of the condition number; Demmel, 1987) and its inversion is numerically unstable. Such a matrix is said to be ill-conditioned. Numerically, ill-conditioning will mean that small changes in the penalty parameter lead to dramatic changes in the condition number. From a numerical point of view one can thus track the domain of the penalty parameter for which the regularized precision matrix is ill-conditioned. When plotting the condition number against the (domain of the) penalty parameter, there is a point of relative stabilization (when working in the $p > n$ situation) that can be characterized by a leveling-off of the acceleration along the curve when plotting the condition number against the (chosen) domain of the penalty parameter. This suggest the following fast heuristic for determining the (minimal) value of the penalty parameter: The value of the penalty parameter for which the spectral condition number starts to stabilize may be termed an acceptable (minimal) value.

The function outputs a graph of the (spectral) matrix condition number over the domain [lambdaMin, lambdaMax]. When norm = "2" the spectral condition number is calculated. It is determined by exact calculation using the spectral decomposition. For most purposes this exact calculation is fast enough, especially when considering rotation equivariant situations (see ridgeP). For such situations the amenities for fast eigenvalue calculation as provided by RSpectra are used internally. When exact computation of the spectral condition number is deemed too costly one may approximate the computationally friendly L1-condition number. This approximation is accessed through the rcond function (Anderson et al. 1999).

When Iaids = TRUE the basic condition number plot is amended/enhanced with two additional plots (over the same domain of the penalty parameter as the basic plot): The approximate loss in digits of accuracy (for the operation of inversion) and an approximation to the second-order derivative of the curvature in the basic plot. These interpretational aids can enhance interpretation of the basic condition number plot and may support choosing a value for the penalty parameter (see Peeters, van de Wiel, & van Wieringen, 2016). When vertical = TRUE a vertical line is added at the constant value. This option can be used to assess if the optimal penalty obtained by, e.g., the routines optPenalty.LOOCV or optPenalty.aLOOCV, has led to a precision estimate that is well-conditioned.

Value

The function returns a graph. If nOutput = TRUE the function also returns an object of class list:

Note

The condition number of a (regularized) covariance matrix is equivalent to the condition number of its corresponding inverse, the (regularized) precision matrix. Please note that the target argument (for Type I ridge estimators) is assumed to be specified in precision terms. In case one is interested in the condition number of a Type I estimator under a covariance target, say $\mathbf{\Gamma}$, then just specify target = solve($\mathbf{\Gamma}$).

Author(s)

Carel F.W. Peeters <[email protected]>

References

Anderson, E, Bai, Z., ..., Sorenson, D. (1999). LAPACK Users' Guide (3rd ed.). Philadelphia, PA, USA: Society for Industrial and Applied Mathematics.

Demmel, J.W. (1987). On condition numbers and the distance to the nearest ill-posed problem. Numerische Mathematik, 51: 251–289.

Peeters, C.F.W., van de Wiel, M.A., & van Wieringen, W.N. (2020). The spectral condition number plot for regularization parameter evaluation. Computational Statistics, 35: 629–646.

See Also

covML, ridgeP, optPenalty.LOOCV, optPenalty.aLOOCV, default.target

Examples

## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]
Cx <- covML(X)

## Assess spectral condition number across grid of penalty parameter
CNplot(Cx, lambdaMin = .0001, lambdaMax = 50, step = 1000)

## Include interpretational aids
CNplot(Cx, lambdaMin = .0001, lambdaMax = 50, step = 1000, Iaids = TRUE)

Description

Function that searches for and visualizes community-structures in graphs.

Usage

Communities(
  P,
  graph = TRUE,
  lay = "layout_with_fr",
  coords = NULL,
  Vsize = 15,
  Vcex = 1,
  Vcolor = "orangered",
  VBcolor = "darkred",
  VLcolor = "black",
  main = ""
)

Arguments

Details

Communities in a network are groups of vertices (modules) that are densely connected within. Community search is performed by the Girvan-Newman algorithm (Newman and Girvan, 2004).

When graph = TRUE the community structure in the graph is visualized. The default layout is according to the Fruchterman-Reingold algorithm (1991). Most layout functions supported by igraph are supported (the function is partly a wrapper around certain igraph functions). The igraph layouts can be invoked by a character that mimicks a call to a igraph layout functions in the lay argument. When using lay = NULL one can specify the placement of vertices with the coords argument. The row dimension of this matrix should equal the number of vertices. The column dimension then should equal 2 (for 2D layouts) or 3 (for 3D layouts). The coords argument can also be viewed as a convenience argument as it enables one, e.g., to layout a graph according to the coordinates of a previous call to Ugraph. If both the the lay and the coords arguments are not NULL, the lay argument takes precedence. Communities are indicated by color markings.

Value

An object of class list:

When graph = TRUE the function also returns a graph.

Author(s)

Carel F.W. Peeters <[email protected]>

References

Csardi, G. and Nepusz, T. (2006). The igraph software package for complex network research. InterJournal, Complex Systems 1695. http://igraph.sf.net

Fruchterman, T.M.J., and Reingold, E.M. (1991). Graph Drawing by Force-Directed Placement. Software: Practice & Experience, 21: 1129-1164.

Newman, M. and Girvan, M. (2004). Finding and evaluating community structure in networks. Physical Review E, 69: 026113.

See Also

Ugraph

Examples

## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]

## Obtain regularized precision under optimal penalty
OPT <- optPenalty.LOOCV(X, lambdaMin = .5, lambdaMax = 30, step = 100)

## Determine support regularized standardized precision under optimal penalty
PC0 <- sparsify(symm(OPT$optPrec), threshold = "localFDR")$sparseParCor

## Search and visualize communities
Commy <- Communities(PC0)

Description

This function is now deprecated. Please use CNplot instead.

Usage

conditionNumberPlot(
  S,
  lambdaMin,
  lambdaMax,
  step,
  type = "Alt",
  target = default.target(S),
  norm = "2",
  digitLoss = FALSE,
  rlDist = FALSE,
  vertical = FALSE,
  value,
  main = TRUE,
  nOutput = FALSE,
  verbose = TRUE
)

Arguments

Details

See CNplot.

Value

The function returns a graph. If nOutput = TRUE the function also returns an object of class list:

Author(s)

Carel F.W. Peeters <[email protected]>

See Also

CNplot

Examples

## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]
Cx <- covML(X)

## Assess spectral condition number across grid of penalty parameter
conditionNumberPlot(Cx, lambdaMin = .0001, lambdaMax = 50, step = 1000)

Description

Function that gives the maximum likelihood estimate of the covariance matrix.

Usage

covML(Y, cor = FALSE)

Arguments

Details

The function gives the maximum likelihood (ML) estimate of the covariance matrix. The input matrix Y assumes that the variables are represented by the columns. Note that when the input data is standardized, the ML covariance matrix of the scaled data is computed. If a correlation matrix is desired, use cor = TRUE.

Value

Function returns the maximum likelihood estimate of the covariance matrix. In case cor = TRUE, the correlation matrix is returned.

Author(s)

Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

See Also

ridgeP

Examples

## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]

## Obtain ML estimate covariance matrix
Cx <- covML(X)

## Obtain correlation matrix
Cx <- covML(X, cor = TRUE)

Description

Function that performs maximum likelihood estimation of the covariance matrix, with various types of assumptions on its structure.

Usage

covMLknown(
  Y,
  covMat = NULL,
  corMat = NULL,
  corType = "none",
  varType = "none",
  nInit = 100
)

Arguments

Details

The function gives the maximum likelihood estimate of the covariance matrix. The input matrix Y assumes that the variables are represented by the columns.

When simultaneously covMat=NULL, corMat=NULL, corType="none" and varType="none" the covML-function is invoked and the regular maximum likelihood estimate of the covariance matrix is returned.

Value

The maximum likelihood estimate of the covariance matrix under the specified assumptions on its structure.

Author(s)

Wessel N. van Wieringen, Carel F.W. Peeters <[email protected]>

See Also

covML

Examples

## Obtain some data
p = 10
n = 100
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:10] = letters[1:10]

## Obtain maximum likelihood estimate covariance matrix
Cx <- covMLknown(X, corType="equi", varType="common")

Description

Simulate data from a p-dimensional (zero-mean) gaussian graphical model (GGM) with a specified (or random) topology and return the sample covariance matrix or matrices. Can also return the original simulated data or underlying precision matrix.

Usage

createS(
  n,
  p,
  topology = "identity",
  dataset = FALSE,
  precision = FALSE,
  nonzero = 0.25,
  m = 1L,
  banded.n = 2L,
  invwishart = FALSE,
  nu = p + 1,
  Plist
)

Arguments

Details

The data is simulated from a zero-mean p-dimensional multivariate gaussian distribution with some precision matrix determined by the argument topology which defines the GGM. If precision is TRUE the population precision matrix is returned. This is useful to see what the actual would-be-used precision matrices are. The available values of topology are described below. Unless otherwise stated the diagonal entries are always one. If m is 2 or greater block diagonal precision matrices are constructed and used.

• "identity": uses the identity matrix (diag(p)) as precision matrix. Corresponds to no conditional dependencies.

• "star": simulate from a star topology. Within each block the first node is selected as the "hub". The off-diagonal entries $(1,j)$ and $(j,1)$ values taper off with the value $1/(j + 1)$.

• "clique": simulate from clique topology where each block is a complete graph with off-diagonal elements equal to nonzero.

• "complete": alias for (and identical to) "clique".

• "chain": simulate from a chain topology where the precision matrix is a tridiagonal matrix with off-diagonal elements (in each block) given by argument nonzero.

• "banded": precision elements (i,j) are given by $1/(|i-j|+1)$ if $|i-j|$ is less than or equal to banded.n and zero otherwise.

• "scale-free": The non-zero pattern of each block is generated by a Barabassi random graph. Non-zero off-diagonal values are given by nonzero. Gives are very "hubby" network.

• "Barabassi": alias for "scale-free".

• "small-world": The non-zero pattern of each block is generated by a 1-dimensional Watts-Strogatz random graph with banded.n starting neighbors and $5\%$ probability of rewiring. Non-zero off-diagonal values are given by nonzero. Gives are very "bandy" network.

• "Watts-Strogatz": alias for "small-world"

• "random-graph": The non-zero pattern of each block is generated by a Erdos-Renyi random graph where each edge is present with probability $1/p$. Non-zero off-diagonal values are given by nonzero.

• "Erdos-Renyi": alias for "random-graph"

When n has length greater than 1, the datasets are generated i.i.d. given the topology and number of blocks.

Arguments invwishart and nu allows for introducing class homogeneity. Large values of nu imply high class homogeneity. nu must be greater than p + 1. More precisely, if invwishart == TRUE then the constructed precision matrix is used as the scale parameter in an inverse Wishart distribution with nu degrees of freedom. Each class covariance is distributed according to this inverse Wishart and independent.

Value

The returned type is dependent on n and covariance. The function generally returns a list of numeric matrices with the same length as n. If covariance is FALSE the simulated datasets with size n[i] by p are given in the i entry of the output. If covariance is TRUE the p by p sample covariances of the datasets are given. When n has length 1 the list structure is dropped and the matrix is returned.

Author(s)

Anders E. Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

Examples

## Generate some simple sample covariance matrices
createS(n = 10, p = 3)
createS(n = c(3, 4, 5), p = 3)
createS(n = c(32, 55), p = 7)

## Generate some datasets and not sample covariance matrices
createS(c(3, 4), p = 6, dataset = TRUE)

## Generate sample covariance matrices from other topologies:
A <- createS(2000, p = 4, topology = "star")
round(solve(A), 3)
B <- createS(2000, p = 4, topology = "banded", banded.n = 2)
round(solve(B), 3)
C <- createS(2000, p = 4, topology = "clique")  # The complete graph (as m = 1)
round(solve(C), 3)
D <- createS(2000, p = 4, topology = "chain")
round(solve(D), 3)

## Generate smaple covariance matrices from block topologies:
C3 <- createS(2000, p = 10, topology = "clique", m = 3)
round(solve(C3), 1)
C5 <- createS(2000, p = 10, topology = "clique", m = 5)
round(solve(C5), 1)

## Can also return the precision matrix to see what happens
## m = 2 blocks, each "banded" with 4 off-diagonal bands
round(createS(1, 12, "banded", m = 2, banded.n = 4, precision = TRUE), 2)

## Simulation using graph-games
round(createS(1, 10, "small-world", precision = TRUE), 2)
round(createS(1, 5, "scale-free", precision = TRUE), 2)
round(createS(1, 5, "random-graph", precision = TRUE), 2)

## Simulation using inverse Wishart distributed class covariance
## Low class homogeneity
createS(n = c(10,10), p = 5, "banded", invwishart = TRUE, nu = 10)
## Extremely high class homogeneity
createS(n = c(10,10), p = 5, "banded", invwishart = TRUE, nu = 1e10)

# The precision argument can again be used to see the actual realised class
# precision matrices used when invwishart = TRUE.

# The Plist argument is used to reuse old precision matrices or
# user-generated ones
P <- createS(n = 1, p = 5, "banded", precision = TRUE)
lapply(createS(n = c(1e5, 1e5), p = 5, Plist = list(P, P+1)), solve)

Description

Function that constructs default or commonly use penalty matrices according to a (factorial) study design. The constructed penalty matrix can be used directly in optPenalty.fused.auto or serve as basis for modification.

Usage

default.penalty(
  G,
  df,
  type = c("Complete", "CartesianEqual", "CartesianUnequal", "TensorProd")
)

Arguments

Details

The type gives a number of common choices for the penalty matrix:

• 'Complete' is the complete penalty graph with equal penalties.

• 'CartesianEqual' corresponds to a penalizing along each "direction" of factors with a common penalty. The choice is named Cartesian as it is the Cartesian graph product of the complete penalty graphs for the individual factors.

• 'CartesianUnequal' corresponds to a penalizing each direction of factors with individual penalties.

• 'TensorProd' correspond to penalizing the "diagonals" only. It is equivalent to the graph tensor products of the complete graphs for each individual factor.

Value

Returns a G by G character matrix which specify the class of penalty graphs to be used. The output is suitable as input for the penalty matrix used in optPenalty.fused.auto.

Author(s)

Anders E. Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

References

Bilgrau, A.E., Peeters, C.F.W., Eriksen, P.S., Boegsted, M., and van Wieringen, W.N. (2020). Targeted Fused Ridge Estimation of Inverse Covariance Matrices from Multiple High-Dimensional Data Classes. Journal of Machine Learning Research, 21(26): 1-52.

See Also

ridgeP.fused, optPenalty.fused, default.target

Examples

# Handling one-way designs
  default.penalty(2)
  default.penalty(4)
  Slist <- vector("list", 6)
  default.penalty(Slist)   # The function uses only the length of the list
  df0 <- expand.grid(Factor = c("lvl1", "lvl2"))
  default.penalty(df0)

  # A more elaborate example
  df1 <- expand.grid(DS = c("DS1", "DS2", "DS3"), ER = c("ER+", "ER-"))

  # Usage (various interface demonstrations)
  default.penalty(6, df1, type = "Complete")
  default.penalty(6, type = "CartesianEqual")  # GIVES WARNING
  default.penalty(6, df1, type = "CartesianEqual")
  default.penalty(Slist, df1, type = "CartesianEqual")
  default.penalty(6, df1, type = "CartesianUnequal")
  default.penalty(df1)

  # A 2 by 2 by 2 design
  df2 <- expand.grid(A = c("A1", "A2"), B = c("B1", "B2"), C = c("C1", "C3"))
  default.penalty(df2)
  default.penalty(df2, type = "CartesianEqual")
  default.penalty(df2, type = "CartesianUnequal")

Description

Function that generates a (data-driven) default target for usage in (type I) ridge shrinkage estimation of the precision matrix (see ridgeP). The target that is generated is to be understood in precision terms. Most options for target generation result in a target that implies a situation of rotation equivariant estimation (see ridgeP).

Usage

default.target(S, type = "DAIE", fraction = 1e-04, const)

Arguments

Details

The function can generate the following default target matrices:

• DAIE: Diagonal matrix with average of inverse nonzero eigenvalues of S as entries;

• DIAES: Diagonal matrix with inverse of average of eigenvalues of S as entries;

• DUPV: Diagonal matrix with unit partial variance as entries (identity matrix);

• DAPV: Diagonal matrix with average of inverse variances of S as entries;

• DCPV: Diagonal matrix with constant partial variance as entries. Allows one to use other constant than DAIE, DIAES, DUPV, DAPV, and in a sense Null;

• DEPV: Diagonal matrix with the inverse variances of S as entries;

• Null: Null matrix.

The targets DUPV, DCPV, and Null are not data-driven in the sense that the input matrix S only provides information on the size of the desired target. The targets DAIE, DIAES, DAPV, and DEPV are data-driven in the sense that the input matrix S provides the information for the diagonal entries. The argument fraction is only used when type = "DAIE". The argument const is only used when type = "DCPV". All types except DEPV and Null lead to rotation equivariant alternative and archetypal Type I ridge estimators. The target Null also leads to a rotation equivariant alternative Type II ridge estimator (see ridgeP). Note that the DIAES, DAPV, and DEPV targets amount to the identity matrix when the sample covariance matrix S is standardized to be the correlation matrix. The same goes, naturally, for the DCPV target when const is specified to be 1.

Value

Function returns a target matrix.

Author(s)

Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

References

van Wieringen, W.N. & Peeters, C.F.W. (2016). Ridge Estimation of Inverse Covariance Matrices from High-Dimensional Data, Computational Statistics & Data Analysis, vol. 103: 284-303. Also available as arXiv:1403.0904v3 [stat.ME].

See Also

ridgeP, covML

Examples

## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]
Cx <- covML(X)

## Obtain default diagonal target matrix
default.target(Cx)

Description

Generates a list of (data-driven) targets to use in fused ridge estimation. Simply a wrapper for default.target.

Usage

default.target.fused(Slist, ns, type = "DAIE", by, ...)

Arguments

Value

A list of $K$ covariance target matrices of the same size.

Author(s)

Anders E. Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

See Also

default.target

Examples

# Make some toy data
ns <- c(3, 4)  # Two classes with sample size 3 and 4
Slist <- createS(ns, p = 3)  # Generate two 3-dimensional covariance matrices
Slist

# Different choices:
default.target.fused(Slist, ns)
default.target.fused(Slist, ns, by = seq_along(Slist)) # The same as before
default.target.fused(Slist, ns, type = "Null")
default.target.fused(Slist, ns, type = "DAPV")
default.target.fused(Slist, ns, type = "DAPV", by = rep(1, length(Slist)))


# Make some (more) toy data
ns <- c(3, 4, 6, 7)  # Two classes with sample size 3 and 4
Slist <- createS(ns, p = 2)  # Generate four 2-dimensional covariance matrices

# Use the same target in class 1 and 2, but another in class 3 and 4:
default.target.fused(Slist, ns, by = c("A", "A", "B", "B"))

Description

Function visualizing the differential graph, i.e., the network of edges that are unique for 2 class-specific graphs over the same vertices

Usage

DiffGraph(
  P1,
  P2,
  lay = "layout_with_fr",
  coords = NULL,
  Vsize = 15,
  Vcex = 1,
  Vcolor = "orangered",
  VBcolor = "darkred",
  VLcolor = "black",
  P1color = "red",
  P2color = "green",
  main = ""
)

Arguments

Details

Say you have 2 class-specific precision matrices that are estimated over the same variables/features. This function visualizes in a single graph the edges that are unique to the respective classes. Hence, it gives the differential graph. Edges unique to P1 are colored according to P1color. Edges unique to P2 are colored according to P2color. Dashed lines indicate negative precision elements while solid lines indicate positive precision elements.

The default layout is according to the Fruchterman-Reingold algorithm (1991). Most layout functions supported by igraph are supported (the function is partly a wrapper around certain igraph functions). The igraph layouts can be invoked by a character that mimicks a call to a igraph layout functions in the lay argument. When using lay = NULL one can specify the placement of vertices with the coords argument. The row dimension of this matrix should equal the number of vertices. The column dimension then should equal 2 (for 2D layouts) or 3 (for 3D layouts). The coords argument can also be viewed as a convenience argument as it enables one, e.g., to layout a graph according to the coordinates of a previous call to Ugraph. If both the the lay and the coords arguments are not NULL, the lay argument takes precedence.

Value

The function returns a graph.

Author(s)

Carel F.W. Peeters <[email protected]>

References

Csardi, G. and Nepusz, T. (2006). The igraph software package for complex network research. InterJournal, Complex Systems 1695. http://igraph.sf.net

Fruchterman, T.M.J., and Reingold, E.M. (1991). Graph Drawing by Force-Directed Placement. Software: Practice & Experience, 21: 1129-1164.

See Also

Ugraph

Examples

## Obtain some (high-dimensional) data, class 1
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]

## Obtain some (high-dimensional) data, class 2
set.seed(123456)
X2 = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X2)[1:25] = letters[1:25]

## Obtain regularized precision under optimal penalty, classes 1 and 2
OPT  <- optPenalty.LOOCV(X, lambdaMin = .5, lambdaMax = 30, step = 100)
OPT2 <- optPenalty.LOOCV(X2, lambdaMin = .5, lambdaMax = 30, step = 100)

## Determine support regularized standardized precision under optimal penalty
PC0  <- sparsify(symm(OPT$optPrec), threshold = "localFDR")$sparseParCor
PC02 <- sparsify(symm(OPT2$optPrec), threshold = "localFDR")$sparseParCor

## Visualize differential graph
DiffGraph(PC0, PC02)

Description

Function that visualizes a (precision) matrix as a heatmap. May be used to assess visually the elements of a single (possibly sparsified precision) matrix. May also be used in assessing the performance of edge selection techniques.

Usage

edgeHeat(
  M,
  lowColor = "blue",
  highColor = "red",
  textsize = 10,
  diag = TRUE,
  legend = TRUE,
  main = ""
)

Arguments

Details

This function utilizes ggplot2 (Wickham, 2009) to visualize a matrix as a heatmap: a false color plot in which the individual matrix entries are represented by colors. lowColor determines the color scale for matrix entries in the negative range. highColor determines the color scale for matrix entries in the positive range. For the colors supported by the arguments lowColor and highColor, see https://www.nceas.ucsb.edu/sites/default/files/2020-04/colorPaletteCheatsheet.pdf. White entries in the plot represent the midscale value of 0. One can opt to set the diagonal entries to the midscale color of white when one is interested in (heatmapping) the off-diagonal elements only. To achieve this, set diag = FALSE. Naturally, the diag argument is only used when the input matrix M is a square matrix.

The intended use of the function is to visualize a, possibly sparsified, precision matrix as a heatmap. The function may also be used, in a graphical modeling setting, to assess the performance of edge selection techniques. However, the function is quite general, in the sense that it can represent any matrix as a heatmap.

Author(s)

Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

References

Wickham, H. (2009). ggplot2: elegant graphics for data analysis. New York: Springer.

See Also

covML, ridgeP, sparsify

Examples

## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]
Cx <- covML(X)

## Obtain regularized precision matrix
P <- ridgeP(Cx, lambda = 10, type = "Alt")

## Obtain sparsified partial correlation matrix
PC0 <- sparsify(P, threshold = "localFDR", FDRcut = .8)$sparseParCor

## Visualize sparsified partial correlation matrix as heatmap
edgeHeat(PC0)

Description

Function that evaluates various numerical properties of a square input matrix. The intended use is to evaluate the various numerical properties of what is assumed to be a covariance matrix. Another use is to evaluate the various numerical properties of a (regularized) precision matrix.

Usage

evaluateS(S, verbose = TRUE)

Arguments

Details

The function evaluates various numerical properties of a covariance or precision input matrix. The function assesses if the input matrix is symmetric, if all its eigenvalues are real, if all its eigenvalues are strictly positive, and if it is a diagonally dominant matrix. In addition, the function calculates the trace, the determinant, and the spectral condition number of the input matrix. See, e.g., Harville (1997) for more details on the mentioned (numerical) matrix properties.

Value

Author(s)

Wessel N. van Wieringen, Carel F.W. Peeters <[email protected]>

References

Harville, D.A.(1997). Matrix algebra from a statistician's perspective. New York: Springer-Verlag.

See Also

covML, ridgeP

Examples

## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]
Cx <- covML(X)

## Evaluate numerical properties covariance matrix
## Obtain, e.g., value trace
Seval <- evaluateS(Cx); Seval
Seval$trace

## Evaluate numerical properties precision matrix after regularization
P <- ridgeP(Cx, lambda = 10, type = 'Alt')
Peval <- evaluateS(P); Peval

Description

Function aiding the visual inspection of the fit of an estimated (possibly regularized) precision matrix vis-a-vis the sample covariance matrix.

Usage

evaluateSfit(
  Phat,
  S,
  diag = FALSE,
  fileType = "pdf",
  nameExt = "",
  dir = getwd()
)

Arguments

Details

The function outputs various visualizations to aid the visual inspection of an estimated and possibly regularized precision matrix vis-a-vis the sample covariance matrix. The inverse of the estimated precision matrix P is taken to represent the estimated covariance matrix. The function then outputs a QQ-plot and a heatmap of the observed covariances against the estimated ones. The heatmap has the estimated covariances as lower-triangular elements and the observed covariances as the upper-triangular elements. The function outputs analogous plots for the estimated and observed correlations. In case the observed covariance matrix S is non-singular also a QQ-plot an a heatmap are generated for the estimated and observed partial correlations.

The function generates files with extension fileType under default output names. These files are stored in the directory dir (default is the working directory). To avoid overwriting of files when working in a single directory one may employ the argument nameExt. By using nameExt the default output names are extended with a character of choice.

Author(s)

Wessel N. van Wieringen, Carel F.W. Peeters <[email protected]>

See Also

ridgeP, covML

Examples

## Not run: 
## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]
Cx <- covML(X)

## Obtain regularized precision matrix
P <- ridgeP(Cx, lambda = 10, type = 'Alt')

## Evaluate visually fit of regularized precision matrix vis-a-vis sample covariance
evaluateSfit(P, Cx, diag = FALSE, fileType = "pdf", nameExt = "test")
## End(Not run)

Description

Function that forms a wrapper around certain rags2ridges functionalities. More specifically, it (automatically) invokes functionalities to get from high-dimensional data to a penalized precision estimate, to the corresponding conditional independence graph and topology summaries.

Usage

fullMontyS(
  Y,
  lambdaMin,
  lambdaMax,
  target = default.target(covML(Y)),
  dir = getwd(),
  fileTypeFig = "pdf",
  FDRcut = 0.9,
  nOutput = TRUE,
  verbose = TRUE
)

Arguments

Details

The wrapper always uses the alternative ridge precision estimator (see ridgeP) with target as the target matrix. The optimal value for the penalty parameter is determined by employing Brent's method to the calculation of a cross-validated negative log-likelihood score (see optPenalty.LOOCVauto). The support of the regularized precision matrix is determined by way of local FDR thresholding (see sparsify). The corresponding conditional independence graph is visualized using Ugraph with type = "fancy". This visualization as well as the calculation of network statistics (see GGMnetworkStats) is based on the standardization of the regularized and sparsified precision matrix to a partial correlation matrix.

Value

The function stores in the specified directory dir a condition number plot (either .pdf or .eps file), a visualization of the network (either .pdf or .eps file), and a file containing network statistics (.txt file). When nOutput = TRUE the function also returns an object of class list:

Note

We consider this to be a preliminary version of an envisioned wrapper than will take better form with subsequent versions of rags2ridges.

Author(s)

Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

See Also

ridgeP, conditionNumberPlot, optPenalty.LOOCVauto, sparsify, Ugraph, GGMnetworkStats

Examples

## Not run: 
## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]

## Employ the wrapper function
theWorks <- fullMontyS(X, lambdaMin = .5, lambdaMax = 30)
## End(Not run)

Description

Function for testing the null hypothesis that all population precision matrices are equal and thus the necessity for the fusion penalty. Note, the test performed is conditional on the supplied penalties and targets.

Usage

fused.test(Ylist, Tlist, lambda, n.permutations = 100, verbose = FALSE, ...)

Arguments

Details

The function computes the observed score statistic $U_obs$ using the fused ridge estimator on the given data. Next, the score statistic is computed a number of times (given by n.permutations) under the null-hypothesis by effectively permuting the class labels of the data.

Value

Returns a list values containing the observed test statistic and the test statistic under the null distribution.

Author(s)

Anders Ellern Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel, N. van Wieringen

References

Bilgrau, A.E., Peeters, C.F.W., Eriksen, P.S., Boegsted, M., and van Wieringen, W.N. (2020). Targeted Fused Ridge Estimation of Inverse Covariance Matrices from Multiple High-Dimensional Data Classes. Journal of Machine Learning Research, 21(26): 1-52.

See Also

ridgeP.fused

Examples

ns <- c(10, 5, 23)
Ylist <- createS(ns, p = 15, topology = "banded", dataset = TRUE)

# Use the identity target matrix for each class
Tlist <- replicate(length(ns), diag(15), simplify = FALSE)

# Do the test
lm <- matrix(10, 3, 3)
diag(lm) <- 1
ft <- fused.test(Ylist, Tlist, lambda = lm,
                 n.permutations = 500)
print(ft)

# Summary spits out a bit more information
summary(ft)

# The returned object can alo be plotted via
hist(ft)
# or via the alias
plot(ft)

# Customization and parameters work a usual:
hist(ft, col = "steelblue", main = "Null distribution", add.extra = FALSE,
     xlab = "Score statistic", freq = FALSE)

Description

Function that serves as a precursor function to the block-independence test (see GGMblockTest). It generates an empirical distribution of the penalty parameter under the null hypothesis of block independence (in the regularized precision matrix).

Usage

GGMblockNullPenalty(
  Y,
  id,
  nPerm = 25,
  lambdaMin,
  lambdaMax,
  lambdaInit = (lambdaMin + lambdaMax)/2,
  target = default.target(covML(Y)),
  type = "Alt",
  ncpus = 1
)

Arguments

Details

This function can be viewed as a precursor to the function for the block-independence test (see GGMblockTest). The mentioned test evaluates the null hypothesis of block-independence against the alternative of block-dependence (presence of non-zero elements in the off-diagonal block) in the precision matrix using high-dimensional data. To accommodate the high-dimensionality the parameters of interest are estimated in a penalized manner (ridge-type penalization, see ridgeP). Penalization involves a degree of freedom (the penalty parameter) which needs to be fixed before testing. This function then generates an empirical distribution of this penalty parameter. Hereto the samples are permutated within block. The resulting permuted data sets represent the null hypothesis. To avoid the dependence on a single permutation, many permuted data sets are generated. For each permutation the optimal penalty parameter is determined by means of cross-validation (see optPenalty.LOOCVauto). The resulting optimal penalty parameters are returned. An estimate of the location (such as the median) is recommended for use in the block-independence test.

Value

A numeric vector, representing the distribution of the (LOOCV optimal) penalty parameter under the null hypothesis of block-independence.

Author(s)

Wessel N. van Wieringen, Carel F.W. Peeters <[email protected]>

See Also

ridgeP, optPenalty.LOOCVauto, default.target, GGMblockTest

Examples

## Obtain some (high-dimensional) data
p = 15
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:15] = letters[1:15]
id <- c(rep(0, 10), rep(1, 5))

## Generate null distribution of the penalty parameter
lambda0dist <- GGMblockNullPenalty(X, id, 5, 0.001, 10)

## Location of null distribution
lambdaNull <- median(lambda0dist)

Description

Function performing a test that evaluates the null hypothesis of block-independence against the alternative of block-dependence (presence of non-zero elements in the off-diagonal block) in the precision matrix using high-dimensional data. The mentioned test is a permutation-based test (see details).

Usage

GGMblockTest(
  Y,
  id,
  nPerm = 1000,
  lambda,
  target = default.target(covML(Y)),
  type = "Alt",
  lowCiThres = 0.1,
  ncpus = 1,
  verbose = TRUE
)

Arguments

Details

The function performs a permutation test for the null hypothesis of block-independence against the alternative of block-dependence (presence of non-zero elements in the off-diagonal block) in the precision matrix using high-dimensional data. In the low-dimensional setting the common test statistic under multivariate normality (cf. Anderson, 2003) is:

$\log( \| \hat{\mathbf{\Sigma}}_a \| ) + \log( \| \hat{\mathbf{\Sigma}}_b \| ) - \log( \| \hat{\mathbf{\Sigma}} \| ),$

where the $\hat{\mathbf{\Sigma}}_a$, $\hat{\mathbf{\Sigma}}_b$, $\hat{\mathbf{\Sigma}}$ are the estimates of the covariance matrix in the sub- and whole group(s), respectively.

To accommodate the high-dimensionality the parameters of interest are estimated in a penalized manner (ridge-type penalization, see ridgeP). Penalization involves a degree of freedom (the penalty parameter: lambda) which needs to be fixed before testing. To decide on the penalty parameter for testing we refer to the GGMblockNullPenalty function. With an informed choice on the penalty parameter at hand, the null hypothesis is evaluated by permutation. Hereto the samples are permutated within block. The resulting permuted data set represents the null hypothesis. Many permuted data sets are generated. For each permutation the test statistic is calculated. The observed test statistic is compared to the null distribution from the permutations.

The function implements an efficient permutation resampling algorithm (see van Wieringen et al., 2008, for details.): If the probability of a p-value being below lowCiThres is smaller than 0.001 (read: the test is unlikely to become significant), the permutation analysis is terminated and a p-value of unity (1) is reported.

When verbose = TRUE also graphical output is generated: A histogram of the null-distribution. Note that, when ncpus is larger than 1, functionalities from snowfall are imported.

Value

An object of class list:

Author(s)

Wessel N. van Wieringen, Carel F.W. Peeters <[email protected]>

References

Anderson, T.W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd Edition. John Wiley.

van Wieringen, W.N., van de Wiel, M.A., and van der Vaart, A.W. (2008). A Test for Partial Differential Expression. Journal of the American Statistical Association 103: 1039-1049.

See Also

ridgeP, optPenalty.LOOCVauto, default.target, GGMblockNullPenalty

Examples

## Obtain some (high-dimensional) data
p = 15
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:15] = letters[1:15]
id <- c(rep(0, 10), rep(1, 5))

## Generate null distribution of the penalty parameter
lambda0dist <- GGMblockNullPenalty(X, id, 5, 0.001, 10)

## Location of null distribution
lambdaNull <- median(lambda0dist)

## Perform test
testRes <- GGMblockTest(X, id, nPerm = 100, lambdaNull)

Description

Function computing the mutual information between two exhaustive and mutually exclusive splits of a set of multivariate normal random variables.

Usage

GGMmutualInfo(S, split1)

Arguments

Value

A numeric, the mutual information between the variates forming split1 and those forming its complement.

Author(s)

Wessel N. van Wieringen, Carel F.W. Peeters <[email protected]>

References

Cover, T.M., Thomas, J.A. (2012), Elements of information theory.

See Also

covML, ridgeP.

Examples

# create a covariance matrix
Sigma <- covML(matrix(rnorm(100), ncol=5))

# impulse response analysis
GGMmutualInfo(Sigma, c(1,2))

Description

Function that calculates various network statistics from a sparse precision matrix. The sparse precision matrix is taken to represent the conditional indepence graph of a Gaussian graphical model.

Usage

GGMnetworkStats(sparseP, as.table = FALSE)

Arguments

Details

The function calculates various network statistics from a sparse matrix. The input matrix P is assumed to be a sparse precision or partial correlation matrix. The sparse matrix is taken to represent a conditional independence graph. In the Gaussian setting, conditional independence corresponds to zero entries in the (standardized) precision matrix. Each node in the graph represents a Gaussian variable, and each undirected edge represents conditional dependence in the sense of a nonzero corresponding precision entry.

The function calculates various measures of centrality: node degree, betweenness centrality, closeness centrality, and eigenvalue centrality. It also calculates the number of positive and the number of negative edges for each node. In addition, for each variate the mutual information (with all other variates), the variance, and the partial variance is represented. It is also indicated if the graph is chordal (i.e., triangulated). For more information on network measures, consult, e.g., Newman (2010).

Value

An object of class list when as.table = FALSE:

When as.table = TRUE the list items above (with the exception of chordal) are represented in tabular form as an object of class matrix.

Author(s)

Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

References

Newman, M.E.J. (2010). "Networks: an introduction", Oxford University Press.

See Also

ridgeP, covML, sparsify, Ugraph

Examples

## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]
Cx <- covML(X)

## Obtain sparsified partial correlation matrix
Pridge   <- ridgeP(Cx, 10, type = "Alt")
PCsparse <- sparsify(Pridge , threshold = "top")$sparseParCor

## Represent the graph and calculate GGM network statistics
Ugraph(PCsparse, "fancy")
## Not run: GGMnetworkStats(PCsparse)

Description

Compute various network statistics from a list sparse precision matrices. The sparse precision matrix is taken to represent the conditional independence graph of a Gaussian graphical model. This function is a simple wrapper for GGMnetworkStats.

Usage

GGMnetworkStats.fused(Plist)

Arguments

Details

For details on the columns see GGMnetworkStats.

Value

A data.frame of the various network statistics for each class. The names of Plist is prefixed to column-names.

Author(s)

Anders E. Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

See Also

GGMnetworkStats

Examples

## Create some "high-dimensional" data
set.seed(1)
p <- 10
ns <- c(5, 6)
Slist <- createS(ns, p)

## Obtain sparsified partial correlation matrix
Plist    <- ridgeP.fused(Slist, ns, lambda = c(5.2, 1.3), verbose = FALSE)
PCsparse <- sparsify.fused(Plist , threshold = "absValue", absValueCut = 0.2)
SPlist <- lapply(PCsparse, "[[", "sparsePrecision") # Get sparse precisions

## Calculate GGM network statistics in each class
## Not run: GGMnetworkStats.fused(SPlist)

Description

Function that calculates, for a specified node pair representing endpoints, path statistics from a sparse precision matrix. The sparse precision matrix is taken to represent the conditional independence graph of a Gaussian graphical model. The contribution to the observed covariance between the specified endpoints is calculated for each (heuristically) determined path between the endpoints.

Usage

GGMpathStats(
  P0,
  node1,
  node2,
  neiExpansions = 2,
  verbose = TRUE,
  graph = TRUE,
  nrPaths = 2,
  lay = "layout_in_circle",
  coords = NULL,
  nodecol = "skyblue",
  Vsize = 15,
  Vcex = 0.6,
  VBcolor = "darkblue",
  VLcolor = "black",
  all.edges = TRUE,
  prune = TRUE,
  legend = TRUE,
  scale = 1,
  Lcex = 0.8,
  PTcex = 2,
  main = ""
)

Arguments

Details

The conditional independence graph (as implied by the sparse precision matrix) is undirected. In undirected graphs origin and destination are interchangeable and are both referred to as 'endpoints' of a path. The function searches for shortest paths between the specified endpoints node1 and node2. It searches for shortest paths that visit nodes only once. The shortest paths between the provided endpoints are determined heuristically by the following procedure. The search is initiated by application of the get.all.shortest.paths-function from the igraph-package, which yields all shortest paths between the nodes. Next, the neighborhoods of the endpoints are defined (excluding the endpoints themselves). Then, the shortest paths are found between: (a) node1 and node Vs in its neighborhood; (b) node Vs in the node1-neighborhood and node Ve in the node2-neighborhood; and (c) node Ve in the node2-neighborhood and node2. These paths are glued and new shortest path candidates are obtained (preserving only novel paths). In additional iterations (specified by neiExpansions) the node1- and node2-neighborhood are expanded by including their neighbors (still excluding the endpoints) and shortest paths are again searched as described above.

The contribution of a particular path to the observed covariance between the specified node pair is calculated in accordance with Theorem 1 of Jones and West (2005). As in Jones and West (2005), paths whose weights have an opposite sign to the marginal covariance (between endnodes of the path) are referred to as 'moderating paths' while paths whose weights have the same sign as the marginal covariance are referred to as 'mediating' paths. Such paths are visualized when graph = TRUE.

All arguments following the graph argument are only (potentially) used when graph = TRUE. When graph = TRUE the conditional independence graph is returned with the paths highlighted that have the highest contribution to the marginal covariance between the specified endpoints. The number of paths highlighted is indicated by nrPaths. The edges of mediating paths are represented in green while the edges of moderating paths are represented in red. When all.edges = TRUE the edges other than those implied by the nrPaths-paths between node1 and node2 are also visualized (in lightgrey). When all.edges = FALSE only the mediating and moderating paths implied by nrPaths are visualized.

The default layout gives a circular placement of the vertices. Most layout functions supported by igraph are supported (the function is partly a wrapper around certain igraph functions). The igraph layouts can be invoked by a character that mimicks a call to a igraph layout functions in the lay argument. When using lay = NULL one can specify the placement of vertices with the coords argument. The row dimension of this matrix should equal the number of (pruned) vertices. The column dimension then should equal 2 (for 2D layouts) or 3 (for 3D layouts). The coords argument can also be viewed as a convenience argument as it enables one, e.g., to layout a graph according to the coordinates of a previous call to Ugraph. If both the the lay and the coords arguments are not NULL, the lay argument takes precedence

The arguments Lcex and PTcex are only used when legend = TRUE. If prune = TRUE the vertices of degree 0 (vertices not implicated by any edge) are removed. For the colors supported by the arguments nodecol, Vcolor, and VBcolor, see https://www.nceas.ucsb.edu/sites/default/files/2020-04/colorPaletteCheatsheet.pdf.

Value

An object of class list:

Note

Eppstein (1998) describes a more sophisticated algorithm for finding the top k shortest paths in a graph.

Author(s)

Wessel N. van Wieringen, Carel F.W. Peeters <[email protected]>

References

Eppstein, D. (1998). Finding the k Shortest Paths. SIAM Journal on computing 28: 652-673.

Jones, B., and West, M. (2005). Covariance Decomposition in Undirected Gaussian Graphical Models. Biometrika 92: 779-786.

See Also

ridgeP, optPenalty.LOOCVauto, sparsify

Examples

## Obtain some (high-dimensional) data
p <- 25
n <- 10
set.seed(333)
X <- matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X) <- letters[1:p]

## Obtain regularized precision under optimal penalty
OPT <- optPenalty.LOOCVauto(X, lambdaMin = .5, lambdaMax = 30)

## Determine support regularized standardized precision under optimal penalty
PC0 <- sparsify(OPT$optPrec, threshold = "localFDR")$sparseParCor

## Obtain information on mediating and moderating paths between nodes 14 and 23
pathStats <- GGMpathStats(PC0, 14, 23, verbose = TRUE, prune = FALSE)
pathStats

Description

A simple wrapper for GGMpathStats.

Usage

GGMpathStats.fused(sparsePlist, ...)

Arguments

Value

A list of path stats.

Note

The function currently fails if no paths are present in one of the groups.

Author(s)

Anders E. Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

See Also

GGMpathStats

Examples

## Obtain some (high-dimensional) data
set.seed(1)
ns <- c(10, 11)
Slist <- createS(ns, p = 7, topology = "banded")
Tlist <- default.target.fused(Slist, ns)

## Obtain regularized precision and sparsify
Plist <- ridgeP.fused(Slist, ns, Tlist, lambda = c(1, 1.6))
sparsePlist <- sparsify.fused(Plist, threshold = "absValue", absValueCut = 0.20)
SPlist <- lapply(sparsePlist, "[[", "sparsePrecision")

## Obtain information on mediating and moderating paths between nodes 14 and 23
res <- GGMpathStats.fused(SPlist, node1 = 3, node2 = 4, graph = FALSE)

Description

Plot a histogram of the null distribution and the observed test statistic in a permutation type "fusion test".

Usage

## S3 method for class 'ptest'
hist(x, add.extra = TRUE, ...)

## S3 method for class 'ptest'
plot(x, add.extra = TRUE, ...)

Arguments

Details

plot.ptest is simply a wrapper for hist.ptest.

Value

Invisibly returns x with extra additions.

Author(s)

Anders Ellern Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

References

Bilgrau, A.E., Peeters, C.F.W., Eriksen, P.S., Boegsted, M., and van Wieringen, W.N. (2020). Targeted Fused Ridge Estimation of Inverse Covariance Matrices from Multiple High-Dimensional Data Classes. Journal of Machine Learning Research, 21(26): 1-52.

See Also

fused.test, print.ptest

Examples

ns <- c(10, 5, 23)
Ylist <- createS(ns, p = 15, topology = "banded", dataset = TRUE)

# Use the identity target matrix for each class
Tlist <- replicate(length(ns), diag(15), simplify = FALSE)

# Do the test
lam <- matrix(10, 3, 3)
diag(lam) <- 1
ft <- fused.test(Ylist, Tlist, lambda = lam, n.permutations = 500)

# The returned object can alo be plotted via
hist(ft)
# or via the alias
plot(ft)

Description

Function to check if the argument submits to the various list-formats used by the fused ridge estimator and related functions are correct. That is, it tests if generic fused list arguments (such as Slist, Tlist, Plist, Ylist) are properly formatted.

Usage

is.Xlist(Xlist, Ylist = FALSE, semi = FALSE)

Arguments

Value

Returns TRUE if all tests are passed, throws error if not.

Author(s)

Anders Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

References

Bilgrau, A.E., Peeters, C.F.W., Eriksen, P.S., Boegsted, M., and van Wieringen, W.N. (2020). Targeted Fused Ridge Estimation of Inverse Covariance Matrices from Multiple High-Dimensional Data Classes. Journal of Machine Learning Research, 21(26): 1-52.

van Wieringen, W.N. & Peeters, C.F.W. (2016). Ridge Estimation of Inverse Covariance Matrices from High-Dimensional Data, Computational Statistics & Data Analysis, vol. 103: 284-303. Also available as arXiv:1403.0904v3 [stat.ME].

See Also

ridgeP.fused, optPenalty.fused

Examples

Slist <- createS(n = c(4, 6, 9), p = 10)
is.Xlist(Slist, semi = TRUE)

Description

Function to test if a matrix is symmetric positive (semi)definite or not.

Usage

isSymmetricPD(M)

isSymmetricPSD(M, tol = 1e-04)

Arguments

Details

Tests positive definiteness by Cholesky decomposition. Tests positive semi-definiteness by checking if all eigenvalues are larger than $-\epsilon|\lambda_1|$ where $\epsilon$ is the tolerance and $\lambda_1$ is the largest eigenvalue.

While isSymmetricPSD returns TRUE if the matrix is symmetric positive definite and FASLE if not. In practice, it tests if all eigenvalues are larger than -tol*|l| where l is the largest eigenvalue. More here.

Value

Returns a logical value. Returns TRUE if the M is symmetric positive (semi)definite and FALSE if not. If M is not even symmetric, the function throws an error.

Author(s)

Anders Ellern Bilgrau Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

See Also

isSymmetric

Examples

A <- matrix(rnorm(25), 5, 5)
## Not run: 
isSymmetricPD(A)

## End(Not run)
B <- symm(A)
isSymmetricPD(B)

C <- crossprod(B)
isSymmetricPD(C)

isSymmetricPSD(C)

Description

Construct a target matrix by combining topology information from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database and pilot data.

Usage

kegg.target(
  Y,
  kegg.id,
  method = "linreg",
  organism = "hsa",
  graph = getKEGGPathway(kegg.id)$graph
)

Arguments

Details

The function estimates the precision matrix based on the topology given by the KEGG database. Requires a connection to the internet.

Value

Returns a target matrix with size depending on the kegg.id.

Note

It is currently nessesary to require("KEGGgraph") (or require("KEGGgraph")) due to a bug in KEGGgraph.

Author(s)

Anders Ellern Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

References

https://www.genome.jp/kegg/

See Also

getKEGGPathway, default.target, and default.target.fused

Examples

## Not run: 
if (require("KEGGgraph")) {
kegg.g <- getKEGGPathway("map04115")$graph

# Create some toy data with the correct names
Y <- createS(n = 10, p = numNodes(kegg.g), dataset = TRUE)
colnames(Y) <- nodes(kegg.g)

T <- kegg.target(Y, "map04115")
print(T[1:10, 1:10])
}

## End(Not run)

Description

Function calculating the Kullback-Leibler divergence between two multivariate normal distributions.

Usage

KLdiv(Mtest, Mref, Stest, Sref, symmetric = FALSE)

Arguments

Details

The Kullback-Leibler (KL) information (Kullback and Leibler, 1951; also known as relative entropy) is a measure of divergence between two probability distributions. Typically, one distribution is taken to represent the ‘true’ distribution and functions as the reference distribution while the other is taken to be an approximation of the true distribution. The criterion then measures the loss of information in approximating the reference distribution. The KL divergence between two $p$-dimensional multivariate normal distributions $\mathcal{N}^{0}_{p}(\boldsymbol{\mu}_{0}, \mathbf{\Sigma}_{0})$ and $\mathcal{N}^{1}_{p}(\boldsymbol{\mu}_{1}, \mathbf{\Sigma}_{1})$ is given as

$\mathrm{I}_{KL}(\mathcal{N}^{0}_{p} \| \mathcal{N}^{1}_{p}) = \frac{1}{2}\left\{\mathrm{tr}(\mathbf{\Omega}_{1}\mathbf{\Sigma}_{0}) + (\boldsymbol{\mu}_{1} - \boldsymbol{\mu}_{0})^{\mathrm{T}} \mathbf{\Omega}_{1}(\boldsymbol{\mu}_{1} - \boldsymbol{\mu}_{0}) - p - \ln|\mathbf{\Sigma}_{0}| + \ln|\mathbf{\Sigma}_{1}| \right\},$

where $\mathbf{\Omega} = \mathbf{\Sigma}^{-1}$. The KL divergence is not a proper metric as $\mathrm{I}_{KL}(\mathcal{N}^{0}_{p} \| \mathcal{N}^{1}_{p}) \neq \mathrm{I}_{KL}(\mathcal{N}^{1}_{p} \| \mathcal{N}^{0}_{p})$. When symmetric = TRUE the function calculates the symmetric KL divergence (also referred to as Jeffreys information), given as

$\mathrm{I}_{KL}(\mathcal{N}^{0}_{p} \| \mathcal{N}^{1}_{p}) + \mathrm{I}_{KL}(\mathcal{N}^{1}_{p} \| \mathcal{N}^{0}_{p}).$

Value

Function returns a numeric representing the (symmetric) Kullback-Leibler divergence.

Author(s)

Wessel N. van Wieringen, Carel F.W. Peeters <[email protected]>

References

Kullback, S. and Leibler, R.A. (1951). On Information and Sufficiency. Annals of Mathematical Statistics 22: 79-86.

See Also

covML, ridgeP

Examples

## Define population
set.seed(333)
p = 25
n = 1000
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]
Cov0  <- covML(X)
mean0 <- colMeans(X)

## Obtain sample from population
samples <- X[sample(nrow(X), 10),]
Cov1  <- covML(samples)
mean1 <- colMeans(samples)

## Regularize singular Cov1
P <- ridgeP(Cov1, 10)
CovR <- solve(P)

## Obtain KL divergence
KLdiv(mean1, mean0, CovR, Cov0)

Description

Function calculating the Kullback-Leibler divergence between two sets of multivariate normal distributions. In other words, it calculates a weigthed mean of Kullback-Leibler divergences between multiple paired normal distributions.

Usage

KLdiv.fused(MtestList, MrefList, StestList, SrefList, ns, symmetric = FALSE)

Arguments

Value

Function returns a numeric representing the (optionally symmetric) fused Kullback-Leibler divergence.

Author(s)

Anders Ellern Bilgrau, Wessel N. van Wieringen, Carel F.W. Peeters <[email protected]>

See Also

KLdiv

Examples

# Create some toy data
n <- c(40, 60, 80)
p <- 10
Stest <- replicate(length(n), diag(p), simplify = FALSE)
Sref <- createS(n, p = p)

KLdiv.fused(StestList = Stest, SrefList = Sref, ns = n, symmetric = FALSE)
KLdiv.fused(StestList = Stest, SrefList = Sref, ns = n, symmetric = TRUE)

Description

Function that evaluates an estimated and possibly regularized precision matrix under various loss functions. The loss functions are formulated in precision terms. This function may be used to estimate the risk (vis-a-vis, say, the true precision matrix) of the various ridge estimators employed.

Usage

loss(E, T, precision = TRUE, type = c("frobenius", "quadratic"))

Arguments

Details

Let $\mathbf{\Omega}$ denote a generic $(p \times p)$ population precision matrix and let $\hat{\mathbf{\Omega}}(\lambda)$ denote a generic ridge estimator of the precision matrix under generic regularization parameter $\lambda$ (see also ridgeP). The function then considers the following loss functions:

1. Squared Frobenius loss, given by:

   $L_{F}[\hat{\mathbf{\Omega}}(\lambda), \mathbf{\Omega}] = \|\hat{\mathbf{\Omega}}(\lambda) - \mathbf{\Omega}\|_{F}^{2};$

2. Quadratic loss, given by:

   $L_{Q}[\hat{\mathbf{\Omega}}(\lambda), \mathbf{\Omega}] = \|\hat{\mathbf{\Omega}}(\lambda) \mathbf{\Omega}^{-1} - \mathbf{I}_{p}\|_{F}^{2}.$

The argument T is considered to be the true precision matrix when precision = TRUE. If precision = FALSE the argument T is considered to represent the true covariance matrix. This statement is needed so that the loss is properly evaluated over the precision, i.e., depending on the value of the logical argument precision inversions are employed where needed.

The function can be employed to assess the risk of a certain ridge precision estimator (see also ridgeP). The risk $\mathcal{R}_{f}$ of the estimator $\hat{\mathbf{\Omega}}(\lambda)$ given a loss function $L_{f}$, with $f \in \{F, Q\}$ can be defined as the expected loss:

$\mathcal{R}_{f}[\hat{\mathbf{\Omega}}(\lambda)] = \mathrm{E}\{L_{f}[\hat{\mathbf{\Omega}}(\lambda), \mathbf{\Omega}]\},$

which can be approximated by the mean or median of losses over repeated simulation runs.

Value

Function returns a numeric representing the loss under the chosen loss function.

Author(s)

Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

References

van Wieringen, W.N. & Peeters, C.F.W. (2016). Ridge Estimation of Inverse Covariance Matrices from High-Dimensional Data, Computational Statistics & Data Analysis, vol. 103: 284-303. Also available as arXiv:1403.0904v3 [stat.ME].

See Also

covML, ridgeP

Examples

## Define population covariance
set.seed(333)
p = 25
n = 1000
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]
Truecov <- covML(X)

## Obtain sample
samples <- X[sample(nrow(X), 10), ]
Cxx <- covML(samples)

## Obtain regularized precision
P <- ridgeP(Cxx, 10, type = "Alt")

## Evaluate estimated precision against population
## precision under Frobenius loss
loss(P, Truecov, precision = FALSE, type = "frobenius")

Description

Calculates the moments of the sample covariance matrix. It assumes that the summands (the outer products of the samples' random data vector) that constitute the sample covariance matrix follow a Wishart-distribution with scale parameter $\mathbf{\Sigma}$ and shape parameter $\nu$. The latter is equal to the number of summands in the sample covariance estimate.

Usage

momentS(Sigma, shape, moment = 1)

Arguments

Value

The $r$-th moment of a sample covariance matrix: $E(\mathbf{S}^r)$.

Author(s)

Wessel N. van Wieringen.

References

Lesac, G., Massam, H. (2004), "All invariant moments of the Wishart distribution", Scandinavian Journal of Statistics, 31(2), 295-318.

Examples

# create scale parameter
Sigma <- matrix(c(1, 0.5, 0, 0.5, 1, 0, 0, 0, 1), byrow=TRUE, ncol=3)

# evaluate expectation of the square of a sample covariance matrix
# that is assumed to Wishart-distributed random variable with the
# above scale parameter Sigma and shape parameter equal to 40.
momentS(Sigma, 40, 2)

Description

Functions that evaluate the (penalized) (fused) likelihood.

Usage

NLL(S, P)

PNLL(S, P, T, lambda)

NLL.fused(Slist, Plist, ns)

PNLL.fused(Slist, Plist, ns, Tlist, lambda)

Arguments

Value

A single number.

Author(s)

Anders Ellern Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

See Also

ridgeP, ridgeP.fused

Examples

ns <- c(4,5)
Slist <- createS(n = ns, p = 5)
Plist <- list(diag(5), diag(2,5))
Tlist <- list(diag(5), diag(5))

NLL(Slist[[1]], Plist[[1]])
PNLL(Slist[[1]], Plist[[1]], Tlist[[1]], lambda = 1)
NLL.fused(Slist, Plist, ns)
PNLL.fused(Slist, Plist, ns, Tlist, lambda = diag(2))

Description

Function that selects the optimal penalty parameter for the ridgeP call by usage of approximate leave-one-out cross-validation. Its output includes (a.o.) the precision matrix under the optimal value of the penalty parameter.

Usage

optPenalty.aLOOCV(
  Y,
  lambdaMin,
  lambdaMax,
  step,
  type = "Alt",
  cor = FALSE,
  target = default.target(covML(Y)),
  output = "light",
  graph = TRUE,
  verbose = TRUE
)

Arguments

Details

The function calculates an approximate leave-one-out cross-validated (aLOOCV) negative log-likelihood score (using a regularized ridge estimator for the precision matrix) for each value of the penalty parameter contained in the search grid. The utilized aLOOCV score was proposed by Lian (2011) and Vujacic et al. (2014). The aLOOCV negative log-likeliho od score is computationally more efficient than its non-approximate counterpart (see optPenalty.LOOCV). For details on the aLOOCV negative log-likelihood score see Lian (2011) and Vujacic et al (2014). For scalar matrix targets (see default.target) the complete solution path of the alternative Type I and II ridge estimators (see ridgeP) depends on only 1 eigendecomposition and 1 matrix inversion, making the determination of the optimal penalty value particularly efficient (see van Wieringen and Peeters, 2015).

The value of the penalty parameter that achieves the lowest aLOOCV negative log-likelihood score is deemed optimal. The penalty parameter must be positive such that lambdaMin must be a positive scalar. The maximum allowable value of lambdaMax depends on the type of ridge estimator employed. For details on the type of ridge estimator one may use (one of: "Alt", "ArchI", "ArchII") see ridgeP. The ouput consists of an object of class list (see below). When output = "light" (default) only the optLambda and optPrec elements of the list are given.

Value

An object of class list:

Note

When cor = TRUE correlation matrices are used in the computation of the approximate (cross-validated) negative log-likelihood score, i.e., the sample covariance matrix is a matrix on the correlation scale. When performing evaluation on the correlation scale the data are assumed to be standardized. If cor = TRUE and one wishes to used the default target specification one may consider using target = default.target(covML(Y, cor = TRUE)). This gives a default target under the assumption of standardized data.

Author(s)

Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

References

Lian, H. (2011). Shrinkage tuning parameter selection in precision matrices estimation. Journal of Statistical Planning and Inference, 141: 2839-2848.

van Wieringen, W.N. & Peeters, C.F.W. (2016). Ridge Estimation of Inverse Covariance Matrices from High-Dimensional Data, Computational Statistics & Data Analysis, vol. 103: 284-303. Also available as arXiv:1403.0904v3 [stat.ME].

Vujacic, I., Abbruzzo, A., and Wit, E.C. (2014). A computationally fast alternative to cross-validation in penalized Gaussian graphical models. arXiv: 1309.6216v2 [stat.ME].

See Also

ridgeP, optPenalty.LOOCV, optPenalty.LOOCVauto,
default.target, covML

Examples

## Obtain some (high-dimensional) data
p = 25
n = 10
set.seed(333)
X = matrix(rnorm(n*p), nrow = n, ncol = p)
colnames(X)[1:25] = letters[1:25]

## Obtain regularized precision under optimal penalty
OPT  <- optPenalty.aLOOCV(X, lambdaMin = .001, lambdaMax = 30, step = 400); OPT
OPT$optLambda	# Optimal penalty
OPT$optPrec	  # Regularized precision under optimal penalty

## Another example with standardized data
X <- scale(X, center = TRUE, scale = TRUE)
OPT  <- optPenalty.aLOOCV(X, lambdaMin = .001, lambdaMax = 30,
                          step = 400, cor = TRUE,
                          target = default.target(covML(X, cor = TRUE))); OPT
OPT$optLambda	# Optimal penalty
OPT$optPrec	  # Regularized precision under optimal penalty

Description

Functions to find the optimal ridge and fusion penalty parameters via leave-one-out cross validation. The functions support leave-one-out cross-validation (LOOCV), $k$-fold CV, and two forms of approximate LOOCV. Depending on the used function, general numerical optimization or a grid-based search is used.

Usage

optPenalty.fused.grid(
  Ylist,
  Tlist,
  lambdas = 10^seq(-5, 5, length.out = 15),
  lambdaFs = lambdas,
  cv.method = c("LOOCV", "aLOOCV", "sLOOCV", "kCV"),
  k = 10,
  verbose = TRUE,
  ...
)

optPenalty.fused.auto(
  Ylist,
  Tlist,
  lambda,
  cv.method = c("LOOCV", "aLOOCV", "sLOOCV", "kCV"),
  k = 10,
  verbose = TRUE,
  lambda.init,
  maxit.ridgeP.fused = 1000,
  optimizer = "optim",
  maxit.optimizer = 1000,
  debug = FALSE,
  optim.control = list(trace = verbose, maxit = maxit.optimizer),
  ...
)

optPenalty.fused(
  Ylist,
  Tlist,
  lambda = default.penalty(Ylist),
  cv.method = c("LOOCV", "aLOOCV", "sLOOCV", "kCV"),
  k = 10,
  grid = FALSE,
  ...
)

Arguments

Details

optPenalty.fused.auto serves a utilizes optim for identifying the optimal fused parameters and works for general classes of penalty graphs.

optPenalty.fused.grid gives a grid-based evaluation of the (approximate) LOOCV loss.

Value

optPenalty.fused.auto returns a list:

optPenalty.fused.LOOCV returns a list:

Author(s)

Anders Ellern Bilgrau, Carel F.W. Peeters <[email protected]>, Wessel N. van Wieringen

References

Bilgrau, A.E., Peeters, C.F.W., Eriksen, P.S., Boegsted, M., and van Wieringen, W.N. (2020). Targeted Fused Ridge Estimation of Inverse Covariance Matrices from Multiple High-Dimensional Data Classes. Journal of Machine Learning Research, 21(26): 1-52.

See Also

See also default.penalty, optPenalty.LOOCV.

Examples

## Not run: 
# Generate some (not so) high-dimensional data witn (not so) many samples
ns <- c(4, 5, 6)
Ylist <- createS(n = ns, p = 6, dataset = TRUE)
Slist <- lapply(Ylist, covML)
Tlist <- default.target.fused(Slist, ns, type = "DIAES")


# Grid-based
lambdas <- 10^seq(-5, 3, length.out = 7)
a <- optPenalty.fused.grid(Ylist, Tlist,
                           lambdas = lambdas,
                           cv.method = "LOOCV", maxit = 1000)
b <- optPenalty.fused.grid(Ylist, Tlist,
                           lambdas = lambdas,
                           cv.method = "aLOOCV", maxit = 1000)
c <- optPenalty.fused.grid(Ylist, Tlist,
                           lambdas = lambdas,
                           cv.method = "sLOOCV", maxit = 1000)
d <- optPenalty.fused.grid(Ylist, Tlist,
                           lambdas = lambdas,
                           cv.method = "kCV", k = 2, maxit = 1000)

# Numerical optimization (uses the default "optim" optimizer with method "BFGS")
aa <- optPenalty.fused.auto(Ylist, Tlist, cv.method = "LOOCV", method = "BFGS")
print(aa)
bb <- optPenalty.fused.auto(Ylist, Tlist, cv.method = "aLOOCV", method = "BFGS")
print(bb)
cc <- optPenalty.fused.auto(Ylist, Tlist, cv.method = "sLOOCV", method = "BFGS")
print(cc)
dd <- optPenalty.fused.auto(Ylist, Tlist, cv.method = "kCV", k=3, method="BFGS")
print(dd)


#
# Plot the results
#

# LOOCV
# Get minimums and plot
amin  <- log(expand.grid(a$lambda, a$lambdaF))[which.min(a$fcvl), ]
aamin <- c(log(aa$lambda[1,1]), log(aa$lambda[1,2]))

# Plot
filled.contour(log(a$lambda), log(a$lambdaF), log(a$fcvl), color = heat.colors,
               plot.axes = {points(amin[1], amin[2], pch = 16);
                            points(aamin[1], aamin[2], pch = 16, col = "purple");
                            axis(1); axis(2)},
               xlab = "lambda", ylab = "lambdaF", main = "LOOCV")

# Approximate LOOCV
# Get minimums and plot
bmin <- log(expand.grid(b$lambda, b$lambdaF))[which.min(b$fcvl), ]
bbmin <- c(log(bb$lambda[1,1]), log(unique(bb$lambda[1,2])))

filled.contour(log(b$lambda), log(b$lambdaF), log(b$fcvl), color = heat.colors,
               plot.axes = {points(bmin[1], bmin[2], pch = 16);
                            points(bbmin[1], bbmin[2], pch = 16, col ="purple");
                            axis(1); axis(2)},
               xlab = "lambda", ylab = "lambdaF", main = "Approximate LOOCV")


#
# Arbitrary penalty graphs
#

# Generate some new high-dimensional data and a 2 by 2 factorial design
ns <- c(6, 5, 3, 2)
df <- expand.grid(Factor1 = LETTERS[1:2], Factor2 = letters[3:4])
Ylist <- createS(n = ns, p = 4, dataset = TRUE)
Tlist <- lapply(lapply(Ylist, covML), default.target, type = "Null")

# Construct penalty matrix
lambda <- default.penalty(df, type = "CartesianUnequal")

# Find optimal parameters,
# Using optim with method "Nelder-Mead" with "special" LOOCV
ans1 <- optPenalty.fused(Ylist, Tlist, lambda = lambda,
                         cv.method = "sLOOCV", verbose = FALSE)
print(ans1$lambda.unique)

# By approximate LOOCV using optim with method "BFGS"
ans2 <- optPenalty.fused(Ylist, Tlist, lambda = lambda,
                         cv.method = "aLOOCV", verbose = FALSE,
                         method = "BFGS")
print(ans2$lambda.unique)

# By LOOCV using nlm
lambda.init <- matrix(1, 4, 4)
lambda.init[cbind(1:4,4:1)] <- 0
ans3 <- optPenalty.fused(Ylist, Tlist, lambda = lambda,
                         lambda.init = lambda.init,
                         cv.method = "LOOCV", verbose = FALSE,
                         optimizer = "nlm")
print(ans3$lambda.unique)

# Quite different results!


#
# Arbitrary penalty graphs with fixed penalties!
#

# Generate some new high-dimensional data and a 2 by 2 factorial design
ns <- c(6, 5, 5, 5)
df <- expand.grid(DS = LETTERS[1:2], ER = letters[3:4])
Ylist <- createS(n = ns, p = 4, dataset = TRUE)
Tlist <- lapply(lapply(Ylist, covML), default.target, type = "Null")

lambda <- default.penalty(df, type = "Tensor")
print(lambda)  # Say we want to penalize the pair (1,2) with strength 2.1;
lambda[2,1] <- lambda[1,2] <- 2.1
print(lambda)

# Specifiying starting values is also possible:
init <- diag(length(ns))
init[2,1] <- init[1,2] <- 2.1

res <- optPenalty.fused(Ylist, Tlist, lambda = lambda, lambda.init = init,
                        cv.method = "aLOOCV", optimizer = "nlm")
print(res)

## End(Not run)