Output of aov() in R varies with changing order of independent variables I need to understand why does the summary output of aov() in R vary when the order of independent variable changes. Like for example:
summary(aov(y~A+B, data)) is different from summary(aov(y~B+A,data)).
 A: As @SmallChess pointed out, 

"The anova and aov functions in R implement a sequential sum of
  squares (type I)." [source]

This procedure has probably been tackled already by different questions in the site. Basically, it works like this [source]

The different types of sums of squares then arise depending on the
  stage of model reduction at which they are carried out. In particular:
  
  
*
  
*Type I, also called "sequential" sum of squares:
  
  
*
  
*SS(A) for factor A.
  
*SS(B | A) for factor B.
  
*SS(AB | B, A) for interaction AB.
  
*This tests the main effect of factor A, followed by the main effect of
  factor B after the main effect of A, followed by the interaction
  effect AB after the main effects.
  
*Because of the sequential nature and the fact that the two main factors are tested in a particular order, this type of sums of squares
  will give different results for unbalanced data depending on which
  main effect is considered first.

So the the effect of B in summary(aov(y~A+B, data)) is tested in the presence of A, while the effect of B in summary(aov(y~B+A, data)) is tested in the absence of A.
A: Very likely your design is unbalanced. Try this function !is.list(replications(y~A+B, data)) to see if this is so. If it is unbalanced, then you might try the lm function from the nlme package. 
See this answer for more in depth information regarding decisions around unbalanced designs. 
https://stats.stackexchange.com/questions/13241/the-order-of-variables-in-anova-matters-doesnt-it
