Study 2: Extension w/ Participant-Generated Stimuli & BS/SB Measures

Data Preparation

In this pre-registered study, we collected data from 361 participants. After filtering these participants with an attention check measure, we obtained analyzable data from 350 participants.

Loading Data/Packages

Before running this chunk, please load “Study 2.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 in the continuous virtuosity condition were instructed to select zero (i.e., the middle option in the 9-point scale) as their response for the virtuosity measure. Additionally, participants in the binary goodness condition, selected “Monica” as their response for the “more good” measure. If a participant failed to select either of these responses, they would be excluded from the analyses.

S2_bin <- Study_2 %>% 
  filter(AttnBinQ == "Monica") %>% 
  select(-c(AttnContV, s4_V1Theft, s4_V1Cheat, s4_V1Drugs, s4_V1Sex,
            s4_V1Sexualizing, s4_V1Health, s4_V1Blame, s4_V1Revenge,
            s4_V1Workethic, s4_V1Criticism))

S2_cont <- Study_2 %>% 
  filter(AttnContV == "0\n") %>% 
  select(-c(AttnBinQ, s4_sbTheft, s4_sbCheat, s4_sbDrugs, s4_sbSex,
            s4_sbSexualization, s4_sbHealth, s4_sbBlame, s4_sbRevenge,
            s4_sbWorkethic, s4_sbCriticism))

Variable Creation

Here, we create the relevant variables necessary to begin analyses for both the binary goodness measure and the continuous virtuosity measure.

# Create a variable per scenario that indicates if each participant chose the hypothesized agent (i.e., 1), or not (i.e., 0)
S2_bin <- S2_bin %>%  
  mutate(theft_hyp = if_else(s4_sbTheft == "Kenna", 1, 0)) %>% 
  mutate(drugs_hyp = if_else(s4_sbDrugs == "Carl", 1, 0)) %>% 
  mutate(cheat_hyp = if_else(s4_sbCheat == "James", 1, 0)) %>% 
  mutate(sex_hyp = if_else(s4_sbSex == "Camila", 1, 0)) %>% 
  mutate(sexualize_hyp = if_else(s4_sbSexualization == "Robert", 1, 0)) %>% 
  mutate(revenge_hyp = if_else(s4_sbRevenge == "Clara", 1, 0)) %>% 
  mutate(workethic_hyp = if_else(s4_sbWorkethic == "Vivian", 1, 0)) %>% 
  mutate(blame_hyp = if_else(s4_sbBlame == "Eric", 1, 0)) %>% 
  mutate(criticism_hyp = if_else(s4_sbCriticism == "Shawn", 1, 0)) %>% 
  mutate(health_hyp = if_else(s4_sbHealth == "Summer", 1, 0))

# Create a variable that averages the binary choice values within participants where the hypothesized non-tempted agent has a value of 1 and the tempted agent has a value of 0
S2_bin$BinaryAverage <- rowMeans(S2_bin[, (c("theft_hyp", "drugs_hyp", "cheat_hyp", "sex_hyp", "sexualize_hyp","revenge_hyp","workethic_hyp", "blame_hyp", "criticism_hyp","health_hyp"))], na.rm = T)

# Create a variable per scenario that indicates if each participant chose the hypothesized agent or not
S2_cont <- S2_cont %>%
  mutate_at(vars(s4_V1Theft, s4_V1Drugs, s4_V1Cheat, s4_V1Sex, s4_V1Sexualizing,                                                                   s4_V1Revenge, s4_V1Workethic, s4_V1Blame, s4_V1Criticism, s4_V1Health), as.numeric)

S2_cont <- S2_cont %>%  
  mutate(theft_hyp = if_else(s4_V1Theft > 0, 1, 0)) %>%
  mutate(drugs_hyp = if_else(s4_V1Drugs > 0, 1, 0)) %>%
  mutate(cheat_hyp = if_else(s4_V1Cheat > 0, 1, 0)) %>%
  mutate(sex_hyp = if_else(s4_V1Sex > 0, 1, 0)) %>%
  mutate(sexualize_hyp = if_else(s4_V1Sexualizing > 0, 1, 0)) %>%
  mutate(revenge_hyp = if_else(s4_V1Revenge > 0, 1, 0)) %>%
  mutate(workethic_hyp = if_else(s4_V1Workethic > 0, 1, 0)) %>%
  mutate(blame_hyp = if_else(s4_V1Blame > 0, 1, 0)) %>%
  mutate(criticism_hyp = if_else(s4_V1Criticism > 0, 1, 0)) %>%
  mutate(health_hyp = if_else(s4_V1Health > 0, 1, 0))

