Study 1
Data preparation
# Creating a smaller data frame with only variables needed for current analyses
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
data1<- data.1 %>%
select(ID, event2_end, event5_end, lec5_08a, lec5_09a, Gender, apref_e2, apref_e5,
Sample, Age, er15e5, er16e5, er33e5, er34e5,
er35e5, er15e2, er16e2, er33e2, er34e2, er35e2, apref_e1, apref_e3, apref_e4, apref_e6)
################# GENDER #########################
table(data1$Gender) # Gender frequences
##
## Male Female
## 175 573
## Genderqueer/non-conforming Transgender
## 5 3
## Other (please specify)
## 4
data1<- data1 %>% #Created new binary gender variable for men vs. women/TNB+
mutate(gender.bin = case_when(
data1$Gender=="Male" ~ -1, #Man
data1$Gender=="Female" ~ 1, #WTNB+
data1$Gender=="Genderqueer/non-conforming" ~ 1,#WTNB+
data1$Gender=="Transgender" ~ 1, #WTNB+
data1$Gender=="Other (please specify)" ~ 1, #WTNB+
))
table(data1$gender.bin) #checking new frequencies, looks correct, 77% of sample is WTNB+
##
## -1 1
## 175 585
prop.table(table(data1$gender.bin))
##
## -1 1
## 0.2302632 0.7697368
data1$gender.bin<- as.factor(data1$gender.bin) # Converting gender to factor
################# SAMPLE #########################
data1<- data1 %>%
mutate(sample.fac = case_when(
data1$Sample=="UMN" ~ 1,
data1$Sample=="WWU" ~ 2,
data1$Sample=="qualtrics" ~ 3))
data1$sample.fac<- as.factor(data1$sample.fac) # Converting sample to factor
################# ASA ENDING #########################
data1<- data1 %>%
mutate(ASA.ending.fac = case_when(
data1$event2_end=="Negative" ~ 1,
data1$event2_end=="Redemption" ~ 2,
data1$event2_end=="Survivor" ~ 3))
data1$ASA.ending.fac<- as.factor(data1$ASA.ending.fac) # Converting ASA ending to factor
table(data1$ASA.ending.fac)
##
## 1 2 3
## 250 252 258
################# CSA ENDING #########################
data1<- data1 %>%
mutate(CSA.ending.fac = case_when(
data1$event5_end=="Negative" ~ 1,
data1$event5_end=="Redemption" ~ 2,
data1$event5_end=="Survivor" ~ 3))
data1$CSA.ending.fac<- as.factor(data1$CSA.ending.fac) # Converting CSA ending to factor
################# Participant SV History #########################
data1$lec5_08a<- as.character(data1$lec5_08a) # First convert to character, then switch back to factor to avoid R error
data1$lec5_09a<- as.character(data1$lec5_09a)
# Replace NAs with 0s (Nas in this case indicate that experience does not apply)
data1$lec5_08a[is.na(data1$lec5_08a)] <- 0
data1$lec5_09a[is.na(data1$lec5_09a)] <- 0
data1<- data1 %>%
mutate(SV = case_when(
data1$lec5_08a=="Happened_to_me" | data1$lec5_09a=="Happened_to_me" ~ 1,
data1$lec5_08a==0 & data1$lec5_09a==0 ~ -1
))
table(data1$SV)
##
## -1 1
## 435 325
data1$SV<- as.factor(data1$SV)
svcheck1<- data1 %>% #creating dataframe with SV variables to check that they're coded correctly
select(lec5_08a, lec5_09a, SV)
Likability Ratings
which(colnames(data1)=="er15e5")
## [1] 11
which(colnames(data1)=="er35e2")
## [1] 20
############ Pulling likability items from full dataset just to make sure nothing was miscoded originally. The means for the new scale and the original are the same and everything seems to be working the same, so looks good.
# Recoding likability items to numeric where 1 = strongly disagree and 5 = strongly agree
data1 <- data1 %>%
mutate(
er15e5 = dplyr::recode(er15e5, "Strongly disagree"=1,
"Disagree"=2,
"Neither agree nor disagree"=3,
"Agree"=4,
"Strongly agree"=5))
data1 <- data1 %>%
mutate(
er16e5 = dplyr::recode(er16e5, "Strongly disagree"=1,
"Disagree"=2,
"Neither agree nor disagree"=3,
"Agree"=4,
"Strongly agree"=5))
data1 <- data1 %>%
mutate(
er15e2 = dplyr::recode(er15e2, "Strongly disagree"=1,
"Disagree"=2,
"Neither agree nor disagree"=3,
"Agree"=4,
"Strongly agree"=5))
data1 <- data1 %>%
mutate(
er33e5 = dplyr::recode(er33e5, "Strongly disagree"=1,
"Disagree"=2,
"Neither agree nor disagree"=3,
"Agree"=4,
"Strongly agree"=5))
data1 <- data1 %>%
mutate(
er34e5 = dplyr::recode(er34e5, "Strongly disagree"=1,
"Disagree"=2,
"Neither agree nor disagree"=3,
"Agree"=4,
"Strongly agree"=5))
data1 <- data1 %>%
mutate(
er35e5 = dplyr::recode(er35e5, "Strongly disagree"=1,
"Disagree"=2,
"Neither agree nor disagree"=3,
"Agree"=4,
"Strongly agree"=5))
data1 <- data1 %>%
mutate(
er16e2 = dplyr::recode(er16e2, "Strongly disagree"=1,
"Disagree"=2,
"Neither agree nor disagree"=3,
"Agree"=4,
"Strongly agree"=5))
data1 <- data1 %>%
mutate(
er33e2 = dplyr::recode(er33e2, "Strongly disagree"=1,
"Disagree"=2,
"Neither agree nor disagree"=3,
"Agree"=4,
"Strongly agree"=5))
data1 <- data1 %>%
mutate(
er34e2 = dplyr::recode(er34e2, "Strongly disagree"=1,
"Disagree"=2,
"Neither agree nor disagree"=3,
"Agree"=4,
"Strongly agree"=5))
data1 <- data1 %>%
mutate(
er35e2 = dplyr::recode(er35e2, "Strongly disagree"=1,
"Disagree"=2,
"Neither agree nor disagree"=3,
"Agree"=4,
"Strongly agree"=5))
data1 <- data1 %>%
rowwise() %>%
mutate(Likability.n1 = mean(c(er15e5, er16e5, er33e5, er34e5,
er35e5, er15e2, er16e2, er33e2, er34e2, er35e2), na.rm = TRUE),
CSA.likability = mean(c(er15e5, er16e5, er33e5, er34e5, er35e5), na.rm = TRUE),
ASA.likability = mean(c(er15e2, er16e2, er33e2, er34e2, er35e2), na.rm = TRUE))%>%
ungroup()
mean(data1$Likability.n1, na.rm=T)
## [1] 3.171053
mean(data1$CSA.likability, na.rm=T)
## [1] 3.193412
mean(data1$ASA.likability, na.rm=T)
## [1] 3.150526
################# LIKABILITY RATINGS #########################
# Averaging ASA and CSA likability ratings, as previous research using this dataset found no mean differences
data1 <- data1 %>%
rowwise() %>%
mutate(Likability = mean(c(apref_e2, apref_e5), na.rm = TRUE))%>%
ungroup()
likability.check<-data1 %>% #creating dataframe likability ratings to check that they look correct
select(apref_e2, apref_e5, Likability)
mean(data1$Likability.n1, na.rm=T)
## [1] 3.171053
mean(data1$CSA.likability, na.rm=T)
## [1] 3.193412
mean(data1$ASA.likability, na.rm=T)
## [1] 3.150526
mean(data1$Likability, na.rm=T)
## [1] 3.171053
# Checking that these variables are the type I want, all factor except likability which is numeric so good
class(data1$gender.bin)
## [1] "factor"
class(data1$sample.fac)
## [1] "factor"
class(data1$ASA.ending.fac)
## [1] "factor"
class(data1$CSA.ending.fac)
## [1] "factor"
class(data1$SV)
## [1] "factor"
class(data1$Likability)
## [1] "numeric"
# Study 1 Likability Reliability
data1 %>%
select(er15e5:er35e2) %>%
psych::alpha( ,check.keys = F, na.rm = TRUE)
##
## Reliability analysis
## Call: psych::alpha(x = ., na.rm = TRUE, check.keys = F)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.84 0.84 0.9 0.35 5.4 0.0088 3.2 0.59 0.29
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.83 0.84 0.86
## Duhachek 0.83 0.84 0.86
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## er15e5 0.82 0.82 0.87 0.34 4.6 0.0100 0.032 0.29
## er16e5 0.82 0.82 0.87 0.34 4.7 0.0099 0.033 0.29
## er33e5 0.84 0.84 0.90 0.36 5.1 0.0094 0.048 0.29
## er34e5 0.83 0.83 0.88 0.35 4.9 0.0095 0.043 0.29
## er35e5 0.84 0.84 0.89 0.36 5.2 0.0092 0.042 0.31
## er15e2 0.82 0.82 0.87 0.34 4.7 0.0098 0.032 0.29
## er16e2 0.82 0.82 0.87 0.34 4.7 0.0099 0.033 0.29
## er33e2 0.83 0.83 0.89 0.35 4.9 0.0097 0.048 0.25
## er34e2 0.83 0.83 0.88 0.35 5.0 0.0095 0.042 0.29
## er35e2 0.84 0.84 0.89 0.36 5.1 0.0092 0.041 0.29
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## er15e5 759 0.71 0.72 0.72 0.62 3.3 0.89
## er16e5 759 0.70 0.71 0.71 0.61 3.2 0.85
## er33e5 759 0.58 0.58 0.50 0.47 3.7 0.87
## er34e5 759 0.64 0.64 0.60 0.53 3.0 0.95
## er35e5 759 0.57 0.57 0.52 0.46 2.8 0.84
## er15e2 760 0.69 0.69 0.69 0.59 3.3 0.94
## er16e2 760 0.70 0.70 0.70 0.60 3.2 0.90
## er33e2 760 0.65 0.64 0.58 0.54 3.6 0.97
## er34e2 760 0.63 0.63 0.58 0.52 3.0 0.93
## er35e2 760 0.58 0.58 0.52 0.46 2.6 0.94
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## er15e5 0.05 0.09 0.43 0.38 0.05 0
## er16e5 0.05 0.08 0.51 0.31 0.04 0
## er33e5 0.03 0.03 0.30 0.48 0.17 0
## er34e5 0.06 0.24 0.43 0.21 0.05 0
## er35e5 0.08 0.25 0.51 0.15 0.01 0
## er15e2 0.05 0.11 0.37 0.40 0.07 0
## er16e2 0.05 0.09 0.48 0.31 0.07 0
## er33e2 0.03 0.08 0.31 0.40 0.17 0
## er34e2 0.06 0.22 0.43 0.25 0.05 0
## er35e2 0.15 0.29 0.42 0.13 0.01 0
mean(data1$Likability, na.rm=T)
## [1] 3.171053
sd(data1$Likability, na.rm=T)
## [1] 0.5873537
# Average Non-SV likability
data1 <- data1 %>%
rowwise() %>%
mutate(Likability.nonSV = mean(c(apref_e1, apref_e3, apref_e4, apref_e6), na.rm = TRUE))%>%
ungroup()
# Difference scores: SV and non-SV likability
data1$DiffScores <- (data1$Likability.nonSV - data1$Likability)
# Positive scores mean that participants rated narrators in the non-SV vignettes as more likable. Negative scores mean that participants rated victims of sexual violence as more likable.
Data visualization
Likability.mean <- mean(data1$Likability)
hist(data1$Likability, main="Histogram of Likability", xlab="Likability", col='hotpink', sub=paste("Skewness:",
round(e1071::skewness(data1$Likability, na.rm=TRUE), 2)))
abline(v = Likability.mean, col = 'orange') # Overlay mean on histogram
qqnorm(data1$Likability, pch = 1, frame = FALSE, main="QQ Plot of Likability")
qqline(data1$Likability, col = "hotpink", lwd = 2)
# Means and sd
mean(data1$Likability, na.rm=T)
## [1] 3.171053
sd(data1$Likability, na.rm=T)
## [1] 0.5873537
library(ggplot2)
ggplot(data1, aes(SV, Likability)) +
geom_boxplot(fill = "lightpink", color = "hotpink") +
labs(x="Assault History (No SV = -1; SV = 1)", y="Likability") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggplot(data1, aes(gender.bin, Likability)) +
geom_boxplot(fill = "lightpink", color = "hotpink") +
labs(x="Gender (Male = -1; WTNB+ = 1)", y="Likability") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggplot(data1, aes(sample.fac, Likability)) +
geom_boxplot(fill = "lightpink", color = "hotpink") +
labs(x="Sample (UMN = 1; WWU = 2; Qualtrics = 3)", y="Likability") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggplot(data1, aes(ASA.ending.fac, Likability)) + # For this study they had same endings for CSA and ASA conditions
geom_boxplot(fill = "lightpink", color = "hotpink") +
labs(x="Story Ending (Negative = 1; Redemptive = 2; Survivor = 3)", y="Likability") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
################### AGE ##########################
data1$Age<- as.numeric(data1$Age)
# Someone entered 188 as their age which is throwing SD for WWU off (should be lower than 10)
data1$Age<-replace(data1$Age, data1$Age==188, NA) #changing the 188 value to NA
hist(data1$Age) #instead of removing the whole person just remove this value
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
scatterplot(Likability ~ Age, data=data1,
xlab="Age", ylab="Likability",
main="Age by Likability", col="hotpink")
#################### NON SV LIKABILITY ################
ggplot(data1, aes(gender.bin, Likability.nonSV)) +
geom_boxplot(fill = "lightpink", color = "hotpink") +
labs(x="Gender (Male = -1; WTNB+ = 1)", y="Non SV Likability") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
#################### DIFFERENCE SCORES ################
Diffscore.mean <- mean(data1$DiffScores)
hist(data1$DiffScores, main="Histogram of Difference Scores", xlab="Difference Scores (Likability.nonSV - Likability", col='hotpink', sub=paste("Skewness:",
round(e1071::skewness(data1$DiffScores, na.rm=TRUE), 2)))
abline(v = Diffscore.mean, col = 'orange') # Overlay mean on histogram
qqnorm(data1$DiffScores, pch = 1, frame = FALSE, main="QQ Plot of Difference SCores")
qqline(data1$DiffScores, col = "hotpink", lwd = 2)
Descriptives
################# Sexual Violence (SV) History ####################
library(expss)
## Loading required package: maditr
##
## To select rows from data: rows(mtcars, am==0)
##
## Attaching package: 'maditr'
## The following objects are masked from 'package:dplyr':
##
## between, coalesce, first, last
## The following object is masked from 'package:base':
##
## sort_by
##
## Attaching package: 'expss'
## The following object is masked from 'package:car':
##
## recode
## The following object is masked from 'package:ggplot2':
##
## vars
## The following objects are masked from 'package:dplyr':
##
## compute, contains, na_if, recode, vars, where
cross_cases(data1, lec5_08a, lec5_09a) # lec8 is sexual assault, lec9 is unwanted sexual experiences
|
|
 lec5_09aÂ
