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

Data Preparation

In this pre-registered study, we collected data from 176 participants. Two of 176 IP addresses were identical, which prompted their removal from the dataset. Additionally, six participants failed the attention checks, leading to an analyzable sample of 168 participants.

Loading Packages

Before running this chunk, please load “Study 1b.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_1b$AttChecks <- rowSums(Study_1b[, c("AttnCheck1_S2", "AttnCheck2_S2", "AttnCheck3_S2")], na.rm = T)

# The exclusion of any sum (within the variable "AttChecks") that did not result in a zero
Study_1b <- Study_1b %>% 
  filter(Study_1b$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_1b$MoralAverage <- rowMeans(Study_1b[, c("AdulteryVirtuous1_S2", "AdulteryHonorable2_S2", "AdulteryRespectable3_S2", "OfficeVirtuous1_S2", "OfficeHonorable2_S2", "OfficeRespectable3_S2", "KidneyVirtuous1_S2", "KidneyHonorable2_S2", "KidneyRespectable3_S2", "AdulteryVirtuous2_S2", "AdulteryHonrable3_S2", "AdulteryRespectable1_S2", "OfficeVirtuous2_S2", "OfficeHonorable3_S2", "OfficeRespectable1_S2", "KidneyVirtuous2_S2", "KidneyHonorable3_S2", "KidneyRespectable1_S2", "AdulteryVirtuous3_S2", "AdulteryHonrable1_S2", "AdulteryRespectable2_S2", "OfficeVirtuous3_S2", "OfficeHonorable1_S2", "OfficeRespectable2_S2", "KidneyVirtuous3_S2", "KidneyHonorable1_S2", "KidneyRespectable2_S2")], na.rm = T)

# Create a variable for virtuosity average in the non-moral condition
Study_1b$NonmoralAverage <- rowMeans(Study_1b[, c("DessertVirtuous1_S2", "DessertHonorable2_S2", "DessertRespectable3_S2", "GymVirtuous1_S2", "GymHonorable2_S2", "GymRespectable3_S2", "SchoolworkVirtuous1_S2", "SchoolworkHonorable2_S2", "SchoolworkRespectable3_S2", "DessertVirtuous2_S2", "DessertHonorable3_S2", "DessertRespectable1_S2", "GymVirtuous2_S2", "GymHonorable3_S2", "GymRespectable1_S2", "SchoolworkVirtuous2_S2", "SchoolworkHonorable3_S2", "SchoolworkRespectable1_S2", "DessertVirtuous3_S2", "DessertHonorable1_S2", "DessertRespectable2_S2", "GymVirtuous3_S2", "GymHonorable1_S2", "GymRespectable2_S2", "SchoolworkVirtuous3_S2", "SchoolworkHonorable1_S2", "SchoolworkRespectable2_S2")], na.rm = T)

# Create a variable for virtuosity average in the Office Supplies scenario
Study_1b$OfficeAverage <- rowMeans(Study_1b[, c("OfficeVirtuous1_S2", "OfficeHonorable2_S2", "OfficeRespectable3_S2", "OfficeVirtuous2_S2", "OfficeHonorable3_S2", "OfficeRespectable1_S2", "OfficeVirtuous3_S2", "OfficeHonorable1_S2", "OfficeRespectable2_S2" )], na.rm = T)

# Create a variable for virtuosity average in the Adultery scenario
Study_1b$AdulteryAverage <- rowMeans(Study_1b[, c("AdulteryVirtuous1_S2", "AdulteryHonorable2_S2", "AdulteryRespectable3_S2", "AdulteryVirtuous2_S2", "AdulteryHonrable3_S2", "AdulteryRespectable1_S2", "AdulteryVirtuous3_S2", "AdulteryHonrable1_S2", "AdulteryRespectable2_S2" )], na.rm = T)

# Create a variable for virtuosity average in the Kidney scenario
Study_1b$KidneyAverage <- rowMeans(Study_1b[, c("KidneyVirtuous1_S2", "KidneyHonorable2_S2", "KidneyRespectable3_S2", "KidneyVirtuous2_S2", "KidneyHonorable3_S2", "KidneyRespectable1_S2", "KidneyVirtuous3_S2", "KidneyHonorable1_S2", "KidneyRespectable2_S2")], na.rm = T)

# Create a variable for virtuosity average in the Dessert scenario
Study_1b$DessertAverage <- rowMeans(Study_1b[, c("DessertVirtuous1_S2", "DessertHonorable2_S2", "DessertRespectable3_S2", "DessertVirtuous2_S2", "DessertHonorable3_S2", "DessertRespectable1_S2", "DessertVirtuous3_S2", "DessertHonorable1_S2", "DessertRespectable2_S2")], na.rm = T)

# Create a variable for virtuosity average in the Gym scenario
Study_1b$GymAverage <- rowMeans(Study_1b[, c("GymVirtuous1_S2", "GymHonorable2_S2", "GymRespectable3_S2", "GymVirtuous2_S2", "GymHonorable3_S2", "GymRespectable1_S2", "GymVirtuous3_S2", "GymHonorable1_S2", "GymRespectable2_S2")], na.rm = T)

# Create a variable for virtuosity average in the Schoolwork scenario
Study_1b$SchoolworkAverage <- rowMeans(Study_1b[, c("SchoolworkVirtuous1_S2", "SchoolworkHonorable2_S2", "SchoolworkRespectable3_S2", "SchoolworkVirtuous2_S2", "SchoolworkHonorable3_S2", "SchoolworkRespectable1_S2","SchoolworkVirtuous3_S2", "SchoolworkHonorable1_S2", "SchoolworkRespectable2_S2")], 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_1b_nonmoral_t <- t.test(Study_1b$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_1b_nonmoral_t)

    One Sample t-test

data:  Study_1b$NonmoralAverage
t = -2.0956, df = 167, p-value = 0.01881
alternative hypothesis: true mean is less than 0
95 percent confidence interval:
       -Inf -0.0303816
sample estimates:
 mean of x 
-0.1441799 

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

    One Sample t-test

data:  Study_1b$NonmoralAverage
t = -2.0956, df = 167, p-value = 0.03762
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.28001104 -0.00834875
sample estimates:
 mean of x 
-0.1441799 

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

Cohen's d (single sample)

d estimate: -0.1616804 (negligible)
Reference mu: 0
95 percent confidence interval:
     lower      upper 
-0.4668144  0.1434535 

Dessert Vignette

# Dessert vignette t-test
Study_1b_dessert_t <- t.test(Study_1b$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_1b_dessert_t)

    One Sample t-test

data:  Study_1b$DessertAverage
t = -2.6352, df = 167, p-value = 0.004601
alternative hypothesis: true mean is less than 0
95 percent confidence interval:
       -Inf -0.0679641
sample estimates:
 mean of x 
-0.1825397 

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

    One Sample t-test

data:  Study_1b$DessertAverage
t = -2.6352, df = 167, p-value = 0.009201
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.31929861 -0.04578076
sample estimates:
 mean of x 
-0.1825397 

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

Cohen's d (single sample)

d estimate: -0.2033077 (small)
Reference mu: 0
95 percent confidence interval:
     lower      upper 
-0.5087303  0.1021149 

Gym Vignette

# Gym vignette t-test
Study_1b_gym_t <- t.test(Study_1b$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_1b_gym_t)

    One Sample t-test

data:  Study_1b$GymAverage
t = -1.9341, df = 167, p-value = 0.0274
alternative hypothesis: true mean is less than 0
95 percent confidence interval:
        -Inf -0.02441968
sample estimates:
 mean of x 
-0.1686508 

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

    One Sample t-test

data:  Study_1b$GymAverage
t = -1.9341, df = 167, p-value = 0.05479
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.340806954  0.003505366
sample estimates:
 mean of x 
-0.1686508 

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

Cohen's d (single sample)

d estimate: -0.1492169 (negligible)
Reference mu: 0
95 percent confidence interval:
     lower      upper 
-0.4542771  0.1558434 

Schoolwork Vignette

# Schoolwork vignette t-test
Study_1b_schoolwork_t <- t.test(Study_1b$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_1b_schoolwork_t)

    One Sample t-test

data:  Study_1b$SchoolworkAverage
t = -0.94962, df = 167, p-value = 0.1718
alternative hypothesis: true mean is less than 0
95 percent confidence interval:
       -Inf 0.06034302
sample estimates:
  mean of x 
-0.08134921 

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

    One Sample t-test

data:  Study_1b$SchoolworkAverage
t = -0.94962, df = 167, p-value = 0.3437
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.25047492  0.08777651
sample estimates:
  mean of x 
-0.08134921 

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

Cohen's d (single sample)

d estimate: -0.07326487 (negligible)
Reference mu: 0
95 percent confidence interval:
     lower      upper 
-0.3780037  0.2314739 

Moral Vignettes

Overall

# Moral t-test using person-level averages across Moral vignettes
Study_1b_moral_t <- t.test(Study_1b$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_1b_moral_t)

    One Sample t-test

data:  Study_1b$MoralAverage
t = 10.334, df = 167, p-value < 0.00000000000000022
alternative hypothesis: true mean is greater than 0
95 percent confidence interval:
 0.8049388       Inf
sample estimates:
mean of x 
0.9583333 

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

    One Sample t-test

data:  Study_1b$MoralAverage
t = 10.334, df = 167, p-value < 0.00000000000000022
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.7752396 1.1414270
sample estimates:
mean of x 
0.9583333 

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

Cohen's d (single sample)

d estimate: 0.7972515 (medium)
Reference mu: 0
95 percent confidence interval:
    lower     upper 
0.4807443 1.1137587 

Office Vignette

# Office vignette t-test
Study_1b_office_t <- t.test(Study_1b$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_1b_office_t)

    One Sample t-test

data:  Study_1b$OfficeAverage
t = 7.0028, df = 167, p-value = 0.00000000002924
alternative hypothesis: true mean is greater than 0
95 percent confidence interval:
 0.592554      Inf
sample estimates:
mean of x 
0.7757937 

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

    One Sample t-test

data:  Study_1b$OfficeAverage
t = 7.0028, df = 167, p-value = 0.00000000005848
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.5570764 0.9945109
sample estimates:
mean of x 
0.7757937 

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

Cohen's d (single sample)

d estimate: 0.5402756 (medium)
Reference mu: 0
95 percent confidence interval:
    lower     upper 
0.2301311 0.8504201 

Adultery Vignette

# Adultery vignette t-test
Study_1b_adultery_t <- t.test(Study_1b$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_1b_adultery_t)

    One Sample t-test

data:  Study_1b$AdulteryAverage
t = 10.843, df = 167, p-value < 0.00000000000000022
alternative hypothesis: true mean is greater than 0
95 percent confidence interval:
 1.229141      Inf
sample estimates:
mean of x 
 1.450397 

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

    One Sample t-test

data:  Study_1b$AdulteryAverage
t = 10.843, df = 167, p-value < 0.00000000000000022
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 1.186302 1.714491
sample estimates:
mean of x 
 1.450397 

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

Cohen's d (single sample)

d estimate: 0.8365266 (large)
Reference mu: 0
95 percent confidence interval:
    lower     upper 
0.5188457 1.1542076 

Kidney Vignette

# Kidney vignette t-test
Study_1b_kidney_t <- t.test(Study_1b$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_1b_kidney_t)

    One Sample t-test

data:  Study_1b$KidneyAverage
t = 6.5382, df = 167, p-value = 0.000000000363
alternative hypothesis: true mean is greater than 0
95 percent confidence interval:
 0.4846732       Inf
sample estimates:
mean of x 
0.6488095 

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

    One Sample t-test

data:  Study_1b$KidneyAverage
t = 6.5382, df = 167, p-value = 0.000000000726
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.4528942 0.8447249
sample estimates:
mean of x 
0.6488095 

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

Cohen's d (single sample)

d estimate: 0.50443 (medium)
Reference mu: 0
95 percent confidence interval:
    lower     upper 
0.1949866 0.8138734 

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_1b_anova <- Study_1b %>%
  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_1b_interaction_aov <- aov(VirtuosityAverage~Condition + Error(as.factor(ResponseId)), data = Study_1b_anova)
summary(Study_1b_interaction_aov)

Error: as.factor(ResponseId)
           Df Sum Sq Mean Sq F value Pr(>F)
Residuals 167  228.8    1.37               

Error: Within
           Df Sum Sq Mean Sq F value              Pr(>F)    
Condition   1  102.1  102.10   117.3 <0.0000000000000002 ***
Residuals 167  145.3    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_1b_interaction_t <- t.test(VirtuosityAverage ~ Condition, data = Study_1b_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_1b_interaction_t

    Paired t-test

data:  VirtuosityAverage by Condition
t = 10.833, df = 167, p-value < 0.00000000000000022
alternative hypothesis: true mean difference is greater than 0
95 percent confidence interval:
 0.9341718       Inf
sample estimates:
mean difference 
       1.102513 

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

    Paired t-test

data:  VirtuosityAverage by Condition
t = 10.833, df = 167, p-value < 0.00000000000000022
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.9015787 1.3034477
sample estimates:
mean difference 
       1.102513 

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

Cohen's d

d estimate: 0.8357596 (large)
95 percent confidence interval:
    lower     upper 
0.6594744 1.0120448 

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

Cohen's d

d estimate: 1.035001 (large)
95 percent confidence interval:
    lower     upper 
0.8021003 1.2679011 

Person-Lvl Analyses

It is first necessary to rearrange data in a way that facilitates conducting various person-level analyses.

# Create person-level variables
Study_1b_plvl <- Study_1b %>% 
  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_1b_plvl$number_of_hyp_fit <- rowSums(Study_1b_plvl[, c("fits_hypothesis_m", "fits_hypothesis_n")], na.rm = T)
Study_1b_plvl <- Study_1b_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_1b_plvl$fits_hypothesis_n) # (0 = "no"; 1 = "yes")

  0   1 
108  60 

Moral Vignettes

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

  0   1 
 43 125 

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_1b_plvl$fits_complex) # (0 = "no"; 1 = "yes")

 No Yes 
124  44 

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

Registered S3 method overwritten by 'mosaic':
  method                           from   
  fortify.SpatialPolygonsDataFrame ggplot2

The 'mosaic' package masks several functions from core packages in order to add 
additional features.  The original behavior of these functions should not be affected by this.

Attaching package: ‘mosaic’

The following object is masked from ‘package:Matrix’:

    mean

The following objects are masked from ‘package:dplyr’:

    count, do, tally

The following object is masked from ‘package:purrr’:

    cross

The following object is masked from ‘package:ggplot2’:

    stat

The following object is masked from ‘package:correlation’:

    cor_test

The following objects are masked from ‘package:sjstats’:

    prop, props

The following objects are masked from ‘package:psych’:

    logit, rescale

The following objects are masked from ‘package:stats’:

    binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test, quantile, sd, t.test, var

The following objects are masked from ‘package:base’:

    max, mean, min, prod, range, sample, sum

Study_1b_random <- Study_1b_plvl %>% 
  dplyr::select(c(ResponseId, MoralAverage, NonmoralAverage))

Study_1b_random <- Study_1b_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_1b <- prop( ~ Interaction == "Pos, Neg, Pos", data = Study_1b_random) # create descriptive proportion of claimed pattern from Berman & Small
prop_Study_1b_dp <- unname(prop_Study_1b[1]) # get descriptive proportion proportion for interaction

nsim <- 1000 # number of random datasets to generate
nullprop_Study_1b <- 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_1b_random$NonmoralAverageR <- shuffle(Study_1b_random$NonmoralAverage) # shuffling 'Non-Moral' values  
  Study_1b_random$MoralAverageR <- shuffle(Study_1b_random$MoralAverage) # shuffling 'Moral' values
  
  ## Creating variable that distinguishes Ps who showed predicted (claimed) pattern from Ps who did not
  Study_1b_R <- Study_1b_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_1b <-prop( ~ InteractionR == "Pos, Neg, Pos", data = Study_1b_R) # create proportion value of claimed direction 
  nullprop_Study_1b[i]<- unname(nullprop_TF_Study_1b[1]) # get claimed pattern proportion for randomly shuffled data and repeat for each nsim dataset
}

