[R] Correlation matrices

[David Wooff]

Patrik Waldmann wrote:
> 
> Hello,
> 
> are there any good methods in R that will test if two correlation matrices (obtained in different ways) are equal? Something better than the Mantel test would be preferable.
> 
> Regards,
> 
> Patrik Waldmann
> 

Not that I know of. I include a function I've written for my own use,
assuming multivariate Normality. I'm sure the experts could improve on
the function. (And please could they or you point out any embarrasing
errors!)

David.

-- 
  David Wooff, Director, Statistics and Mathematics Consultancy Unit,
  Department of Mathematical Sciences, University of Durham.
  Science Laboratories, South Road, Durham, DH1 3LE, UK.
  Tel. 0191 374 4531, Fax 0191 374 7388.


sphericity.daw<-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

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) <- "Sphericity test statistic"
parameter <- 0.5*(p+2)*(p-1)
names(parameter) <- "df"
rval<-list(data.name="",sighat=sighat,statistic=statistic,parameter=parameter,p
value=1-pchisq(statistic,parameter))
}
else
{
cc<--n*(p+ldet-tr)
statistic <- cc
names(statistic) <- "Covariance equality test statistic"
parameter <- 0.5*(p+1)*p
names(parameter) <- "df"
rval<-list(data.name="",statistic=statistic,parameter=parameter,p.value=1-pchis
q(statistic,parameter))
}
class(rval) <- "htest"
return(rval)
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._