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


# 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

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

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


Zlist <- Array2List(Z)
conlist <- lapply(Zlist, function(x) apply(x, 2, function(a) pic(a, phy)))

relist.l <- lapply( conlist, function(x) {
	re.lorec <- lorec(cor(x), 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))
	Lhat + Shat
})

luo <- relist.l[[1]]


# Set up plotting parameters
mycols <- rainbow(100)
mycols <- heat.colors(100)
mycols <- tim.colors(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()