Study 1a: 1st W-Ss Replication of Berman & Small (2018)’s Study 2

Data Preparation

In this pre-registered study, we collected data from 174 participant. After filtering these participants with three attention check measures, we obtained analyzable data from 166 participants.

Loading Packages

Before running this chunk, please load “Study 1a.csv” into the R environment.

# packages should be loaded in the following order to avoid function conflicts
options(scipen = 99) # turns off scientific notation
library(psych) # for describing data
library(effsize) # for mean difference effect sizes
library(sjstats) # for eta-squared effect sizes
library(correlation) # for cleaner correlation test output
library(rmcorr) # for repeated-measures correlation tests
library(cocor) # for comparing dependent correlation coefficients
library(tidyverse) # for data manipulation and plotting

Attention Check

Participants were instructed to select zero (i.e., the middle option in the 9-point scale) as their response for three measures that followed a modified scenario. If a participant failed to select “0” as their response for any of the three measures, they would be excluded from the analyses.

# Experimental task attention check variable
Study_1a$AttChecks <- rowSums(Study_1a[, c("AttnCheck1", "AttnCheck2", "AttnCheck3")], na.rm = T)

# The exclusion of any sum (within the variable "AttChecks") that did not result in a zero
Study_1a <- Study_1a %>% 
  filter(Study_1a$AttChecks == 0)

Variable Creation

The original design provided by Berman & Small required virtuosity, honorability, and respectability measures to be averaged across for all participants. Moreover, since the study had a between-subjects design, it included an overall “Virtuosity Average”. For this replication, the within-subjects design does not allow to have an overall average, so there will be a moral and a non-moral average for each participant. Additionally, there will be a virtuosity average for each vignette.

# Create a variable for virtuosity average in the moral condition
Study_1a$MoralAverage <- rowMeans(Study_1a[, c("AdulteryVirtuous1", "AdulteryHonorable2", "AdulteryRespectable3", "OfficeVirtuous1", "OfficeHonorable2", "OfficeRespectable3", "KidneyVirtuous1", "KidneyHonorable2", "KidneyRespectable3", "AdulteryVirtuous2", "AdulteryHonorable3", "AdulteryRespectable1", "OfficeVirtuous2", "OfficeHonorable3", "OfficeRespectable1", "KidneyVirtuous2", "KidneyHonorable3", "KidneyRespectable1", "AdulteryVirtuous3", "AdulteryHonorable1", "AdulteryRespectable2", "OfficeVirtuous3", "OfficeHonorable1", "OfficeRespectable2", "KidneyVirtuous3", "KidneyHonorable1", "KidneyRespectable2" )], na.rm = T)

# Create a variable for virtuosity average in the non-moral condition
Study_1a$NonmoralAverage <- rowMeans(Study_1a[, c("DessertVirtuous1", "DessertHonorable2", "DessertRespectable3", "GymVirtuous1", "GymHonorable2", "GymRespectable3", "SchoolworkVirtuous1", "SchoolworkHonorable2", "SchoolworkRespctable3", "DessertVirtuous2", "DessertHonorable3", "DessertRespectable1", "GymVirtuous2", "GymHonorable3", "GymRespectable1", "SchoolworkVirtuous2", "SchoolworkHonorable3", "SchoolworkRespectable1", "DessertVirtuous3", "DessertHonorable1", "DessertRespectable2", "GymVirtuous3", "GymHonorable1", "GymRespectable2", "SchoolworkVirtuous3", "SchoolworkHonorable1", "SchoolworkRespectable2")], na.rm = T)

# Create a variable for virtuosity average in the Office Supplies scenario
Study_1a$OfficeAverage <- rowMeans(Study_1a[, c("OfficeVirtuous1", "OfficeHonorable2", "OfficeRespectable3", "OfficeVirtuous2", "OfficeHonorable3", "OfficeRespectable1", "OfficeVirtuous3", "OfficeHonorable1", "OfficeRespectable2" )], na.rm = T)

# Create a variable for virtuosity average in the Adultery scenario
Study_1a$AdulteryAverage <- rowMeans(Study_1a[, c("AdulteryVirtuous1", "AdulteryHonorable2", "AdulteryRespectable3", "AdulteryVirtuous2", "AdulteryHonorable3", "AdulteryRespectable1", "AdulteryVirtuous3", "AdulteryHonorable1", "AdulteryRespectable2" )], na.rm = T)

# Create a variable for virtuosity average in the Kidney scenario
Study_1a$KidneyAverage <- rowMeans(Study_1a[, c("KidneyVirtuous1", "KidneyHonorable2", "KidneyRespectable3", "KidneyVirtuous2", "KidneyHonorable3", "KidneyRespectable1", "KidneyVirtuous3", "KidneyHonorable1", "KidneyRespectable2")], na.rm = T)

# Create a variable for virtuosity average in the Dessert scenario
Study_1a$DessertAverage <- rowMeans(Study_1a[, c("DessertVirtuous1", "DessertHonorable2", "DessertRespectable3", "DessertVirtuous2", "DessertHonorable3", "DessertRespectable1", "DessertVirtuous3", "DessertHonorable1", "DessertRespectable2")], na.rm = T)

# Create a variable for virtuosity average in the Gym scenario
Study_1a$GymAverage <- rowMeans(Study_1a[, c("GymVirtuous1", "GymHonorable2", "GymRespectable3", "GymVirtuous2", "GymHonorable3", "GymRespectable1", "GymVirtuous3", "GymHonorable1", "GymRespectable2")], na.rm = T)

# Create a variable for virtuosity average in the Schoolwork scenario
Study_1a$SchoolworkAverage <- rowMeans(Study_1a[, c("SchoolworkVirtuous1", "SchoolworkHonorable2", "SchoolworkRespctable3", "SchoolworkVirtuous2", "SchoolworkHonorable3", "SchoolworkRespectable1","SchoolworkVirtuous3", "SchoolworkHonorable1", "SchoolworkRespectable2")], na.rm = T)

B&S (2018)’s Group-Lvl Analyses

