Why do the results of a MANOVA change when the order of the predictor variables is changed? So, for example, using the Iris data and treating iris species as the predictor variable and sepal length, sepal width, petal length, and petal width as the dependent variables we get MANOVA output that looks like this:
set.seed(2)
# Creating a matrix of the 4 dependent variables (DVs)
Y <- as.matrix(iris[,c(1:4)]) 

# MANOVA looking at the effect of species on DVs
summary(manova(Y ~ iris$Species))
    #                 Df  Pillai    approx F  num Df  den Df    Pr(>F)    
    # iris$Species    2   1.1919    53.466    8       290       < 2.2e-16 ***

That seems to make sense. Species has a significant effect on our DVs. Now what if we add another predictor variable (a random one which we shouldn’t expect to have an effect on the DVs)?
# Creating a random dummy variable to be used as a predictor variable
iris$random.dummy <- sample(x = c(0,1), size = 150, replace = TRUE)

# MANOVA looking at the effect of species + our random dummy on DVs
summary(manova(Y ~ iris$Species + iris$random.dummy))
    #                     Df  Pillai    approx F  num Df  den Df  Pr(>F)    
    # iris$Species        2   1.19339   53.263    8       288     <2e-16 ***
    # iris$random.dummy   1   0.03784   1.406     4       143     0.2349  

That also seems to make sense. Species is significant still, but our random dummy variable is not. Now what if we simply switch the order of those variables?
# Switching the order of our two predictor variables in the formula
summary(manova(Y ~ iris$random.dummy + iris$Species))
    #                     Df  Pillai    approx F  num Df  den Df  Pr(>F)  
    # iris$random.dummy   1   0.13031   5.357     4       143     0.0004764 ***
    # iris$Species        2   1.19526   53.470    8       288     < 2.2e-16 ***

Now, the Pillai’s trace and approximate F-values change and our random dummy variable has become significant. 
So my questions are these.
Why do the results of a MANOVA change when the order of the predictor variables is changed?
and
What does this mean for those of us trying to use and interpret a MANOVA?
 A: I read this some weeks ago:

MANOVA uses a regression approach to analysis of variance. More particularly, the program employs a hierarchical model. There is an important consequence for the user: if a MANOVA execution involves more than 1 factor variable, and if there are disproportionate number of cases in the cells formed by the cross-classification of the factors, then consideration must be given to the order in which factor variables are specified. Disproportionality of subclass numbers confounds the main effects and the researcher must choose the order in which the confounded effects should be eliminated. When using MANOVA, this choice is accomplished by the order in which factor variables are specified. When using standard ordering, variables early in the specification have the effects of later variables removed, e.g. the first listed effect will be tested with all other main effects eliminated. The general rule is that each test eliminates effects listed before it on the test name specifications and ignores effects listed afterward. For a standard two-way analysis, the interaction term is not affected by the order of factor variables; more generally, for a standard n-way analysis, the n-th order interaction term and that term only, is unaffected. The problem exists for both univariate and multivariate analysis. 

