CODEWITHSUNDEEP
r
Klaus
Klaus

Reputation: 2056

Hotellings statistic

I expected I could proof, if a sample of p-variate normal random vectors has a theoretical mean value, by the hotellings test. But a crosscheck with ks.test, if the distribution from the HottelingsT2 function matches the distribution of the test statistic which is used by HottelingsT2-Test failed. It means that the simulated experiments has not the mean value 0, but obviously they have. So there should be something wrong in context. Are there some mistakes?

require(mvtnorm)
require(ICSNP)
subject<-50
treatment<-4
V<-matrix(c(644.03100226056, 184.319025225855, 572.5312199559, 143.106678641056, 184.319025225855, 73.5310268006399, 230.838267981476, 130.977532385651, 572.5312199559, 230.838267981476, 736.378779002912, 429.445506266528, 143.106678641056, 130.977532385651, 429.445506266528, 435.124191935888),treatment,treatment)

experiment<-list()
R<-3000
seed<-split(1:(R*subject),1:R)
for(i in 1:R){
  e<-c()
  for(j in 1:subject){
    set.seed(seed[[i]][j]) 
    e<-c(e,rmvnorm(mean=rep(0,treatment),sigma=V,n=1,method="chol"))
   }
   experiment<-c(experiment,list(matrix(e,subject,treatment,byrow=T)))
 }

 p.values<-c()
 for(e in experiment){
   fit<-lm(e~1)
   p.values<-c(p.values,HotellingsT2(e, mu=rep(0,treatment))[["p.value"]])
 }

 ks.test(p.values, punif,alternative = "two.sided")

Upvotes: 2

Views: 267

Answers (2)

Jon
Jon

Reputation: 763

Hong Ooi is correct about this being a problem with set.seed. I ran your code as it was posted and got the following results:

> ks.test(p.values, punif,alternative = "two.sided")

    One-sample Kolmogorov-Smirnov test

data:  p.values
D = 0.0615, p-value = 2.729e-10
alternative hypothesis: two-sided

But if you change your code such that:

... everything the same before here ...
experiment <- list()
R <- 3000 # experiment
set.seed(42) # set new seed
for (i in 1:R) { # for each of 3000 experiments
  e <- c() # empty vector
  for (j in 1:subject){ # for each of 50 subjects
    e <- c(e,rmvnorm(mean=rep(0,treatment),sigma=V,n=1,method="chol"))
  }
  experiment <- c(experiment,list(matrix(e,subject,treatment,byrow=T)))
}
... everything the same after here ...

Then you get the following:

> ks.test(p.values, punif,alternative = "two.sided")

    One-sample Kolmogorov-Smirnov test

data:  p.values
D = 0.0122, p-value = 0.7613
alternative hypothesis: two-sided

Essentially by continuing to set random seeds anew at every iteration, even though you are careful to choose a different value, you are still removing the independence of the successive draws.

Upvotes: 4

Hong Ooi
Hong Ooi

Reputation: 57696

I haven't checked the code, but I wouldn't be surprised if this was the same problem as described in Klaus' other post: Using Kolmogorov Smirnov Test in R. Basically, don't put set.seed in the middle of the loop: set it once, at the top of the code, and leave it alone afterwards.

Upvotes: 2

Related Questions