[R] TukeyHSD model thing

[Chuck Cleland]

On 3/6/2010 4:38 PM, casperyc wrote:
> Hi,
> 
> I am trying to reproduce a tukey test in R
> 
> ==========================
> x=c(145,40,40,120,180,
> 	140,155,90,160,95,
> 	195,150,205,110,160,
> 	45,40,195,65,145,
> 	195,230,115,235,225,
> 	120,55,50,80,45
> 	)
> y2=c(
> 	rep(as.character(1),5),
> 	rep(as.character(2),5),
> 	rep(as.character(3),5),
> 	rep(as.character(4),5),
> 	rep(as.character(5),5),
> 	rep(as.character(6),5)
> 	)
> 
> crd2=data.frame(x,y2)
> 
> model1=aov(x~y2,data=crd2)
> TukeyHSD(model1)
> 
> ==========================
> 
> The result in the 'diff' of the means are the same as those did using SAS,
> (which is in my tutorial sheet, i got a MAC, so cant use SAS)
> however, the 95% confiden limits are quite different.
> 
> ===========================
> 
> 2-1   23  -73.817441 119.817441 0.975518208
> 3-1   59  -37.817441 155.817441 0.435116628
> 4-1   -7 -103.817441  89.817441 0.999912318
> 5-1   95   -1.817441 191.817441 0.056613465
> 6-1  -35 -131.817441  61.817441 0.869242006
> 3-2   36  -60.817441 132.817441 0.855531189
> 4-2  -30 -126.817441  66.817441 0.926612938
> 5-2   72  -24.817441 168.817441 0.232896275
> 6-2  -58 -154.817441  38.817441 0.453535553
> 4-3  -66 -162.817441  30.817441 0.316718292
> 5-3   36  -60.817441 132.817441 0.855531189
> 6-3  -94 -190.817441   2.817441 0.060579795
> 5-4  102    5.182559 198.817441 0.034819938
> 6-4  -28 -124.817441  68.817441 0.944203446
> 6-5 -130 -226.817441 -33.182559 0.004294761
> 
> ===========================
> 
> in the SAS output, it's
> (in slightly different order, you can just check one of the set)
> ===========================
> 5 - 3 36.00 -28.63 100.63
> 5 - 2 72.00 7.37 136.63 ***
> 5 - 1 95.00 30.37 159.63 ***
> 5 - 4 102.00 37.37 166.63 ***
> 5 - 6 130.00 65.37 194.63 ***
> 3 - 5 -36.00 -100.63 28.63
> 3 - 2 36.00 -28.63 100.63
> 3 - 1 59.00 -5.63 123.63
> 3 - 4 66.00 1.37 130.63 ***
> 3 - 6 94.00 29.37 158.63 ***
> 2 - 5 -72.00 -136.63 -7.37 ***
> 2 - 3 -36.00 -100.63 28.63
> 2 - 1 23.00 -41.63 87.63
> 2 - 4 30.00 -34.63 94.63
> 2 - 6 58.00 -6.63 122.63
> 1 - 5 -95.00 -159.63 -30.37 ***
> 1 - 3 -59.00 -123.63 5.63
> 1 - 2 -23.00 -87.63 41.63
> 1 - 4 7.00 -57.63 71.63
> 1 - 6 35.00 -29.63 99.63
> 4 - 5 -102.00 -166.63 -37.37 ***
> 4 - 3 -66.00 -130.63 -1.37 ***
> 4 - 2 -30.00 -94.63 34.63
> 4 - 1 -7.00 -71.63 57.63
> 4 - 6 28.00 -36.63 92.63
> 6 - 5 -130.00 -194.63 -65.37 ***
> 6 - 3 -94.00 -158.63 -29.37 ***
> 6 - 2 -58.00 -122.63 6.63
> 6 - 1 -35.00 -99.63 29.63
> 6 - 4 -28.00 -92.63 36.63
> ===========================
> 
> say, betweet treatment 5 and 3
> 
> R
> 5-3   36       -60.817441 132.817441
> SAS
> 5-3   36       -28.63       100.63
> 
> i am wondering if i have done something wrong in R?

  It looks like the confidence intervals for your SAS output are not
family-wise intervals.  For example, you can get intervals that match
your SAS output when you don't adjust for multiple comparisons.

> confint(lm(x ~ as.factor(y2), data=crd2))
                    2.5 %    97.5 %
(Intercept)     59.302005 150.69800
as.factor(y2)2 -41.626725  87.62672
as.factor(y2)3  -5.626725 123.62672
as.factor(y2)4 -71.626725  57.62672
as.factor(y2)5  30.373275 159.62672
as.factor(y2)6 -99.626725  29.62672

> full documents are in the attachment,
> can someone suggest to me the relevent R codes
> that does the same sort of thing?
> (tukeyHSD,fisherLSD,and anova table )
> 
> Thanks!
> 
> casper http://n4.nabble.com/file/n1583109/SAS.pdf SAS.pdf 
> http://n4.nabble.com/file/n1583109/R.pdf R.pdf 
> http://n4.nabble.com/file/n1583109/ws1.pdf ws1.pdf 
> http://n4.nabble.com/file/n1583109/ws1sols.pdf ws1sols.pdf

-- 
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894