Skip to content

stats/debugging question hotelling t-sq

2 messages · toby909 at gmail.com, Klaus Nordhausen

toby909 at gmail.com
# 
Hi

I spent hours looking over my formula. Somehow I cant find the reason
why it gives me different answer.

help appreciated.




x = as.matrix(read.table("http://www.niehs.nih.gov/research/atniehs/core/microarrays/docs/heinloth.txt",1))
x = t(x)    #now rows are subjects, cols are genes
x = x[order(rownames(x)),]    #order by treatment group oxygen,
ultra-violet, gamma radiation
y = x[26:45,1:10]
x = x[2:25,1:10]
p = ncol(x); p
nx = nrow(x); nx
ny = nrow(y); ny
n = nx+ny; n

# (t(x)-colMeans(x)) %*% t(t(x)-colMeans(x))

T2 = nx*ny/n * t(colMeans(x)-colMeans(y)) %*% solve( (
(nx-1)*cov(x)+(ny-1)*cov(y) )/( n-2 ) ) %*% (colMeans(x)-colMeans(y));
T2

library(ICSNP)
HotellingsT2(y,x)





http://en.wikipedia.org/wiki/Hotelling's_T-square_distribution

http://finzi.psych.upenn.edu/R/Rhelp02a/archive/60962.html
Klaus Nordhausen
# 
Hi!

I think the results agree:

using simulated data:

set.seed(1)
library(mvtnorm)
x<-rmvnorm(44,rep(0,10))
y = x[(26:45)-1,1:10]
x = x[(2:25)-1,1:10]
p = ncol(x); p
nx = nrow(x); nx
ny = nrow(y); ny
n = nx+ny; n

# (t(x)-colMeans(x)) %*% t(t(x)-colMeans(x))

T2 = nx*ny/n * t(colMeans(x)-colMeans(y)) %*% solve( (
(nx-1)*cov(x)+(ny-1)*cov(y) )/( n-2 ) ) %*% (colMeans(x)-colMeans(y));
T2

library(ICSNP)
HotellingsT2(y,x)

note that the default here returns the test statistic value that is F distributed. So using 

HotellingsT2(y,x,test="chi")

gives the same.

Or if you transform your T2 to be F distributed

(n-p-1) / ((n-2)*p) * T2

Best wishes,

Klaus



-------- Original-Nachricht --------