Simple Effects

For all t-tests, we conduct both one-tailed and two-tailed tests. Our pre-registration specified directional hypotheses which therefore ought to be analyzed with one-tailed tests. However, in order for readers to know any estimate’s uncertainty in both directions, we also conduct two-tailed tests to get relevant CIs.

Non-Moral Vignettes

Overall

# Non-Moral t-test using person-level averages across Non-Moral vignettes
Study_1a_nonmoral_t <- t.test(Study_1a$NonmoralAverage, mu = 0, alternative = "less") 
## 'greater' is for one-tailed p-value because of the tempted > non-tempted hypothesis (negative numbers mean tempted > non-tempted)
print(Study_1a_nonmoral_t)

    One Sample t-test

data:  Study_1a$NonmoralAverage
t = -0.65805, df = 165, p-value = 0.2557
alternative hypothesis: true mean is less than 0
95 percent confidence interval:
       -Inf 0.09726628
sample estimates:
  mean of x 
-0.06425703 

Study_1a_nonmoral_t_2side <- t.test(Study_1a$NonmoralAverage, mu = 0, alternative = "two.sided") 
## 'two.sided' is to get uncertainty estimates in both directions
print(Study_1a_nonmoral_t_2side)

    One Sample t-test

data:  Study_1a$NonmoralAverage
t = -0.65805, df = 165, p-value = 0.5114
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.2570574  0.1285433
sample estimates:
  mean of x 
-0.06425703 

# returns cohen's d
effsize::cohen.d(Study_1a$NonmoralAverage,
                 f = NA,
                 mu = 0,
                 data = data)

Cohen's d (single sample)

d estimate: -0.05107448 (negligible)
Reference mu: 0
95 percent confidence interval:
     lower      upper 
-0.3576178  0.2554689 

Dessert Vignette

# Dessert vignette t-test
Study_1a_dessert_t <- t.test(Study_1a$DessertAverage, mu = 0, alternative = "less") 
## 'greater' is for one-tailed p-value because of the tempted > non-tempted hypothesis (negative numbers mean tempted > non-tempted)
print(Study_1a_dessert_t)

    One Sample t-test

data:  Study_1a$DessertAverage
t = -1.001, df = 165, p-value = 0.1591
alternative hypothesis: true mean is less than 0
95 percent confidence interval:
       -Inf 0.07206212
sample estimates:
 mean of x 
-0.1104418 

Study_1a_dessert_t_2side <- t.test(Study_1a$DessertAverage, mu = 0, alternative = "two.sided") 
## 'two.sided' is to get uncertainty estimates in both directions
print(Study_1a_dessert_t_2side)

    One Sample t-test

data:  Study_1a$DessertAverage
t = -1.001, df = 165, p-value = 0.3183
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.3282854  0.1074018
sample estimates:
 mean of x 
-0.1104418 

# returns cohen's d
effsize::cohen.d(Study_1a$DessertAverage,
                 f = NA,
                 mu = 0,
                 data = data)

Cohen's d (single sample)

d estimate: -0.07769261 (negligible)
Reference mu: 0
95 percent confidence interval:
     lower      upper 
-0.3843016  0.2289164 

Gym Vignette

# Gym vignette t-test
Study_1a_gym_t <- t.test(Study_1a$GymAverage, mu = 0, alternative = "less") 
## 'greater' is for one-tailed p-value because of the tempted > non-tempted hypothesis (negative numbers mean tempted > non-tempted)
print(Study_1a_gym_t)

    One Sample t-test

data:  Study_1a$GymAverage
t = -0.41237, df = 165, p-value = 0.3403
alternative hypothesis: true mean is less than 0
95 percent confidence interval:
      -Inf 0.1330292
sample estimates:
  mean of x 
-0.04417671 

Study_1a_gym_t_2side <- t.test(Study_1a$GymAverage, mu = 0, alternative = "two.sided") 
## 'two.sided' is to get uncertainty estimates in both directions
print(Study_1a_gym_t_2side)

    One Sample t-test

data:  Study_1a$GymAverage
t = -0.41237, df = 165, p-value = 0.6806
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.2556964  0.1673430
sample estimates:
  mean of x 
-0.04417671 

# returns cohen's d
effsize::cohen.d(Study_1a$GymAverage,
                 f = NA,
                 mu = 0,
                 data = data)

Cohen's d (single sample)

d estimate: -0.03200616 (negligible)
Reference mu: 0
95 percent confidence interval:
     lower      upper 
-0.3385192  0.2745069 

Schoolwork Vignette

# Schoolwork vignette t-test
Study_1a_schoolwork_t <- t.test(Study_1a$SchoolworkAverage, mu = 0, alternative = "less") 
## 'greater' is for one-tailed p-value because of the tempted > non-tempted hypothesis (negative numbers mean tempted > non-tempted)
print(Study_1a_schoolwork_t)

    One Sample t-test

data:  Study_1a$SchoolworkAverage
t = -0.32359, df = 165, p-value = 0.3733
alternative hypothesis: true mean is less than 0
95 percent confidence interval:
      -Inf 0.1568773
sample estimates:
  mean of x 
-0.03815261 

Study_1a_schoolwork_t_2side <- t.test(Study_1a$SchoolworkAverage, mu = 0, alternative = "two.sided") 
## 'two.sided' is to get uncertainty estimates in both directions
print(Study_1a_schoolwork_t_2side)

    One Sample t-test

data:  Study_1a$SchoolworkAverage
t = -0.32359, df = 165, p-value = 0.7467
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.2709477  0.1946425
sample estimates:
  mean of x 
-0.03815261 

# returns cohen's d
effsize::cohen.d(Study_1a$SchoolworkAverage,
                 f = NA,
                 mu = 0,
                 data = data)

Cohen's d (single sample)

d estimate: -0.02511548 (negligible)
Reference mu: 0
95 percent confidence interval:
    lower     upper 
-0.331621  0.281390 

Moral Vignettes

Overall

