# ANOVA in R: How MSresiduals are estimated when using Error(interaction)?

I have checked the existing answers, but I found no answer to my question. It is about choosing the right ANOVA model in R and how the MSresiduals are estimated.

I have an experiment with two factors: MANAGEMENT and REGION. I have three MANAGEMENT types and three REGIONS. That means nine sites under study but no repetition, like this:

REGION: A, B and C (those regions are not clonal copies of each other)
MANAGEMENT: N, P, ST
the dependent variable is called PHt

I have set the aov in R in two possible ways.

Model1: aov(PHt~REGION+MANAGEMENT+REGION*MANAGEMENT, data=data)

Model2: aov(PHt~REGION*MANAGEMENT+Error(subject/(REGION*MANAGEMENT)),data=data)

In respect to these two models I have three questions:

1. Is it correct that the second model should be used if we aim to include the interaction term into the Error term?
2. How does the model 2 estimate the MSresiduals? It seems to use different MSresiduals for each variable. Why? (See output for both models below.)
3. Since I have unbalanced data, is R using an approximation of a Wald test?

Here are results of the two models.

Model 1
Model <- aov(PHt~REGION+MANAGEMENT+REGION*MANAGEMENT, data=data)
summary(model)
Df      Sum Sq        Mean Sq    F value       Pr(>F)
REGION               2      0.05          0.0229       0.907        0.4039
MANAGEMENT           2      0.18          0.0884       3.494        0.0306 *
REGION:MANAGEMENT    4      3.34          0.8351      33.021        <2e-16 ***
Residuals           1744    44.10         0.0253

Model 2
aov(PHt~REGION*MANAGEMENT+Error(subject/(REGION*MANAGEMENT)),data=data)

SUMMARY:
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 195  4.714 0.02494

Error: subject_f:REGION
Df Sum Sq Mean Sq F value Pr(>F)
REGION      2  0.066 0.03281   1.406  0.246
Residuals 388  8.822 0.02334

Error: subject_f:MANAGEMENT
Df Sum Sq Mean Sq F value Pr(>F)
MANAGEMENT   2  0.205 0.10260   4.413 0.0127 *
Residuals  388  8.788 0.02325
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: subject_f:REGION:MANAGEMENT
Df Sum Sq Mean Sq F value Pr(>F)
REGION:MANAGEMENT   4  3.451  0.8627   31.99 <2e-16 ***
Residuals         776 20.389  0.0270
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


To start with, PHt is your dependent variable, not independent variable. In addition, your model 1 can simply be written as aov(PHt~REGION*MANAGEMENT,data=data) because the * symbol denotes both main effects and interaction. If you want interaction only, it should be REGION:MANAGEMENT. See ?formula for more detail.
A3: The fact that the design is unbalanced really suggests you should be using mixed model (unless it is very close to being balanced, then perhaps you can stay with aov) Google the R package nlme for fitting mixed model.
A2: The MS for Residuals in each error stratum (subject_f:REGION, then subject_f:MANAGEMENT, then subject_f:REGION:MANAGEMENT) is trying to estimate the variability attributed to that specific stratum. So the first MSResidual, 0.02334, is saying subject_f:REGION adds $\sigma^2\approx 0.02334^2$ to the total variability in the data, while subject_f:MANAGEMENT adds approximately $0.02325^2$ to the total variation in data. Since they are different factors that add to the total variability, it makes sense that they are different.
A1: I am actually not very sure what you are asking. I guess the surface answer is yes, if you want to include the interaction term in the Error() term, then it is correct to 'include the interaction term in the Error() term'. However, is this really what you want? It is not immediately clear to me that MANAGEMENT should appear in the Error() term. I am not saying it can't be the case, but are you sure that's what you want? Also, what is the interpretation of subject/MANAGEMENT*REGION?