covgpenal: covgpenal
| covgpenal | R Documentation |
covgpenal
Description
This function computes the empirical penalized Gaussian copula covariance matrix with the Gaussian log-likelihood plus a lasso-type penalty as objective function. Model selection is done by choosing omega such that BIC is maximal.
Usage
covgpenal(
S,
n,
omegas,
derpenal = function(t, omega) {
derSCAD(t, omega, 3.7)
},
nsteps = 1
)
Arguments
S |
The sample matrix of normal scores covariances. |
n |
The sample size. |
omegas |
The candidate values for the tuning parameter in |
derpenal |
The derivative of the penalty function to be used (default = scad with parameter |
nsteps |
The number of weighted covariance graphical lasso iterations (default = 1). |
Details
The aim is to solve/compute
\widehat{\boldsymbol{\Sigma}}_{\text{LT},n} \in \text{arg min}_{\boldsymbol{\Sigma} > 0} \left \{
\ln \left | \boldsymbol{\Sigma} \right | + \text{tr} \left (\boldsymbol{\Sigma}^{-1} \widehat{\boldsymbol{\Sigma}}_{n} \right ) + P_{\text{LT}}\left (\boldsymbol{\Sigma},\omega_{n} \right ) \right \},
where the penalty function P_{\text{LT}} is of lasso-type:
P_{\text{LT}} \left (\boldsymbol{\Sigma},\omega_{n} \right ) = \sum_{ij} p_{\omega_{n}} \left (\Delta_{ij} \left |\sigma_{ij} \right | \right ),
for a certain penalty function p_{\omega_{n}} with penalty parameter \omega_{n}, and \sigma_{ij} the
(i,j)'th entry of \boldsymbol{\Sigma} with \Delta_{ij} = 1 if i \neq j and \Delta_{ij} = 0 if i = j (in order to not shrink the variances).
The matrix \widehat{\boldsymbol{\Sigma}}_{n} is the matrix of sample normal scores covariances.
In case p_{\omega_{n}}(t) = \omega_{n} t is the lasso penalty, the implementation for the
(weighted) covariance graphical lasso is available in the R package ‘covglasso’ (see the manual for further explanations). For general penalty functions,
we perform a local linear approximation to the penalty function and iteratively do (nsteps, default = 1) weighted covariance graphical lasso optimizations.
The default for the penalty function is the scad (derpenal = derivative of scad penalty), i.e.,
p_{\omega_{n},\text{scad}}^{\prime}(t) = \omega_{n} \left [1 \left (t \leq \omega_{n} \right ) + \frac{\max \left (a \omega_{n} - t,0 \right )}{\omega_{n} (a-1)} 1 \left (t > \omega_{n} \right ) \right ],
with a = 3.7 by default.
For tuning \omega_{n}, we maximize (over a grid of candidate values) the BIC criterion
\text{BIC} \left (\widehat{\boldsymbol{\Sigma}}_{\omega_{n}} \right ) = -n \left [\ln \left |\widehat{\boldsymbol{\Sigma}}_{\omega_{n}} \right | + \text{tr} \left (\widehat{\boldsymbol{\Sigma}}_{\omega_{n}}^{-1} \widehat{\boldsymbol{\Sigma}}_{n} \right ) \right ] - \ln(n) \text{df} \left (\widehat{\boldsymbol{\Sigma}}_{\omega_{n}} \right ),
where \widehat{\boldsymbol{\Sigma}}_{\omega_{n}} is the estimated candidate covariance matrix using \omega_{n}
and df (degrees of freedom) equals the number of non-zero entries in \widehat{\boldsymbol{\Sigma}}_{\omega_{n}}, not taking the elements under the diagonal into account.
Value
A list with elements "est" containing the (lasso-type) penalized matrix of sample normal scores rank correlations (output as provided by the function “covglasso.R”), and "omega" containing the optimal tuning parameter.
References
De Keyser, S. & Gijbels, I. (2024). High-dimensional copula-based Wasserstein dependence. doi: https://doi.org/10.48550/arXiv.2404.07141.
Fop, M. (2021). covglasso: sparse covariance matrix estimation, R package version 1.0.3. url: https://CRAN.R-project.org/package=covglasso.
Wang, H. (2014). Coordinate descent algorithm for covariance graphical lasso. Statistics and Computing 24:521-529. doi: https://doi.org/10.1007/s11222-013-9385-5.
See Also
grouplasso for group lasso estimation of the normal scores rank correlation matrix.
Examples
q = 10
dim = c(5,5)
n = 100
# AR(1) correlation matrix with correlation 0.5
R = 0.5^(abs(matrix(1:q-1,nrow = q, ncol = q, byrow = TRUE) - (1:q-1)))
# Sparsity on off-diagonal blocks
R0 = createR0(R,dim)
# Sample from multivariate normal distribution
sample = mvtnorm::rmvnorm(n,rep(0,q),R0,method = "chol")
# Normal scores
scores = matrix(0,n,q)
for(j in 1:q){scores[,j] = qnorm((n/(n+1)) * ecdf(sample[,j])(sample[,j]))}
# Sample matrix of normal scores covariances
Sigma_est = cov(scores) * ((n-1)/n)
# Candidate tuning parameters
omega = seq(0.01, 0.6, length = 50)
Sigma_est_penal = covgpenal(Sigma_est,n,omega)