gf_histogram( ~ nullprop_Study_1b, fill = ~ (nullprop_Study_1b >= prop_Study_1b_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_1b >= prop_Study_1b_dp) # suggests that empirical interaction proportion is very likely to have occurred by chance

prop_TRUE 
    0.672 

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 <- 60  # number of persons matching hypothesized / group-level pattern
n <- 60+108 # 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.9999957

Moral Vignettes

k <- 125  # number of persons matching hypothesized / group-level pattern
n <- 125+43 # 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.000000004627076

Interaction

k <- 44  # number of persons matching hypothesized / group-level pattern
n <- 44+124 # 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

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 <- 60  # number of persons matching hypothesized / group-level pattern
n <- 60+108 # 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.2466803 0.3233083 0.4048897

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

[1] 0.000006394378

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

[1] 1

Moral Vignettes

k <- 125  # number of persons matching hypothesized / group-level pattern
n <- 125+43 # 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.6542208 0.7305764 0.7984644

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

Interaction

k <- 44  # number of persons matching hypothesized / group-level pattern
n <- 44+124 # 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.1545315 0.2230576 0.2998444

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

[1] 0.000000000002327472

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

[1] 1

Primary Plots

# set breaks for histograms / groups for histogram values
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(S1b_nonmoral <- ggplot(Study_1b_plvl, aes(x = NonmoralAverage)) +
  geom_histogram(breaks = breaks_seezero,
                 color = "grey30", fill = c("#8b1616", "#8b1616", "#8b1616", "#8b1616",
                                           "#8b1616", "#8b1616", "#8b1616", "#8b1616",
                                           "black",
                                           "#8b1616", "#8b1616", "#8b1616", "#8b1616",
                                           "#8b1616", "#8b1616", "#8b1616", "#8b1616")) +
  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,75.1), breaks = c(0, 15, 30, 45, 60, 75)) +
  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("S1b_nonmoral_plot.png")

