1 Loading Libraries

#install.packages("afex")
#install.packages("emmeans")
#install.packages("ggbeeswarm")

library(psych) # for the describe() command
library(ggplot2) # to visualize our results
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
library(expss) # for the cross_cases() command
## Loading required package: maditr
## 
## Use magrittr pipe '%>%' to chain several operations:
##              mtcars %>%
##                  let(mpg_hp = mpg/hp) %>%
##                  take(mean(mpg_hp), by = am)
## 
## 
## Attaching package: 'expss'
## The following object is masked from 'package:ggplot2':
## 
##     vars
library(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':
## 
##     logit
library(afex) # to run the ANOVA 
## 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':
## 
##     lmer
library(ggbeeswarm) # to run plot results
library(emmeans) # for posthoc tests
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'

2 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)

3 State Your Hypothesis

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 (at least 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.

One-Way: We predict that there will be a high significant difference generalized anxiety disorder scores based on participants sexual orientation (asexual, bisexual, gay/lesbian)

Two-Way: We predict that there will be a significant difference in Rosenberg self-esteem scores based on participants’ mental health disorder status (no diagnosis, one diagnosis, multiple diagnoses) and their level of intolerance of uncertainty (low, moderate, high). We also predict that there will be a significant interaction between mental health disorder status and intolerance of uncertainty.

4 Check Your Variables

# you only need to check the variables you're using in the current analysis
# even if 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':    1225 obs. of  8 variables:
##  $ X                 : int  1 321 401 469 520 1390 1422 2183 2247 2482 ...
##  $ sexual_orientation: chr  "Heterosexual/Straight" "Heterosexual/Straight" "Heterosexual/Straight" "Heterosexual/Straight" ...
##  $ mhealth           : chr  "none or NA" "none or NA" "obsessive compulsive disorder" "depression" ...
##  $ big5_open         : num  5.33 4 6 5 3.67 ...
##  $ iou               : num  3.19 2.48 2.81 2.59 2.22 ...
##  $ rse               : num  2.3 3.8 3.1 3 2.6 3 1.3 3 3.2 1.8 ...
##  $ gad               : num  1.86 1 2.14 1.71 1.14 ...
##  $ 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$sexual_orientation <- as.factor(d$sexual_orientation) 
d$mhealth <- as.factor(d$mhealth)
d$row_id <- as.factor(d$row_id)

# we're going to recode our race variable into two groups: poc and white
# in doing so, we are creating a new variable "poc" that has 2 levels
d$sexual_orientation <- as.factor(d$sexual_orientation)
d$mhealth <- as.factor(d$mhealth)
d$row_id <- as.factor(d$row_id)

# Recode to create poc vs asexual
d$poc[d$sexual_orientation == "I use another term"] <- "poc"
d$poc[d$sexual_orientation == "Bi"] <- "poc"
d$poc[d$sexual_orientation == "Gay/Lesbian"] <- "poc"
d$poc[d$sexual_orientation == "Heterosexual/Straight"] <- "poc"
d$poc[d$sexual_orientation == "other"] <- "poc"
d$poc[d$sexual_orientation == "prefer_not"] <- NA
d$poc[d$sexual_orientation == "Asexual"] <- "Asexual"
d$poc <- as.factor(d$poc)

str(d)
## 'data.frame':    1225 obs. of  9 variables:
##  $ X                 : int  1 321 401 469 520 1390 1422 2183 2247 2482 ...
##  $ sexual_orientation: Factor w/ 6 levels "Asexual","Bi",..: 4 4 4 4 6 4 4 4 4 6 ...
##  $ mhealth           : Factor w/ 8 levels "anxiety disorder",..: 5 5 6 3 5 5 5 5 5 5 ...
##  $ big5_open         : num  5.33 4 6 5 3.67 ...
##  $ iou               : num  3.19 2.48 2.81 2.59 2.22 ...
##  $ rse               : num  2.3 3.8 3.1 3 2.6 3 1.3 3 3.2 1.8 ...
##  $ gad               : num  1.86 1 2.14 1.71 1.14 ...
##  $ row_id            : Factor w/ 1225 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ poc               : Factor w/ 2 levels "Asexual","poc": 2 2 2 2 NA 2 2 2 2 NA ...
describe(d$gad)
##    vars    n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 1225 2.05 0.91   1.71    1.96 0.85   1   4     3 0.67    -0.74 0.03
describeBy(d$gad, group = d$sexual_orientation)
## 
##  Descriptive statistics by group 
## group: Asexual
##    vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 35  2.4 1.08   2.14    2.37 1.48   1   4     3 0.13    -1.56 0.18
## ------------------------------------------------------------ 
## group: Bi
##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 155 2.57 0.88   2.57    2.57 1.06   1   4     3 0.03     -1.1 0.07
## ------------------------------------------------------------ 
## group: Gay/Lesbian
##    vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 51  2.6 0.92   2.57     2.6 1.27   1   4     3 0.12    -1.25 0.13
## ------------------------------------------------------------ 
## group: Heterosexual/Straight
##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 866 1.87 0.83   1.57    1.76 0.64   1   4     3 0.94    -0.14 0.03
## ------------------------------------------------------------ 
## group: I use another term
##    vars  n mean   sd median trimmed  mad min max range  skew kurtosis   se
## X1    1 34 2.89 0.97   3.14    2.96 1.27   1   4     3 -0.46    -1.17 0.17
## ------------------------------------------------------------ 
## group: Prefer not to say
##    vars  n mean   sd median trimmed  mad min max range skew kurtosis  se
## X1    1 84 2.14 0.89   1.86    2.07 1.06   1   4     3  0.5    -1.01 0.1
hist(d$gad)

cross_cases(d, sexual_orientation, poc)
 poc 
 Asexual   poc 
 sexual_orientation 
   Asexual  35
   Bi  155
   Gay/Lesbian  51
   Heterosexual/Straight  866
   I use another term  34
   Prefer not to say 
   #Total cases  35 1106

5 Check Your Assumptions

5.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 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 2 error)
  • Homogeneity of variance should be assured (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.

5.1.1 Check levels of IVs

# One-Way
table(d$sexual_orientation)
## 
##               Asexual                    Bi           Gay/Lesbian 
##                    35                   155                    51 
## Heterosexual/Straight    I use another term     Prefer not to say 
##                   866                    34                    84

5.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

# One-Way
leveneTest(gad~ sexual_orientation, data = d)
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value   Pr(>F)   
## group    5  3.5018 0.003791 **
##       1219                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.1.3 Check for outliers using Cook’s distance and Residuals VS Leverage plot

5.1.3.1 Run a Regression to get these 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))

