In previous study, we found that bravery, selflessness, exposure to
danger, and helpfulness were important predictors of heroism in 6
different occupations: Soldiers, Nurses, Underwater welders,
Psychiatrists, Firefighters, and Police officers. This would suggest a
‘universal’ frame for heroes, defined by “risk + altruism”. These two
components might be sufficient (though maybe not necessary) across
occupations.
In the present study, we are interested in possible exceptions for
this frame. Mapping a large and representative subset of occupations on
these heroism components will be a valuable resources for our future
studies. Moreover, mapping the occupations and assessing how some
occupations might escape the typical and previously observed ‘risk +
altruism’ frame of heroism, will provide insights into the essence of
heroism.
Using a recent dataset where a large representative British sample
evaluated the social value and occupational prestige of more than 500
occupations (Newlands & Lutz, 2024), we sub-sampled a representative
subset of occupations. To do so, we categorised each occupation into
their respective ISCO major group, and quintilised each major group on
their social value. We then picked one occupation per quintile per major
group. This resulted in a selection of 53 occupations (see full details
on the OSF webpage of the project, https://osf.io/fj7pw/?view_only=2751ad848d4647348c677d09a960942e).
A representative sample of 430 UK residents will indicate to what
extent each occupation in a list of 54 items are qualified by the
following attributes: Heroic, Brave, Selfless, Helpful, Exposed to
danger (all evaluations using 7-points Likert scales from ‘1 - Not at
all’ to ‘7 - Very much’).
Participants will then report their general attitude toward each of
the 54 target occupations (using a 7-points Likert scale ranging from ‘1
- Very negative’ to ‘7 - Very positive’).
Using mixed models, we will assess how each component of heroism
influence heroism when accounting for occupation type as a random
effect. We will plot the correlations between each key component and
assess any atypical data points.
Data preparation
We load required packages.
if(!require("dplyr")) install.packages("dplyr")
if(!require("MASS")) install.packages("MASS")
if(!require("Matrix")) install.packages("Matrix")
if(!require("ggplot2")) install.packages("ggplot2")
if(!require("robustbase")) install.packages("robustbase") # For Minimum Covariance distance
if(!require("ggrepel")) install.packages("ggrepel")
if(!require("ppcor")) install.packages("ppcor") # For partial correlations
if(!require("corrplot")) install.packages("corrplot")
if(!require("plotly")) install.packages("plotly") # Interactive plots <3
if(!require("performance")) install.packages("performance")
if(!require("tidyr")) install.packages("tidyr")
if(!require("partR2")) install.packages("partR2")
if(!require("report")) install.packages("report")
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("interactions")) install.packages("interactions")
if(!require("gridExtra")) install.packages("gridExtra") # Plot layout
if(!require("scatterplot3d")) install.packages("scatterplot3d") # Really cool 3D scatterplots
Some data wrangling is required to end up with a nice df following
long format: an occupation column, and one column per main variable.
As always, there are participants who appears in the dataset, but
failed attention checks and/or returned their surveys.
CorrSet <- read.csv("~/Downloads/50+shades+of+heroism_July+2,+2025_09.30.csv", comment.char="#")
#head(CorrSet)
CorrSet <- CorrSet[ -c(1,2),]
CorrSet<- subset(CorrSet, CorrSet$Gender != "")
Test <- read.csv("~/Downloads/prolific_export_684feb7d4ef7e1096f146ef9(2).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)]
#nrow(CorrSet)
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 %>%
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
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
mutate(
Occupation = factor(
as.integer(OccCode),
levels = 1:54,
labels = occupation_labels
)
) %>%
# 6. rename columns to exactly what you want
transmute(
Participant = ResponseId,
Occupation,
Heroes = Hero,
Bravery = Brav,
Selflessness = Selfle,
Danger = Dang,
Help = Helpf,
Attitude = Attitude
)
DF<-long_df
Main Analyses
We registered mixed models predicting Heroism using each variable
independtly.
Bravery
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
)
model <- lmerTest::lmer(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: Heroes_num ~ zBravery + zAttitude + (zBravery | Occupation) +
## (1 | Participant)
## Data: DF
##
## REML criterion at convergence: 67373.5
##
## 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.80167 0.8954
## Occupation (Intercept) 0.32929 0.5738
## zBravery 0.07216 0.2686 0.83
## Residual 1.03202 1.0159
## Number of obs: 22790, groups: Participant, 430; Occupation, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.705e+00 9.028e-02 8.477e+01 29.96 <2e-16 ***
## zBravery 5.219e-01 3.849e-02 5.445e+01 13.56 <2e-16 ***
## zAttitude 5.727e-01 1.076e-02 2.065e+04 53.21 <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
## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## -----------------------------------------------------------------------------------
## 67389.465 | 67389.472 | 67453.738 | 0.673 | 0.292 | 0.538 | 1.004 | 1.016
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zBravery and zAttitude (formula: 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 2.71 (95% CI [2.53, 2.88], t(22782) =
## 29.96, p < .001). Within this model:
##
## - The effect of zBravery is statistically significant and positive (beta =
## 0.52, 95% CI [0.45, 0.60], 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.57, 95% CI [0.55, 0.59], 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.
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.538
## Unadjusted ICC: 0.381
We can use ICC to see which occupation is associated to large or
small Bravery -> Heroism slopes:
slopes_df <- ranef(model)$Occupation %>%
tibble::rownames_to_column("Occupation") %>%
rename(slope = zBravery)
p <- 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(p, tooltip = "text")
Danger
model <- lmerTest::lmer(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: Heroes_num ~ zDanger + zAttitude + (zDanger | Occupation) + (1 |
## Participant)
## Data: DF
##
## REML criterion at convergence: 68336.9
##
## 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.73102 0.8550
## Occupation (Intercept) 0.39108 0.6254
## zDanger 0.05874 0.2424 0.48
## Residual 1.07828 1.0384
## Number of obs: 22790, groups: Participant, 430; Occupation, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.763e+00 9.592e-02 7.688e+01 28.80 <2e-16 ***
## zDanger 4.615e-01 3.578e-02 5.497e+01 12.90 <2e-16 ***
## zAttitude 6.753e-01 1.042e-02 2.256e+04 64.83 <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
## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## -----------------------------------------------------------------------------------
## 68352.894 | 68352.900 | 68417.166 | 0.654 | 0.275 | 0.523 | 1.027 | 1.038
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zDanger and zAttitude (formula: 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 2.76 (95% CI [2.57, 2.95], t(22782) =
## 28.80, p < .001). Within this model:
##
## - The effect of zDanger is statistically significant and positive (beta = 0.46,
## 95% CI [0.39, 0.53], 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.68, 95% CI [0.65, 0.70], 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.
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.523
## Unadjusted ICC: 0.379
We can use ICC to see which occupation is associated to large or
small Bravery -> Heroism slopes:
slopes_df <- ranef(model)$Occupation %>%
tibble::rownames_to_column("Occupation") %>%
rename(slope = zDanger)
p <- 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(p, tooltip = "text")
Selflessness
model <- lmerTest::lmer(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:
## Heroes_num ~ zSelflessness + zAttitude + (zSelflessness | Occupation) +
## (1 | Participant)
## Data: DF
##
## REML criterion at convergence: 68326.7
##
## 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.83929 0.9161
## Occupation (Intercept) 0.50401 0.7099
## zSelflessness 0.06129 0.2476 0.84
## Residual 1.07551 1.0371
## Number of obs: 22790, groups: Participant, 430; Occupation, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.722e+00 1.074e-01 7.431e+01 25.34 < 2e-16 ***
## zSelflessness 3.878e-01 3.588e-02 5.641e+01 10.81 2.31e-15 ***
## zAttitude 5.985e-01 1.139e-02 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
## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## -----------------------------------------------------------------------------------
## 68342.707 | 68342.714 | 68406.980 | 0.673 | 0.247 | 0.566 | 1.025 | 1.037
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zSelflessness and zAttitude (formula: 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 2.72
## (95% CI [2.51, 2.93], t(22782) = 25.34, p < .001). Within this model:
##
## - The effect of zSelflessness is statistically significant and positive (beta =
## 0.39, 95% CI [0.32, 0.46], 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.60, 95% CI [0.58, 0.62], 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.
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.566
## Unadjusted ICC: 0.427
slopes_df <- ranef(model)$Occupation %>%
tibble::rownames_to_column("Occupation") %>%
rename(slope = zSelflessness)
p <- 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(p, tooltip = "text")
Helpfulness
model <- lmerTest::lmer(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: Heroes_num ~ zHelp + zAttitude + (zHelp | Occupation) + (1 |
## Participant)
## Data: DF
##
## REML criterion at convergence: 67490.1
##
## 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.75098 0.8666
## Occupation (Intercept) 0.40601 0.6372
## zHelp 0.09119 0.3020 0.90
## Residual 1.03874 1.0192
## Number of obs: 22790, groups: Participant, 430; Occupation, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.700e+00 9.741e-02 7.686e+01 27.72 <2e-16 ***
## zHelp 5.327e-01 4.325e-02 5.440e+01 12.31 <2e-16 ***
## zAttitude 5.232e-01 1.133e-02 1.953e+04 46.20 <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
## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## -----------------------------------------------------------------------------------
## 67506.132 | 67506.139 | 67570.405 | 0.676 | 0.287 | 0.546 | 1.008 | 1.019
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zHelp and zAttitude (formula: 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 2.70 (95% CI [2.51, 2.89], t(22782) = 27.72, p <
## .001). Within this model:
##
## - The effect of zHelp is statistically significant and positive (beta = 0.53,
## 95% CI [0.45, 0.62], 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.52, 95% CI [0.50, 0.55], 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.
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.546
## Unadjusted ICC: 0.389
slopes_df <- ranef(model)$Occupation %>%
tibble::rownames_to_column("Occupation") %>%
rename(slope = zHelp)
p <- 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 Selflessness",
y = "Occupation"
)
# Convert to an interactive Plotly plot
ggplotly(p, tooltip = "text")
Interactions
We might explore interactions (see additional comments of the
registration).
There are two possible models of interest: - Interaction between
character evaluations: Bravery x Selflessness - Interaction between
situation evaluations: Danger x Helpfulness
Character interactions: Bravery x Selflessness
m <- lmerTest::lmer(
Heroes_num ~ zBravery * zSelflessness + zAttitude +
(1 + zBravery * zSelflessness | Occupation) +
(1 | Participant),
data = DF,
REML = TRUE
)
summary(m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Heroes_num ~ zBravery * zSelflessness + zAttitude + (1 + zBravery *
## zSelflessness | Occupation) + (1 | Participant)
## Data: DF
##
## REML criterion at convergence: 66237.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.0466 -0.5944 -0.0555 0.5521 6.8729
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Participant (Intercept) 0.762982 0.87349
## Occupation (Intercept) 0.259734 0.50964
## zBravery 0.024244 0.15570 0.83
## zSelflessness 0.019032 0.13796 0.56 0.21
## zBravery:zSelflessness 0.002855 0.05343 -0.66 -0.36 -0.33
## Residual 0.980165 0.99003
## Number of obs: 22790, groups: Participant, 430; Occupation, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.640e+00 8.216e-02 9.163e+01 32.14 < 2e-16 ***
## zBravery 4.794e-01 2.422e-02 5.170e+01 19.80 < 2e-16 ***
## zSelflessness 2.951e-01 2.220e-02 5.935e+01 13.29 < 2e-16 ***
## zAttitude 4.744e-01 1.124e-02 1.732e+04 42.23 < 2e-16 ***
## zBravery:zSelflessness 1.563e-01 1.096e-02 3.567e+01 14.26 2.65e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) zBrvry zSlfls zAtttd
## zBravery 0.629
## zSelflssnss 0.409 0.102
## zAttitude 0.000 -0.127 -0.180
## zBrvry:zSlf -0.397 -0.174 -0.159 0.042
## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## -----------------------------------------------------------------------------------
## 66271.772 | 66271.799 | 66408.351 | 0.691 | 0.357 | 0.519 | 0.978 | 0.990
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zBravery, zSelflessness and zAttitude (formula:
## Heroes_num ~ zBravery * zSelflessness + zAttitude). The model included zBravery
## as random effects (formula: list(~1 + zBravery * zSelflessness | Occupation, ~1
## | Participant)). The model's total explanatory power is substantial
## (conditional R2 = 0.69) and the part related to the fixed effects alone
## (marginal R2) is of 0.36. The model's intercept, corresponding to zBravery = 0,
## zSelflessness = 0 and zAttitude = 0, is at 2.64 (95% CI [2.48, 2.80], t(22773)
## = 32.14, p < .001). Within this model:
##
## - The effect of zBravery is statistically significant and positive (beta =
## 0.48, 95% CI [0.43, 0.53], t(22773) = 19.80, p < .001; Std. beta = 0.24, 95% CI
## [0.22, 0.26])
## - The effect of zSelflessness is statistically significant and positive (beta =
## 0.30, 95% CI [0.25, 0.34], t(22773) = 13.29, p < .001; Std. beta = 0.15, 95% CI
## [0.13, 0.17])
## - The effect of zAttitude is statistically significant and positive (beta =
## 0.47, 95% CI [0.45, 0.50], t(22773) = 42.22, p < .001; Std. beta = 0.24, 95% CI
## [0.23, 0.25])
## - The effect of zBravery × zSelflessness is statistically significant and
## positive (beta = 0.16, 95% CI [0.13, 0.18], t(22773) = 14.26, p < .001; Std.
## beta = 0.08, 95% CI [0.07, 0.09])
##
## 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.
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.519
## Unadjusted ICC: 0.334
interact_plot(
model = m,
pred = zBravery, # x-axis
modx = zSelflessness, # moderator
plot.points = TRUE, # raw data in the background
interval = TRUE, # 95% CI
int.width = 0.95,
x.label = "Standardized Bravery",
modx.labels = c("−1 SD", "Mean", "+1 SD"),
y.label = "Predicted Heroes_num",
legend.main = "Selflessness"
)
re_df <- ranef(m)$Occupation %>%
tibble::rownames_to_column("Occupation") %>%
select(Occupation, zBravery, zSelflessness, `zBravery:zSelflessness`)
# Create a long-format list of separate dataframes for each predictor
slope_list <- re_df %>%
pivot_longer(-Occupation, names_to = "Predictor", values_to = "Slope") %>%
group_split(Predictor)
plot_list <- lapply(slope_list, function(df) {
predictor_name <- unique(df$Predictor)
gg <- ggplot(df, aes(
x = Slope,
y = reorder(Occupation, Slope),
text = paste0("Occupation: ", Occupation,
"<br>Slope: ", round(Slope, 3))
)) +
geom_point() +
labs(
title = paste("Random Slopes for", predictor_name),
x = "Slope",
y = "Occupation"
) +
theme_minimal()
ggplotly(gg, tooltip = "text")
})
plot_list[[1]]
You can hover over the plot to see the interesting
cases
Situation Interaction: Danger x Helpfulness
m <- lmerTest::lmer(
Heroes_num ~ zDanger * zHelp + zAttitude +
(1 + zDanger * zHelp | Occupation) +
(1 | Participant),
data = DF,
REML = TRUE
)
summary(m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Heroes_num ~ zDanger * zHelp + zAttitude + (1 + zDanger * zHelp |
## Occupation) + (1 | Participant)
## Data: DF
##
## REML criterion at convergence: 66212.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7928 -0.5763 -0.0426 0.5548 6.7488
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Participant (Intercept) 0.658608 0.81155
## Occupation (Intercept) 0.210067 0.45833
## zDanger 0.011989 0.10950 0.60
## zHelp 0.053218 0.23069 0.92 0.43
## zDanger:zHelp 0.008673 0.09313 -0.76 -0.24 -0.81
## Residual 0.982649 0.99129
## Number of obs: 22790, groups: Participant, 430; Occupation, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.665e+00 7.486e-02 9.342e+01 35.60 < 2e-16 ***
## zDanger 3.997e-01 1.972e-02 5.493e+01 20.27 < 2e-16 ***
## zHelp 4.908e-01 3.425e-02 5.408e+01 14.33 < 2e-16 ***
## zAttitude 4.681e-01 1.102e-02 1.210e+04 42.49 < 2e-16 ***
## zDanger:zHelp 1.386e-01 1.627e-02 4.734e+01 8.52 4.07e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) zDangr zHelp zAtttd
## zDanger 0.399
## zHelp 0.722 0.251
## zAttitude 0.002 -0.059 -0.157
## zDangr:zHlp -0.523 -0.147 -0.572 -0.027
## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## -----------------------------------------------------------------------------------
## 66246.616 | 66246.643 | 66383.195 | 0.689 | 0.397 | 0.485 | 0.979 | 0.991
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zDanger, zHelp and zAttitude (formula: Heroes_num ~
## zDanger * zHelp + zAttitude). The model included zDanger as random effects
## (formula: list(~1 + zDanger * zHelp | Occupation, ~1 | Participant)). The
## model's total explanatory power is substantial (conditional R2 = 0.69) and the
## part related to the fixed effects alone (marginal R2) is of 0.40. The model's
## intercept, corresponding to zDanger = 0, zHelp = 0 and zAttitude = 0, is at
## 2.66 (95% CI [2.52, 2.81], t(22773) = 35.60, p < .001). Within this model:
##
## - The effect of zDanger is statistically significant and positive (beta = 0.40,
## 95% CI [0.36, 0.44], t(22773) = 20.27, p < .001; Std. beta = 0.20, 95% CI
## [0.18, 0.22])
## - The effect of zHelp is statistically significant and positive (beta = 0.49,
## 95% CI [0.42, 0.56], t(22773) = 14.33, p < .001; Std. beta = 0.24, 95% CI
## [0.21, 0.28])
## - The effect of zAttitude is statistically significant and positive (beta =
## 0.47, 95% CI [0.45, 0.49], t(22773) = 42.49, p < .001; Std. beta = 0.23, 95% CI
## [0.22, 0.24])
## - The effect of zDanger × zHelp is statistically significant and positive (beta
## = 0.14, 95% CI [0.11, 0.17], t(22773) = 8.52, p < .001; Std. beta = 0.07, 95%
## CI [0.05, 0.08])
##
## 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.
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.485
## Unadjusted ICC: 0.292
interact_plot(
model = m,
pred = zDanger, # x-axis
modx = zHelp, # moderator
plot.points = TRUE, # raw data in the background
interval = TRUE, # 95% CI
int.width = 0.95,
x.label = "Standardized Bravery",
modx.labels = c("−1 SD", "Mean", "+1 SD"),
y.label = "Predicted Heroes_num",
legend.main = "Selflessness"
)
re_df <- ranef(m)$Occupation %>%
tibble::rownames_to_column("Occupation") %>%
select(Occupation, zDanger, zHelp, `zDanger:zHelp`)
# Create a long-format list of separate dataframes for each predictor
slope_list <- re_df %>%
pivot_longer(-Occupation, names_to = "Predictor", values_to = "Slope") %>%
group_split(Predictor)
plot_list <- lapply(slope_list, function(df) {
predictor_name <- unique(df$Predictor)
gg <- ggplot(df, aes(
x = Slope,
y = reorder(Occupation, Slope),
text = paste0("Occupation: ", Occupation,
"<br>Slope: ", round(Slope, 3))
)) +
geom_point() +
labs(
title = paste("Random Slopes for", predictor_name),
x = "Slope",
y = "Occupation"
) +
theme_minimal()
ggplotly(gg, tooltip = "text")
})
plot_list[[1]]
Comparing components directly
In a single model, we predict heroism using all components:
m <- lmerTest::lmer(
Heroes_num ~ zDanger + zHelp + zBravery + zSelflessness + zAttitude +
(1 | Occupation) +
(1 | Participant),
data = DF,
REML = TRUE
)
summary(m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Heroes_num ~ zDanger + zHelp + zBravery + zSelflessness + zAttitude +
## (1 | Occupation) + (1 | Participant)
## Data: DF
##
## REML criterion at convergence: 66075.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.8944 -0.6033 -0.0246 0.5695 6.9479
##
## Random effects:
## Groups Name Variance Std.Dev.
## Participant (Intercept) 0.6495 0.8059
## Occupation (Intercept) 0.2676 0.5173
## Residual 0.9816 0.9907
## Number of obs: 22790, groups: Participant, 430; Occupation, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.846e+00 8.126e-02 8.544e+01 35.02 <2e-16 ***
## zDanger 2.736e-01 1.221e-02 2.215e+04 22.41 <2e-16 ***
## zHelp 3.752e-01 1.219e-02 2.272e+04 30.77 <2e-16 ***
## zBravery 2.943e-01 1.101e-02 2.261e+04 26.74 <2e-16 ***
## zSelflessness 2.123e-01 1.114e-02 2.267e+04 19.05 <2e-16 ***
## zAttitude 3.601e-01 1.182e-02 2.263e+04 30.47 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) zDangr zHelp zBrvry zSlfls
## zDanger 0.000
## zHelp 0.000 -0.100
## zBravery 0.000 -0.274 -0.113
## zSelflssnss 0.000 -0.010 -0.162 -0.187
## zAttitude 0.000 -0.020 -0.313 -0.187 -0.293
## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## -----------------------------------------------------------------------------------
## 66093.512 | 66093.519 | 66165.818 | 0.703 | 0.426 | 0.483 | 0.980 | 0.991
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
## to predict Heroes_num with zDanger, zHelp, zBravery, zSelflessness and
## zAttitude (formula: Heroes_num ~ zDanger + zHelp + zBravery + zSelflessness +
## zAttitude). The model included Occupation as random effects (formula: list(~1 |
## Occupation, ~1 | Participant)). The model's total explanatory power is
## substantial (conditional R2 = 0.70) and the part related to the fixed effects
## alone (marginal R2) is of 0.43. The model's intercept, corresponding to zDanger
## = 0, zHelp = 0, zBravery = 0, zSelflessness = 0 and zAttitude = 0, is at 2.85
## (95% CI [2.69, 3.01], t(22781) = 35.03, p < .001). Within this model:
##
## - The effect of zDanger is statistically significant and positive (beta = 0.27,
## 95% CI [0.25, 0.30], t(22781) = 22.40, p < .001; Std. beta = 0.14, 95% CI
## [0.12, 0.15])
## - The effect of zHelp is statistically significant and positive (beta = 0.38,
## 95% CI [0.35, 0.40], t(22781) = 30.77, p < .001; Std. beta = 0.19, 95% CI
## [0.18, 0.20])
## - The effect of zBravery is statistically significant and positive (beta =
## 0.29, 95% CI [0.27, 0.32], t(22781) = 26.74, p < .001; Std. beta = 0.15, 95% CI
## [0.14, 0.16])
## - The effect of zSelflessness is statistically significant and positive (beta =
## 0.21, 95% CI [0.19, 0.23], t(22781) = 19.05, p < .001; Std. beta = 0.11, 95% CI
## [0.09, 0.12])
## - The effect of zAttitude is statistically significant and positive (beta =
## 0.36, 95% CI [0.34, 0.38], t(22781) = 30.47, p < .001; Std. beta = 0.18, 95% CI
## [0.17, 0.19])
##
## 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.
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.483
## Unadjusted ICC: 0.277
# 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)
## [1] "Parameter" "Coefficient" "SE" "CI" "CI_low"
## [6] "CI_high" "t" "df_error" "p" "Effects"
# 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"
)
Fixed‐effect estimates from mixed model
| (Intercept) |
0.00 |
0.00 |
22781 |
-0.08 |
0.08 |
| zDanger |
0.14 |
22.40 |
22781 |
0.12 |
0.15 |
| zHelp |
0.19 |
30.77 |
22781 |
0.18 |
0.20 |
| zBravery |
0.15 |
26.74 |
22781 |
0.14 |
0.16 |
| zSelflessness |
0.11 |
19.05 |
22781 |
0.09 |
0.12 |
| zAttitude |
0.18 |
30.47 |
22781 |
0.17 |
0.19 |
Main correlations
Of importance: observe the discrepant cases where help and danger or
bravery and selflessness do not work together.
Not accounting for attitude.
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
)
# 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")
## 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 <- 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"
Selflessness:Heroism
paste0("Correlation = ", round(cor(OccupDF2$mean_zSelfl, OccupDF2$mean_zHeroism, use = "complete.obs"), 2))
## [1] "Correlation = 0.87"
Danger:Heroism
paste0("Correlation = ", round(cor(OccupDF2$mean_zDang, OccupDF2$mean_zHeroism, use = "complete.obs"), 2))
## [1] "Correlation = 0.69"
Helpfulness:Heroism
paste0("Correlation = ", round(cor(OccupDF2$mean_zHelp, OccupDF2$mean_zHeroism, use = "complete.obs"), 2))
## [1] "Correlation = 0.86"
Selflessness:Bravery
paste0("Correlation = ", round(cor(OccupDF2$mean_zBrav, OccupDF2$mean_zSelfl, use = "complete.obs"), 2))
## [1] "Correlation = 0.79"
Danger:Bravery
paste0("Correlation = ", round(cor(OccupDF2$mean_zDang, OccupDF2$mean_zBrav, use = "complete.obs"), 2))
## [1] "Correlation = 0.86"
Helpfulness:Bravery
paste0("Correlation = ", round(cor(OccupDF2$mean_zBrav, OccupDF2$mean_zHelp, use = "complete.obs"), 2))
## [1] "Correlation = 0.77"
Danger:Selflessness
paste0("Correlation = ", round(cor(OccupDF2$mean_zDang, OccupDF2$mean_zSelfl, use = "complete.obs"), 2))
## [1] "Correlation = 0.44"
Helpfulness:Selflessness
paste0("Correlation = ", round(cor(OccupDF2$mean_zHelp, OccupDF2$mean_zSelfl, use = "complete.obs"), 2))
## [1] "Correlation = 0.91"
Helpfulness:Danger
paste0("Correlation = ", round(cor(OccupDF2$mean_zHelp, OccupDF2$mean_zDang, use = "complete.obs"), 2))
## [1] "Correlation = 0.48"
Residual approach
Scatterplots for each pair of variable when accounting
for attitude.
In the previous approach, we used raw scores. However, to account for
a possible Halo effect, we can use residual
scores after regressing our variables on Attitude.
Note than all the plots below use residual scores after partialing
out Attitude. So we are using attitude-free outcomes. This can results
in some weird results for occupations that have ceiling
(e.g. firefighters) or floor (e.g., Organised Crime Boss) attitudes.
# Load required packages
library(dplyr)
library(ggplot2)
library(plotly)
library(RColorBrewer)
# 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"))
}
# 4) Create your 10 residual-based interactive plots
p1_resid <- make_resid_plot(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")
p2_resid <- make_resid_plot(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_plot(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_plot(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_plot(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_plot(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_plot(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_plot(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_plot(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_plot(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")
Bravery:Heroism
paste0("Correlation = ", round(cor(OccupDF3$resid_zBrav, OccupDF3$resid_zHeroism, use = "complete.obs"), 2))
## [1] "Correlation = 0.81"
Selflessness:Heroism
paste0("Correlation = ", round(cor(OccupDF3$resid_zSelfl, OccupDF3$resid_zHeroism, use = "complete.obs"), 2))
## [1] "Correlation = 0.52"
Danger:Heroism
paste0("Correlation = ", round(cor(OccupDF3$resid_zDang, OccupDF3$resid_zHeroism, use = "complete.obs"), 2))
## [1] "Correlation = 0.72"
Helpfulness:Heroism
paste0("Correlation = ", round(cor(OccupDF3$resid_zHelp, OccupDF3$resid_zHeroism, use = "complete.obs"), 2))
## [1] "Correlation = 0.47"
Selflessness:Bravery
paste0("Correlation = ", round(cor(OccupDF3$resid_zBrav, OccupDF3$resid_zSelfl, use = "complete.obs"), 2))
## [1] "Correlation = 0.39"
Danger:Bravery
paste0("Correlation = ", round(cor(OccupDF3$resid_zBrav, OccupDF3$resid_zDang, use = "complete.obs"), 2))
## [1] "Correlation = 0.94"
Helpfulness:Bravery
paste0("Correlation = ", round(cor(OccupDF3$resid_zBrav, OccupDF3$resid_zHelp, use = "complete.obs"), 2))
## [1] "Correlation = 0.33"
Danger:Selflessness
paste0("Correlation = ", round(cor(OccupDF3$resid_zDang, OccupDF3$resid_zSelfl, use = "complete.obs"), 2))
## [1] "Correlation = 0.3"
Helpfulness:Selflessness
paste0("Correlation = ", round(cor(OccupDF3$resid_zHelp, OccupDF3$resid_zSelfl, use = "complete.obs"), 2))
## [1] "Correlation = 0.33"
Helpfulness:Danger
paste0("Correlation = ", round(cor(OccupDF3$resid_zHelp, OccupDF3$resid_zDang, use = "complete.obs"), 2))
## [1] "Correlation = 0.35"
3D viz
- Finally, 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")
)
)
Conclusion
Do our components of heroism ‘work’ across a
representative set of occupations?
==> YES. Particularly character evaluations. But the best
predictors are Bravery + Helpfulness
Any exceptions to these rules?
==> Very few ‘ambivalent cases’ when using raw scores – but
controlling for attitudes does bring out some cases: Mafia boss defies
our model on the danger dimension, so does all ‘villain’ occupation.
Some Interesting jobs: Charity worker, caretaker, home maker
(not particularly brave, but helpful), criminals, machine
operator/roofers (brave, but not helpful).
For any question relating to this project or this script, please
contact j.moneger@kent.ac.uk