Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

trouble with cbind in manova call in r

Tags:

r

I'm trying to do a multivariate ANOVA with the manova function in R. My problem is that I'm trying to find a way to pass the list of dependent variables without typing them all in manually, as there are many and they have horrible names. My data are in a data frame where "unit" is the dependent variable (factor), and the rest of the columns are various numeric response variables. e.g.

  unit C_pct  Cln C_N_mol Cnmolln C_P_mol N_P_mol
1    C 48.22 3.88   53.92    3.99 3104.75   68.42
2    C 49.91 3.91   56.32    4.03 3454.53   62.04
3    C 50.75 3.93   56.96    4.04 3922.01   69.16
4   SH 50.72 3.93   46.58    3.84 2590.16   57.12
5   SH 51.06 3.93   43.27    3.77 2326.04   53.97
6   SH 48.62 3.88   40.97    3.71 2357.16   59.67

If I write the manova call as

fit <- manova(cbind(C_pct, Cln) ~ unit, data = plots) 

it works fine, but I'd like to be able to pass a long list of columns without naming them one by one, something like

fit <- manova(cbind(colnames(plots[5:32])) ~ unit, data = plots) 

or

fit <- manove(cbind(plots[,5:32]) ~ unit, data = plots)

I get the error

"Error in model.frame.default(formula = as.matrix(cbind(colnames(plots[5:32]))) ~  : 
  variable lengths differ (found for 'unit')

I'm sure it's because I'm using cbind wrong, but can't figure it out. Any help is appreciated! Sorry if the formatting is rough, this is my first question posted.

EDIT: Both ways (all 3, actually) work. thanks all!

like image 874
Iceberg Slim Avatar asked Jun 18 '13 17:06

Iceberg Slim


People also ask

How to interpret MANOVA test in R?

Interpretation of MANOVAIf the global multivariate test is significant, we conclude that the corresponding effect (treatment) is significant. In that case, the next question is to determine if the treatment affects only the weight, only the height or both.

What does Manova do in R?

MANOVA in R uses Pillai's Trace test for the calculations, which is then converted to an F-statistic when we want to check the significance of the group mean differences. You can use other tests, such as Wilk's Lambda, Roy's Largest Root, or Hotelling-Lawley's test, but Pillai's Trace test is the most powerful one.


2 Answers

manova, like most R modelling functions, builds its formula out of the names of the variables found in the dataset. However, when you pass it the colnames, you're technically passing the strings that represent the names of those variables. Hence the function doesn't know what to do with them, and chokes.

You can actually get around this. The LHS of the formula only has to resolve to a matrix; the use of cbind(C_pct, Cln, ...) is a way of obtaining a matrix by evaluating the names of its arguments C_pct, Cln, etc in the environment of your data frame. But if you provide a matrix to start with, then no evaluation is necessary.

fit <- manova(as.matrix(plots[, 5:32]) ~ unit, data=plots)

Some notes. The as.matrix is necessary because getting columns from a data frame like this, returns a data frame. manova won't like this, so we coerce the data frame to a matrix. Second, this works assuming you don't have an actual variable called plots inside your data frame plots. This is because, if R doesn't find a name inside your data frame, it then looks in the environment of the caller, in this case the global environment.

You can also create the matrix before fitting the model, with

plots$response <- as.matrix(plots[, 5:32])
fit <- manova(response ~ unit, data=plots)
like image 180
Hong Ooi Avatar answered Sep 22 '22 21:09

Hong Ooi


You can build your formula as string and cast it to a formula:

responses <- paste( colnames( plots )[2:6], collapse=",")
myformula <- as.formula( paste0( "cbind(", responses , ")~ unit" ) )
manova( myformula, data = plots ) 

Call:
   manova(myformula, data = plots)

Terms:
                     unit Residuals
resp 1                0.4       6.8
resp 2                  0         0
resp 3              220.6      21.0
resp 4                0.1       0.0
resp 5          1715135.8  377938.1
Deg. of Freedom         1         4

Residual standard error: 1.3051760.027080132.293640.04966555307.3834
Estimated effects may be unbalanced
like image 26
Beasterfield Avatar answered Sep 19 '22 21:09

Beasterfield