Running a One-Way & Two-Way ANOVA Lab
Charlotte Wiser
2026-06-15
1 Loading Libraries
#install.packages("afex")
#install.packages("emmeans")
#install.packages("ggbeeswarm")
#install.packages("expss")
library(psych) # for the describe() command
library(ggplot2) # to visualize our results
library(expss) # for the cross_cases() command
library(car) # for the leveneTest() command
library(afex) # to run the ANOVA
library(ggbeeswarm) # to run plot results
library(emmeans) # for posthoc tests2 Importing Data
# For HW, import the project dataset you cleaned previously this will be the dataset you'll use throughout the rest of the semester
d <- read.csv(file="Data/projectdata.csv", header=T)
# new code! this adds a column with a number for each row. It will make it easier if we need to drop outliers later
d$row_id<- 1:nrow(d)Note: For your HW, you will choose to run EITHER a one-way ANOVA (a single IV with 3 or more levels) OR a two-way/factorial ANOVA (two IVs with 2 or 3 levels each). You will need to specify your hypothesis and customize your code based on the choice you make. We will run BOTH versions of the test in the lab for illustrative purposes.
3 One-Way ANOVA
3.1 State Your Hypothesis
We predict that there will be a significant difference in people’s openness to trying something new based on people’s level of mental health disorders (eating disorder, an anxiety disorder, no mental health disorder).
3.2 Check Your Variables
# you only need to check the variables you're using in the current analysis
str(d)## 'data.frame': 684 obs. of 8 variables:
## $ X : int 520 2814 3146 3295 717 6056 4753 5365 2044 1965 ...
## $ relationship_status: chr "Single, never married" "Single, never married" "Prefer not to say" "Single, never married" ...
## $ mhealth : chr "none or NA" "none or NA" "none or NA" "none or NA" ...
## $ pswq : num 2.71 1.43 1.86 1.79 2.36 ...
## $ mfq_26 : num 2.7 4.55 4.8 3.8 4.5 4 5.8 4.2 4.5 5.25 ...
## $ rse : num 2.6 3.1 3.7 3 3 3 4 3.8 2.5 4 ...
## $ phq : num 1.56 1.44 1.11 1.33 1.44 ...
## $ row_id : int 1 2 3 4 5 6 7 8 9 10 ...# make our categorical variables of interest "factors"
# because we'll use our newly created row ID variable for this analysis, so make sure it's coded as a factor, too.
d$mhealth <- as.factor(d$mhealth)
d$row_id <- as.factor(d$row_id)
d<- subset(d, mhealth != "depression") # use subset() to remove all participants from the additional level
table(d$mhealth, useNA = "always") # verify that now there are ZERO participants in the additional level##
## anxiety disorder bipolar
## 80 3
## depression eating disorders
## 0 18
## none or NA obsessive compulsive disorder
## 528 16
## other ptsd
## 17 12
## <NA>
## 0 d$mhealth <- droplevels(d$mhealth) # use droplevels() to drop the empty factor
table(d$mhealth, useNA = "always") # verify that now the entire factor level is removed ##
## anxiety disorder bipolar
## 80 3
## eating disorders none or NA
## 18 528
## obsessive compulsive disorder other
## 16 17
## ptsd <NA>
## 12 0d<- subset(d, mhealth != "bipolar") # use subset() to remove all participants from the additional level
table(d$mhealth, useNA = "always") # verify that now there are ZERO participants in the additional level##
## anxiety disorder bipolar
## 80 0
## eating disorders none or NA
## 18 528
## obsessive compulsive disorder other
## 16 17
## ptsd <NA>
## 12 0 d$mhealth <- droplevels(d$mhealth) # use droplevels() to drop the empty factor
table(d$mhealth, useNA = "always") # verify that now the entire factor level is removed ##
## anxiety disorder eating disorders
## 80 18
## none or NA obsessive compulsive disorder
## 528 16
## other ptsd
## 17 12
## <NA>
## 0d<- subset(d, mhealth != "obsessive compulsive disorder") # use subset() to remove all participants from the additional level
table(d$mhealth, useNA = "always") # verify that now there are ZERO participants in the additional level##
## anxiety disorder eating disorders
## 80 18
## none or NA obsessive compulsive disorder
## 528 0
## other ptsd
## 17 12
## <NA>
## 0 d$mhealth <- droplevels(d$mhealth) # use droplevels() to drop the empty factor
table(d$mhealth, useNA = "always") # verify that now the entire factor level is removed ##
## anxiety disorder eating disorders none or NA other
## 80 18 528 17
## ptsd <NA>
## 12 0d<- subset(d, mhealth != "ptsd") # use subset() to remove all participants from the additional level
table(d$mhealth, useNA = "always") # verify that now there are ZERO participants in the additional level##
## anxiety disorder eating disorders none or NA other
## 80 18 528 17
## ptsd <NA>
## 0 0 d$mhealth <- droplevels(d$mhealth) # use droplevels() to drop the empty factor
table(d$mhealth, useNA = "always") # verify that now the entire factor level is removed ##
## anxiety disorder eating disorders none or NA other
## 80 18 528 17
## <NA>
## 0d<- subset(d, mhealth != "other") # use subset() to remove all participants from the additional level
table(d$mhealth, useNA = "always") # verify that now there are ZERO participants in the additional level##
## anxiety disorder eating disorders none or NA other
## 80 18 528 0
## <NA>
## 0 d$mhealth <- droplevels(d$mhealth) # use droplevels() to drop the empty factor
table(d$mhealth, useNA = "always") # verify that now the entire factor level is removed ##
## anxiety disorder eating disorders none or NA <NA>
## 80 18 528 0# check that our categorical variables of interest are now factors
str(d)## 'data.frame': 626 obs. of 8 variables:
## $ X : int 520 2814 3146 3295 717 6056 4753 5365 2044 1965 ...
## $ relationship_status: chr "Single, never married" "Single, never married" "Prefer not to say" "Single, never married" ...
## $ mhealth : Factor w/ 3 levels "anxiety disorder",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ pswq : num 2.71 1.43 1.86 1.79 2.36 ...
## $ mfq_26 : num 2.7 4.55 4.8 3.8 4.5 4 5.8 4.2 4.5 5.25 ...
## $ rse : num 2.6 3.1 3.7 3 3 3 4 3.8 2.5 4 ...
## $ phq : num 1.56 1.44 1.11 1.33 1.44 ...
## $ row_id : Factor w/ 684 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...# check our DV skew and kurtosis
describe(d$mfq_26)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 626 4.27 0.64 4.3 4.28 0.59 1.8 5.8 4 -0.28 0.12 0.03# we'll use the describeBy() command to view our DV's skew and kurtosis across our IVs' levels
describeBy(d$mfq_26, group = d$mhealth)##
## Descriptive statistics by group
## group: anxiety disorder
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 80 3.96 0.71 3.88 3.94 0.74 2.4 5.5 3.1 0.17 -0.56 0.08
## ------------------------------------------------------------
## group: eating disorders
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 18 4.02 0.65 4.07 4.03 0.63 2.8 5.05 2.25 -0.31 -1.05 0.15
## ------------------------------------------------------------
## group: none or NA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 528 4.32 0.61 4.35 4.33 0.59 1.8 5.8 4 -0.3 0.39 0.03# also use histograms to examine your continuous variable
hist(d$mfq_26)
3.3 Check Your Assumptions
3.3.1 ANOVA Assumptions
- DV should be normally distributed across levels of the IV (we checked previously using “describeBy” function)
- All levels of the IVs should have an equal number of cases and there should be no empty cells. Cells with low numbers decrease the power of the test (which increases chance of Type II error)
- Homogeneity of variance should be confirmed (using Levene’s Test)
- Outliers should be identified and removed – we will actually remove them this time!
- If you have confirmed everything above, the sampling distribution should be normal.
3.3.2 Check levels of IVs
table(d$mhealth)##
## anxiety disorder eating disorders none or NA
## 80 18 528# REMEMBER your test's level of POWER is determined by your SMALLEST subsample3.3.3 Check for outliers using Cook’s distance and Residuals VS Leverage plot
3.3.3.1 Run a Regression to get both outlier plots
# use this commented out section below ONLY IF if you need to remove outliers
# to drop a single outlier, use this code:
#d <- subset(d, row_id!=c(1108)) # NOTE: Participant 1108 is not truly problaemtic in the Lab, but we are going to removed them for demonstration.
# to drop multiple outliers, use this code:
# d <- subset(d, row_id!=c(1108) & row_id!=c(602))
# use the lm() command to run the regression
# formula is y~x, where y is our DV, x is our IV
reg_model <- lm(mfq_26~mhealth, data = d) 3.3.3.2 Check for outliers
# Cook's distance
plot(reg_model, 4)
# Residuals VS Leverage
plot(reg_model, 5)
# IF you find outliers, go back up to line 114 or 117 and remove it/them, then re-run the reg_model code3.3.4 Check homogeneity of variance
# use the leveneTest() command from the car package to test homogeneity of variance
# uses the 'formula' setup: formula is y~x, where y is our DV and x is our IV
leveneTest(mfq_26~mhealth, data = d)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 2.443 0.08773 .
## 623
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.3.5 Issues with My Data
Our cell sizes are very unbalanced between the mental disorder group levels. A small sample size for one of the levels of our variable (eating disorder, n = 18) limits our power and increases our Type II error rate.
Levene’s test was not significant for our three-level mental health disorder variable with the One-Way ANOVA.
3.4 Run a One-Way ANOVA
aov_model <- aov_ez(data = d,
id = "X",
between = c("mhealth"),
dv = "mfq_26",
anova_table = list(es = "pes"))## Contrasts set to contr.sum for the following variables: mhealth3.5 View One-Way Output
nice(aov_model)## Anova Table (Type 3 tests)
##
## Response: mfq_26
## Effect df MSE F pes p.value
## 1 mhealth 2, 623 0.39 13.23 *** .041 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1ANOVA Effect Size [partial eta-squared] cutoffs from Cohen (1988): * η^2 < 0.01 indicates a trivial effect * η^2 >= 0.01 indicates a small effect * η^2 >= 0.06 indicates a medium effect * η^2 >= 0.14 indicates a large effect
3.6 Visualize One-Way Results
afex_plot(aov_model, x = "mhealth")
3.7 Run One-Way Posthoc Tests
Remember: We ONLY run posthoc IF the ANOVA test is SIGNIFICANT! E.g., only run the posthoc tests on pet type if there is a main effect for pet type
emmeans(aov_model, specs="mhealth")## mhealth emmean SE df lower.CL upper.CL
## anxiety disorder 3.96 0.0697 623 3.82 4.09
## eating disorders 4.02 0.1470 623 3.73 4.31
## none or NA 4.32 0.0271 623 4.27 4.37
##
## Confidence level used: 0.95pairs(emmeans(aov_model, specs="mhealth", adjust="tukey"))## contrast estimate SE df t.ratio p.value
## anxiety disorder - eating disorders -0.0585 0.1630 623 -0.360 0.9310
## anxiety disorder - none or NA -0.3624 0.0747 623 -4.849 <0.0001
## eating disorders - none or NA -0.3039 0.1490 623 -2.035 0.1048
##
## P value adjustment: tukey method for comparing a family of 3 estimates3.8 Write Up One-Way ANOVA Results
To test our hypothesis that there will be a significant difference in people’s openness to trying something new based on people’s level of mental health disorders (eating disorder, an anxiety disorder, no mental health disorder), we used a one-way ANOVA. Our data was unbalanced, with many more people have no mental health disorder participating in our survey (n = 528) than who have an anxiety disorder (n = 80) or an eating disorder (n = 18). This significantly reduces the power of our test and increases the chances of a Type II error. We removed no outliers following visual analysis of Cook’s distance and Residual VS Leverage plots. A non-significant Levene’s test (p = .08) also indicates that our data meets the assumption of homogeneity of variance.
We found a significant effect of pet type, F(2, 623) = 13.23, p <.001 , ηp2 = .041 (small effect; Cohen, 1988). Posthoc tests using Tukey’s HSD adjustment revealed that participants who had an eating disorder (M = 4.02, SE = .15) reported more openness to trying something new than those who had an anxiety disorder (M = 3.96, SE = .07) but less openness to trying something new than those who had no mental health disorder (M = 4.32, SE = .03); participants who own had no mental health disorder reported the highest amount of openness to trying something new overall (see Figure 1 for a comparison).
References
Cohen J. (1988). Statistical Power Analysis for the Behavioral Sciences. New York, NY: Routledge Academic.