# Commits

committed e064ae1

renamed files for clarity

• Participants
• Parent commits f7e2b2b
• Branches master

# File functions.r

`-#' regularization_by_thresholding`
`-#' `
`-#' Regularizes a correlation matrix by thresholding`
`-#' `
`-#' @param cor_matrix A correlation matrix`
`-#' `
`-#' @param n The number of independent observations that were used to construct `
`-#' the `
`-#' `
`-#' @return The regularized correlation matrix`
`-#' `
`-#' @references`
`-#' Bickel, P. J. & Levina, E. Covariance regularization by thresholding. `
`-#' Ann. Statist. 36, 2577–2604 (2008). http://dx.doi.org/10.1214/08-AOS600 `
`-#' `
`-regularization_by_thresholding <- function ( cor_matrix, n ) {`
`-	# 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( cor_matrix )`
`-	`
`-	regularized <- cor_matrix * ( abs(cor_matrix) > sqrt(log(p)/n) )`
`-	`
`-	return( regularized )`
`-}`
`-`
`-`
`-#' regularization_by_convex_minimization`
`-#' `
`-#' Regularizes a correlation matrix by convex minimization`
`-#' `
`-#' @param cor_matrix A correlation matrix`
`-#' `
`-#' @param n The number of independent observations that were used to construct `
`-#' the `
`-#' `
`-#' @return The regularized correlation matrix`
`-#' `
`-#' @references`
`-#' Luo, X. High Dimensional Low Rank and Sparse Covariance Matrix Estimation via `
`-#' Convex Minimization. arXiv.org (2011). http://arxiv.org/abs/1111.1133`
`-#' `
`-regularization_by_convex_minimization <- function ( cor_matrix, n ) {`
`-`
`-	p <- ncol( cor_matrix )`
`-	`
`-	re.lorec <- lorec( cor_matrix, 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))`
`-	regularized <- Lhat + Shat`
`-	return( regularized )`
`-}`

# File regularization.r

`-source("functions.r")`
`-require(picante)`
`-require(ape) `
`-require(geiger)`
`-require(lorec)`
`-require(fields)`
`-`
`-set.seed(123456)`
`-`
`-`
`-# A couple functions to facilitate simulation`
`-`
`-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 ) `
`-}`
`-`
`-plot_matrix <- function(m, ... ) {`
`-	# Basic plotting`
`-	nr <- nrow(m)`
`-	nc <- ncol(m)`
`-	image(1:nc, 1:nr, t(m[nr:1, ]), axes=F,xlab="", ylab="", ... )`
`-}`
`-`
`-`
`-# 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`
`-`
`-n <- nrow( contrasts )`
`-`
`-# Regularization by thresholding`
`-`
`-bickel <- regularization_by_thresholding( contrastcor, n )`
`-`
`-# Regularization by convex Minimization`
`-`
`-luo <- regularization_by_convex_minimization( contrastcor, n )`
`-`
`-`
`-# 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()`

# File regularization_functions.r

`+#' regularization_by_thresholding`
`+#' `
`+#' Regularizes a correlation matrix by thresholding`
`+#' `
`+#' @param cor_matrix A correlation matrix`
`+#' `
`+#' @param n The number of independent observations that were used to construct `
`+#' the `
`+#' `
`+#' @return The regularized correlation matrix`
`+#' `
`+#' @references`
`+#' Bickel, P. J. & Levina, E. Covariance regularization by thresholding. `
`+#' Ann. Statist. 36, 2577–2604 (2008). http://dx.doi.org/10.1214/08-AOS600 `
`+#' `
`+regularization_by_thresholding <- function ( cor_matrix, n ) {`
`+	# 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( cor_matrix )`
`+	`
`+	regularized <- cor_matrix * ( abs(cor_matrix) > sqrt(log(p)/n) )`
`+	`
`+	return( regularized )`
`+}`
`+`
`+`
`+#' regularization_by_convex_minimization`
`+#' `
`+#' Regularizes a correlation matrix by convex minimization`
`+#' `
`+#' @param cor_matrix A correlation matrix`
`+#' `
`+#' @param n The number of independent observations that were used to construct `
`+#' the `
`+#' `
`+#' @return The regularized correlation matrix`
`+#' `
`+#' @references`
`+#' Luo, X. High Dimensional Low Rank and Sparse Covariance Matrix Estimation via `
`+#' Convex Minimization. arXiv.org (2011). http://arxiv.org/abs/1111.1133`
`+#' `
`+regularization_by_convex_minimization <- function ( cor_matrix, n ) {`
`+`
`+	p <- ncol( cor_matrix )`
`+	`
`+	re.lorec <- lorec( cor_matrix, 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))`
`+	regularized <- Lhat + Shat`
`+	return( regularized )`
`+}`

# File simulation.r

`+source("regularization_functions.r")`
`+require(picante)`
`+require(ape) `
`+require(geiger)`
`+require(lorec)`
`+require(fields)`
`+`
`+set.seed(123456)`
`+`
`+`
`+# A couple functions to facilitate simulation`
`+`
`+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 ) `
`+}`
`+`
`+plot_matrix <- function(m, ... ) {`
`+	# Basic plotting`
`+	nr <- nrow(m)`
`+	nc <- ncol(m)`
`+	image(1:nc, 1:nr, t(m[nr:1, ]), axes=F,xlab="", ylab="", ... )`
`+}`
`+`
`+`
`+# 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`
`+`
`+n <- nrow( contrasts )`
`+`
`+# Regularization by thresholding`
`+`
`+bickel <- regularization_by_thresholding( contrastcor, n )`
`+`
`+# Regularization by convex Minimization`
`+`
`+luo <- regularization_by_convex_minimization( contrastcor, n )`
`+`
`+`
`+# 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()`