#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 tests
# 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.
We predict that there will be a significant difference in people’s level of Satisfaction with Life based on their political party (Apolitical, Democrat, Republican, Independent).
# you only need to check the variables you're using in the current analysis
str(d)
## 'data.frame': 3146 obs. of 8 variables:
## $ ResponseID: chr "R_BJN3bQqi1zUMid3" "R_2TGbiBXmAtxywsD" "R_12G7bIqN2wB2N65" "R_39pldNoon8CePfP" ...
## $ gender : chr "f" "m" "m" "f" ...
## $ party_rc : chr "democrat" "independent" "apolitical" "apolitical" ...
## $ swb : num 4.33 4.17 1.83 5.17 3.67 ...
## $ belong : num 2.8 4.2 3.6 4 3.4 4.2 3.9 3.6 2.9 2.5 ...
## $ socmeduse : int 47 23 34 35 37 13 37 43 37 29 ...
## $ stress : num 3.3 3.3 4 3.2 3.1 3.5 3.3 2.4 2.9 2.7 ...
## $ 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$party_rc <- as.factor(d$party_rc)
d$row_id <- as.factor(d$row_id)
# check that our categorical variables of interest are now factors
str(d)
## 'data.frame': 3146 obs. of 8 variables:
## $ ResponseID: chr "R_BJN3bQqi1zUMid3" "R_2TGbiBXmAtxywsD" "R_12G7bIqN2wB2N65" "R_39pldNoon8CePfP" ...
## $ gender : chr "f" "m" "m" "f" ...
## $ party_rc : Factor w/ 4 levels "apolitical","democrat",..: 2 3 1 1 1 1 3 2 2 1 ...
## $ swb : num 4.33 4.17 1.83 5.17 3.67 ...
## $ belong : num 2.8 4.2 3.6 4 3.4 4.2 3.9 3.6 2.9 2.5 ...
## $ socmeduse : int 47 23 34 35 37 13 37 43 37 29 ...
## $ stress : num 3.3 3.3 4 3.2 3.1 3.5 3.3 2.4 2.9 2.7 ...
## $ row_id : Factor w/ 3146 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
# check our DV skew and kurtosis
describe(d$swb)
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3146 4.48 1.32 4.67 4.53 1.48 1 7 6 -0.36 -0.46 0.02
# we'll use the describeBy() command to view our DV's skew and kurtosis across our IVs' levels
describeBy(d$swb, group = d$party_rc)
##
## Descriptive statistics by group
## group: apolitical
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 439 4.29 1.29 4.33 4.32 1.48 1 7 6 -0.16 -0.64 0.06
## ------------------------------------------------------------
## group: democrat
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1597 4.37 1.33 4.5 4.42 1.48 1 7 6 -0.32 -0.51 0.03
## ------------------------------------------------------------
## group: independent
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 327 4.39 1.38 4.33 4.46 1.48 1 7 6 -0.38 -0.4 0.08
## ------------------------------------------------------------
## group: republican
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 783 4.84 1.23 5 4.92 1.24 1 7 6 -0.53 -0.15 0.04
# also use histograms to examine your continuous variable
hist(d$swb)
table(d$party_rc)
##
## apolitical democrat independent republican
## 439 1597 327 783
# REMEMBparty_rc# REMEMBER your test's level of POWER is determined by your SMALLEST subsample
# 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 problematic 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))
# use the lm() command to run the regression
# formula is y~x, where y is our DV, x is our IV
reg_model <- lm(swb~party_rc, data = d )
# 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 code
# 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(swb~party_rc, data = d)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 3.8918 0.008662 **
## 3142
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Our cell sizes are highly unbalanced between the levels of our political party variable. While the democrat group has a large sample size (n = 1597), the independent group is notably smaller (n = 327), which can limit our statistical power and increase Type II error rates when evaluating differences between specific groups.
Levene’s Test was statistically significant (F(3, 3142) = 3.8918, p = 0.008662), indicating that the assumption of equal variances across our four groups has been violated. However, for the purposes of this analysis and the class, we are acknowledging this limitation and proceeding with the One-Way ANOVA.
aov_model <- aov_ez(data = d,
id = "ResponseID",
between = c("party_rc"),
dv = "swb" ,
anova_table = list(es = "pes"))
## Contrasts set to contr.sum for the following variables: party_rc
nice(aov_model)
## Anova Table (Type 3 tests)
##
## Response: swb
## Effect df MSE F pes p.value
## 1 party_rc 3, 3142 1.70 27.30 *** .025 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
ANOVA 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
afex_plot(aov_model, x = "party_rc")
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="party_rc")
## party_rc emmean SE df lower.CL upper.CL
## apolitical 4.29 0.0623 3142 4.17 4.41
## democrat 4.37 0.0327 3142 4.30 4.43
## independent 4.39 0.0722 3142 4.25 4.53
## republican 4.84 0.0466 3142 4.75 4.93
##
## Confidence level used: 0.95
pairs(emmeans(aov_model, specs="party_rc", adjust="tukey"))
## contrast estimate SE df t.ratio p.value
## apolitical - democrat -0.0790 0.0703 3142 -1.123 0.6755
## apolitical - independent -0.1009 0.0953 3142 -1.058 0.7151
## apolitical - republican -0.5499 0.0778 3142 -7.066 <0.0001
## democrat - independent -0.0219 0.0792 3142 -0.276 0.9926
## democrat - republican -0.4709 0.0569 3142 -8.270 <0.0001
## independent - republican -0.4490 0.0859 3142 -5.225 <0.0001
##
## P value adjustment: tukey method for comparing a family of 4 estimates
To test our hypothesis that there will be a significant difference in people’s level of subjective well-being based on political affiliation (apolitical, democrat, independent, republican), we used a one-way ANOVA. Our data was somewhat unbalanced, with many more people identifying as a democrat (n = 1357) or republican (n = 1007) participating in our survey than people who identify as independent (n = 422) or apolitical (n = 360). This significantly reduces the power of our test and increases the chances of a Type II error. We continued with our analysis for the purpose of this class. A significant Levene’s test (p=.008) indicated that our data failed the assumption of homogeneity of variance. This suggests that there is an increased chance of Type I error. We continued with our analysis for the purpose of this class.
We found a significant effect of political affiliation, F(3, 3142) = 27.30, p<.001, np2= .042 (small effect; Cohen, 1988). Post-hoc tests using Tukey’s HSD adjustment revealed that participants who identify as republican (M = 4.84, SE = .05) reported significantly higher subjective well-being than those who identify as democrat (M = 4.37, SE = .03), independent (M = 4.39, SE = .07), or apolitical (M = 4.29, SE = .06).Posthoc tests using Tukey’s HSD adjustment revealed that participants who identify as a democrat (M = 4.37, SE = .03) reported less subjective well-being than those who identify as a republican (M = 4.84, SE = .05) but similar subjective well-being to those who identify as an independent (M = 4.39, SE = .07) or apolitical (M = 4.29, SE = .06); participants who identify as a republican reported the highest amount of subjective well-being overall (see Figure 1 for comparison).
References
Cohen J. (1988). Statistical Power Analysis for the Behavioral Sciences. New York, NY: Routledge Academic.