8 Workshop 8: Factorial ANOVA: Mixed Design Factorial ANOVA

Aims:

Practice importing and wrangling data ready for analysis.
Practice running a mixed design factorial ANOVA, including interpreting the model and running any follow-up analyses.
Practice evaluating model assumptions.
Practice generating plots and writing-up the analysis.

Introduction

This week’s workshop is going to bring together several things you have covered so far during PS2010, including right back to the week one workshop. Today you are going to analyse a data set starting from a “messy” data set right the way through to producing some fancy data visualisations to show off your findings.

All of the steps covered in today’s workshop could be helpful for your lab report analysis and write-up.

There are quite a few steps today and so don’t worry if you cannot complete them all in the allotted time, but instead save your script and use your independent study time to complete it.

Today’s Research Scenario

There are lots of tasks where our performance can be negatively impacted by distractions. There are some tasks whereby we should solely focus on what we are doing, for example, when completing this workshop or driving a car.

If you are distracted completing this workshop, the worst that could happen is that you make a mistake with your code. Distractions when driving can have much more serious consequences. Today’s research question is: How do distractions impact our ability to detect different types of hazards when driving.

Sixty experienced drivers (minimum 1 year driving experience) were recruited to take part in a driving simulation task, whereby they had to identify different hazards. Participants were first randomly allocated to one of two conditions.

Independent variable 1 (distraction yes/no):
- Yes: participants were required to answer questions by the researcher via a hands free headset whilst in the simulator.
- No: participants were left to drive uninterrupted in the simulator.

Question 1: Was this variable independent or repeated measures?

During the driving simulation task, participants had to identify hazards (e.g., a car pulling out in front, or a pedestrian stepping out) by pressing a button. There were three different types of hazard that participants saw.

Independent variable 2 (hazard type: car, bike, pedestrian)
- Car: a car was the cause of the hazard.
- Bike: a bicycle was the cause of the hazard
- Pedestrian: a pedestrian was the cause of the hazard.

Question 2: Was this variable independent or repeated measures?

The dependent variable in this study was the number of hazards correctly identified. The driving task was split into two blocks, and there were a maximum of 10 hazards for each hazard type in each block (total 30 hazards per block).

Workshop Outline

The workshop is split into different parts:

Exercises 1-2: Import and Prepare Data for Analysis This will involve setting up RStudio, importing the data, and any data wrangling tasks to ensure the data is in the correct format prior to analysis.

Exercise 3: Calculate Statistics This will involve running initial descriptive statistics and the initial ANOVA model.

Exercises 4-6 Carrying out any additional analyses (e.g., examine the interaction or interpret main effects) depending on the outcome of the initial ANOVA model.

Exercise 7 Evaluate the model. This will be checking any necessary assumptions.

Exercise 8 Generate publication worthy figures and write-up the results.

Let’s begin…

8.1 Exercise 1: Prepare RStudio and Import the Data

You’ll need these packages:

library(tidyverse)
library(afex)
library(car)
library(emmeans)
library(effectsize)

Import the data set from the csv file called hazard.csv.

View the data to familiarise yourself with it:

id: the identifier.
age: participant age in years.
exp_code: was a code used to log in to the simulator experiment system for that participant.
years_driving: how many years they have been driving.
distraction: which distraction condition they were randomly assigned to.
block1_car and block2_car: the number of car hazards identified in blocks 1 and 2 respectively.
block1_bike and block2_bike: the number of bike hazards identified in blocks 1 and 2 respectively.
block1_pedestrian and block2_pedestrian: the number of pedestrian hazards identified in blocks 1 and 2 respectively.

8.2 Exercise 2: Data Wrangling

The data set is a bit busier than what you have seen before, and requires some work before it is ready to analyse. This is a normal part of data analysis.

You need to calculate the total number of hazards correctly identified, because these are split over two blocks at the moment.

Question 3: Which R function could help with this?

There are a few variables in the data set we don’t really need, you could get rid of them to clean up the data set if you want.