# Create a variable per scenario that indicates if each participant chose the non-hypothesized agent
## this would rule out that any one vignette that does not have a predicted majority pattern is due to an opposite majority pattern
## for instance, it could be the case that a non-majority shows the predicted pattern but that the modal or second-modal pattern is judging both agents as equally virtuous (and not the tempted agent as more virtuous)
S2_cont <- S2_cont %>%  
  mutate(theft_opp = if_else(s4_V1Theft < 0, 1, 0)) %>%
  mutate(drugs_opp = if_else(s4_V1Drugs < 0, 1, 0)) %>%
  mutate(cheat_opp = if_else(s4_V1Cheat < 0, 1, 0)) %>%
  mutate(sex_opp = if_else(s4_V1Sex < 0, 1, 0)) %>%
  mutate(sexualize_opp = if_else(s4_V1Sexualizing < 0, 1, 0)) %>%
  mutate(revenge_opp = if_else(s4_V1Revenge < 0, 1, 0)) %>%
  mutate(workethic_opp = if_else(s4_V1Workethic < 0, 1, 0)) %>%
  mutate(blame_opp = if_else(s4_V1Blame < 0, 1, 0)) %>%
  mutate(criticism_opp = if_else(s4_V1Criticism < 0, 1, 0)) %>%
  mutate(health_opp = if_else(s4_V1Health < 0, 1, 0))


# Create a variable that averages the continuous values within participants where positive numbers indicate choosing of non-tempted agent
S2_cont$ContinuousAverage <- rowMeans(S2_cont[, (c("s4_V1Theft", "s4_V1Drugs", "s4_V1Cheat", "s4_V1Sex", "s4_V1Sexualizing",                                                    "s4_V1Revenge", "s4_V1Workethic", "s4_V1Blame", "s4_V1Criticism", "s4_V1Health"))],
                                      na.rm = T)

Group-Lvl Analyses

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.

Binary Goodness

# t-test using person-level averages across binary measure
S2_bin_t <- t.test(S2_bin$BinaryAverage, mu = 0.5, alternative = "greater") # mu = 0.5 because of binary nature of DV
## 'greater' is for one-tailed p-value because of the non-tempted > tempted hypothesis (numbers > 0.5 mean non-tempted > tempted was chosen more often)
print(S2_bin_t)

    One Sample t-test

data:  S2_bin$BinaryAverage
t = 10.465, df = 176, p-value < 0.00000000000000022
alternative hypothesis: true mean is greater than 0.5
95 percent confidence interval:
 0.6869496       Inf
sample estimates:
mean of x 
0.7220339 

S2_bin_t_2side <- t.test(S2_bin$BinaryAverage, mu = 0.5, alternative = "two.sided") 
## 'two.sided' is to get uncertainty estimates in both directions
print(S2_bin_t_2side)

    One Sample t-test

data:  S2_bin$BinaryAverage
t = 10.465, df = 176, p-value < 0.00000000000000022
alternative hypothesis: true mean is not equal to 0.5
95 percent confidence interval:
 0.6801604 0.7639073
sample estimates:
mean of x 
0.7220339 

# returns cohen's d
effsize::cohen.d(S2_bin$BinaryAverage,
                 f = NA,
                 mu = 0.5,
                 data = data)

Cohen's d (single sample)

d estimate: 0.7865722 (medium)
Reference mu: 0.5
95 percent confidence interval:
    lower     upper 
0.4786338 1.0945106 

Continuous Virtuosity

# t-test using person-level averages across binary measure
S2_cont_t <- t.test(S2_cont$ContinuousAverage, mu = 0, alternative = "greater") # mu = 0 because of continuous nature of DV
## 'greater' is for one-tailed p-value because of the non-tempted > tempted hypothesis (pos numbers mean non-tempted > tempted)
print(S2_cont_t)

    One Sample t-test

data:  S2_cont$ContinuousAverage
t = 9.9387, df = 172, p-value < 0.00000000000000022
alternative hypothesis: true mean is greater than 0
95 percent confidence interval:
 0.7656622       Inf
sample estimates:
mean of x 
0.9184971 

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

    One Sample t-test

data:  S2_cont$ContinuousAverage
t = 9.9387, df = 172, p-value < 0.00000000000000022
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.7360804 1.1009138
sample estimates:
mean of x 
0.9184971 

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

Cohen's d (single sample)

d estimate: 0.7556224 (medium)
Reference mu: 0
95 percent confidence interval:
   lower    upper 
0.444958 1.066287 

Person-Lvl Analyses

# create variables for whether hypothesis fit for each participant
S2_bin <- S2_bin %>%
  mutate(fits_hypothesis = if_else(BinaryAverage > 0.50, 1, 0))

S2_cont <- S2_cont %>%
  mutate(fits_hypothesis = if_else(ContinuousAverage > 0, 1, 0))

Descriptive Pervasiveness

Binary Goodness

