1

How to create a covariance matrix by eigen-values, so that mauchly's test indicates to reject the null hypothesis. Is there a pattern of a covariance-matrix for a experimental design which produces r.v., so that the assumption of sphericity has been violated? If I calculated a arbitrary covariance-matrix by given eigen-values, and for that matrix I produce an experiment by the following source-code, then mauchlys test rejects the null hypothesis sometimes. How to simulate normally distributed data which violates the assumption of sphericity?

require(corpcor)
require(reshape2)
require(mvtnorm)

getCovariance<-function(lambda){
   j<-1
   Sigma<--1
   while(any(Sigma<=0)){    
      set.seed(j)
  X<-matrix(rnorm(50*length(lambda)),50,length(lambda))  
  R<-cor(X)           
  P<-eigen(R)$vector 
  Sigma<-t(P)%*%diag(lambda)%*%P   
      j<-j+1
   }             
   Sigma
 }

 subject<-50
 treatment<-4
 set.seed(1)
 lambda<-sample(1:subject,treatment)
 cov<-getCovariance(lambda)
 exp<-rmvnorm(mean=rep(0,treatment),sigma=cov,n=subject,method="chol")
 fit<-lm(exp~1)
 mauchly.test(fit,X=~1)
Klaus
  • 1,946
  • 3
  • 19
  • 34
  • I would think that you would need to draw from a non-Normal distribution or to create a non-linear dependence. Also note that `cov` and `exp` are essential R functions and creating a data object with those names is going to create confusion in the human readers of your code. There's even a fortune to the effect: `fortunes::fortune("dog")` – IRTFM Jun 21 '13 at 18:11
  • I dont see how the choice of distribution, in getCovariance, makes the dependencies more virtuos. Or did you mean, to use a other distribution instead of rmvnorm? But this is not my intention. I know about the programming conventions, but in this case I have a other point of view. – Klaus Jun 21 '13 at 20:17
  • The second. A distribution where you would get a covariance structure that was not ... Normal. "Another point of view"? Can you be any more vague? – IRTFM Jun 21 '13 at 20:20
  • You mean I should use an other generating function for my experimental design. But this creates wrong data. I will use the data to simulate Anova's and the distribution of a Anova teststatistic. therefore have I have to use rmvnorm because the theory is based on quadratic forms in k-variate normal distributed r.v. – Klaus Jun 21 '13 at 20:29
  • With k-variate normal distributed r.v. I ment arbitrary normal distributed r.v., I ment not independent r.v. that follows from the assumption of the linear model, where the distribution of the teststatistic of a classic anova is based on. – Klaus Jun 21 '13 at 20:42
  • You don't have any groups , i.e. `treatment` or `lambda` in your models, so maybe you should review what you are actually supposed to be testing. I would have expected that `lambda` would be on the rhs if you were looking for group differences in variances. – IRTFM Jun 21 '13 at 21:05
  • I do have groups, my design is the same if I would calculate `library(car) a<-Anova(fit,type=3,idata=data.frame(treatment=factor(1:treatment)),idesign=~treatment) summary(a,multivariate=FALSE)` there you see the same mauchly test values. I think that mauchly's test is not very useful for a sample size like 50. Because if mauchlys test rejects the alternativ. It means there is a transformation to uncorrelated r.v. and in this case it would possible to use a standard anova `anova(lm(response~treatment,data.frame(response=c(exp),treatment=factor(rep(1:4,each=subject)))))`. – Klaus Jun 22 '13 at 09:19
  • But if you look closer to the distribution of the teststatistic of this anova you will see it doesnt work, because of the correlation in the data. Also 50 is not a small number for anova studies, i.g. in 'Girden, E.R. (1992). ANOVA: repeated measures. Newbury Park, CA: Sage.' are many experiments with a smaller number of subjects. – Klaus Jun 22 '13 at 09:20

0 Answers0