[R] Hotelling's T^2 Test (here is the code)

[Peter B. Mandeville]

At 09:28 AM 16/11/98 +0100, you wrote:
># Peter B. Mandeville kindly offered to send me code for Hotelling's T^2
># Test. Unfortunately there seems to be no route to his machine.
># So i'm trying to reach him via the Mailing List.
>
>------------------------------------------------------------------------
>
>Sir,
>
>this morning i recieved your message about the availability of the code
>for Hotelling's Test. I hurried to find out whether the book you cited is
>available at german universities/libraries. This seems not to be the case
>yet. 
>
>So the transmission of the code would be appreciated a lot. Many, many
>thanks for your offer. 
>
>Could i ask you for a little comment on Steve Selvin's publication? If
>it's not available publicly in Germany within the next time, i probably
>will try to get it via a bookstore at my own cost.
>
>Regards,
>
>Ulrich Flenker
>______________
>
>Institute of Biochemistry
>German Sports University Cologne
>
>
>-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
.-.-
>r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
>Send "info", "help", or "[un]subscribe"
>(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
>_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
._._
>

Here is the code
Peter B.

Steve Selvin
1998 Modern Applied Biostatistical Methods Using S-Plus
Monographs in Epidemiology and Biostatistics Volume 28
Oxford University Press
pages 54-57

SLIGHTLY MODIFIED TO WORK IN R, GETS THE SAME RESULTS

> 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




-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._