install.packages(c("WebPower", "esc", "appRiori", "shiny", "ggpubr"),
repos = list(CRAN="http://cran.rstudio.com/"))
## Installing packages into '/Users/perso204/Library/R/x86_64/4.5/library'
## (as 'lib' is unspecified)
##
## The downloaded binary packages are in
## /var/folders/96/xz1fg2r56mb2kymyn8r37rqc0000gq/T//RtmpuHcIMm/downloaded_packages
library(WebPower)
## Loading required package: MASS
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: lavaan
## This is lavaan 0.6-20
## lavaan is FREE software! Please report any bugs.
## Loading required package: parallel
## Loading required package: PearsonDS
library(esc)
########################## Hypothesis 1 ########################
## Person et al. 2025 Study 1
# SV M = 3.17; SD = .59
# Car accident M = 3.33; SD = .65
# We are taking means from the car accident vignette to
#most closely match our current study design, but this effect size was
#very similar to the non-SV vignette average what was published in Person et al. (2025)
esc_mean_sd(grp2m = 3.17, # mean of sexual violence
grp2sd = .59, # standard deviation of sexual violence vignettes
grp2n = 760, # sample size for grp 2
grp1m = 3.33, # mean of car accident vignettes
grp1sd = .65, # standard deviation of car accident vignettes
grp1n = 760, # sample size for grp 1
es.type = "d")
##
## Effect Size Calculation for Meta Analysis
##
## Conversion: mean and sd to effect size d
## Effect Size: 0.2578
## Standard Error: 0.0515
## Variance: 0.0027
## Lower CI: 0.1568
## Upper CI: 0.3587
## Weight: 376.8700
## Person et al. (2025) Study 2
# SV M = 3,20; SD = .56
# Car accident M = 3.42; SD = ..69
esc_mean_sd(grp2m = 3.20, # mean of sexual violence
grp2sd = .56, # standard deviation of sexual violence vignettes
grp2n = 1112, # sample size for grp 2
grp1m = 3.42, # mean of car accident vignettes
grp1sd = .69, # standard deviation of car accident vignettes
grp1n = 1112, # sample size for grp 1
es.type = "d")
##
## Effect Size Calculation for Meta Analysis
##
## Conversion: mean and sd to effect size d
## Effect Size: 0.3501
## Standard Error: 0.0427
## Variance: 0.0018
## Lower CI: 0.2664
## Upper CI: 0.4339
## Weight: 547.6094
## Average effect size across two studies = .31
(.35 + .26)/2
## [1] 0.305
### cohens d of .31 = f of .155
(effectsize<- 0.155*sqrt(3/(1-0.40))) #effect size with correlation
## [1] 0.3465905
#between repeated measures,
#effect.size.RP = effect.size.G*sqrt(#(repeatedMeasure)/(1-corr)), 0.55*sqrt(3/(1-0.40)
#WebPower
wp.rmanova(n = NULL,
ng = 1, #number of groups
nm = 3, #number of measurements
f = effectsize,
alpha = 0.05,
power = .80,
nscor = 1, #default, assumes that assumption of sphericity is met
type = 1) # type=1 specifies that we want power for the repeated effect,
## Repeated-measures ANOVA analysis
##
## n f ng nm nscor alpha power
## 81.71282 0.3465905 1 3 1 0.05 0.8
##
## NOTE: Power analysis for within-effect test
## URL: http://psychstat.org/rmanova
#whereas 0 (default) is for between effect
#################### Hypothesis 2#############################
# Comparison between no-trauma and less-stigmatized trauma, Harter et al. (2009)
# Means and SDs came from men (because current study is
#only recruiting men) on the likability outcome variable
esc_mean_sd(grp2m = 3.08, # mean of sexual violence
grp2sd = .62, # standard deviation of sexual violence vignettes
grp2n = 31, # sample size for grp 2
grp1m = 3.24, # mean of car accident vignettes
grp1sd = .86, # standard deviation of car accident vignettes
grp1n = 27, # sample size for grp 1
es.type = "d")
##
## Effect Size Calculation for Meta Analysis
##
## Conversion: mean and sd to effect size d
## Effect Size: 0.2159
## Standard Error: 0.2640
## Variance: 0.0697
## Lower CI: -0.3016
## Upper CI: 0.7333
## Weight: 14.3478
# d=.22 is f = .11
(effectsize.2<- 0.11*sqrt(3/(1-0.40))) #effect size with correlation
## [1] 0.2459675
#between repeated measures,
#effect.size.RP = effect.size.G*sqrt(#(repeatedMeasure)/(1-corr)), 0.55*sqrt(3/(1-0.40)
#WebPower
wp.rmanova(n = NULL,
ng = 1, #number of groups
nm = 3, #number of measurements
f = effectsize.2,
alpha = 0.05,
power = .80,
nscor = 1, #default, assumes that assumption of sphericity is met
type = 1) # type=1 specifies that we want power for the repeated effect,
## Repeated-measures ANOVA analysis
##
## n f ng nm nscor alpha power
## 160.7537 0.2459675 1 3 1 0.05 0.8
##
## NOTE: Power analysis for within-effect test
## URL: http://psychstat.org/rmanova
#whereas 0 (default) is for between effect
library(readr)
SASD_data <- read_csv("~/Downloads/SASD_data.csv")
## Rows: 192 Columns: 125
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): Q_RelevantIDLastStartDate, Year_6_TEXT, SexualAttraction_7_TEXT, ...
## dbl (114): Status, Progress, Duration (in seconds), Finished, Q_RecaptchaSco...
## lgl (3): Q_RelevantIDDuplicate, Q_DuplicateRespondent, Gender_6_TEXT
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# In line with our eligibility criteria, we will exclude participants who do not identify as men and indicate that they are only or mostly attracted to men prior to data analysis. If participants respond to at least three of the five dating evaluation items, we will use their mean score in analyses. Outliers on the dependent variable will be recoded to three standard deviations from the mean. Participant responses that score >75 on the Qualtrics Relevant ID Duplicate variable will be reviewed as scores above this threshold typically suggest a duplicate response (Fraud detection, 2023). Participant responses will be removed prior to data analysis if they complete the survey in under two minutes and score > 90 on the Qualtrics Relevant ID Duplicate variable.
library(dplyr)
##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:MASS':
##
## select
##
## The following objects are masked from 'package:stats':
##
## filter, lag
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
SASD_data.1<- SASD_data %>% # filtering out 1 person who didn't consent
filter(consent == 1)
SASD_data.2<- SASD_data.1 %>%
filter(debrief_4 == "why did you need me to put my x500?" | is.na(debrief_4)) # removed 4 people who asked to have their data deleted after debriefing
SASD_data.3<- SASD_data.2 %>%
filter(`Duration (in seconds)`> 120) # removing 15 people who completed survey in under 2 minutes
SASD_data.4<- SASD_data.3 %>%
filter(Q_RelevantIDDuplicateScore<75) # removing 6 duplicate responders
SASD_data.5<- SASD_data.4 %>%
filter(Gender==7) # removing 49 participants that did not identify as men (which was part of eligibility requirement)
SASD_data.6<- SASD_data.5 %>%
filter(SexualAttraction!= 4 & SexualAttraction!=5) # removing 4 people who said they were only or mostly attracted to women
SASD_data.7<- SASD_data.6 %>%
filter(Q_RecaptchaScore>.5) # forgot to preregister this but removed 2 people who were likely bots
SASD<- SASD_data.7
#Demographics
#Variable coding is Freshman 1, Sophomore 2, Junior 3, Senior 4, Graduate student 5, Other 6
table(SASD$Year)
##
## 1 2 3 4 6
## 49 25 21 13 3
prop.table(table(SASD$Year))
##
## 1 2 3 4 6
## 0.44144144 0.22522523 0.18918919 0.11711712 0.02702703
# Someone said they were born in 2005, so recoding that person to age instead
SASD$Age_1[SASD$Age_1 == 2005] <- 20
mean(SASD$Age_1, na.rm=T)
## [1] 20.06306
sd(SASD$Age_1, na.rm=T)
## [1] 4.454373
table(SASD$Sex) #So 3 transgender participants
##
## 1 2
## 108 3
prop.table(table(SASD$Sex))
##
## 1 2
## 0.97297297 0.02702703
# Variable coding is 1 Straight or heterosexual, 2 Gay, 4 Lesbian, 6 Bisexual, 3 Queer, 7 Pansexual, 8 Asexual, 9 I am not sure, 10 I #don't know what this question means, 11 Decline to answer
table(SASD$SexualOrientation)
##
## 1 3 6 7 11
## 101 1 6 1 2
prop.table(table(SASD$SexualOrientation))
##
## 1 3 6 7 11
## 0.909909910 0.009009009 0.054054054 0.009009009 0.018018018
# 1=only attracted to women, 2=mostly attracted to women, 3 =equally attracted to men and women, 6 =not sure, 7= other
table(SASD$SexualAttraction)
##
## 1 2 3 6 7
## 88 13 6 3 1
prop.table(table(SASD$SexualAttraction))
##
## 1 2 3 6 7
## 0.792792793 0.117117117 0.054054054 0.027027027 0.009009009
# 1=single and looking, 2= single and not looking, 3=short term, 4=committed relationship, 6=married
table(SASD$Rel_status)
##
## 1 2 3 4 6
## 36 23 8 42 2
prop.table(table(SASD$Rel_status))
##
## 1 2 3 4 6
## 0.32432432 0.20720721 0.07207207 0.37837838 0.01801802
# 1=white, 2=Black, 3=Asian, 4=nativehawaian, 5=nativeamerican, 6=multiracial, 7=other, 8=hispanic
table(SASD$Race)
##
## 1 2 3 5 6 7 8
## 72 8 19 1 5 3 2
prop.table(table(SASD$Race))
##
## 1 2 3 5 6 7
## 0.654545455 0.072727273 0.172727273 0.009090909 0.045454545 0.027272727
## 8
## 0.018181818
# 1= strong democrat to 7 = strong republican
table(SASD$PoliticalIdentitiy)
##
## 1 2 3 4 5 6 7
## 23 20 14 22 13 12 6
prop.table(SASD$PoliticalIdentitiy)
## [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [51] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [76] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [101] NA NA NA NA NA NA NA NA NA NA NA
SASD$ID<- 1:111 #giving everyone a participant ID
library(dplyr)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpack
library(stringr)
data_long_items <- SASD %>%
pivot_longer(
# all the profile-condition-item columns
cols = matches("^(Ashley|Madison|Emma)(VT|N|C)_\\d+$"),
names_to = c("profile", "condition", "item"),
names_pattern = "^(Ashley|Madison|Emma)(VT|N|C)_(\\d+)$",
values_to = "response"
)
data <- data_long_items %>%
mutate(item = paste0("Q", item)) %>%
pivot_wider(
id_cols = c(ID, profile, condition),
names_from = item,
values_from = response
) %>%
# keep only blocks the person actually saw
filter(if_any(starts_with("Q"), ~ !is.na(.)))
#creating mean score for RAS measure
data <- data %>%
rowwise() %>%
mutate(RAS = mean(c_across(Q1:Q5), na.rm = TRUE))%>%
ungroup()
# pulling IDs and conditions from original dataset to make sure everything looks correct
conds<- select(SASD, ID, ConditionN, ConditionC, ConditionVT)
#everything looks correct
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ purrr 1.0.4
## ✔ ggplot2 4.0.1 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ dplyr::select() masks MASS::select()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggpubr)
library(rstatix)
##
## Attaching package: 'rstatix'
##
## The following objects are masked from 'package:esc':
##
## cohens_d, eta_squared
##
## The following object is masked from 'package:MASS':
##
## select
##
## The following object is masked from 'package:stats':
##
## filter
# making sure conditions were successfully randomized
rand<- table(data$condition, data$profile)
rand
##
## Ashley Emma Madison
## C 39 31 40
## N 34 41 36
## VT 37 38 35
chi_square_test <- chisq.test(rand)
chi_square_test # p>.05so conditions were successfully randomized
##
## Pearson's Chi-squared test
##
## data: rand
## X-squared = 2.1541, df = 4, p-value = 0.7074
# means and sds by condition
data %>%
group_by(condition) %>%
get_summary_stats(RAS, type = "mean_sd")
# means and sds by profile, these were randomized across conditions
data %>%
group_by(profile) %>%
get_summary_stats(RAS, type = "mean_sd")
boxplot(RAS ~ profile, data = data,
main = "Profile Photo Ratings",
xlab = "Profile Photo",
ylab = "Dating Evaluations")
boxplot(RAS ~ condition, data = data,
main = "Condition Ratings",
xlab = "Condition",
ylab = "Dating Evaluations")
hist(data$RAS, main="RAS", xlab="RAS", col='hotpink', sub=paste("Skewness:",
round(e1071::skewness(data$RAS, na.rm=TRUE), 2)))
qqnorm(data$RAS, pch = 1, frame = FALSE, main="QQ Plot of RAS")
qqline(data$RAS, col = "hotpink", lwd = 2)
data %>%
group_by(condition) %>%
shapiro_test(RAS)
# none of the groups are normally distributed when using Shapiro wilk test
ggqqplot(data, "RAS", facet.by = "condition")
sum(is.na(data$Q1))
## [1] 0
sum(is.na(data$Q2))
## [1] 0
sum(is.na(data$Q3))
## [1] 0
sum(is.na(data$Q4))
## [1] 0
sum(is.na(data$Q5))
## [1] 0
# no missing data
# we preregistered that we'd recode outliers that were abs(3) mean, but those values aren't possible.
mean<- mean(data$RAS, na.rm=T)
sd<- sd(data$RAS, na.rm=T)
sd3<-3*sd
mean-sd3
## [1] 0.1706184
mean+sd3
## [1] 8.098868
# no extreme outliers
data %>%
group_by(condition) %>%
identify_outliers(RAS)
ras.aov <- anova_test(data = data, dv = RAS, wid = ID, within = condition)
get_anova_table(ras.aov)
# not significant. this makes sense based on the visual plots as well