install.packages(c('tidyverse', 'MASS',
               'magrittr', 'robustHD',
               'tibble', 'psych',
               'knitr', 'tidyr', 'ggplot2',
               'lavaan', 'semPlot', 'car',
               'lmtest', 'ggpubr',
               'FSA', 'rstatix',
               'rcompanion', 'coin', 'lm.beta', 'Hmisc', 'nnet', 
               'gtsummary',
               'ggeffects', 'effects', 'ggeffects', 'foreign', 
               'nnet', 'reshape2', 'naniar', 'dplyr', 
               'expss', 'writexl', 'datawizard'), repos = list(CRAN="http://cran.rstudio.com/"))
## Installing packages into '/Users/perso204/Library/R/x86_64/4.4/library'
## (as 'lib' is unspecified)
## 
## The downloaded binary packages are in
##  /var/folders/96/xz1fg2r56mb2kymyn8r37rqc0000gq/T//Rtmpc9tCcr/downloaded_packages

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
contrasts(data1$ASA.ending.fac) = contr.sum(3)

################# 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 drop variable use NULL: let(mtcars, am = NULL) %>% head()
## 
## 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
## 
## Use 'expss_output_rnotebook()' to display tables inside R Notebooks.
##  To return to the console output, use 'expss_output_default()'.
## 
## 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),
    IQR = IQR(Age, na.rm = TRUE)
  )
## # A tibble: 3 × 6
##   Sample    count  mean    sd median   IQR
##   <chr>     <int> <dbl> <dbl>  <dbl> <dbl>
## 1 UMN         290  20.1  2.52     20   2  
## 2 WWU         280  20.1  4.42     19   2  
## 3 qualtrics   190  50.9 15.8      53  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
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

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

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
# 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, SV.num, gender.num, Age, DiffScores)

(S1.corr<- rcorr(as.matrix(Study1.cor)))
##                Likability.win SV.num gender.num   Age DiffScores
## Likability.win           1.00   0.09       0.06 -0.09      -0.48
## SV.num                   0.09   1.00       0.30  0.00      -0.15
## gender.num               0.06   0.30       1.00  0.02      -0.10
## Age                     -0.09   0.00       0.02  1.00       0.08
## DiffScores              -0.48  -0.15      -0.10  0.08       1.00
## 
## n
##                Likability.win SV.num gender.num Age DiffScores
## Likability.win            760    760        760 759        760
## SV.num                    760    760        760 759        760
## gender.num                760    760        760 759        760
## Age                       759    759        759 759        759
## DiffScores                760    760        760 759        760
## 
## P
##                Likability.win SV.num gender.num Age    DiffScores
## Likability.win                0.0137 0.1282     0.0096 0.0000    
## SV.num         0.0137                0.0000     0.9744 0.0000    
## gender.num     0.1282         0.0000            0.5591 0.0063    
## Age            0.0096         0.9744 0.5591            0.0369    
## DiffScores     0.0000         0.0000 0.0063     0.0369
# 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

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)
## 
## Call:
## lm(formula = Likability.n1 ~ ASA.ending.fac + Age + gender.bin + 
##     SV + gender.bin * SV, data = data1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.74327 -0.32023 -0.00399  0.32065  2.32953 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.872662   0.060008  47.871  < 2e-16 ***
## ASA.ending.fac2  0.507450   0.047266  10.736  < 2e-16 ***
## ASA.ending.fac3  0.554533   0.046853  11.835  < 2e-16 ***
## Age             -0.003815   0.001214  -3.143  0.00174 ** 
## gender.bin1     -0.005245   0.053590  -0.098  0.92206    
## SV1              0.024134   0.108894   0.222  0.82466    
## gender.bin1:SV1  0.105746   0.117470   0.900  0.36830    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5275 on 752 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2001, Adjusted R-squared:  0.1938 
## F-statistic: 31.36 on 6 and 752 DF,  p-value: < 2.2e-16
confint(Study1.reg)
##                        2.5 %       97.5 %
## (Intercept)      2.754857874  2.990465132
## ASA.ending.fac2  0.414661718  0.600239097
## ASA.ending.fac3  0.462553405  0.646511620
## Age             -0.006197907 -0.001431822
## gender.bin1     -0.110448730  0.099959428
## SV1             -0.189637386  0.237906176
## gender.bin1:SV1 -0.124861202  0.336353491
library(lm.beta)
Study1.stan<- lm.beta(Study1.reg)
Study1.stan
## 
## Call:
## lm(formula = Likability.n1 ~ ASA.ending.fac + Age + gender.bin + 
##     SV + gender.bin * SV, data = data1)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac2 ASA.ending.fac3             Age     gender.bin1 
##              NA     0.406680432     0.447452238    -0.102603940    -0.003762989 
##             SV1 gender.bin1:SV1 
##     0.020335015     0.087860131
vif(Study1.reg, type=c("predictor"))
## GVIFs computed for predictors
##                    GVIF Df GVIF^(1/(2*Df)) Interacts With
## ASA.ending.fac 1.007915  2        1.001973           --  
## Age            1.002168  1        1.001083           --  
## gender.bin     1.009578  3        1.001590             SV
## SV             1.009578  3        1.001590     gender.bin
##                              Other Predictors
## ASA.ending.fac            Age, gender.bin, SV
## Age            ASA.ending.fac, gender.bin, SV
## gender.bin                ASA.ending.fac, Age
## SV                        ASA.ending.fac, Age