Question 4: Which R function could help with this?

The data has a repeated measures variable and this is in wide format currently. You will need long format data.

Question 5: Which R function could help with this?

Calculate the Total Scores

mydata <- mydata %>%
  mutate(car = block1_car + block2_car,
         bike = block1_bike + block2_bike,
         pedestrian = block1_pedestrian + block2_pedestrian)

This code will create three new variables: car, bike, and pedestrian to contain the sum of the two blocks. You used the mutate() function in week 1’s workshop.

view(mydata)

Can you spot the new variables that you created?

Select Relevant Variables

Use names(mydata) for a list of the variables in the data set. Use the code below to select the variables you will need for your ANOVA analysis. You do not need everything in the original mydata data set.

For example, we really don’t need the original block1_ or block2_ scores now we have the totals from above, and similarly, the exp_code is not needed for our analysis.

mydata_clean <- mydata %>%
  select(NULL, NULL, NULL, NULL) # add as many as necessary

Change each NULL to a variable you want to keep (each separated by a comma). You used the select() function in week 1’s workshop. Note: this code creates a new object called mydata_clean so that we don’t overwrite the original mydata object.

Use the hint below if needed.

👀 Click for a hint

mydata_clean <- mydata %>%
  select(id, age, years_driving, distraction, car, bike, pedestrian)

view(mydata_clean)

Do you notice that there are fewer variables in this data set? That looks a bit tidier!

Convert from Wide to Long format

longd <- mydata_clean %>%
  pivot_longer(
    cols = car:pedestrian,
    names_to = "hazard",
    values_to = "acc")

This will create two new labels: hazard for type of hazard, and acc for number of hazards correctly identified (accuracy). From now on you should use the longd object.

Re-code Any Factors

Last bit of wrangling. Make sure you have any factors coded as a Factor and not a character.

longd$distraction <- factor(longd$distraction)
longd$hazard <- factor(longd$hazard)

Now you are ready for the real analysis.

8.3 Exercise 3: Descriptive Stats and the Model

Calculate descriptive statistics just like you have done in previous weeks.

Use names(longd) if you need a reminder of the variable names.

desc <- longd %>%
  group_by(NULL, NULL) %>%
    summarise(mean_acc = mean(NULL),
              sd_acc = sd(NULL))
view(desc)

You could also generate a box plot if you like (use code from previous weeks) but we will generate some visualisations later on.

Now to run the initial ANOVA model…

mod <- aov_ez(id = "NULL",
              dv = "NULL",
              within = "NULL",
              between = "NULL",
              type = 3,
              include_aov = TRUE,
              data = NULL)

Take extra care to ensure you enter the correct variable for the within and between lines of code. Look back over the study scenario section at the start of the workshop if you need.

Sphericity Check

As one of the independent variables is repeated measures we need to run an extra check here.

aov_ez() will automatically check Sphericity for you when you have a repeated measures variable with three or more levels, and it will essentially produce two ANOVA models: a corrected model and an uncorrected model. You need to decide which one to use.

If the assumption of Sphericity is violated, you should use the corrected model. If the assumption of Sphericity is met, you can use the uncorrected model.

To check, produce the model output but this time using summary(mod) as this gives us some extra information.

summary(mod)

Check the output and look for Mauchly Tests for Sphericity:

Mauchly Tests for Sphericity

                   Test statistic p-value
hazard                    0.96285 0.33993
distraction:hazard        0.96285 0.33993

To report the assumption based on the output above, i’d write something like: > Mauchly’s test indicated that the assumption of sphericity was not violated for the main effect of hazard, w = .96, p = .340 or the distraction × hazard interaction, W = .96 , p = .340.

As these p-values were not-significant this means that Sphericity was NOT violated and the assumption was met. This now means we can report the ANOVA model without any correction.

If the assumption is violated (either p-value is less than 0.05 and therefore significant) then you should use and report the corrected model. To find this just use:

nice(mod)

