Posted by fernandohrosa on Oct 9, 2011 in R, Statistics Tags: multivariate analysis, multivariate testing, R | 9 comments

This is a piece of code I implemented in 2004, which was supposed to be part of an R-package in multivariate testing (to be named, rather creatively, mvttests).

Time has flown, I haven’t still got around to implementing the said package, but people keep asking me for the varcomp function, so here it is, for posterity:

Download varcomp.R

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 | varcomp <- function(covmat,n) { if (is.list(covmat)) { if (length(covmat) < 2) stop("covmat must be a list with at least 2 elements") ps <- as.vector(sapply(covmat,dim)) if (sum(ps[1] == ps) != length(ps)) stop("all covariance matrices must have the same dimension") p <- ps[1] q <- length(covmat) if (length(n) == 1) Ng <- rep(n,q) else if (length(n) == q) Ng <- n else stop("n must be equal length(covmat) or 1") DNAME <- deparse(substitute(covmat)) } else stop("covmat must be a list") ng <- Ng - 1 Ag <- lapply(1:length(covmat),function(i,mat,n) { n[i] * mat[[i]] },mat=covmat,n=ng) A <- matrix(colSums(matrix(unlist(Ag),ncol=p^2,byrow=T)),ncol=p) detAg <- sapply(Ag,det) detA <- det(A) V1 <- prod(detAg^(ng/2))/(detA^(sum(ng)/2)) kg <- ng/sum(ng) l1 <- prod((1/kg)^kg)^(p*sum(ng)/2) * V1 rho <- 1 - (sum(1/ng) - 1/sum(ng))*(2*p^2+3*p-1)/(6*(p+1)*(q-1)) w2 <- p*(p+1) * ((p-1)*(p+2) * (sum(1/ng^2) - 1/(sum(ng)^2)) - 6*(q-1)*(1-rho)^2) / (48*rho^2) f <- 0.5 * (q-1)*p*(p+1) STATISTIC <- -2*rho*log(l1) PVAL <- 1 - (pchisq(STATISTIC,f) + w2*(pchisq(STATISTIC,f+4) - pchisq(STATISTIC,f))) names(STATISTIC) <- "corrected lambda*" names(f) <- "df" RVAL <- structure(list(statistic = STATISTIC, parameter = f,p.value = PVAL, data.name = DNAME, method = "Equality of Covariances Matrices Test"),class="htest") return(RVAL) } |

Hi Fernando,

I was wondering what test for equality of covariance matrices your varcomp function is based on? I cant really figure it out myself.

Sara

Good question… I have to trace it sometime but I think it was the one described by Anderson (2003): http://www.amazon.com/gp/product/0471360910/002-3606257-8048037?ie=UTF8&tag=feferraznet-20&linkCode=xm2&camp=1789&creativeASIN=0471360910

And also what is the parameter n?

It’s the number of data points you used to calculate each of the matrices…

Hi Fernando,

I am trying to implement the above function but keep receiving this error message:

Error in ps[1] == ps : comparison of these types is not implemented

The code I am trying to use is this:

mn.test <- t(cbind(Shelter$Toptosub, Shelter$Dbot,

Shelter$Rtopratio, Shelter$Rbotratio))

covmat <- as.list(cov(mn.test))

varcomp(covmat, 133) # this is my overall n. I have also tried

# it with length(covmat) and 1.

Thanks for any help and also for sharing this function!

Hi Evan,

Sorry for taking so long to reply. My WordPress install went berserk and I stopped receiving e-mail notices upon new comments.

As to your question: you have only one covariance matrix to compare, that’s why it’s not working. The function is supposed to work on two or more (potentially several) covariance matrices. cov(mn.test) is a single covariance matrix, you’ve got nothing to compare it against!

Fernando,

This function has been implemented correctly. Do you have R code example to work on it after this function in order to understand the steps of covariance testing?

Regards

You just call the function with two arguments, the first is a list of covariance matrices, and n is a vector of the number of observations used to calculate each of the covariance matrices. If n is of length 1, it is assumed that all matrices were calculated from data sets with the same length.

Hi. Thanks. I changed the code to use the logarithm of the determinant in the computations because for larger dimensions and n it wouldnt compute anymore. Below is the change. Seems to work. Ivo.

————– snippet ————–

ng <- Ng – 1

Ag <- lapply(1:length(covmat),function(i,mat,n) { n[i] * mat[[i]] },mat=covmat,n=ng)

A <- matrix(colSums(matrix(unlist(Ag),ncol=p^2,byrow=T)),ncol=p)

log.detAg <- sapply(Ag,function(x) determinant(x, logarithm=TRUE)$modulus)

log.detA <- determinant(A, logarithm=TRUE)$modulus

log.V1 <- as.numeric( sum((ng/2)*log.detAg) – (sum(ng)/2)*log.detA )

kg <- ng/sum(ng)

log.l1 <- (p*sum(ng)/2) * sum( kg*log(1/kg)) + log.V1

rho <- 1 – (sum(1/ng) – 1/sum(ng))*(2*p^2+3*p-1)/(6*(p+1)*(q-1))

w2 <- p*(p+1) * ((p-1)*(p+2) * (sum(1/ng^2) – 1/(sum(ng)^2)) – 6*(q-1)*(1-rho)^2) / (48*rho^2)

f <- 0.5 * (q-1)*p*(p+1)

STATISTIC <- -2*rho*log.l1