Assumption Checks

plot(Study1.reg)

# The residual vs fitted plot shows if residuals have non-linear pattern. There appears to be a linear pattern as the line is doesn't follow other patter (e.g., U curve) though there is slight gap in the residuals in middle. 
# The QQ plot shows if the residuals are normally distributed, they look very normally distributed
# The Scale-location plot shows if the residuals are equally spread along a range of predictors (e.g., homoscedasticity). You should see horizontal line with randomly spread points about it. The line is horizontal, but again gap of data in middle of plot
# Residuals vs leverage plot is looking at outliers, no issues there

Exploratory Analyses

SV

SV.1.reg<- lm(Likability.win ~ ASA.ending.fac + Age + SV, data=data1)
summary(SV.1.reg)
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + Age + SV, data = data1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.74113 -0.31697 -0.00683  0.32816  2.30627 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.871092   0.049875  57.566  < 2e-16 ***
## ASA.ending.fac2  0.503064   0.046874  10.732  < 2e-16 ***
## ASA.ending.fac3  0.552328   0.046513  11.875  < 2e-16 ***
## Age             -0.003740   0.001205  -3.104  0.00198 ** 
## SV1              0.117621   0.038498   3.055  0.00233 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.524 on 754 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1986, Adjusted R-squared:  0.1943 
## F-statistic:  46.7 on 4 and 754 DF,  p-value: < 2.2e-16
confint(SV.1.reg)
##                        2.5 %       97.5 %
## (Intercept)      2.773182012  2.969001202
## ASA.ending.fac2  0.411046047  0.595082924
## ASA.ending.fac3  0.461017867  0.643638548
## Age             -0.006105728 -0.001375086
## SV1              0.042045010  0.193196952
SV.1.reg.stan<- lm.beta(SV.1.reg)
SV.1.reg.stan
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + Age + SV, data = data1)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac2 ASA.ending.fac3             Age             SV1 
##              NA      0.40566134      0.44843260     -0.10122413      0.09971791

Gender

Gender.1.reg<- lm(Likability.win ~ ASA.ending.fac + Age + gender.bin, data=data1)
summary(Gender.1.reg)
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + Age + gender.bin, 
##     data = data1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.77593 -0.31898 -0.00612  0.33452  2.29729 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.881439   0.057644  49.987  < 2e-16 ***
## ASA.ending.fac2  0.495113   0.047085  10.515  < 2e-16 ***
## ASA.ending.fac3  0.549154   0.046775  11.740  < 2e-16 ***
## Age             -0.003766   0.001211  -3.109  0.00195 ** 
## gender.bin1      0.057559   0.045430   1.267  0.20555    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5267 on 754 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1904, Adjusted R-squared:  0.1861 
## F-statistic: 44.32 on 4 and 754 DF,  p-value: < 2.2e-16
confint(Gender.1.reg)
##                       2.5 %       97.5 %
## (Intercept)      2.76827805  2.994600917
## ASA.ending.fac2  0.40267971  0.587546536
## ASA.ending.fac3  0.45732905  0.640978895
## Age             -0.00614409 -0.001388306
## gender.bin1     -0.03162440  0.146742378
Gender.1.reg.stan<- lm.beta(Gender.1.reg)
Gender.1.reg.stan
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + Age + gender.bin, 
##     data = data1)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac2 ASA.ending.fac3             Age     gender.bin1 
##              NA      0.39924951      0.44585546     -0.10192209      0.04155371

Sample

Conducting analyses (minus interaction term) on the separate samples.

UMN.1 <- data1 %>% 
dplyr::filter(Sample== "UMN")

WWU.1 <- data1 %>% 
dplyr::filter(Sample== "WWU")

qualtrics.1 <- data1 %>% 
dplyr::filter(Sample== "qualtrics")


