corpcor (version 1.0.0)

cov.shrink: Shrinkage Estimates of Covariance and (Partial) Correlation

Description

The functions cov.shrink, cor.shrink, and pcor.shrink implement a shrinkage approach to estimate covariance and (partial) correlation matrices. The advantages of using this approach in comparison with the standard empirical estimates (cov and cor) are that 1) the shrinkage estimates are always positive definite, 2) well conditioned (so that the inverse always exists), and 3) exhibit (sometimes dramatically) better mean squared error. Furthermore, they are inexpensive to compute and do not need the specification of any parameters (the shrinkage intensity is analytically estimated from the data).

Usage

cov.shrink(x, lambda, verbose=TRUE)
cor.shrink(x, lambda, verbose=TRUE)
pcor.shrink(x, lambda, verbose=TRUE)

Arguments

x

a data matrix

lambda

the shrinkage intensity (range 0-1). It is is not specified (the default) lambda is chosen such that the resulting shrinkage estimate has minimal MSE (see below for details).

verbose

report progress while computing (default: TRUE)

Value

• cov.shrink returns a matrix with the shrinkage estimate of the covariance matrix. cor.shrink returns the corresponding standardized matrix. pcor.shrink returns the partical correlations computed from the output of cov.shrink.

Details

Let $S$ be the usual empirical unbiased estimator of the true covariance matrix $\Sigma$ with $E(S)=\Sigma$, and the target $T$ an unbiased estimator of a structured covariance matrix (with $E(T)=\Theta$. Then a shrinkage estimator $S^{*}$ can be constructed by setting $$S^{*} = \lambda T + (1-\lambda) S .$$ If the shrinkage intensity $\lambda$ is large, then the resulting estimate will be close the structured estimate $T$, otherwise the unstructured estimate $S$ dominates. Ledoit and Wolf (2003) suggest an analytic formula to estimate $\lambda$ from the observed data such that the mean squared error (MSE) of $S^{*}$ is minimal. In cov.shrink the target $T$ is the diagonal matrix $diag(S)$. The above equation then simplifies to a shrinkage estimator of the offdiagonal elements only, i.e. with $s^{*}_{ii} = s_{ii}$ and $s^{*}_{ij} = (1-\lambda) s_{ij}$ for $i \neq j$. It then can be shown (see Schaefer and Strimmer 2005) that an estimator for the optimal intensity is given by $$\lambda^{*} = \sum_{i \neq j} Var(s_{ij}) / \sum_{i \neq j} s_{ij}^2 .$$ This shrinkage estimator are especially useful when computing sparse large-scale covariance matrices in a 'small n, large p' setting. These situations are often encountered in problems bioinformatics and statistical genomics (see Schaefer and Strimmer 2005 for examples).

References

Ledoit, O., and Wolf. M. (2003). Improved estimation of the covariance matrix of stock returns with an application to portfolio selection. J. Emp. Finance 10:503-621.

Schaefer, J., and Strimmer, K. (2005). A shrinkage approach to large-scale covariance estimation and implications for functional genomics. Submitted to SAGMB.

See Also

cor2pcor

Examples

# load corpcor library
library("corpcor")

# small n, large p
p <- 100
n <- 20

# generate random pxp covariance matrix
sigma <- matrix(rnorm(p*p),ncol=p)
sigma <- crossprod(sigma)+ diag(rep(0.1, p))

# simulate multinormal data of sample size n  
sigsvd <- svd(sigma)
Y <- t(sigsvd$v %*% (t(sigsvd$u) * sqrt(sigsvd$d)))
X <- matrix(rnorm(n * ncol(sigma)), nrow = n) %*% Y


# estimate covariance matrix
s1 <- cov(X)
s2 <- cov.shrink(X)

# squared error
sum((s1-sigma)^2)
sum((s2-sigma)^2)

# varcov produces the same results as cov
vc <- varcov(X)
sum(abs(vc$S-s1))

# compare positive definiteness
is.positive.definite(s1)
is.positive.definite(s2)
is.positive.definite(sigma)

# compare ranks and condition
rank.condition(s1)
rank.condition(s2)
rank.condition(sigma)

# compare eigenvalues
e1 <- eigen(s1, symmetric=TRUE)$values
e2 <- eigen(s2, symmetric=TRUE)$values
e3 <- eigen(sigma, symmetric=TRUE)$values
m <-max(e1, e2, e3)
yl <- c(0, m)

par(mfrow=c(1,3))
plot(e1,  main="empirical")
plot(e2,  ylim=yl, main="shrinkage")
plot(e3,  ylim=yl, main="true")
par(mfrow=c(1,1))