David,
2.3.1 is a bit old to be reporting bugs -- we do ask people to check that their
problem is still present in a contemporary version of R. However, your data do
still give the same output in R 2.7.2 (which is not current, but was current
less than a year ago).
I've tidied up the code to remove all the weird characters:
booth<-read.table(tmp<-textConnection("COUNT REWGRP COMMIT SATIS STAY
1 1 16 19 18
2 1 18 15 17
3 1 18 14 14
4 1 16 20 10
5 1 15 13 17
6 1 12 15 11
7 2 16 20 13
8 2 18 14 16
9 2 13 10 14
10 2 17 13 19
11 2 14 18 15
12 2 19 16 18
13 3 20 18 16
14 3 18 15 19
15 3 13 14 17
16 3 12 16 15
17 3 16 17 18
18 3 14 19 15
"),header=TRUE)
fit<-manova(cbind(COMMIT,SATIS,STAY)~REWGRP,data=booth)
Now, as to the question of whether this is a bug. You don't give the SAS
answers that you are happy with, just the R answers. This makes it a lot more
difficult.
It's possible that there is a bug in manova(), but another possibility,
since you are concerned about degrees of freedom, and based on the last three
letters of the name of your predictor variable, is that you wanted
> fit2<-manova(cbind(COMMIT,SATIS,STAY)~factor(REWGRP),data=booth)
> summary(fit2, test="Pillai")
Df Pillai approx F num Df den Df Pr(>F)
factor(REWGRP) 2 0.28342 0.77049 6 28 0.5995
Residuals 15 > summary(fit2, test="Roy")
Df Roy approx F num Df den Df Pr(>F)
factor(REWGRP) 2 0.31963 1.49159 3 14 0.2599
Residuals 15 > summary(fit2, test="Hotelling")
Df Hotelling-Lawley approx F num Df den Df Pr(>F)
factor(REWGRP) 2 0.36261 0.72521 6 24 0.6336
Residuals 15
Googling for the variable names and SAS, MANOVA found some programs in which
REWGRP was specified as a CLASS variable, ie, a factor.
Also
http://my.safaribooksonline.com/9781590474174/ch11lev1sec3
has what might be the output of this code. The test statistics all match the
ones in R, but the p-values are slightly different except for Wilks' lambda.
So, it looks as though at least you need to specify that the variable is a
factor. I will have to leave the question of whether the p-values are correct
to someone with more knowledge of MANOVA. It does seem from the documentation
that agreement with SAS is intended at least for the Pillai trace and Roy's
largest root.
We do appreciate bug reports, but it shouldn't be necessary to do all this
work to find out what you think the correct answer is.
-thomas
On Mon, 16 Mar 2009 dvdbooth at cs.com wrote:
>
> Hi.? There appears to be a bug in R function manova.? My friend and I both
ran it the same way as shown below (his run) with the shown data set. His
results are shown below. we both got the same results.? I was running with R
2.3.1. I'm not sure what version he used.
> Thanks very much,
> David Booth
> Kent State University
>
>
>
>
>
>
>
> -----Original Message-----
> From: dvdbooth at cs.com
> To: kberk at ilstu.edu
> Sent: Sun, 15 Mar 2009 7:01 pm
> Subject: Re: MANOVA Data
>
>
>
>
>
>
>
>
>
>
>
> Ken,
>
> Did you notice that Wilks, Roy, etc p-values are all the same?? Pillai is
almost the SAS result.? Can't figure it out.? I'll submit a bug report.
What's Velleman going to talk about?? Thanks for looking at the R.
>
> Best,
>
> Dave
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> -----Original Message-----
>
> From: Ken Berk <kberk at ilstu.edu>
>
> To: dvdbooth at cs.com
>
> Sent: Sun, 15 Mar 2009 3:45 pm
>
> Subject: Re: Fwd: MANOVA Data
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> At 08:07 PM 3/5/2009, you wrote:
>
>
>
> Hi Ken,
>
>
> I've run the attached data set ( a one way MANOVA ex. from the SAS
manual
> chapter on MANOVA) in both SAS and R and I don't get the same
> results.? Do you have any suggestions about how I can find out
> what's going on?
>
>
> Thanks,
>
>
> Dave
>
>
>
>
>
>
>
> -----Original Message-----
>
>
> From: dvdbooth at cs.com
>
>
> To: DvdBooth at aol.com
>
>
> Sent: Thu, 5 Mar 2009 5:06 pm
>
>
> Subject: MANOVA Data
>
>
>
>
>
>
>
>
>
>
> Email message sent from CompuServe - visit us today at
> http://www.cs.com
>
>
>
>
>
> Email message sent from CompuServe - visit us today at
> http://www.cs.com
>
>
>
>
>
>
>
>
> Hello, David
>
>
>
>
> My R results are clearly crap, as shown below.
>
>
>
>
> The degrees of freedom are clearly wrong, as is apparent when looking at
> the univariate anovas.
>
>
>
>
> SAS gives the correct answers.
>
>
>
>
> I don't know what to do about R.
>
>
>
>
> Ken
>
>
>
>
>
>
>
> COUNT??? REWGRP??? COMMIT???
> SATIS??? STAY
>
>
> 1????
> 1????????
> 16???????
> 19?????? 18
>
>
> 2????
> 1????????
> 18???????
> 15?????? 17
>
>
> 3????
> 1????????
> 18???????
> 14?????? 14
>
>
> 4????
> 1????????
> 16???????
> 20?????? 10
>
>
> 5????
> 1????????
> 15???????
> 13?????? 17
>
>
> 6????
> 1????????
> 12???????
> 15?????? 11
>
>
> 7????
> 2????????
> 16???????
> 20?????? 13
>
>
> 8????
> 2????????
> 18???????
> 14?????? 16
>
>
> 9????
> 2????????
> 13???????
> 10?????? 14
>
>
> 10??? 2????????
> 17???????
> 13?????? 19
>
>
> 11??? 2????????
> 14???????
> 18?????? 15
>
>
> 12??? 2????????
> 19???????
> 16?????? 18
>
>
> 13??? 3????????
> 20???????
> 18?????? 16
>
>
> 14??? 3????????
> 18???????
> 15?????? 19
>
>
> 15??? 3????????
> 13???????
> 14?????? 17
>
>
> 16??? 3????????
> 12???????
> 16?????? 15
>
>
> 17??? 3????????
> 16???????
> 17?????? 18
>
>
> 18??? 3????????
> 14???????
> 19?????? 15
>
>
>
>
>> attach(booth)
>
>
>> Y <- cbind(COMMIT, SATIS, STAY)
>
>
>> fit <- manova(Y ~ REWGRP)
>
>
>> summary(fit, test="Pillai")
>
>
> ????????? Df? Pillai
> approx F num Df den Df Pr(>F)
>
>
> REWGRP???? 1 0.22731?
> 1.37283????? 3???? 14
> 0.2918
>
>
> Residuals
> 16?????????????????????????????????????
>
>
>
>> summary(fit, test="Wilks")
>
>
> ????????? Df??
> Wilks approx F num Df den Df Pr(>F)
>
>
> REWGRP???? 1 0.77269?
> 1.37283????? 3???? 14
> 0.2918
>
>
> Residuals
> 16?????????????????????????????????????
>
>
>
>> summary(fit, test="Hotelling-Lawley")
>
>
> ????????? Df
> Hotelling-Lawley approx F num Df den Df Pr(>F)
>
>
> REWGRP????
> 1????????? 0.29418?
> 1.37283????? 3???? 14
> 0.2918
>
>
> Residuals
> 16??????????????????????????????????????????????
>
>
>
>> summary(fit, test="Roy")
>
>
> ?????????
> Df???? Roy approx F num Df den Df Pr(>F)
>
>
> REWGRP???? 1 0.29418?
> 1.37283????? 3???? 14
> 0.2918
>
>
> Residuals
> 16?????????????????????????????????????
>
>
>
>> summary(fit)
>
>
> ????????? Df? Pillai
> approx F num Df den Df Pr(>F)
>
>
> REWGRP???? 1 0.22731?
> 1.37283????? 3???? 14
> 0.2918
>
>
> Residuals
> 16?????????????????????????????????????
>
>
>
>> summary.aov(fit)
>
>
> ?Response COMMIT :
>
>
> ???????????
> Df? Sum Sq Mean Sq F value Pr(>F)
>
>
> REWGRP?????? 1??
> 0.333?? 0.333? 0.0532 0.8204
>
>
> Residuals?? 16 100.167??
> 6.260??????????????
>
>
>
>
>
> ?Response SATIS :
>
>
> ???????????
> Df? Sum Sq Mean Sq F value Pr(>F)
>
>
> REWGRP?????? 1??
> 0.750?? 0.750? 0.0945 0.7625
>
>
> Residuals?? 16 127.028??
> 7.939??????????????
>
>
>
>
>
> ?Response STAY :
>
>
> ??????????? Df Sum
> Sq Mean Sq F value Pr(>F)
>
>
> REWGRP?????? 1 14.083? 14.083?
> 2.3013 0.1488
>
>
> Residuals?? 16 97.917??
> 6.120??????????????
>
>
>
>
>
>>
>
>
>
>
>
>
>
>
>
>
>
> Email message sent from CompuServe - visit us today at http://www.cs.com
>
>
>
>
>
>
>
> ________________________________________________________________________
> Email message sent from CompuServe - visit us today at http://www.cs.com
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
Thomas Lumley Assoc. Professor, Biostatistics
tlumley at u.washington.edu University of Washington, Seattle