UMN.1.reg<- lm(Likability.win ~ ASA.ending.fac + gender.bin + SV, data=UMN.1)
summary(UMN.1.reg)
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + gender.bin + SV, 
##     data = UMN.1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6490 -0.2720  0.0406  0.3086  1.3535 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.70081    0.06303  42.851  < 2e-16 ***
## ASA.ending.fac2  0.57788    0.06771   8.534 8.50e-16 ***
## ASA.ending.fac3  0.54821    0.06740   8.134 1.29e-14 ***
## gender.bin1      0.11039    0.06424   1.718  0.08680 .  
## SV1              0.18712    0.06197   3.019  0.00276 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4712 on 285 degrees of freedom
## Multiple R-squared:  0.2801, Adjusted R-squared:  0.2699 
## F-statistic: 27.72 on 4 and 285 DF,  p-value: < 2.2e-16
confint(UMN.1.reg)
##                       2.5 %    97.5 %
## (Intercept)      2.57674717 2.8248663
## ASA.ending.fac2  0.44460007 0.7111557
## ASA.ending.fac3  0.41555327 0.6808708
## gender.bin1     -0.01605175 0.2368225
## SV1              0.06514278 0.3091068
UMN.1.reg.stan<- lm.beta(UMN.1.reg)
UMN.1.reg.stan
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + gender.bin + SV, 
##     data = UMN.1)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac2 ASA.ending.fac3     gender.bin1             SV1 
##              NA      0.49126522      0.46856843      0.09275466      0.16300183
WWU.1.reg<- lm(Likability.win ~ ASA.ending.fac + gender.bin + SV, data=WWU.1)
summary(WWU.1.reg)
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + gender.bin + SV, 
##     data = WWU.1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.63625 -0.29884 -0.03625  0.28138  1.20116 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.66116    0.07572  35.145  < 2e-16 ***
## ASA.ending.fac2  0.51359    0.06686   7.682 2.77e-13 ***
## ASA.ending.fac3  0.55100    0.06535   8.432 1.94e-15 ***
## gender.bin1      0.15961    0.07439   2.145   0.0328 *  
## SV1              0.06448    0.05520   1.168   0.2438    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4448 on 275 degrees of freedom
## Multiple R-squared:  0.2597, Adjusted R-squared:  0.249 
## F-statistic: 24.12 on 4 and 275 DF,  p-value: < 2.2e-16
confint(WWU.1.reg)
##                       2.5 %    97.5 %
## (Intercept)      2.51209832 2.8102285
## ASA.ending.fac2  0.38197699 0.6452108
## ASA.ending.fac3  0.42235229 0.6796445
## gender.bin1      0.01315629 0.3060543
## SV1             -0.04419849 0.1731558
WWU.1.reg.stan<- lm.beta(WWU.1.reg)
WWU.1.reg.stan
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + gender.bin + SV, 
##     data = WWU.1)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac2 ASA.ending.fac3     gender.bin1             SV1 
##              NA      0.47212870      0.51530559      0.11542834      0.06292323
qualtrics.1.reg<-  lm(Likability.n1 ~ ASA.ending.fac + gender.bin + SV, data=qualtrics.1)
summary(qualtrics.1.reg)
## 
## Call:
## lm(formula = Likability.n1 ~ ASA.ending.fac + gender.bin + SV, 
##     data = qualtrics.1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.66370 -0.41949 -0.03013  0.45942  2.01792 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.98208    0.13007  22.926  < 2e-16 ***
## ASA.ending.fac2  0.38180    0.12003   3.181  0.00172 ** 
## ASA.ending.fac3  0.55579    0.12106   4.591 8.13e-06 ***
## gender.bin1     -0.31838    0.12438  -2.560  0.01128 *  
## SV1              0.07688    0.10159   0.757  0.45018    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6766 on 185 degrees of freedom
## Multiple R-squared:  0.1319, Adjusted R-squared:  0.1131 
## F-statistic: 7.025 on 4 and 185 DF,  p-value: 2.752e-05
confint(qualtrics.1.reg)
##                      2.5 %      97.5 %
## (Intercept)      2.7254634  3.23870366
## ASA.ending.fac2  0.1449988  0.61859725
## ASA.ending.fac3  0.3169568  0.79462262
## gender.bin1     -0.5637754 -0.07298855
## SV1             -0.1235528  0.27730817
qualtrics.1.reg.stan<- lm.beta(qualtrics.1.reg)
qualtrics.1.reg.stan
## 
## Call:
## lm(formula = Likability.n1 ~ ASA.ending.fac + gender.bin + SV, 
##     data = qualtrics.1)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac2 ASA.ending.fac3     gender.bin1             SV1 
##              NA      0.25276945      0.36365582     -0.17945207      0.05313716
# Looking into qualtrics sample more closely as gender differences were odd
qualtrics.S1<- subset(data1, data1$Sample=="qualtrics")

