Friday, August 13, 2010

R::Hotelling's T test

> ab <- scan("c:/r/ttest.dat")
Read 189 items
> dtemp <- t(matrix(ab,3))
> k <- ncol(dtemp)
> d1 <- dtemp[1:25,1:3]
> d2 <- dtemp[26:63,1:3]
> xbar1 <- apply(d1,2,mean)
> xbar2 <- apply(d2,2,mean)
> dbar <- xbar2-xbar1
> round(cbind(xbar1,xbar2,dbar),3)
[,1] [,2] [,3]
[1,] 70.28 69.76316 -0.5168423
[2,] 174.40 177.44737 3.0473687
[3,] 47.12 49.94737 2.8273687
> v <- ((n1-1)*var(d1)+(n2-1)*var(d2))/(n1+n2-2)
> round(v,3)
[,1] [,2] [,3]
[1,] 6.916531 33.413546 1.224366
[2,] 33.413546 457.596635 6.372045
[3,] 1.224366 6.372045 37.484176
> t2 <- n1*n2*dbar%*%solve(v)%*%dbar/(n1+n2)
> f <- (n1+n2-k-1)*t2/((n1+n2-2)*k)
> pvalue <- 1-pf(f,ncol(dtemp),n1+n2-k-1)
> cbind(f,pvalue)
[,1] [,2]
[1,] 1.797264 0.1575505
>

THE DATA FILE (ttest.dat)

73 165 54
70 175 59
67 157 44
70 168 40
68 180 59
70 165 39
70 200 40
68 175 46
72 178 47
70 170 46
69 167 49
71 145 41
74 172 48
68 162 43
77 210 49
70 170 48
71 180 44
71 170 39
68 153 52
74 226 57
70 180 39
71 174 53
68 178 43
70 140 45
67 200 54
65 150 54
67 173 49
68 173 45
74 185 41
69 167 44
76 265 53
70 160 52
70 172 58
69 165 59
63 155 57
68 164 51
69 190 43
70 160 52
72 200 59
72 176 58
70 180 55
72 160 59
71 168 54
63 137 48
70 175 52
68 185 49
71 228 40
72 180 48
74 200 52
72 191 44
68 178 40
72 190 52
66 160 39
69 165 49
70 165 52
69 145 52
69 168 40
74 205 51
69 185 59
68 182 49
70 188 51
72 188 41
70 165 47




No comments:

Post a Comment