If you run this now to compare, you’ll notice that the degrees of freedon are slightly adjusted (they’re not whole numbers) and some other values, such as p-values, might have changed slightly. If you see Sphericity correction method: GG this meand a Greenhouse-Geisser correction has been applied and you should include that information in your write-up.

As for today’s scenario the assumption was met, continue wit the uncorrected model.

Interpret the ANOVA model and decide what should happen next…

The first main effect (distraction) is

The second main effect (hazard type) is

The distraction by hazard type interaction is

8.4 Exercise 4: Interpreting the Main effect of Distraction

As there are only two levels, you can just use the mean values to interpret the direction. You can also use a box plot too.

m_distract <- longd %>%
  group_by(NULL) %>%
    summarise(mean_acc = mean(NULL),
              sd_acc = sd(NULL))
view(m_distract)

ggplot(longd, aes(x = distraction, y = acc)) +
  geom_boxplot() +
  theme_classic()

Question 9: Which statement is the correct interpretation of this main effect?

8.5 Exercise 5: Interpreting the Main Effect of Hazard Type

As there are more than two levels for this main effect, we cannot simply look at the means just yet, as it won’t be clear which conditions are significantly different. All the main tell us is that at least one of them are significantly different, but which one?

Here you will need to follow the same steps as a one-way ANOVA (see Workshops 4 & 5 for a refresher). As the main effect was significant, in this scenario you will run post-hoc multiple comparisons to see how they differ. Work through the code below to figure out how to best interpret this main effect using:

pairwise comparisons
box plot for visualising data
means and standard deviation descriptive statistics

emmeans(mod,
        pairwise ~ hazard,
        adjust = "holm")

ggplot(longd, aes(x = hazard, y = acc)) +
  geom_boxplot() +
  theme_classic()

m_hazard <- longd %>%
  group_by(hazard) %>%
  summarise(mean_acc = mean(acc),
            sd_acc = sd(acc))
view(m_hazard)

You have used the emmeans() function before in week 5’s workshop. Interpret the output.

Question 10: Which was the only pairwise comparison that was NOT significant?

Question 11: Using all available information, which statement is NOT a correct interpretation of this main effect?

Question 12: Using all available information, which statement is the correct interpretation of this main effect?

8.6 Exercise 6: Effect Sizes

So far you have the initial ANOVA model which gave you:

F, df, and p for each main effect and the interaction.

Descriptive statistics for main effect of distraction which gave you:

Mean and standard deviation of identified hazard for the two groups (yes, no for presence of distraction)

Descriptive statistics and pairwise comparisons for main effect of hazard type which gave you:

Mean and standard deviation of identified hazards for the three conditions (car, bike, pedestrian for hazard type) and t, df, and p from pairwise comparisons between the three.

But what about effect sizes?

For the ANOVA model, you should look at partial eta squared (\(η_P^2\)). This is super quick and easy to check from the ANOVA model mod using this code:

eta_squared(mod)

You can report the full F-statement in APA format like this:

F(df1, df2) = value, p = value, \(η_P^2\) = value

Anything underlined needs to be replaced by values from your output so far.

For the pairwise comparisons, it is a bit tricky as it depends if the variable for the main effect is an independent (between) or repeated (within) measures. If it is independent measures then follow the code from last week’s workshop (see cohens_d).

As we have a repeated design we can use the same code as the t-test workshop earlier in term. **One annoying this here is we need to go back and use the original mydata data as this is a rare occasion that wide data is preferred (how annoying, right?!)

cohens_d(mydata$car, mydata$bike, paired = TRUE)
cohens_d(mydata$car, mydata$pedestrian, paired = TRUE)
cohens_d(mydata$pedestrian, mydata$bike, paired = TRUE)

These need to match the three comparisons carried out in exercise 4.

Question 13: Which pairwise comparison had the largest effect size?

8.7 Exercise 7: Evaluating Assumptions

We Sphericity was already checker earlier on. As this is a mixed design, we now have two more assumptions to check.

