[R] manova question
Ranjan Maitra
maitra at iastate.edu
Wed Mar 23 18:32:47 CET 2011
Dear John,
Thanks very much! Truly a duh moment...
But I am trying to now use this to test H_0: C\beta = 0. C is a 2x3
matrix, and \beta is 3 x 4.
I have been trying to use the linearHypothesis() function in the car
package. My reading of the help file is that my hypothesis.matrix is
this C. The RHS is by default the zero matrix. Is this correct?
# So I try the following:
morel <- read.table(file = "http://www.public.iastate.edu/~maitra/stat501/datasets/morel.dat", col.names=c("studentgroup", "aptitude", "mathematics", "language", "generalknowledge"))
morel$studentgroup <- as.factor(morel$studentgroup)
library(car)
fit.lm <- lm(cbind(aptitude, mathematics, language, generalknowledge)~studentgroup , data = morel)
C <- matrix(c(0, 1, 0, 0, 0, 1), ncol = 3, by = T)
newfit <- linearHypothesis(model = fit.lm, hypothesis.matrix = C)
print(newfit)
# But the print(newfit) here gives me the same output as:
summary(Manova(fit.lm))
# What am I doing wrong here?
# Many thanks and best wishes,
# Ranjan
On Tue, 22 Mar 2011 20:21:32 -0500 John Fox <jfox at mcmaster.ca> wrote:
> Dear Ranjan,
>
> Looking at your data, studentgroup is a numeric variable, so your call to lm() produces a multivariate regression, not a one-way MANOVA. Here's a correction:
>
> > morel$studentgroup <- as.factor(morel$studentgroup)
> > Manova( lm( cbind(aptitude, mathematics, language, generalknowledge)
> + ~ studentgroup , data = morel), test="Wilks")
>
> Type II MANOVA Tests: Wilks test statistic
> Df test stat approx F num Df den Df Pr(>F)
> studentgroup 2 0.54345 6.7736 8 152 1.384e-07 ***
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Best,
> John
>
> --------------------------------
> John Fox
> Senator William McMaster
> Professor of Social Statistics
> Department of Sociology
> McMaster University
> Hamilton, Ontario, Canada
> http://socserv.mcmaster.ca/jfox
>
>
>
>
> > -----Original Message-----
> > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
> > On Behalf Of Ranjan Maitra
> > Sent: March-22-11 8:36 PM
> > To: 'R-help'
> > Subject: Re: [R] manova question
> >
> > Dear John, Peter and others,
> >
> > So, I now have a query at an even more elementary level and that is
> > regarding my results from anova.mlm() not matching the car package's
> > Manova(). Specifically, I have been trying the following out with regard
> > to a simple one-way MANOVA setup. So, I try out the following using R:
> >
> > ******* R code *******
> >
> > morel <- read.table(file =
> > "http://www.public.iastate.edu/~maitra/stat501/datasets/morel.dat",
> > col.names = c("studentgroup", "aptitude", "mathematics", "language",
> > "generalknowledge"))
> >
> > morel[,1] <- as.factor(morel[,1])
> > fit <- anova.mlm(as.matrix(morel[,-1]) ~ morel[,1])
> >
> > summary(fit, test="Wilks")
> >
> >
> >
> > *** providing the output ***
> >
> >
> > Df Wilks approx F num Df den Df Pr(>F)
> > morel[, 1] 2 0.54345 6.7736 8 152 1.384e-07 ***
> > Residuals 79
> >
> > *** end of output
> >
> >
> >
> > The above is correct, also by doing the calculations "by hand".
> >
> > Then, I use the car package, following the help function on Anova() and
> > do the following:
> >
> >
> > ******* R code ********
> >
> > morel <- read.table(file =
> > "http://www.public.iastate.edu/~maitra/stat501/datasets/morel.dat",
> > col.names=c("studentgroup", "aptitude", "mathematics", "language",
> > "generalknowledge"))
> >
> > library(car)
> > fit1 <- Manova( lm( cbind(aptitude, mathematics, language,
> > generalknowledge) ~ studentgroup , data = morel)) summary(fit1, test =
> > "Wilks")
> >
> >
> > ****** providing the output *****
> >
> >
> > Type II MANOVA Tests:
> >
> > Sum of squares and products for error:
> > aptitude mathematics language generalknowledge
> > aptitude 78506.68 13976.5677 11041.9434 4330.1304
> > mathematics 13976.57 16040.3996 3979.9528 -416.4845
> > language 11041.94 3979.9528 6035.6132 -372.8491
> > generalknowledge 4330.13 -416.4845 -372.8491 7097.9562
> >
> > ------------------------------------------
> >
> > Term: studentgroup
> >
> > Sum of squares and products for the hypothesis:
> > aptitude mathematics language generalknowledge
> > aptitude 1129.7271 996.0542 237.54441 -880.4353
> > mathematics 996.0542 878.1980 209.43741 -776.2594
> > language 237.5444 209.4374 49.94777 -185.1266
> > generalknowledge -880.4353 -776.2594 -185.12655 686.1536
> >
> > Multivariate Test: studentgroup
> > Df test stat approx F num Df den Df Pr(>F)
> > studentgroup 1 0.8620544 3.080378 4 77 0.020805 *
> > ---
> > Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> >
> >
> > ***** end of output.
> >
> >
> >
> > Which is very different from the previous results. So what am I doing
> > wrong here? Same issues arise with the other tests also (Pillai, Roy,
> > Hotelling-Lawley, etc).
> >
> > Many thanks and best wishes,
> > Ranjan
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > On Sun, 20 Mar 2011 19:29:41 -0500 John Fox <jfox at mcmaster.ca> wrote:
> >
> > > Dear Peter and Ranjan,
> > >
> > > In addition to Anova(), linearHypothesis() in the car package handles
> > > multivariate linear models, including those for repeated measures.
> > >
> > > Best,
> > > John
> > >
> > > --------------------------------
> > > John Fox
> > > Senator William McMaster
> > > Professor of Social Statistics
> > > Department of Sociology
> > > McMaster University
> > > Hamilton, Ontario, Canada
> > > http://socserv.mcmaster.ca/jfox
> > >
> > >
> > >
> > >
> > > > -----Original Message-----
> > > > From: r-help-bounces at r-project.org
> > > > [mailto:r-help-bounces at r-project.org]
> > > > On Behalf Of peter dalgaard
> > > > Sent: March-20-11 6:50 PM
> > > > To: Ranjan Maitra
> > > > Cc: R-help
> > > > Subject: Re: [R] manova question
> > > >
> > > >
> > > > On Mar 20, 2011, at 21:05 , Ranjan Maitra wrote:
> > > >
> > > > > Dear friends,
> > > > >
> > > > > Sorry for this somewhat generically titled posting but I had a
> > > > > question with using contrasts in a manova context. So here is my
> > > > question:
> > > > >
> > > > > Suppose I am interested in doing inference on \beta in the case of
> > > > > the model given by:
> > > > >
> > > > > Y = X %*% \beta + e
> > > > >
> > > > > where Y is a n x p matrix of observations, X is a n x m design
> > > > > matrix, \beta is m x p matrix of parameters, and e is a
> > > > > normally-distributed random matrix with mean zero and independent
> > > > > rows, each having dispersion matrix given by \Sigma. Then, I know
> > > > > (I think) how to perform MANOVA. Specifically, I use:
> > > > >
> > > > > fit <- manova(Y ~ X)
> > > > >
> > > > > and
> > > > >
> > > > > summary(fit) will allow me to perform appropriate inference on
> > beta.
> > > > >
> > > > > Now, suppose I am interested in doing inference on C %*% \beta %*%
> > > > > M (say testing whether this is equal to zero) with C and M being q
> > > > > x m and p x r matrices, respectively (with q, r both being no more
> > > > > than p), then can this be done using the manova object from the
> > above? How?
> > > > > If not, is there an efficient way to do this?
> > > >
> > > > Check out anova.mlm(), it does most of this sort of testing. Not
> > > > quite the "C %*% ..." bit because the linear model code is not
> > > > really built to handle linear constraints, but rather compare nested
> > > > models, each specified using a set of betas. (So you usually test
> > > > whether a subset of betas is zero).
> > > >
> > > > Also check out the "car" package. Its Anova() function does some
> > > > similar stuff.
> > > >
> > > > If noone has done so already, I wouldn't think it to be very hard to
> > > > implement the general case. Most of the bits are there already.
> > > >
> > > > --
> > > > Peter Dalgaard
> > > > Center for Statistics, Copenhagen Business School Solbjerg Plads 3,
> > > > 2000 Frederiksberg, Denmark
> > > > Phone: (+45)38153501
> > > > Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
> > > >
> > > > ______________________________________________
> > > > R-help at 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.
> > >
> >
> > ______________________________________________
> > R-help at 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.
>
More information about the R-help
mailing list