library(afex) # to run the ANOVA and plot results
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(emmeans) # for posthoc testsANOVA Homework
Loading Libraries
Importing Data
d <- read.csv(file="Data/mydata.csv", header=T)
# new code! this adds a column with a number for each row. it acts as an identifier for our participations
d$row_id <- 1:nrow(d)State Your Hypothesis - PART OF YOUR WRITEUP
Note: You can chose to run either a one-way ANOVA (a single IV with at least 3 levels) or a two-way/factorial ANOVA (at least two IVs) for the homework. You will need to specify your hypothesis and customize your code based on the choice you make. I will run both versions of the test here for illustrative purposes.
One-way ANOVA: People between the ages of 18 and 25 will have higher need to belong scores than people of all other age groups.
State your hypotheses. Remember, your DV will be a continuous variable. For your IV, you need either one categorical variable with three levels, or two categorical variables with at least two levels each.
Check Your Assumptions
ANOVA Assumptions
- Independence of observations (confirmed by data report)
- All levels of the IVs should have equal number of cases (ideally; in the real world, this varies) and there should be no empty cells. Cells with low numbers increase chance of Type II error. (we will check this below)
- Homogeneity of variance should be assured (we will check this below)
- Outliers should be identified and removed (we will check this below)
- DV should be normally distributed for each level of the IV (we will check this below)
Check levels of IVs
# for a one-way ANOVA
# COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
table(d$age, useNA = "always")
1 between 18 and 25 2 between 26 and 35 3 between 36 and 45 4 over 45
1941 111 37 17
<NA>
0 d2 <- d
# for a two-way ANOVA
# COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
# cross_cases(d2, pet, mhealth)
# if you need to recode
# to drop levels from your variable
# this subsets the data and says that any participant who is coded as 'LEVEL BAD' should be removed
# if you don't need this for the homework, comment it out (add a # at the beginning of the line)
# d <- subset(d, pet != "bird")
# d <- subset(d, pet != "cat and dog")
# d <- subset(d, pet != "fish")
# d <- subset(d, pet != "multiple types of pet")
# d <- subset(d, pet != "other")
#
# # to combine levels
# # this says that where any participant is coded as 'LEVEL BAD' it should be replaced by 'LEVEL GOOD'
# # you can repeat this as needed, changing 'LEVEL BAD' if you have multiple levels that you want to combine into a single level
# # if you don't need this for the homework, comment it out (add a # at the beginning of the line)
d2$age_rc[d2$age == "1 between 18 and 25"] <- "emerging"
d2$age_rc[d2$age == "2 between 26 and 35"] <- "young"
d2$age_rc[d2$age == "3 between 36 and 45"] <- "middle"
d2$age_rc[d2$age == "4 over 45"] <- "middle"
#
# # preview your changes and make sure everything is correct
table(d2$age_rc, useNA = "always")
emerging middle young <NA>
1941 54 111 0 # # or
#
# # check your variable types
str(d)'data.frame': 2106 obs. of 7 variables:
$ gender : chr "f" "m" "m" "f" ...
$ age : chr "1 between 18 and 25" "1 between 18 and 25" "1 between 18 and 25" "1 between 18 and 25" ...
$ moa_role: num 3 2.67 2.5 2 2.67 ...
$ belong : num 2.8 4.2 3.6 4 3.4 4.2 3.9 3.6 2.9 2.5 ...
$ swb : num 4.33 4.17 1.83 5.17 3.67 ...
$ mindful : num 2.4 1.8 2.2 2.2 3.2 ...
$ row_id : int 1 2 3 4 5 6 7 8 9 10 ...str(d2)'data.frame': 2106 obs. of 8 variables:
$ gender : chr "f" "m" "m" "f" ...
$ age : chr "1 between 18 and 25" "1 between 18 and 25" "1 between 18 and 25" "1 between 18 and 25" ...
$ moa_role: num 3 2.67 2.5 2 2.67 ...
$ belong : num 2.8 4.2 3.6 4 3.4 4.2 3.9 3.6 2.9 2.5 ...
$ swb : num 4.33 4.17 1.83 5.17 3.67 ...
$ mindful : num 2.4 1.8 2.2 2.2 3.2 ...
$ row_id : int 1 2 3 4 5 6 7 8 9 10 ...
$ age_rc : chr "emerging" "emerging" "emerging" "emerging" ...#
# # make sure that your IV is recognized as a factor by R
d$age <- as.factor(d$age)
d2$age_rc <- as.factor(d2$age_rc)Check homogeneity of variance
# use the leveneTest() command from the car package to test homogeneity of variance
# uses the 'formula' setup: formula is y~x1*x2, where y is our DV and x1 is our first IV and x2 is our second IV
# for a one-way ANOVA
# COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
leveneTest(belong~age_rc, data = d2)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 2 2.7078 0.06691 .
2103
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# # for a two-way ANOVA
# # COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
# leveneTest(phq~pet_rc*mhealth_rc, data = d2)Check for outliers using Cook’s distance and Residuals vs Leverage plot
Run a Regression
# # use this commented out section only if you need to remove outliers
# # to drop a single outlier, remove the # at the beginning of the line and use this code:
# # d <- subset(d, row_id!=c(1108))
#
# # to drop multiple outliers, remove the # at the beginning of the line and 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~x1*x2 + c, where y is our DV, x1 is our first IV, x2 is our second IV, and c is our covariate
# for a one-way ANOVA
# COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
reg_model <- lm(belong ~ age_rc, data = d2) #for one-way
# # for a two-way ANOVA
# # COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
# reg_model2 <- lm(phq ~ pet_rc*mhealth_rc, data = d2) #for two-wayCheck for outliers
# for a one-way ANOVA
# COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
# Cook's distance (.5 is the cutoff)
plot(reg_model, 4)
# Residuals vs Leverage (we want the red line to be close to the dotted line)
plot(reg_model, 5)
# # for a two-way ANOVA
# # COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
# # Cook's distance
# plot(reg_model2, 4)
#
# # Residuals vs Leverage
# plot(reg_model2, 5)Check Your Variables
# we'll use the describeBy() command to view skew and kurtosis across our IVs and make sure the DV is normally distributed across all of the levels
describeBy(d2$belong, group = d2$age_rc)
Descriptive statistics by group
group: emerging
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 1941 3.24 0.6 3.3 3.26 0.59 1.3 5 3.7 -0.3 -0.01 0.01
------------------------------------------------------------
group: middle
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 54 2.8 0.53 2.8 2.8 0.59 1.7 4 2.3 -0.09 -0.61 0.07
------------------------------------------------------------
group: young
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 111 2.88 0.66 2.8 2.88 0.74 1.4 4.3 2.9 -0.02 -0.82 0.06Issues with My Data - PART OF YOUR WRITEUP
Before running the one-way ANOVA test, I confirmed the independence of observations and checked cell sizes. To ensure my data was appropriate for an ANOVA, I combined the age groups “between 35 and 45” and “over 45” to one category, making my independent variable 3 categories. I confirmed homogeneity of variance using Levine’s test (Pr(>F) = .067) and checked for outliers. This test revealed one outlier, which may skew my data. I will need to be careful in my analyses. I also checked the normality of my dependent variables by all levels of my independent variable and determined that they are normal (all values fall between -2 and +2).
Run an ANOVA
# for a one-way ANOVA
# COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
aov_model <- aov_ez(data = d2,
id = "row_id",
between = c("age_rc"),
dv = "belong",
anova_table = list(es = "pes"))Contrasts set to contr.sum for the following variables: age_rc# # for a two-way ANOVA
# # COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
# aov_model2 <- aov_ez(data = d2,
# id = "row_id",
# between = c("pet_rc","mhealth_rc"),
# dv = "phq",
# anova_table = list(es = "pes"))View Output
Effect size cutoffs from Cohen (1988):
- η2 = 0.01 indicates a small effect
- η2 = 0.06 indicates a medium effect
- η2 = 0.14 indicates a large effect
# for a one-way ANOVA
# COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
nice(aov_model)Anova Table (Type 3 tests)
Response: belong
Effect df MSE F pes p.value
1 age_rc 2, 2103 0.36 32.19 *** .030 <.001
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1# # for a two-way ANOVA
# # COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
# nice(aov_model2)Visualize Results
# for a one-way ANOVA
# COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
afex_plot(aov_model, x = "age_rc")
# # for a two-way ANOVA
# # COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
# afex_plot(aov_model2, x = "pet_rc", trace = "mhealth_rc")
# afex_plot(aov_model2, x = "mhealth_rc", trace = "pet_rc")Run Posthoc Tests
Only run posthocs if the test is significant! E.g., only run the posthoc tests on gender if there is a main effect for gender.
# for a one-way ANOVA
# COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
emmeans(aov_model, specs="age_rc", adjust="tukey")Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons age_rc emmean SE df lower.CL upper.CL
emerging 3.24 0.0136 2103 3.21 3.28
middle 2.80 0.0817 2103 2.61 3.00
young 2.88 0.0570 2103 2.74 3.02
Confidence level used: 0.95
Conf-level adjustment: sidak method for 3 estimates pairs(emmeans(aov_model, specs="age_rc", adjust="tukey")) contrast estimate SE df t.ratio p.value
emerging - middle 0.4395 0.0828 2103 5.309 <.0001
emerging - young 0.3639 0.0586 2103 6.215 <.0001
middle - young -0.0756 0.0996 2103 -0.759 0.7281
P value adjustment: tukey method for comparing a family of 3 estimates # # for a two-way ANOVA
# # COMMENT THIS SECTION OUR FROM THE HW IF YOU DO NOT USE IT
# emmeans(aov_model2, specs="pet_rc", adjust="tukey")
# pairs(emmeans(aov_model2, specs="pet_rc", adjust="tukey"))
#
# emmeans(aov_model2, specs="mhealth_rc", adjust="tukey")
# pairs(emmeans(aov_model2, specs="mhealth_rc", adjust="tukey"))
#
# emmeans(aov_model2, specs="pet_rc", by="mhealth_rc", adjust="sidak")
# pairs(emmeans(aov_model2, specs="pet_rc", by="mhealth_rc", adjust="sidak"))
#
# emmeans(aov_model2, specs="mhealth_rc", by="pet_rc", adjust="sidak")
# pairs(emmeans(aov_model2, specs="mhealth_rc", by="pet_rc", adjust="sidak"))Write Up Results
One-Way ANOVA
Write up your results. Again, make sure to maintain an appropriate tone, and follow APA guidelines for reporting statistical results. I recommend following the below outline:
- Briefly restate your hypothesis
- Describe any issues with your data (you can copy/paste from above, just make sure everything flows).
- Report your results. Make sure to include your F-value, degrees of freedom, p-value, and effect size. Since we are showing our means and standard deviations for the levels of our IV in our plot, you do NOT have to report them in the text (normally, you would report it like this: (M = #.##, SD = .##), repeating for each level).
- Specify where the differences occur. For instance, one level of the IV may be significantly different than the other two, or there may be multiple differences.
- Describe your interaction, if you ran a two-way ANOVA and there is one.
- Interpret your effect size (trivial, small, medium, or large) and include the citation.
- Make sure to include a reference to Figure 1 (created using the code below)
My hypothesis was that emerging adults would have higher need for belonging scores than other age groups. My data was almost perfect for an ANOVA, but there was one outlier that may have affected results. I found a significant main effect of age on need for belonging, F(2,2103) = 32.19, p<.001, ηp2 = .03. Emerging adults were statistically different from young adults (p<.0001) and middle adults (p<.0001). There was no statistically significant difference between young and middle aged adults (p=.7281). My effect size was .3, which is a medium/large effect size (Cohen 1988).
- We did not find a significant main effect for pet ownership (p = .400)
- We did find a significant main effect of mental health diagnosis on depression score, F(1, 1197) = 120.95, p<.001, ηp2 = .09
- We did not find a significant interaction for pet ownership and mental health diagnosis (p = .800)
References
Cohen J. (1988). Statistical Power Analysis for the Behavioral Sciences. New York, NY: Routledge Academic.