"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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.