cor.test(qualtrics.S1$Likability, qualtrics.S1$Age)
## 
##  Pearson's product-moment correlation
## 
## data:  qualtrics.S1$Likability and qualtrics.S1$Age
## t = -0.15261, df = 188, p-value = 0.8789
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1532400  0.1314321
## sample estimates:
##         cor 
## -0.01112942
cor.test(UMN.1$Likability, UMN.1$Age)
## 
##  Pearson's product-moment correlation
## 
## data:  UMN.1$Likability and UMN.1$Age
## t = -2.1285, df = 288, p-value = 0.03415
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.236239194 -0.009400932
## sample estimates:
##        cor 
## -0.1244458
cor.test(WWU.1$Likability, WWU.1$Age)
## 
##  Pearson's product-moment correlation
## 
## data:  WWU.1$Likability and WWU.1$Age
## t = -1.0845, df = 277, p-value = 0.2791
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.18107322  0.05281069
## sample estimates:
##         cor 
## -0.06502423
table(qualtrics.S1$Gender, qualtrics.S1$SV)
##                             
##                              -1  1
##   Male                       30  9
##   Female                     75 73
##   Genderqueer/non-conforming  0  0
##   Transgender                 2  0
##   Other (please specify)      1  0
group_by(qualtrics.1, Gender) %>%
  summarise(
    count = n(),
    mean = mean(Likability.win, na.rm = TRUE),
    sd = sd(Likability.win, na.rm = TRUE),
    median = median(Likability.win, na.rm = TRUE),
    IQR = IQR(Likability.win, na.rm = TRUE)
  ) # so mean of men slightly higher at bivariate level
## # A tibble: 4 × 6
##   Gender                 count  mean     sd median   IQR
##   <fct>                  <int> <dbl>  <dbl>  <dbl> <dbl>
## 1 Male                      39  3.30  0.792   3.1   1.05
## 2 Female                   148  3.02  0.680   3     0.9 
## 3 Transgender                2  2.75  0.354   2.75  0.25
## 4 Other (please specify)     1  3.7  NA       3.7   0
ggplot(qualtrics.1, aes(gender.bin, Likability.win)) + 
  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())

### Difference Scores

Linear Regression

S1.reg.diffScores<- lm(Diffscore.win ~ Age + gender.bin + SV + gender.bin*SV, data=data1)
summary(S1.reg.diffScores)
## 
## Call:
## lm(formula = Diffscore.win ~ Age + gender.bin + SV + gender.bin * 
##     SV, data = data1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9684 -0.2282 -0.0282  0.2166  0.9572 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.196371   0.036574   5.369 1.05e-07 ***
## Age              0.001805   0.000818   2.207   0.0276 *  
## gender.bin1     -0.052470   0.036058  -1.455   0.1460    
## SV1             -0.077163   0.073340  -1.052   0.2931    
## gender.bin1:SV1 -0.027682   0.079041  -0.350   0.7263    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3555 on 754 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.03558,    Adjusted R-squared:  0.03046 
## F-statistic: 6.953 on 4 and 754 DF,  p-value: 1.69e-05
confint(S1.reg.diffScores)
##                         2.5 %      97.5 %
## (Intercept)      0.1245714738 0.268169471
## Age              0.0001991282 0.003410839
## gender.bin1     -0.1232556149 0.018315704
## SV1             -0.2211366402 0.066810956
## gender.bin1:SV1 -0.1828482160 0.127483684
Study1.diff.stan<- lm.beta(S1.reg.diffScores)
Study1.diff.stan
## 
## Call:
## lm(formula = Diffscore.win ~ Age + gender.bin + SV + gender.bin * 
##     SV, data = data1)
## 
## Standardized Coefficients::
##     (Intercept)             Age     gender.bin1             SV1 gender.bin1:SV1 
##              NA      0.07899133     -0.06125596     -0.10578828     -0.03742400

Mixed model

S1.mixeddata<- data1 %>%
  select(ID, Age, gender.bin, SV, Likability, Likability.nonSV)


#for (row_i in 1:nrow(S1.mixeddata)) {     # loop thru each rwo in S1.Mixeddata
 # for (vig in 0:1) {                      # loop twice , 0 and 1
    #print(vig)
    #print(row_i)
    #id  = S1.mixeddata[row_i,1]
    #print(id)
    #age = S1.mixeddata[row_i,2]
    #print(age)
    #gen = S1.mixeddata[row_i,3]
    #print(gen)
    #sv  = S1.mixeddata[row_i,4]
    #print(sv)
   # if (vig == 0) {
        #like = S1.mixeddata[row_i,6]
   # } else {
       # like = S1.mixeddata[row_i,5]
   # }
   # print(like)
    # new_df[nrow(new_df) + 1,] = c(id, age, gen, sv, vig, like)
  #}

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 frequences
## 
##                       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)

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.199011
sd(data2$Likability, na.rm=T)
## [1] 0.5581802
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  723 223
   1  18 148
   #Total cases  741 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      723  3.17 0.563    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           395  3.11 0.580    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      392  3.24 0.524    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     374  3.31 0.558    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     376  3.33 0.533    3.3 0.7