# Moral t-test using person-level averages across Moral vignettes
Study_1a_moral_t <- t.test(Study_1a$MoralAverage, mu = 0, alternative = "greater") 
## 'less' is for one-tailed p-value because of the non-tempted > tempted hypothesis (positive numbers mean non-tempted > tempted)
print(Study_1a_moral_t)

    One Sample t-test

data:  Study_1a$MoralAverage
t = 10.276, df = 165, p-value < 0.00000000000000022
alternative hypothesis: true mean is greater than 0
95 percent confidence interval:
 0.9861768       Inf
sample estimates:
mean of x 
 1.175368 

Study_1a_moral_t_2side <- t.test(Study_1a$MoralAverage, mu = 0, alternative = "two.sided") 
## 'two.sided' is to get uncertainty estimates in both directions
print(Study_1a_moral_t_2side)

    One Sample t-test

data:  Study_1a$MoralAverage
t = 10.276, df = 165, p-value < 0.00000000000000022
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.9495421 1.4011941
sample estimates:
mean of x 
 1.175368 

# returns cohen's d
effsize::cohen.d(Study_1a$MoralAverage,
                 f = NA,
                 mu = 0,
                 data = data)

Cohen's d (single sample)

d estimate: 0.797611 (medium)
Reference mu: 0
95 percent confidence interval:
   lower    upper 
0.479164 1.116058 

Office Vignette

# Office vignette t-test
Study_1a_office_t <- t.test(Study_1a$OfficeAverage, mu = 0, alternative = "greater") 
## 'less' is for one-tailed p-value because of the non-tempted > tempted hypothesis (positive numbers mean non-tempted > tempted)
print(Study_1a_office_t)

    One Sample t-test

data:  Study_1a$OfficeAverage
t = 9.568, df = 165, p-value < 0.00000000000000022
alternative hypothesis: true mean is greater than 0
95 percent confidence interval:
 1.091196      Inf
sample estimates:
mean of x 
 1.319277 

Study_1a_office_t_2side <- t.test(Study_1a$OfficeAverage, mu = 0, alternative = "two.sided") 
## 'two.sided' is to get uncertainty estimates in both directions
print(Study_1a_office_t_2side)

    One Sample t-test

data:  Study_1a$OfficeAverage
t = 9.568, df = 165, p-value < 0.00000000000000022
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 1.047031 1.591523
sample estimates:
mean of x 
 1.319277 

# returns cohen's d
effsize::cohen.d(Study_1a$OfficeAverage,
                 f = NA,
                 mu = 0,
                 data = data)

Cohen's d (single sample)

d estimate: 0.742618 (medium)
Reference mu: 0
95 percent confidence interval:
    lower     upper 
0.4257365 1.0594994 

Adultery Vignette

# Adultery vignette t-test
Study_1a_adultery_t <- t.test(Study_1a$AdulteryAverage, mu = 0, alternative = "greater") 
## 'less' is for one-tailed p-value because of the non-tempted > tempted hypothesis (positive numbers mean non-tempted > tempted)
print(Study_1a_adultery_t)

    One Sample t-test

data:  Study_1a$AdulteryAverage
t = 9.8834, df = 165, p-value < 0.00000000000000022
alternative hypothesis: true mean is greater than 0
95 percent confidence interval:
 1.141947      Inf
sample estimates:
mean of x 
 1.371486 

Study_1a_adultery_t_2side <- t.test(Study_1a$AdulteryAverage, mu = 0, alternative = "two.sided") 
## 'two.sided' is to get uncertainty estimates in both directions
print(Study_1a_adultery_t_2side)

    One Sample t-test

data:  Study_1a$AdulteryAverage
t = 9.8834, df = 165, p-value < 0.00000000000000022
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 1.097499 1.645473
sample estimates:
mean of x 
 1.371486 

# returns cohen's d
effsize::cohen.d(Study_1a$AdulteryAverage,
                 f = NA,
                 mu = 0,
                 data = data)

Cohen's d (single sample)

d estimate: 0.7671017 (medium)
Reference mu: 0
95 percent confidence interval:
    lower     upper 
0.4495362 1.0846673 

Kidney Vignette

# Kidney vignette t-test
Study_1a_kidney_t <- t.test(Study_1a$KidneyAverage, mu = 0, alternative = "greater") 
## 'less' is for one-tailed p-value because of the non-tempted > tempted hypothesis (positive numbers mean non-tempted > tempted)
print(Study_1a_kidney_t)

    One Sample t-test

data:  Study_1a$KidneyAverage
t = 6.0929, df = 165, p-value = 0.000000003781
alternative hypothesis: true mean is greater than 0
95 percent confidence interval:
 0.6085589       Inf
sample estimates:
mean of x 
0.8353414 

Study_1a_kidney_t_2side <- t.test(Study_1a$KidneyAverage, mu = 0, alternative = "two.sided") 
## 'two.sided' is to get uncertainty estimates in both directions
print(Study_1a_kidney_t_2side)

    One Sample t-test

data:  Study_1a$KidneyAverage
t = 6.0929, df = 165, p-value = 0.000000007561
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.5646451 1.1060376
sample estimates:
mean of x 
0.8353414 

# returns cohen's d
effsize::cohen.d(Study_1a$KidneyAverage,
                 f = NA,
                 mu = 0,
                 data = data)

Cohen's d (single sample)

d estimate: 0.4729039 (small)
Reference mu: 0
95 percent confidence interval:
    lower     upper 
0.1621561 0.7836518 

Interaction

Repeated-Measures ANOVA

The purpose of this test is to determine whether there is a difference in mean virtuosity depending on condition (i.e., non-moral vs moral).

# Prepare data by using pivot_longer on a new dataset
Study_1a_anova <- Study_1a %>%
  dplyr::select(c(ResponseId, MoralAverage, NonmoralAverage))%>% 
  group_by(ResponseId) %>%
  pivot_longer(cols = (c(MoralAverage,NonmoralAverage)),
               names_to = "Condition",
               values_to = "VirtuosityAverage")

