Heroes Factory

Jean Monéger & Roger Giner-Sorolla

2025-08-14

Visit project site

Visit Project Site

Visit project site

▲ Top

R session Environment

Loading packages:

if(!require("dplyr")) install.packages("dplyr")
if(!require("tidyr")) install.packages("tidyr")
if(!require("stringr")) install.packages("stringr")
if(!require("ggplot2")) install.packages("ggplot2")
if(!require("emmeans")) install.packages("emmeans")
if(!require("data.table")) install.packages("data.table")
if(!require("PerformanceAnalytics")) install.packages("PerformanceAnalytics")
if(!require("interactions")) install.packages("interactions")
if(!require("car")) install.packages("car")
if(!require("effectsize")) install.packages("effectsize")
if(!require("RColorBrewer")) install.packages("RColorBrewer")
if(!require("report")) install.packages("report")
if(!require("ordinal")) install.packages("ordinal")
if(!require("robustbase")) install.packages("robustbase") # For Minimum Covariance distance
if(!require("olsrr")) install.packages("olsrr")
if(!require("knitr")) install.packages("knitr")
if(!require("kableExtra")) install.packages("kableExtra")
if(!require("gt")) install.packages("gt")
if(!require("lavaan")) install.packages("lavaan")
if(!require("ppcor")) install.packages("ppcor") # For partial correlations
if(!require("sjPlot")) install.packages("sjPlot")
if(!require("MASS")) install.packages("MASS")
if(!require("Matrix")) install.packages("Matrix")
if(!require("ggrepel")) install.packages("ggrepel")
if(!require("corrplot")) install.packages("corrplot")
if(!require("plotly")) install.packages("plotly") # Interactive plots <3
if(!require("performance")) install.packages("performance") 
if(!require("partR2")) install.packages("partR2") 
if(!require("multcomp")) install.packages("multcomp") 
if(!require("parameters")) install.packages("parameters") 
if(!require("lme4")) install.packages("lme4") 
if(!require("lmerTest")) install.packages("lmerTest") 
if(!require("MuMIn")) install.packages("MuMIn") 
if(!require("gridExtra")) install.packages("gridExtra") # Plot layout
if(!require("scatterplot3d")) install.packages("scatterplot3d") # Really cool 3D scatterplots

Session info:

sessionInfo()

## R version 4.2.2 (2022-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS 14.7.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] scatterplot3d_0.3-44       gridExtra_2.3             
##  [3] MuMIn_1.47.5               lmerTest_3.1-3            
##  [5] lme4_1.1-35.1              parameters_0.25.0         
##  [7] multcomp_1.4-28            TH.data_1.1-3             
##  [9] survival_3.8-3             mvtnorm_1.2-4             
## [11] partR2_0.9.2               performance_0.13.0        
## [13] plotly_4.10.4              corrplot_0.95             
## [15] ggrepel_0.9.6              Matrix_1.6-5              
## [17] sjPlot_2.8.17              ppcor_1.1                 
## [19] MASS_7.3-60.0.1            lavaan_0.6-19             
## [21] gt_1.0.0                   kableExtra_1.4.0          
## [23] knitr_1.50                 olsrr_0.6.1               
## [25] robustbase_0.99-2          ordinal_2023.12-4.1       
## [27] report_0.6.1               RColorBrewer_1.1-3        
## [29] effectsize_1.0.0           car_3.1-3                 
## [31] carData_3.0-5              interactions_1.2.0        
## [33] PerformanceAnalytics_2.0.8 xts_0.14.1                
## [35] zoo_1.8-14                 data.table_1.17.0         
## [37] emmeans_1.11.1             ggplot2_3.5.2             
## [39] stringr_1.5.1              tidyr_1.3.1               
## [41] dplyr_1.1.4               
## 
## loaded via a namespace (and not attached):
##  [1] minqa_1.2.6         sjlabelled_1.2.0    estimability_1.5.1 
##  [4] rstudioapi_0.17.1   listenv_0.9.1       furrr_0.3.1        
##  [7] farver_2.1.2        xml2_1.3.8          codetools_0.2-20   
## [10] splines_4.2.2       mnormt_2.1.1        cachem_1.1.0       
## [13] sjmisc_2.8.10       Formula_1.2-5       jsonlite_2.0.0     
## [16] nloptr_2.0.3        ggeffects_2.2.1     broom_1.0.8        
## [19] broom.mixed_0.2.9.6 compiler_4.2.2      httr_1.4.7         
## [22] sjstats_0.19.0      backports_1.5.0     fastmap_1.2.0      
## [25] lazyeval_0.2.2      cli_3.6.5           htmltools_0.5.8.1  
## [28] tools_4.2.2         coda_0.19-4.1       gtable_0.3.6       
## [31] glue_1.8.0          Rcpp_1.0.14         jquerylib_0.1.4    
## [34] vctrs_0.6.5         svglite_2.1.3       nlme_3.1-164       
## [37] insight_1.2.0       xfun_0.52           globals_0.18.0     
## [40] lifecycle_1.0.4     goftest_1.2-3       future_1.40.0      
## [43] DEoptimR_1.1-3-1    scales_1.4.0        parallel_4.2.2     
## [46] sandwich_3.1-1      yaml_2.3.10         pander_0.6.6       
## [49] sass_0.4.10         stringi_1.8.7       ucminf_1.2.1       
## [52] bayestestR_0.15.3   nortest_1.0-4       boot_1.3-31        
## [55] rlang_1.1.6         pkgconfig_2.0.3     systemfonts_1.1.0  
## [58] evaluate_1.0.3      lattice_0.22-7      purrr_1.0.4        
## [61] htmlwidgets_1.6.4   tidyselect_1.2.1    parallelly_1.44.0  
## [64] magrittr_2.0.3      R6_2.6.1            generics_0.1.4     
## [67] pillar_1.10.2       withr_3.0.2         jtools_2.3.0       
## [70] datawizard_1.0.2    abind_1.4-8         tibble_3.2.1       
## [73] rmarkdown_2.29      grid_4.2.2          pbivnorm_0.6.0     
## [76] forcats_1.0.0       digest_0.6.37       xtable_1.8-4       
## [79] numDeriv_2016.8-1.1 stats4_4.2.2        viridisLite_0.4.2  
## [82] bslib_0.9.0         quadprog_1.5-8

Introduction

Certain occupations, like healthcare workers, military personnel, and other first responders, are often labeled as “heroes,” especially during crises. While this label brings praise, it has been established that heroism perception might justify workers’ exploitation, burden workers with high expectations, and limit their ability to define their working conditions. We report a set of three pre-registered studies testing danger and altruism in different definitions, as contributing elements of heroism perception in entire occupations. Studies 1 and 2 explored whether framing occupations as high in risk or altruism causally influenced perceived heroism. An additional correlational study emphasized the centrality of risk and altruism perception in explaining heroism in a representative set of 53 occupations. Taken together, these studies found that risk and altruism are major, relatively independent components of heroism perception, that can be manipulated to increase or decrease heroism perception in occupations that are not already subject to widespread discourse about their heroism.

Perception of physical threat contributes to the perception of heroism

Altruistic motivations of the group’s member contribute to the perception of heroism

Occupations differ in the manner they are depicted in the media. Whereas representations of psychiatrists might emphasize their intellectual curiosity and over-analytical interactions, fictional and news portrayals might paint occupations such as firefighters as heroes: brave and selfless individuals, ready to put their life at risk for the sake of others (Parkinson et al., 2023). Interestingly, this heroic framing might change over time, as illustrated by the fluctuations and construction of the image of New-York firefighters through time (Hochbruck, 2019), or the evolution of soviet labour heroes throughout the history of USSR (Neutatz & Tibilova, 2024). For example, in a crisis situation such as the recent covid-19 pandemic, entire occupations, such as frontline healthcare workers, might be typecast as heroes (Franco, 2016; Mayor et al., 2012; Wagner-Egger et al., 2011), but also less socially valued occupations such as gig workers (Cameron & Fekete, 2020; Cameron et al., 2022; Friedland & Balkin, 2023).

Framing any occupation with a moral role, be it hero, victim or indeed villain, carries consequences, as moral roles have been studied as a central feature of moral social cognition (Bergstrand & Jasper, 2018; Klapp, 1954; Hartman et al., 2022). Moral roles such as heroes have been suggested to be rhetoric tools to support social control by creating norms for emulation, directing collective responses through propaganda, and affirming shared social norms (Brink et al., 2022; Klapp, 1954; Neutatz & Tibilova, 2024; Parkinson et al., 2023). Moreover, the heroization of professions has been shown to carry ironic negative consequences for the workers. It was reported that “heroes” are more likely to be exploited (Stanley & Kay, 2024). For instance, veterans described as heroes would be more likely to be funneled in low wage positions (Stanley et al., 2023). Perception of healthcare workers as heroes, rather than victims, also appear to elicit sympathy, but reduce perceptions of inequalities and buffer intents to support the workers (Yuan et al., 2021).

Workers appear to be quite aware of possible ironic consequences of being framed as heroes, and calls from professionals to tune down the hero rhetoric have been heard nurses (Cox, 2020; Hennekam et al., 2020; Pangborn et al., 2023), but also from firefighters (Cutcher & Dyer, 2020). Tangherlini (2000) noted a tradition of resisting the hero rhetoric by sharing anti-heroic anecdotes among paramedics and other first-responders, thus underlining workers’ willingness to step away from the heroic frame. In particular, Hennekam et al. (2020) noted in their interviews with nurses during the covid-19 pandemic, that most healthcare workers were skeptical and pessimistic regarding the hero frame. It has been suggested that the framing might serve the purpose of masking the responsibility of the government (Mohammed et al., 2021). Moreover, accolades such as ritually applauding the workers every night might taste bitter coming from a population that also flouts regulations aiming at containing the virus (Hennekam et al., 2020; Pangborn et al., 2023:, cf. analysis of the so-called hero contract, Sumner & Kinsella, 2020).

For some occupations, such as fire and rescue workers, hero status seems natural. In others, such as policing and the military, individual heroes have always been recognized, but collective heroism is only claimed in some places and times (de Camargo & Whiley, 2022; Giner-Sorolla et al., 2020). And in dangerous times, such as the COVID-19 pandemic, the hero badge has been awarded not just to normally valued members of society such as healthcare workers, but to less likely key worker occupations such as delivery drivers (Cameron & Fekete, 2020; Cameron et al., 2022; Friedland & Balkin, 2023). But labelling an occupation as “heroes,” however well-intended, might lead to a typecasting process that justifies exploitation. As such, it is urgent to understand the appraisals that facilitate awarding this moral role.

What exactly is a hero? Several definitions co-exist in the literature. “To act in a prosocial manner despite personal risk” (Franco et al., 2011, p. 99), “A hero is a person who knowingly, and voluntarily, acts for the good of one or more people at significant risk to the self, without being motivated by reward” (Riches, 2018, p.585), or “Heroism consists of actions undertaken to help others, despite the possibility that they may result in the helper’s death or injury” (Becker & Eagly, 2004, p. 163). All these definitions share two common elements: heroes put themselves at risk, and heroes are altruistic – in particular if they have a positive influence on outcomes and are well-intended (Harvey et al., 2009). Indeed, one study of popular prototypes of the hero indicated that brave and selfless were among the most cited traits of heroes among laypeople (Kinsella et al., 2015).

If risk and altruism are the building blocks of heroism, then it should be possible to deploy narratives manipulating these two aspects to foster the perception of occupations as heroic – or, if heroism narratives are having a negative effect, to promote a more mundane appraisal. However, from the previous literature it is not clear how best a collective can effectively be associated with risk and altruism. Which is more important for heroic groups: that their motivations and moral character are brave and selfless; or that their actual behavior reflects these traits– that is, actually exposed to danger and actually helping people? Although the perception of moral characteristics and situation characteristics might overlap, it is unclear which elements would be the most effective in turning occupations such as underwater welders into heroes; or conversely, making firefighters appear to be less heroic. Moral psychology often emphasizes that intentions can be inferred from outcomes under certain conditions (Aiyer et al., 2024; Baron & Hershey, 1988; Gino et al., 2009; Knobe, 2003; Nakamura, 2018), whereas it is unlikely that outcomes are inferred solely from intentions. Moreover, while it is possible that laypeople do not see selflessness as exceptional, given possible self-enhancement strategy whereby individuals would believe they are themselves more selfless than average (Sedikides et al., 2025; Xiao et al., 2021), describing altruistic outcomes – rather than characters – of entire occupations might create a more undeniable perception of exceptionality of the workers. However, to our knowledge, no direct comparisons of the effectiveness of framing altruistic intentions vs. altruistic outcomes on collective perceptions have been made in the literature.

Another, more theoretical question arises in determining whether danger and altruism create heroes additively or interactively. If both traits are absolutely required to gain the label of hero, it would make sense to find an interaction effect layered on the main effects of each trait, such that occupations with one trait but the other do not proportionately become “halfway heroes.” But it is also possible that the traits, instead of being necessary, are each partially sufficient, predicting an additive effect without interaction; or, less likely, the situation where each trait is sufficient to create hero perceptions would be marked by two main effects and an interaction in the opposite direction, such that the combination of the two increases heroism perceptions less than their simple additive effect would predict.

On the basis that risk and altruism have been consistently identified as central elements of heroism, we investigated 1) the causal link between risk perception and altruism perception on heroism perception of collectives, 2) the relationship between these two elements (i.e., whether they interact together or have independent effects on heroism), 3) the relative importance of character elements – and particularly altruistic motives – over situational elements – and particularly altruistic outcomes – in predicting heroism.

In particular, Study 1 assessed the centrality and causal effects of exposure to risk and altruistic motives on heroism perceptions of occupations, Study 2 contrasted altruistic motives and altruistic outcomes to assess their causal effects on occupations framed as exposed to risk, and study 3 investigated the universality of heroism perceptions being driven by risk and altruism across a representative subsample of 53 occupations. Throughout all the studies, we investigated whether risk and altruism interacted vs. had independent effects on heroism perception. All studies were pre-registered and conducted using representative UK samples.

Study 1

We assessed how framing the occupations as exposed to physical (vs no) risk and altruistic (vs self-centered) motivations would influence people typecasting them as heroic. As well, we used manipulation checks as correlational predictors of heroism judgments: evaluations that the workers are brave, selfless, exposed to danger, and helpful. Additionally, we measured attitude toward target occupations to account for possible Halo effects in correlational analyses.

We focused on five occupations viewed with varying degrees of heroism: firefighters, nurses, police officers, psychiatrists, and underwater welders. The latter occupation was chosen because it is lesser known than the other occupations. We reasoned that most participants would not be familiar with the work environments, motivations, and character of underwater welders. To manipulate perception of risk and altruism, study 1 used narrative vignettes emphasizing exposure to physical risk (vs. a boredom condition emphasizing the lack of physical risk) and the selfless motives of workers (vs. their motives of self-improvement). The design of the study was between-participant: each participant was randomly assigned to one target occupation.

