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? Many thanks in advance for all your help, and best wishes, Ranjan
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
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.
>
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.
>