# Repeated-Meausures ANOVA
Study_1a_interaction_aov <- aov(VirtuosityAverage~Condition + Error(as.factor(ResponseId)), data = Study_1a_anova)
summary(Study_1a_interaction_aov)

Error: as.factor(ResponseId)
           Df Sum Sq Mean Sq F value Pr(>F)
Residuals 165  475.9   2.884               

Error: Within
           Df Sum Sq Mean Sq F value              Pr(>F)    
Condition   1  127.5  127.54   146.6 <0.0000000000000002 ***
Residuals 165  143.6    0.87                                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Paired-Samples t-Test

Replicates results of the repeated-measures ANOVA (Berman & Small used ANOVAs rather than paired-samples t-tests for the interaction test). These tests are statistically identical, but we conduct both in case researchers have preferences for one over the other.

Study_1a_interaction_t <- t.test(VirtuosityAverage ~ Condition, data = Study_1a_anova, alternative = "greater", paired = TRUE)
## 'less' is for one-tailed p-value because of the non-moral > moral hypothesis (positive numbers mean non-moral > moral)
Study_1a_interaction_t

    Paired t-test

data:  VirtuosityAverage by Condition
t = 12.107, df = 165, p-value < 0.00000000000000022
alternative hypothesis: true mean difference is greater than 0
95 percent confidence interval:
 1.070262      Inf
sample estimates:
mean difference 
       1.239625 

Study_1a_interaction_t_2side <- t.test(VirtuosityAverage ~ Condition, data = Study_1a_anova, alternative = "two.sided", paired = TRUE) 
## 'two.sided' is to get uncertainty estimates in both directions
print(Study_1a_interaction_t_2side)

    Paired t-test

data:  VirtuosityAverage by Condition
t = 12.107, df = 165, p-value < 0.00000000000000022
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 1.037467 1.441784
sample estimates:
mean difference 
       1.239625 

# returns dz effect size and 95% CIs
effsize::cohen.d(VirtuosityAverage ~ Condition | Subject(ResponseId),
                 data = Study_1a_anova %>% droplevels(), 
                 paired = T,
                 within = F)

Cohen's d

d estimate: 0.9397003 (large)
95 percent confidence interval:
    lower     upper 
0.7563845 1.1230161 

# returns d-av effect size and 95% CIs
effsize::cohen.d(VirtuosityAverage ~ Condition | Subject(ResponseId),
                 data = Study_1a_anova %>% droplevels(), 
                 paired = T,
                 within = T)

Cohen's d

d estimate: 0.8981868 (large)
95 percent confidence interval:
    lower     upper 
0.7253034 1.0710701 

Person-Lvl Analyses

# Create person-level variables
Study_1a_plvl <- Study_1a %>% 
  mutate(fits_hypothesis_n = if_else(NonmoralAverage < 0, 1, 0)) %>% 
  mutate(fits_hypothesis_dessert = if_else(DessertAverage < 0, 1, 0)) %>% 
  mutate(fits_hypothesis_gym = if_else(GymAverage < 0, 1, 0)) %>% 
  mutate(fits_hypothesis_schoolwork = if_else(SchoolworkAverage < 0, 1, 0)) %>% 
  mutate(fits_hypothesis_m = if_else(MoralAverage > 0, 1, 0)) %>%
  mutate(fits_hypothesis_office = if_else(OfficeAverage > 0, 1, 0)) %>% 
  mutate(fits_hypothesis_adultery = if_else(AdulteryAverage > 0, 1, 0)) %>% 
  mutate(fits_hypothesis_kidney = if_else(KidneyAverage > 0, 1, 0))

# Create variable that indicates whether participants fit the full complex hypothesis or not 
Study_1a_plvl$number_of_hyp_fit <- rowSums(Study_1a_plvl[, c("fits_hypothesis_m", "fits_hypothesis_n")], na.rm = T)
Study_1a_plvl <- Study_1a_plvl %>% 
   mutate(fits_complex = if_else(number_of_hyp_fit == 2, "Yes", "No"))

Descriptive Pervasiveness

Non-Moral Vignettes

# How many participants fit hypothesis for Non-Moral vignettes
table(Study_1a_plvl$fits_hypothesis_n) # (0 = "no"; 1 = "yes")

 0  1 
97 69 

Moral Vignettes

# How many participants fit hypothesis for Moral vignettes
table(Study_1a_plvl$fits_hypothesis_m) # (0 = "no"; 1 = "yes")

  0   1 
 37 129 

Interaction

# How many participants fit complex interaction hypothesis (i.e., in Non-Moral vignettes: tempted > non-tempted; in Moral vignettes: non-tempted > tempted)
table(Study_1a_plvl$fits_complex) # (0 = "no"; 1 = "yes")

 No Yes 
115  51 

Randomization Test

This test is only possible for the complex interaction hypothesis, as randomly shuffling within-condition data would yield the same result across shufflings.

library(mosaic) # for randomly shuffling data

Study_1a_random <- Study_1a_plvl_all %>% 
  dplyr::select(c(ResponseId, MoralAverage, NonmoralAverage))

Study_1a_random <- Study_1a_random %>% 
  mutate("Interaction" = case_when(
    (MoralAverage == 0 & NonmoralAverage == 0) ~ "Zero, Zero, Zero",
    (MoralAverage == 0 & NonmoralAverage < 0) ~ "Zero, Neg, Pos",
    (MoralAverage == 0 & NonmoralAverage > 0) ~ "Zero, Pos, Neg",
    (MoralAverage < 0 & NonmoralAverage == 0) ~ "Neg, Zero, Neg",
    (MoralAverage < 0 & NonmoralAverage < 0 & MoralAverage == NonmoralAverage) ~ "Neg, Neg, Zero",
    (MoralAverage < 0 & NonmoralAverage > 0) ~ "Neg, Pos, Neg",
    (MoralAverage < 0 & NonmoralAverage < 0 & MoralAverage > NonmoralAverage) ~ "Neg, Neg, Pos",
    (MoralAverage < 0 & NonmoralAverage < 0 & MoralAverage < NonmoralAverage) ~ "Neg, Neg, Neg",
    (MoralAverage > 0 & NonmoralAverage == 0) ~ "Pos, Zero, Pos",
    (MoralAverage > 0 & NonmoralAverage < 0) ~ "Pos, Neg, Pos", # this is the predicted pattern
    (MoralAverage > 0 & NonmoralAverage > 0 & MoralAverage == NonmoralAverage) ~ "Pos, Pos, Zero",
    (MoralAverage > 0 & NonmoralAverage > 0 & MoralAverage < NonmoralAverage) ~ "Pos, Pos, Neg",
    (MoralAverage > 0 & NonmoralAverage > 0 & MoralAverage > NonmoralAverage) ~ "Pos, Pos, Pos"))