# How many participants fit hypothesis for binary measure
table(S2_bin$fits_hypothesis) # 0 = no, 1 = yes

  0   1 
 49 128 

Continuous Virtuosity

# How many participants fit hypothesis for continuous
table(S2_cont$fits_hypothesis) # 0 = no, 1 = yes

  0   1 
 43 130 

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))
}

Binary Goodness

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

[1] 0.00000005534935

Continuous Virtuosity

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

[1] 0.0000000007686286

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)
  
}

Binary Goodness

k <- 128  # number of persons matching hypothesized / group-level pattern
n <- 128+49 # 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.6328624 0.7085935 0.7769709

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

Continuous Virtuosity

k <- 130  # number of persons matching hypothesized / group-level pattern
n <- 130+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.6637135 0.7383632 0.8045365

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

Primary Plots

# re-organize data to plot vignette-by-vignette variability

## binary
S2_vignettes_bin <- S2_bin %>% 
  select(c(ResponseId, 
           theft_hyp, drugs_hyp, cheat_hyp, sex_hyp, sexualize_hyp, blame_hyp, workethic_hyp, revenge_hyp,
           criticism_hyp, health_hyp)) %>% 
  group_by(ResponseId) %>%
  pivot_longer(cols = (c(theft_hyp, drugs_hyp, cheat_hyp, sex_hyp, sexualize_hyp, blame_hyp, workethic_hyp, revenge_hyp,
                         criticism_hyp, health_hyp)),
               names_to = "vignette_name",
               values_to = "fits_hypothesis")

S2_vignettes_bin <- S2_vignettes_bin %>% 
  mutate(fits_hyp = recode(fits_hypothesis, "1" = "Yes", "0" = "No"))
  
bin_labs <- c("Blame", "Cheating", "Criticism", "Drugs", "Health", "Revenge", "Sex", "Sexualization", "Theft", 
              "Work Ethic")
names(bin_labs) <- c("blame_hyp","cheat_hyp", "criticism_hyp", "drugs_hyp", "health_hyp", "revenge_hyp", "sex_hyp",
                      "sexualize_hyp", "theft_hyp", "workethic_hyp")

# continuous
S2_vignettes_cont <- S2_cont %>% 
  select(c(ResponseId, 
           theft_hyp, drugs_hyp, cheat_hyp, sex_hyp, sexualize_hyp, blame_hyp, workethic_hyp, revenge_hyp,
           criticism_hyp, health_hyp,
           theft_opp, drugs_opp, cheat_opp, sex_opp, sexualize_opp, blame_opp, workethic_opp, revenge_opp,
           criticism_opp, health_opp)) %>%
  mutate(theft = case_when(
    theft_hyp == 1 ~ "Yes",
    theft_opp == 1 ~ "No - Opposite",
    (theft_hyp == 0 & theft_opp == 0) ~ "No - Equal")) %>%
  mutate(drugs = case_when(
    drugs_hyp == 1 ~ "Yes",
    drugs_opp == 1 ~ "No - Opposite",
    (drugs_hyp == 0 & drugs_opp == 0) ~ "No - Equal")) %>%
  mutate(cheat = case_when(
    cheat_hyp == 1 ~ "Yes",
    cheat_opp == 1 ~ "No - Opposite",
    (cheat_hyp == 0 & cheat_opp == 0) ~ "No - Equal")) %>%
  mutate(sex = case_when(
    sex_hyp == 1 ~ "Yes",
    sex_opp == 1 ~ "No - Opposite",
    (sex_hyp == 0 & sex_opp == 0) ~ "No - Equal")) %>%
  mutate(sexualize = case_when(
    sexualize_hyp == 1 ~ "Yes",
    sexualize_opp == 1 ~ "No - Opposite",
    (sexualize_hyp == 0 & sexualize_opp == 0) ~ "No - Equal")) %>%
  mutate(blame = case_when(
    blame_hyp == 1 ~ "Yes",
    blame_opp == 1 ~ "No - Opposite",
    (blame_hyp == 0 & blame_opp == 0) ~ "No - Equal")) %>%
  mutate(workethic = case_when(
    workethic_hyp == 1 ~ "Yes",
    workethic_opp == 1 ~ "No - Opposite",
    (workethic_hyp == 0 & workethic_opp == 0) ~ "No - Equal")) %>%
  mutate(revenge = case_when(
    revenge_hyp == 1 ~ "Yes",
    revenge_opp == 1 ~ "No - Opposite",
    (revenge_hyp == 0 & revenge_opp == 0) ~ "No - Equal")) %>%
  mutate(criticism = case_when(
    criticism_hyp == 1 ~ "Yes",
    criticism_opp == 1 ~ "No - Opposite",
    (criticism_hyp == 0 & criticism_opp == 0) ~ "No - Equal")) %>%
  mutate(health = case_when(
    health_hyp == 1 ~ "Yes",
    health_opp == 1 ~ "No - Opposite",
    (health_hyp == 0 & health_opp == 0) ~ "No - Equal")) %>%
  group_by(ResponseId) %>%
  pivot_longer(cols = (c(theft, drugs, cheat, sex, sexualize, blame, workethic, revenge,
                         criticism, health)),
               names_to = "vignette_name",
               values_to = "fits_hyp")

