Two-Way ANOVA using R

Looking at the bottom plot, there appears to be an interaction between partner.status and fcategory. The difference in means in the two partner status levels is small when F-score category is high but larger when the F-score category is medium or low.

Next, we calculate our two-way ANOVA. To use type-III sum of squares in R, we cannot use the base R aov function. Instead, we fit the model using the lm function and then pipe the results into the Anova function from the car package. However, when using lm we have to carry out one extra step. By default, lm treats the factor levels as dummy variables. To correctly calculate the Type-III sum of squares, we need to use a different coding of the contrasts. We tell lm to use the appropriate sum to zero contrasts by adding an argument, contrasts = list(fcategory = contr.sum, partner.status = contr.sum).

mod <- lm(
  conformity ~ fcategory*partner.status, 
  data = Moore, 
  contrasts = list(fcategory = contr.sum, partner.status = contr.sum)
)
Anova(mod, type = "III")

## Anova Table (Type III tests)
## 
## Response: conformity
##                          Sum Sq Df  F value    Pr(>F)    
## (Intercept)              5752.8  1 274.3592 < 2.2e-16 ***
## fcategory                  36.0  2   0.8589  0.431492    
## partner.status            239.6  1  11.4250  0.001657 ** 
## fcategory:partner.status  175.5  2   4.1846  0.022572 *  
## Residuals                 817.8 39                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The output table gives:

• Sum Sq: The sum of squared errors for each model term
• Df: The degree of freedom
• F value: The F statistic for each variable
• Pr(>F): The p-value associated with each term (\(\lt .05\) indicates significance)

The first thing to investigate is the significance of the interaction. If it is non-significant, we can proceed to looking at the main effects. However, if the interaction is significant, the main effects will not be very helpful, as we will need to explore when each factor is significant given levels of the other factor. We indeed find that \(p = 0.023\) for partner.status*fcategory. This tells us that the effect of partner status will depend on levels of F-category, or vice versa. We will want to determine when each factor is significant.

The next step is to determine the nature of the interaction. For example, in the second interaction plot above we saw that the effect of partner.status appeared weak among the factory = high group, but it was larger for the other two fcategory levels. Perhaps partner.status has a significant effect on conformity, but only when fcategory is not high. To test this, we’ll need to perform something akin to a one-way ANOVA for partner.status within each level of fcategory. However, the correct F-test will utilize all of the information from the full two-way factorial model. To get this test, it is necessary to use the emmeans package (https://www.rdocumentation.org/packages/emmeans/versions/1.4.1). The following syntax calculates the means for each level of partner.status separately within each level of fcategory.

em_out_category <- emmeans(mod,  ~ partner.status | fcategory)
print(em_out_category)

## fcategory = low:
##  partner.status emmean   SE df lower.CL upper.CL
##  low              8.90 1.45 39     5.97     11.8
##  high            17.40 2.05 39    13.26     21.5
## 
## fcategory = medium:
##  partner.status emmean   SE df lower.CL upper.CL
##  low              7.25 2.29 39     2.62     11.9
##  high            14.27 1.38 39    11.48     17.1
## 
## fcategory = high:
##  partner.status emmean   SE df lower.CL upper.CL
##  low             12.62 1.62 39     9.35     15.9
##  high            11.86 1.73 39     8.36     15.4
## 
## Confidence level used: 0.95

Next, we pipe this object into pairs and use the joint=TRUE option in the test function to get the joint test of significance (i.e. the ANOVA F-test).

em_out_category %>%
  pairs() %>%
  test(joint = TRUE)

##  fcategory df1 df2 F.ratio p.value
##  low         1  39  11.486 0.0016 
##  medium      1  39   6.899 0.0123 
##  high        1  39   0.105 0.7477

The test function is akin to running multiple one-way ANOVAs, except we are using the error term from the full two-way factorial model. The results give:

• df1: The degrees of freedom between groups
• df2: The degrees of freedom within groups
• F.ratio: The F statistic for partner.status within each level of fcategory
• p.value: The probability of a type-I error for each category

We find that partner status is significant within the low (\(p = 0.002\)) and medium (\(p = 0.012\)) F-score category groups. The effect of partner status is not significant when F-category is high. This corroborates our theory from the interaction plots.

Because partner.status only has two levels, we do not need to do any further pairwise comparisons. If there were more than two levels, we could re-run the prior code chunk but remove the joint = TRUE option.

Now we will test if F-score category is significant in different levels of partner status. To do so, we reverse the variable order in the emmeans function:

em_out_status <- emmeans(mod,  ~ fcategory | partner.status)

em_out_status %>%
  pairs() %>%
  test(joint = TRUE)

##  partner.status df1 df2 F.ratio p.value note
##  low              2  39   2.323 0.1114   d  
##  high             2  39   2.138 0.1315   d  
## 
## d: df1 reduced due to linear dependence

We do not find a significant result, so no follow-up tests will be necessary. Note that, if fcategory had been significant in at least one level of partner.status, the following code would tell us which pairwise difference led to the significant joint F-test.

pairs(em_out_category)

## fcategory = low:
##  contrast   estimate   SE df t.ratio p.value
##  low - high   -8.500 2.51 39 -3.389  0.0016 
## 
## fcategory = medium:
##  contrast   estimate   SE df t.ratio p.value
##  low - high   -7.023 2.67 39 -2.627  0.0123 
## 
## fcategory = high:
##  contrast   estimate   SE df t.ratio p.value
##  low - high    0.768 2.37 39  0.324  0.7477