prop_Study_1a <- prop( ~ Interaction == "Pos, Neg, Pos", data = Study_1a_random) # create descriptive proportion of claimed pattern from Berman & Small
prop_Study_1a_dp <- unname(prop_Study_1a[1]) # get descriptive proportion proportion for interaction

nsim <- 1000 # number of random datasets to generate
nullprop_Study_1a <- rep(NA, nsim) # creating empty bin

for (i in 1:nsim) { # for loop to do random shuffling across nsim datasets and calculate proportions in those randomly shuffled datasets
  #permutation, size equal to number of observations for each condition
  Study_1a_random$NonmoralAverageR <- shuffle(Study_1a_random$NonmoralAverage) # shuffling 'Non-Moral' values  
  Study_1a_random$MoralAverageR <- shuffle(Study_1a_random$MoralAverage) # shuffling 'Moral' values
  
  ## Creating variable that distinguishes Ps who showed predicted (claimed) pattern from Ps who did not
  Study_1a_R <- Study_1a_random %>%
  mutate("InteractionR" = case_when(
    (MoralAverageR == 0 & NonmoralAverageR == 0) ~ "Zero, Zero, Zero",
    (MoralAverageR == 0 & NonmoralAverageR < 0) ~ "Zero, Neg, Pos",
    (MoralAverageR == 0 & NonmoralAverageR > 0) ~ "Zero, Pos, Neg",
    (MoralAverageR < 0 & NonmoralAverageR == 0) ~ "Neg, Zero, Neg",
    (MoralAverageR < 0 & NonmoralAverageR < 0 & MoralAverageR == NonmoralAverageR) ~ "Neg, Neg, Zero",
    (MoralAverageR < 0 & NonmoralAverageR > 0) ~ "Neg, Pos, Neg",
    (MoralAverageR < 0 & NonmoralAverageR < 0 & MoralAverageR > NonmoralAverageR) ~ "Neg, Neg, Pos",
    (MoralAverageR < 0 & NonmoralAverageR < 0 & MoralAverageR < NonmoralAverageR) ~ "Neg, Neg, Neg",
    (MoralAverageR > 0 & NonmoralAverageR == 0) ~ "Pos, Zero, Pos",
    (MoralAverageR > 0 & NonmoralAverageR < 0) ~ "Pos, Neg, Pos", # this is the predicted pattern
    (MoralAverageR > 0 & NonmoralAverageR > 0 & MoralAverageR == NonmoralAverageR) ~ "Pos, Pos, Zero",
    (MoralAverageR > 0 & NonmoralAverageR > 0 & MoralAverageR < NonmoralAverageR) ~ "Pos, Pos, Neg",
    (MoralAverageR > 0 & NonmoralAverageR > 0 & MoralAverageR > NonmoralAverageR) ~ "Pos, Pos, Pos"))

  
  nullprop_TF_Study_1a <-prop( ~ InteractionR == "Pos, Neg, Pos", data = Study_1a_R) # create proportion value of claimed direction 
  nullprop_Study_1a[i]<- unname(nullprop_TF_Study_1a[1]) # get claimed pattern proportion for randomly shuffled data and repeat for each nsim dataset
}

gf_histogram( ~ nullprop_Study_1a, fill = ~ (nullprop_Study_1a >= prop_Study_1a_dp)) # create histogram of how many randomly shuffled datasets lead to a prop value similar to or greater than empirical prop

#c(chance)-value for empirical proportion (similar to Grice's c-values computed in his OOM software)
prop( ~ nullprop_Study_1a >= prop_Study_1a_dp) # suggests that empirical interaction proportion is very likely to have occurred by chance

prop_TRUE 
    0.871 

Frequentist Prevalence

# CREATE NECESSARY FREQUENTIST FUNCTION (code adapted from Donhauser et al.'s (2018) MATLAB code)

prevalence_test <- function(observed, alpha_ind=0.05, beta_ind=1, alpha_group=0.05, gamma_0=0.5) {
  
  # Inputs:
  # observed    = vector of p-values OR vector with two values, e.g., 'c(positive cases, total cases)'
  # alpha_ind   = alpha threshold for person-level tests, if 'observed' = vector of p-values (typically set to .05)
  # beta_ind    = sensitivity of person-level tests (probably set to 1.00)
  # alpha_group = see "outputs"
  # gamma_0     = null gamma value for empirical prevalence to be tested against (defaults to majority null)
  
  # Outputs:
  # p_null      = p-value for majority null or whatever null you specify in the gamma_0 argument (e.g., to specify a global null set "gamma_0 = 0") 
  # gamma_0     = highest gamma_0 value that can be rejected at threshold of 'alpha_group' given positive cases. Note this output is not the same as the input argument
  if(sum(observed > 1) == 0) {
    pvals <- observed
  } else {
    pvals <- c(rep(0, observed[1]), rep(1, observed[2]-observed[1]))
  }
  
  n <- length(pvals)
  k_obs <- sum(pvals < alpha_ind)
  kvals <- 0:n
  dgamma <- 0.001
  gammavals <- seq(0, 1, dgamma)
  
  p_k_gamma <- matrix(0, nrow=length(gammavals), ncol=length(kvals))
  p_kk_gamma <- matrix(0, nrow=length(gammavals), ncol=length(kvals))
  
  for(iGam in 1:length(gammavals)) {
    p_pos <- gammavals[iGam] * beta_ind + (1 - gammavals[iGam]) * alpha_ind
    p_neg <- gammavals[iGam] * (1 - beta_ind) + (1 - gammavals[iGam]) * (1 - alpha_ind)
    for(iK in 1:length(kvals)) {
      k <- kvals[iK]
      tmp <- choose(n, k) * p_pos^k * p_neg^(n - k)
      p_k_gamma[iGam, iK] <- tmp
      p_kk_gamma[iGam, iK] <- 1 - sum(p_k_gamma[iGam, 1:iK]) + tmp
    }
  }
  
  iK <- which(kvals == k_obs)
  p_null <- p_kk_gamma[which(gammavals == gamma_0), iK]
  
  index <- sum(p_kk_gamma[, iK] < alpha_group)
  if(index != 0) {
    gamma_0 <- gammavals[index]
  } else {
    gamma_0 <- 0
  }
  
  return(list(p_null=p_null, gamma_0=gamma_0, p_kk_gamma=p_kk_gamma, gammavals=gammavals, p_k_gamma=p_k_gamma))
}

