# Sphericity Test for Covariance Matrices in R (sphericity.test)

October 9, 2011
By

(This article was first published on fernandohrosa.com.br - en » R, and kindly contributed to R-bloggers)

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 sphericity.test function, so here it is, for posterity:

 ```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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 ``` ```sphericity.test<-function(n,s1,s2=NULL,estsigma=TRUE){ #### Performs a hypothesis test that a covariance matrix is of specified #### form. Test is of the form H0: S1=sigma^2*S2. n is the number of #### observations on which the sample covariance matrix is based. #### If the input parameter estsigma is TRUE: #### Perform test of the hypothesis that S1=sigma^2 S2, for unknown sigma. #### If S2 not specified, assumed that S2=I. Reference is Basilevsky, #### Statistical Factor Analysis and Related Methods, page 191. #### If the input parameter estsigma is FALSE: #### Perform test of the hypothesis that S1=S2. If S2 not specified, #### assumed that S2=I. Reference is Seber, Multivariate Observations, #### sec 3.5.4 #### Only the lower triangle+diagonal is required at entry, and the upper #### triangle is ignored. #### DAW July 2000 dname <- paste(substitute(s1)) p<-nrow(s1) for (i in 1:(p-1)){for (j in ((i+1):p)){ s1[i,j]<-s1[j,i] s2[i,j]<-s2[j,i] }} if (!is.null(s2)){ b<-eigen(s2,symmetric=T,only.values=F) r<-b\$vectors %*% diag(1/sqrt(b\$values)) s<-t(r) %*% s1 %*% r } else { s<-s1 }   d<-eigen(s,symmetric=T,only.values=T)\$values ldet<-sum(log(d)) tr<-sum(d)   if (estsigma==TRUE){ sighat<-tr/p cc<--(n-(2*p^2+p+2)/(6*p))*(ldet-p*log(tr/p)) statistic <- cc sighat<-sighat names(statistic) <- "L statistic" parameter <- 0.5*(p+2)*(p-1) names(parameter) <- "df" rval<-list(data.name=dname,sighat=sighat,statistic=statistic,parameter=parameter,p.value=1-pchisq(statistic,parameter),method="Sphericity test") } else { cc<--n*(p+ldet-tr) statistic <- cc names(statistic) <- "L statistic" parameter <- 0.5*(p+1)*p names(parameter) <- "df" rval<-list(data.name="",statistic=statistic,parameter=parameter,p.value=1-pchisq(statistic,parameter),method="Covariance equality test statistic") } class(rval) <- "htest" return(rval) } pw <- function(q,n) { pdf <- function(w) { 1/2 * (n-2) * w^((n-3)/2) } integrate(pdf,0,q) }   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) }```

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...