We predicted that risk and altruism manipulations will predict heroism perception across all five occupations. We explored the possible interaction between these two hypothetical components of collective heroism perception. The study was pre-registered (https://osf.io/ysvxb).

Protocol

As registered, we used a representative UK sample of 1360 participants using Prolific. Two participants failed the comprehension checks, two participants timed-out despite completing the survey (see Prolific policies), and three participants returned their surveys despite completing it. To account for this attrition, we sampled seven more participants to result in a complete sample of 1360 participants (692 women, 653 men, 11 others, and 4 preferred not to say, Mage = 46.55, SD = 15.3). Participants were recruited using the representative sampling procedure offered by Prolific. Prolific representative sampling stratifies the participants according to sex (Female, Male), ethnicity (Asian, Black, Other, White) and age (18-24, 25-34, 35-44, 45-54, 55-100) to match the characteristics of the target population using census data from the UK Office of National Statistics. Participants received £9 per hour.

After providing informed consent, entering their Prolific ID, and completing a commitment check (Geisen, 2022; Peer et al., 2024), participants were randomly assigned (between-subjects) to read a vignette describing one of the target occupations. Each vignette comprised three paragraphs: the first introduced the occupation; the second described the level of risk associated with the job (physical risk vs. boredom); and the third outlined the workers’ primary motivation (altruistic vs. self-improvement). We used deceptive information to manipulate risk and motivations. Specifically, we reported false results from a survey indicating that “83% of XXX thought that their life was at risk in the past 12 months” (vs “83% of XXX reported being bored most of the time”), and “74% of XXX identified”Helping people” as their primary motivation” (vs “74% of XXX identified”self-improvement” as their primary motivation”). These two paragraphs were presented in randomized order. Conditions were randomly assigned.

Participants then completed three comprehension check questions. Those who failed twice were redirected to a page instructing them to return their submission. Participants who passed continued to the main dependent measures.

Material

Heroism evaluation. Participants evaluated to what extent they agreed with the target occupation being described as “Heroes” using a 7-point Likert scale from “Strongly disagree” to “Strongly agree”.

Character attribution evaluation. Participants will evaluate the applicability of two adjectives to describe the character of the target workers using two 7 points bipolar scales. One will be ranging from cowardly to brave, and the other will be ranging from Selfish to selfless.

Occupation descriptors evaluation. Participants then rated to what extent they believed the target occupation is objectively associated with physical danger, and to what extent the profession is helpful. Both questions used 7-point Likert scales from “Not at all” to “Extremely”.

General attitude. Participants reported their general attitude toward each target occupation using a 7-point-scale from “Very negative” to “Very positive”.

Comprehension checks. Participants needed to respond to three multiple choice questions asking them about the content of the vignette they just read (see material for details). They had two chances of getting all comprehension checks correct - after which they were urged to return their survey. The comprehension checks were displayed right beneath the descriptions and consisted of easy questions about the content of the vignette.

Credibility check. Participants answered the following question: “In your opinion, how believable was the information that you read at the beginning of this study, about the target occupation’s motivation and working conditions?” using a 7-point scale from “1- Very unbelievable” to “7 - Very believable”.

Results

The following code (that you can unfold by clicking on ‘show code’) show how the data was wrangled from the raw qualtrics data to the clean usable data – users can run it if they wish to reproduce our data wrangling procedure by using the raw data file as an input and changing eval = FALSE to eval = TRUE.

Set <- read.csv("Raw_Qualtrics_Data_HeroFactory_April2025.csv", comment.char="#")
Set <- subset(Set, Set$Q233 != "No, I won't")
Set <- Set[-c(1:2),]
Set <- subset(Set, Set$FailedComp == "No") # Keep only participants who succeeded the completion checks

# Two participants timed out but completed the survey. Let's remove them to not threaten representativity and stick closely to registration:

Set <- subset(Set, !(Set$Prol_ID %in% c("5cb6f38ffdc7fa0013f809a3", "5cf84c0f4b639a0016a45a54")))  

QualCheck <- Set[, c(21:26)] # Quality evaluation is Good. There are 8 participant flagged as duplicatated, but I do not trust this opaque quality check

Demographics <- Set[, 38:48] # Keeping demographics aside
# First, add a unique row identifier if one does not exist already:
job_cols <- c("Job_match_1", "Job_match_2", "Job_match_3", 
              "Job_match_4", "Job_match_5", "Job_match_7", "Job_match_6")

Demographics$JOB <- apply(Demographics[, job_cols], 1, function(x) {
  # Remove empty strings
  reported_jobs <- x[x != ""]
  # If no jobs reported, assign NA; else paste the jobs together
  if (length(reported_jobs) == 0) NA else paste(reported_jobs, collapse = ", ")
})

Set$Credibility[Set$Credibility == "Very unbelievable"] <- 1
Set$Credibility[Set$Credibility == "Quite unbelievable"] <- 2
Set$Credibility[Set$Credibility == "Somewhat unbelievable"] <- 3
Set$Credibility[Set$Credibility == "Neutral"] <- 4
Set$Credibility[Set$Credibility == "Somewhat believable"] <- 5
Set$Credibility[Set$Credibility == "Quite believable"] <- 6
Set$Credibility[Set$Credibility == "Very believable"] <- 7

Set <- Set %>%
  dplyr::mutate(across(
    where(is.character),
    ~ dplyr::case_match(.,
                        "83% of the workers reported feeling bored most of the time" ~ "Bored",
                        "83% of the workers reported that their life was at risk in the past 12 months" ~ "Risk",
                        "74% of the workers identified \"helping people\" as their primary motivation" ~ "Help",
                        "74% of the workers identified \"self-improvement\" as their primary motivation" ~ "Self",
                        .default = .
    )
  ))


# I need to code the condition. It can be retrieved using non-empty heroism ratings. e.g., NRH is Nurses, Risk, Help. Note that it can also be re-computed from succeeded Comprehension checks (AC___1, AC____2, and AC_____3)

# I concatenate Comprehension checks (AC1, AC2, and AC3)
Set_long <- Set %>%
  # Gather all columns matching the pattern
  pivot_longer(
    cols = matches("^AC.*[123]$"),
    names_to = c("base", "suffix"),
    names_pattern = "^(AC.*?)([123])$",
    values_to = "value"
  ) %>%
  # Now spread the suffix into separate columns.
  pivot_wider(
    names_from = suffix,
    values_from = value,
    names_prefix = "AC"
  )



# Identify all columns ending with "h_1" and use it to fill the column "CONDITION"
h1_cols <- grep("_h_1$", names(Set_long), value = TRUE)
# Creating the condition
Set_long_filtered <- Set_long %>%
  rowwise() %>%  # process each row individually
  mutate(
    # Count how many of the h1_cols are non-empty ("")
    non_empty_count = sum(c_across(all_of(h1_cols)) != ""),
    
    # Determine the Condition based on non-empty count
    Condition_2 = if (non_empty_count == 1) {
      # Get the name of the column that is non-empty
      selected_col <- h1_cols[which(c_across(all_of(h1_cols)) != "")]
      # Remove the "_h_1" suffix to get the condition (e.g., "NRS" from "NRS_h_1")
      str_remove(selected_col, "_h_1")
    } else {
      "Error"
    }
  ) %>%
  ungroup() %>%  # exit rowwise mode
  dplyr::select(-non_empty_count) %>% distinct(ResponseId, .keep_all=T) # remove the helper column if no longer needed



# Now I'll neatly create the main data frame. I create a helper function that extract all non empty values and apply it on all rows ending with our codes: 
# - _h_1 is heroism rating, _m_1_1 is first item of MC1, m_1_2, is second item of mc2, _at is attitude etc.


# Helper function: Given a vector, return the only non-empty element (or NA if none/few found)
extract_value <- function(x) {
  non_empty <- x[x != ""]
  if(length(non_empty) == 1) return(non_empty)
  else return(NA)
}

# Create the final data frame
final_df <- Set_long_filtered %>%
  dplyr::mutate(
    # Extract the non-empty value from each group of columns:
    Heroism = apply(dplyr::select(., ends_with("_h_1")), 1, extract_value),
    Danger        = apply(dplyr::select(., ends_with("m1_1")), 1, extract_value),
    Helpfulness        = apply(dplyr::select(., ends_with("m1_2")), 1, extract_value),
    Selfless        = apply(dplyr::select(., ends_with("m2_1")), 1, extract_value),
    Brave        = apply(dplyr::select(., ends_with("m2_2")), 1, extract_value),
    Attitude     = apply(dplyr::select(., matches("(_at)$")), 1, extract_value)
  ) %>%
  # Select only the required columns in the final data frame
  dplyr::select(ResponseId, Condition, Condition_2, Heroism, Danger, Helpfulness, Selfless, Brave, Attitude, Credibility, Gender, Age)

# Recoding values:


final_df <- final_df %>%
  # 1. Split Condition into Job, Risk, and Help
  mutate(
    Job  = if_else(str_sub(Condition_2, 1, 2) == "Ps", "Ps", str_sub(Condition_2, 1, 1)),
    Risk = if_else(str_sub(Condition_2, 1, 2) == "Ps", str_sub(Condition_2, 3, 3), str_sub(Condition_2, 2, 2)),
    Help = if_else(str_sub(Condition_2, 1, 2) == "Ps", str_sub(Condition_2, 4, 4), str_sub(Condition_2, 3, 3))
  ) %>%
  # 2. Dummy code Risk and Help: R -> 0.5, B -> -0.5; H -> 0.5, S -> -0.5
  mutate(
    Risk_dummy = case_when(
      Risk == "R" ~ 0.5,
      Risk == "B" ~ -0.5,
      TRUE ~ NA_real_
    ),
    Help_dummy = case_when(
      Help == "H" ~ 0.5,
      Help == "S" ~ -0.5,
      TRUE ~ NA_real_
    )
  ) %>%
  # 3. Recode rating items (HeroismScore, MC1_1, MC1_2, MC2_1, MC2_2)
  mutate_at(vars(Heroism, Danger, Helpfulness, Selfless, Brave), 
            ~ case_when(
              . == "1 - Strongly disagree" ~ "1",
              . == "7 - Strongly agree"    ~ "7",
              TRUE ~ .
            )) %>%
  # Optionally convert these recoded columns to numeric
  mutate_at(vars(Heroism, Danger, Helpfulness, Selfless, Brave), as.numeric) %>%
  # 4. Recode Attitude values
  mutate(
    Attitude = case_when(
      Attitude == "Very negative"    ~ 1,
      Attitude == "Quite negative" ~ 2,
      Attitude == "Somewhat negative" ~ 3,
      Attitude == "Neutral"           ~ 4,
      Attitude == "Somewhat positive" ~ 5,
      Attitude == "Quite positive" ~ 6,
      Attitude == "Very positive" ~ 7,
      TRUE ~ NA_real_
    )
  ) 

final_df$Part_Job <- Demographics$JOB
Set<-final_df

#write.csv(Demographics, "Demog_Hero_factory_2025.csv", row.names = F)
#write.csv(Set, "DF_HeroFactory_April2025.csv", row.names = F)

Demographics

Set <- read.csv("DF_HeroFactory_April2025.csv", comment.char="#") # This is the clean Data /
# If you want to start with the raw qualtrics data, use the chunk above with eval = T

Demographics <- read.csv("Demog_Hero_factory_2025.csv") # This is the clean demographic 


Demographics$flag_inconsistent <- (Demographics$Job_match_6 == "None of the above") & 
  (Demographics$Job_match_1 != "" | Demographics$Job_match_2 != "" | Demographics$Job_match_3 != "" | 
   Demographics$Job_match_4 != "" | Demographics$Job_match_5 != "")




paste0("1360 participants took part in the study. Mean age in the sample is ", mean(as.numeric(Set$Age)), ", SD = ", sd(as.numeric(Set$Age)))

## [1] "1360 participants took part in the study. Mean age in the sample is 46.5455882352941, SD = 15.359732853326"

## Gender

Set %>% group_by(Gender) %>% summarise(N=n()) %>%
  ggplot(aes(x=Gender,y=N,fill=Gender))+
  geom_bar(stat = 'identity',color='black')+
  scale_y_continuous(labels = scales::comma_format(accuracy = 2))+
  geom_text(aes(label=N),vjust=-0.25,fontface='bold')+
  theme_bw()+
  theme(axis.text = element_text(color='black',face='bold'),
        axis.title = element_text(color='black',face='bold'),
        legend.text = element_text(color='black',face='bold'),
        legend.title = element_text(color='black',face='bold')) +
  ggtitle("Gender distribution")

## Occupations
#colnames(Set)
jobs <- unlist(Demographics[-which(Demographics$flag_inconsistent == T), 4:10])           # Make a long list of all jobs that were named

jobs <- jobs[jobs != ""]     # Remove empty strings

job_df <- as.data.frame(table(jobs))
colnames(job_df) <- c("Job", "Count")

ggplot(job_df, aes(x = Job, y = Count, fill = Job)) +
  geom_bar(stat = 'identity',color='black')+
  scale_y_continuous(labels = scales::comma_format(accuracy = 2))+
  geom_text(aes(label=Count),vjust=-0.25,fontface='bold')+
  theme_bw()+
  theme(axis.text = element_text(color='black',face='bold'),
        axis.title = element_text(color='black',face='bold'),
        legend.text = element_text(color='black',face='bold'),
        legend.title = element_text(color='black',face='bold')) +
  ggtitle("Job distribution")

Descriptives

⚠️ NOTE: In our report, the following codes are used:

• F: Firefighter condition

• N: Nurse condition

• P: Police officer condition

• Ps: Psychiatrists condition

• W: Welders condition

• B: Boredom condition

• R: Risk condition

• H: Help condition

• S: Self-improvement condition

# 1. Create a summary dataframe for facet annotations (mean and SD)
df_summary <- Set %>%
  group_by(Job) %>%
  summarize(
    mean_score = mean(Heroism, na.rm = TRUE),
    sd_score   = sd(Heroism, na.rm = TRUE),
    .groups = "drop"
  )

# 2. Create the ggplot using the long format data
ggplot(Set, aes(x = Heroism)) +
  geom_histogram(aes(fill = after_stat(count)),
                 binwidth = 1,
                 color = "black", show.legend = FALSE) +
  facet_grid( ~ Job, scales = "free") +
  scale_fill_gradientn(
    colours = brewer.pal(9, "YlOrBr"),
    name = "Count"
  ) +
  labs(
    title = "Histograms of Variable by Occupation",
    x = "Score",
    y = "Count"
  ) +
  # Annotate each facet with the mean and standard deviation
  geom_text(data = df_summary,
            aes(x = 7, y = Inf,
                label = paste0("Mean = ", round(mean_score, 2),
                               "\nSD = ", round(sd_score, 2))),
            vjust = 1.5, hjust = 1.1, size = 3) +
  theme_classic() +
  theme(panel.grid.major.y = element_line(linewidth = 0.5),
        panel.grid.minor.y = element_line(linewidth = 0.5))

Main reported analyses - as reported

Distribution of credibility scores within each occupation indicated that globally, participants thought the descriptions were credible, with the lowest credibility being in the nurse condition, M = 4.9, SD = 1.69 (see SOM, for details including influential cases analyses, model diagnostics, and exploratory analyses).

#Fi <- subset(Set, Set$Job == "F") # Firef subset
#N <- subset(Set, Set$Job == "N") # Nurses subset
#P <- subset(Set, Set$Job == "P") # Police subset
#Ps <- subset(Set, Set$Job == "Ps") # Psych subset
#W <- subset(Set, Set$Job == "W") # Welders subset

#paste0("Firefighter condition mean credibility:")
#mean(Fi$Credibility)
#paste0("Firefighter condition sd credibility:")
#sd(Fi$Credibility)
#
#paste0("Nurses condition mean credibility:")
#mean(N$Credibility)
#paste0("Nurses condition sd credibility:")
#sd(N$Credibility)
#
#paste0("Police officer condition mean credibility:")
#mean(P$Credibility)
#paste0("Police officer condition sd credibility:")
#sd(P$Credibility)
#
#paste0("Psychiatrists condition mean credibility:")
#mean(Ps$Credibility)
#paste0("Psychiatrists condition sd credibility:")
#sd(Ps$Credibility)
#
#paste0("Welders condition mean credibility:")
#mean(W$Credibility)
#paste0("Welders condition sd credibility:")
#sd(W$Credibility)



credibility_data <- tibble::tribble(
  ~Condition,       ~Mean,      ~SD,
  "Firefighter",    5.274725,   1.595485,
  "Nurse",          4.900369,   1.686888,
  "Police officer", 5.080882,   1.328402,
  "Psychiatrist",   4.992593,   1.532771,
  "Welder",         5.240876,   1.422094
)

# Create gt table
credibility_data %>%
  mutate(
    Mean = round(Mean, 2),
    SD   = round(SD, 2)
  ) %>%
  gt() %>%
  tab_header(
    title = "Credibility Ratings by Occupation"
  ) %>%
  cols_label(
    Condition = "Occupation",
    Mean      = "Mean Credibility",
    SD        = "SD Credibility"
  ) %>%
  fmt_number(
    columns = c(Mean, SD),
    decimals = 2
  ) %>%
  tab_options(
    table.border.top.width = px(2),
    table.border.bottom.width = px(2),
    table.font.size = px(14)
  )

Credibility Ratings by Occupation
Occupation	Mean Credibility	SD Credibility
Firefighter	5.27	1.60
Nurse	4.90	1.69
Police officer	5.08	1.33
Psychiatrist	4.99	1.53
Welder	5.24	1.42

---

Finer grained Credibility analyses: Credibility by Occupations x experimental condition

Set$Credibility <- as.numeric(Set$Credibility)
df_long <- Set %>%
  pivot_longer(
    cols = c(Heroism, Danger, Helpfulness, Selfless, Brave, Attitude, Credibility),
    names_to = "Variable",
    values_to = "Score"
  )

# 2. Compute summary statistics by Job and Variable
df_long2<- subset(df_long, df_long$Variable != "Credibility")

df_summary <- df_long %>%
  group_by(Job, Variable) %>%
  summarize(
    mean_score = mean(Score, na.rm = TRUE),
    sd_score   = sd(Score, na.rm = TRUE),
    .groups = "drop"
  )
df_summary2 <- df_long2 %>%
  group_by(Job, Variable) %>%
  summarize(
    mean_score = mean(Score, na.rm = TRUE),
    sd_score   = sd(Score, na.rm = TRUE),
    .groups = "drop"
  )
# 3. Create the ggplot using the long format data
ggplot(df_long2, aes(x = Score)) +
  geom_histogram(aes(fill = after_stat(count)),
                 binwidth = 1,
                 color = "black", show.legend = FALSE) +
  facet_grid(Variable ~ Job, scales = "free") +
  scale_fill_gradientn(
    colours = brewer.pal(9, "YlOrBr"),
    name = "Count"
  ) +
  labs(
    title = "Histograms of Variable by Occupation",
    x = "Score",
    y = "Count"
  ) +
  # Annotate each facet with the mean and standard deviation
  geom_text(data = df_summary2,
            aes(x = 7, y = Inf,
                label = paste0("Mean = ", round(mean_score, 2),
                               "\nSD = ", round(sd_score, 2))),
            vjust = 1.5, hjust = 1.1, size = 3) +
  theme_classic() +
  theme(
    panel.grid.major.y = element_line(linewidth = 0.5),
    panel.grid.minor.y = element_line(linewidth = 0.5)
  )

Using Welch t-test, our manipulations did successfully influence character-based manipulation checks. Participants in the boredom condition were associated with lower perceived bravery (M = 5.62, SD = 1.28) than participants in the risk condition (M = 5.98, SD = 1.26), and this difference reached significance in a Welch two sample t-test, t(1357.02) = -5.21, p < .001, Cohen’s d = 0.28, 95% CI [0.18, 0.39].

t.test(Set$Brave ~Set$Risk)

## 
##  Welch Two Sample t-test
## 
## data:  Set$Brave by Set$Risk
## t = -5.206, df = 1357, p-value = 2.227e-07
## alternative hypothesis: true difference in means between group B and group R is not equal to 0
## 95 percent confidence interval:
##  -0.4931945 -0.2232317
## sample estimates:
## mean in group B mean in group R 
##        5.622781        5.980994

cohens_d(Set$Brave ~Set$Risk)

my_sum <- Set %>%
  group_by(Risk) %>%
  summarise( 
    n=n(),
    mean=mean(Brave),
    sd=sd(Brave)
  ) %>%
  mutate( se=sd/sqrt(n))
 
# Standard deviation
ggplot(my_sum) +
  geom_bar( aes(x=Risk, y=mean), stat="identity", fill="forestgreen", alpha=0.5) +
  geom_errorbar( aes(x=Risk, ymin=mean-se, ymax=mean+se), width=0.4, colour="orange", alpha=0.9, size=1.5) +
  ggtitle("Brave ~ Risk type (B = Boredom; R = High risk); bars are SE")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

by(Set$Brave, Set$Risk, sd)

## Set$Risk: B
## [1] 1.278198
## ------------------------------------------------------------ 
## Set$Risk: R
## [1] 1.259083

sd(Set$Brave)

## [1] 1.280748

Similarly, participants in the altruistic motivation condition rated target occupation as more selfless (M = 5.67, SD = 1.24) than participants in the self-improvement condition (M = 5.33, SD = 1.40), and this difference also reached significance, t(1336.56) = 4.80, p < .001, d = 0.26, 95% CI [0.15, 0.37].

t.test(Set$Selfless ~ Set$Help)

## 
##  Welch Two Sample t-test
## 
## data:  Set$Selfless by Set$Help
## t = 4.7993, df = 1336.6, p-value = 1.771e-06
## alternative hypothesis: true difference in means between group H and group S is not equal to 0
## 95 percent confidence interval:
##  0.2034586 0.4847767
## sample estimates:
## mean in group H mean in group S 
##        5.669118        5.325000

by(Set$Selfless, Set$Help, mean)

## Set$Help: H
## [1] 5.669118
## ------------------------------------------------------------ 
## Set$Help: S
## [1] 5.325

by(Set$Selfless, Set$Help, sd)

## Set$Help: H
## [1] 1.235551
## ------------------------------------------------------------ 
## Set$Help: S
## [1] 1.403328

my_sum <- Set %>%
  group_by(Help) %>%
  summarise( 
    n=n(),
    mean=mean(Brave),
    sd=sd(Brave)
  ) %>%
  mutate( se=sd/sqrt(n))
 
# Standard deviation
ggplot(my_sum) +
  geom_bar( aes(x=Help, y=mean), stat="identity", fill="forestgreen", alpha=0.5) +
  geom_errorbar( aes(x=Help, ymin=mean-se, ymax=mean+se), width=0.4, colour="orange", alpha=0.9, size=1.5) +
  ggtitle("Selfless ~ Help type (H = Helping; S = Self improve); bars are SE")

For objective assessments also, our manipulations were successful: participants in the risk condition evaluated the target occupation as significantly more exposed to risk (M = 6.28, SD = 0.95) than participants in the boredom condition (M = 5.11, SD = 1.38), t(1190.54) = 18.15, p < .001; Cohen’s d = 0.99, 95% CI [0.87, 1.10].

t.test(Set$Danger ~Set$Risk)

## 
##  Welch Two Sample t-test
## 
## data:  Set$Danger by Set$Risk
## t = -18.15, df = 1190.5, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group B and group R is not equal to 0
## 95 percent confidence interval:
##  -1.294546 -1.041971
## sample estimates:
## mean in group B mean in group R 
##        5.113905        6.282164

by(Set$Danger, Set$Risk, mean)

## Set$Risk: B
## [1] 5.113905
## ------------------------------------------------------------ 
## Set$Risk: R
## [1] 6.282164

by(Set$Danger, Set$Risk, sd)

## Set$Risk: B
## [1] 1.384694
## ------------------------------------------------------------ 
## Set$Risk: R
## [1] 0.9454696

my_sum <- Set %>%
  group_by(Risk) %>%
  summarise( 
    n=n(),
    mean=mean(Danger),
    sd=sd(Danger)
  ) %>%
  mutate( se=sd/sqrt(n))
 
# Standard deviation
ggplot(my_sum) +
  geom_bar( aes(x=Risk, y=mean), stat="identity", fill="forestgreen", alpha=0.5) +
  geom_errorbar( aes(x=Risk, ymin=mean-se, ymax=mean+se), width=0.4, colour="orange", alpha=0.9, size=1.5) +
  ggtitle("Danger ~ Risk type (B = Boredom; R = High risk); bars are SE")

Participants in the altrustic motivation condition indicated that the target occupation was more helping others (M = 6.11, SD = 0.98) than participants in the self-improvement motivation condition (M = 5.80, SD = 1.24), t(1289.51) = 5.01, p < .001; Cohen’s d = 0.27, 95% CI [0.16, 0.38].

t.test(Set$Helpfulness ~Set$Help)

## 
##  Welch Two Sample t-test
## 
## data:  Set$Helpfulness by Set$Help
## t = 5.0127, df = 1289.5, p-value = 6.113e-07
## alternative hypothesis: true difference in means between group H and group S is not equal to 0
## 95 percent confidence interval:
##  0.1852758 0.4235477
## sample estimates:
## mean in group H mean in group S 
##        6.107353        5.802941

by(Set$Helpfulness, Set$Help, mean)

## Set$Help: H
## [1] 6.107353
## ------------------------------------------------------------ 
## Set$Help: S
## [1] 5.802941

by(Set$Helpfulness, Set$Help, sd)

## Set$Help: H
## [1] 0.9822903
## ------------------------------------------------------------ 
## Set$Help: S
## [1] 1.242115

my_sum <- Set %>%
  group_by(Help) %>%
  summarise( 
    n=n(),
    mean=mean(Helpfulness),
    sd=sd(Helpfulness)
  ) %>%
  mutate( se=sd/sqrt(n))
 
# Standard deviation
ggplot(my_sum) +
  geom_bar( aes(x=Help, y=mean), stat="identity", fill="forestgreen", alpha=0.5) +
  geom_errorbar( aes(x=Help, ymin=mean-se, ymax=mean+se), width=0.4, colour="orange", alpha=0.9, size=1.5) +
  ggtitle("Helpfulness ~ Motivation type (H = Helping; S = Self-improve); bars are SE")

As registered, we evaluated the effects of our manipulations using Ordinary Least Square (OLS) linear models. To assess the normative influence of occupational category on heroism perception, we compared a model that included occupation as a moderator to one that did not. The model with the occupation covariate showed a significantly better fit (RSS = 2378.1) than the model without it (RSS = 2906.4), F(16, 1340) = 18.603, p < .001, indicating that the effects were significantly driven by normative evaluations of each occupation.

mod <-lm(Heroism ~ Risk_dummy * Help_dummy, data = Set)

mod_cov<-lm(Heroism ~ Risk_dummy * Help_dummy * Job, data = Set)

anova(mod, mod_cov)

Importantly, risk and helping motivation had independent effects on perceived heroism across both models. In our main model (excluding occupation), the risk manipulation significantly increased heroism ratings, t(1356) = 6.02, p < .001, std. b = 0.16, 95% CI [0.11, 0.21]. Altruistic motivation also had a statistically significant, albeit appreciably smaller, effect, t(1356) = 2.32, p = .021, b = 0.06, 95% CI [0.0095, 0.11]. However, their interaction was not significant, std.b = -0.11, 95% CI [-0.42, 0.21], t(1356) = -0.66, p = .508, b = -0.02, 95% CI [-0.07, 0.03].

summary(mod)

## 
## Call:
## lm(formula = Heroism ~ Risk_dummy * Help_dummy, data = Set)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3246 -0.8994  0.1006  1.1006  2.3373 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             5.0199     0.0397 126.448  < 2e-16 ***
## Risk_dummy              0.4777     0.0794   6.017 2.29e-09 ***
## Help_dummy              0.1841     0.0794   2.319   0.0205 *  
## Risk_dummy:Help_dummy  -0.1051     0.1588  -0.662   0.5082    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.464 on 1356 degrees of freedom
## Multiple R-squared:  0.03004,    Adjusted R-squared:  0.0279 
## F-statistic:    14 on 3 and 1356 DF,  p-value: 5.444e-09

report(mod)

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Risk_dummy and Help_dummy (formula: Heroism ~ Risk_dummy * Help_dummy). The
## model explains a statistically significant and weak proportion of variance (R2
## = 0.03, F(3, 1356) = 14.00, p < .001, adj. R2 = 0.03). The model's intercept,
## corresponding to Risk_dummy = 0 and Help_dummy = 0, is at 5.02 (95% CI [4.94,
## 5.10], t(1356) = 126.45, p < .001). Within this model:
## 
##   - The effect of Risk dummy is statistically significant and positive (beta =
## 0.48, 95% CI [0.32, 0.63], t(1356) = 6.02, p < .001; Std. beta = 0.16, 95% CI
## [0.11, 0.21])
##   - The effect of Help dummy is statistically significant and positive (beta =
## 0.18, 95% CI [0.03, 0.34], t(1356) = 2.32, p = 0.021; Std. beta = 0.06, 95% CI
## [9.45e-03, 0.11])
##   - The effect of Risk dummy × Help dummy is statistically non-significant and
## negative (beta = -0.11, 95% CI [-0.42, 0.21], t(1356) = -0.66, p = 0.508; Std.
## beta = -0.02, 95% CI [-0.07, 0.03])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

FOOTNOTE: about the model using job as a moderator

anova(mod_cov)

As shown in table 1, manipulations of risk (vs boredom) was effective for firefighters, psychiatrists and underwater welders; but not for nurses and police officers. On the other hand, our manipulation of altruistic motivation failed to predict heroism perception, with the only exception of underwater welders. Across all occupations, our predictors failed to interact to predict perceived heroism (all ps > .1). Table 2 displays the standardized effect sizes when comparing means between experimental conditions within each target occupations.

Nurses <- subset(Set, Set$Job == "N")
Pol <- subset(Set, Set$Job == "P")
Firef <- subset(Set, Set$Job == "F")
Psych <- subset(Set, Set$Job == "Ps")
Weld <- subset(Set, Set$Job == "W")

paste0("NURSES")

## [1] "NURSES"

report(lm(Heroism ~ Risk_dummy * Help_dummy, data = Nurses))

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Risk_dummy and Help_dummy (formula: Heroism ~ Risk_dummy * Help_dummy). The
## model explains a statistically not significant and very weak proportion of
## variance (R2 = 0.01, F(3, 267) = 1.15, p = 0.330, adj. R2 = 1.63e-03). The
## model's intercept, corresponding to Risk_dummy = 0 and Help_dummy = 0, is at
## 5.48 (95% CI [5.31, 5.66], t(267) = 62.64, p < .001). Within this model:
## 
##   - The effect of Risk dummy is statistically non-significant and positive (beta
## = 0.20, 95% CI [-0.15, 0.54], t(267) = 1.12, p = 0.264; Std. beta = 0.07, 95%
## CI [-0.05, 0.19])
##   - The effect of Help dummy is statistically non-significant and positive (beta
## = 0.05, 95% CI [-0.30, 0.39], t(267) = 0.27, p = 0.789; Std. beta = 0.02, 95%
## CI [-0.10, 0.14])
##   - The effect of Risk dummy × Help dummy is statistically non-significant and
## negative (beta = -0.51, 95% CI [-1.20, 0.18], t(267) = -1.46, p = 0.146; Std.
## beta = -0.09, 95% CI [-0.21, 0.03])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

paste0("Police")

## [1] "Police"

report(lm(Heroism ~ Risk_dummy * Help_dummy, data = Pol))

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Risk_dummy and Help_dummy (formula: Heroism ~ Risk_dummy * Help_dummy). The
## model explains a statistically not significant and very weak proportion of
## variance (R2 = 0.01, F(3, 268) = 0.94, p = 0.421, adj. R2 = -6.52e-04). The
## model's intercept, corresponding to Risk_dummy = 0 and Help_dummy = 0, is at
## 4.52 (95% CI [4.36, 4.68], t(268) = 55.82, p < .001). Within this model:
## 
##   - The effect of Risk dummy is statistically non-significant and positive (beta
## = 0.24, 95% CI [-0.08, 0.56], t(268) = 1.49, p = 0.137; Std. beta = 0.09, 95%
## CI [-0.03, 0.21])
##   - The effect of Help dummy is statistically non-significant and positive (beta
## = 0.12, 95% CI [-0.19, 0.44], t(268) = 0.77, p = 0.442; Std. beta = 0.05, 95%
## CI [-0.07, 0.17])
##   - The effect of Risk dummy × Help dummy is statistically non-significant and
## positive (beta = 0.01, 95% CI [-0.63, 0.65], t(268) = 0.03, p = 0.973; Std.
## beta = 2.07e-03, 95% CI [-0.12, 0.12])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

paste0("FIREF")

## [1] "FIREF"

report(lm(Heroism ~ Risk_dummy * Help_dummy, data = Firef))

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Risk_dummy and Help_dummy (formula: Heroism ~ Risk_dummy * Help_dummy). The
## model explains a statistically not significant and very weak proportion of
## variance (R2 = 0.02, F(3, 269) = 1.39, p = 0.246, adj. R2 = 4.30e-03). The
## model's intercept, corresponding to Risk_dummy = 0 and Help_dummy = 0, is at
## 5.94 (95% CI [5.80, 6.08], t(269) = 85.00, p < .001). Within this model:
## 
##   - The effect of Risk dummy is statistically significant and positive (beta =
## 0.29, 95% CI [0.01, 0.56], t(269) = 2.04, p = 0.042; Std. beta = 0.12, 95% CI
## [4.40e-03, 0.24])
##   - The effect of Help dummy is statistically non-significant and positive (beta
## = 5.86e-03, 95% CI [-0.27, 0.28], t(269) = 0.04, p = 0.967; Std. beta =
## 2.52e-03, 95% CI [-0.12, 0.12])
##   - The effect of Risk dummy × Help dummy is statistically non-significant and
## positive (beta = 0.02, 95% CI [-0.53, 0.57], t(269) = 0.06, p = 0.950; Std.
## beta = 3.84e-03, 95% CI [-0.12, 0.12])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

paste0("PSYCH")

## [1] "PSYCH"

report(lm(Heroism ~ Risk_dummy * Help_dummy, data = Psych))

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Risk_dummy and Help_dummy (formula: Heroism ~ Risk_dummy * Help_dummy). The
## model explains a statistically significant and weak proportion of variance (R2
## = 0.05, F(3, 266) = 4.79, p = 0.003, adj. R2 = 0.04). The model's intercept,
## corresponding to Risk_dummy = 0 and Help_dummy = 0, is at 4.47 (95% CI [4.30,
## 4.64], t(266) = 52.27, p < .001). Within this model:
## 
##   - The effect of Risk dummy is statistically significant and positive (beta =
## 0.52, 95% CI [0.19, 0.86], t(266) = 3.05, p = 0.002; Std. beta = 0.18, 95% CI
## [0.06, 0.30])
##   - The effect of Help dummy is statistically significant and positive (beta =
## 0.36, 95% CI [0.02, 0.70], t(266) = 2.11, p = 0.035; Std. beta = 0.13, 95% CI
## [7.95e-03, 0.24])
##   - The effect of Risk dummy × Help dummy is statistically non-significant and
## negative (beta = -0.29, 95% CI [-0.97, 0.38], t(266) = -0.86, p = 0.392; Std.
## beta = -0.05, 95% CI [-0.17, 0.07])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

paste0("WELDERS")

## [1] "WELDERS"

report(lm(Heroism ~ Risk_dummy * Help_dummy, data = Weld))

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Risk_dummy and Help_dummy (formula: Heroism ~ Risk_dummy * Help_dummy). The
## model explains a statistically significant and moderate proportion of variance
## (R2 = 0.19, F(3, 270) = 20.65, p < .001, adj. R2 = 0.18). The model's
## intercept, corresponding to Risk_dummy = 0 and Help_dummy = 0, is at 4.68 (95%
## CI [4.53, 4.84], t(270) = 59.29, p < .001). Within this model:
## 
##   - The effect of Risk dummy is statistically significant and positive (beta =
## 1.18, 95% CI [0.87, 1.49], t(270) = 7.46, p < .001; Std. beta = 0.41, 95% CI
## [0.30, 0.52])
##   - The effect of Help dummy is statistically significant and positive (beta =
## 0.40, 95% CI [0.08, 0.71], t(270) = 2.50, p = 0.013; Std. beta = 0.14, 95% CI
## [0.03, 0.25])
##   - The effect of Risk dummy × Help dummy is statistically non-significant and
## positive (beta = 0.18, 95% CI [-0.44, 0.80], t(270) = 0.57, p = 0.568; Std.
## beta = 0.03, 95% CI [-0.08, 0.14])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

# Regression results data frame
regression_results <- data.frame(
  Group    = rep(c("Nurses", "Police", "Firefighters", "Psychiatrists", "Welders"), each = 3),
  Effect   = rep(c("Risk dummy", "Help dummy", "Risk dummy × Help dummy"), times = 5),
  Beta     = c(0.20, 0.05, -0.51,
               0.24, 0.12, 0.01,
               0.29, 0.00586, 0.02,
               0.52, 0.36, -0.29,
               1.18, 0.40, 0.18),
  CI       = c("[-0.15, 0.54]", "[-0.30, 0.39]", "[-1.20, 0.18]",
               "[-0.08, 0.56]", "[-0.19, 0.44]", "[-0.63, 0.65]",
               "[0.01, 0.56]", "[-0.27, 0.28]", "[-0.53, 0.57]",
               "[0.19, 0.86]", "[0.02, 0.70]", "[-0.97, 0.38]",
               "[0.87, 1.49]", "[0.08, 0.71]", "[-0.44, 0.80]"),
  t_value  = c(1.12, 0.27, -1.46,
               1.49, 0.77, 0.03,
               2.04, 0.04, 0.06,
               3.05, 2.11, -0.86,
               7.46, 2.50, 0.57),
  p_value  = c(0.264, 0.789, 0.146,
               0.137, 0.442, 0.973,
               0.042, 0.967, 0.950,
               0.002, 0.035, 0.392,
               "<0.001", 0.013, 0.568),
  Std_Beta = c(0.07, 0.02, -0.09,
               0.09, 0.05, 0.00207,
               0.12, 0.00252, 0.00384,
               0.18, 0.13, -0.05,
               0.41, 0.14, 0.03),
  Std_CI   = c("[-0.05, 0.19]", "[-0.10, 0.14]", "[-0.21, 0.03]",
               "[-0.03, 0.21]", "[-0.07, 0.17]", "[-0.12, 0.12]",
               "[4.40e-03, 0.24]", "[-0.12, 0.12]", "[-0.12, 0.12]",
               "[0.06, 0.30]", "[7.95e-03, 0.24]", "[-0.17, 0.07]",
               "[0.30, 0.52]", "[0.03, 0.25]", "[-0.08, 0.14]"),
  stringsAsFactors = FALSE
)

# Create the effect sizes data frame
effect_sizes <- data.frame(
  Group    = rep(c("Nurses", "Police", "Firefighters", "Psychiatrists", "Welders"), each = 2),
  Effect   = rep(c("Risk", "Motivation"), times = 5),
  Cohen_d  = c(0.14, 0.03,
               0.18, 0.09,
               0.25, 0.00593,
               0.37, 0.25,
               0.89, 0.27),
  CI       = c("[-0.10, 0.37]", "[-0.20, 0.27]",
               "[-0.06, 0.42]", "[-0.14, 0.33]",
               "[0.01, 0.49]", "[-0.23, 0.24]",
               "[0.13, 0.61]", "[0.01, 0.49]",
               "[0.64, 1.14]", "[0.03, 0.51]"),
  stringsAsFactors = FALSE
)

regression_results_gt <- regression_results %>%
  gt(groupname_col = "Group") %>%  # automatically groups rows by the Group column
  fmt_number(
    columns = c("Beta", "t_value", "Std_Beta"),
    decimals = 2
  ) %>%
  tab_header(title = "Regression Results Summary")



## Create the gt table for effect sizes:
effect_sizes_gt <- effect_sizes %>%
  gt(groupname_col = "Group") %>%
  # Format numeric columns to display 2 decimals
  fmt_number(
    columns = c("Cohen_d"), # add any other numeric columns you want to format
    decimals = 2
  ) %>%
  tab_header(
    title = "Effect Sizes (Cohen's d)"
  ) %>%
  cols_label(
    Effect  = "Effect",
    Cohen_d = "Cohen's d",
    CI      = "95% CI"
  )

# Display the tables in an HTML document (in R Markdown, simply putting the table object in a code chunk will render it)
regression_results_gt

Regression Results Summary
Effect	Beta	CI	t_value	p_value	Std_Beta	Std_CI
Nurses
Risk dummy	0.20	[-0.15, 0.54]	1.12	0.264	0.07	[-0.05, 0.19]
Help dummy	0.05	[-0.30, 0.39]	0.27	0.789	0.02	[-0.10, 0.14]
Risk dummy × Help dummy	−0.51	[-1.20, 0.18]	−1.46	0.146	−0.09	[-0.21, 0.03]
Police
Risk dummy	0.24	[-0.08, 0.56]	1.49	0.137	0.09	[-0.03, 0.21]
Help dummy	0.12	[-0.19, 0.44]	0.77	0.442	0.05	[-0.07, 0.17]
Risk dummy × Help dummy	0.01	[-0.63, 0.65]	0.03	0.973	0.00	[-0.12, 0.12]
Firefighters
Risk dummy	0.29	[0.01, 0.56]	2.04	0.042	0.12	[4.40e-03, 0.24]
Help dummy	0.01	[-0.27, 0.28]	0.04	0.967	0.00	[-0.12, 0.12]
Risk dummy × Help dummy	0.02	[-0.53, 0.57]	0.06	0.95	0.00	[-0.12, 0.12]
Psychiatrists
Risk dummy	0.52	[0.19, 0.86]	3.05	0.002	0.18	[0.06, 0.30]
Help dummy	0.36	[0.02, 0.70]	2.11	0.035	0.13	[7.95e-03, 0.24]
Risk dummy × Help dummy	−0.29	[-0.97, 0.38]	−0.86	0.392	−0.05	[-0.17, 0.07]
Welders
Risk dummy	1.18	[0.87, 1.49]	7.46	<0.001	0.41	[0.30, 0.52]
Help dummy	0.40	[0.08, 0.71]	2.50	0.013	0.14	[0.03, 0.25]
Risk dummy × Help dummy	0.18	[-0.44, 0.80]	0.57	0.568	0.03	[-0.08, 0.14]

As predicted, both the character attributions (perceived bravery and perceived selflessness) and the occupations’ descriptors (perceived danger and perceived helpfulness) were strong predictors of perceived heroism (see correlation matrix, Table 3).

Character attribution main model

Set$Selfless_scale <- scale(Set$Selfless)
Set$Brave_scale <- scale(Set$Brave)

Set$Danger_scale <- scale(Set$Danger)
Set$Helpful_scale <- scale(Set$Helpfulness)
Set$Attitude_scale <- scale(Set$Attitude)

summary(mod<-lm(Heroism ~ Selfless_scale * Brave_scale , data = Set))

## 
## Call:
## lm(formula = Heroism ~ Selfless_scale * Brave_scale, data = Set)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.2958 -0.4015  0.2531  0.7042  4.1631 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 4.88313    0.03448  141.62   <2e-16 ***
## Selfless_scale              0.57453    0.04157   13.82   <2e-16 ***
## Brave_scale                 0.58440    0.04483   13.04   <2e-16 ***
## Selfless_scale:Brave_scale  0.20734    0.02281    9.09   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.141 on 1356 degrees of freedom
## Multiple R-squared:  0.4105, Adjusted R-squared:  0.4092 
## F-statistic: 314.8 on 3 and 1356 DF,  p-value: < 2.2e-16

Situation attribution model

Set$Danger_scale <- scale(Set$Danger)
Set$Helpfulness_scale <- scale(Set$Helpfulness)

summary(mod<-lm(Heroism ~ Helpfulness_scale * Danger_scale , data = Set))

## 
## Call:
## lm(formula = Heroism ~ Helpfulness_scale * Danger_scale, data = Set)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0765 -0.5704  0.1585  0.8965  3.6948 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     5.00226    0.03344 149.575   <2e-16 ***
## Helpfulness_scale               0.74007    0.03610  20.502   <2e-16 ***
## Danger_scale                    0.35615    0.03533  10.082   <2e-16 ***
## Helpfulness_scale:Danger_scale  0.04330    0.02647   1.636    0.102    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.156 on 1356 degrees of freedom
## Multiple R-squared:  0.3952, Adjusted R-squared:  0.3938 
## F-statistic: 295.3 on 3 and 1356 DF,  p-value: < 2.2e-16

PerformanceAnalytics::chart.Correlation(Set[, c("Heroism", "Selfless", "Brave", "Danger", "Helpfulness", "Attitude")], method = "spearman")

## Additional Unreported analyses

Diagnostics

Below, we present some plots showing QQ-plot (normality of the residuals), fitted vs predicted residuals (heteroscedasticity) and Cook’s distance of predicted values. These diagnostics flags major concerns about normality of the residuals and the influence of some observations (Cook’s d plot). This further warrants some robust analyses (see the robust model). In the mean time, we also present the OLS results.

mod <-lm(Heroism ~ Risk_dummy * Help_dummy, data = Set)

plot(mod)

Outliers analyses and robustness checks

Note on Outlier management: As registered, We compared the outputs from our models to the outputs of a robust model using a smoothed Huber function to down-weight extreme residuals. Robust models used the default parameters from the lmrob command in the robustbase package – that is a set of parameters resulting in a 95% efficiency (assuming normal residual distribution and no contamination). This model is robust to the presence of extreme values and less sensitive to deviations from normality assumption

Cook’s distance plot does show that there are some influential cases:

ols_plot_cooksd_bar(mod, type = 1)

(Note that in all Cook’s d barplots, the threshold for identifying influential cases is 4/N = 0.003 – this is an usual way of finding influential cases based on Cook’s distance but… it’s a rule of thumb).

This warrant the use of robust models, or at the very least, a model comparison with robust models.

summary(modrob<-lmrob(Set$Heroism ~ Set$Risk_dummy * Set$Help_dummy))

## 
## Call:
## lmrob(formula = Set$Heroism ~ Set$Risk_dummy * Set$Help_dummy)
##  \--> method = "MM"
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.40872 -0.96122  0.03878  1.03878  2.25711 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    5.10065    0.04178 122.071  < 2e-16 ***
## Set$Risk_dummy                 0.49720    0.08155   6.097 1.41e-09 ***
## Set$Help_dummy                 0.16863    0.08137   2.072   0.0384 *  
## Set$Risk_dummy:Set$Help_dummy -0.09939    0.16268  -0.611   0.5413    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 1.462 
## Multiple R-squared:  0.03149,    Adjusted R-squared:  0.02935 
## Convergence in 10 IRWLS iterations
## 
## Robustness weights: 
##  93 weights are ~= 1. The remaining 1267 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3432  0.8792  0.9545  0.9088  0.9852  0.9972 
## Algorithmic parameters: 
##        tuning.chi                bb        tuning.psi        refine.tol 
##         1.548e+00         5.000e-01         4.685e+00         1.000e-07 
##           rel.tol         scale.tol         solve.tol          zero.tol 
##         1.000e-07         1.000e-10         1.000e-07         1.000e-10 
##       eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
##         7.353e-05         1.819e-12         5.000e-01         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##            500             50              2              1            200 
##    maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
##            200              0           1000              0           2000 
##                   psi           subsampling                   cov 
##            "bisquare"         "nonsingular"         ".vcov.avar1" 
## compute.outlier.stats 
##                  "SM" 
## seed : int(0)

paste0("Weights applied to residuals - a value of zero would mean that the observation was discarded, a value of 1 means no re-weighting")

## [1] "Weights applied to residuals - a value of zero would mean that the observation was discarded, a value of 1 means no re-weighting"

plot(modrob$rweights)

==> It does not influence our inferences.

Attitude as a covariate (explaining Attitude-Free Heroism)

When accounting for attitude, the effect of Helpfulness became non significant. In other words, the part of variance of heroism that was successfully influenced by our helpfulness manipulation was shared with attitude.

Moreover, attitude accounted for Heroism ratings above and beyond our manipulation.

summary(mod_cov2<-lm(Heroism ~ Risk_dummy * Help_dummy + Attitude, data = Set))

## 
## Call:
## lm(formula = Heroism ~ Risk_dummy * Help_dummy + Attitude, data = Set)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1585 -0.4928 -0.0421  0.8415  3.6071 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.196735   0.148368   1.326    0.185    
## Risk_dummy             0.276350   0.059317   4.659 3.49e-06 ***
## Help_dummy            -0.008258   0.059290  -0.139    0.889    
## Attitude               0.831791   0.025076  33.171  < 2e-16 ***
## Risk_dummy:Help_dummy  0.012139   0.118064   0.103    0.918    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.088 on 1355 degrees of freedom
## Multiple R-squared:  0.4647, Adjusted R-squared:  0.4631 
## F-statistic: 294.1 on 4 and 1355 DF,  p-value: < 2.2e-16

report(mod_cov2)

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Risk_dummy, Help_dummy and Attitude (formula: Heroism ~ Risk_dummy * Help_dummy
## + Attitude). The model explains a statistically significant and substantial
## proportion of variance (R2 = 0.46, F(4, 1355) = 294.08, p < .001, adj. R2 =
## 0.46). The model's intercept, corresponding to Risk_dummy = 0, Help_dummy = 0
## and Attitude = 0, is at 0.20 (95% CI [-0.09, 0.49], t(1355) = 1.33, p = 0.185).
## Within this model:
## 
##   - The effect of Risk dummy is statistically significant and positive (beta =
## 0.28, 95% CI [0.16, 0.39], t(1355) = 4.66, p < .001; Std. beta = 0.09, 95% CI
## [0.05, 0.13])
##   - The effect of Help dummy is statistically non-significant and negative (beta
## = -8.26e-03, 95% CI [-0.12, 0.11], t(1355) = -0.14, p = 0.889; Std. beta =
## -2.77e-03, 95% CI [-0.04, 0.04])
##   - The effect of Attitude is statistically significant and positive (beta =
## 0.83, 95% CI [0.78, 0.88], t(1355) = 33.17, p < .001; Std. beta = 0.67, 95% CI
## [0.63, 0.71])
##   - The effect of Risk dummy × Help dummy is statistically non-significant and
## positive (beta = 0.01, 95% CI [-0.22, 0.24], t(1355) = 0.10, p = 0.918; Std.
## beta = 2.05e-03, 95% CI [-0.04, 0.04])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

Discussion

Study 1 was a relative success. We successfully identified the predicted correlational structure between perceived bravery, selflessness, danger exposure, helpfulness and perceived heroism. Our manipulations of risk and altruism successfully predicted heroism overall. In particular, the risk manipulation produced a medium effect size on Heroism ratings (average Cohen’s d = 0.33, 95%CI [0.22, 0.33]) – indicating that boredom vs risk exposure was a strong manipulation. Notably, firefighters, typically heroized group of workers – for whom being exposed to risk is a given – saw their heroism level decrease as a consequence of our boredom description. Nevertheless, regarding altruism, our manipulation, albeit successful in influencing heroism ratings, was associated with a small overall effect size (Cohen’s d = 0.12, 95% CI [0.02, 0.23], see Table 2 for effect sizes within each occupation). Considering our manipulation of altruism, it appeared retrospectively that the two conditions were not mutually exclusive. That is, “Most X are motivated to help people” and “Most X are motivated by self-improvement” are not inconsistent motives as one can be motivated to both help people and self-improve. A stronger manipulation of motivation would create a contrast between the two conditions where the alternative condition is negated (e.g., explaining that “Most X are motivated to help people, but not self-improve” and vice-versa). In addition, it is possible, and suggested in a model directly confronting perceived bravery, selflessness, risk exposure and helpfulness (see Exploratory analyses in our SOM), that heroism ratings are more about the occupation effectively helping people, rather than a character attribution of altruism. Perhaps a group would be more heroic in helping people – despite not being driven to do so – rather than being ineffective in helping people – despite a strong desire to help. Study 2 addressed directly the centrality of altruistic motivation vs altruistic outcomes by manipulating them orthogonally.

Study 2

We revised our manipulation to include a framing of the occupations as being either factually helpful or factually not very helpful. In study 2, we manipulated the presentation of altruistic outcome (i.e., being actually helpful) in combination with our previous manipulation of altruistic motivation (i.e., being selfless). In particular, we asked participants to consider that target occupations are helpful in addition to being selfless, or helpful despite not being selfless, or unhelpful despite being selfless, or finally unhelpful in addition to not being selfless. Inferring intentions from outcomes is well-documented in moral psychology (e.g., Baron & Hershey, 1988; Knobe, 2003), but not the reverse. It is thus possible that framing occupations by their observable, altruistic outcomes might be more effective in changing the perception of the occupation, rather than trying to manipulate the perception of altruistic motivations. Study 2 directly tested whether heroism perception is about altruistic motivation (i.e., wanting to help), or altruistic outcome (i.e., actually helping). To our knowledge, no prior studies have directly compared the effects of framing altruistic intentions versus altruistic outcomes on collective perceptions. We registered this study on OSF (https://osf.io/96awk).

Protocol

We registered a sample of 816 participants UK representative participants from Prolific. Two participants failed the comprehension checks, three participants returned their survey despite completing it, and two participants were timed-out by Prolific despite finishing the survey – we collected seven additional participants to compensate for this attrition, resulting in a complete representative sample of 816 UK residents (420 women, 394 men, and 2 others, Mage = 46.39, SD = 15.4). Participants received £9 per hour.

The protocol of study 2 was similar to Study 1 with the following changes: only the occupations that were sensitive to our previously tested exposure to danger manipulation were tested (i.e., firefighters, underwater welders, and psychiatrists). Each occupation was framed as exposed to danger (high danger condition from study 2). However, we manipulated frames of target occupations as having altruistic motivations (i.e., 82% of XXX reported helping others as their motivation vs. Only 8% of XXX reported helping others as their motivation, see Material online for exact scenarios) and/or altruistic outcomes (i.e., effectively being helpful, ‘89 % of the missions of the workers effectively help people’ vs ‘Only 11% of the missions of the workers effectively help people’). The two manipulated elements were presented to participants in a random order (see SOM for the specific vignettes used). In addition, we shortened the comprehension check to only include one question that participants had two chances of getting correct. As in the previous studies, after reading the descriptions of the occupations, participants reported perceived heroism of the target occupation. They then reported perceived bravery, perceived selflessness (manipulation check of the altruistic motivation manipulation), perceived exposure to danger, and perceived helpfulness (manipulation check of the altruistic outcome manipulation) of the target occupations. Finally, they completed the attitude measure, completed demographics, were debriefed and thanked.

Results

Set <- read.csv("Raw_Data_GW1_ter.csv", comment.char="#")
Set <- subset(Set, Set$Q233 != "No, I won't")
Set <- Set[-c(1:2),]
Set <- subset(Set, Set$FailedComp == "No") # Keep only participants who succeeded the completion checks
Set <- subset(Set, Set$Age != "") # Keep only participants who succeeded the completion checks

# Remove this weird participant that completed the survey yet returned the survey to preserve representativity
ProlExport<- read.csv("~/Downloads/prolific_export_68541b7296d11f4b2e3bc37e(2).csv")
Ret <- subset(ProlExport, ProlExport$Status == "RETURNED")
# intersect(Set$Prol_ID, Ret$Participant.id)

# Remove Timed Out yet Complete Questionnaires to preserve representativity
ProlExport<- read.csv("Prolific_Demographic_GW1_ter.csv")
TO <- subset(ProlExport, ProlExport$Status == "TIMED-OUT")
# intersect(Set$Prol_ID, TO$Participant.id)

# Also, one participant reported a mistake, we can correct this:
idx <- which(
  Set$Q523 == "I made a mistake in the first question when asked about whether firefighters are heroes and I said I strongly disagree while I meant to say strongly agree."
)

Set$F_Selfish_Help_h_1[idx] <- "7 - Strongly agree"


Set <- subset(Set, 
              !(Set$PROLIFIC_PID %in% TO$Participant.id | 
                Set$PROLIFIC_PID %in% Ret$Participant.id))


#table(Set$SubCond, Set$Condition)
#colnames(Set)
Demographics <- Set[,c(37, 39:46)]# First, add a unique row identifier if one does not exist already:
job_cols <- c("Job_match_1", "Job_match_2", "Job_match_3", 
              "Job_match_4", "Job_match_5", "Job_match_7", "Job_match_6")

Demographics$JOB <- apply(Demographics[, job_cols], 1, function(x) {
  # Remove empty strings
  reported_jobs <- x[x != ""]
  # If no jobs reported, assign NA; else paste the jobs together
  if (length(reported_jobs) == 0) NA else paste(reported_jobs, collapse = ", ")
})

Set$Credibility[Set$Credibility == "Very unbelievable"] <- 1
Set$Credibility[Set$Credibility == "Quite unbelievable"] <- 2
Set$Credibility[Set$Credibility == "Somewhat unbelievable"] <- 3
Set$Credibility[Set$Credibility == "Neutral"] <- 4
Set$Credibility[Set$Credibility == "Somewhat believable"] <- 5
Set$Credibility[Set$Credibility == "Quite believable"] <- 6
Set$Credibility[Set$Credibility == "Very believable"] <- 7



# 1. Identify all columns whose names start with "AC" (attention checks)
ac_cols <- grep("^AC", names(Set), value = TRUE)

# 2. Pivot those AC columns into long form so each respondent has exactly one row
Set_long <- Set %>%
  mutate(row_id = row_number()) %>%        # create a unique ID for each row
  pivot_longer(
    cols      = all_of(ac_cols),           # take all attention-check columns
    names_to  = "ac_var",                  # put their names into a column called ac_var
    values_to = "ac_resp"                  # put their values into a column called ac_resp
  ) %>%
  filter(ac_resp != "")                    # keep only the non-empty entry per respondent



# Identify all columns ending with "h_1" and use it to fill the column "CONDITION"
h1_cols <- grep("_h_1$", names(Set_long), value = TRUE)
# Creating the condition
Set_long_filtered <- Set_long %>%
  rowwise() %>%  # process each row individually
  mutate(
    # Count how many of the h1_cols are non-empty ("")
    non_empty_count = sum(c_across(all_of(h1_cols)) != ""),
    
    # Determine the Condition based on non-empty count
    Condition_2 = if (non_empty_count == 1) {
      # Get the name of the column that is non-empty
      selected_col <- h1_cols[which(c_across(all_of(h1_cols)) != "")]
      # Remove the "_h_1" suffix to get the condition (e.g., "NRS" from "NRS_h_1")
      str_remove(selected_col, "_h_1")
    } else {
      "Error"
    }
  ) %>%
  ungroup() %>%  # exit rowwise mode
  dplyr::select(-non_empty_count) %>% distinct(ResponseId, .keep_all=T) # remove the helper column if no longer needed



### RESPONSES TO ATTENTION CHECKS MATCH THE CONDITION

# Now I'll neatly create the main data frame. I create a helper function that extract all non empty values and apply it on all rows ending with our codes: 
# - _h_1 is heroism rating, _m_1_1 is first item of MC1, m_1_2, is second item of mc2, _at is attitude etc.


# Helper function: Given a vector, return the only non-empty element (or NA if none/few found)
extract_value <- function(x) {
  non_empty <- x[x != ""]
  if(length(non_empty) == 1) return(non_empty)
  else return(NA)
}

# Create the final data frame
final_df <- Set_long_filtered %>%
  mutate(
    # Extract the non-empty value from each group of columns:
    Heroism = apply(dplyr::select(., ends_with("_h_1")), 1, extract_value),
    Selfless        = apply(dplyr::select(., ends_with("m_1")), 1, extract_value),
    Helpfulness        = apply(dplyr::select(., ends_with("m_2")), 1, extract_value),
    Brave        = apply(dplyr::select(., ends_with("m_3")), 1, extract_value),
    Danger        = apply(dplyr::select(., ends_with("m_4")), 1, extract_value),
    Attitude     = apply(dplyr::select(., matches("(_at)$")), 1, extract_value)
  ) %>%
  # Select only the required columns in the final data frame
  dplyr::select(ResponseId, Condition, Condition_2, Heroism, Danger, Helpfulness, Brave, Selfless, Attitude, Credibility, Gender, Age)

# Recoding values:
final_df<-subset(final_df, final_df$Condition_2 != "Error")
final_df <- final_df %>%
  # 1. Split Condition into Job, Motivation, and Help
  mutate(
    Job  = str_sub(Condition_2, 1, 1),
    Motivation = if_else(grepl("Selfless", Condition_2, ignore.case = TRUE), "Selfless", "Selfish"),
    Help       = if_else(grepl("Unhelp",     Condition_2, ignore.case = TRUE), "Unhelp", "Help")  ) %>%
  # 2. Dummy code Motivation and Help: R -> 0.5, B -> -0.5; H -> 0.5, S -> -0.5
  mutate(
    Motivation_dummy = case_when(
      Motivation == "Selfless" ~ 0.5,
      Motivation == "Selfish" ~ -0.5,
      TRUE ~ NA_real_
    ),
    Help_dummy = case_when(
      Help == "Help" ~ 0.5,
      Help == "Unhelp" ~ -0.5,
      TRUE ~ NA_real_
    )
  ) %>%
  # 3. Recode rating items (HeroismScore, MC1_1, MC1_2, MC2_1, MC2_2)
  mutate_at(vars(Heroism, Danger, Helpfulness, Brave, Selfless), 
            ~ case_when(
              . == "1 - Strongly disagree" ~ "1",
              . == "7 - Strongly agree"    ~ "7",
              TRUE ~ .
            )) %>%
  # Optionally convert these recoded columns to numeric
  mutate_at(vars(Heroism, Danger, Helpfulness,Brave, Selfless), as.numeric) %>%
  # 4. Recode Attitude values
  mutate(
    Attitude = case_when(
      Attitude == "Very negative"    ~ 1,
      Attitude == "Quite negative" ~ 2,
      Attitude == "Somewhat negative" ~ 3,
      Attitude == "Neutral"           ~ 4,
      Attitude == "Somewhat positive" ~ 5,
      Attitude == "Quite positive" ~ 6,
      Attitude == "Very positive" ~ 7,
      TRUE ~ NA_real_
    )
  ) 

#final_df$Part_Job <- Demographics$JOB
Set<-final_df

Set$Job[Set$Job == "F"] <- "Firefighters"
Set$Job[Set$Job == "W"] <- "Welders"
Set$Job[Set$Job == "P"] <- "Psychiatrists"

Demographics

Demographics$flag_inconsistent <- (Demographics$Job_match_6 == "None of the above") & 
  (Demographics$Job_match_1 != "" | Demographics$Job_match_2 != "" | Demographics$Job_match_3 != "" | 
   Demographics$Job_match_4 != "" | Demographics$Job_match_5 != "")




paste0("816 participants took part in the study. Mean age in the sample is ", mean(as.numeric(Set$Age)), ", SD = ", sd(as.numeric(Set$Age)))

## [1] "816 participants took part in the study. Mean age in the sample is 46.389705882353, SD = 15.4076737127492"

## Gender

Set %>% group_by(Gender) %>% summarise(N=n()) %>%
  ggplot(aes(x=Gender,y=N,fill=Gender))+
  geom_bar(stat = 'identity',color='black')+
  scale_y_continuous(labels = scales::comma_format(accuracy = 2))+
  geom_text(aes(label=N),vjust=-0.25,fontface='bold')+
  theme_bw()+
  theme(axis.text = element_text(color='black',face='bold'),
        axis.title = element_text(color='black',face='bold'),
        legend.text = element_text(color='black',face='bold'),
        legend.title = element_text(color='black',face='bold')) +
  ggtitle("Gender distribution")

## Occupations
#colnames(Set)
#jobs <- unlist(Demographics[-which(Demographics$flag_inconsistent == T), 4:10])           # Make a long list of all jobs that were named

jobs <- unlist(Demographics[, 3:9])           # Make a long list of all jobs that were named


jobs <- jobs[jobs != ""]     # Remove empty strings

job_df <- as.data.frame(table(jobs))
colnames(job_df) <- c("Job", "Count")

ggplot(job_df, aes(x = Job, y = Count, fill = Job)) +
  geom_bar(stat = 'identity',color='black')+
  scale_y_continuous(labels = scales::comma_format(accuracy = 2))+
  geom_text(aes(label=Count),vjust=-0.25,fontface='bold')+
  theme_bw()+
  theme(axis.text = element_text(color='black',face='bold'),
        axis.title = element_text(color='black',face='bold'),
        legend.text = element_text(color='black',face='bold'),
        legend.title = element_text(color='black',face='bold')) +
  ggtitle("Job distribution")

Descriptives

# 1. Create a summary dataframe for facet annotations (mean and SD)
df_summary <- Set %>%
  group_by(Job) %>%
  summarize(
    mean_score = mean(Heroism, na.rm = TRUE),
    sd_score   = sd(Heroism, na.rm = TRUE),
    .groups = "drop"
  )

# 2. Create the ggplot using the long format data
ggplot(Set, aes(x = Heroism)) +
  geom_histogram(aes(fill = after_stat(count)),
                 binwidth = 1,
                 color = "black", show.legend = FALSE) +
  facet_grid( ~ Job, scales = "free") +
  scale_fill_gradientn(
    colours = brewer.pal(9, "YlOrBr"),
    name = "Count"
  ) +
  labs(
    title = "Histograms of Heroism by Occupation",
    x = "Score",
    y = "Count"
  ) +
  # Annotate each facet with the mean and standard deviation
  geom_text(data = df_summary,
            aes(x = 7, y = Inf,
                label = paste0("Mean = ", round(mean_score, 2),
                               "\nSD = ", round(sd_score, 2))),
            vjust = 1.5, hjust = 1.1, size = 3) +
  theme_classic() +
  theme(panel.grid.major.y = element_line(linewidth = 0.5),
        panel.grid.minor.y = element_line(linewidth = 0.5))

Set$Credibility <- as.numeric(Set$Credibility)
mean(Set$Credibility)

## [1] 5.362745

sd(Set$Credibility)

## [1] 1.395271

df_long <- Set %>%
  pivot_longer(
    cols = c(Heroism, Danger, Helpfulness, Attitude, Credibility),
    names_to = "Variable",
    values_to = "Score"
  )

# 2. Compute summary statistics by Job and Variable
df_long2<- subset(df_long, df_long$Variable != "Credibility")

df_summary <- df_long %>%
  group_by(Job, Variable) %>%
  summarize(
    mean_score = mean(Score, na.rm = TRUE),
    sd_score   = sd(Score, na.rm = TRUE),
    .groups = "drop"
  )
df_summary2 <- df_long2 %>%
  group_by(Job, Variable) %>%
  summarize(
    mean_score = mean(Score, na.rm = TRUE),
    sd_score   = sd(Score, na.rm = TRUE),
    .groups = "drop"
  )
# 3. Create the ggplot using the long format data
ggplot(df_long2, aes(x = Score)) +
  geom_histogram(aes(fill = after_stat(count)),
                 binwidth = 1,
                 color = "black", show.legend = FALSE) +
  facet_grid(Variable ~ Job, scales = "free") +
  scale_fill_gradientn(
    colours = brewer.pal(9, "YlOrBr"),
    name = "Count"
  ) +
  labs(
    title = "Histograms of Variable by Occupation",
    x = "Score",
    y = "Count"
  ) +
  # Annotate each facet with the mean and standard deviation
  geom_text(data = df_summary2,
            aes(x = 7, y = Inf,
                label = paste0("Mean = ", round(mean_score, 2),
                               "\nSD = ", round(sd_score, 2))),
            vjust = 1.5, hjust = 1.1, size = 3) +
  theme_classic() +
  theme(
    panel.grid.major.y = element_line(linewidth = 0.5),
    panel.grid.minor.y = element_line(linewidth = 0.5)
  )

Main reported analyses - as reported

The average credibility of our manipulations was M = 5.36, SD = 1.40, indicating a generally credible study. The effects of our experimental conditions on credibility of the manipulation are detailed online (see SOM).

Cred <- subset(df_long, df_long$Variable == "Credibility")

mean(Cred$Score)

## [1] 5.362745

sd(Cred$Score)

## [1] 1.395271

df_summary <- Cred %>%
  group_by(Job, Help, Motivation) %>%
  summarize(
    mean_score = mean(Score, na.rm = TRUE),
    sd_score   = sd(Score, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(Cred, aes(x = Score)) +
  geom_histogram(aes(fill = after_stat(count)),
                 binwidth = 1,
                 color = "black", show.legend = FALSE) +
  facet_grid(
    Help + Motivation ~ Job,
    scales = "free"
  ) +
  scale_fill_gradientn(
    colours = RColorBrewer::brewer.pal(9, "YlOrBr"),
    name = "Count"
  ) +
  labs(
    title = "Frequency of Credibility Scores by Job, Help & Motivation",
    x = "Score", y = "Count"
  ) +
  geom_text(data = df_summary,
            aes(x = 7, y = Inf,
                label = paste0("Mean = ", round(mean_score, 2),
                               "\nSD = ", round(sd_score, 2))),
            vjust = 1.5, hjust = 1.1, size = 3) +
  theme_classic() +
  theme(
    panel.grid.major.y = element_line(linewidth = 0.5),
    panel.grid.minor.y = element_line(linewidth = 0.5)
  )

Manipulation checks were assessed using Welch t-test. Results indicated that workers were perceived as significantly more selfless (M = 5.28, SD = 1.44) in the altruistic motivation condition in comparison to the non-altruistic motivation condition (M = 4.85, SD = 1.54), t(810.46) = 4.10, p < .001; Cohen’s d = 0.29, 95% CI [0.15, 0.42].

Set$Motivation <- factor(Set$Motivation, levels = c("Selfless", "Selfish"))

by(Set$Selfless, Set$Help, mean)

## Set$Help: Help
## [1] 5.293399
## ------------------------------------------------------------ 
## Set$Help: Unhelp
## [1] 4.835381

by(Set$Selfless, Set$Help, sd)

## Set$Help: Help
## [1] 1.460922
## ------------------------------------------------------------ 
## Set$Help: Unhelp
## [1] 1.526426

t.test(Set$Selfless ~Set$Motivation)

## 
##  Welch Two Sample t-test
## 
## data:  Set$Selfless by Set$Motivation
## t = 4.0986, df = 810.46, p-value = 4.575e-05
## alternative hypothesis: true difference in means between group Selfless and group Selfish is not equal to 0
## 95 percent confidence interval:
##  0.2236695 0.6348156
## sample estimates:
## mean in group Selfless  mean in group Selfish 
##               5.280098               4.850856

cohens_d(Set$Selfless ~ Set$Motivation)

my_sum <- Set %>%
  group_by(Motivation) %>%
  summarise( 
    n=n(),
    mean=mean(Selfless),
    sd=sd(Selfless)
  ) %>%
  mutate( se=sd/sqrt(n))
 
# Standard deviation
ggplot(my_sum) +
  geom_bar( aes(x=Motivation, y=mean), stat="identity", fill="forestgreen", alpha=0.5) +
  geom_errorbar( aes(x=Motivation, ymin=mean-se, ymax=mean+se), width=0.4, colour="orange", alpha=0.9, size=1.5) +
  ggtitle("Selfless ~ Motivation type; bars are SE")

To a lesser extent, workers were also perceived as more helpful in the altruistic motivation condition (M = 5.80, SD = 1.38) than in the non-altruistic motivation condition (M = 5.66, SD = 1.40), but this effect was considerably small and did not reach significance, t(813.94) = 1.47, p = 0.141; Cohen’s d = 0.10, 95% CI [-0.03, 0.24].

by(Set$Helpfulness, Set$Motivation, mean)

## Set$Motivation: Selfless
## [1] 5.80344
## ------------------------------------------------------------ 
## Set$Motivation: Selfish
## [1] 5.660147

by(Set$Helpfulness, Set$Motivation, sd)

## Set$Motivation: Selfless
## [1] 1.380971
## ------------------------------------------------------------ 
## Set$Motivation: Selfish
## [1] 1.40007

t.test(Set$Helpfulness ~Set$Motivation)

## 
##  Welch Two Sample t-test
## 
## data:  Set$Helpfulness by Set$Motivation
## t = 1.4718, df = 813.94, p-value = 0.1415
## alternative hypothesis: true difference in means between group Selfless and group Selfish is not equal to 0
## 95 percent confidence interval:
##  -0.04780704  0.33439325
## sample estimates:
## mean in group Selfless  mean in group Selfish 
##               5.803440               5.660147

cohens_d(Set$Helpfulness ~Set$Motivation)

Similarly, framing the workers as effectively helping people increased perceived helpfulness perception (M = 6.36, SD = 0.85) in comparison to framing the workers as not being actually helping people (M = 5.10, SD = 1.54), t(810.46) = 4.10, p < .001; Cohen’s d = 0.29, 95% CI [0.15, 0.42].

by(Set$Helpfulness, Set$Help, mean)

## Set$Help: Help
## [1] 6.361858
## ------------------------------------------------------------ 
## Set$Help: Unhelp
## [1] 5.09828

by(Set$Helpfulness, Set$Help, sd)

## Set$Help: Help
## [1] 0.8495169
## ------------------------------------------------------------ 
## Set$Help: Unhelp
## [1] 1.536156

t.test(Set$Helpfulness ~Set$Help)

## 
##  Welch Two Sample t-test
## 
## data:  Set$Helpfulness by Set$Help
## t = 14.53, df = 632.43, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Help and group Unhelp is not equal to 0
## 95 percent confidence interval:
##  1.092808 1.434348
## sample estimates:
##   mean in group Help mean in group Unhelp 
##             6.361858             5.098280

cohens_d(t.test(Set$Helpfulness ~Set$Help))

my_sum <- Set %>%
  group_by(Help) %>%
  summarise( 
    n=n(),
    mean=mean(Helpfulness),
    sd=sd(Helpfulness)
  ) %>%
  mutate( se=sd/sqrt(n))
 
# Standard deviation
ggplot(my_sum) +
  geom_bar( aes(x=Help, y=mean), stat="identity", fill="forestgreen", alpha=0.5) +
  geom_errorbar( aes(x=Help, ymin=mean-se, ymax=mean+se), width=0.4, colour="orange", alpha=0.9, size=1.5) +
  ggtitle("Helpfulness ~ Motivation type (Helpful vs Unhelpful); bars are SE")

Our manipulation of helpfulness also lead participants to evaluate workers as more selfless in the helpful condition (M = 5.29, SD = 1.46) in comparison to the non-helpful condition (M = 4.84, SD = 1.53), t(812.07) = 4.38, p < .001; Cohen’s d = 0.31, 95%CI [0.17, 0.44].

by(Set$Selfless, Set$Help, mean)

## Set$Help: Help
## [1] 5.293399
## ------------------------------------------------------------ 
## Set$Help: Unhelp
## [1] 4.835381

by(Set$Selfless, Set$Help, sd)

## Set$Help: Help
## [1] 1.460922
## ------------------------------------------------------------ 
## Set$Help: Unhelp
## [1] 1.526426

t.test(Set$Selfless ~Set$Help)

## 
##  Welch Two Sample t-test
## 
## data:  Set$Selfless by Set$Help
## t = 4.3784, df = 812.07, p-value = 1.352e-05
## alternative hypothesis: true difference in means between group Help and group Unhelp is not equal to 0
## 95 percent confidence interval:
##  0.2526814 0.6633540
## sample estimates:
##   mean in group Help mean in group Unhelp 
##             5.293399             4.835381

cohens_d(t.test(Set$Selfless ~Set$Help))

As registered, and similar to previous studies, we compared two models assessing the normative effect of occupation on heroism ratings. As in previous studies, occupations significantly accounted for heroism: models including Occupation as a moderator provided a significantly better fit than model omitting occupation (RSS = 1348.8 vs RSS = 1731.8), F(8, 804) = 28.54, p < .001.

mod <-lm(Heroism ~ Help_dummy * Motivation_dummy, data = Set)

mod_cov<-lm(Heroism ~ Help_dummy* Motivation_dummy * Job, data = Set)

anova(mod, mod_cov)

In the registered model not accounting for occupations, altruistic motivations failed to predict heroism perception, t(812) = 1.87, p = 0.061; std. b= 0.06, 95% CI [-2.94e-03, 0.13]. However, our manipulation of altruistic outcomes successfully predicted perceived heroism, t(812) = 6.61, p < .001; std. b= 0.23, 95% CI [0.16, 0.29]4. The two predictors did not interact, t(812) = 1.45, p = .147; std. b = 0.05, 95% CI [-0.02, 0.12].

summary(mod)

## 
## Call:
## lm(formula = Heroism ~ Help_dummy * Motivation_dummy, data = Set)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3122 -0.8276  0.1724  1.2157  2.2157 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  5.14401    0.05112 100.618  < 2e-16 ***
## Help_dummy                   0.67613    0.10225   6.613 6.84e-11 ***
## Motivation_dummy             0.19152    0.10225   1.873   0.0614 .  
## Help_dummy:Motivation_dummy  0.29649    0.20450   1.450   0.1475    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.46 on 812 degrees of freedom
## Multiple R-squared:  0.05724,    Adjusted R-squared:  0.05376 
## F-statistic: 16.43 on 3 and 812 DF,  p-value: 2.243e-10

report(mod)

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Help_dummy and Motivation_dummy (formula: Heroism ~ Help_dummy *
## Motivation_dummy). The model explains a statistically significant and weak
## proportion of variance (R2 = 0.06, F(3, 812) = 16.43, p < .001, adj. R2 =
## 0.05). The model's intercept, corresponding to Help_dummy = 0 and
## Motivation_dummy = 0, is at 5.14 (95% CI [5.04, 5.24], t(812) = 100.62, p <
## .001). Within this model:
## 
##   - The effect of Help dummy is statistically significant and positive (beta =
## 0.68, 95% CI [0.48, 0.88], t(812) = 6.61, p < .001; Std. beta = 0.23, 95% CI
## [0.16, 0.29])
##   - The effect of Motivation dummy is statistically non-significant and positive
## (beta = 0.19, 95% CI [-9.18e-03, 0.39], t(812) = 1.87, p = 0.061; Std. beta =
## 0.06, 95% CI [-2.94e-03, 0.13])
##   - The effect of Help dummy × Motivation dummy is statistically non-significant
## and positive (beta = 0.30, 95% CI [-0.10, 0.70], t(812) = 1.45, p = 0.147; Std.
## beta = 0.05, 95% CI [-0.02, 0.12])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

As displayed in table 4, when decomposing the effects within each occupation, we observed that for firefighters, there was no effect of the altruistic motivation, nor any effect of altruistic outcome (ps > .05). However, regarding psychiatrists, whereas the altruistic motivation condition was non significant (p = .11), there was a main effect of altruistic outcomes (p = .025). Regarding underwater welders, we found no effect associated with framing altruistic motivations (p = .51) but a very large effect associated with framing altruistic outcome on heroic perception (p < .001). See table 5 for standardized effect sizes when comparing means between experimental conditions within each target occupations.

Firef <- subset(Set, Set$Job == "Firefighters") 
Psych <- subset(Set, Set$Job == "Psychiatrists")
Weld <- subset(Set, Set$Job == "Welders")
paste0("FIREFIGHTER")

## [1] "FIREFIGHTER"

report(lm(Heroism ~ Motivation_dummy * Help_dummy, data = Firef))

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Motivation_dummy and Help_dummy (formula: Heroism ~ Motivation_dummy *
## Help_dummy). The model explains a statistically not significant and weak
## proportion of variance (R2 = 0.03, F(3, 268) = 2.39, p = 0.069, adj. R2 =
## 0.02). The model's intercept, corresponding to Motivation_dummy = 0 and
## Help_dummy = 0, is at 6.02 (95% CI [5.88, 6.16], t(268) = 83.71, p < .001).
## Within this model:
## 
##   - The effect of Motivation dummy is statistically non-significant and positive
## (beta = 0.16, 95% CI [-0.13, 0.44], t(268) = 1.08, p = 0.281; Std. beta = 0.07,
## 95% CI [-0.05, 0.18])
##   - The effect of Help dummy is statistically non-significant and positive (beta
## = 0.18, 95% CI [-0.10, 0.47], t(268) = 1.27, p = 0.205; Std. beta = 0.08, 95%
## CI [-0.04, 0.20])
##   - The effect of Motivation dummy × Help dummy is statistically significant and
## positive (beta = 0.60, 95% CI [0.03, 1.17], t(268) = 2.09, p = 0.038; Std. beta
## = 0.13, 95% CI [7.25e-03, 0.25])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

paste0("PSYCH")

## [1] "PSYCH"

report(lm(Heroism ~ Motivation_dummy * Help_dummy, data = Psych))

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Motivation_dummy and Help_dummy (formula: Heroism ~ Motivation_dummy *
## Help_dummy). The model explains a statistically significant and weak proportion
## of variance (R2 = 0.03, F(3, 271) = 2.87, p = 0.037, adj. R2 = 0.02). The
## model's intercept, corresponding to Motivation_dummy = 0 and Help_dummy = 0, is
## at 4.71 (95% CI [4.55, 4.88], t(271) = 56.72, p < .001). Within this model:
## 
##   - The effect of Motivation dummy is statistically non-significant and positive
## (beta = 0.27, 95% CI [-0.06, 0.59], t(271) = 1.60, p = 0.111; Std. beta = 0.09,
## 95% CI [-0.02, 0.21])
##   - The effect of Help dummy is statistically significant and positive (beta =
## 0.38, 95% CI [0.05, 0.70], t(271) = 2.26, p = 0.025; Std. beta = 0.13, 95% CI
## [0.02, 0.25])
##   - The effect of Motivation dummy × Help dummy is statistically non-significant
## and positive (beta = 0.35, 95% CI [-0.31, 1.00], t(271) = 1.04, p = 0.300; Std.
## beta = 0.06, 95% CI [-0.06, 0.18])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

by(Psych$Heroism, Psych$Help, mean)

## Psych$Help: Help
## [1] 4.897059
## ------------------------------------------------------------ 
## Psych$Help: Unhelp
## [1] 4.52518

by(Psych$Heroism, Psych$Help, sd)

## Psych$Help: Help
## [1] 1.362348
## ------------------------------------------------------------ 
## Psych$Help: Unhelp
## [1] 1.400471

paste0("WELDERS")

## [1] "WELDERS"

report(lm(Heroism ~ Motivation_dummy * Help_dummy, data = Weld))

## We fitted a linear model (estimated using OLS) to predict Heroism with
## Motivation_dummy and Help_dummy (formula: Heroism ~ Motivation_dummy *
## Help_dummy). The model explains a statistically significant and moderate
## proportion of variance (R2 = 0.25, F(3, 265) = 29.00, p < .001, adj. R2 =
## 0.24). The model's intercept, corresponding to Motivation_dummy = 0 and
## Help_dummy = 0, is at 4.69 (95% CI [4.54, 4.85], t(265) = 58.55, p < .001).
## Within this model:
## 
##   - The effect of Motivation dummy is statistically non-significant and positive
## (beta = 0.11, 95% CI [-0.21, 0.42], t(265) = 0.67, p = 0.506; Std. beta = 0.04,
## 95% CI [-0.07, 0.14])
##   - The effect of Help dummy is statistically significant and positive (beta =
## 1.49, 95% CI [1.17, 1.80], t(265) = 9.29, p < .001; Std. beta = 0.50, 95% CI
## [0.39, 0.60])
##   - The effect of Motivation dummy × Help dummy is statistically non-significant
## and positive (beta = 7.46e-03, 95% CI [-0.62, 0.64], t(265) = 0.02, p = 0.981;
## Std. beta = 1.24e-03, 95% CI [-0.10, 0.11])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

by(Weld$Heroism, Weld$Help, mean)

## Weld$Help: Help
## [1] 5.437956
## ------------------------------------------------------------ 
## Weld$Help: Unhelp
## [1] 3.94697

by(Weld$Heroism, Weld$Help, sd)

## Weld$Help: Help
## [1] 1.253591
## ------------------------------------------------------------ 
## Weld$Help: Unhelp
## [1] 1.366527

# A list of our models
mods <- list(
  Firefighter = lm(Heroism ~ Motivation_dummy * Help_dummy, data = Firef),
  Psychiatrist = lm(Heroism ~ Motivation_dummy * Help_dummy, data = Psych),
  Welder       = lm(Heroism ~ Motivation_dummy * Help_dummy, data = Weld)
)

# tidy + bind + format
table1 <- bind_rows(
  lapply(names(mods), function(occ) {
    broom::tidy(mods[[occ]], conf.int = TRUE) %>%
      mutate(Occupation = occ)
  }),
  .id = NULL
) %>%
    filter(term != "(Intercept)") %>%        # drop intercepts
  mutate(
    CI = paste0("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"), # Format a CI string
    Beta = round(estimate, 2),
    t     = round(statistic, 2),
    p     = sprintf("%.3f", p.value)
  ) %>%
  dplyr::select(
    Occupation,
    Predictor = term,
    Beta,
    `95% CI` = CI,
    t,
    `p‐value` = p
  )


# render
kable(
  table1,
  caption = "Table 1. Regression coefficients by occupation",
  align   = "c"
)

Table 1. Regression coefficients by occupation

Occupation	Predictor	Beta	95% CI	t	p‐value
Firefighter	Motivation_dummy	0.16	[-0.13, 0.44]	1.08	0.281
Firefighter	Help_dummy	0.18	[-0.1, 0.47]	1.27	0.205
Firefighter	Motivation_dummy:Help_dummy	0.60	[0.03, 1.17]	2.09	0.038
Psychiatrist	Motivation_dummy	0.27	[-0.06, 0.59]	1.60	0.111
Psychiatrist	Help_dummy	0.38	[0.05, 0.7]	2.26	0.025
Psychiatrist	Motivation_dummy:Help_dummy	0.35	[-0.31, 1]	1.04	0.300
Welder	Motivation_dummy	0.11	[-0.21, 0.42]	0.67	0.506
Welder	Help_dummy	1.49	[1.17, 1.8]	9.29	0.000
Welder	Motivation_dummy:Help_dummy	0.01	[-0.62, 0.64]	0.02	0.981

d_F_Help <- cohens_d(Heroism ~ Help,       data = Firef)
d_F_Mot  <- cohens_d(Heroism ~ Motivation, data = Firef)
d_P_Help <- cohens_d(Heroism ~ Help,       data = Psych)
d_P_Mot  <- cohens_d(Heroism ~ Motivation, data = Psych)
d_W_Help <- cohens_d(Heroism ~ Help,       data = Weld)
d_W_Mot  <- cohens_d(Heroism ~ Motivation, data = Weld)

# assemble into a single data.frame
table2 <- bind_rows(
  data.frame(
    Occupation = "Firefighters",
    Predictor  = "Helpfulness",
    d          = round(d_F_Help$Cohens_d, 2),
    CI         = paste0("[", round(d_F_Help$CI_low, 2), ", ", round(d_F_Help$CI_high, 2), "]")
  ),
  data.frame(
    Occupation = "Firefighters",
    Predictor  = "Motivation",
    d          = round(d_F_Mot$Cohens_d, 2),
    CI         = paste0("[", round(d_F_Mot$CI_low, 2), ", ", round(d_F_Mot$CI_high, 2), "]")
  ),
  data.frame(
    Occupation = "Psychiatrists",
    Predictor  = "Helpfulness",
    d          = round(d_P_Help$Cohens_d, 2),
    CI         = paste0("[", round(d_P_Help$CI_low, 2), ", ", round(d_P_Help$CI_high, 2), "]")
  ),
  data.frame(
    Occupation = "Psychiatrists",
    Predictor  = "Motivation",
    d          = round(d_P_Mot$Cohens_d, 2),
    CI         = paste0("[", round(d_P_Mot$CI_low, 2), ", ", round(d_P_Mot$CI_high, 2), "]")
  ),
  data.frame(
    Occupation = "Welders",
    Predictor  = "Helpfulness",
    d          = round(d_W_Help$Cohens_d, 2),
    CI         = paste0("[", round(d_W_Help$CI_low, 2), ", ", round(d_W_Help$CI_high, 2), "]")
  ),
  data.frame(
    Occupation = "Welders",
    Predictor  = "Motivation",
    d          = round(d_W_Mot$Cohens_d, 2),
    CI         = paste0("[", round(d_W_Mot$CI_low, 2), ", ", round(d_W_Mot$CI_high, 2), "]")
  )
) %>%
  rename(`Cohen's d` = d)

# render
kable(
  table2,
  caption = "Table 2. Cohen’s d by occupation and predictor",
  align   = "c"
)

Table 2. Cohen’s d by occupation and predictor

Occupation	Predictor	Cohen’s d	CI
Firefighters	Helpfulness	0.15	[-0.08, 0.39]
Firefighters	Motivation	0.13	[-0.11, 0.37]
Psychiatrists	Helpfulness	0.27	[0.03, 0.51]
Psychiatrists	Motivation	0.19	[-0.05, 0.42]
Welders	Helpfulness	1.14	[0.88, 1.39]
Welders	Motivation	0.09	[-0.15, 0.33]

Once again, we confirmed that perceived heroism strongly correlated with all the measures of interest: perceived bravery, perceived selflessness, perceived exposure to danger, and perceived helpfulness (see Table 6).

PerformanceAnalytics::chart.Correlation(Set[, c("Heroism", "Danger", "Helpfulness", "Attitude")], method = "spearman")

## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties
## Warning in cor.test.default(as.numeric(x), as.numeric(y), method = method):
## Cannot compute exact p-value with ties

Additional unreported analyses

Diagnostics

mod <-lm(Heroism ~ Help_dummy * Motivation_dummy, data = Set)
plot(mod)

Outliers analyses

As registered: we’ll check for outliers by comparing robust vs non robust models.

Cook’s distance plot does show that there are some influential cases:

ols_plot_cooksd_bar(mod, type = 1)

To account for potential influential cases, we run Robust regression (huber weighted, MM-estimates, tuned to achieve 50% breakdown and 95% Gaussian efficiency).

summary(modrob<-lmrob(Set$Heroism ~ Set$Help_dummy*Set$Motivation_dummy))

## 
## Call:
## lmrob(formula = Set$Heroism ~ Set$Help_dummy * Set$Motivation_dummy)
##  \--> method = "MM"
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.38635 -0.91267  0.08733  1.17555  2.17555 
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                          5.21183    0.05526  94.310  < 2e-16 ***
## Set$Help_dummy                       0.68653    0.10750   6.386 2.86e-10 ***
## Set$Motivation_dummy                 0.21285    0.10766   1.977   0.0484 *  
## Set$Help_dummy:Set$Motivation_dummy  0.24927    0.21505   1.159   0.2467    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 1.51 
## Multiple R-squared:  0.05734,    Adjusted R-squared:  0.05385 
## Convergence in 10 IRWLS iterations
## 
## Robustness weights: 
##  51 weights are ~= 1. The remaining 765 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3790  0.8848  0.9455  0.9147  0.9850  0.9988 
## Algorithmic parameters: 
##        tuning.chi                bb        tuning.psi        refine.tol 
##         1.548e+00         5.000e-01         4.685e+00         1.000e-07 
##           rel.tol         scale.tol         solve.tol          zero.tol 
##         1.000e-07         1.000e-10         1.000e-07         1.000e-10 
##       eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
##         1.225e-04         1.819e-12         5.000e-01         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##            500             50              2              1            200 
##    maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
##            200              0           1000              0           2000 
##                   psi           subsampling                   cov 
##            "bisquare"         "nonsingular"         ".vcov.avar1" 
## compute.outlier.stats 
##                  "SM" 
## seed : int(0)

paste0("Weights applied to residuals - a value of zero would mean that the observation was discarded, a value of 1 means no re-weighting")

## [1] "Weights applied to residuals - a value of zero would mean that the observation was discarded, a value of 1 means no re-weighting"

plot(modrob$rweights)

==> It does not influence our inferences.

Discussion

This study assessed if the altruistic component of heroism was more related to selfless motivations or to being factually helping people. Indeed, it is possible some workers could be motivated to help people, despite their occupations being only relatively helpful – and conversely, some occupations might be occupied by workers who are not necessarily driven by the idea of helping others, but regardless of their motivation, their work helps a lot. We observed that, at least for underwater welders and psychiatrists, our framing of altruistic outcomes produced medium-to-large effect sizes, whereas our framing of altruistic motivations produced small and non-significant effect sizes.

Confirming correlational observations from previous studies, heroism perceptions were significantly more driven by being helpful, above and beyond, being selfless (at least in the subset of occupations tested in study 3). In contrast, being framed as selfless appeared to have a weak effect on perceived heroism (replicating study 2 small effect sizes).

Study 3

Across two experimental study, we observed overvall evidence for causal links between risk (study 1) and altruism (study 2) framing of occupations on heroism perception. Nevertheless, when decomposing the effects between each target occupation, our manipulations appeared to be effective only in subsamples of occupations. This raises the question of the universality of a risk + altruism model of collective heroism perception.

In a wide-scope correlational and descriptive study, we conducted a large assessment of a representative sample of 53 occupations among 430 participants. The aim of this study was to test the predictive power of judgments of bravery, actual danger, selflessness, and actual helping upon judgments of the whole occupation as heroic. Data, materials, codes are publicly available online (see https://osf.io/fj7pw/?view_only=2751ad848d4647348c677d09a960942e). The study was pre-registered (https://osf.io/n6mp3).

Protocol

Participants were recruited online using the representative sampling procedure provided by Prolific. We registered a representative sample of 430 UK residents. Three participants failed the attention checks, and one participant was timed-out by the Prolific platform despite completing the survey. We thus recruited four additional participants to reach our target sample size of 430 participants (225 women, 204 men, and one participant identifying as other, mean age = 47.07, SD = 15.82). After completing the informed consent and a commitment check to ensure attention to the study (see Geisen, 2022; Peer et al., 2024), participants rated 53 target occupations on perceived heroism, bravery, selflessness, exposure to danger, and helpfulness. The occupations were presented in a new random order for each evaluated feature, and the order of the evaluated features themselves was also randomized, except for attitude, which was measured at the end of the survey.

Target occupations. Our sample of occupations was determined using recent data about 580 occupations evaluated by UK residents (Newlands & Lutz, 2024), including ratings of social value for each occupation. Our goal was to include occupations that participants would recognize, spanning all major categories of the International Standard Classification of Occupations (ISCO) groups and capturing high and low social value ratings within each group. For each major ISCO group (e.g., Armed forces; Managers, etc.), we categorized each occupation in quintiles based on their social value. We then selected one to two occupations per quintile per ISCO group using our best judgment. Our selection aimed to obtain familiar occupations where possible, as well as to include occupations that we were particularly interested in as heroes (i.e., soldiers, nurses, firefighters, police officers, but also underwater welders6 and psychiatrists)7. Details on the protocol to select our final sample of occupations can be found on our OSF webpage.

Heroism. A single item (“To what extent do you believe the following occupations to be heroes?”) scaled from ‘1 - Not at all a hero’ to ‘7 - Very much a hero’.

Bravery. A single item (“To what extent do you believe the following occupations to be cowardly or brave?”) scaled from ‘1 - Cowardly’ to ‘7 - Brave’.

Selflessness. A single item (“To what extent do you believe the following occupations to be selfish or selfless?”) scaled from ‘1 - Selfish’ to ‘7 - Selfless’.

Actual helpfulness. A single item (“To what extent do you believe the following occupations to actually help other people?”) scaled from ‘1 - Not at all helpful’ to ‘7 - Very much helpful’.

Exposure to danger. A single item (“To what extent do you believe the following occupations to be exposed to physical danger?”) scaled from ‘1 - Not at all exposed to danger’ to ‘7 - Very much exposed to danger’.

Attitude. A single item (“What is your attitude (from ‘1 - Very negative’ to ‘7 - Very positive’) toward the following occupations?”) .

Attention checks. The survey included five items instructing participants to select a specific response (e.g., ‘Select 3’, ‘Tick 5’) instead of choosing an occupation. Participants who failed two or more of these attention checks were excluded from the survey.

Results

CorrSet <- read.csv("Raw_Data_50Shades.csv", comment.char="#")

#head(CorrSet)
CorrSet <- CorrSet[ -c(1,2),]
CorrSet<- subset(CorrSet, CorrSet$Gender != "")

Test <- read.csv("Prolific_Demographics_50Shades.csv")
Test <- subset(Test, Test$Status == "AWAITING REVIEW")

CorrSet <- CorrSet[-which(CorrSet$ProlID %in%setdiff(CorrSet$ProlID, Test$Participant.id)),]

Demog<- CorrSet[, c(357:361)]

occupation_labels <- c(
  "Army sergeant", 
  "Admiral", 
  "Soldier", 
  "Organised crime boss",
  "Shopping centre manager", 
  "Political party leader", 
  "Construction project manager",
  "Foreign diplomat", 
  "Street musician", 
  "Investment banker", 
  "Journalist",
  "Bridge engineer", 
  "Psychiatrist", 
  "Wedding planner", 
  "Nurse",
  "Fitness instructor", 
  "Computer help desk operator", 
  "Building surveyor",
  "Air traffic controller", 
  "Debt collector", 
  "Social media community manager",
  "Charity worker", 
  "Secretary", 
  "Postal worker", 
  "Fashion model", 
  "Bodyguard",
  "Bus conductor", 
  "Caretaker", 
  "Firefighter", 
  "Police officer", 
  "Cannabis farmer",
  "Tree feller", 
  "Fisher", 
  "Beekeeper", 
  "Farmer", 
  "Cigar maker", 
  "Goldsmith",
  "Factory worker", 
  "Underwater welder", 
  "Roofer", 
  "Machine operator",
  "Uber driver", 
  "Miner", 
  "Ambulance driver", 
  "Leaflet distributor",
  "Food delivery courier", 
  "Window cleaner", 
  "Garden labourer", 
  "Tree planter",
  "Unemployed", 
  "Retiree", 
  "Student", 
  "Homemaker", 
  "Attention check"
)

# 2. reshape
long_df <- CorrSet %>%
  dplyr::select(ResponseId,
         starts_with("Hero_"),    # Heroes
         starts_with("Brav_"),    # Bravery
         starts_with("Selfle_"),  # Selflessness
         starts_with("Dang_"),    # Danger
         starts_with("Helpf_"),   # Help
         starts_with("Attitude_") # Attitude
  ) %>%
  pivot_longer(
    cols      = -ResponseId,
    names_to  = c("Scale", "OccCode"),
    names_sep = "_",
    values_to = "Score"
  ) %>%
  # 3. drop attention-check
  dplyr::filter(OccCode != "54") %>%
  # 4. spread the Scale back out
  pivot_wider(
    names_from  = Scale,
    values_from = Score
  ) %>%
  # 5. map OccCode → Occupation label and cast numeric
  dplyr::mutate(
    Occupation = factor(
      as.integer(OccCode),
      levels = 1:54,
      labels = occupation_labels
    )
  ) %>%
  # 6. rename columns to exactly what you want
  dplyr::transmute(
    Participant   = ResponseId,
    Occupation,
    Heroes        = Hero,
    Bravery       = Brav,
    Selflessness  = Selfle,
    Danger        = Dang,
    Help          = Helpf,
    Attitude      = Attitude
  )
DF<-long_df

Demographics

table(Demog$Gender)

## 
##                 Female                   Male Other (Please specify) 
##                    225                    204                      1

paste0("430 participants took part in the study. Mean age in the sample is ", mean(as.numeric(Demog$Q19)), ", SD = ", sd(as.numeric(Demog$Q19)))

## [1] "430 participants took part in the study. Mean age in the sample is 47.0697674418605, SD = 15.8206665300891"

## Gender

Demog %>% group_by(Gender) %>% summarise(N=n()) %>%
  ggplot(aes(x=Gender,y=N,fill=Gender))+
  geom_bar(stat = 'identity',color='black')+
  scale_y_continuous(labels = scales::comma_format(accuracy = 2))+
  geom_text(aes(label=N),vjust=-0.25,fontface='bold')+
  theme_bw()+
  theme(axis.text = element_text(color='black',face='bold'),
        axis.title = element_text(color='black',face='bold'),
        legend.text = element_text(color='black',face='bold'),
        legend.title = element_text(color='black',face='bold')) +
  ggtitle("Gender distribution")

Main reported analyses - as registered

As registered, we evaluated the associations between each dimension of heroism and heroism by conducting independent mixed models (one per predictor) accounting for attitude as a covariate and using random intercept for each occupation level to model random slopes for each main predictor in each model. That way, it was possible to evaluate how each dimension predicted heroism differently across occupations by plotting the random effects in caterpillar plots. All variables were standardized prior to analysis. Across the four mixed-effects models, perceived bravery, selflessness, exposure to danger, and helpfulness each significantly predicted perceived heroism above and beyond attitudes toward occupations (all ps < .001). All models showed substantial between-occupation variability (ICCs = .52–.57), and random slopes indicated that while most occupations followed the average pattern, certain groups (e.g., soldiers and firefighters) were especially heroized for bravery, danger, or selflessness, whereas others (e.g., unemployed individuals, cigar makers) showed weaker associations. Details of these independent models are provided in our SOM. As registered, we explored the interactions between the character evaluations (i.e., bravery and selflessness) on the one hand, and the situation evaluations (i.e., exposure to danger and helpfulness) on the other hand. In both models, the interactions were significant, but of a smaller magnitude than the main independent effects of the predictors (see Table 7).

DF <- DF %>%
  dplyr::mutate(
    Heroes_num = dplyr::recode(
      Heroes,
      "Not at all a hero" = 1,
      "Very little"        = 2,
      "A little"           = 3,
      "Somewhat"           = 4,
      "Quite a bit"        = 5,
      "A lot"              = 6,
      "Very much a hero"   = 7,
      .default = NA_real_   # just in case there’s a typo or new level
    )
  )

DF <- DF %>%
  dplyr::mutate(
    Bravery_num = dplyr::recode(
      Bravery,
      "1 - Cowardly" = 1,
      "2" = 2,
      "3" = 3,
      "4" = 4,
      "5" = 5,
      "6" = 6,
      "7 - Brave"   = 7,
      .default = NA_real_   # just in case there’s a typo or new level
    )
  )

DF <- DF %>%
  dplyr::mutate(
    Selflessness_num = dplyr::recode(
      Selflessness,
      "1 - Selfish" = 1,
      "2" = 2,
      "3" = 3,
      "4" = 4,
      "5" = 5,
      "6" = 6,
      "7 - Selfless"   = 7,
      .default = NA_real_   # just in case there’s a typo or new level
    )
  )

DF <- DF %>%
  dplyr::mutate(
    Danger_num = dplyr::recode(
      Danger,
      "Not at all exposed to danger" = 1,
      "Very little" = 2,
      "A little" = 3,
      "Somewhat" = 4,
      "Quite a bit" = 5,
      "A lot" = 6,
      "Very much exposed to danger"   = 7,
      .default = NA_real_   # just in case there’s a typo or new level
    )
  )

DF <- DF %>%
  dplyr::mutate(
    Help_num = dplyr::recode(
      Help,
      "Not at all helpful" = 1,
      "Very little" = 2,
      "A little" = 3,
      "Somewhat" = 4,
      "Quite a bit" = 5,
      "A lot" = 6,
      "Very much helpful"   = 7,
      .default = NA_real_   # just in case there’s a typo or new level
    )
  )


DF <- DF %>%
  dplyr::mutate(
    Attitude_num = dplyr::recode(
      Attitude,
      "Very negative" = 1,
      "Quite negative" = 2,
      "Slightly negative" = 3,
      "Neutral" = 4,
      "Slightly positive" = 5,
      "Quite positive" = 6,
      "Very positive"   = 7,
      .default = NA_real_   # just in case there’s a typo or new level
    )
  )


DF <- DF %>%
  mutate(
    zHeroes = as.numeric(scale(Heroes_num)), 
    zBravery     = as.numeric(scale(Bravery_num)),
    zSelflessness= as.numeric(scale(Selflessness_num)),
    zDanger      = as.numeric(scale(Danger_num)),
    zHelp        = as.numeric(scale(Help_num)),
    zAttitude    = as.numeric(scale(Attitude_num)),
    diff_bs      = zBravery - zSelflessness,
    diff_dh      = zDanger   - zHelp
  )


df2 <- DF %>%
  mutate(
    zHeroes = as.numeric(scale(Heroes_num)), 
    zBravery     = as.numeric(scale(Bravery_num)),
    zSelflessness= as.numeric(scale(Selflessness_num)),
    zDanger      = as.numeric(scale(Danger_num)),
    zHelp        = as.numeric(scale(Help_num)),
    diff_bs      = zBravery - zSelflessness,
    diff_dh      = zDanger   - zHelp
  )


OccupDF <- df2 %>%
  group_by(Occupation) %>%
  summarise(
    Heroes_Mean    = mean(Heroes_num),
    Brav_Mean      = mean(Bravery_num),
    Selfl_Mean     = mean(Selflessness_num),
    Danger_Mean    = mean(Danger_num),
    Help_Mean      = mean(Help_num),
    Attitude_Mean  = mean(Attitude_num),
    mean_zBrav     = mean(zBravery),
    mean_zSelfl    = mean(zSelflessness),
    mean_zDang     = mean(zDanger),
    mean_zHelp     = mean(zHelp),
    mean_zHeroism  = mean(zHeroes),
    mean_zAttitude = mean(zAttitude),
    diff_bs        = mean(diff_bs),
    diff_dh        = mean(diff_dh)
  ) 



OccupDF2 <- OccupDF %>%
  mutate(
    # original ones
    ambiv_brav_self   = mean_zBrav    * mean_zSelfl   < 0,  # Bravery vs Selflessness
    ambiv_dang_help   = mean_zDang    * mean_zHelp    < 0,  # Danger vs Helpfulness
    
    # vs Heroism
    ambiv_brav_hero   = mean_zBrav    * mean_zHeroism < 0,  # Bravery vs Heroism
    ambiv_self_hero   = mean_zSelfl   * mean_zHeroism < 0,  # Selflessness vs Heroism
    ambiv_dang_hero   = mean_zDang    * mean_zHeroism < 0,  # Danger vs Heroism
    ambiv_help_hero   = mean_zHelp    * mean_zHeroism < 0,  # Helpfulness vs Heroism
    
    # the remaining cross-combos
    ambiv_brav_dang   = mean_zBrav    * mean_zDang    < 0,  # Bravery vs Danger
    ambiv_brav_help   = mean_zBrav    * mean_zHelp    < 0,  # Bravery vs Helpfulness
    ambiv_self_dang   = mean_zSelfl   * mean_zDang    < 0,  # Selflessness vs Danger
    ambiv_self_help   = mean_zSelfl   * mean_zHelp    < 0   # Selflessness vs Helpfulness
  )


model <- lmerTest::lmer(scale(Heroes_num) ~ zBravery + zAttitude + 
                (zBravery | Occupation) +
                (1 | Participant), data = DF)
summary(model)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(Heroes_num) ~ zBravery + zAttitude + (zBravery | Occupation) +  
##     (1 | Participant)
##    Data: DF
## 
## REML criterion at convergence: 35618.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6069 -0.6032 -0.0500  0.5628  6.7687 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev. Corr
##  Participant (Intercept) 0.19898  0.4461       
##  Occupation  (Intercept) 0.08172  0.2859       
##              zBravery    0.01791  0.1338   0.83
##  Residual                0.25614  0.5061       
## Number of obs: 22790, groups:  Participant, 430; Occupation, 53
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept) -7.026e-02  4.497e-02  8.480e+01  -1.562    0.122    
## zBravery     2.600e-01  1.918e-02  5.446e+01  13.560   <2e-16 ***
## zAttitude    2.853e-01  5.362e-03  2.065e+04  53.208   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) zBrvry
## zBravery   0.695       
## zAttitude  0.003 -0.106

model_performance(model)

report(model)

## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zBravery and zAttitude (formula: scale(Heroes_num) ~
## zBravery + zAttitude). The model included zBravery as random effects (formula:
## list(~zBravery | Occupation, ~1 | Participant)). The model's total explanatory
## power is substantial (conditional R2 = 0.67) and the part related to the fixed
## effects alone (marginal R2) is of 0.29. The model's intercept, corresponding to
## zBravery = 0 and zAttitude = 0, is at -0.07 (95% CI [-0.16, 0.02], t(22782) =
## -1.56, p = 0.118). Within this model:
## 
##   - The effect of zBravery is statistically significant and positive (beta =
## 0.26, 95% CI [0.22, 0.30], t(22782) = 13.56, p < .001; Std. beta = 0.26, 95% CI
## [0.22, 0.30])
##   - The effect of zAttitude is statistically significant and positive (beta =
## 0.29, 95% CI [0.27, 0.30], t(22782) = 53.21, p < .001; Std. beta = 0.29, 95% CI
## [0.27, 0.30])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

icc(model) 

Bravery Model

#colnames(DF)
model <- lmerTest::lmer(scale(Heroes_num) ~ zBravery + zAttitude + 
                (zBravery | Occupation) +
                (1 | Participant), data = DF)
summary(model)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(Heroes_num) ~ zBravery + zAttitude + (zBravery | Occupation) +  
##     (1 | Participant)
##    Data: DF
## 
## REML criterion at convergence: 35618.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6069 -0.6032 -0.0500  0.5628  6.7687 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev. Corr
##  Participant (Intercept) 0.19898  0.4461       
##  Occupation  (Intercept) 0.08172  0.2859       
##              zBravery    0.01791  0.1338   0.83
##  Residual                0.25614  0.5061       
## Number of obs: 22790, groups:  Participant, 430; Occupation, 53
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept) -7.026e-02  4.497e-02  8.480e+01  -1.562    0.122    
## zBravery     2.600e-01  1.918e-02  5.446e+01  13.560   <2e-16 ***
## zAttitude    2.853e-01  5.362e-03  2.065e+04  53.208   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) zBrvry
## zBravery   0.695       
## zAttitude  0.003 -0.106

model_performance(model)

report(model)

## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zBravery and zAttitude (formula: scale(Heroes_num) ~
## zBravery + zAttitude). The model included zBravery as random effects (formula:
## list(~zBravery | Occupation, ~1 | Participant)). The model's total explanatory
## power is substantial (conditional R2 = 0.67) and the part related to the fixed
## effects alone (marginal R2) is of 0.29. The model's intercept, corresponding to
## zBravery = 0 and zAttitude = 0, is at -0.07 (95% CI [-0.16, 0.02], t(22782) =
## -1.56, p = 0.118). Within this model:
## 
##   - The effect of zBravery is statistically significant and positive (beta =
## 0.26, 95% CI [0.22, 0.30], t(22782) = 13.56, p < .001; Std. beta = 0.26, 95% CI
## [0.22, 0.30])
##   - The effect of zAttitude is statistically significant and positive (beta =
## 0.29, 95% CI [0.27, 0.30], t(22782) = 53.21, p < .001; Std. beta = 0.29, 95% CI
## [0.27, 0.30])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

icc(model) 

slopes_df <- ranef(model)$Occupation %>%
  tibble::rownames_to_column("Occupation") %>%
  rename(slope = zBravery)

p1 <- ggplot(slopes_df, aes(x = slope, 
                           y = reorder(Occupation, slope),
                           text = paste0("Occupation: ", Occupation,
                                         "<br>Slope: ", round(slope, 3)))) +
  geom_point() +
  labs(
    title = "Bravery → Heroism Slopes by Occupation",
    x     = "Random Slope for Bravery",
    y     = "Occupation"
  )
# Convert to an interactive Plotly plot
ggplotly(p1, tooltip = "text")

Selflessness Model

model <- lmerTest::lmer(scale(Heroes_num) ~ zSelflessness + zAttitude + 
                (zSelflessness | Occupation) +
                (1 | Participant), data = DF)
summary(model)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(Heroes_num) ~ zSelflessness + zAttitude + (zSelflessness |  
##     Occupation) + (1 | Participant)
##    Data: DF
## 
## REML criterion at convergence: 36572.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5083 -0.6080 -0.0450  0.5617  6.7037 
## 
## Random effects:
##  Groups      Name          Variance Std.Dev. Corr
##  Participant (Intercept)   0.20831  0.4564       
##  Occupation  (Intercept)   0.12509  0.3537       
##              zSelflessness 0.01521  0.1233   0.84
##  Residual                  0.26694  0.5167       
## Number of obs: 22790, groups:  Participant, 430; Occupation, 53
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   -6.205e-02  5.350e-02  7.431e+01   -1.16     0.25    
## zSelflessness  1.932e-01  1.787e-02  5.641e+01   10.81 2.31e-15 ***
## zAttitude      2.982e-01  5.673e-03  2.192e+04   52.55  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) zSlfls
## zSelflssnss  0.721       
## zAttitude    0.000 -0.138

model_performance(model)

report(model)

## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zSelflessness and zAttitude (formula:
## scale(Heroes_num) ~ zSelflessness + zAttitude). The model included
## zSelflessness as random effects (formula: list(~zSelflessness | Occupation, ~1
## | Participant)). The model's total explanatory power is substantial
## (conditional R2 = 0.67) and the part related to the fixed effects alone
## (marginal R2) is of 0.25. The model's intercept, corresponding to zSelflessness
## = 0 and zAttitude = 0, is at -0.06 (95% CI [-0.17, 0.04], t(22782) = -1.16, p =
## 0.246). Within this model:
## 
##   - The effect of zSelflessness is statistically significant and positive (beta =
## 0.19, 95% CI [0.16, 0.23], t(22782) = 10.81, p < .001; Std. beta = 0.19, 95% CI
## [0.16, 0.23])
##   - The effect of zAttitude is statistically significant and positive (beta =
## 0.30, 95% CI [0.29, 0.31], t(22782) = 52.55, p < .001; Std. beta = 0.30, 95% CI
## [0.29, 0.31])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

icc(model) 

slopes_df <- ranef(model)$Occupation %>%
  tibble::rownames_to_column("Occupation") %>%
  rename(slope = zSelflessness)

p3 <- ggplot(slopes_df, aes(x = slope, 
                           y = reorder(Occupation, slope),
                           text = paste0("Occupation: ", Occupation,
                                         "<br>Slope: ", round(slope, 3)))) +
  geom_point() +
  labs(
    title = "Selflessness → Heroism Slopes by Occupation",
    x     = "Random Slope for Selflessness",
    y     = "Occupation"
  )
# Convert to an interactive Plotly plot
ggplotly(p3, tooltip = "text")

Risk exposure Model

model <- lmerTest::lmer(scale(Heroes_num) ~ zDanger + zAttitude + 
                (zDanger | Occupation) +
                (1 | Participant), data = DF)
summary(model)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(Heroes_num) ~ zDanger + zAttitude + (zDanger | Occupation) +  
##     (1 | Participant)
##    Data: DF
## 
## REML criterion at convergence: 36582.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1912 -0.5952 -0.0401  0.5664  6.7052 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev. Corr
##  Participant (Intercept) 0.18144  0.4260       
##  Occupation  (Intercept) 0.09706  0.3116       
##              zDanger     0.01458  0.1207   0.48
##  Residual                0.26762  0.5173       
## Number of obs: 22790, groups:  Participant, 430; Occupation, 53
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept) -4.163e-02  4.779e-02  7.688e+01  -0.871    0.386    
## zDanger      2.299e-01  1.782e-02  5.497e+01  12.899   <2e-16 ***
## zAttitude    3.364e-01  5.189e-03  2.256e+04  64.830   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) zDangr
## zDanger    0.406       
## zAttitude  0.003 -0.071

model_performance(model)

report(model)

## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zDanger and zAttitude (formula: scale(Heroes_num) ~
## zDanger + zAttitude). The model included zDanger as random effects (formula:
## list(~zDanger | Occupation, ~1 | Participant)). The model's total explanatory
## power is substantial (conditional R2 = 0.65) and the part related to the fixed
## effects alone (marginal R2) is of 0.28. The model's intercept, corresponding to
## zDanger = 0 and zAttitude = 0, is at -0.04 (95% CI [-0.14, 0.05], t(22782) =
## -0.87, p = 0.384). Within this model:
## 
##   - The effect of zDanger is statistically significant and positive (beta = 0.23,
## 95% CI [0.19, 0.26], t(22782) = 12.90, p < .001; Std. beta = 0.23, 95% CI
## [0.19, 0.26])
##   - The effect of zAttitude is statistically significant and positive (beta =
## 0.34, 95% CI [0.33, 0.35], t(22782) = 64.83, p < .001; Std. beta = 0.34, 95% CI
## [0.33, 0.35])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

icc(model) 

slopes_df <- ranef(model)$Occupation %>%
  tibble::rownames_to_column("Occupation") %>%
  rename(slope = zDanger)

p2 <- ggplot(slopes_df, aes(x = slope, 
                           y = reorder(Occupation, slope),
                           text = paste0("Occupation: ", Occupation,
                                         "<br>Slope: ", round(slope, 3)))) +
  geom_point() +
  labs(
    title = "Danger → Heroism Slopes by Occupation",
    x     = "Random Slope for Danger",
    y     = "Occupation"
  )
# Convert to an interactive Plotly plot
ggplotly(p2, tooltip = "text")

Helpfulness exposure Model

model <- lmerTest::lmer(scale(Heroes_num) ~ zHelp + zAttitude + 
                (zHelp | Occupation) +
                (1 | Participant), data = DF)
summary(model)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(Heroes_num) ~ zHelp + zAttitude + (zHelp | Occupation) +  
##     (1 | Participant)
##    Data: DF
## 
## REML criterion at convergence: 35735.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4467 -0.5950 -0.0365  0.5644  6.7451 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev. Corr
##  Participant (Intercept) 0.18639  0.4317       
##  Occupation  (Intercept) 0.10077  0.3174       
##              zHelp       0.02263  0.1504   0.90
##  Residual                0.25781  0.5078       
## Number of obs: 22790, groups:  Participant, 430; Occupation, 53
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept) -7.269e-02  4.853e-02  7.686e+01  -1.498    0.138    
## zHelp        2.654e-01  2.155e-02  5.440e+01  12.315   <2e-16 ***
## zAttitude    2.607e-01  5.642e-03  1.953e+04  46.198   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) zHelp 
## zHelp      0.774       
## zAttitude  0.001 -0.133

model_performance(model)

report(model)

## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zHelp and zAttitude (formula: scale(Heroes_num) ~
## zHelp + zAttitude). The model included zHelp as random effects (formula:
## list(~zHelp | Occupation, ~1 | Participant)). The model's total explanatory
## power is substantial (conditional R2 = 0.68) and the part related to the fixed
## effects alone (marginal R2) is of 0.29. The model's intercept, corresponding to
## zHelp = 0 and zAttitude = 0, is at -0.07 (95% CI [-0.17, 0.02], t(22782) =
## -1.50, p = 0.134). Within this model:
## 
##   - The effect of zHelp is statistically significant and positive (beta = 0.27,
## 95% CI [0.22, 0.31], t(22782) = 12.32, p < .001; Std. beta = 0.27, 95% CI
## [0.22, 0.31])
##   - The effect of zAttitude is statistically significant and positive (beta =
## 0.26, 95% CI [0.25, 0.27], t(22782) = 46.20, p < .001; Std. beta = 0.26, 95% CI
## [0.25, 0.27])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

icc(model) 

slopes_df <- ranef(model)$Occupation %>%
  tibble::rownames_to_column("Occupation") %>%
  rename(slope = zHelp)

p4 <- ggplot(slopes_df, aes(x = slope, 
                           y = reorder(Occupation, slope),
                           text = paste0("Occupation: ", Occupation,
                                         "<br>Slope: ", round(slope, 3)))) +
  geom_point() +
  labs(
    title = "Helpfulness → Heroism Slopes by Occupation",
    x     = "Random Slope for Helpfulness",
    y     = "Occupation"
  )
# Convert to an interactive Plotly plot
ggplotly(p4, tooltip = "text")

We deviated from the registration and computed a general mixed model estimating the unique contribution of each predictor (danger, help, bravery, selflessness) to perceived heroism, while controlling for random variation across both occupations and participants and using attitude as a covariate to account for possible halo effects. Consistent with the independent models, each predictor were significant (see Table 8). Interestingly, and consistent with findings from study 2, being perceived as actually helpful appears to predict heroism to a larger extent than perceiving workers as being selflessly motivated, as indicated by their non-overlapping 95% confidence intervals. Indeed, when conducting an exploratory, non-registered, linear hypothesis testing, a simultaneous Wald test indicated that the slope for Helpfulness was significantly greater than the slope for Selflessness, delta beta = 0.081, z =9.15, p<.001. Excluding helpfulness from the full model resulted in a drop of 6.7% of the explained variance, whereas excluding selflessness from the full model resulted in a more modest drop of 2.3% as indicated by semi-partial R2.

m <- lmerTest::lmer(
  scale(Heroes_num) ~ zDanger + zHelp + zBravery + zSelflessness + zAttitude +
          (1 | Occupation) +
                (1 | Participant),
  data   = DF,
  REML   = TRUE
)

res <- partR2(
  m, 
  data = DF,
  partvars = c("zBravery", "zDanger", "zHelp", "zSelflessness", "zAttitude"),
  R2_type = "marginal", 
  nboot = 3
)

r2_table <- res$R2 %>% dplyr::filter(term %in% c("zBravery", "zDanger", "zHelp", "zSelflessness", "zAttitude"))

# Create the gt table
r2_table %>%
  gt() %>%
  tab_header(
    title = md("**Part R² for Each Predictor**"),
    subtitle = "Marginal R² with 95% CI"
  ) %>%
  cols_label(
    term = "Predictor",
    estimate = "Part R²",
    CI_lower = "95% CI (Lower)",
    CI_upper = "95% CI (Upper)",
    ndf = "Num. DF"
  ) %>%
  fmt_number(
    columns = c(estimate, CI_lower, CI_upper),
    decimals = 3
  ) %>%
  tab_options(
    table.font.size = px(14),
    column_labels.font.weight = "bold"
  )

summary(glht(m, linfct = "zHelp - zSelflessness = 0"))


summary(m)
model_performance(m)
report(m)
icc(m)



# 1. Get a table of fixed effects, standardized (by refitting the model with all variables
#    on their original scales but extracting standardized slopes), plus CIs, t, and df
tab <- model_parameters(
  m,
  effects     = "fixed",        # only fixed effects
  standardize = "refit",        # compute beta on refit model → "standardized β"
  ci          = 0.95            # 95% confidence intervals
)
names(tab)
# 2. Tidy up: select and rename the columns you asked for
tab2 <- tab %>%
  select(Parameter,
         `Std. beta` = Coefficient,
         `t value` = t,
         df = df_error,
         `Low CI` = CI_low,
         `High CI` = CI_high) 

# 3. Print a neat markdown table
knitr::kable(
  tab2,
  col.names = c("Predictor", "Std. beta", "t value", "df", "Low CI", "High CI"),
  digits = c(NA, 2, 2, 0, 2, 2),
  caption = "Fixed‐effect estimates from mixed model"
)

Across all occupations, each dimension of heroism correlated strongly with heroism when partialing out attitude, with correlations ranging from r = .32 (Bravery) to r = .79 (helpfulness), see Figure 1A-1D. Further visualisations can be found in the SOM.

# 1) Residualize each z-scored variable on mean_zAttitude
OccupDF3 <- OccupDF2 %>%
  mutate(
    resid_zBrav    = resid(lm(mean_zBrav    ~ mean_zAttitude, data = .)),
    resid_zSelfl   = resid(lm(mean_zSelfl   ~ mean_zAttitude, data = .)),
    resid_zDang    = resid(lm(mean_zDang    ~ mean_zAttitude, data = .)),
    resid_zHelp    = resid(lm(mean_zHelp    ~ mean_zAttitude, data = .)),
    resid_zHeroism = resid(lm(mean_zHeroism ~ mean_zAttitude, data = .))
  ) %>%
  # 2) Compute ambivalence flags on the residuals
  mutate(
    ambiv_brav_self   = resid_zBrav    * resid_zSelfl   < 0,
    ambiv_dang_help   = resid_zDang    * resid_zHelp    < 0,
    ambiv_brav_hero   = resid_zBrav    * resid_zHeroism < 0,
    ambiv_self_hero   = resid_zSelfl   * resid_zHeroism < 0,
    ambiv_dang_hero   = resid_zDang    * resid_zHeroism < 0,
    ambiv_help_hero   = resid_zHelp    * resid_zHeroism < 0,
    ambiv_brav_dang   = resid_zBrav    * resid_zDang    < 0,
    ambiv_brav_help   = resid_zBrav    * resid_zHelp    < 0,
    ambiv_self_dang   = resid_zSelfl   * resid_zDang    < 0,
    ambiv_self_help   = resid_zSelfl   * resid_zHelp    < 0
  )

# 3) Re-use the helper function to make interactive plots
make_resid_plot <- function(df, xvar, yvar, ambiv_var, title, xlab, ylab) {
  p <- ggplot(df, aes_string(
         x = xvar,
         y = yvar,
         text = "paste0('Occupation: ', Occupation,
                       '<br>Heroism (mean): ', round(Heroes_Mean,2))"
       )) +
       geom_point(aes_string(color = ambiv_var), size = 3) +
       scale_color_manual(values = c("FALSE" = "grey50", "TRUE" = "firebrick")) +
       geom_hline(yintercept = 0, linetype = "dashed") +
       geom_vline(xintercept = 0, linetype = "dashed") +
       labs(title = title,
            x     = xlab,
            y     = ylab,
            color = "Ambivalent?") +
       theme_minimal()
  ggplotly(p, tooltip = c("text","x","y"))
}

library(ggrepel)   # install.packages("ggrepel") if you haven’t already
make_resid_plot2 <- function(df, xvar, yvar, ambiv_var, title, xlab, ylab) {
  p <- ggplot(df, aes_string(
       x = xvar,
       y = yvar
     )) +
  geom_point(aes_string(color = ambiv_var), size = 3) +
  # add this layer to put the occupation name next to each point
  geom_text_repel(aes_string(label = "Occupation"),
                  size = 3,            # adjust text size
                  max.overlaps = 30,   # let it try to place up to 30 labels
                  box.padding = 0.3,
                  point.padding = 0.5) +
  scale_color_manual(values = c("FALSE" = "grey50", "TRUE" = "firebrick")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(title = title,
       x     = xlab,
       y     = ylab,
       color = "Ambivalent?") +
  theme_minimal()
}



# 4) Create your 10 residual-based interactive plots
p1_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zBrav", yvar = "resid_zHeroism",
                 ambiv_var = "ambiv_brav_hero",
                 title     = "Bravery vs Heroism (residuals)",
                 xlab      = "Residual z-Bravery",
                 ylab      = "Residual z-Heroism")

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

p2_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zSelfl", yvar = "resid_zHeroism",
                 ambiv_var = "ambiv_self_hero",
                 title     = "Selflessness vs Heroism (residuals)",
                 xlab      = "Residual z-Selflessness",
                 ylab      = "Residual z-Heroism")

p3_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zDang", yvar = "resid_zHeroism",
                 ambiv_var = "ambiv_dang_hero",
                 title     = "Danger vs Heroism (residuals)",
                 xlab      = "Residual z-Danger",
                 ylab      = "Residual z-Heroism")

p4_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zHelp", yvar = "resid_zHeroism",
                 ambiv_var = "ambiv_help_hero",
                 title     = "Helpfulness vs Heroism (residuals)",
                 xlab      = "Residual z-Helpfulness",
                 ylab      = "Residual z-Heroism")

p5_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zSelfl", yvar = "resid_zBrav",
                 ambiv_var = "ambiv_brav_self",
                 title     = "Selflessness vs Bravery (residuals)",
                 xlab      = "Residual z-Selflessness",
                 ylab      = "Residual z-Bravery")

p6_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zDang", yvar = "resid_zBrav",
                 ambiv_var = "ambiv_brav_dang",
                 title     = "Danger vs Bravery (residuals)",
                 xlab      = "Residual z-Danger",
                 ylab      = "Residual z-Bravery")

