Source

sicb2013 / regularization.r

Full commit
source("functions.r")
require(picante)
require(ape) 
require(geiger)
require(lorec)
require(fields)

set.seed(123456)

# Build the true matrix
trueCovariance = test_cov_matrix()

# Run the phylogenetic analyses
tree_text <- "(((Species_A:10.0,Species_B:10.0)I:40.0,(Species_C:40.0,Species_D:40.0)J:10.0)K:10.0,(Species_E:45.0,(Species_F:30.0,(Species_G:15.0,Species_H:15.0)O:15.0)N:15.0)M:15.0)L:0.0;"
phy <- read.tree( text=tree_text )
Z <- sim.char(phy, trueCovariance, nsims = 1, model = "brownian", root.state = 1)
W <- Z[,,1]	# Grab one simulation
contrasts <- apply( W, 2, function(a) pic(a, phy) )


# The covariance matrix estimated directly from the contrasts

contrastcor <- cor.table(contrasts)$r

p <- ncol( contrasts )
n <- nrow( contrasts )

# Regularization by thresholding
# Bickel, P. J. & Levina, E. Covariance regularization by thresholding. 
# Ann. Statist. 36, 2577–2604 (2008). http://dx.doi.org/10.1214/08-AOS600

bickel <- contrastcor * ( abs(contrastcor) > sqrt(log(p)/n) )

# Regularization by convex Minimization
# Luo, X. High Dimensional Low Rank and Sparse Covariance Matrix Estimation via 
# Convex Minimization. arXiv.org (2011). http://arxiv.org/abs/1111.1133


re.lorec <- lorec( contrastcor, diag(1, 100), diag(1, 100), sqrt(p/n), 
	sqrt(log(p)/n) )
re.lorec.eig <- eigen( re.lorec$L )
# threholding both
sel <- re.lorec.eig$values > sqrt( p/n )
if ( sum(sel) > 0 ) {
	V <- re.lorec.eig$vectors
	V <- V * ( abs(V)>sqrt(1/p) )
	V <- V[ ,1:sum(sel) ]
	Lhat <- V%*%diag( re.lorec.eig$values[1:sum(sel)] )%*%t(V)
} else {
	Lhat <- matrix(0, p, p)
}
Shat <- re.lorec$S
Shat <- Shat*(abs(Shat)>sqrt(1/p))
luo <- Lhat + Shat


# Set up plotting parameters
# Color plots
# mycols <- tim.colors(100)
# Greyscale plots
mycols <- gray(1:100 / 100)

mybreaks <- seq(-1, 1, length.out=101)


# Plot the results

pdf(file="regularization.pdf", width=7, height=3.5)
layout(rbind(c(1,2,3,4),c(0,5,5,0)), heights=c(.75,.5))

m <- 1
par( mar=c(m,m,m,m) )
par( oma=c(0,0,0,0) )

plot_matrix( trueCovariance, main="(a) True", breaks=mybreaks, col=mycols )
plot_matrix( contrastcor, main="(b) Contrasts", breaks=mybreaks, col=mycols )
plot_matrix( luo, main="(c) Minimization", breaks=mybreaks, col=mycols )
plot_matrix( bickel, main="(d) Thresholding", breaks=mybreaks, col=mycols )

frame()
par( mar=c(0,0,0,0))
image.plot( legend.only=TRUE, zlim=range(-1:1), col=mycols, legend.width=3, 
	horizontal=TRUE)

dev.off()