# 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~x1*x2 + c, where y is our DV, x1 is our first IV, x2 is our second IV.
reg_model <- lm(gad~sexual_orientation, data = d) #for One-Way
 #

5.1.3.2 Check for outliers (One-Way)

# Cook's distance
plot(reg_model, 4)

# Residuals VS Leverage
plot(reg_model, 5)

5.1.3.3 Check for outliers (Two-Way)

5.2 Issues with My Data

Our cell sizes are very unbalanced between the GAD type group levels. A small sample size for one of the levels of our variable limits our power and increases our Type 2 error rate.

Levene’s test was significant for our level GAD type variable with the One-Way ANOVA. We are ignoring this and continuing with the analysis anyway for this class.

We identified and removed a single outlier for the One-Way ANOVA.

[UPDATE this section in your HW.]

6 Run an ANOVA

# One-Way
aov_model <- aov_ez(id = "row_id", dv = "gad", between = "sexual_orientation", data = d)
## Contrasts set to contr.sum for the following variables: sexual_orientation

7 View Output

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

8 Visualize Results

# One-Way
afex_plot(aov_model, x = "sexual_orientation")

9 Run Posthoc Tests (One-Way)

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="sexual_orientation", adjust="sidak")
##  sexual_orientation    emmean     SE   df lower.CL upper.CL
##  Asexual                 2.40 0.1450 1218     2.01     2.78
##  Bi                      2.57 0.0689 1218     2.39     2.75
##  Gay/Lesbian             2.60 0.1200 1218     2.28     2.92
##  Heterosexual/Straight   1.87 0.0291 1218     1.79     1.94
##  I use another term      2.89 0.1470 1218     2.50     3.28
##  Prefer not to say       2.14 0.0941 1218     1.90     2.39
## 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 6 estimates
pairs(emmeans(aov_model, specs="sexual_orientation", adjust="sidak"))
##  contrast                                     estimate     SE   df t.ratio
##  Asexual - Bi                                  -0.1774 0.1610 1218  -1.105
##  Asexual - (Gay/Lesbian)                       -0.2035 0.1880 1218  -1.081
##  Asexual - (Heterosexual/Straight)              0.5300 0.1480 1218   3.584
##  Asexual - I use another term                  -0.4948 0.2070 1218  -2.396
##  Asexual - Prefer not to say                    0.2513 0.1730 1218   1.454
##  Bi - (Gay/Lesbian)                            -0.0262 0.1380 1218  -0.189
##  Bi - (Heterosexual/Straight)                   0.7074 0.0748 1218   9.457
##  Bi - I use another term                       -0.3175 0.1620 1218  -1.955
##  Bi - Prefer not to say                         0.4287 0.1170 1218   3.675
##  (Gay/Lesbian) - (Heterosexual/Straight)        0.7336 0.1240 1218   5.936
##  (Gay/Lesbian) - I use another term            -0.2913 0.1900 1218  -1.534
##  (Gay/Lesbian) - Prefer not to say              0.4549 0.1530 1218   2.981
##  (Heterosexual/Straight) - I use another term  -1.0249 0.1500 1218  -6.835
##  (Heterosexual/Straight) - Prefer not to say   -0.2787 0.0985 1218  -2.828
##  I use another term - Prefer not to say         0.7462 0.1750 1218   4.273
##  p.value
##   0.8794
##   0.8890
##   0.0047
##   0.1583
##   0.6936
##   1.0000
##   <.0001
##   0.3693
##   0.0034
##   <.0001
##   0.6422
##   0.0347
##   <.0001
##   0.0538
##   0.0003
## 
## P value adjustment: tukey method for comparing a family of 6 estimates