mean(data2$Likability)
## [1] 3.199011
median(data2$Likability)
## [1] 3.2
sd(data2$Likability)
## [1] 0.5581802
################# SV by Gender #############

cross_cases(data2, SV, gender.bin)
 gender.bin 
 -1   1 
 SV 
   -1  343 380
   1  52 337
   #Total cases  395 717 
cross_cases(data2, lec8_bin, gender.bin) # sexual assault
 gender.bin 
 -1   1 
 lec8_bin 
   0  381 565
   1  14 152
   #Total cases  395 717 
cross_cases(data2, lec9_bin, gender.bin) # other unwanted or uncomfortable sexual experiences
 gender.bin 
 -1   1 
 lec9_bin 
   0  345 396
   1  50 321
   #Total cases  395 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 
##      343       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       343  3.11 0.579    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 ######################

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      392  19.8  2.14     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 109429   54715    1301 <2e-16 ***
## Residuals   1109  46636      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.552296 -22.683746 -20.4208455 0.0000000
## WWU-mturk -21.654762 -22.791639 -20.5178849 0.0000000
## WWU-UMN    -0.102466  -1.195175   0.9902435 0.9736552
################# 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      723 0.234 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           395 0.239 0.543  0.200  0.65
## 2 1            717 0.172 0.515  0.150  0.6

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.1961668
mean(data2$Diffscore.win)
## [1] 0.1930356

Missingness

table(data2$lec5_08a)
## 
##                0 Happened to me\n   Happened_to_me 
##              946               69               97
table(data2$lec5_09a)
## 
##                0 Happened to me\n   Happened_to_me 
##              741              160              211
table(data2$Gender)
## 
##                       Male                     Female 
##                        395                        706 
## Genderqueer/non-conforming                Transgender 
##                          7                          4 
##     Other (please specify) 
##                          0
table(data2$SV)
## 
##  -1   1 
## 723 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  46  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  43  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, SV.num, gender.num, Age, Diffscore.win)

(S2.corr<- rcorr(as.matrix(Study2.cor)))
##                Likability.win SV.num gender.num   Age Diffscore.win
## Likability.win           1.00   0.06       0.11 -0.08         -0.65
## SV.num                   0.06   1.00       0.34 -0.07         -0.10
## gender.num               0.11   0.34       1.00 -0.06         -0.06
## Age                     -0.08  -0.07      -0.06  1.00          0.17
## Diffscore.win           -0.65  -0.10      -0.06  0.17          1.00
## 
## n= 1112 
## 
## 
## P
##                Likability.win SV.num gender.num Age    Diffscore.win
## Likability.win                0.0364 0.0002     0.0070 0.0000       
## SV.num         0.0364                0.0000     0.0255 0.0007       
## gender.num     0.0002         0.0000            0.0633 0.0346       
## Age            0.0070         0.0255 0.0633            0.0000       
## Diffscore.win  0.0000         0.0007 0.0346     0.0000

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)
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     Age + gender.bin + SV + gender.bin * SV, data = data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.61964 -0.30202 -0.01634  0.33880  1.37265 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.173321   0.042724  74.274  < 2e-16 ***
## ASA.ending.fac1 -0.194026   0.020811  -9.323  < 2e-16 ***
## ASA.ending.fac2  0.081744   0.020827   3.925 9.22e-05 ***
## CSA.ending.fac1 -0.233713   0.020846 -11.211  < 2e-16 ***
## CSA.ending.fac2  0.102049   0.020939   4.874 1.26e-06 ***
## Age             -0.003086   0.001249  -2.471   0.0136 *  
## gender.bin1      0.154116   0.036639   4.206 2.81e-05 ***
## SV1              0.099205   0.073178   1.356   0.1755    
## gender.bin1:SV1 -0.075436   0.082073  -0.919   0.3582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.49 on 1103 degrees of freedom
## Multiple R-squared:  0.1757, Adjusted R-squared:  0.1697 
## F-statistic: 29.39 on 8 and 1103 DF,  p-value: < 2.2e-16
confint(Study2.reg)
##                        2.5 %       97.5 %
## (Intercept)      3.089490571  3.257151084
## ASA.ending.fac1 -0.234859137 -0.153193551
## ASA.ending.fac2  0.040877977  0.122609412
## CSA.ending.fac1 -0.274615303 -0.192811093
## CSA.ending.fac2  0.060964424  0.143132593
## Age             -0.005536288 -0.000635381
## gender.bin1      0.082226508  0.226004704
## SV1             -0.044377460  0.242788291
## gender.bin1:SV1 -0.236471920  0.085599977
library(lm.beta)
Study2.stan<- lm.beta(Study2.reg)
Study2.stan
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     Age + gender.bin + SV + gender.bin * SV, data = data2)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac1 ASA.ending.fac2 CSA.ending.fac1 CSA.ending.fac2 
##              NA     -0.29508260      0.12431890     -0.35735348      0.15455243 
##             Age     gender.bin1             SV1 gender.bin1:SV1 
##     -0.06801751      0.13722948      0.08802824     -0.06450405
library(car)
vif(Study2.reg, type=c("predictor")) # variance inflation factors are all low
## GVIFs computed for predictors
##                    GVIF Df GVIF^(1/(2*Df)) Interacts With
## ASA.ending.fac 1.003630  2        1.000906           --  
## CSA.ending.fac 1.013392  2        1.003331           --  
## Age            1.014002  1        1.006977           --  
## gender.bin     1.028972  3        1.004771             SV
## SV             1.028972  3        1.004771     gender.bin
##                                              Other Predictors
## ASA.ending.fac            CSA.ending.fac, Age, gender.bin, SV
## CSA.ending.fac            ASA.ending.fac, Age, gender.bin, SV
## Age            ASA.ending.fac, CSA.ending.fac, gender.bin, SV
## gender.bin                ASA.ending.fac, CSA.ending.fac, Age
## SV                        ASA.ending.fac, CSA.ending.fac, Age

