This tutorial is going to take the theory learned in our Two-Way ANOVA tutorial and walk through how to apply it using R. We will be using the `Moore`

dataset from the `carData`

package.
This data frame consists of subjects in a “social-psychological experiment who were faced with manipulated disagreement from a partner of either of low or high status. The subjects could either conform to the partner’s judgment or stick with their own judgment.” (John Fox, Sanford Weisberg and Brad Price (2018). carData: Companion to Applied Regression Data Sets. R package version 3.0-2. https://CRAN.R-project.org/package=carData).

Start by loading the packages we’ll use.

```
library(tidyverse)
library(carData)
library(car)
library(emmeans)
Moore <- Moore %>%
mutate(
partner.status = factor(
partner.status,
levels = c("low", "high")
),
fcategory = factor(
fcategory,
levels = c("low", "medium", "high")
)
)
```

Depending on the package used to perform two-way ANOVA, R will default to either Type-I or Type-II sum of squares. To be consistent with output from SPSS, we’ll use Type-III sum of squares in this tutorial. An in-depth summary of the types of sum of squares can be found in our previous post on the three different sum of squares.

We are interested in exploring an individuals conformity based on their partner’s status:

- Low Status
- High Status

We are also interested in exploring whether their F-score category (a measure of authoritarianism) affects outcomes or interacts with partner status. Our three treatment levels are:

- Low F-score
- Medium F-score
- High F-score

Our sample size is \(N = 45\). First, let’s inspect the data for outliers or funky distributions. The following boxplot shows the distribution of scores on the conformity variable within each combination of `partner.status`

and `fcategory`

.

```
Moore %>%
ggplot(aes(x = fcategory, y = conformity, color = partner.status)) +
geom_boxplot() +
xlab("F-Score Category") +
ylab("Conformity") +
scale_color_discrete(name = "Partner Status")
```

We can also get a sense of whether an interaction is present by looking at an interaction plot. An interaction plot shows the means for the outcome within each level of one of the factors, with separate lines for the other factor. Parallel lines indicate that no interaction is present, because the mean differences in the first factor are the same regardless of the level of the other factor. Non-parallel lines mean that an interaction is likely present. In other words, the mean differences on the first factor depend on the the level of the second factor. The `ggline`

function from the `ggpubr`

package (https://www.rdocumentation.org/packages/ggpubr/versions/0.2.3) makes nice interaction plots:

```
Moore %>%
ggpubr::ggline(
x = "fcategory",
y = "conformity",
color = "partner.status",
add = c("mean_se", "dotplot"),
palette = c("#00AFBB", "#E7B800"),
xlab = "F-Score Category",
ylab = "Conformity",
legend.title = "Partner Status"
)
```

```
Moore %>%
ggpubr::ggline(
x = "partner.status",
y = "conformity",
color = "fcategory",
add = c("mean_se", "dotplot"),
palette = c("#00AFBB", "#E7B800", "red"),
xlab = "Partner Status",
ylab = "Conformity",
legend.title = "F-Score Category")
```

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
```

Still have questions? Contact us!