Non-Moral Vignettes

k <- 69  # number of persons matching hypothesized / group-level pattern
n <- 69+97 # total N
nonmoral_freqPrevalence <- prevalence_test(c(k, n)) # specify positive cases out of total cases for prev test function
nonmoral_freqPrevalence$p_null # print p-value testing against majority null

[1] 0.9981335

# Originally printed: 0.9981335

Moral Vignettes

k <- 129  # number of persons matching hypothesized / group-level pattern
n <- 129+37 # total N
moral_freqPrevalence <- prevalence_test(c(k, n)) # specify positive cases out of total cases for prev test function
moral_freqPrevalence$p_null # print p-value testing against majority null

[1] 0.00000000001652683

# Originally printed: 0.00000000001652683

Interaction

k <- 51  # number of persons matching hypothesized / group-level pattern
n <- 51+115 # total N
primary_freqPrevalence <- prevalence_test(c(k, n)) # specify positive cases out of total cases for prev test function
primary_freqPrevalence$p_null # print p-value testing against majority null

[1] 1

# Originally printed: 0.9981335

Bayesian Prevalence

# CREATE NECESSARY BAYES FUNCTIONS (code from Ince et al., 2021)
library(nleqslv)

bayesprev_map <- function(k, n, a=0.05, b=1) {
  # Bayesian maximum a posteriori estimate of population prevalence gamma
  # under a uniform prior
  # 
  # Args:
  #  k: number of participants/tests significant out of 
  #  n: total number of participants/tests
  #  a: alpha value of within-participant test (default=0.05)
  #  b: sensitivity/beta of within-participant test (default=1)
  
  gm <- (k/n -a)/(b-a)
  if(gm <0) gm <- 0
  if(gm>1) gm <- 1
  return(gm)
} 

bayesprev_posterior <- function(x, k, n, a=0.05, b=1) {
  # Bayesian posterior of population prevalence gamma under a uniform prior
  #
  # Args:
  # x : values of gamma at which to evaluate the posterior density
  # k : number of participants significant out of 
  # n : total number of participants
  # a : alpha value of within-participant test (default=0.05)
  # b : sensitivity/beta of within-participant test (default=1)
  
  m1 <-  k + 1
  m2 <- n - k + 1
  theta <- a + (b-a)*x
  post <- (b -a)*dbeta(theta,m1, m2)
  post <- post/(pbeta(b, m1, m2) - pbeta(a, m1, m2))
  return(post)
}


bayesprev_bound <- function(p, k, n, a=0.05, b=1) {
  # Bayesian lower bound of population prevalence gamma under a uniform prior
  #
  # Args:
  #  p : density the lower bound should bound (e.g. 0.95)
  #  k : number of participants significant out of 
  #  n : total number of participants
  #  a : alpha value of within-participant test (default=0.05)
  #  b : sensitivity/beta of within-participant test (default=1)
  
  m1 <-  k + 1
  m2 <- n - k + 1
  th_c <- qbeta( p*pbeta(a, m1, m2) + (1-p)*pbeta(b, m1, m2), m1, m2 )
  g_c <- (th_c -a)/(b-a)
  return(g_c)
}


bayesprev_hpdi <- function(p, k, n, a=0.05, b=1) {
  # Bayesian highest posterior density interval of population prevalence gamma
  # under a uniform prior
  #
  # Args:
  #  p : HPDI to return (e.g. 0.95 for 95%)
  #  k : number of participants significant out of 
  #  n : total number of participants
  #  a : alpha value of within-participant test (default=0.05)
  #  b : sensitivity/beta of within-participant test (default=1)
  
  m1 <- k+1
  m2 <- n-k+1
  
  if(m1 ==1) {
    endpts <- c(a, qbeta( (1 -p)*pbeta(a, m1, m2) + p*pbeta(b, m1, m2), m1, m2 ) ) 
    return((endpts -a)/(b-a))
  }
  
  if(m2 ==1) {
    endpts <- c( qbeta( p*pbeta(a, m1, m2) + (1- p)* pbeta(b, m1, m2), m1, m2 ) , b) 
    return( (endpts-a)/(b-a))
  }
  
  if(k<= n*a) {
    endpts <- c(a, qbeta( (1 -p)*pbeta(a, m1, m2) + p*pbeta(b, m1, m2), m1, m2 ) ) 
    return((endpts -a)/(b-a))
  }
  
  if(k>= n*b) {
    endpts <- c( qbeta( p*pbeta(a, m1, m2) + (1- p)* pbeta(b, m1, m2), m1, m2 ) , b) 
    return( (endpts-a)/(b-a))
  }
  
  
  g <- function(x, m1, m2, a, b, p ) {
    y <- numeric(2)
    y[1] <-  pbeta(x[2], m1, m2) - pbeta(x[1], m1, m2) - p*(pbeta(b, m1, m2) - pbeta(a, m1, m2))
    y[2] <- log(dbeta(x[2], m1, m2)) - log(dbeta(x[1], m1, m2))
    return(y)
  }
  
  x_init <- numeric(2)
  
  p1 <- (1-p)/2 
  p2 <- (1 +p)/2
  
  x_init[1] <- qbeta( (1 -p1)*pbeta(a, m1, m2) + p1* pbeta(b, m1, m2), m1, m2 )
  x_init[2] <- qbeta( (1 -p2)*pbeta(a, m1, m2) + p2* pbeta(b, m1, m2), m1, m2 )
  
  opt <- nleqslv(x_init, g, method ="Newton", control=list(maxit=1000), m1=m1, m2=m2, a=a, b=b, p=p)
  
  if (opt$termcd ==1)  print("convergence achieved") 
  if (opt$termcd != 1)  print("failed to converge") 
  
  temp <- opt$x
  if (temp[1] <a) {
    temp[1] <- a
    temp[2] <- qbeta( (1 -p)*pbeta(a, m1, m2) + p* pbeta(b, m1, m2), m1, m2 )
  }
  if (temp[2] > b) {
    temp[1] <- qbeta( p*pbeta(a, m1, m2) + (1-p)* pbeta(b, m1, m2), m1, m2 )
    temp[2] <- b
  }
  endpts <- (temp -a)/(b-a)
  return(endpts) 
}