8.7.1 Checking Normality

You have done this before. Same code, providing your ANOVA model is called mod.

# testing normality, looking at normal distribution of residuals from the model.
residuals <- residuals(mod) #pulls residuals from the model
qqPlot(residuals) #produces a qq-plot of residuals
hist(residuals, breaks = 20, main = "Histogram of Residuals", xlab = "Residuals", col = "lightblue") #histogram of residuals
shapiro.test(residuals) #Shapiro test for residuals

8.7.2 Checking Homogeneity of Variance

You should use Levene’s test. You’ve also done this before but for a one-way independent ANOVA.

leveneTest(NULL1 ~ NULL2, longd)

Change NULL1 to the dependent variable.
Change NULL2 to the independent (between subjects) measures independent variable.
The third item is the object which contains the data, in this case longd again.

Remember, a non-significant p-value (>.05) means the assumption was met.

8.8 Exercise 8: Data Visualisation

When considering data visualisation options, you should consider “what are the significant and important findings I want to show my reader?”. Given you have two significant main effects and a non-significant interaction, if you only present a figure showing the interaction, you’re not teling the full story.

Use the code below to generate three plots. Take a look at them and interpret what they show.

8.8.1 Plot 1

📊 Click for Code

p1 <- ggplot(longd, aes(x = distraction, y = acc)) +
  stat_summary(fun = mean, geom = "bar", 
               position = position_dodge(), 
               color = "black",
               fill = c("slategray4", "slategray2")) +
  stat_summary(fun.data = mean_se, geom = "errorbar", 
               position = position_dodge(0.9), width = 0.25) +
  labs(
    title = NULL,
    x = "Distraction",
    y = "Mean Number of Hazards Identified") +
  scale_x_discrete(labels = c("yes" = "Yes", "no" = "No")) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold"))

8.8.2 Plot 2

📊 Click for Code

p2 <- ggplot(longd, aes(x = hazard, y = acc)) +
  stat_summary(fun = mean, geom = "bar", 
               position = position_dodge(), 
               color = "black",
               fill = c("slategray4", "slategray2", "darkslategray3")) +
  stat_summary(fun.data = mean_se, geom = "errorbar", 
               position = position_dodge(0.9), width = 0.25) +
  labs(
    title = NULL,
    x = "Hazard Type",
    y = "Mean Hazards Identified") +
  scale_x_discrete(labels = c("car" = "Car", "bike" = "Bike", "pedestrian" = "Pedestrian")) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold"))

8.8.3 Plot 3

📊 Click for Code

p3 <- ggplot(longd, aes(x = hazard, y = acc, fill = distraction)) +
  stat_summary(fun = mean, geom = "bar", position = position_dodge(), color = "black") +
  stat_summary(fun.data = mean_se, geom = "errorbar", 
               position = position_dodge(0.9), width = 0.25) +
  labs(
    title = NULL,   # leave this as NULL because APA figures do not use titles.
    x = "Hazard",
    y = "Mean Number of Hazards Identified",
    fill = "Distraction") +
  scale_fill_manual(values = c("yes" = "slategray4", "no" = "slategray2"), 
                    labels = c("Yes", "No")) +
  scale_x_discrete(labels = c("car" = "Car", "bike" = "Bike", "pedestrian" = "Pedestrian")) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.title = element_text(face = "bold"),
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold"))

Once you have run all three, use the code below to present all three at once. You’ll need a new package.

install.packages("gridExtra")
library(gridExtra)

grid.arrange(p1,p2,p3, ncol = 3)

NOTE: You will likely need to resize the output panel (bottom right panel) by making it wider so all three plots are visible.

Question 14: Which plot tells us about the main effect of distraction?

Question 15: Which plot tells us about the main effect of hazard?

Question 16: Which plot tells us about the interaction?

Well Done. You have reached the end of the workshop

7 Workshop 7: Factorial ANOVA: Different Designs and Assumptions

9 Workshop 9: Non-Parametric Statistics