Automatic search for penalty parameter of ridge precision estimator with known chordal support

Automic search for the optimal ridge penalty parameter for the ridge estimator of the precision matrix with known chordal support. Optimal in the sense that it yields the maximum cross-validated likelihood. The search employs the Brent algorithm as implemented in the optim function.

optPenaltyPchordal(
  Y,
  lambdaMin,
  lambdaMax,
  lambdaInit = (lambdaMin + lambdaMax)/2,
  zeros,
  cliques = list(),
  separators = list(),
  target = default.target(covML(Y)),
  type = "Alt"
)

Arguments

Y: Data matrix. Variables assumed to be represented by columns.
lambdaMin: A numeric giving the minimum value for the penalty parameter.
lambdaMax: A numeric giving the maximum value for the penalty parameter.
lambdaInit: A numeric giving the initial value for the penalty parameter.
zeros: A matrix with indices of entries of the precision matrix that are constrained to zero. The matrix comprises two columns, each row corresponding to an entry of the precision matrix. The first column contains the row indices and the second the column indices. The specified conditional independence graph implied by the zero-structure of the precision should be undirected and decomposable. If not, it is symmetrized and triangulated.
cliques: A list-object containing the node indices per clique as obtained from the support4ridgeP-function.
separators: A list-object containing the node indices per separator as obtained from the support4ridgeP-function.
target: A target matrix (in precision terms) for Type I ridge estimators.
type: A character indicating the type of ridge estimator to be used. Must be one of: Alt, ArchI, ArchII.

Value

A numeric with the LOOCV optimal choice for the ridge penalty parameter.

Details

See the function optim for details on the implementation of the Brent algorithm.

References

Miok, V., Wilting, S.M., Van Wieringen, W.N. (2016), "Ridge estimation of the VAR(1) model and its time series chain graph from multivariate time-course omics data", Biometrical Journal, 59(1), 172-191.

Van Wieringen, W.N. and Peeters, C.F.W. (2016), "Ridge Estimation of Inverse Covariance Matrices from High-Dimensional Data", Computational Statistics and Data Analysis, 103, 284-303.

Author

Wessel N. van Wieringen.

Examples


# generate data
p <- 8
n <- 100
set.seed(333)
Y <- matrix(rnorm(n*p), nrow = n, ncol = p)

# define zero structure
S <- covML(Y)
S[1:3, 6:8] <- 0
S[6:8, 1:3] <- 0
zeros <- which(S==0, arr.ind=TRUE)

# obtain (triangulated) support info
supportP <- support4ridgeP(nNodes=p, zeros=zeros)

# determine optimal penalty parameter
if (FALSE) {
optLambda <- optPenaltyPchordal(Y, 10^(-10), 10, 0.1, zeros=supportP$zeros,
  cliques=supportP$cliques, separators=supportP$separators)
}
optLambda <- 0.1

# estimate precision matrix with known (triangulated) support
Phat <- ridgePchordal(S, optLambda, zeros=supportP$zeros,
  cliques=supportP$cliques, separators=supportP$separators)
#>   
#> Progress report .... 
#> -> ---------------------------------------------------------------------- 
#> -> optimization over                   : 27 out of 36 unique 
#> ->                                       parameters (75%) 
#> -> cond. number of initial estimate    : 1.95 (if >> 100, consider 
#> ->                                       larger values of lambda) 
#> -> # graph components                  : 1 (optimization per component) 
#> -> optimization per component          :  
#> 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |......................................................................| 100%
#> 
#> -> estimation done ... 
#> -> formatting output ... 
#> -> overall summary .... 
#> -> initial pen. log-likelihood         : -7.45199157 
#> -> optimized pen. log-likelihood       : -7.45199081 
#> -> optimization                        : converged (most likely) 
#> ->                                       for all components 
#> -> ----------------------------------------------------------------------