Saving 14.1 x 9.04 in image

Moral Vignettes

print(S1b_moral <- ggplot(Study_1b_plvl, aes(x = MoralAverage)) +
  geom_histogram(breaks = breaks_seezero,
                 color = "grey30", fill = c("#8b1616", "#8b1616", "#8b1616", "#8b1616",
                                           "#8b1616", "#8b1616", "#8b1616", "#8b1616",
                                           "black",
                                           "#8b1616", "#8b1616", "#8b1616", "#8b1616",
                                           "#8b1616", "#8b1616", "#8b1616", "#8b1616")) +
  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("S1b_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(S1b_int <- ggplot(Study_1b_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#8b1616") +
  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("S1b_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_1b_plvl <- Study_1b_plvl %>%
  mutate(fits_nonmoral = if_else(fits_hypothesis_n == 1, "Yes", "No"))

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

Education

Non-Moral Vignettes

print(S1b_nonmoral_educ <- ggplot(Study_1b_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Education_S2) +
  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(S1b_moral_educ <- ggplot(Study_1b_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Education_S2) +
  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(S1b_int_educ <- ggplot(Study_1b_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Education_S2) +
  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(S1b_nonmoral_pol <- ggplot(Study_1b_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Political_Party_S2) +
  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(S1b_moral_pol <- ggplot(Study_1b_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Political_Party_S2) +
  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(S1b_int_pol <- ggplot(Study_1b_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Political_Party_S2) +
  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(S1b_nonmoral_relig <- ggplot(Study_1b_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Religiosity_S2) +
  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(S1b_moral_relig <- ggplot(Study_1b_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Religiosity_S2) +
  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(S1b_int_relig <- ggplot(Study_1b_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Religiosity_S2) +
  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(S1b_nonmoral_urb <- ggplot(Study_1b_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~UrbanRural_S2) +
  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(S1b_moral_urb <- ggplot(Study_1b_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~UrbanRural_S2) +
  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(S1b_int_urb <- ggplot(Study_1b_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~UrbanRural_S2) +
  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(S1b_nonmoral_gender <- ggplot(Study_1b_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Non-Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Gender_S2) +
  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(S1b_moral_gender <- ggplot(Study_1b_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Moral Vignette Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Gender_S2) +
  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(S1b_int_gender <- ggplot(Study_1b_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#8b1616") +
  xlab("\nMatched Predicted Interaction Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Gender_S2) +
  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_1b_plvl <- Study_1b_plvl %>%
  mutate(Race = coalesce(Race_S2_1, Race_S2_2, Race_S2_3, Race_S2_4, Race_S2_5, Race_S2_6, Race_S2_7, Race_S2_8))

Non-Moral Vignettes

print(S1b_nonmoral_race <- ggplot(Study_1b_plvl, aes(x = fits_nonmoral)) +
  geom_bar(fill = "#8b1616") +
  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(S1b_moral_race <- ggplot(Study_1b_plvl, aes(x = fits_moral)) +
  geom_bar(fill = "#8b1616") +
  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(S1b_int_race <- ggplot(Study_1b_plvl, aes(x = fits_complex)) +
  geom_bar(fill = "#8b1616") +
  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)))