p7_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zHelp", yvar = "resid_zBrav",
                 ambiv_var = "ambiv_brav_help",
                 title     = "Helpfulness vs Bravery (residuals)",
                 xlab      = "Residual z-Helpfulness",
                 ylab      = "Residual z-Bravery")

p8_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zDang", yvar = "resid_zSelfl",
                 ambiv_var = "ambiv_self_dang",
                 title     = "Danger vs Selflessness (residuals)",
                 xlab      = "Residual z-Danger",
                 ylab      = "Residual z-Selflessness")

p9_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zHelp", yvar = "resid_zSelfl",
                 ambiv_var = "ambiv_self_help",
                 title     = "Helpfulness vs Selflessness (residuals)",
                 xlab      = "Residual z-Helpfulness",
                 ylab      = "Residual z-Selflessness")

p10_resid <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zHelp", yvar = "resid_zDang",
                 ambiv_var = "ambiv_dang_help",
                 title     = "Helpfulness vs Danger (residuals)",
                 xlab      = "Residual z-Helpfulness",
                 ylab      = "Residual z-Danger")

p1_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zBrav", yvar = "resid_zHeroism",
                 ambiv_var = "ambiv_brav_hero",
                 title     = "Bravery vs Heroism (residuals)",
                 xlab      = "Bravery (standardized and partialled out for attitude)",
                 ylab      = "Heroism (standardized and partialled out for attitude)")