10 Write Up Results

10.1 One-Way ANOVA

To test our hypothesis that there would be a significant difference in generalized anxiety disorder (GAD) scores based on participants’ sexual orientation, we conducted a one-way ANOVA. Our data was unbalanced, with substantially more participants identifying as heterosexual/straight (n = 866) compared to bisexual (n = 155), gay/lesbian (n = 51), asexual (n = 35), or other identities. This imbalance reduced the statistical power for detecting differences in smaller groups and increased the risk of Type II errors. One outlier was identified and removed based on Cook’s Distance and Residuals vs. Leverage plots. Additionally, Levene’s test was significant (p < .001), indicating a violation of the assumption of homogeneity of variance. This violation increases the risk of Type I error, but we continued with our analysis for educational purposes..

The one-way ANOVA revealed a significant effect of sexual orientation on GAD scores, F(5, 1219) = 31.56, p < .001, partial η² = .115, indicating a large effect size (Cohen, 1988). Post-hoc comparisons with Sidak correction showed that bisexual (M = 2.57, SE = 0.07) and gay/lesbian (M = 2.60, SE = 0.12) participants reported significantly higher GAD scores than heterosexual/straight participants (M = 1.87, SE = 0.03), p < .001. Participants who identified with another term (M = 2.89, SE = 0.15) also reported significantly higher GAD scores than heterosexual/straight individuals, p < .001. Additionally, those who preferred not to say their sexual orientation (M = 2.14, SE = 0.09) scored significantly higher than heterosexuals, p = .029. No significant differences were found between bisexual and gay/lesbian participants (p = 1.000), or between asexual individuals (M = 2.40, SE = 0.15) and other non-heterosexual groups.

References

Cohen J. (1988). Statistical Power Analysis for the Behavioral Sciences. New York, NY: Routledge Academic.