Running a One-Way or Two-Way/Factorial ANOVA
Logan Williams
2024-04-21
1 Loading Libraries
library(psych) # for the describe() command
library(ggplot2) # to visualize our results##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphalibrary(expss) # for the cross_cases() command## Loading required package: maditr##
## To get total summary skip 'by' argument: take_all(mtcars, mean)##
## Attaching package: 'expss'## The following object is masked from 'package:ggplot2':
##
## varslibrary(car) # for the leveneTest() command## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:expss':
##
## recode## The following object is masked from 'package:psych':
##
## logitlibrary(afex) # to run the ANOVA and plot results## Loading required package: lme4## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:expss':
##
## dummy## ************
## Welcome to afex. For support visit: http://afex.singmann.science/## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - Get and set global package options with: afex_options()
## - Set sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************##
## Attaching package: 'afex'## The following object is masked from 'package:lme4':
##
## lmerlibrary(emmeans) # for posthoc tests
library(ggbeeswarm) # one of the error messages said this was needed2 Importing Data
# import the dataset you cleaned previously
# this will be the dataset you'll use throughout the rest of the semester
# use ARC data
d <- read.csv(file="HWData/updated_eammi2_data.csv", header=T)
# for the HW, you may or may not need to use the code below this comment
# check to see if you have a variable called 'X' or 'ResponseId' in your data
names(d)## [1] "income" "usdream" "idea" "swb" "efficacy" "support"# if you do have 'X' or 'ResponseId' (aka your ID variable) then DELETE THE CODE BELOW
# if you don't have those variables, keep the code and run it
d$row_id <- 1:nrow(d)3 State Your Hypothesis
Note: You can chose to run either a one-way ANOVA (a single IV with more than 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.
Two-Way: I predict that there will be significant effects of belief in attainability of American Dream and household income on General Self-efficacy. I predict that those who believe the American Dream is achievable will have higher self-efficacy scores than those who do not believe it is achievable. I also believe that belief in attainability of the American Dream will interact with household income.
4 Check Your Variables
# you only need to check the variables you're using in the current analysis
# although you checked them previously, it's always a good idea to look them over again and be sure that everything is correct
str(d)## 'data.frame': 3141 obs. of 7 variables:
## $ income : chr "1 low" "1 low" "rather not say" "rather not say" ...
## $ usdream : chr "american dream is important and achievable for me" "american dream is important and achievable for me" "american dream is not important and maybe not achievable for me" "american dream is not important and maybe not achievable for me" ...
## $ idea : num 3.75 3.88 3.75 3.75 3.5 ...
## $ swb : num 4.33 4.17 1.83 5.17 3.67 ...
## $ efficacy: num 3.4 3.4 2.2 2.8 3 2.4 2.3 3 3 3.7 ...
## $ support : num 6 6.75 5.17 5.58 6 ...
## $ row_id : int 1 2 3 4 5 6 7 8 9 10 ...# make our categorical variables factors
d$row_id <- as.factor(d$row_id)
#we'll actually use our ID variable for this analysis, so make sure it's coded as a factor
d$usdream <- as.factor(d$usdream)
d$income <- as.factor(d$income)
# check your categorical variables
table(d$usdream) ##
## american dream is important and achievable for me
## 1454
## american dream is important but maybe not achievable for me
## 341
## american dream is not important and maybe not achievable for me
## 585
## american dream is not important but is achievable for me
## 182
## not sure if american dream important
## 579table(d$income) ##
## 1 low 2 middle 3 high rather not say
## 876 879 534 852# also use histograms to examine your continuous variable
hist(d$idea)
hist(d$swb)
hist(d$efficacy)
hist(d$support)
5 Check Your Assumptions
5.1 ANOVA Assumptions
Assumptions checked below:
- DV should be normally distributed across levels of the IV
- All levels of the IVs should have equal number of cases and there should be no empty cells. Cells with low numbers decrease the power of the test (increase change of Type II error)
- Homogeneity of variance should be assured
- Outliers should be identified and removed
If you have confirmed everything else…
- The sampling distribution should be normal. (For a demonstration of what the sampling distribution is, go here.)
5.1.1 Check levels of IVs
# check your categorical variables and make sure they have decent cell sizes
# they should have at least 5 participants in each cell
# but larger numbers are always better
table(d$usdream)##
## american dream is important and achievable for me
## 1454
## american dream is important but maybe not achievable for me
## 341
## american dream is not important and maybe not achievable for me
## 585
## american dream is not important but is achievable for me
## 182
## not sure if american dream important
## 579table(d$income)##
## 1 low 2 middle 3 high rather not say
## 876 879 534 852cross_cases(d, usdream, income)| income | ||||
|---|---|---|---|---|
| 1 low | 2 middle | 3 high | rather not say | |
| usdream | ||||
| american dream is important and achievable for me | 379 | 456 | 289 | 330 |
| american dream is important but maybe not achievable for me | 111 | 81 | 59 | 90 |
| american dream is not important and maybe not achievable for me | 186 | 153 | 63 | 183 |
| american dream is not important but is achievable for me | 56 | 46 | 43 | 37 |
| not sure if american dream important | 144 | 143 | 80 | 212 |
| #Total cases | 876 | 879 | 534 | 852 |
# we're going to recode our race/ethnicity variable into two groups: poc and white
# you may or may not need to combine and/or drop groups for the HW
# watch variable name and name of levels of variable
table(d$usdream)##
## american dream is important and achievable for me
## 1454
## american dream is important but maybe not achievable for me
## 341
## american dream is not important and maybe not achievable for me
## 585
## american dream is not important but is achievable for me
## 182
## not sure if american dream important
## 579d$usdream2[d$usdream == "american dream is important and achievable for me"] <- "Achievable"
d$usdream2[d$usdream == "american dream is important but maybe not achievable for me"] <- "Not Achievable"
d$usdream2[d$usdream == "american dream is not important and maybe not achievable for me"] <- "Not Achievable"
d$usdream2[d$usdream == "american dream is not important but is achievable for me"] <- "Achievable"
table(d$usdream2)##
## Achievable Not Achievable
## 1636 926d$usdream2 <- as.factor(d$usdream2)
cross_cases(d, usdream2, income)| income | ||||
|---|---|---|---|---|
| 1 low | 2 middle | 3 high | rather not say | |
| usdream2 | ||||
| Achievable | 435 | 502 | 332 | 367 |
| Not Achievable | 297 | 234 | 122 | 273 |
| #Total cases | 732 | 736 | 454 | 640 |
# our number of small nb participants is going to hurt us for the two-way anova, but it should be okay for the one-way anova
# so we'll create a new dataframe for the two-way analysis and call it d2
d2 <- subset(d, income != "rather not say")
# to double-check any changes we made
table(d2$usdream2)##
## Achievable Not Achievable
## 1269 653# we'll use the describeBy() command to view skew and kurtosis across our IVs
describeBy(d2$efficacy, group = d2$usdream2)##
## Descriptive statistics by group
## group: Achievable
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1269 3.24 0.41 3.2 3.24 0.44 1.7 4 2.3 -0.1 -0.29 0.01
## ------------------------------------------------------------
## group: Not Achievable
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 653 3.05 0.45 3 3.05 0.44 1.1 4 2.9 -0.25 0.97 0.02describeBy(d2$efficacy, group = d2$income)##
## Descriptive statistics by group
## group: 1 low
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 876 3.15 0.44 3.1 3.16 0.44 1.1 4 2.9 -0.3 0.43 0.01
## ------------------------------------------------------------
## group: 2 middle
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 879 3.15 0.43 3.1 3.15 0.44 1.5 4 2.5 -0.1 -0.02 0.01
## ------------------------------------------------------------
## group: 3 high
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 534 3.18 0.44 3.1 3.18 0.44 1.2 4 2.8 -0.16 0.76 0.02
## ------------------------------------------------------------
## group: rather not say
## NULL5.1.2 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
leveneTest(efficacy~income*usdream2, data = d2)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.4117 0.8409
## 19165.1.3 Check for outliers using Cook’s distance and Residuals vs Leverage plot
5.1.3.1 Run a Regression
# use the lm() command to run the regression
# formula is y~x1*x2, where y is our DV, x1 is our first IV and x2 is our second IV
reg_model2 <- lm(efficacy~income*usdream2, data = d2) #for two-way5.1.3.2 Check for outliers (Two-Way)
# Cook's distance
plot(reg_model2, 4)
# Residuals vs Leverage
plot(reg_model2, 5)
5.2 Issues with My Data
I removed the “rather not say” income group, as this group likely contains a mix of participants in different household income categories who likely selected that option because of some extraneous factor that will be difficult to interpret for the sake of this study. This income group consisted of 640 participants, leaving us with a significantly smaller sample size. Running a Levene’s test yielded a non-significant result (p=0.8409), suggesting there is homogeneity of variance in our sample. Our cell sizes are somewhat unbalanced, though it does not appear this will have a large impact on the results.
6 Run an ANOVA
aov_model2 <- aov_ez(data = d2,
id = "row_id",
between = c("income","usdream2"),
dv = "efficacy",
anova_table = list(es = "pes"))## Warning: Missing values for 367 ID(s), which were removed before analysis:
## 5, 16, 18, 21, 36, 45, 48, 54, 56, 57, ... [showing first 10 only]
## Below the first few rows (in wide format) of the removed cases with missing data.
## row_id income usdream2 .
## # 3 5 2 middle <NA> 3.0
## # 8 16 1 low <NA> 2.9
## # 10 18 1 low <NA> 2.5
## # 12 21 1 low <NA> 3.1
## # 17 36 1 low <NA> 3.7
## # 23 45 2 middle <NA> 2.8## Contrasts set to contr.sum for the following variables: income, usdream27 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
nice(aov_model2)## Anova Table (Type 3 tests)
##
## Response: efficacy
## Effect df MSE F pes p.value
## 1 income 2, 1916 0.18 0.70 <.001 .499
## 2 usdream2 1, 1916 0.18 96.48 *** .048 <.001
## 3 income:usdream2 2, 1916 0.18 4.06 * .004 .017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 18 Visualize Results
afex_plot(aov_model2, x = "income")## NOTE: Results may be misleading due to involvement in interactions
afex_plot(aov_model2, x = "usdream2")## NOTE: Results may be misleading due to involvement in interactions
afex_plot(aov_model2, x = "income", trace = "usdream2")
afex_plot(aov_model2, x = "usdream2", trace = "income")
9 Run Posthoc Tests (Two-Way)
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.
emmeans(aov_model2, specs="usdream2", adjust="tukey")## NOTE: Results may be misleading due to involvement in interactions## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons## usdream2 emmean SE df lower.CL upper.CL
## Achievable 3.25 0.0121 1916 3.22 3.27
## Not Achievable 3.03 0.0179 1916 2.99 3.07
##
## Results are averaged over the levels of: income
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 2 estimatespairs(emmeans(aov_model2, specs="usdream2", adjust="tukey"))## NOTE: Results may be misleading due to involvement in interactions## contrast estimate SE df t.ratio p.value
## Achievable - Not Achievable 0.212 0.0216 1916 9.823 <.0001
##
## Results are averaged over the levels of: incomeemmeans(aov_model2, specs="usdream2", by="income", adjust="sidak")## income = 1 low:
## usdream2 emmean SE df lower.CL upper.CL
## Achievable 3.24 0.0204 1916 3.19 3.28
## Not Achievable 3.07 0.0247 1916 3.02 3.13
##
## income = 2 middle:
## usdream2 emmean SE df lower.CL upper.CL
## Achievable 3.21 0.0190 1916 3.17 3.26
## Not Achievable 3.05 0.0278 1916 2.99 3.11
##
## income = 3 high:
## usdream2 emmean SE df lower.CL upper.CL
## Achievable 3.29 0.0234 1916 3.23 3.34
## Not Achievable 2.98 0.0385 1916 2.89 3.06
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 2 estimatespairs(emmeans(aov_model2, specs="usdream2", by="income", adjust="sidak"))## income = 1 low:
## contrast estimate SE df t.ratio p.value
## Achievable - Not Achievable 0.163 0.0320 1916 5.083 <.0001
##
## income = 2 middle:
## contrast estimate SE df t.ratio p.value
## Achievable - Not Achievable 0.165 0.0337 1916 4.893 <.0001
##
## income = 3 high:
## contrast estimate SE df t.ratio p.value
## Achievable - Not Achievable 0.308 0.0451 1916 6.841 <.0001emmeans(aov_model2, specs="income", by="usdream2", adjust="sidak")## usdream2 = Achievable:
## income emmean SE df lower.CL upper.CL
## 1 low 3.24 0.0204 1916 3.19 3.29
## 2 middle 3.21 0.0190 1916 3.17 3.26
## 3 high 3.29 0.0234 1916 3.23 3.34
##
## usdream2 = Not Achievable:
## income emmean SE df lower.CL upper.CL
## 1 low 3.07 0.0247 1916 3.02 3.13
## 2 middle 3.05 0.0278 1916 2.98 3.12
## 3 high 2.98 0.0385 1916 2.88 3.07
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimatespairs(emmeans(aov_model2, specs="income", by="usdream2", adjust="sidak"))## usdream2 = Achievable:
## contrast estimate SE df t.ratio p.value
## 1 low - 2 middle 0.0237 0.0279 1916 0.850 0.6719
## 1 low - 3 high -0.0480 0.0310 1916 -1.548 0.2689
## 2 middle - 3 high -0.0717 0.0301 1916 -2.382 0.0456
##
## usdream2 = Not Achievable:
## contrast estimate SE df t.ratio p.value
## 1 low - 2 middle 0.0257 0.0372 1916 0.691 0.7689
## 1 low - 3 high 0.0974 0.0458 1916 2.128 0.0846
## 2 middle - 3 high 0.0717 0.0475 1916 1.508 0.2873
##
## P value adjustment: tukey method for comparing a family of 3 estimates10 Write Up Results
10.1 Two-Way ANOVA
To test my hypothesis that income and belief in the attainability of the American Dream would impact self-efficacy and would interact significantly with each other, we used a two-way/factorial ANOVA. Our data met all of the assumptions of the test.
As predicted, there was a significant main effect for belief in attainability of American Dream, F(1,1916) = 96.48, p < .001, ηp2 .048 (small effect size; Cohen, 1988). Participants who believed the American Dream was achievable (M= 3.25, SE= 0.01) reported significantly higher self efficacy scores than participants who did not believe the American Dream was achievable (M= 3.03, SE= 0.01, p<0.001). I believe this a result of self-sufficient ideals surrounding the American Dream. Those who believe that they can achieve the American Dream probably are aware of the effort and resources required to do so and thus need to believe in their own capabilities overcome these obstacles. Those who do not believe they can achieve the American Dream perhaps have a lower self-efficacy score due to outside factors that result in a more negative future outlook in life, caused by something such as mental illness, low self-esteem, or perhaps an external locus of control in combination with an unfavorable environment.
Unlike predicted, I found no main effect for household income. I predicted that household income would affect self-efficacy because I believed the stress and restrictions a lack of financial resources places on lower income people would result in lower self-efficacy scores for lower income participants. As well, I believed the freedom and ability to do more things that comes with having abundant financial resources would result in higher self-efficacy scores. This logic was not reflected in the data, with no significant differences in self-efficacy scores among the three income categories.
Lastly, we found a significant interaction between income and belief in attainability of the american dream (see Figure 2), F(2,1916) = 4.06, p = .017, ηp2 = .004 (trivial effect size; Cohen, 1988). Specifically comparing across income groups among participants who believed in the attainability of the American Dream, participants in the middle household income group (M=3.21, SE 0.02) had significantly lower self-efficacy scores than participants in the higher household income category (M=3.29, SE=0.02, p=0.0456). In contrast, participants in the middle and high household income groups who believed they maybe could not achieve the American Dream did not have significantly different self-efficacy scores.
This was unexpected to me, and I believe is the result of an idea Çakar discusses in the literature. Finding little difference in self-efficacy scores across income groups in his study, he discussed the idea that people with less financial resources gain a sense of self-efficacy from overcoming their struggles (Çakar, 2012). I believe that this could be the reason why there is no difference between the self-efficacy scores of low income and high income participants who believed they could achieve the American Dream, yet middle income participants had lower scores than their high income group peers. Lower income participants likely gain a sense of self-efficacy from overcoming their struggles, and the fact that they believe they can achieve the American Dream despite their lack of resources demonstrates an optimistic outlook of their future that contributes to a greater sense of self-efficacy. In contrast, middle income participants who believe in the American Dream likely do not gain this sense of self-efficacy from overcoming the obstacles that come with lower income and must find it from another source. When looking for other sources, I believe that middle income people who believe in the American Dream are more likely to compare their lives and circumstances to higher income individuals, falling short and leading to a lower self-efficacy score.
References
Cohen J. (1988). Statistical Power Analysis for the Behavioral Sciences. New York, NY: Routledge Academic.
Çakar, F. S. (2012). The Relationship between the Self-Efficacy and Life Satisfaction of Young Adults. International Education Studies, 5(6), 123-130.
SIDE NOTE Where does the “Ç” character fall for alphabetic purposes in comparison to “C”?