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 --------
> Datum: Sun, 16 Mar 2008 04:09:00 -0700
> Von: [EMAIL PROTECTED]
> An: r-help@r-project.org
> Betreff: [R] stats/debugging question hotelling t-sq
> 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
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Psst! Geheimtipp: Online Games kostenlos spielen bei den GMX Free Games!
http://games.entertainment.gmx.net/de/entertainment/games/free
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