bayesprev_posteriorprob <- function(x, k, n, a=0.05, b=1) {
  # Bayesian posterior probability in favour of the population prevalence gamma being greater than x
  #
  # Args:
  # x : values of gamma at which to evaluate the posterior probability
  # k : number of participants significant out of 
  # n : total number of participants
  # a : alpha value of within-participant test (default=0.05)
  # b : sensitivity/beta of within-participant test (default=1)
  
  theta <- a + (b-a)*x
  m1 <-  k + 1
  m2 <- n - k + 1
  p <- (pbeta(b,m1,m2)-pbeta(theta,m1,m2)) / (pbeta(b,m1,m2)-pbeta(a,m1,m2))
  return(p)
  
}

bayesprev_posteriorlogodds <- function(x, k, n, a=0.05, b=1) {
  # Bayesian posterior log-odds in favour of the population prevalence gamma being greater than x
  #
  # Args:
  # x : log-odds threshold
  # k : number of participants significant out of 
  # n : total number of participants
  # a : alpha value of within-participant test (default=0.05)
  # b : sensitivity/beta of within-participant test (default=1)
  
  theta <- a + (b-a)*x
  m1 <-  k + 1
  m2 <- n - k + 1
  p <- (pbeta(b,m1,m2)-pbeta(theta,m1,m2)) / (pbeta(b,m1,m2)-pbeta(a,m1,m2))
  lo<- log(p / (1 - p))
  return(lo)
  
}

Non-Moral Vignettes

k <- 69  # number of persons matching hypothesized / group-level pattern
n <- 69+97 # total N
x <- 0.50 # minority cut-off of population prevalence gamma for Bayesian posterior probability
x_0 <- 0.00 # global null of population prevalence gamma for Bayesian posterior probability
map = bayesprev_map(k, n) # get maximum a posteriori (MAP) estimate (i.e., the likeliest person-level prevalence population value)
int = bayesprev_hpdi(0.96, k, n) # get 96% HPDIs

[1] "convergence achieved"

i1 = int[1] # get lower 96% HPDI
i2 = int[2] # get upper 96% HPDI
print(c(i1, map ,i2)) # print lower 96 HPDI, map estimate, and upper 96 HPDI

[1] 0.3045692 0.3849081 0.4682054

bayesprev_posteriorprob(x, k, n) # Bayesian posterior probability in favor of population prevalence being greater than x

[1] 0.002422369

bayesprev_posteriorprob(x_0, k, n) # Bayesian posterior probability in favor of population prevalence being greater than x_0

[1] 1

# ORIGINALLY PRINTED: 0.3045692 0.3849081 0.4682054

Moral Vignettes

k <- 129  # number of persons matching hypothesized / group-level pattern
n <- 129+37 # total N
x <- 0.50 # minority cut-off of population prevalence gamma for Bayesian posterior probability
x_0 <- 0.00 # global null of population prevalence gamma for Bayesian posterior probability
map = bayesprev_map(k, n) # get maximum a posteriori (MAP) estimate (i.e., the likeliest person-level prevalence population value)
int = bayesprev_hpdi(0.96, k, n) # get 96% HPDIs

[1] "convergence achieved"

i1 = int[1] # get lower 96% HPDI
i2 = int[2] # get upper 96% HPDI
print(c(i1, map ,i2)) # print lower 96 HPDI, map estimate, and upper 96 HPDI

[1] 0.6912801 0.7653773 0.8297420

bayesprev_posteriorprob(x, k, n) # Bayesian posterior probability in favor of population prevalence being greater than x

[1] 1

bayesprev_posteriorprob(x_0, k, n) # Bayesian posterior probability in favor of population prevalence being greater than x_0

[1] 1

# ORIGINALLY PRINTED: 0.6912801 0.7653773 0.8297420

Interaction

k <- 51  # number of persons matching hypothesized / group-level pattern
n <- 51+115 # total N
x <- 0.50 # minority cut-off of population prevalence gamma for Bayesian posterior probability
x_0 <- 0.00 # global null of population prevalence gamma for Bayesian posterior probability
map = bayesprev_map(k, n) # get maximum a posteriori (MAP) estimate (i.e., the likeliest person-level prevalence population value)
int = bayesprev_hpdi(0.96, k, n) # get 96% HPDIs

[1] "convergence achieved"

i1 = int[1] # get lower 96% HPDI
i2 = int[2] # get upper 96% HPDI
print(c(i1, map ,i2)) # print lower 96 HPDI, map estimate, and upper 96 HPDI

[1] 0.1975061 0.2707673 0.3507936

bayesprev_posteriorprob(x, k, n) # Bayesian posterior probability in favor of population prevalence being greater than x

[1] 0.000000007825034