Assumption Checks

plot(Study2.reg)

library(car)
ncvTest(Study2.reg) #this function tests whether data meets assumption of homoscedasticity, where the null hypothesis is that it does meet this assumption. So we reject the null hypothesis, meaning that there is non-constant variance as you move along the regression line
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 2.900798, Df = 1, p = 0.088536

Exploratory Analyses

SV

SV.2.reg<- lm(Likability.win ~ ASA.ending.fac + CSA.ending.fac + Age + SV, data=data2)
summary(SV.2.reg)
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     Age + SV, data = data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.54015 -0.30877 -0.00484  0.32979  1.44167 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.257043   0.038378  84.868  < 2e-16 ***
## ASA.ending.fac1 -0.193081   0.020966  -9.209  < 2e-16 ***
## ASA.ending.fac2  0.080049   0.020969   3.817 0.000142 ***
## CSA.ending.fac1 -0.226218   0.020926 -10.810  < 2e-16 ***
## CSA.ending.fac2  0.100297   0.021087   4.756 2.23e-06 ***
## Age             -0.003178   0.001253  -2.537 0.011324 *  
## SV1              0.085459   0.031178   2.741 0.006225 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4937 on 1105 degrees of freedom
## Multiple R-squared:  0.1616, Adjusted R-squared:  0.1571 
## F-statistic: 35.51 on 6 and 1105 DF,  p-value: < 2.2e-16
confint(SV.2.reg)
##                        2.5 %        97.5 %
## (Intercept)      3.181741350  3.3323436702
## ASA.ending.fac1 -0.234219234 -0.1519421011
## ASA.ending.fac2  0.038904730  0.1211929696
## CSA.ending.fac1 -0.267278056 -0.1851581962
## CSA.ending.fac2  0.058921923  0.1416717065
## Age             -0.005635949 -0.0007199353
## SV1              0.024283068  0.1466343736
SV.2.reg.stan<- lm.beta(SV.2.reg)
SV.2.reg.stan
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     Age + SV, data = data2)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac1 ASA.ending.fac2 CSA.ending.fac1 CSA.ending.fac2 
##              NA     -0.29364438      0.12174132     -0.34589332      0.15189949 
##             Age             SV1 
##     -0.07004774      0.07583034

Gender

Gender.2.reg<- lm(Likability.win ~ ASA.ending.fac + CSA.ending.fac + Age + gender.bin, data=data2)
summary(Gender.2.reg)
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     Age + gender.bin, data = data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.63051 -0.30862 -0.01518  0.33588  1.35878 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.185846   0.041903  76.028  < 2e-16 ***
## ASA.ending.fac1 -0.193918   0.020811  -9.318  < 2e-16 ***
## ASA.ending.fac2  0.082371   0.020821   3.956 8.10e-05 ***
## CSA.ending.fac1 -0.232090   0.020819 -11.148  < 2e-16 ***
## CSA.ending.fac2  0.101096   0.020930   4.830 1.56e-06 ***
## Age             -0.003062   0.001243  -2.464   0.0139 *  
## gender.bin1      0.152074   0.030911   4.920 9.98e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.49 on 1105 degrees of freedom
## Multiple R-squared:  0.174,  Adjusted R-squared:  0.1695 
## F-statistic:  38.8 on 6 and 1105 DF,  p-value: < 2.2e-16
confint(Gender.2.reg)
##                        2.5 %        97.5 %
## (Intercept)      3.103626577  3.2680652562
## ASA.ending.fac1 -0.234752718 -0.1530840795
## ASA.ending.fac2  0.041518626  0.1232231278
## CSA.ending.fac1 -0.272939432 -0.1912400758
## CSA.ending.fac2  0.060029491  0.1421625412
## Age             -0.005500201 -0.0006239028
## gender.bin1      0.091422481  0.2127259758
Gender.2.reg.stan<- lm.beta(Gender.2.reg)
Gender.2.reg.stan
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     Age + gender.bin, data = data2)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac1 ASA.ending.fac2 CSA.ending.fac1 CSA.ending.fac2 
##              NA      -0.2949184       0.1252727      -0.3548712       0.1531099 
##             Age     gender.bin1 
##      -0.0674933       0.1354118