|
|
|
 0Â
|
 Happened_to_meÂ
|
|
 lec5_08aÂ
|
|
   0Â
|
435
|
169
|
|
   Happened_to_meÂ
|
20
|
136
|
|
   #Total casesÂ
|
455
|
305
|
################# SV History by Likability #############
group_by(data1, SV) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 2 × 6
## SV count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -1 435 3.13 0.620 3.1 0.8
## 2 1 325 3.23 0.536 3.2 0.7
################# Gender by Likability #############
group_by(data1, gender.bin) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 2 × 6
## gender.bin count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -1 175 3.11 0.641 3 0.900
## 2 1 585 3.19 0.570 3.2 0.8
################# Sample by Likability #############
group_by(data1, Sample) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 3 × 6
## Sample count mean sd median IQR
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 UMN 290 3.21 0.553 3.25 0.7
## 2 WWU 280 3.19 0.513 3.2 0.8
## 3 qualtrics 190 3.07 0.718 3 0.975
################# Ending by Likability #############
group_by(data1, event2_end) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 3 × 6
## event2_end count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Negative 250 2.82 0.519 2.8 0.500
## 2 Redemption 252 3.32 0.557 3.4 0.7
## 3 Survivor 258 3.37 0.521 3.4 0.7
mean(data1$Likability, na.rm=T)
## [1] 3.171053
sd(data1$Likability, na.rm=T)
## [1] 0.5873537
median(data1$Likability)
## [1] 3.2
################# Sample by Gender #############
cross_cases(data1, Sample, gender.bin)
|
|
 gender.binÂ
|
|
|
 -1Â
|
 1Â
|
|
 SampleÂ
|
|
   UMNÂ
|
90
|
200
|
|
   WWUÂ
|
46
|
234
|
|
   qualtricsÂ
|
39
|
151
|
|
   #Total casesÂ
|
175
|
585
|
################# SV by Gender #############
cross_cases(data1, SV, gender.bin)
|
|
 gender.binÂ
|
|
|
 -1Â
|
 1Â
|
|
 SVÂ
|
|
   -1Â
|
147
|
288
|
|
   1Â
|
28
|
297
|
|
   #Total casesÂ
|
175
|
585
|
cross_cases(data1, lec5_08a, gender.bin) # sexual assault
|
|
 gender.binÂ
|
|
|
 -1Â
|
 1Â
|
|
 lec5_08aÂ
|
|
   0Â
|
162
|
442
|
|
   Happened_to_meÂ
|
13
|
143
|
|
   #Total casesÂ
|
175
|
585
|
cross_cases(data1, lec5_09a, gender.bin) # other unwanted or uncomfortable sexual experiences
|
|
 gender.binÂ
|
|
|
 -1Â
|
 1Â
|
|
 lec5_09aÂ
|
|
   0Â
|
152
|
303
|
|
   Happened_to_meÂ
|
23
|
282
|
|
   #Total casesÂ
