R/cortest.mat.R

In psych: Procedures for Psychological, Psychometric, and Personality Research

"cortest.mat" <- 
function(R1,R2=NULL,n1=NULL,n2 = NULL) {
cl <- match.call()
p <- dim(R1)[2]
if(dim(R1)[1] != p) { n1 <- dim(R1)[1]
                    R1 <- cor(R1,use="pairwise")
                    warning ("R1 matrix was not square, correlations found") 
                    }
 if(!is.matrix(R1)) R1 <- as.matrix(R1)  #in case R1 is a data.frame
 
 if (is.null(n1)) {warning("n1 not specified, set to 100")
            n1 <-100}
 if (is.null(R2)) {message("Bartlett's test of is R = I") 
                  detR1 <- det(R1)
                chi2  <- -log(detR1) *(n1 -1 - (2*p + 5)/6)
                  df <- p * (p-1)/2
                  pval <- pchisq(chi2,df,lower.tail=FALSE)
                  n.obs <- n1
                } else {
                       if(dim(R2)[1] != dim(R2)[2] ) {n2 <- dim(R2)[1]
                         R2 <- cor(R2,use="pairwise") 
                    warning ("R2 matrix was not square, correlations found") }
                    
if (p != dim(R2)[2]) stop("correlation matrices R1 and R2 must be of the same size!")
 if(!is.matrix(R2)) R2 <- as.matrix(R2)  #in case R1 is a data.frame
R1.inv <- solve(R1)    #inverse of R
R2.inv <- solve(R2)    #inverse of R

R.inv.2 <- R1.inv %*% R2 #inverse of R1 times R2
R.inv.1 <- R2.inv %*% R1 #inverse of R2 times R1

E1 <- .5*(sum((diag(R.inv.2))) -log(det(R.inv.2)) - p ) #likelihood

E2 <- .5*(sum((diag(R.inv.1))) -log(det(R.inv.1)) - p ) #likelihood

df1 <- p * (p-1)/2
df <- 2*df1
 
if (is.null(n2)) {n2 <- n1}
n.obs <- min(n1,n2)
chi21 <- E1 * (n1-1-(2*p - 5)/6)
chi22 <- E2 * (n2-1-(2*p - 5)/6)
chi2 <- chi21 + chi22} 
 results <- list(chi2 =chi2,prob =pchisq(chi2,df,lower.tail=FALSE), df= df,n.obs=n.obs,Call =cl)
 class(results) <- c("psych","cortest")
 return(results)
}