Sample

Conducting analyses (minus interaction term) on each of the separate samples.

UMN.2 <- data2 %>% 
dplyr::filter(Sample== "UMN")

WWU.2 <- data2 %>% 
dplyr::filter(Sample== "WWU")

mturk.2 <- data2 %>% 
dplyr::filter(Sample== "mturk")


UMN.2.reg<- lm(Likability.win ~ ASA.ending.fac + CSA.ending.fac + gender.bin + SV, data=UMN.2)
summary(UMN.2.reg)
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     gender.bin + SV, data = UMN.2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.56040 -0.27159 -0.02679  0.30082  1.16806 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.15321    0.04034  78.163  < 2e-16 ***
## ASA.ending.fac1 -0.17985    0.03318  -5.420 1.05e-07 ***
## ASA.ending.fac2  0.07462    0.03290   2.268  0.02388 *  
## CSA.ending.fac1 -0.23826    0.03272  -7.281 1.88e-12 ***
## CSA.ending.fac2  0.10411    0.03319   3.137  0.00184 ** 
## gender.bin1      0.06408    0.05251   1.220  0.22305    
## SV1              0.14254    0.05198   2.742  0.00639 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4594 on 385 degrees of freedom
## Multiple R-squared:  0.192,  Adjusted R-squared:  0.1795 
## F-statistic: 15.25 on 6 and 385 DF,  p-value: 1.077e-15
confint(UMN.2.reg)
##                        2.5 %     97.5 %
## (Intercept)      3.073892720  3.2325261
## ASA.ending.fac1 -0.245090766 -0.1146037
## ASA.ending.fac2  0.009934028  0.1393046
## CSA.ending.fac1 -0.302605132 -0.1739247
## CSA.ending.fac2  0.038862855  0.1693658
## gender.bin1     -0.039156782  0.1673211
## SV1              0.040343459  0.2447348
UMN.2.reg.stan<- lm.beta(UMN.2.reg)
UMN.2.reg.stan
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     gender.bin + SV, data = UMN.2)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac1 ASA.ending.fac2 CSA.ending.fac1 CSA.ending.fac2 
##              NA     -0.28964719      0.12136015     -0.38751183      0.16643338 
##     gender.bin1             SV1 
##      0.06001693      0.13419574
WWU.2.reg<- lm(Likability.win ~ ASA.ending.fac + CSA.ending.fac + gender.bin + SV, data=WWU.2)
summary(WWU.2.reg)
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     gender.bin + SV, data = WWU.2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.32878 -0.25313 -0.00215  0.29364  1.12268 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.04050    0.04192  72.534  < 2e-16 ***
## ASA.ending.fac1 -0.16391    0.03119  -5.256 2.47e-07 ***
## ASA.ending.fac2  0.08125    0.03079   2.638 0.008673 ** 
## CSA.ending.fac1 -0.21382    0.03110  -6.875 2.58e-11 ***
## CSA.ending.fac2  0.10330    0.03108   3.323 0.000977 ***
## gender.bin1      0.23987    0.05234   4.583 6.25e-06 ***
## SV1             -0.02531    0.04755  -0.532 0.594796    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4276 on 377 degrees of freedom
## Multiple R-squared:  0.1938, Adjusted R-squared:  0.1809 
## F-statistic:  15.1 on 6 and 377 DF,  p-value: 1.656e-15
confint(WWU.2.reg)
##                       2.5 %      97.5 %
## (Intercept)      2.95807417  3.12291935
## ASA.ending.fac1 -0.22523073 -0.10258774
## ASA.ending.fac2  0.02069863  0.14179389
## CSA.ending.fac1 -0.27497024 -0.15266577
## CSA.ending.fac2  0.04218137  0.16442053
## gender.bin1      0.13695221  0.34278172
## SV1             -0.11880668  0.06818137
WWU.2.reg.stan<- lm.beta(WWU.2.reg)
WWU.2.reg.stan
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     gender.bin + SV, data = WWU.2)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac1 ASA.ending.fac2 CSA.ending.fac1 CSA.ending.fac2 
##              NA     -0.28253349      0.14114420     -0.37213295      0.17906949 
##     gender.bin1             SV1 
##      0.22792668     -0.02657743
mturk.2.reg<-  lm(Likability.win ~ ASA.ending.fac + CSA.ending.fac + gender.bin + SV, data=mturk.2)
summary(mturk.2.reg)
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     gender.bin + SV, data = mturk.2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.54552 -0.38468  0.01376  0.44032  1.33560 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.09392    0.04879  63.409  < 2e-16 ***
## ASA.ending.fac1 -0.23593    0.04493  -5.251 2.72e-07 ***
## ASA.ending.fac2  0.08694    0.04565   1.905   0.0577 .  
## CSA.ending.fac1 -0.24090    0.04573  -5.268 2.50e-07 ***
## CSA.ending.fac2  0.10068    0.04538   2.219   0.0272 *  
## gender.bin1      0.12445    0.06697   1.858   0.0640 .  
## SV1             -0.01663    0.07598  -0.219   0.8269    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.585 on 329 degrees of freedom
## Multiple R-squared:  0.1573, Adjusted R-squared:  0.1419 
## F-statistic: 10.23 on 6 and 329 DF,  p-value: 2.173e-10
confint(mturk.2.reg)
##                        2.5 %     97.5 %
## (Intercept)      2.997930721  3.1899032
## ASA.ending.fac1 -0.324311069 -0.1475460
## ASA.ending.fac2 -0.002858632  0.1767301
## CSA.ending.fac1 -0.330862646 -0.1509343
## CSA.ending.fac2  0.011407620  0.1899584
## gender.bin1     -0.007302459  0.2561972
## SV1             -0.166104131  0.1328459
mturk.2.reg.stan<- lm.beta(mturk.2.reg)
mturk.2.reg.stan
## 
## Call:
## lm(formula = Likability.win ~ ASA.ending.fac + CSA.ending.fac + 
##     gender.bin + SV, data = mturk.2)
## 
## Standardized Coefficients::
##     (Intercept) ASA.ending.fac1 ASA.ending.fac2 CSA.ending.fac1 CSA.ending.fac2 
##              NA     -0.30812792      0.11129584     -0.31330943      0.13006722 
##     gender.bin1             SV1 
##      0.09833255     -0.01150811
cor.test(mturk.2$Likability, mturk.2$Age)
## 
##  Pearson's product-moment correlation
## 
## data:  mturk.2$Likability and mturk.2$Age
## t = -1.1102, df = 334, p-value = 0.2677
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.16655052  0.04666016
## sample estimates:
##         cor 
## -0.06063681
cor.test(UMN.2$Likability, UMN.2$Age)
## 
##  Pearson's product-moment correlation
## 
## data:  UMN.2$Likability and UMN.2$Age
## t = -0.52736, df = 390, p-value = 0.5982
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.12541123  0.07254558
## sample estimates:
##         cor 
## -0.02669453
cor.test(WWU.2$Likability, WWU.2$Age)
## 
##  Pearson's product-moment correlation
## 
## data:  WWU.2$Likability and WWU.2$Age
## t = -1.3116, df = 382, p-value = 0.1905
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.16591891  0.03334453
## sample estimates:
##         cor 
## -0.06695478

