# manova question

6 messages
Open this post in threaded view
|
Report Content as Inappropriate

## manova question

 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: 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: [hidden email]  Priv: [hidden email] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: manova question

Open this post in threaded view
|
Report Content as Inappropriate

## Re: manova question

Open this post in threaded view
|
Report Content as Inappropriate

## Re: manova question

 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: [hidden email] [mailto:[hidden email]] > 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 <[hidden email]> 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: [hidden email] > > > [mailto:[hidden email]] > > > 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: [hidden email]  Priv: [hidden email] > > > > > > ______________________________________________ > > > [hidden email] 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. > > > > ______________________________________________ > [hidden email] 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. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
 In reply to this post by Ranjan Maitra 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 <[hidden email]> 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: [hidden email] [mailto:[hidden email]] > > 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 <[hidden email]> 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: [hidden email] > > > > [mailto:[hidden email]] > > > > 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: [hidden email]  Priv: [hidden email] > > > > > > > > ______________________________________________ > > > > [hidden email] 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. > > > > > > > ______________________________________________ > > [hidden email] 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. > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.