p2_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zSelfl", yvar = "resid_zHeroism",
                 ambiv_var = "ambiv_self_hero",
                 title     = "Selflessness vs Heroism (residuals)",
xlab      = "Selflessness (standardized and partialled out for attitude)",
ylab      = "Heroism (standardized and partialled out for attitude)")

p3_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zDang", yvar = "resid_zHeroism",
                 ambiv_var = "ambiv_dang_hero",
                 title     = "Danger vs Heroism (residuals)",
xlab      = "Exposure to danger (standardized and partialled out for attitude)",
ylab      = "Heroism (standardized and partialled out for attitude)")
p4_resid  <- make_resid_plot2(OccupDF3,
                 xvar      = "resid_zHelp", yvar = "resid_zHeroism",
                 ambiv_var = "ambiv_help_hero",
                 title     = "Helpfulness vs Heroism (residuals)",
xlab      = "Helpfulness (standardized and partialled out for attitude)",
ylab      = "Heroism (standardized and partialled out for attitude)")

Bravery:Heroism

paste0("Correlation = ", round(cor(OccupDF3$resid_zBrav, OccupDF3$resid_zHeroism, use = "complete.obs"), 2))

## [1] "Correlation = 0.81"

p1_resid

Selflessness:Heroism

paste0("Correlation = ", round(cor(OccupDF3$resid_zSelfl, OccupDF3$resid_zHeroism, use = "complete.obs"), 2))

