# This function tests the null hypothesis of no divergence between two correlation matrices using matrix correlation as the similarity index and parametric bootstrap as the testing procedure. H0 is that the matrix correlation (MR) is equal to 1.
# Input is the two correlation matrices ("R1","R2") and sample size for both ("N1","N2"). The number of bootstrap iterations ("n.boot") can also be specified; the default is 1000.
# The analysis is run by typing "bootMR1(R1,R2,N1,N2)".
# Output is the observed correlation ("r.obs"), the repeatabilities for each of the matrices ("t1" and "t2"), and the two p-values, one p based on R1 ("p1") and one p based on R2 ("p2"). The repeatabilities are in fact based on n.boot*2 iterations.

# A summary of the code follows:
#	X1 and X2 are n-by-p matrices of numbers drawn randomly from a multivariate normal distribution.
#	The function "rmvnorm" is constrained by the "sigma" term to impose the correlation structure of observed matrix R onto the X1 or X2 matrix.
#	The matrix correlation is then computed for this pair of matrices.
#	This procedure is repeated 1000 times, generating a distribution of matrix correlation values between matrices of the same structure and differing only in random, normally distributed noise.
#	The matrix correlation between the two observed correlation matrices is then compared to the distribution.
#	The probability of achieving the matrix correlation between the observed correlation matrices if those matrices were drawn from the same population is given by the proportion of the bootstrapped distribution having a lower matrix correlation value.

# This code was generously provided my Annat Haber (Ph.D. candidate, CEB, University of Chicago) in August 2010.



#### auxiliary functions
require(mvtnorm)

matcor <- function(R1, R2) {cor(R1[lower.tri(R1)], R2[lower.tri(R2)])}

boot1 <- function(N, R) {
		X1 <- rmvnorm(N, sigma=R)
		X2 <-rmvnorm(N, sigma=R)
		rb <- matcor(cor(X1), cor(X2))
		tb1 <- matcor(cor(X1),R)
		tb2 <- matcor(cor(X2),R)
		c(rb,tb1,tb2)}

#### the actual bootH0MR1 function which wraps up the previous ones
bootMR1 <- function(R1,R2,N1,N2, n.boot=1000) {
	r.obs <- matcor(R1,R2)
	Nb <- rep(N1, n.boot)
	Drt <- sapply(Nb, boot1, R=R1)
	br <- sort(Drt[1,])
	p1 <- length(br[br<=r.obs])/(length(br)+1)
	t1 <- mean(Drt[2:3,])
	Nb <- rep(N2, n.boot)
	Drt <- sapply(Nb, boot1, R=R2)
	br <- sort(Drt[1,])
	p2 <- length(br[br<=r.obs])/(length(br)+1)
	t2 <- mean(Drt[2:3,])	
	list(r.obs=r.obs, t1=t1, t2=t2, p1=p1, p2=p2)}