0
$\begingroup$

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
$\endgroup$

1 Answer 1

0
$\begingroup$

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.

Now I will try to answer your 3 questions in reverse order.

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?

I personally would suggest that you consult your local statistical expert. But if you really must do it yourself, then perhaps read up on mixed model?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.