Commits

Casey Dunn  committed e229d13

initialized repo

  • Participants

Comments (0)

Files changed (3)

+# sicb2013
+
+This repository includes the source code for the analyses presented 
+at the [EDEN](http://edenrcn.com) sponsored symposium ["Understanding First 
+Order Phenotypes: Transcriptomics for Emerging Model Systems"](http://www.sicb.org/meetings/2013/symposia/phenotypes.php) 
+at the 2013 SICB conference, described in 
+[abstract S4-2.1](http://www.sicb.org/meetings/2013/SICB%202013%20abstracts.pdf).
+
+This study is a collaborative project by Brown University faculty members:
+
+- [Casey Dunn](http://www.brown.edu/Faculty/Dunn_Lab/)
+
+- [Xi Luo](http://www.stat.brown.edu/FacultyDisplay.aspx?id=1317070681)
+
+- [Zhijin Wu](http://www.stat.brown.edu/FacultyDisplay.aspx?id=1128605312)
+
+All collaborators contributed to the code presented here.
+
+Additional detail will be provided in a forthcoming publication.
+
+Many thanks to [Joe Felsenstein](http://www.gs.washington.edu/faculty/felsenstein.htm) 
+for helpful discussions about multivariate comparative analyses.
+
+## License
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the [GNU General Public License as published by
+the Free Software Foundation](http://www.gnu.org/licenses/), either 
+version 3 of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+test_cov_matrix <- function() {
+	# Builds a simple character covariance matrix for simulating data
+	
+	G <- 100
+	
+	trueCovariance = matrix(0,G,G)
+	trueCovariance[1:10,1:10] = 0.95
+	trueCovariance[11:80,11:80] = 0.3
+	trueCovariance[81:100,81:100] = 0.7
+	diag(trueCovariance) = 1
+
+	return( trueCovariance ) 
+}
+
+
+Array2List <- function(A) {
+	s <- dim(A)
+	R <- vector("list", s[3])
+
+	for (j in 1:s[3]) {
+		R[[j]] <- A[,,j]
+	}
+	return(R)
+}
+
+
+plot_matrix <- function(m, ... ) {
+	
+	nr <- nrow(m)
+	nc <- ncol(m)
+	image(1:nc, 1:nr, t(m[nr:1, ]), axes=F,xlab="", ylab="", ... )
+}

File regularization.r

+source("functions.r")
+require(picante)
+require(ape) 
+require(geiger)
+require(lorec)
+require(fields)
+
+set.seed(123456)
+
+
+# Set up plotting parameters
+mycols <- rainbow(100)
+mycols <- heat.colors(100)
+mycols <- tim.colors(100)
+mybreaks <- seq(-1, 1, length.out=101)
+
+
+# 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) )
+
+p <- ncol( contrasts )
+n <- nrow( contrasts )
+
+# 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
+
+Zlist <- Array2List(Z)
+conlist <- lapply(Zlist, function(x) apply(x, 2, function(a) pic(a, phy)))
+
+relist.s <- lapply(conlist, function(x) {
+	cor(x)
+})
+
+bickel <- relist.s[[1]]*(abs(relist.s[[1]])>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
+
+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]]
+
+
+# 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()