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)
}

flattr this!

To leave a comment for the author, please follow the link and comment on his blog: fernandohrosa.com.br - en » R.

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



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.