Commits

Casey Dunn  committed e064ae1

renamed files for clarity

  • Participants
  • Parent commits f7e2b2b

Comments (0)

Files changed (4)

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()