cont_labs <- c("Blame", "Cheating", "Criticism", "Drugs", "Health", "Revenge", "Sex", "Sexualization", "Theft", 
              "Work Ethic")
names(cont_labs) <- c("blame","cheat", "criticism", "drugs", "health", "revenge", "sex",
                      "sexualize", "theft", "workethic")

Binary Goodness

print(S2_bin_vigplot <- ggplot(S2_vignettes_bin, 
                               aes(x = as.character(fits_hyp), 
                                  fill = as.character(fits_hyp))) +
  geom_bar() +
  scale_fill_manual(values = c("#8b1616", "#CEA335"), name = "Fits Hypothesis")+
  scale_y_continuous(limits = c(-0.1,200.1), breaks = c(0, 25, 50, 75, 100, 125, 150, 175, 200)) +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~ vignette_name, ncol = 5, labeller = labeller(vignette_name = bin_labs)) +
  theme(axis.title.x = element_blank(), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_blank(), 
              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("S2_bin_vigplot.png")

Saving 14.1 x 9.04 in image

Continuous Virtuosity

print(S2_cont_vigplot <- ggplot(S2_vignettes_cont, 
                               aes(x = as.character(fits_hyp), 
                                  fill = as.character(fits_hyp))) +
  geom_bar() +
  scale_fill_manual(values = c("grey30","#8b1616", "#CEA335"), name = "Fits Hypothesis")+
  scale_y_continuous(limits = c(-0.1,150.1), breaks = c(0, 25, 50, 75, 100, 125, 150)) +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~ vignette_name, ncol = 5, labeller = labeller(vignette_name = cont_labs)) +
  theme(axis.title.x = element_blank(), 
              axis.title.y = element_text(size = 18),
              axis.text.x = element_blank(), 
              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("S2_cont_vigplot.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
S2_bin <- S2_bin %>%
  mutate(fits_hyp = if_else(fits_hypothesis == 1, "Yes", "No"))

# create yes/no variable for simple hypotheses
S2_cont <- S2_cont %>%
  mutate(fits_hyp = if_else(fits_hypothesis == 1, "Yes", "No"))

Education

Binary Goodness

print(S2_bin_educ <- ggplot(S2_bin, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted 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)))

Continuous Virtuosity

print(S2_cont_educ <- ggplot(S2_cont, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted 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

Binary Goodness

print(S2_bin_pol <- ggplot(S2_bin, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~PoliticalParty) +
  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)))

Continuous Virtuosity

print(S2_cont_pol <- ggplot(S2_cont, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~PoliticalParty) +
  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

Binary Goodness

print(S2_bin_relig <- ggplot(S2_bin, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Religious) +
  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)))

Continuous Virtuosity

print(S2_cont_relig <- ggplot(S2_cont, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~Religious) +
  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

Binary Goodness

print(S2_bin_urb <- ggplot(S2_bin, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~LivingArea) +
  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)))

Continuous Virtuosity

print(S2_cont_urban <- ggplot(S2_cont, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted Pattern") +
  ylab("Number of Participants\n") +
  theme_classic() +
  facet_wrap(~LivingArea) +
  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

Binary Goodness

print(S2_bin_gender <- ggplot(S2_bin, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted 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)))

Continuous Virtuosity

print(S2_cont_gender <- ggplot(S2_cont, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted 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
S2_bin <- S2_bin %>%
  mutate(Race = coalesce(RaceEthnicity_1, RaceEthnicity_2, RaceEthnicity_3, RaceEthnicity_4, 
                         RaceEthnicity_5, RaceEthnicity_6, RaceEthnicity_7, RaceEthnicity_8))

S2_cont <- S2_cont %>%
  mutate(Race = coalesce(RaceEthnicity_1, RaceEthnicity_2, RaceEthnicity_3, RaceEthnicity_4, 
                         RaceEthnicity_5, RaceEthnicity_6, RaceEthnicity_7, RaceEthnicity_8))

Binary Goodness

print(S2_bin_race <- ggplot(S2_bin, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted 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)))

Continuous Virtuosity

print(S2_cont_race <- ggplot(S2_cont, aes(x = fits_hyp)) +
  geom_bar(fill = "grey30") +
  xlab("\nMatched Predicted 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)))