## [1] "Correlation = 0.52"

p2_resid

Danger:Heroism

paste0("Correlation = ", round(cor(OccupDF3$resid_zDang, OccupDF3$resid_zHeroism, use = "complete.obs"), 2))

## [1] "Correlation = 0.72"

p3_resid

Helpfulness:Heroism

paste0("Correlation = ", round(cor(OccupDF3$resid_zHelp, OccupDF3$resid_zHeroism, use = "complete.obs"), 2))

## [1] "Correlation = 0.47"

p4_resid

Selflessness:Bravery

paste0("Correlation = ", round(cor(OccupDF3$resid_zBrav, OccupDF3$resid_zSelfl, use = "complete.obs"), 2))

## [1] "Correlation = 0.39"

p5_resid

Danger:Bravery

paste0("Correlation = ", round(cor(OccupDF3$resid_zBrav, OccupDF3$resid_zDang, use = "complete.obs"), 2))

## [1] "Correlation = 0.94"

p6_resid

Helpfulness:Bravery

paste0("Correlation = ", round(cor(OccupDF3$resid_zBrav, OccupDF3$resid_zHelp, use = "complete.obs"), 2))

## [1] "Correlation = 0.33"

p7_resid

Danger:Selflessness

paste0("Correlation = ", round(cor(OccupDF3$resid_zDang, OccupDF3$resid_zSelfl, use = "complete.obs"), 2))