|
175
|
585
|
data1<- data1 %>% #Created new interaction term for descriptives that includes SV history x gender
mutate(SVxGen = case_when(
data1$gender.bin==-1 & data1$SV==-1 ~ "M_no", # Men with no SV history
data1$gender.bin==-1 & data1$SV==1 ~ "M_yes", # Men with SV history
data1$gender.bin==1 & data1$SV==-1 ~ "WTNB_no", #WTNB+ with no SV history
data1$gender.bin==1 & data1$SV==1 ~ "WTNB_yes", #WTNB+ with SV history
))
table(data1$SVxGen) # Frequencies by gender and SV history
##
## M_no M_yes WTNB_no WTNB_yes
## 147 28 288 297
group_by(data1, SVxGen) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 4 × 6
## SVxGen count mean sd median IQR
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 M_no 147 3.11 0.669 3 0.900
## 2 M_yes 28 3.12 0.477 3.05 0.725
## 3 WTNB_no 288 3.13 0.594 3.15 0.8
## 4 WTNB_yes 297 3.24 0.541 3.2 0.7
################## AGE ##########################
group_by(data1, Sample) %>%
summarise(
count = n(),
mean = mean(Age, na.rm = TRUE),
sd = sd(Age, na.rm = TRUE),
median = median(Age, na.rm = TRUE),
range = range(Age, na.rm=TRUE),
IQR = IQR(Age, na.rm = TRUE)
)
## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `summarise()` has grouped output by 'Sample'. You can override using the
## `.groups` argument.
## # A tibble: 6 × 7
## # Groups: Sample [3]
## Sample count mean sd median range IQR
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 UMN 290 20.1 2.52 20 18 2
## 2 UMN 290 20.1 2.52 20 42 2
## 3 WWU 280 20.1 4.42 19 17 2
## 4 WWU 280 20.1 4.42 19 67 2
## 5 qualtrics 190 50.9 15.8 53 18 25.8
## 6 qualtrics 190 50.9 15.8 53 82 25.8
S1.samplexage<- aov(Age ~ Sample, data=data1)
summary(S1.samplexage)
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 2 134905 67453 939 <2e-16 ***
## Residuals 756 54306 72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingness
TukeyHSD(S1.samplexage)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Age ~ Sample, data = data1)
##
## $Sample
## diff lwr upr p adj
## UMN-qualtrics -30.74464610 -32.602307 -28.886985 0.0000000
## WWU-qualtrics -30.80703641 -32.679137 -28.934936 0.0000000
## WWU-UMN -0.06239031 -1.731469 1.606688 0.9957608
library(rstatix)
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
eta_squared(S1.samplexage)
## Sample
## 0.7129885
mean(data1$Age, na.rm=T)
## [1] 27.81818
sd(data1$Age, na.rm=T)
## [1] 15.79932
range(data1$Age, na.rm=T)
## [1] 17 82
################# SV History by Difference Scores #############
group_by(data1, SV) %>%
summarise(
count = n(),
mean = mean(DiffScores, na.rm = TRUE),
sd = sd(DiffScores, na.rm = TRUE),
median = median(DiffScores, na.rm = TRUE),
IQR = IQR(DiffScores, na.rm = TRUE)
)
## # A tibble: 2 × 6
## SV count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -1 435 0.221 0.400 0.200 0.450
## 2 1 325 0.0949 0.397 0.0500 0.4
################# Gender by Difference Scores #############
group_by(data1, gender.bin) %>%
summarise(
count = n(),
mean = mean(DiffScores, na.rm = TRUE),
sd = sd(DiffScores, na.rm = TRUE),
median = median(DiffScores, na.rm = TRUE),
IQR = IQR(DiffScores, na.rm = TRUE)
)
## # A tibble: 2 × 6
## gender.bin count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -1 175 0.24 0.396 0.200 0.45
## 2 1 585 0.145 0.403 0.100 0.4
mean(data1$Likability.nonSV)
## [1] 3.338032
sd(data1$Likability.nonSV)
## [1] 0.5314687
Outliers
Likabilitymin <- mean(data1$Likability, na.rm=T) - (3*(sd(data1$Likability, na.rm=T)))
Likabilitymax <- mean(data1$Likability, na.rm=T) + (3*(sd(data1$Likability, na.rm=T)))
data1$Likability[which(data1$Likability < Likabilitymin | data1$Likability > Likabilitymax)] # 6 outliers
## [1] 1.3 1.4 5.0 1.2 5.0 1.0 1.4
library(datawizard)
data1$Likability.win<- winsorize(data1$Likability, method="zscore",
threshold=3, robust=TRUE) # Robust= TRUE means that values are winsorized based on their Median Absolute Deviation (MAD)
mean(data1$Likability.win)
## [1] 3.172056
mean(data1$Likability)
## [1] 3.171053
Diffscoremin <- mean(data1$DiffScores, na.rm=T) - (3*(sd(data1$DiffScores, na.rm=T)))
Diffscoremax <- mean(data1$DiffScores, na.rm=T) + (3*(sd(data1$DiffScores, na.rm=T)))
data1$DiffScores[which(data1$DiffScores < Diffscoremin | data1$DiffScores > Diffscoremax)] # 12 outliers
## [1] 1.60 -1.40 -1.30 -1.30 1.55 -1.25 1.40 1.75 1.45 1.60 2.05 2.00
data1$Diffscore.win<- winsorize(data1$DiffScores, method="zscore",
threshold=3, robust=TRUE)
mean(data1$DiffScores)
## [1] 0.1669792
mean(data1$Diffscore.win)
## [1] 0.1619322
Lik.nonSV.min <- mean(data1$Likability.nonSV, na.rm=T) - (3*(sd(data1$Likability.nonSV, na.rm=T)))
Lik.nonSV.max <- mean(data1$Likability.nonSV, na.rm=T) + (3*(sd(data1$Likability.nonSV, na.rm=T)))
data1$Likability.nonSV[which(data1$Likability.nonSV < Lik.nonSV.min | data1$Likability.nonSV > Lik.nonSV.max)] # 5 outliers
## [1] 1.6 1.5 1.4 5.0 1.6
data1$Likability.nonSV.win<- winsorize(data1$Likability.nonSV, method="zscore",
threshold=3, robust=TRUE)
data1$CSA.likability.win<- winsorize(data1$CSA.likability, method="zscore",
threshold=3, robust=TRUE)
data1$ASA.likability.win<- winsorize(data1$ASA.likability, method="zscore",
threshold=3, robust=TRUE)
Missingness
table(data1$lec5_08a)
##
## 0 Happened_to_me
## 604 156
table(data1$lec5_09a)
##
## 0 Happened_to_me
## 455 305
table(data1$Gender)
##
## Male Female
## 175 573
## Genderqueer/non-conforming Transgender
## 5 3
## Other (please specify)
## 4
table(data1$apref_e2)
##
## 1 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5
## 6 5 6 9 10 35 44 60 76 119 69 89 77 66 36 26 16 5 4 2
table(data1$apref_e5)
##
## 1 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5
## 5 2 2 11 16 16 27 56 80 128 98 85 87 60 41 18 12 8 3 4
sum(is.na(data1$lec5_08a))
## [1] 0
sum(is.na(data1$lec5_09a))
## [1] 0
sum(is.na(data1$Gender))
## [1] 0
sum(is.na(data1$lec5_08a))
## [1] 0
sum(is.na(data1$apref_e2))
## [1] 0
sum(is.na(data1$apref_e5))
## [1] 1
sum(is.na(data1$DiffScores))
## [1] 0
sum(is.na(data1$Age))
## [1] 1
sum(is.na(data1$Likability.nonSV))
## [1] 0
# 1 missing value for csa likability rating. NAs in LEC variables indicate that a given experience did NOT occur
Bivariate Correlations
data1$gender.num<- as.numeric(data1$gender.bin)
data1$SV.num<- as.numeric(data1$SV)
library(Hmisc)
## Registered S3 methods overwritten by 'Hmisc':
## method from
## [.labelled expss
## print.labelled expss
## as.data.frame.labelled expss
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
##
## src, summarize
## The following objects are masked from 'package:base':
##
## format.pval, units
Study1.cor<- data1 %>%
select(Likability.win, Likability.nonSV.win, SV.num, gender.num, Age)
(S1.corr<- rcorr(as.matrix(Study1.cor)))
## Likability.win Likability.nonSV.win SV.num gender.num
## Likability.win 1.00 0.74 0.09 0.06
## Likability.nonSV.win 0.74 1.00 -0.02 -0.01
## SV.num 0.09 -0.02 1.00 0.30
## gender.num 0.06 -0.01 0.30 1.00
## Age -0.09 -0.05 0.00 0.02
## Age
## Likability.win -0.09
## Likability.nonSV.win -0.05
## SV.num 0.00
## gender.num 0.02
## Age 1.00
##
## n
## Likability.win Likability.nonSV.win SV.num gender.num Age
## Likability.win 760 760 760 760 759
## Likability.nonSV.win 760 760 760 760 759
## SV.num 760 760 760 760 759
## gender.num 760 760 760 760 759
## Age 759 759 759 759 759
##
## P
## Likability.win Likability.nonSV.win SV.num gender.num
## Likability.win 0.0000 0.0137 0.1282
## Likability.nonSV.win 0.0000 0.6090 0.7171
## SV.num 0.0137 0.6090 0.0000
## gender.num 0.1282 0.7171 0.0000
## Age 0.0096 0.2042 0.9744 0.5591
## Age
## Likability.win 0.0096
## Likability.nonSV.win 0.2042
## SV.num 0.9744
## gender.num 0.5591
## Age
print(S1.corr$r, digits=3) #correlations with 3 decimal places
## Likability.win Likability.nonSV.win SV.num gender.num
## Likability.win 1.0000 0.7416 0.08937 0.0552
## Likability.nonSV.win 0.7416 1.0000 -0.01858 -0.0132
## SV.num 0.0894 -0.0186 1.00000 0.2959
## gender.num 0.0552 -0.0132 0.29587 1.0000
## Age -0.0940 -0.0461 0.00117 0.0212
## Age
## Likability.win -0.09398
## Likability.nonSV.win -0.04614
## SV.num 0.00117
## gender.num 0.02123
## Age 1.00000
# Paired samples t test comparing SV vignette ratings to non-SV ratings (e.g., car accident)
t.test(data1$Likability, data1$Likability.nonSV, paired=T, alternative="two.sided")
##
## Paired t-test
##
## data: data1$Likability and data1$Likability.nonSV
## t = -11.41, df = 759, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.1957068 -0.1382515
## sample estimates:
## mean difference
## -0.1669792
#data1$DiffScores <- (data1$Likability.nonSV - data1$Likability)
mean(data1$DiffScores, na.rm=T) / sd(data1$DiffScores, na.rm=T) # Cohen's d = 0.413901
## [1] 0.413901
ASA vs CSA Likability
vig.long.1<- data1 %>%
select(ID, apref_e1, apref_e2, apref_e3, apref_e4, apref_e5, apref_e6)
# apref_e1 = friend car accident, e2=ASA, e3=car accident, e4=hurricane survivor, e5=CSA, e6=childhood leukemia
# Convert to long format
vig.long.1<- vig.long.1 %>%
gather(key="Vignette", value="Likability", apref_e1, apref_e2, apref_e3, apref_e4, apref_e5, apref_e6) %>%
convert_as_factor(ID, Vignette)
# Summary statistics
vig.long.1 %>%
group_by(Vignette) %>%
get_summary_stats(Likability, type="mean_sd")
## # A tibble: 6 × 5
## Vignette variable n mean sd
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 apref_e1 Likability 759 3.17 0.633
## 2 apref_e2 Likability 760 3.15 0.66
## 3 apref_e3 Likability 760 3.33 0.648
## 4 apref_e4 Likability 760 3.37 0.603
## 5 apref_e5 Likability 759 3.19 0.617
## 6 apref_e6 Likability 759 3.48 0.641
# Repeated measures Anova
aov.1<- anova_test(data=vig.long.1, dv=Likability, wid=ID, within=Vignette)
## Warning: NA detected in rows: 86,3353,4489.
## Removing this rows before the analysis.
get_anova_table(aov.1)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Vignette 4.71 3557.06 81.946 2.9e-77 * 0.035
# Pairwise coparisons
pwc.1<- vig.long.1 %>%
pairwise_t_test(
Likability~Vignette, paired=TRUE,
p.adjust.method="bonferroni"
)
pwc.1
## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Likability apref_e1 apref_e2 759 760 0.878 758 3.8 e- 1 1 e+ 0
## 2 Likability apref_e1 apref_e3 759 760 -7.77 758 2.57e-14 3.86e-13
## 3 Likability apref_e1 apref_e4 759 760 -9.69 758 5.28e-21 7.92e-20
## 4 Likability apref_e1 apref_e5 759 759 -1.06 757 2.9 e- 1 1 e+ 0
## 5 Likability apref_e1 apref_e6 759 759 -14.1 757 2.36e-40 3.54e-39
## 6 Likability apref_e2 apref_e3 760 760 -8.10 759 2.11e-15 3.16e-14
## 7 Likability apref_e2 apref_e4 760 760 -10.5 759 3.98e-24 5.97e-23
## 8 Likability apref_e2 apref_e5 760 759 -2.23 758 2.6 e- 2 3.87e- 1
## 9 Likability apref_e2 apref_e6 760 759 -14.8 758 8.72e-44 1.31e-42
## 10 Likability apref_e3 apref_e4 760 760 -2.31 759 2.1 e- 2 3.21e- 1
## 11 Likability apref_e3 apref_e5 760 759 6.50 758 1.43e-10 2.14e- 9
## 12 Likability apref_e3 apref_e6 760 759 -7.73 758 3.5 e-14 5.25e-13
## 13 Likability apref_e4 apref_e5 760 759 8.70 758 2.06e-17 3.09e-16
## 14 Likability apref_e4 apref_e6 760 759 -5.67 758 2.03e- 8 3.04e- 7
## 15 Likability apref_e5 apref_e6 759 759 -14.0 757 7.44e-40 1.12e-38
## # ℹ 1 more variable: p.adj.signif <chr>
Multiple Regression
# Dummy coding: On 5/31/24 we decided to use dummy coding instead of deviation coding using UMN as reference group
contrasts(data1$ASA.ending.fac) = contr.treatment(3)
# I only controlled for ASA ending because for Study 1 participants received same ending for both vignettes
#Study1.reg<- lm(Likability.n1 ~ ASA.ending.fac + Age + gender.bin + SV + gender.bin*SV, data=data1)
#summary(Study1.reg)
#confint(Study1.reg)
library(lm.beta)
#Study1.stan<- lm.beta(Study1.reg)
#Study1.stan
# On 7/23/24 we decided to not control for ending because they are randomized across participants and to keep results consistent with the mixed effect models. We note this deviation from our preregistered data analysis plan in the manuscript.
Study1.reg.lim<- lm(Likability.n1 ~ Age + gender.bin + SV + gender.bin*SV, data=data1)
summary(Study1.reg.lim)
##
## Call:
## lm(formula = Likability.n1 ~ Age + gender.bin + SV + gender.bin *
## SV, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.03740 -0.37575 -0.03493 0.39394 1.98300
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.206400 0.060029 53.415 < 2e-16 ***
## Age -0.003574 0.001343 -2.662 0.00794 **
## gender.bin1 0.027554 0.059181 0.466 0.64165
## SV1 0.021216 0.120371 0.176 0.86014
## gender.bin1:SV1 0.084900 0.129729 0.654 0.51302
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5834 on 754 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.01871, Adjusted R-squared: 0.01351
## F-statistic: 3.595 on 4 and 754 DF, p-value: 0.006504
round(confint(Study1.reg.lim), 3)
## 2.5 % 97.5 %
## (Intercept) 3.089 3.324
## Age -0.006 -0.001
## gender.bin1 -0.089 0.144
## SV1 -0.215 0.258
## gender.bin1:SV1 -0.170 0.340
lm.beta(Study1.reg.lim)
##
## Call:
## lm(formula = Likability.n1 ~ Age + gender.bin + SV + gender.bin *
## SV, data = data1)
##
## Standardized Coefficients::
## (Intercept) Age gender.bin1 SV1 gender.bin1:SV1
## NA -0.09611576 0.01976963 0.01787646 0.07054016
vif(Study1.reg.lim, type=c("predictor"))
## GVIFs computed for predictors
## GVIF Df GVIF^(1/(2*Df)) Interacts With Other Predictors
## Age 1.001929 1 1.000964 -- gender.bin, SV
## gender.bin 1.001929 3 1.000321 SV Age
## SV 1.001929 3 1.000321 gender.bin Age
Assumption Checks
plot(Study1.reg.lim)
Exploratory Analyses
Mixed model
S1.mixeddata<- data1 %>%
select(ID, Age, gender.bin, SV, Likability.win, Likability.nonSV.win)
mixed1_df<-NULL # 1.) construct empty (null) dataframe
for (row_i in 1:nrow(S1.mixeddata)) { # 2.) loop thru each row in S1.Mixeddata
for (vig in 0:1) { # 3.) loop twice , 0 and 1 (for each row of S1.Mixed data, we wil generate 2 rows in mixed1_df)
id = S1.mixeddata[[row_i,1]]
age = S1.mixeddata[[row_i,2]]
gen = S1.mixeddata[[row_i,3]]
sv = S1.mixeddata[[row_i,4]]
if (vig == 0) {
like = S1.mixeddata[[row_i,6]]
} else {
like = S1.mixeddata[[row_i,5]]
}
mixed1_temp_df <- data.frame(ID=id,
Age=age,
gender.bin=gen,
SV=sv,
vig=vig, # 0 = non-SV likability ratings; 1 = SV likability ratings
Likability=like) # create a temporary dataframe with the latest values, and append temp frame onto mixed1_df
rbind(mixed1_df,mixed1_temp_df)->mixed1_df
}
}
# Computing intraclass correlation coefficient (ICC)
mixed1_df$ID<- as.factor(mixed1_df$ID) #making ID a factor variable
library(ICC)
(ICC.1 <- ICCbare(mixed1_df$ID, mixed1_df$Likability, data = mixed1_df)) # ICC = .70
## [1] 0.6994779
################## Visualizing Data ###################
x1<- ggplot(mixed1_df, aes(x = vig, y = Likability, colour=gender.bin, group=gender.bin)) + scale_x_continuous(breaks=c(0, 1), name="Vignette", labels=c("Non-SV Vignettes", "SV Vignettes")) + scale_color_hue(name= "Gender", labels=c("Cisgender Men", "Women/TNB+")) + ggtitle("Figure 1") +
stat_summary(fun="mean",geom="line")
x2<- x1 + ylim(1, 5)
x2
# + theme(text=element_text(size=12, family="serif")) ## code to change font to times new roman
a1<- ggplot(mixed1_df, aes(x = vig, y = Likability, colour=SV,group=SV)) + scale_x_continuous(breaks=c(0, 1), name="Vignette", labels=c("Non-SV Vignettes", "SV Vignettes"))+ scale_color_hue(name= "SV History", labels=c("Non-victim", "Victim")) + ggtitle("Figure 2") +
stat_summary(fun="mean",geom="line")
a2<- a1 + ylim(1, 5)
a2
b1<- ggplot(mixed1_df, aes(x = vig, y = Likability)) + scale_x_continuous(breaks=c(0, 1), name="Vignette")+ ggtitle("Likability Scores Across Vignettes by SV History and Gender") +
stat_summary(fun="mean",geom="line")
b1 + facet_grid(gender.bin~SV, labeller=label_both) + ylim(1, 5)
#mixed1_df$Age.3 <- cut_number(mixed1_df$Age,n=3) #create groups, originally plotted at +-1 SD M but realized that doesn't really make sense because age is very right skewed. Cut_number cuts a numeric vector into intervals containing equal number of points
#levels(mixed1_df$Age.3) <- c("Low","Med","High") # name groups
mixed1_df.c <- na.omit(mixed1_df) # creating complete dataset because one value missing for age
mixed1_df.c$Age.bin<- cut(mixed1_df.c$Age, breaks=c(0, 43.62, 100), labels=c("Mean Age", "+1 SD")) #creating age bins that are mean, 1 sd above mean
d2<- ggplot(mixed1_df.c, aes(x = vig, y = Likability, colour=Age.bin,group=Age.bin)) + scale_x_continuous(breaks=c(0, 1), name="Vignette", labels=c("Non-SV Vignettes", "SV Vignettes")) + labs(color = "Age") + ggtitle("Figure 3") +
stat_summary(fun="mean",geom="line")
d3<- d2 + ylim(1, 5)
d3
################ Random Effect Structure #########################
#mod.1<- lmer(Likability ~ 1 + vig + (1 + vig | ID), data=mixed1_df, REML=TRUE) # Random intercept, random slope, model does not converge
mod.2<- lme4::lmer(Likability ~ 1 + vig + (1 | ID), data=mixed1_df, REML=TRUE)# Random intercept, fixed slope
#mod.3<- lmer(Likability ~ 1 + vig + (-1 | ID), data=mixed1_df, REML=TRUE) # Random slope, fixed intercept, does not converge
mod.4<- lm(Likability ~ 1 + vig, data=mixed1_df) # Fixed intercept, fixed slope
AIC(mod.2) # So the random intercept model better fits the data than the model with fixed intercept and slope (so regular regression)
## [1] 1949.723
AIC(mod.4)
## [1] 2532.102
######################### Error Structure ######################
# Independent error structure
ind.error <- nlme::lme(Likability ~ 1 + vig, # Fixed portion of the model
random = ~ 1 | ID, # Random portion of the model
data = mixed1_df)
nlme::ACF(ind.error) # If independent error structure is plausible, we expect autocorrelations to be about 0 (this does not seem to be the case)
## lag ACF
## 1 0 1.0000000
## 2 1 -0.7376897
# Compound Symmetry
comp.error <- nlme::lme(Likability ~ 1 + vig, # Fixed portion of the model
random = ~ 1 | ID, # Random portion of the model
data = mixed1_df,
correlation=nlme::corCompSymm(.3),
control=lme4::lmerControl(check.conv.singular = lme4::.makeCC(action = "ignore", tol = 1e-3)))
# FIRST ORDER AUTOREGRESSIVE STRUCTURE
auto.reg <- nlme::lme(Likability ~ 1 + vig, # Fixed portion of the model
random = ~ 1 | ID, # Random portion of the model
data = mixed1_df,
correlation=nlme::corAR1())
# Unstructured
unstructured <- nlme::lme(Likability ~ 1 + vig,
random = ~ 1 | ID,
data = mixed1_df,
correlation=nlme::corSymm(),
weights = nlme::varIdent(form = ~ 1 | vig)) # unstructured error variances provide the most flexibility but contain more parameters
AIC(ind.error)
## [1] 1949.723
AIC(comp.error)
## [1] 1951.723
AIC(auto.reg)
## [1] 1951.723
AIC(unstructured) # So the model with the unstructured error variance provides the best fit to the data
## [1] 1935.838
########################### FULL MODEL ###################
full.mod.1 <- nlme::lme(Likability ~ vig*Age + vig*gender.bin + vig*SV,
random = ~ 1 | ID,
data = mixed1_df.c,
correlation=nlme::corSymm(),
weights = nlme::varIdent(form = ~ 1 | vig))
summary(full.mod.1)
## Linear mixed-effects model fit by REML
## Data: mixed1_df.c
## AIC BIC logLik
## 1960.594 2024.433 -968.2971
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.4695287 0.2395552
##
## Correlation Structure: General
## Formula: ~1 | ID
## Parameter estimate(s):
## Correlation:
## 1
## 2 0.096
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | vig
## Parameter estimates:
## 0 1
## 1.000000 1.419897
## Fixed effects: Likability ~ vig * Age + vig * gender.bin + vig * SV
## Value Std.Error DF t-value p-value
## (Intercept) 3.396567 0.05214705 755 65.13441 0.0000
## vig -0.202920 0.03925738 755 -5.16895 0.0000
## Age -0.001532 0.00121208 755 -1.26425 0.2065
## gender.bin1 -0.009269 0.04755767 755 -0.19489 0.8455
## SV1 -0.017746 0.04048846 755 -0.43830 0.6613
## vig:Age -0.001970 0.00091248 755 -2.15909 0.0312
## vig:gender.bin1 0.055672 0.03580241 755 1.55499 0.1204
## vig:SV1 0.110142 0.03048056 755 3.61350 0.0003
## Correlation:
## (Intr) vig Age gndr.1 SV1 vig:Ag vg:g.1
## vig -0.237
## Age -0.633 0.150
## gender.bin1 -0.590 0.140 -0.022
## SV1 -0.128 0.030 0.005 -0.295
## vig:Age 0.150 -0.633 -0.237 0.005 -0.001
## vig:gender.bin1 0.140 -0.590 0.005 -0.237 0.070 -0.022
## vig:SV1 0.030 -0.128 -0.001 0.070 -0.237 0.005 -0.295
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.6711127063 -0.4208341120 0.0004137708 0.4215568136 3.0979926028
##
## Number of Observations: 1518
## Number of Groups: 759
AIC(full.mod.1) # AIC
## [1] 1960.594
r2glmm::r2beta(full.mod.1) # R2 square
## Effect Rsq upper.CL lower.CL
## 1 Model 0.049 0.075 0.033
## 2 vig 0.007 0.017 0.001
## 8 vig:SV1 0.003 0.012 0.000
## 3 Age 0.001 0.008 0.000
## 6 vig:Age 0.001 0.007 0.000
## 7 vig:gender.bin1 0.001 0.006 0.000
## 5 SV1 0.000 0.004 0.000
## 4 gender.bin1 0.000 0.003 0.000
nlme::intervals(full.mod.1, which="fixed") # Confidence intervals
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 3.294196764 3.396567203 3.4989376427
## vig -0.279986152 -0.202919551 -0.1258529504
## Age -0.003911823 -0.001532375 0.0008470739
## gender.bin1 -0.102629682 -0.009268696 0.0840922895
## SV1 -0.097229585 -0.017746246 0.0617370924
## vig:Age -0.003761425 -0.001970127 -0.0001788282
## vig:gender.bin1 -0.014611643 0.055672453 0.1259565492
## vig:SV1 0.050304873 0.110141587 0.1699783014
# For more informations, see https://stats.oarc.ucla.edu/r/seminars/interactions-r/#s4b
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
emvigSV <- emmeans(full.mod.1, ~ vig*SV) #calculate predicted means for SV by vig
contrast(emvigSV, "revpairwise",by="SV") # simple slopes
## SV = -1:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.23 0.0199 755 -11.557 <.0001
##
## SV = 1:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.12 0.0266 755 -4.509 <.0001
##
## Results are averaged over the levels of: gender.bin
## Degrees-of-freedom method: containment
emvigGen<- emmeans(full.mod.1, ~vig*gender.bin) #calculate predicted means for gender by vig
contrast(emvigGen, "revpairwise",by="gender.bin")
## gender.bin = -1:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.203 0.0317 755 -6.384 <.0001
##
## gender.bin = 1:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.147 0.0164 755 -8.950 <.0001
##
## Results are averaged over the levels of: SV
## Degrees-of-freedom method: containment
mylist.1<- list(Age=c(28.82, 43.62), vig=c(1, 0))
emcontcat.1<- emmeans(full.mod.1, ~vig*Age, at=mylist.1)
contrast(emcontcat.1, "pairwise", by="Age")
## Age = 28.8:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.177 0.0179 755 -9.892 <.0001
##
## Age = 43.6:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.206 0.0231 755 -8.925 <.0001
##
## Results are averaged over the levels of: gender.bin, SV
## Degrees-of-freedom method: containment
Study 2
Data preparation
data2<- data.2 %>%
select(ID, event2_end, event5_end, lec5_08a, lec5_09a, Gender, apref_e2, apref_e5, Sample, Age, apref_e1, apref_e3, apref_e4, apref_e6)
################# GENDER #########################
table(data2$Gender) # Gender frequencies
##
## Male Female
## 395 706
## Genderqueer/non-conforming Transgender
## 7 4
## Other (please specify)
## 0
data2<- data2 %>% #Created new binary gender variable for men vs. women/TNB+
mutate(gender.bin = case_when(
data2$Gender=="Male" ~ -1, #Man
data2$Gender=="Female" ~ 1, #WTNB+
data2$Gender=="Genderqueer/non-conforming" ~ 1,#WTNB+
data2$Gender=="Transgender" ~ 1, #WTNB+
data2$Gender=="Other (please specify)" ~ 1, #WTNB+
))
table(data2$gender.bin) #checking new frequencies, looks correct, __of sample is WTNB+
##
## -1 1
## 395 717
prop.table(table(data2$gender.bin))
##
## -1 1
## 0.3552158 0.6447842
data2$gender.bin<- as.factor(data2$gender.bin) # Converting gender to factor
################# SAMPLE #########################
data2<- data2 %>%
mutate(sample.fac = case_when(
data2$Sample=="UMN" ~ 1,
data2$Sample=="WWU" ~ 2,
data2$Sample=="mturk" ~ 3))
data2$sample.fac<- as.factor(data2$sample.fac) # Converting sample to factor
################# ASA ENDING #########################
data2<- data2 %>%
mutate(ASA.ending.fac = case_when(
data2$event2_end=="Negative" ~ 1,
data2$event2_end=="Redemption" ~ 2,
data2$event2_end=="Survivor" ~ 3))
data2$ASA.ending.fac<- as.factor(data2$ASA.ending.fac) # Converting ASA ending to factor
contrasts(data2$ASA.ending.fac) = contr.sum(3)
################# CSA ENDING #########################
data2<- data2 %>%
mutate(CSA.ending.fac = case_when(
data2$event5_end=="Negative" ~ 1,
data2$event5_end=="Redemption" ~ 2,
data2$event5_end=="Survivor" ~ 3))
data2$CSA.ending.fac<- as.factor(data2$CSA.ending.fac) # Converting CSA ending to factor
contrasts(data2$CSA.ending.fac) = contr.sum(3)
################# Participant SV History #########################
data2$lec5_08a<- as.character(data2$lec5_08a) # First convert to character, then switch back to factor to avoid R error
data2$lec5_09a<- as.character(data2$lec5_09a)
# Replace NAs with 0s (Nas in this case indicate that experience does not apply)
data2$lec5_08a[is.na(data2$lec5_08a)] <- 0
data2$lec5_09a[is.na(data2$lec5_09a)] <- 0
# There were some oddities with spacing (some responses Happened_to_me some happened to me) so creating new binary variables
data2$lec8_bin <- ifelse(data2$lec5_08a==0, 0, 1) # so if 9a is 0 (did not happen to me) then 0 else 1 (happened to me)
data2$lec9_bin <- ifelse(data2$lec5_09a==0, 0, 1) # so if 9a is 0 (did not happen to me) then 0 else 1 (happened to me)
data2<- data2 %>%
mutate(SV = case_when(
data2$lec8_bin==1 | data2$lec9_bin==1 ~ 1, #so if experienced SA or other unwanted sexual experience, then 1
data2$lec8_bin==0 & data2$lec9_bin==0 ~ -1 # if did not experience SA AND did not experience unwanted experience then 0
))
svcheck<- data2 %>% #creating dataframe with SV variables to check that they're coded correctly
select(lec5_08a, lec5_09a, lec8_bin, lec9_bin, SV)
data2$SV<- as.factor(data2$SV)
table(data2$SV)
##
## -1 1
## 723 389
data2$Age<- as.numeric(data2$Age)
data2 <- data2 %>% #removing 13 year old from analyses but this person is still included in total N and demographic information
dplyr::filter(Age!=13)
Likability Ratings
# Averaging ASA and CSA likability ratings, as previous research using this dataset found no mean differences
data2 <- data2 %>%
rowwise() %>%
mutate(Likability = mean(c(apref_e2, apref_e5), na.rm = TRUE))%>%
ungroup()
# Checking that these variables are the type I want, all factor except likability which is numeric so good
class(data2$gender.bin)
## [1] "factor"
class(data2$sample.fac)
## [1] "factor"
class(data2$ASA.ending.fac)
## [1] "factor"
class(data2$CSA.ending.fac)
## [1] "factor"
class(data2$SV)
## [1] "factor"
class(data2$Likability)
## [1] "numeric"
# Average Non-SV likability
data2 <- data2 %>%
rowwise() %>%
mutate(Likability.nonSV = mean(c(apref_e1, apref_e3, apref_e4, apref_e6), na.rm = TRUE))%>%
ungroup()
# Difference scores: SV and non-SV likability
data2$DiffScores <- (data2$Likability.nonSV - data2$Likability)
# Positive scores mean that participants rated narrators in the non-SV vignettes as more likable. Negative scores mean that participants rated victims of sexual violence as more likable.
Data visualization
Likability.mean <- mean(data2$Likability)
hist(data2$Likability, main="Study 2 Histogram of Likability", xlab="Likability", col='hotpink', sub=paste("Skewness:",
round(e1071::skewness(data2$Likability, na.rm=TRUE), 2)))
abline(v = Likability.mean, col = 'orange') # Overlay mean on histogram
qqnorm(data2$Likability, pch = 1, frame = FALSE, main="QQ Plot of Likability")
qqline(data2$Likability, col = "hotpink", lwd = 2)
# Means and sd
mean(data2$Likability, na.rm=T)
## [1] 3.19811
sd(data2$Likability, na.rm=T)
## [1] 0.5576221
library(ggplot2)
ggplot(data2, aes(SV, Likability)) +
geom_boxplot(fill = "lightpink", color = "hotpink") +
labs(x="Assault History (No SV = -1; SV = 1)", y="Likability") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggplot(data2, aes(gender.bin, Likability)) +
geom_boxplot(fill = "lightpink", color = "hotpink") +
labs(x="Gender (Male = -1; WTNB+ = 1)", y="Likability") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggplot(data2, aes(sample.fac, Likability)) +
geom_boxplot(fill = "lightpink", color = "hotpink") +
labs(x="Sample (UMN = 1; WWU = 2; MTurk = 3)", y="Likability") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggplot(data2, aes(ASA.ending.fac, Likability)) + # For this study they had same endings for CSA and ASA conditions
geom_boxplot(fill = "lightpink", color = "hotpink") +
labs(x="Story Ending (Negative = 1; Redemptive = 2; Survivor = 3)", y="Likability") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
################### AGE ##########################
hist(data2$Age)
scatterplot(Likability ~ Age, data=data2,
xlab="Age", ylab="Likability",
main="Age by Likability", col="hotpink")
#################### DIFFERENCE SCORES ################
Diffscore2.mean <- mean(data2$DiffScores)
hist(data2$DiffScores, main="Histogram of Difference Scores", xlab="Difference Scores (Likability.nonSV - Likability", col='hotpink', sub=paste("Skewness:",
round(e1071::skewness(data2$DiffScores, na.rm=TRUE), 2)))
abline(v = Diffscore2.mean, col = 'orange') # Overlay mean on histogram
qqnorm(data2$DiffScores, pch = 1, frame = FALSE, main="QQ Plot of Difference SCores")
qqline(data2$DiffScores, col = "hotpink", lwd = 2)
Descriptives
################# Sexual Violence (SV) History ####################
cross_cases(data2, lec8_bin, lec9_bin) # lec8 is sexual assault, lec9 is unwanted sexual experiences
|
|
 lec9_binÂ
