Analysis

Analysis Priors

prior_m1_interaction_gen <- c(
  prior(normal(0, 1), class = "Intercept"),
  prior(normal(-2.85, 0.98), class = "b", coef = "Age"),
  prior(normal(3.00, 0.98), class = "b", coef = "dominance_Sum"),
  prior(normal(-0.02, 0.98), class = "b", coef = "dominance_Sum:Gender2"),
  prior(normal(-3.02, 1.00), class = "b", coef = "Gender2"),
  prior(normal(0.05, 0.98), class = "b", coef = "Gender2:leadership_Sum"),
  prior(normal(0.05, 0.99), class = "b", coef = "Gender2:prestige_Sum"),
  prior(normal(-1.91, 0.98), class = "b", coef = "leadership_Sum"),
  prior(normal(0.09, 0.98), class = "b", coef = "prestige_Sum"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum"),
  prior(normal(0, 1), class = "b", coef = "Gender2:PNI_Sum")
)

prior_m1_interaction_gen_no_pni <- c(
  prior(normal(-2.85, 0.98), class = "b", coef = "Age"),
  prior(normal(3.00, 0.98), class = "b", coef = "dominance_Sum"),
  prior(normal(-0.02, 0.98), class = "b", coef = "dominance_Sum:Gender2"),
  prior(normal(-3.02, 1.00), class = "b", coef = "Gender2"),
  prior(normal(0.05, 0.98), class = "b", coef = "Gender2:leadership_Sum"),
  prior(normal(0.05, 0.99), class = "b", coef = "Gender2:prestige_Sum"),
  prior(normal(-1.91, 0.98), class = "b", coef = "leadership_Sum"),
  prior(normal(0.09, 0.98), class = "b", coef = "prestige_Sum")
)

prior_m1 <- c(
  prior(normal(0, 1), class = "Intercept"),
  prior(normal(-2.86, 0.99), class = "b", coef = "Age"),
  prior(normal(3.00, 0.99), class = "b", coef = "dominance_Sum"),
  prior(normal(-3.02, 0.99), class = "b", coef = "Gender2"),
  prior(normal(-1.91, 0.99), class = "b", coef = "leadership_Sum"),
  prior(normal(0.09, 0.99), class = "b", coef = "prestige_Sum"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum")
)

prior_m2 <- c(
  # Risk
  prior(normal(-0.121402797, 0.317918624), class = "Intercept", resp = "riskSumz"),
  prior(normal(0, 1), class = "sigma", resp = "riskSumz"),
  prior(normal(0.39426886, 0.097544417), class = "b", coef = "dominance_Sum", resp = "riskSumz"),
  prior(normal(0.168132287, 0.097451296), class = "b", coef = "prestige_Sum", resp = "riskSumz"),
  prior(normal(-0.009328966, 0.09515019), class = "b", coef = "leadership_Sum", resp = "riskSumz"),
  prior(normal(0.000831313, 0.010064451), class = "b", coef = "Age", resp = "riskSumz"),
  prior(normal(0.200477395, 0.179578595), class = "b", coef = "Gender2", resp = "riskSumz"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "riskSumz"),
  # Risk Perception
  prior(normal(0.541755102, 0.343784152), class = "Intercept", resp = "riskPerceptionSumz"),
  prior(normal(0, 1), class = "sigma", resp = "riskPerceptionSumz"),
  prior(normal(-0.270748615, 0.108267928), class = "b", coef = "dominance_Sum", resp = "riskPerceptionSumz"),
  prior(normal(0.210735032, 0.108341126), class = "b", coef = "prestige_Sum", resp = "riskPerceptionSumz"),
  prior(normal(0.031158527, 0.105772838), class = "b", coef = "leadership_Sum", resp = "riskPerceptionSumz"),
  prior(normal(-0.017714536, 0.010901932), class = "b", coef = "Age", resp = "riskPerceptionSumz"),
  prior(normal(-0.187998634, 0.194431114), class = "b", coef = "Gender2", resp = "riskPerceptionSumz"),
  prior(normal(-0.0, 1), class = "b", coef = "PNI_Sum", resp = "riskPerceptionSumz"),
  # Risk Benefit
  prior(normal(-0.373633372, 0.33954029), class = "Intercept", resp = "riskBenefitSumz"),
  prior(normal(0, 1), class = "sigma", resp = "riskBenefitSumz"),
  prior(normal(0.248770622, 0.104869086), class = "b", coef = "dominance_Sum", resp = "riskBenefitSumz"),
  prior(normal(0.142274426, 0.104666143), class = "b", coef = "prestige_Sum", resp = "riskBenefitSumz"),
  prior(normal(-0.025816283, 0.102197521), class = "b", coef = "leadership_Sum", resp = "riskBenefitSumz"),
  prior(normal(0.007281249, 0.010703955), class = "b", coef = "Age", resp = "riskBenefitSumz"),
  prior(normal(0.297908079, 0.194394947), class = "b", coef = "Gender2", resp = "riskBenefitSumz"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "riskBenefitSumz")
)


prior_m2_interaction_gender <- c(
  prior(normal(0.211359466039411, 0.122792310103859), class = "Intercept", resp = "riskSumz"),
  prior(normal(-0.233748907829191, 0.131188637761139), class = "Intercept", resp = "riskPerceptionSumz"),
  prior(normal(0.249317104065427, 0.133157381895898), class = "Intercept", resp = "riskBenefitSumz"),
  prior(normal(0.653445259241636, 0.149635216859058), class = "b", coef = "dominance_Sum", resp = "riskSumz"),
  prior(normal(-0.495103273346896, 0.177473117281467), class = "b", coef = "Gender2", resp = "riskSumz"),
  prior(normal(0.0990458133895065, 0.140168789301184), class = "b", coef = "prestige_Sum", resp = "riskSumz"),
  prior(normal(-0.0285312341213451, 0.134246520912695), class = "b", coef = "leadership_Sum", resp = "riskSumz"),
  prior(normal(0.0162970878783793, 0.0908982500875085), class = "b", coef = "Age", resp = "riskSumz"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "riskSumz"),
  prior(normal(-0.47831593537673, 0.190465828006423), class = "b", coef = "dominance_Sum:Gender2", resp = "riskSumz"),
  prior(normal(0.095717435377145, 0.188962898776915), class = "b", coef = "Gender2:prestige_Sum", resp = "riskSumz"),
  prior(normal(0.0358118424670447, 0.182068107768659), class = "b", coef = "Gender2:leadership_Sum", resp = "riskSumz"),
  prior(normal(0, 1), class = "b", coef = "Gender2:PNI_Sum", resp = "riskSumz"),
  prior(normal(-0.226410040336677, 0.163334625155952), class = "b", coef = "dominance_Sum", resp = "riskPerceptionSumz"),
  prior(normal(0.426961907181259, 0.191507665064676), class = "b", coef = "Gender2", resp = "riskPerceptionSumz"),
  prior(normal(0.305553399549671, 0.152964977686876), class = "b", coef = "prestige_Sum", resp = "riskPerceptionSumz"),
  prior(normal(-0.191711084170981, 0.147449563411409), class = "b", coef = "leadership_Sum", resp = "riskPerceptionSumz"),
  prior(normal(-0.170813512027718, 0.0989265773273746), class = "b", coef = "Age", resp = "riskPerceptionSumz"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "riskPerceptionSumz"),
  prior(normal(-0.055455232784238, 0.206139416132584), class = "b", coef = "dominance_Sum:Gender2", resp = "riskPerceptionSumz"),
  prior(normal(-0.192454377436851, 0.20943763365537), class = "b", coef = "Gender2:prestige_Sum", resp = "riskPerceptionSumz"),
  prior(normal(0.431973601346632, 0.200374884521192), class = "b", coef = "Gender2:leadership_Sum", resp = "riskPerceptionSumz"),
  prior(normal(0, 1), class = "b", coef = "Gender2:PNI_Sum", resp = "riskPerceptionSumz"),
  prior(normal(0.384554988587524, 0.162183757747619), class = "b", coef = "dominance_Sum", resp = "riskBenefitSumz"),
  prior(normal(-0.595963440527305, 0.190748625267443), class = "b", coef = "Gender2", resp = "riskBenefitSumz"),
  prior(normal(0.0215234988730571, 0.151796284179672), class = "b", coef = "prestige_Sum", resp = "riskBenefitSumz"),
  prior(normal(-0.028365512364167, 0.145804876707747), class = "b", coef = "leadership_Sum", resp = "riskBenefitSumz"),
  prior(normal(0.0737186614972512, 0.0976640041773247), class = "b", coef = "Age", resp = "riskBenefitSumz"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "riskBenefitSumz"),
  prior(normal(-0.26549479719041, 0.208575601040728), class = "b", coef = "dominance_Sum:Gender2", resp = "riskBenefitSumz"),
  prior(normal(0.217074362051788, 0.206053971873005), class = "b", coef = "Gender2:prestige_Sum", resp = "riskBenefitSumz"),
  prior(normal(-0.0022324287109488, 0.197145691092869), class = "b", coef = "Gender2:leadership_Sum", resp = "riskBenefitSumz"),
  prior(normal(0, 1), class = "b", coef = "Gender2:PNI_Sum", resp = "riskBenefitSumz")
)

prior_m3 <- c(
  prior(normal(3.61135311921322, 0.398349216512441), class = "Intercept", resp = "ethicalPreferencez"),
  prior(normal(8.60471748613202, 0.565601362688251), class = "Intercept", resp = "financialPreferencez"),
  prior(normal(10.024592015549, 0.877847943744065), class = "Intercept", resp = "socialPreferencez"),
  prior(normal(5.62829428110632, 0.486128699758622), class = "Intercept", resp = "healthAndSafetyPreferencez"),
  prior(normal(1.66534684425098, 0.395641365631446), class = "Intercept", resp = "recreationalPreferencez"),
  prior(normal(0, 1), class = "sigma", resp = "ethicalPreferencez"),
  prior(normal(0, 1), class = "sigma", resp = "financialPreferencez"),
  prior(normal(0, 1), class = "sigma", resp = "socialPreferencez"),
  prior(normal(0, 1), class = "sigma", resp = "healthAndSafetyPreferencez"),
  prior(normal(0, 1), class = "sigma", resp = "recreationalPreferencez"),
  prior(normal(1.01820236866574, 0.230324886400241), class = "b", coef = "dominance_Sum", resp = "ethicalPreferencez"),
  prior(normal(-0.0448056408845481, 0.227661535353379), class = "b", coef = "prestige_Sum", resp = "ethicalPreferencez"),
  prior(normal(0.0644149801605582, 0.220754430664067), class = "b", coef = "leadership_Sum", resp = "ethicalPreferencez"),
  prior(normal(-0.0807925790683574, 0.342054140420214), class = "b", coef = "Gender2", resp = "ethicalPreferencez"),
  prior(normal(0.257535398706787, 0.213176289264519), class = "b", coef = "Age", resp = "ethicalPreferencez"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "ethicalPreferencez"),
  prior(normal(0.770897209832422, 0.323196910570939), class = "b", coef = "dominance_Sum", resp = "financialPreferencez"),
  prior(normal(0.218932574933387, 0.312178873831866), class = "b", coef = "prestige_Sum", resp = "financialPreferencez"),
  prior(normal(-0.0527566662392789, 0.30751451857796), class = "b", coef = "leadership_Sum", resp = "financialPreferencez"),
  prior(normal(-0.122615013279046, 0.465191008236103), class = "b", coef = "Gender2", resp = "financialPreferencez"),
  prior(normal(-0.0756269171084713, 0.30176572695622), class = "b", coef = "Age", resp = "financialPreferencez"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "financialPreferencez"),
  prior(normal(1.45607644842891, 0.521662114182799), class = "b", coef = "dominance_Sum", resp = "socialPreferencez"),
  prior(normal(0.490922123303897, 0.513950094547241), class = "b", coef = "prestige_Sum", resp = "socialPreferencez"),
  prior(normal(0.178147050305328, 0.516402464149945), class = "b", coef = "leadership_Sum", resp = "socialPreferencez"),
  prior(normal(-0.540788557994253, 0.708620378346552), class = "b", coef = "Gender2", resp = "socialPreferencez"),
  prior(normal(-0.130724948888778, 0.499035382854254), class = "b", coef = "Age", resp = "socialPreferencez"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "socialPreferencez"),
  prior(normal(1.09239737499994, 0.281490924431951), class = "b", coef = "dominance_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(-0.00413639061471685, 0.278414215040186), class = "b", coef = "prestige_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(-0.0730575993095214, 0.277181969750254), class = "b", coef = "leadership_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(0.0371780082169458, 0.409955298320128), class = "b", coef = "Gender2", resp = "healthAndSafetyPreferencez"),
  prior(normal(0.290053188796595, 0.26225865354914), class = "b", coef = "Age", resp = "healthAndSafetyPreferencez"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(1.20303099592023, 0.229344906845703), class = "b", coef = "dominance_Sum", resp = "recreationalPreferencez"),
  prior(normal(-0.288711309452481, 0.223487787997842), class = "b", coef = "prestige_Sum", resp = "recreationalPreferencez"),
  prior(normal(-0.103324574146601, 0.220311845696074), class = "b", coef = "leadership_Sum", resp = "recreationalPreferencez"),
  prior(normal(-1.12995776527766, 0.341594265540998), class = "b", coef = "Gender2", resp = "recreationalPreferencez"),
  prior(normal(0.443085112592513, 0.209303734303984), class = "b", coef = "Age", resp = "recreationalPreferencez"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "recreationalPreferencez")
)

prior_m3_interaction_gender <- c(
  prior(normal(3.60556936655962, 0.396580813334836), class = "Intercept", resp = "ethicalPreferencez"),
  prior(normal(8.60016846571163, 0.55989827635477), class = "Intercept", resp = "financialPreferencez"),
  prior(normal(9.98313800859403, 0.870783002602423), class = "Intercept", resp = "socialPreferencez"),
  prior(normal(5.5993299179946, 0.494545763665189), class = "Intercept", resp = "healthAndSafetyPreferencez"),
  prior(normal(1.67501584950533, 0.395261603649963), class = "Intercept", resp = "recreationalPreferencez"),
  prior(normal(0, 1), class = "sigma", resp = "ethicalPreferencez"),
  prior(normal(0, 1), class = "sigma", resp = "financialPreferencez"),
  prior(normal(0, 1), class = "sigma", resp = "socialPreferencez"),
  prior(normal(0, 1), class = "sigma", resp = "healthAndSafetyPreferencez"),
  prior(normal(0, 1), class = "sigma", resp = "recreationalPreferencez"),
  prior(normal(1.15219692310003, 0.280496121652513), class = "b", coef = "dominance_Sum", resp = "ethicalPreferencez"),
  prior(normal(-0.0904502715842434, 0.337483849502847), class = "b", coef = "Gender2", resp = "ethicalPreferencez"),
  prior(normal(-0.22598791634437, 0.273713173932402), class = "b", coef = "prestige_Sum", resp = "ethicalPreferencez"),
  prior(normal(0.0308646711551821, 0.258538873662647), class = "b", coef = "leadership_Sum", resp = "ethicalPreferencez"),
  prior(normal(0.270330442109019, 0.212089881486683), class = "b", coef = "Age", resp = "ethicalPreferencez"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "ethicalPreferencez"),
  prior(normal(-0.301180433615289, 0.340641941157536), class = "b", coef = "dominance_Sum:Gender2", resp = "ethicalPreferencez"),
  prior(normal(0.402692272506763, 0.342934073772113), class = "b", coef = "Gender2:prestige_Sum", resp = "ethicalPreferencez"),
  prior(normal(0, 1), class = "b", coef = "Gender2:leadership_Sum", resp = "ethicalPreferencez"),
  prior(normal(0.0502639776850436, 0.329985134655418), class = "b", coef = "Gender2:PNI_Sum", resp = "ethicalPreferencez"),
  prior(normal(0.865622123607589, 0.372240213681623), class = "b", coef = "dominance_Sum", resp = "financialPreferencez"),
  prior(normal(-0.127788500245137, 0.463896897597577), class = "b", coef = "Gender2", resp = "financialPreferencez"),
  prior(normal(0.0311011331505854, 0.36633048162564), class = "b", coef = "prestige_Sum", resp = "financialPreferencez"),
  prior(normal(-0.00757281760064025, 0.35387521728573), class = "b", coef = "leadership_Sum", resp = "financialPreferencez"),
  prior(normal(-0.0576019017149339, 0.297779255969678), class = "b", coef = "Age", resp = "financialPreferencez"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "financialPreferencez"),
  prior(normal(-0.284327282526985, 0.443947922626661), class = "b", coef = "dominance_Sum:Gender2", resp = "financialPreferencez"),
  prior(normal(0.473660365039343, 0.456992923507399), class = "b", coef = "Gender2:prestige_Sum", resp = "financialPreferencez"),
  prior(normal(-0.107204442773651, 0.44522509325961), class = "b", coef = "Gender2:leadership_Sum", resp = "financialPreferencez"),
  prior(normal(0, 1), class = "b", coef = "Gender2:PNI_Sum", resp = "financialPreferencez"),
  prior(normal(1.81050537949088, 0.57441374280994), class = "b", coef = "dominance_Sum", resp = "socialPreferencez"),
  prior(normal(-0.551234838177188, 0.710533759227118), class = "b", coef = "Gender2", resp = "socialPreferencez"),
  prior(normal(0.576829094824566, 0.570410077343782), class = "b", coef = "prestige_Sum", resp = "socialPreferencez"),
  prior(normal(0.143240540516643, 0.557318279915116), class = "b", coef = "leadership_Sum", resp = "socialPreferencez"),
  prior(normal(-0.121074350979943, 0.50181515854942), class = "b", coef = "Age", resp = "socialPreferencez"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "socialPreferencez"),
  prior(normal(-0.869028894989823, 0.675225455650546), class = "b", coef = "dominance_Sum:Gender2", resp = "socialPreferencez"),
  prior(normal(-0.0294477702432216, 0.680282522915085), class = "b", coef = "Gender2:prestige_Sum", resp = "socialPreferencez"),
  prior(normal(0.064613231262165, 0.682913631130467), class = "b", coef = "Gender2:leadership_Sum", resp = "socialPreferencez"),
  prior(normal(0, 1), class = "b", coef = "Gender2:PNI_Sum", resp = "socialPreferencez"),
  prior(normal(1.08508992869999, 0.34337702887508), class = "b", coef = "dominance_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(0.0278355144042715, 0.422659766789246), class = "b", coef = "Gender2", resp = "healthAndSafetyPreferencez"),
  prior(normal(0.0973594822885823, 0.326608098927675), class = "b", coef = "prestige_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(0.00875269191510724, 0.318469718048124), class = "b", coef = "leadership_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(0.29128133142697, 0.268017521856451), class = "b", coef = "Age", resp = "healthAndSafetyPreferencez"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(-0.0999481706002996, 0.401775133599389), class = "b", coef = "dominance_Sum:Gender2", resp = "healthAndSafetyPreferencez"),
  prior(normal(-0.10403479994876, 0.406001141965322), class = "b", coef = "Gender2:prestige_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(-0.168768978826444, 0.40639358037353), class = "b", coef = "Gender2:leadership_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(0, 1), class = "b", coef = "Gender2:PNI_Sum", resp = "healthAndSafetyPreferencez"),
  prior(normal(1.21601692610744, 0.273519339038594), class = "b", coef = "dominance_Sum", resp = "recreationalPreferencez"),
  prior(normal(-1.14450749448218, 0.349856438123967), class = "b", coef = "Gender2", resp = "recreationalPreferencez"),
  prior(normal(-0.520251439573338, 0.268087046058269), class = "b", coef = "prestige_Sum", resp = "recreationalPreferencez"),
  prior(normal(0.0170805281214146, 0.261190196429715), class = "b", coef = "leadership_Sum", resp = "recreationalPreferencez"),
  prior(normal(0.45864492029498, 0.210442768414193), class = "b", coef = "Age", resp = "recreationalPreferencez"),
  prior(normal(0, 1), class = "b", coef = "PNI_Sum", resp = "recreationalPreferencez"),
  prior(normal(-0.0844462213487249, 0.337748960519003), class = "b", coef = "dominance_Sum:Gender2", resp = "recreationalPreferencez"),
  prior(normal(0.521288662601058, 0.33914921337773), class = "b", coef = "Gender2:prestige_Sum", resp = "recreationalPreferencez"),
  prior(normal(-0.246293286419541, 0.346071807160438), class = "b", coef = "Gender2:leadership_Sum", resp = "recreationalPreferencez"),
  prior(normal(0, 1), class = "b", coef = "Gender2:PNI_Sum", resp = "recreationalPreferencez")
)

Demographics Section

Demographic table

Experiment_2_demographics$Gender <- as.factor(Experiment_2_demographics$Gender)
demo_table_j <- d2
label(demo_table_j$Ethnic_Origin) <- "Ethnic Origin"
table1(~ Gender + Age + Education + Ethnicity + Ethnic_Origin, data = demo_table_j)
Overall
(N=289)
Gender
Female 124 (42.9%)
Male 155 (53.6%)
Gender Non-Binary 8 (2.8%)
Prefer not to respond 2 (0.7%)
Age
Mean (SD) 29.3 (9.84)
Median [Min, Max] 26.0 [18.0, 78.0]
Education
Primary School 5 (1.7%)
GCSEs or Equivalent 18 (6.2%)
A-Levels or Equivalent 65 (22.5%)
University Undergraduate Program 130 (45.0%)
University Post-Graduate Program 62 (21.5%)
Doctoral Degree 4 (1.4%)
Prefer not to respond 5 (1.7%)
Ethnicity
White 222 (76.8%)
Mixed or Multi-ethnic 7 (2.4%)
Asian or Asian Scottish or Asian British 5 (1.7%)
African 51 (17.6%)
Other ethnicity 3 (1.0%)
Prefer not to respond 1 (0.3%)
Ethnic Origin
English 16 (5.5%)
European 200 (69.2%)
Latin American 6 (2.1%)
Asian 7 (2.4%)
African 50 (17.3%)
Other 10 (3.5%)

Gender

Experiment_2_demographics$Gender <- as.factor(Experiment_2_demographics$Gender)
ggplot(d2, aes(x = Gender, fill = Gender)) +
  geom_histogram(stat = "count") +
  labs(x = "Gender2") +
  scale_x_discrete(labels = c("Female", "Male", "Gender \nNon-Binary", "Prefer not \nto respond"), guide = "prism_offset") +
  scale_y_continuous(breaks = seq(0, 160, 10), guide = "prism_offset") +
  theme(legend.position = "none")
## Warning: Ignoring unknown parameters: binwidth, bins, pad

Age

Experiment_2_demographics_Gender$Gender <- as.factor(Experiment_2_demographics_Gender$Gender)
d2 <- Experiment_2_demographics_Gender %>%
  mutate_at(vars(locfunc(Experiment_2_demographics_Gender, "Gender")), ~ as.factor(recode(., "1" = "Female", "2" = "Male")))
age_plot <- ggplot(d2, aes(x = Age, fill = Gender)) +
  geom_bar(data = subset(d2, Gender == "Female")) +
  geom_bar(data = subset(d2, Gender == "Male"), aes(y = ..count.. * (-1))) +
  scale_y_continuous(breaks = seq(-30, 30, 1), labels = abs(seq(-30, 30, 1))) +
  scale_x_continuous(breaks = seq(20, 80, 5)) +
  ylab("Number of Participants") +
  xlab("Age of Participants (In years)") +
  geom_hline(yintercept = 0) +
  coord_flip()
ggplotly(age_plot)
## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

Ethnicity

Experiment_2_demographics$Ethnicity <- as.factor(Experiment_2_demographics$Ethnicity)
ggplot(Experiment_2_demographics, aes(x = Ethnicity, fill = Ethnicity)) +
  geom_histogram(stat = "count") +
  scale_x_discrete(labels = c("White ", "Mixed \nor \nMulti-ethnic ", "Asian \nor \nAsian Scottish \nor \nAsian British", "African", "Caribbean \nor \nBlack", "Arab ", "Other ethnicity", "Prefer not \nto respond"), guide = "prism_offset") +
  scale_y_continuous(breaks = seq(0, 250, 20), guide = "prism_offset") +
  theme(legend.position = "none")
## Warning: Ignoring unknown parameters: binwidth, bins, pad

Ethnic Origin

Experiment_2_demographics$Ethnic_Origin <- as.factor(Experiment_2_demographics$Ethnic_Origin)
ggplot(Experiment_2_demographics, aes(x = Ethnic_Origin, fill = Ethnic_Origin)) +
  geom_histogram(stat = "count") +
  scale_x_discrete(labels = c("Scottish", "English", "European", "Latin \nAmerican", "Asian", "Arab", "African", "Other", "Prefer not \nto respond"), guide = "prism_offset") +
  scale_y_continuous(breaks = seq(0, 250, 20), guide = "prism_offset") +
  theme(legend.position = "none")
## Warning: Ignoring unknown parameters: binwidth, bins, pad

Education

Experiment_2_demographics$Education <- as.factor(Experiment_2_demographics$Education)
ggplot(Experiment_2_demographics, aes(x = Education, fill = Education)) +
  geom_histogram(stat = "count") +
  scale_x_discrete(labels = c("Primary School ", "GCSEs \nor \nEquivalent", "A-Levels \nor \nEquivalent", "University \nUndergraduate \nProgram", "University \nPost-Graduate \nProgram", "Doctoral \nDegree", "Prefer not \nto respond"), guide = "prism_offset") +
  scale_y_continuous(breaks = seq(0, 250, 20), guide = "prism_offset") +
  theme(legend.position = "none")
## Warning: Ignoring unknown parameters: binwidth, bins, pad

Analysis

General Correlation

Experiment_2_demographics_Gender$riskBenefitSum_z <- scale(Experiment_2_demographics_Gender$riskBenefitSum)
correlation_df <- Experiment_2_demographics_Gender[, columns_names$x]

correlation_df$dominance_Sum_z <- Experiment_2_demographics$dominance_Sum_z
correlation_df$prestige_Sum_z <- Experiment_2_demographics$prestige_Sum_z
correlation_df$leadership_Sum_z <- Experiment_2_demographics$leadership_Sum_z
correlation_df$PNI_Sum <- as.numeric(scale(Experiment_2_demographics_Gender$PNI_Sum))

correlation_df <- correlation_df %>% rename(
  "Ethical Preference" = "ethicalPreference_z",
  "Financial Preference" = "financialPreference_z",
  "Health and Safety Preference" = "healthAndSafetyPreference_z",
  "Social Preference" = "socialPreference_z",
  "Recreational Preference" = "recreationalPreference_z",
  "Dominance" = "dominance_Sum",
  "Prestige" = "prestige_Sum",
  "Leadership" = "leadership_Sum",
  "UMS Affiliation" = "UMSAffiliationSum_z",
  "UMS Intimacy" = "UMSIntimacySum_z",
  "UMS Sum" = "UMSSum_z",
  "B-PNI" = "PNI_Sum",
  "Risk Perception" = "riskPerceptionSum_z",
  "Risk Benefit" = "riskBenefitSum_z",
  "Risk Sum" = "riskSum_z",
  "General Expected Benefits" = "generalExpectedBenefits_z",
  "General Risk Preference" = "generalRiskPreference_z"
)
corr_1 <- correlation(correlation_df, bayesian = TRUE, method = "auto")
saveRDS(corr_1, "corr_1.rds")
corr_1 <- readRDS("corr_1.rds")
print(summary(corr_1))
## # Correlation Matrix (auto-method)
## 
## Parameter                    | General Risk Preference | General Expected Benefits | Risk Sum | Risk Benefit | Risk Perception |   B-PNI | UMS Sum | UMS Intimacy | UMS Affiliation | Leadership | Prestige | Dominance | Social Preference | Recreational Preference | Health and Safety Preference | Financial Preference
## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
## Ethical Preference           |                 0.68*** |                   0.67*** |  0.42*** |      0.68*** |        -0.30*** |  0.18** |   -0.05 |        -0.02 |          -0.14* |      -0.10 |     0.02 |   0.33*** |              0.08 |                 0.28*** |                      0.56*** |              0.38***
## Financial Preference         |                 0.60*** |                   0.68*** |  0.37*** |      0.68*** |           -0.09 |    0.05 |    0.02 |         0.06 |           -0.08 |       0.10 |     0.06 |     0.14* |           0.23*** |                 0.27*** |                      0.25*** |                     
## Health and Safety Preference |                 0.71*** |                   0.74*** |  0.44*** |      0.74*** |        -0.24*** |  0.15** |   -0.02 |         0.02 |          -0.12* |   6.14e-03 |    -0.07 |   0.27*** |           0.28*** |                 0.50*** |                              |                     
## Recreational Preference      |                 0.68*** |                   0.70*** |  0.43*** |      0.70*** |        -0.23*** |   0.13* |    0.05 |         0.09 |           -0.07 |      0.12* |    -0.01 |   0.21*** |           0.38*** |                         |                              |                     
## Social Preference            |                 0.43*** |                   0.56*** |  0.22*** |      0.56*** |            0.08 | 0.27*** | 0.28*** |      0.27*** |         0.24*** |    0.32*** |  0.22*** |      0.09 |                   |                         |                              |                     
## Dominance                    |                 0.33*** |                   0.30*** |  0.35*** |      0.30*** |        -0.19*** | 0.47*** |   0.11* |        0.13* |            0.01 |    0.29*** |  0.30*** |           |                   |                         |                              |                     
## Prestige                     |                    0.03 |                      0.06 |   0.17** |         0.06 |            0.05 | 0.45*** | 0.66*** |      0.62*** |         0.55*** |    0.46*** |          |           |                   |                         |                              |                     
## Leadership                   |                    0.08 |                     0.13* |    0.14* |        0.13* |            0.07 | 0.29*** | 0.42*** |      0.40*** |         0.35*** |            |          |           |                   |                         |                              |                     
## UMS Affiliation              |                  -0.12* |                     -0.06 |    -0.09 |        -0.06 |         0.19*** | 0.34*** | 0.74*** |      0.56*** |                 |            |          |           |                   |                         |                              |                     
## UMS Intimacy                 |                    0.09 |                     0.12* |  0.21*** |        0.12* |            0.03 | 0.27*** | 0.97*** |              |                 |            |          |           |                   |                         |                              |                     
## UMS Sum                      |                    0.04 |                      0.08 |   0.14** |         0.08 |            0.07 | 0.31*** |         |              |                 |            |          |           |                   |                         |                              |                     
## B-PNI                        |                  0.17** |                   0.22*** |  0.26*** |      0.22*** |            0.04 |         |         |              |                 |            |          |           |                   |                         |                              |                     
## Risk Perception              |                -0.58*** |                  -0.23*** | -0.39*** |     -0.23*** |                 |         |         |              |                 |            |          |           |                   |                         |                              |                     
## Risk Benefit                 |                 0.92*** |                   1.00*** |  0.56*** |              |                 |         |         |              |                 |            |          |           |                   |                         |                              |                     
## Risk Sum                     |                 0.62*** |                   0.56*** |          |              |                 |         |         |              |                 |            |          |           |                   |                         |                              |                     
## General Expected Benefits    |                 0.92*** |                           |          |              |                 |         |         |              |                 |            |          |           |                   |                         |                              |
ggcorrplot(corr_1, type = "lower", lab = TRUE, insig = "blank", show.diag = TRUE, sig.level = 0.05) + scale_x_discrete(labels = 1:17) + theme(axis.text.x = element_text(angle = 0, hjust = .5))


apa_table(correlation_table_1, landscape = TRUE, row.names = FALSE, placement = "ht", note = "* denotes significance level", caption = "General Correlaiton Matrix | Experiment 2")
(#tab:General Correlation Summary Table)
General Correlaiton Matrix | Experiment 2
Parameter X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17
Ethical Preference 0.68*** 0.67*** 0.42*** 0.68*** -0.30*** 0.18** -0.05 -0.02 -0.14* -0.10 0.02 0.33*** 0.08 0.28*** 0.56*** 0.38*** 1
Financial Preference 0.60*** 0.68*** 0.37*** 0.68*** -0.09 0.05 0.02 0.06 -0.08 0.10 0.06 0.14* 0.23*** 0.27*** 0.25*** 1.00
Health and Safety Preference 0.71*** 0.74*** 0.44*** 0.74*** -0.24*** 0.15** -0.02 0.02 -0.12* 0.01 -0.07 0.27*** 0.28*** 0.50*** 1.00
Recreational Preference 0.68*** 0.70*** 0.43*** 0.70*** -0.23*** 0.13* 0.05 0.09 -0.07 0.12* -0.01 0.21*** 0.38*** 1.00
Social Preference 0.43*** 0.56*** 0.22*** 0.56*** 0.08 0.27*** 0.28*** 0.27*** 0.24*** 0.32*** 0.22*** 0.09 1.00
Dominance 0.33*** 0.30*** 0.35*** 0.30*** -0.19*** 0.47*** 0.11* 0.13* 0.01 0.29*** 0.30*** 1.00
Prestige 0.03 0.06 0.17** 0.06 0.05 0.45*** 0.66*** 0.62*** 0.55*** 0.46*** 1.00
Leadership 0.08 0.13* 0.14* 0.13* 0.07 0.29*** 0.42*** 0.40*** 0.35*** 1.00
UMS Affiliation -0.12* -0.06 -0.09 -0.06 0.19*** 0.34*** 0.74*** 0.56*** 1.00
UMS Intimacy 0.09 0.12* 0.21*** 0.12* 0.03 0.27*** 0.97*** 1.00
UMS Sum 0.04 0.08 0.14** 0.08 0.07 0.31*** 1.00
B-PNI 0.17** 0.22*** 0.26*** 0.22*** 0.04 1.00
Risk Perception -0.58 -0.23*** -0.39*** -0.23*** 1.00
Risk Benefit 0.92*** 1.00*** 0.56*** 1.00
Risk Sum 0.62*** 0.56*** 1.00
General Expected Benefits 0.92*** 1.00
General Risk Preference 1.00

Note. * denotes significance level

 

B-PNI distribution

## B-PNI distribution

B_PNI_1 <- ggplot(Experiment_2_demographics, aes(x = Age, y = PNI_Sum)) +
  geom_point(size = 0.7, alpha = 0.8, position = "jitter") +
  geom_smooth(method = "lm", se = FALSE, size = 2, alpha = 0.8)
ggplotly(B_PNI_1)
## brms SEM attempt
riska <- lm(riskBenefitSum ~ riskSum, data = Experiment_2_demographics)
riskb <- lm(riskBenefitSum ~ riskPerceptionSum, data = Experiment_2_demographics)
Experiment_2_demographics$generalRiskPreference <- (Experiment_2_J_Analysis[, "riskBenefitSum"] * riska$coefficients[2]) + (Experiment_2_J_Analysis[, "riskPerceptionSum"] * riskb$coefficients[2])

Experiment_2_demographics_Gender$generalExpectedBenefits <- (Experiment_2_J_Analysis[, "riskBenefitSum"] * riska$coefficients[2])
Experiment_2_demographics_Gender$generalPerceievedRisk <- (Experiment_2_J_Analysis[, "riskPerceptionSum"] * riskb$coefficients[2])
Experiment_2_demographics_Gender$riskBenefitSum_z <- scale(Experiment_2_demographics_Gender$riskBenefitSum)
Experiment_2_demographics_Gender$generalRiskPreference_z <- scale(Experiment_2_demographics$generalRiskPreference)
Experiment_2_demographics_Gender$generalExpectedBenefits_z <- scale(Experiment_2_demographics$generalExpectedBenefits)
Experiment_2_demographics_Gender$generalPerceievedRisk_z <- scale(Experiment_2_demographics$generalPerceievedRisk)

Experiment_2_demographics$UMSSum_z <- scale(Experiment_2_demographics$UMSSum)
Experiment_2_demographics$UMSAffiliationSum_z <- scale(Experiment_2_demographics$UMSAffiliationSum)
Experiment_2_demographics$UMSIntimacySum_z <- scale(Experiment_2_demographics$UMSIntimacySum)


Experiment_2_demographics_Gender <- Experiment_2_demographics[!(Experiment_2_demographics$Gender == "3" | Experiment_2_demographics$Gender == "6"), ]

m1 Interaction model

m1_interaction <- brm(generalRiskPreference_z ~ dominance_Sum * Gender + prestige_Sum * Gender + leadership_Sum * Gender + PNI_Sum * Gender + Age, data = Experiment_2_demographics_Gender, warmup = 1000, iter = 10000, prior = prior_m1_interaction_gen, save_pars = save_pars(all = T), cores = parallel::detectCores(), backend = "cmdstanr")
saveRDS(m1_interaction, "m1_interaction.rds")
m1_interaction_no_pni <- brm(generalRiskPreference_z ~ dominance_Sum * Gender + prestige_Sum * Gender + leadership_Sum * Gender + Age, data = Experiment_2_demographics_Gender, warmup = 1000, iter = 10000, prior = prior_m1_interaction_gen_no_pni, save_pars = save_pars(all = T), cores = parallel::detectCores(), backend = "cmdstanr")
saveRDS(m1_interaction_no_pni, "m1_interaction_no_pni.rds")
# Model without PNI is favored over model with PNI
m1 <- brm(generalRiskPreference_z ~ dominance_Sum + prestige_Sum + leadership_Sum + PNI_Sum + Gender + Age, data = Experiment_2_demographics_Gender, warmup = 1000, iter = 10000, , cores = parallel::detectCores(), backend = "cmdstanr", prior = prior_m1, save_pars = save_pars(all = T))
saveRDS(m1, "m1.rds")
m1_fixef <- MutateHDI::mutate_each_hdi_no_save(brms::fixef(m1))


kable(m1_fixef[
  sign_match(m1_fixef[, 4]) == sign_match(m1_fixef[, 5]),
  c("Parameter", "Estimate", "CI", "CI Low", "CI High")
], format = "html", booktabs = T, escape = F, longtable = F, digits = 2) %>%
  kable_styling(full_width = T) %>%
  remove_column(1)
Parameter Estimate CI CI Low CI High
Intercept 0.59 0.95 0.22 0.95
Dominance 0.28 0.95 0.16 0.41
Gender 0.23 0.95 0.02 0.45
Age -0.03 0.95 -0.04 -0.01

Model comparison

m1_comparison <- loo(m1, m1_interaction)
m1_comparison

bayes_factor(m1, m1_interaction)

# m1_interaction over m1

m2 Additive Model with DoPL and DOSPERT + PNI

m2 <- brm(mvbind(riskSum_z, riskBenefitSum_z, riskPerceptionSum_z) ~ dominance_Sum + prestige_Sum + leadership_Sum + PNI_Sum + Gender + Age, data = Experiment_2_demographics_Gender, warmup = 1000, iter = 10000, prior = prior_m2, save_pars = save_pars(all = T), cores = parallel::detectCores(), backend = "cmdstanr")
saveRDS(m2, "m2.rds")
m2_fixef <- MutateHDI::mutate_each_hdi_no_save(brms::fixef(m2))


kable(m2_fixef[
  sign_match(m2_fixef[, 4]) == sign_match(m2_fixef[, 5]),
  c("Parameter", "Estimate", "CI", "CI Low", "CI High")
], format = "html", booktabs = T, escape = F, longtable = F, digits = 2) %>%
  kable_styling(full_width = T) %>%
  remove_column(1)
Parameter Estimate CI CI Low CI High
Risk Benefit * Intercept 0.39 0.95 0.06 0.71
Risk * Dominance 0.3 0.95 0.2 0.4
Risk * Gender 0.25 0.95 0.06 0.43
Risk Benefit * Dominance 0.24 0.95 0.14 0.35
Risk Benefit * Age -0.02 0.95 -0.03 -0.01
Risk Perception * Dominance -0.27 0.95 -0.38 -0.16
Risk Perception * Gender -0.29 0.95 -0.48 -0.09

Additive Model with DoPL and DOSPERT + PNI Gender Interaction

m2_interaction_gender <- brm(mvbind(riskSum_z, riskBenefitSum_z, riskPerceptionSum_z) ~ dominance_Sum * Gender + prestige_Sum * Gender + leadership_Sum * Gender + PNI_Sum * Gender + Age, data = Experiment_2_demographics_Gender, warmup = 1000, iter = 10000, prior = prior_m2_interaction_gender, save_pars = save_pars(all = T), cores = parallel::detectCores(), backend = "cmdstanr")
saveRDS(m2_interaction_gender, "m2_interaction_gender.rds")
m2_interaction_gender_fixef <- MutateHDI::mutate_each_hdi_no_save(brms::fixef(m2_interaction_gender))


kable(m2_interaction_gender_fixef[
  sign_match(m2_interaction_gender_fixef[, 4]) == sign_match(m2_interaction_gender_fixef[, 5]),
  c("Parameter", "Estimate", "CI", "CI Low", "CI High")
], format = "html", booktabs = T, escape = F, longtable = F, digits = 2) %>%
  kable_styling(full_width = T) %>%
  remove_column(1)
Parameter Estimate CI CI Low CI High
Risk * Intercept 0.42 0.95 0.05 0.81
Risk Benefit * Intercept 0.88 0.95 0.5 1.26
Risk * Dominance 0.5 0.95 0.35 0.65
Risk * Dominance : Gender -0.26 0.95 -0.45 -0.06
Risk Benefit * Dominance 0.32 0.95 0.16 0.48
Risk Benefit * Age -0.02 0.95 -0.04 -0.01
Risk Perception * Dominance -0.3 0.95 -0.47 -0.14
Risk Perception * Leadership : Gender 0.27 0.95 0.07 0.47

Model Comparison

m2_comparison <- loo(m2, m2_interaction_gender)
m2_comparison
bayes_factor(m2_interaction_gender, m2)

# m2 over m2_interaction_gender

DOSPERT and DoPL and PNI

m3 <- brm(mvbind(ethicalPreference_z, financialPreference_z, healthAndSafetyPreference_z, recreationalPreference_z, socialPreference_z) ~ dominance_Sum + prestige_Sum + leadership_Sum + PNI_Sum + Gender + Age, data = Experiment_2_demographics_Gender, warmup = 1000, iter = 10000, prior = prior_m3, save_pars = save_pars(all = T), cores = parallel::detectCores(), backend = "cmdstanr")
saveRDS(m3, "m3.rds")
m3_fixef <- MutateHDI::mutate_each_hdi_no_save(brms::fixef(m3))


kable(m3_fixef[
  sign_match(m3_fixef[, 4]) == sign_match(m3_fixef[, 5]),
  c("Parameter", "Estimate", "CI", "CI Low", "CI High")
], format = "html", booktabs = T, escape = F, longtable = F, digits = 2) %>%
  kable_styling(full_width = T) %>%
  remove_column(1)
Parameter Estimate CI CI Low CI High
Ethical Preference * Intercept 0.56 0.95 0.19 0.93
Health and Safety Preference * Intercept 0.63 0.95 0.24 1.02
Recreational Preference * Intercept 0.84 0.95 0.46 1.23
Social Preference * Intercept 0.76 0.95 0.4 1.13
Ethical Preference * Dominance 0.4 0.95 0.28 0.53
Ethical Preference * Leadership -0.19 0.95 -0.3 -0.07
Ethical Preference * Gender 0.23 0.95 0.02 0.43
Ethical Preference * Age -0.02 0.95 -0.03 -0.01
Financial Preference * Dominance 0.18 0.95 0.05 0.32
Health and Safety Preference * Dominance 0.39 0.95 0.26 0.52
Health and Safety Preference * Prestige -0.2 0.95 -0.33 -0.07
Recreational Preference * Dominance 0.28 0.95 0.16 0.41
Recreational Preference * Prestige -0.18 0.95 -0.31 -0.05
Recreational Preference * Age -0.03 0.95 -0.04 -0.01
Social Preference * Leadership 0.24 0.95 0.12 0.36
Social Preference * Gender -0.49 0.95 -0.71 -0.27

m3 Interaction Gender DOSPERT, DoPL, and PNI

m3_interaction_gender <- brm(mvbind(ethicalPreference_z, financialPreference_z, healthAndSafetyPreference_z, recreationalPreference_z, socialPreference_z) ~ dominance_Sum * Gender + prestige_Sum * Gender + leadership_Sum * Gender + PNI_Sum * Gender + Age, data = Experiment_2_demographics_Gender, warmup = 1000, iter = 10000, prior = prior_m3_interaction_gender, save_pars = save_pars(all = T), cores = parallel::detectCores(), backend = "cmdstanr")

saveRDS(m3_interaction_gender, "m3_interaction_gender.rds")
m3_interaction_gender_fixef <- MutateHDI::mutate_each_hdi_no_save(brms::fixef(m3_interaction_gender))


kable(m3_interaction_gender_fixef[
  sign_match(m3_interaction_gender_fixef[, 4]) == sign_match(m3_interaction_gender_fixef[, 5]),
  c("Parameter", "Estimate", "CI", "CI Low", "CI High")
], format = "html", booktabs = T, escape = F, longtable = F, digits = 2) %>%
  kable_styling(full_width = T) %>%
  remove_column(1)
Parameter Estimate CI CI Low CI High
Ethical Preference * Intercept 0.59 0.95 0.22 0.97
Financial Preference * Intercept 0.45 0.95 0.04 0.88
Health and Safety Preference * Intercept 0.63 0.95 0.23 1.03
Recreational Preference * Intercept 0.86 0.95 0.47 1.26
Social Preference * Intercept 0.79 0.95 0.4 1.16
Ethical Preference * Dominance 0.43 0.95 0.25 0.61
Ethical Preference * Gender 0.23 0.95 0.02 0.43
Ethical Preference * Age -0.02 0.95 -0.03 -0.01
Financial Preference * Dominance 0.28 0.95 0.07 0.49
Health and Safety Preference * Dominance 0.39 0.95 0.2 0.59
Health and Safety Preference * Prestige -0.36 0.95 -0.55 -0.17
Health and Safety Preference * Prestige : Gender 0.26 0.95 0.02 0.5
Recreational Preference * Dominance 0.35 0.95 0.16 0.54
Recreational Preference * Prestige -0.31 0.95 -0.49 -0.12
Recreational Preference * Age -0.03 0.95 -0.04 -0.01
Recreational Preference * Gender: PNI -0.31 0.95 -0.58 -0.04
Social Preference * Gender -0.49 0.95 -0.71 -0.27
Social Preference * Leadership 0.22 0.95 0.04 0.41

Model Comparison

m3_comparison <- loo(m3, m3_interaction_gender)
m3_comparison
bayes_factor(m3_interaction_gender, m3)

# m3 over m3_interaction_gender

Mediation

Mediation Model attempt

mediation_model.1 <- bf(riskSum_z ~ riskBenefitSum_z + riskPerceptionSum_z + PNI_Sum)
mediation_model.2 <- bf(riskBenefitSum_z ~ riskSum_z + riskPerceptionSum_z + PNI_Sum)
mediation_model_1 <- brm(mediation_model.1 + mediation_model.2 + set_rescor(FALSE), warmup = 1000, iter = 10000, data = mediation_dataset, backend = "cmdstanr", cores = parallel::detectCores(), save_pars = save_pars(all = TRUE))

mediation_model.3 <- bf(riskSum_z ~ riskBenefitSum_z + riskPerceptionSum_z + dominance_Sum)
mediation_model.4 <- bf(riskBenefitSum_z ~ riskSum_z + riskPerceptionSum_z + dominance_Sum)
mediation_model_2 <- brm(mediation_model.3 + mediation_model.4 + set_rescor(FALSE), warmup = 1000, iter = 10000, data = mediation_dataset, backend = "cmdstanr", cores = parallel::detectCores(), save_pars = save_pars(all = TRUE))

mediation_model.5 <- bf(riskSum_z ~ riskBenefitSum_z + riskPerceptionSum_z + prestige_Sum)
mediation_model.6 <- bf(riskBenefitSum_z ~ riskSum_z + riskPerceptionSum_z + prestige_Sum)
mediation_model_3 <- brm(mediation_model.5 + mediation_model.6 + set_rescor(FALSE), warmup = 1000, iter = 10000, data = mediation_dataset, backend = "cmdstanr", cores = parallel::detectCores(), save_pars = save_pars(all = TRUE))

mediation_model.7 <- bf(riskSum_z ~ riskBenefitSum_z + riskPerceptionSum_z + leadership_Sum)
mediation_model.8 <- bf(riskBenefitSum_z ~ riskSum_z + riskPerceptionSum_z + leadership_Sum)
mediation_model_4 <- brm(mediation_model.7 + mediation_model.8 + set_rescor(FALSE), warmup = 1000, iter = 10000, data = mediation_dataset, backend = "cmdstanr", cores = parallel::detectCores(), save_pars = save_pars(all = TRUE))
mediation_loo <- brms::loo(mediation_model_1, mediation_model_2, mediation_model_3, mediation_model_4)
mediation_comparison <- bayesfactor_models(mediation_model_1, mediation_model_2, mediation_model_3, mediation_model_4, denominator = mediation_model_4)
saveRDS(mediation_loo, "mediation_loo.rds")
saveRDS(mediation_comparison, "mediation_comparison.rds")
# I think this indicates that model 4 with dominance is the strongest predictor
mediation_loo
## Output of model 'mediation_model_1':
## 
## Computed from 36000 by 279 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -668.6 23.8
## p_loo        11.7  1.6
## looic      1337.3 47.5
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'mediation_model_2':
## 
## Computed from 36000 by 279 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -667.7 23.3
## p_loo        11.5  1.5
## looic      1335.4 46.6
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'mediation_model_3':
## 
## Computed from 36000 by 279 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -670.6 23.5
## p_loo        11.7  1.5
## looic      1341.1 47.0
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'mediation_model_4':
## 
## Computed from 36000 by 279 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -673.7 24.1
## p_loo        11.7  1.5
## looic      1347.4 48.2
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##                   elpd_diff se_diff
## mediation_model_2  0.0       0.0   
## mediation_model_1 -0.9       3.1   
## mediation_model_3 -2.8       4.4   
## mediation_model_4 -6.0       3.0
mediation_comparison
## Bayes Factors for Model Comparison
## 
##     Model     BF
## [1]       143.64
## [2]       353.65
## [3]        21.68
## 
## * Against Denominator: [4]
## *   Bayes Factor Type: marginal likelihoods (bridgesampling)

m_model_1 <- bf(riskSum_z ~ riskBenefitSum_z + riskPerceptionSum_z + riskPerceptionSum_z + PNI_Sum)
m_model_2 <- bf(riskBenefitSum_z ~ PNI_Sum)

Mediation_comparison_1 <- brm(m_model_1 + m_model_2 + set_rescor(FALSE), warmup = 1000, iter = 10000, data = mediation_dataset, backend = "cmdstanr", save_pars = save_pars(all = TRUE))

summary(Mediation_comparison_1)
##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: riskSum_z ~ riskBenefitSum_z + riskPerceptionSum_z + riskPerceptionSum_z + PNI_Sum 
##          riskBenefitSum_z ~ PNI_Sum 
##    Data: mediation_dataset (Number of observations: 279) 
##   Draws: 4 chains, each with iter = 10000; warmup = 1000; thin = 1;
##          total post-warmup draws = 36000
## 
## Population-Level Effects: 
##                              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## riskSumz_Intercept              -0.01      0.05    -0.10     0.08 1.00    64647
## riskBenefitSumz_Intercept       -0.00      0.06    -0.12     0.11 1.00    63908
## riskSumz_riskBenefitSum_z        0.46      0.05     0.36     0.55 1.00    56122
## riskSumz_riskPerceptionSum_z    -0.30      0.05    -0.39    -0.21 1.00    59026
## riskSumz_PNI_Sum                 0.17      0.05     0.08     0.26 1.00    63480
## riskBenefitSumz_PNI_Sum          0.22      0.06     0.11     0.34 1.00    64607
##                              Tail_ESS
## riskSumz_Intercept              25739
## riskBenefitSumz_Intercept       26697
## riskSumz_riskBenefitSum_z       30146
## riskSumz_riskPerceptionSum_z    28832
## riskSumz_PNI_Sum                29119
## riskBenefitSumz_PNI_Sum         25675
## 
## Family Specific Parameters: 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sigma_riskSumz            0.76      0.03     0.70     0.83 1.00    64866
## sigma_riskBenefitSumz     0.98      0.04     0.90     1.07 1.00    63886
##                       Tail_ESS
## sigma_riskSumz           27469
## sigma_riskBenefitSumz    27507
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Mediation_comparison_1_fixef <- MutateHDI::mutate_each_hdi_no_save(brms::fixef(Mediation_comparison_1))


kable(Mediation_comparison_1_fixef[
  sign_match(Mediation_comparison_1_fixef[, 4]) == sign_match(Mediation_comparison_1_fixef[, 5]),
  c("Parameter", "Estimate", "CI", "CI Low", "CI High")
], format = "html", booktabs = T, escape = F, longtable = F, digits = 2) %>%
  kable_styling(full_width = T) %>%
  remove_column(1)
Parameter Estimate CI CI Low CI High
Risk * Risk Benefit 0.46 0.95 0.36 0.55
Risk * Risk Perception -0.3 0.95 -0.39 -0.21
Risk * PNI 0.17 0.95 0.08 0.26
Risk Benefit * PNI 0.22 0.95 0.11 0.34 

m_model_3 <- bf(riskSum_z ~ riskBenefitSum_z + riskPerceptionSum_z + PNI_Sum + dominance_Sum)
m_model_4 <- bf(riskBenefitSum_z ~ PNI_Sum + dominance_Sum)

Mediation_comparison_2 <- brm(m_model_3 + m_model_4 + set_rescor(FALSE), warmup = 1000, iter = 10000, data = mediation_dataset, backend = "cmdstanr", save_pars = save_pars(all = TRUE))
saveRDS(Mediation_comparison_2, "Mediation_comparison_2.rds")
summary(Mediation_comparison_2)
##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: riskSum_z ~ riskBenefitSum_z + riskPerceptionSum_z + PNI_Sum + dominance_Sum 
##          riskBenefitSum_z ~ PNI_Sum + dominance_Sum 
##    Data: mediation_dataset (Number of observations: 279) 
##   Draws: 4 chains, each with iter = 10000; warmup = 1000; thin = 1;
##          total post-warmup draws = 36000
## 
## Population-Level Effects: 
##                               Estimate Est.Error l-95% CI u-95% CI Rhat
## riskSumz_Intercept               -0.01      0.05    -0.10     0.08 1.00
## riskBenefitSumz_Intercept        -0.00      0.06    -0.11     0.11 1.00
## riskSumz_riskBenefitSum_z         0.44      0.05     0.34     0.54 1.00
## riskSumz_riskPerceptionSum_z     -0.28      0.05    -0.37    -0.19 1.00
## riskSumz_PNI_Sum                  0.12      0.05     0.01     0.22 1.00
## riskSumz_dominance_Sum            0.11      0.05     0.01     0.22 1.00
## riskBenefitSumz_PNI_Sum           0.10      0.07    -0.03     0.23 1.00
## riskBenefitSumz_dominance_Sum     0.26      0.06     0.13     0.38 1.00
##                               Bulk_ESS Tail_ESS
## riskSumz_Intercept               71590    26117
## riskBenefitSumz_Intercept        72103    27881
## riskSumz_riskBenefitSum_z        61552    27485
## riskSumz_riskPerceptionSum_z     63021    29332
## riskSumz_PNI_Sum                 52292    30237
## riskSumz_dominance_Sum           53907    30708
## riskBenefitSumz_PNI_Sum          55044    30846
## riskBenefitSumz_dominance_Sum    54338    29457
## 
## Family Specific Parameters: 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sigma_riskSumz            0.76      0.03     0.70     0.83 1.00    68389
## sigma_riskBenefitSumz     0.96      0.04     0.88     1.04 1.00    64282
##                       Tail_ESS
## sigma_riskSumz           26952
## sigma_riskBenefitSumz    26142
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Mediation_comparison_2_fixef <- MutateHDI::mutate_each_hdi_no_save(brms::fixef(Mediation_comparison_2))


kable(Mediation_comparison_2_fixef[
  sign_match(Mediation_comparison_2_fixef[, 4]) == sign_match(Mediation_comparison_2_fixef[, 5]),
  c("Parameter", "Estimate", "CI", "CI Low", "CI High")
], format = "html", booktabs = T, escape = F, longtable = F, digits = 2) %>%
  kable_styling(full_width = T) %>%
  remove_column(1)
Parameter Estimate CI CI Low CI High
Risk * Risk Benefit 0.44 0.95 0.34 0.54
Risk * Risk Perception -0.28 0.95 -0.37 -0.19
Risk * PNI 0.12 0.95 0.01 0.22
Risk * Dominance 0.11 0.95 0.01 0.22
Risk Benefit * Dominance 0.26 0.95 0.13 0.38 

m_model_5 <- bf(riskSum_z ~ riskBenefitSum_z + riskPerceptionSum_z + PNI_Sum + dominance_Sum + leadership_Sum + prestige_Sum)
m_model_6 <- bf(riskBenefitSum_z ~ PNI_Sum + dominance_Sum + leadership_Sum + prestige_Sum)

Mediation_comparison_3 <- brm(m_model_5 + m_model_6 + set_rescor(FALSE), warmup = 1000, iter = 10000, data = mediation_dataset, backend = "cmdstanr", save_pars = save_pars(all = TRUE), cores = parallel::detectCores())
model_blavaan_test <- "
riskSum_z ~ riskBenefitSum_z + riskPerceptionSum_z + PNI_Sum + dominance_Sum
riskBenefitSum_z ~ PNI_Sum + dominance_Sum
"

fit2 <- sem(model_blavaan_test, data = mediation_dataset)

graph_sem(fit2)