Difference Scores

Linear Regression

S2.reg.diffScores<- lm(Diffscore.win ~ Age + gender.bin + SV + gender.bin*SV, data=data2)
summary(S2.reg.diffScores)
## 
## Call:
## lm(formula = Diffscore.win ~ Age + gender.bin + SV + gender.bin * 
##     SV, data = data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.49130 -0.32268 -0.02686  0.32314  1.40991 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.055287   0.043583   1.269    0.205    
## Age              0.006754   0.001277   5.291 1.47e-07 ***
## gender.bin1     -0.011513   0.037309  -0.309    0.758    
## SV1             -0.018479   0.074658  -0.248    0.805    
## gender.bin1:SV1 -0.085946   0.083758  -1.026    0.305    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5009 on 1107 degrees of freedom
## Multiple R-squared:  0.03748,    Adjusted R-squared:  0.034 
## F-statistic: 10.78 on 4 and 1107 DF,  p-value: 1.431e-08
confint(S2.reg.diffScores)
##                        2.5 %      97.5 %
## (Intercept)     -0.030228436 0.140802006
## Age              0.004249208 0.009258511
## gender.bin1     -0.084717192 0.061690224
## SV1             -0.164967479 0.128008487
## gender.bin1:SV1 -0.250287991 0.078396267
Study2.diff.stan<- lm.beta(S2.reg.diffScores)
Study2.diff.stan
## 
## Call:
## lm(formula = Diffscore.win ~ Age + gender.bin + SV + gender.bin * 
##     SV, data = data2)
## 
## Standardized Coefficients::
##     (Intercept)             Age     gender.bin1             SV1 gender.bin1:SV1 
##              NA      0.15706173     -0.01081628     -0.01730004     -0.07753606