bayesprev_posteriorprob(x_0, k, n) # Bayesian posterior probability in favor of population prevalence being greater than x_0

[1] 1

# ORIGINALLY PRINTED: 0.1975061 0.2707673 0.3507936

Primary Plots

# set breaks for histograms
breaks_seezero <- c(-4, -3.5, -3, -2.5, -2, -1.5, -1, -0.5, -0.01, 0.01, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4)

Non-Moral Vignettes

print(S1a_nonmoral <- ggplot(Study_1a_plvl, aes(x = NonmoralAverage)) +
  geom_histogram(breaks = breaks_seezero,
                 color = "grey30", fill = c("#CEA335", "#CEA335", "#CEA335", "#CEA335",
                                           "#CEA335", "#CEA335", "#CEA335", "#CEA335",
                                           "black",
                                           "#CEA335", "#CEA335", "#CEA335", "#CEA335",
                                           "#CEA335", "#CEA335", "#CEA335", "#CEA335")) +
  scale_x_continuous(limits = c(-4.1,4.1), breaks = c(-4,-3,-2,-1,0,1,2,3,4)) +
  scale_y_continuous(limits = c(-0.1,50.1), breaks = c(0, 10, 20, 30, 40, 50)) +
  xlab("\nGlobal Virtuosity Judgment") +
  ylab("Number of Participants\n") +
  theme_classic() +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 16),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

ggsave("S1a_nonmoral_plot.png")

Saving 14.1 x 9.04 in image

Moral Vignettes

print(S1a_moral <- ggplot(Study_1a_plvl, aes(x = MoralAverage)) +
  geom_histogram(breaks = breaks_seezero,
                 color = "grey30", fill = c("#CEA335", "#CEA335", "#CEA335", "#CEA335",
                                           "#CEA335", "#CEA335", "#CEA335", "#CEA335",
                                           "black",
                                           "#CEA335", "#CEA335", "#CEA335", "#CEA335",
                                           "#CEA335", "#CEA335", "#CEA335", "#CEA335")) +
  scale_x_continuous(limits = c(-4.1,4.1), breaks = c(-4,-3,-2,-1,0,1,2,3,4)) +
  scale_y_continuous(limits = c(-0.1,50.1), breaks = c(0, 10, 20, 30, 40, 50)) +
  xlab("\nGlobal Virtuosity Judgment") +
  ylab("Number of Participants\n") +
  theme_classic() +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 16),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

ggsave("S1a_moral_plot.png")

Saving 14.1 x 9.04 in image

Interaction

Berman & Small predicted an interaction such that tempted agents would be judged as more virtuous in non-moral contexts, but less virtuous in moral contexts. Therefore, this interaction plot shows how many participants’ responses mirrored that nuanced prediction.

print(S1a_int <- ggplot(Study_1a_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#CEA335") +
  scale_y_continuous(limits = c(-0.1,151.1), breaks = c(0, 25, 50, 75, 100, 125, 150)) +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 16),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

    
ggsave("S1a_int_plot.png")

Saving 14.1 x 9.04 in image

Robustness Plots

Here, we investigate the number of participants matching predicted effects by demographics in order to catalog the generality (or lack thereof) of these effects across levels of various the demographics that we collected.

# create yes/no variable for simple hypotheses
Study_1a_plvl <- Study_1a_plvl %>%
  mutate(fits_nonmoral = if_else(fits_hypothesis_n == 1, "Yes", "No"))

# create yes/no variable for simple hypotheses
Study_1a_plvl <- Study_1a_plvl %>%
  mutate(fits_moral = if_else(fits_hypothesis_m == 1, "Yes", "No"))

Education

Non-Moral Vignettes

print(S1a_nonmoral_educ <- ggplot(Study_1a_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Education) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Moral Vignettes

print(S1a_moral_educ <- ggplot(Study_1a_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Education) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Interaction

print(S1a_int_educ <- ggplot(Study_1a_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Education) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Political Party

Non-Moral Vignettes

print(S1a_nonmoral_pol <- ggplot(Study_1a_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Political_Party) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Moral Vignettes

print(S1a_moral_pol <- ggplot(Study_1a_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Political_Party) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Interaction

print(S1a_int_pol <- ggplot(Study_1a_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Political_Party) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Religiosity

Non-Moral Vignettes

print(S1a_nonmoral_relig <- ggplot(Study_1a_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Religiosity) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Moral Vignettes

print(S1a_moral_relig <- ggplot(Study_1a_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Religiosity) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Interaction

print(S1a_int_relig <- ggplot(Study_1a_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Religiosity) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Urban vs Rural

Non-Moral Vignettes

print(S1a_nonmoral_urb <- ggplot(Study_1a_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~UrbanRural) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Moral Vignettes

print(S1a_moral_urb <- ggplot(Study_1a_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~UrbanRural) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Interaction

print(S1a_int_urb <- ggplot(Study_1a_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~UrbanRural) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Gender

Non-Moral Vignettes

print(S1a_nonmoral_gender <- ggplot(Study_1a_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Gender) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Moral Vignettes

print(S1a_moral_gender <- ggplot(Study_1a_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Gender) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Interaction

print(S1a_int_gender <- ggplot(Study_1a_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Gender) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Race

# make race one variable
Study_1a_plvl <- Study_1a_plvl %>%
  mutate(Race = coalesce(Race_1, Race_2, Race_3, Race_4, Race_5, Race_6, Race_7, Race_8))

Non-Moral Vignettes

print(S1a_nonmoral_race <- ggplot(Study_1a_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Race) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Moral Vignettes

print(S1a_moral_race <- ggplot(Study_1a_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Race) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))

Interaction

print(S1a_int_race <- ggplot(Study_1a_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#CEA335") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Race) +
  theme(axis.title.x = element_text(size = 18), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_text(color = "black", size = 16), 
              axis.text.y = element_text(color = "black", size = 16),
              strip.text.x = element_text(color = "black", size = 12),
              legend.position = "right",
              legend.title = element_text(color = "black", size = 18),
              legend.text = element_text(color = "black", size = 16)))