|
|
|
 0Â
|
 1Â
|
|
 lec8_binÂ
|
|
   0Â
|
722
|
223
|
|
   1Â
|
18
|
148
|
|
   #Total casesÂ
|
740
|
371
|
################# SV History by Likability #############
group_by(data2, SV) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 2 × 6
## SV count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -1 722 3.17 0.562 3.1 0.6
## 2 1 389 3.25 0.547 3.3 0.7
################# Gender by Likability #############
group_by(data2, gender.bin) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 2 × 6
## gender.bin count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -1 394 3.11 0.579 3.1 0.7
## 2 1 717 3.25 0.540 3.2 0.7
################# Sample by Likability #############
group_by(data2, Sample) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 3 × 6
## Sample count mean sd median IQR
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 UMN 391 3.24 0.523 3.2 0.700
## 2 WWU 384 3.20 0.487 3.2 0.600
## 3 mturk 336 3.15 0.661 3.1 0.825
################# Ending by Likability #############
group_by(data2, event2_end) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 3 × 6
## event2_end count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Redemption 369 3.27 0.583 3.3 0.7
## 2 Negative 369 3.01 0.479 3 0.600
## 3 Survivor 373 3.31 0.557 3.3 0.7
group_by(data2, event5_end) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 3 × 6
## event5_end count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Redemption 361 3.30 0.539 3.3 0.7
## 2 Negative 375 2.97 0.532 3 0.600
## 3 Survivor 375 3.32 0.531 3.3 0.7
mean(data2$Likability)
## [1] 3.19811
median(data2$Likability)
## [1] 3.2
sd(data2$Likability)
## [1] 0.5576221
################# SV by Gender #############
cross_cases(data2, SV, gender.bin)
|
|
 gender.binÂ