## [1] "Correlation = 0.3"

p8_resid

Helpfulness:Selflessness

paste0("Correlation = ", round(cor(OccupDF3$resid_zHelp, OccupDF3$resid_zSelfl, use = "complete.obs"), 2))

## [1] "Correlation = 0.33"

p9_resid

Helpfulness:Danger

paste0("Correlation = ", round(cor(OccupDF3$resid_zHelp, OccupDF3$resid_zDang, use = "complete.obs"), 2))

## [1] "Correlation = 0.35"

p10_resid

Additional analyses

Using raw scores

… and not accounting for attitudes

# 2) Define all 10 ggplot objects and convert to Plotly
# Helper to avoid repetition
make_plot <- function(df, xvar, yvar, ambiv_var, title, xlab, ylab, highlight = "firebrick") {
  p <- ggplot(df, aes_string(x = xvar, y = yvar,
                             text = "paste0('Occupation: ', Occupation,
                                           '<br>Heroism (mean): ', round(Heroes_Mean,2))")) +
       geom_point(aes_string(color = ambiv_var), size = 3) +
       scale_color_manual(values = c("FALSE" = "grey50", "TRUE" = highlight)) +
       geom_hline(yintercept = 0, linetype = "dashed") +
       geom_vline(xintercept = 0, linetype = "dashed") +
       labs(title = title,
            x     = xlab,
            y     = ylab,
            color = "Ambivalent?") +
       theme_minimal()
  ggplotly(p, tooltip = c("text","x","y"))
}

# 10 interactive plots
p1  <- make_plot(OccupDF2,
                 xvar     = "mean_zBrav", yvar     = "mean_zHeroism",
                 ambiv_var= "ambiv_brav_hero",
                 title    = "Bravery vs Heroism (z-scores)",
                 xlab     = "Mean z-Bravery", ylab = "Mean z-Heroism")

p2  <- make_plot(OccupDF2,
                 xvar     = "mean_zSelfl", yvar     = "mean_zHeroism",
                 ambiv_var= "ambiv_self_hero",
                 title    = "Selflessness vs Heroism (z-scores)",
                 xlab     = "Mean z-Selflessness", ylab = "Mean z-Heroism")

p3  <- make_plot(OccupDF2,
                 xvar     = "mean_zDang", yvar     = "mean_zHeroism",
                 ambiv_var= "ambiv_dang_hero",
                 title    = "Danger vs Heroism (z-scores)",
                 xlab     = "Mean z-Danger", ylab = "Mean z-Heroism")

p4  <- make_plot(OccupDF2,
                 xvar     = "mean_zHelp", yvar     = "mean_zHeroism",
                 ambiv_var= "ambiv_help_hero",
                 title    = "Helpfulness vs Heroism (z-scores)",
                 xlab     = "Mean z-Helpfulness", ylab = "Mean z-Heroism")

p5  <- make_plot(OccupDF2,
                 xvar     = "mean_zSelfl", yvar     = "mean_zBrav",
                 ambiv_var= "ambiv_brav_self",
                 title    = "Selflessness vs Bravery (z-scores)",
                 xlab     = "Mean z-Selflessness", ylab = "Mean z-Bravery")

p6  <- make_plot(OccupDF2,
                 xvar     = "mean_zDang", yvar     = "mean_zBrav",
                 ambiv_var= "ambiv_brav_dang",
                 title    = "Danger vs Bravery (z-scores)",
                 xlab     = "Mean z-Danger", ylab = "Mean z-Bravery")

p7  <- make_plot(OccupDF2,
                 xvar     = "mean_zHelp", yvar     = "mean_zBrav",
                 ambiv_var= "ambiv_brav_help",
                 title    = "Helpfulness vs Bravery (z-scores)",
                 xlab     = "Mean z-Helpfulness", ylab = "Mean z-Bravery")

p8  <- make_plot(OccupDF2,
                 xvar     = "mean_zDang", yvar     = "mean_zSelfl",
                 ambiv_var= "ambiv_self_dang",
                 title    = "Danger vs Selflessness (z-scores)",
                 xlab     = "Mean z-Danger", ylab = "Mean z-Selflessness")

p9  <- make_plot(OccupDF2,
                 xvar     = "mean_zHelp", yvar     = "mean_zSelfl",
                 ambiv_var= "ambiv_self_help",
                 title    = "Helpfulness vs Selflessness (z-scores)",
                 xlab     = "Mean z-Helpfulness", ylab = "Mean z-Selflessness")

p10 <- make_plot(OccupDF2,
                 xvar     = "mean_zHelp", yvar     = "mean_zDang",
                 ambiv_var= "ambiv_dang_help",
                 title    = "Helpfulness vs Danger (z-scores)",
                 xlab     = "Mean z-Helpfulness", ylab = "Mean z-Danger")

Brave:Heroism

paste0("Correlation = ", round(cor(OccupDF2$mean_zBrav, OccupDF2$mean_zHeroism, use = "complete.obs"), 2))

## [1] "Correlation = 0.92"

p1

p1

Selflessness:Heroism

paste0("Correlation = ", round(cor(OccupDF2$mean_zSelfl, OccupDF2$mean_zHeroism, use = "complete.obs"), 2))

## [1] "Correlation = 0.87"

p2

Danger:Heroism

paste0("Correlation = ", round(cor(OccupDF2$mean_zDang, OccupDF2$mean_zHeroism, use = "complete.obs"), 2))

## [1] "Correlation = 0.69"

p3

Helpfulness:Heroism

paste0("Correlation = ", round(cor(OccupDF2$mean_zHelp, OccupDF2$mean_zHeroism, use = "complete.obs"), 2))

## [1] "Correlation = 0.86"

p4

Selflessness:Bravery

paste0("Correlation = ", round(cor(OccupDF2$mean_zBrav, OccupDF2$mean_zSelfl, use = "complete.obs"), 2))

## [1] "Correlation = 0.79"

p5

Danger:Bravery

paste0("Correlation = ", round(cor(OccupDF2$mean_zDang, OccupDF2$mean_zBrav, use = "complete.obs"), 2))

## [1] "Correlation = 0.86"

p6

Helpfulness:Bravery

paste0("Correlation = ", round(cor(OccupDF2$mean_zBrav, OccupDF2$mean_zHelp, use = "complete.obs"), 2))

## [1] "Correlation = 0.77"

p7

Danger:Selflessness

paste0("Correlation = ", round(cor(OccupDF2$mean_zDang, OccupDF2$mean_zSelfl, use = "complete.obs"), 2))

## [1] "Correlation = 0.44"

p8

Helpfulness:Selflessness

paste0("Correlation = ", round(cor(OccupDF2$mean_zHelp, OccupDF2$mean_zSelfl, use = "complete.obs"), 2))

## [1] "Correlation = 0.91"

p9

Helpfulness:Danger

paste0("Correlation = ", round(cor(OccupDF2$mean_zHelp, OccupDF2$mean_zDang, use = "complete.obs"), 2))

## [1] "Correlation = 0.48"

p10

3D plots

we use MCD to spot cases that deviate from the bulk of the data - cases that break the mold of heroism:

X3 <- select(OccupDF2, Heroes_Mean, Brav_Mean, Selfl_Mean)
mcd3 <- covMcd(X3, alpha = 0.10)

## Warning in covMcd(X3, alpha = 0.1): subsample size h < n/2 may be too small

rdist3 <- mahalanobis(X3, center = mcd3$center, cov = mcd3$cov)

# add to df
OccupDF2$rdist3 <- rdist3


plot_ly(
  data = OccupDF2,
  x    = ~Heroes_Mean,
  y    = ~Brav_Mean,
  z    = ~Selfl_Mean,
  color = ~rdist3,                # continuous colour scale
  colors = c("steelblue","gold","firebrick"),
  text = ~Occupation,
  hovertemplate = paste(
    "Job: %{text}<br>",
    "Heroes: %{x:.2f}<br>",
    "Bravery: %{y:.2f}<br>",
    "Selflessness: %{z:.2f}<br>",
    "MCD‐dist: %{marker.color:.1f}<extra></extra>"
  ),
  type = "scatter3d",
  mode = "markers"
) %>%
  layout(
    scene = list(
      xaxis = list(title = "Heroes_Mean"),
      yaxis = list(title = "Brav_Mean"),
      zaxis = list(title = "Selfl_Mean")
    )
  )

X3 <- select(OccupDF2, Heroes_Mean, Danger_Mean, Help_Mean)
mcd3 <- covMcd(X3, alpha = 0.10)

## Warning in covMcd(X3, alpha = 0.1): subsample size h < n/2 may be too small

rdist3 <- mahalanobis(X3, center = mcd3$center, cov = mcd3$cov)

# add to df
OccupDF2$rdist3 <- rdist3



plot_ly(
  data = OccupDF2,
  x    = ~Heroes_Mean,
  y    = ~Danger_Mean,
  z    = ~Help_Mean,
  color = ~rdist3,                # continuous colour scale
  colors = c("steelblue","gold","firebrick"),
  text = ~Occupation,
  hovertemplate = paste(
    "Job: %{text}<br>",
    "Heroes: %{x:.2f}<br>",
    "Danger: %{y:.2f}<br>",
    "Helpfulness: %{z:.2f}<br>",
    "MCD‐dist: %{marker.color:.1f}<extra></extra>"
  ),
  type = "scatter3d",
  mode = "markers"
) %>%
  layout(
    scene = list(
      xaxis = list(title = "Heroes_Mean"),
      yaxis = list(title = "Exposure to danger"),
      zaxis = list(title = "Helpfulness")
    )
  )

---