|
|
|
 -1Â
|
 1Â
|
|
 SVÂ
|
|
   -1Â
|
342
|
380
|
|
   1Â
|
52
|
337
|
|
   #Total casesÂ
|
394
|
717
|
cross_cases(data2, lec8_bin, gender.bin) # sexual assault
|
|
 gender.binÂ
|
|
|
 -1Â
|
 1Â
|
|
 lec8_binÂ
|
|
   0Â
|
380
|
565
|
|
   1Â
|
14
|
152
|
|
   #Total casesÂ
|
394
|
717
|
cross_cases(data2, lec9_bin, gender.bin) # other unwanted or uncomfortable sexual experiences
|
|
 gender.binÂ
|
|
|
 -1Â
|
 1Â
|
|
 lec9_binÂ
|
|
   0Â
|
344
|
396
|
|
   1Â
|
50
|
321
|
|
   #Total casesÂ
|
394
|
717
|
data2<- data2 %>% #Created new interaction term for descriptives that includes SV history x gender
mutate(SVxGen = case_when(
data2$gender.bin==-1 & data2$SV==-1 ~ "M_no", # Men with no SV history
data2$gender.bin==-1 & data2$SV==1 ~ "M_yes", # Men with SV history
data2$gender.bin==1 & data2$SV==-1 ~ "WTNB_no", #WTNB+ with no SV history
data2$gender.bin==1 & data2$SV==1 ~ "WTNB_yes", #WTNB+ with SV history
))
table(data2$SVxGen) # Frequencies by gender and SV history
##
## M_no M_yes WTNB_no WTNB_yes
## 342 52 380 337
group_by(data2, SVxGen) %>%
summarise(
count = n(),
mean = mean(Likability, na.rm = TRUE),
sd = sd(Likability, na.rm = TRUE),
median = median(Likability, na.rm = TRUE),
IQR = IQR(Likability, na.rm = TRUE)
)
## # A tibble: 4 × 6
## SVxGen count mean sd median IQR
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 M_no 342 3.10 0.577 3.1 0.7
## 2 M_yes 52 3.16 0.591 3.2 0.65
## 3 WTNB_no 380 3.23 0.541 3.2 0.7
## 4 WTNB_yes 337 3.26 0.539 3.3 0.7
################## AGE ######################
mean(data2$Age, na.rm=T)
## [1] 26.32043
range(data2$Age, na.rm=T)
## [1] 18 71
sd(data2$Age, na.rm=T)
## [1] 11.85072
group_by(data2, Sample) %>%
summarise(
count = n(),
mean = mean(Age, na.rm = TRUE),
sd = sd(Age, na.rm = TRUE),
median = median(Age, na.rm = TRUE),
IQR = IQR(Age, na.rm = TRUE)
)
## # A tibble: 3 × 6
## Sample count mean sd median IQR
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 UMN 391 19.8 2.11 19 2
## 2 WWU 384 19.7 2.33 19 2
## 3 mturk 336 41.4 11.3 39 17
S2.samplexage<- aov(Age ~ Sample, data=data2)
summary(S2.samplexage)
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 2 109299 54649 1300 <2e-16 ***
## Residuals 1108 46589 42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(S2.samplexage)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Age ~ Sample, data = data2)
##
## $Sample
## diff lwr upr p adj
## UMN-mturk -21.5348237 -22.666885 -20.4027619 0.0000000
## WWU-mturk -21.6547619 -22.791583 -20.5179413 0.0000000
## WWU-UMN -0.1199382 -1.213285 0.9734083 0.9641243
eta_squared(S2.samplexage)
## Sample
## 0.7011353
################# SV History by Difference Scores #############
group_by(data2, SV) %>%
summarise(
count = n(),
mean = mean(DiffScores, na.rm = TRUE),
sd = sd(DiffScores, na.rm = TRUE),
median = median(DiffScores, na.rm = TRUE),
IQR = IQR(DiffScores, na.rm = TRUE)
)
## # A tibble: 2 × 6
## SV count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -1 722 0.235 0.516 0.200 0.65
## 2 1 389 0.125 0.536 0.100 0.7
################# Gender by Difference Scores #############
group_by(data2, gender.bin) %>%
summarise(
count = n(),
mean = mean(DiffScores, na.rm = TRUE),
sd = sd(DiffScores, na.rm = TRUE),
median = median(DiffScores, na.rm = TRUE),
IQR = IQR(DiffScores, na.rm = TRUE)
)
## # A tibble: 2 × 6
## gender.bin count mean sd median IQR
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -1 394 0.241 0.542 0.200 0.65
## 2 1 717 0.172 0.515 0.150 0.6
mean(data2$Likability.nonSV)
## [1] 3.394993
sd(data2$Likability.nonSV)
## [1] 0.4560459
Outliers
Likabilitymin2 <- mean(data2$Likability, na.rm=T) - (3*(sd(data2$Likability, na.rm=T)))
Likabilitymax2 <- mean(data2$Likability, na.rm=T) + (3*(sd(data2$Likability, na.rm=T)))
data2$Likability[which(data2$Likability < Likabilitymin2 | data2$Likability > Likabilitymax2)] # 8 outliers
## [1] 1.2 1.0 1.0 1.5 4.9 5.0 1.2 1.0
data2$Likability.win<- winsorize(data2$Likability, method="zscore",
threshold=3, robust=TRUE) # Robust= TRUE means that values are winsorized based on their Median Absolute Deviation (MAD)
Diffscoremin2 <- mean(data2$DiffScores, na.rm=T) - (3*(sd(data2$DiffScores, na.rm=T)))
Diffscoremax2 <- mean(data2$DiffScores, na.rm=T) + (3*(sd(data2$DiffScores, na.rm=T)))
data2$DiffScores[which(data2$DiffScores < Diffscoremin2 | data2$DiffScores > Diffscoremax2)] #7 outliers
## [1] -1.50 1.95 1.80 1.80 2.30 -1.95 1.95
data2$Diffscore.win<- winsorize(data2$DiffScores, method="zscore",
threshold=3, robust=TRUE)
mean(data2$DiffScores)
## [1] 0.1968834
mean(data2$Diffscore.win)
## [1] 0.1937494
Lik.nonSV.min.2 <- mean(data2$Likability.nonSV, na.rm=T) - (3*(sd(data2$Likability.nonSV, na.rm=T)))
Lik.nonSV.max.2 <- mean(data2$Likability.nonSV, na.rm=T) + (3*(sd(data2$Likability.nonSV, na.rm=T)))
data2$Likability.nonSV[which(data2$Likability.nonSV < Lik.nonSV.min.2 | data2$Likability.nonSV > Lik.nonSV.max.2)] # 8 outliers
## [1] 4.85 1.75 1.00 1.75 1.90 4.90 1.75 1.90
data2$Likability.nonSV.win<- winsorize(data2$Likability.nonSV, method="zscore",
threshold=3, robust=TRUE)
data2$CSA.likability.win<- winsorize(data2$apref_e5, method="zscore",
threshold=3, robust=TRUE)
data2$ASA.likability.win<- winsorize(data2$apref_e2, method="zscore",
threshold=3, robust=TRUE)
Missingness
table(data2$lec5_08a)
##
## 0 Happened to me\n Happened_to_me
## 945 69 97
table(data2$lec5_09a)
##
## 0 Happened to me\n Happened_to_me
## 740 160 211
table(data2$Gender)
##
## Male Female
## 394 706
## Genderqueer/non-conforming Transgender
## 7 4
## Other (please specify)
## 0
table(data2$SV)
##
## -1 1
## 722 389
table(data2$apref_e2)
##
## 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8
## 4 2 9 15 19 17 41 57 81 109 144 130 130 109 79 74 45 24 14 4
## 5
## 4
table(data2$apref_e5)
##
## 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8
## 6 1 5 4 9 23 38 38 101 97 159 168 107 95 84 80 42 25 15 7
## 5
## 7
sum(is.na(data2$lec5_08a))
## [1] 0
sum(is.na(data2$lec5_09a))
## [1] 0
sum(is.na(data2$Gender))
## [1] 0
sum(is.na(data2$SV))
## [1] 0
sum(is.na(data2$apref_e2))
## [1] 0
sum(is.na(data2$apref_e5))
## [1] 0
sum(is.na(data2$Diffscore.win))
## [1] 0
Bivariate Correlations
data2$gender.num<- as.numeric(data2$gender.bin)
data2$SV.num<- as.numeric(data2$SV)
Study2.cor<- data2 %>%
select(Likability.win, Likability.nonSV.win, SV.num, gender.num, Age)
(S2.corr<- rcorr(as.matrix(Study2.cor)))
## Likability.win Likability.nonSV.win SV.num gender.num
## Likability.win 1.00 0.46 0.06 0.11
## Likability.nonSV.win 0.46 1.00 -0.04 0.07
## SV.num 0.06 -0.04 1.00 0.34
## gender.num 0.11 0.07 0.34 1.00
## Age -0.08 0.10 -0.07 -0.06
## Age
## Likability.win -0.08
## Likability.nonSV.win 0.10
## SV.num -0.07
## gender.num -0.06
## Age 1.00
##
## n= 1111
##
##
## P
## Likability.win Likability.nonSV.win SV.num gender.num
## Likability.win 0.0000 0.0327 0.0002
## Likability.nonSV.win 0.0000 0.2309 0.0236
## SV.num 0.0327 0.2309 0.0000
## gender.num 0.0002 0.0236 0.0000
## Age 0.0084 0.0013 0.0238 0.0569
## Age
## Likability.win 0.0084
## Likability.nonSV.win 0.0013
## SV.num 0.0238
## gender.num 0.0569
## Age
print(S2.corr$r, digits=3) #correlations with 3 decimal places
## Likability.win Likability.nonSV.win SV.num gender.num
## Likability.win 1.0000 0.4613 0.0641 0.1122
## Likability.nonSV.win 0.4613 1.0000 -0.0360 0.0679
## SV.num 0.0641 -0.0360 1.0000 0.3390
## gender.num 0.1122 0.0679 0.3390 1.0000
## Age -0.0791 0.0965 -0.0678 -0.0571
## Age
## Likability.win -0.0791
## Likability.nonSV.win 0.0965
## SV.num -0.0678
## gender.num -0.0571
## Age 1.0000
t.test(data2$Likability, data2$Likability.nonSV, paired=T, alternative="two.sided")
##
## Paired t-test
##
## data: data2$Likability and data2$Likability.nonSV
## t = -12.496, df = 1110, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.2277975 -0.1659694
## sample estimates:
## mean difference
## -0.1968834
mean(data2$DiffScores, na.rm=T) / sd(data2$DiffScores, na.rm=T) # Cohen's d = 0.373319
## [1] 0.374902
ASA vs. CSA Likability
vig.long.2<- data2 %>%
select(ID, apref_e1, apref_e2, apref_e3, apref_e4, apref_e5, apref_e6)
# apref_e1 = friend car accident, e2=ASA, e3=car accident, e4=hurricane survivor, e5=CSA, e6=childhood leukemia
# Convert to long format
vig.long.2<- vig.long.2 %>%
gather(key="Vignette", value="Likability", apref_e1, apref_e2, apref_e3, apref_e4, apref_e5, apref_e6) %>%
convert_as_factor(ID, Vignette)
# Summary statistics
vig.long.2 %>%
group_by(Vignette) %>%
get_summary_stats(Likability, type="mean_sd")
## # A tibble: 6 × 5
## Vignette variable n mean sd
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 apref_e1 Likability 1111 3.20 0.68
## 2 apref_e2 Likability 1111 3.18 0.681
## 3 apref_e3 Likability 1111 3.42 0.691
## 4 apref_e4 Likability 1111 3.44 0.614
## 5 apref_e5 Likability 1111 3.22 0.657
## 6 apref_e6 Likability 1111 3.52 0.719
# Repeated measures Anova
aov.2<- anova_test(data=vig.long.2, dv=Likability, wid=ID, within=Vignette)
get_anova_table(aov.2)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Vignette 4.82 5355.05 71.438 3.14e-70 * 0.037
# Pairwise coparisons using Bonferroni correction
pwc.2<- vig.long.2 %>%
pairwise_t_test(
Likability~Vignette, paired=TRUE,
p.adjust.method="bonferroni"
)
pwc.2
## # A tibble: 15 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Likability apref_e1 apref_e2 1111 1111 1.08 1110 2.8 e- 1 1 e+ 0
## 2 Likability apref_e1 apref_e3 1111 1111 -8.48 1110 6.93e-17 1.04e-15
## 3 Likability apref_e1 apref_e4 1111 1111 -9.58 1110 5.86e-21 8.79e-20
## 4 Likability apref_e1 apref_e5 1111 1111 -0.566 1110 5.71e- 1 1 e+ 0
## 5 Likability apref_e1 apref_e6 1111 1111 -12.8 1110 5.05e-35 7.58e-34
## 6 Likability apref_e2 apref_e3 1111 1111 -9.36 1110 4.23e-20 6.35e-19
## 7 Likability apref_e2 apref_e4 1111 1111 -11.5 1110 7.12e-29 1.07e-27
## 8 Likability apref_e2 apref_e5 1111 1111 -1.81 1110 7.1 e- 2 1 e+ 0
## 9 Likability apref_e2 apref_e6 1111 1111 -13.0 1110 5.32e-36 7.98e-35
## 10 Likability apref_e3 apref_e4 1111 1111 -0.543 1110 5.87e- 1 1 e+ 0
## 11 Likability apref_e3 apref_e5 1111 1111 8.07 1110 1.85e-15 2.78e-14
## 12 Likability apref_e3 apref_e6 1111 1111 -3.69 1110 2.37e- 4 4 e- 3
## 13 Likability apref_e4 apref_e5 1111 1111 9.87 1110 4.48e-22 6.72e-21
## 14 Likability apref_e4 apref_e6 1111 1111 -3.34 1110 8.7 e- 4 1.3 e- 2
## 15 Likability apref_e5 apref_e6 1111 1111 -12.0 1110 2.46e-31 3.69e-30
## # ℹ 1 more variable: p.adj.signif <chr>
Multiple Regression
# Dummy coding
contrasts(data1$ASA.ending.fac) = contr.treatment(3)
contrasts(data1$CSA.ending.fac) = contr.treatment(3)
#Study2.reg<- lm(Likability.win ~ ASA.ending.fac + CSA.ending.fac + Age + gender.bin + SV + gender.bin*SV, data=data2)
#summary(Study2.reg)
#confint(Study2.reg)
library(lm.beta)
#Study2.stan<- lm.beta(Study2.reg)
#Study2.stan
# On 7/23/24 we decided to not control for ending because they are randomized across participants and to keep results consistent with the mixed effect models. We note this deviation from our preregistered data analysis plan in the manuscript.
Study2.reg.lim<- lm(Likability.win ~ Age + gender.bin + SV + gender.bin*SV, data=data2)
summary(Study2.reg.lim)
##
## Call:
## lm(formula = Likability.win ~ Age + gender.bin + SV + gender.bin *
## SV, data = data2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.41848 -0.35493 -0.02496 0.34342 1.44052
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.203156 0.046468 68.933 < 2e-16 ***
## Age -0.003288 0.001359 -2.419 0.01573 *
## gender.bin1 0.117536 0.039734 2.958 0.00316 **
## SV1 0.052516 0.079463 0.661 0.50882
## gender.bin1:SV1 -0.029880 0.089143 -0.335 0.73754
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5331 on 1106 degrees of freedom
## Multiple R-squared: 0.01856, Adjusted R-squared: 0.01501
## F-statistic: 5.228 on 4 and 1106 DF, p-value: 0.000357
round(confint(Study2.reg.lim),3)
## 2.5 % 97.5 %
## (Intercept) 3.112 3.294
## Age -0.006 -0.001
## gender.bin1 0.040 0.195
## SV1 -0.103 0.208
## gender.bin1:SV1 -0.205 0.145
lm.beta(Study2.reg.lim)
##
## Call:
## lm(formula = Likability.win ~ Age + gender.bin + SV + gender.bin *
## SV, data = data2)
##
## Standardized Coefficients::
## (Intercept) Age gender.bin1 SV1 gender.bin1:SV1
## NA -0.07254539 0.10473484 0.04666048 -0.02558475
library(car)
vif(Study2.reg.lim, type=c("predictor")) # variance inflation factors are all low
## GVIFs computed for predictors
## GVIF Df GVIF^(1/(2*Df)) Interacts With Other Predictors
## Age 1.013594 1 1.006774 -- gender.bin, SV
## gender.bin 1.013594 3 1.002253 SV Age
## SV 1.013594 3 1.002253 gender.bin Age
Assumption Checks
plot(Study2.reg.lim)
Exploratory Analyses
Mixed Model
S2.mixeddata<- data2 %>%
select(ID, Age, gender.bin, SV, Likability.win, Likability.nonSV.win)
mixed2_df<-NULL # 1.) construct empty (null) dataframe
for (row_i in 1:nrow(S2.mixeddata)) { # 2.) loop thru each row in S2.Mixeddata
for (vig in 0:1) { # 3.) loop twice , 0 and 1 (for each row of S2.Mixed data, we wil generate 2 rows in mixed2_df)
id = S2.mixeddata[[row_i,1]]
age = S2.mixeddata[[row_i,2]]
gen = S2.mixeddata[[row_i,3]]
sv = S2.mixeddata[[row_i,4]]
if (vig == 0) {
like = S2.mixeddata[[row_i,6]]
} else {
like = S2.mixeddata[[row_i,5]]
}
mixed2_temp_df <- data.frame(ID=id,
Age=age,
gender.bin=gen,
SV=sv,
vig=vig, # 0 = non-SV likability ratings; 1 = SV likability ratings
Likability=like) # create a temporary dataframe with the latest values, and append temp frame onto mixed2_df
rbind(mixed2_df,mixed2_temp_df)->mixed2_df
}
}
# Computing intraclass correlation coefficient (ICC)
mixed2_df$ID<- as.factor(mixed2_df$ID) #making ID a factor variable
(ICC.2 <- ICCbare(mixed2_df$ID, mixed2_df$Likability, data = mixed2_df)) # ICC = .40
## [1] 0.3996031
################## Visualizing Data ###################
x3<- ggplot(mixed2_df, aes(x = vig, y = Likability, colour=gender.bin, group=gender.bin)) + scale_x_continuous(breaks=c(0, 1), name="Vignette", labels=c("Non-SV Vignettes", "SV Vignettes")) + scale_color_hue(name= "Gender", labels=c("Cisgender Men", "Women/TNB+")) + ggtitle("Figure 4") +
stat_summary(fun="mean",geom="line")
x4<- x3 + ylim(1, 5)
x4
# + theme(text=element_text(size=12, family="serif")) ## code to change font to times new roman
a3<- ggplot(mixed2_df, aes(x = vig, y = Likability, colour=SV,group=SV)) + scale_x_continuous(breaks=c(0, 1), name="Vignette", labels=c("Non-SV Vignettes", "SV Vignettes"))+ scale_color_hue(name= "SV History", labels=c("Non-victim","Victim")) + ggtitle("Figure 5") +
stat_summary(fun="mean",geom="line")
a4<- a3 + ylim(1, 5)
a4
b2<- ggplot(mixed2_df, aes(x = vig, y = Likability)) + scale_x_continuous(breaks=c(0, 1), name="Vignette")+ ggtitle("Likability Scores Across Vignettes by SV History and Gender") +
stat_summary(fun="mean",geom="line")
b2 + facet_grid(gender.bin~SV, labeller=label_both) + ylim(1, 5)
####### Age ########
mixed2_df$Age.bin<- cut(mixed2_df$Age, breaks=c(0, 38.17, 100), labels=c("Mean Age", "+1 SD")) #creating age bins that are mean, 1 sd above mean, and 2 sds above mean
e2<- ggplot(mixed2_df, aes(x = vig, y = Likability, colour=Age.bin,group=Age.bin)) + scale_x_continuous(breaks=c(0, 1), name="Vignette", labels=c("Non-SV Vignettes", "SV Vignettes"))+ ggtitle("Figure 6") + labs(color = "Age") +
stat_summary(fun="mean",geom="line")
e3<- e2 + ylim(1, 5)
e3
############ Fixed Effect Structure ###############
#mod.5<- lmer(Likability ~ 1 + vig + (1 + vig | ID), data=mixed2_df, REML=TRUE) # Random intercept, random slope, model does not converge
mod.6<- lme4::lmer(Likability ~ 1 + vig + (1 | ID), data=mixed2_df, REML=TRUE)# Random intercept, fixed slope
#mod.7<- lmer(Likability ~ 1 + vig + (-1 | ID), data=mixed2_df, REML=TRUE) # Random slope, fixed intercept, does not converge
mod.8<- lm(Likability ~ 1 + vig, data=mixed2_df)
AIC(mod.6) # So model with the random intercept and fixed slope best fit the data
## [1] 2937.848
AIC(mod.8)
## [1] 3178.551
######################### Error Structure ######################
# Independent error structure
ind.error.2 <- nlme::lme(Likability ~ 1 + vig, # Fixed portion of the model
random = ~ 1 | ID, # Random portion of the model
data = mixed2_df)
nlme::ACF(ind.error.2) # If independent error structure is plausible, we expect autocorrelations to be about 0 (this does not seem to be the case)
## lag ACF
## 1 0 1.0000000
## 2 1 -0.4537313
# Compound Symmetry
comp.error.2 <- nlme::lme(Likability ~ 1 + vig, # Fixed portion of the model
random = ~ 1 | ID, # Random portion of the model
data = mixed2_df,
correlation=nlme::corCompSymm(.3),
control=lme4::lmerControl(check.conv.singular = lme4::.makeCC(action = "ignore", tol = 1e-3)))
# FIRST ORDER AUTOREGRESSIVE STRUCTURE
auto.reg.2 <- nlme::lme(Likability ~ 1 + vig, # Fixed portion of the model
random = ~ 1 | ID, # Random portion of the model
data = mixed2_df,
correlation=nlme::corAR1())
# Unstructured
unstructured.2 <- nlme::lme(Likability ~ 1 + vig,
random = ~ 1 | ID,
data = mixed2_df,
correlation=nlme::corSymm(),
weights = nlme::varIdent(form = ~ 1 | vig)) # unstructured error variances provide the most flexibility but contain more parameters
AIC(ind.error.2)
## [1] 2937.848
AIC(comp.error.2)
## [1] 2939.848
AIC(auto.reg.2)
## [1] 2939.848
AIC(unstructured.2) # So the model with the unstructured error variance provides the best fit to the data
## [1] 2895.242
########################### FULL MODEL ###################
full.mod.2 <- nlme::lme(Likability ~ vig*Age + vig*gender.bin + vig*SV,
random = ~ 1 | ID,
data = mixed2_df,
correlation=nlme::corSymm(),
weights = nlme::varIdent(form = ~ 1 | vig))
summary(full.mod.2)
## Linear mixed-effects model fit by REML
## Data: mixed2_df
## AIC BIC logLik
## 2897.726 2966.157 -1436.863
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.3134218 0.3144443
##
## Correlation Structure: General
## Formula: ~1 | ID
## Parameter estimate(s):
## Correlation:
## 1
## 2 0.097
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | vig
## Parameter estimates:
## 0 1
## 1.000000 1.370463
## Fixed effects: Likability ~ vig * Age + vig * gender.bin + vig * SV
## Value Std.Error DF t-value p-value
## (Intercept) 3.262916 0.03835940 1107 85.06169 0.0000
## vig -0.057696 0.04390247 1107 -1.31419 0.1891
## Age 0.003691 0.00112781 1107 3.27290 0.0011
## gender.bin1 0.088156 0.02961433 1107 2.97679 0.0030
## SV1 -0.057482 0.02972054 1107 -1.93409 0.0534
## vig:Age -0.006940 0.00129078 1107 -5.37633 0.0000
## vig:gender.bin1 0.023437 0.03389371 1107 0.69149 0.4894
## vig:SV1 0.086200 0.03401526 1107 2.53416 0.0114
## Correlation:
## (Intr) vig Age gndr.1 SV1 vig:Ag vg:g.1
## vig -0.380
## Age -0.806 0.306
## gender.bin1 -0.435 0.165 0.036
## SV1 -0.144 0.055 0.052 -0.336
## vig:Age 0.306 -0.806 -0.380 -0.014 -0.020
## vig:gender.bin1 0.165 -0.435 -0.014 -0.380 0.128 0.036
## vig:SV1 0.055 -0.144 -0.020 0.128 -0.380 0.052 -0.336
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.07492635 -0.51601748 -0.01151473 0.51399456 2.83986410
##
## Number of Observations: 2222
## Number of Groups: 1111
0.3278391^2 #Variance around the intercept
## [1] 0.1074785
AIC(full.mod.2)
## [1] 2897.726
r2glmm::r2beta(full.mod.2) # R2 square
## Effect Rsq upper.CL lower.CL
## 1 Model 0.063 0.086 0.047
## 6 vig:Age 0.008 0.017 0.002
## 3 Age 0.005 0.012 0.001
## 4 gender.bin1 0.004 0.011 0.000
## 8 vig:SV1 0.002 0.007 0.000
## 5 SV1 0.002 0.007 0.000
## 2 vig 0.000 0.004 0.000
## 7 vig:gender.bin1 0.000 0.003 0.000
nlme::intervals(full.mod.2, which="fixed") # Confidence intervals
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 3.187650537 3.262915878 3.3381812188
## vig -0.143837634 -0.057696183 0.0284452669
## Age 0.001478317 0.003691197 0.0059040770
## gender.bin1 0.030049055 0.088155607 0.1462621580
## SV1 -0.115797236 -0.057482292 0.0008326514
## vig:Age -0.009472301 -0.006939652 -0.0044070029
## vig:gender.bin1 -0.043066043 0.023437111 0.0899402659
## vig:SV1 0.019458374 0.086200035 0.1529416949
emvigSV.2 <- emmeans(full.mod.2, ~ vig*SV) #calculate predicted means for SV by vig
emvigSV.2
## vig SV emmean SE df lower.CL upper.CL
## 0 -1 3.40 0.0166 1107 3.37 3.44
## 1 -1 3.18 0.0199 1107 3.14 3.21
## 0 1 3.35 0.0250 1107 3.30 3.40
## 1 1 3.20 0.0300 1107 3.15 3.26
##
## Results are averaged over the levels of: gender.bin
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
contrast(emvigSV.2, "revpairwise",by="SV") # simple slopes
## SV = -1:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.229 0.0189 1107 -12.066 <.0001
##
## SV = 1:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.142 0.0286 1107 -4.978 <.0001
##
## Results are averaged over the levels of: gender.bin
## Degrees-of-freedom method: containment
emvigGen.2<- emmeans(full.mod.2, ~vig*gender.bin) #calculate predicted means for gender by vig
contrast(emvigGen.2, "revpairwise",by="gender.bin")
## gender.bin = -1:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.197 0.0285 1107 -6.923 <.0001
##
## gender.bin = 1:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.174 0.0190 1107 -9.140 <.0001
##
## Results are averaged over the levels of: SV
## Degrees-of-freedom method: containment
emtrends(full.mod.2, ~ vig, var="Age") #obtain simple slopes of Age by vignette
## vig Age.trend SE df lower.CL upper.CL
## 0 0.00369 0.00113 1107 0.00148 0.005904
## 1 -0.00325 0.00135 1107 -0.00590 -0.000593
##
## Results are averaged over the levels of: gender.bin, SV
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
mylist.2<- list(Age=c(26.32, 38.17), vig=c(1, 0)) # Plotting at mean and one SD above
emcontcat.2<- emmeans(full.mod.2, ~vig*Age, at=mylist.2)
contrast(emcontcat.2, "pairwise", by="Age")
## Age = 26.3:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.186 0.0173 1107 -10.719 <.0001
##
## Age = 38.2:
## contrast estimate SE df t.ratio p.value
## vig1 - vig0 -0.268 0.0232 1107 -11.565 <.0001
##
## Results are averaged over the levels of: gender.bin, SV
## Degrees-of-freedom method: containment
