---
title: "Cultural Normativity Index: Study 2"
format:
pdf:
prefer-html: true
mainfont: Times
mainfontoptions:
- Numbers=OldStyle
- Ligatures=TeX
monofontoptions: Scale=0.9
fontsize: 11.2pt
code-block-font-size: 9pt
toc: true
number-sections: true
number-depth: 3
section-divs: true
toc-depth: 3
secnumdepth: 3
lof: true
lot: true
execute:
warning: false
message: false
timeout: 300
code-block-wrap: true
code-overflow: wrap
include-in-header:
text: |
\usepackage[T1]{fontenc}
\usepackage{enumitem}
documentclass: article
geometry: margin=1in
papersize: letter
keep-tex: true
crossref:
fig-prefix: "Figure"
lof-title: "List of Figures"
extra_dependencies: ["pdflscape", "booktabs", "longtable"]
editor: source
---
\setlistdepth{6}
```{r eval=FALSE, echo=FALSE}
#Notes: For code chunks that have `DONT-RUN` in their labels involve computations that take longer, require parallel processing or higher GPUs. These will be typically followed by code for loading objects that have been saved after running these chunks.
```
# Data set up
The data, originally in the .sav format, is identified by labels for each set of items. Scale scores for individual and community (here, country) level are computed for each measure.
```{r Libraries}
options(timeout = max(300, getOption("timeout")))
#if (!require("pacman")) install.packages("pacman")
#pacman::p_install() #for installing libraries
pacman::p_load(tidyverse, dplyr, rio, expss, scales, glue,
rlang, kableExtra, stringr, tibble, psych, here, tools,
haven, parallel, tictoc, papaja, #wrangling and setup
patchwork,gridExtra, ggplot2, grid, ggtext, ggcorrplot, #plots
lme4, lmerTest, broom.mixed, tidyr, #modeling - mlm
emmeans, marginaleffects, modelsummary,
pbkrtest, car #model estimates
)
options(modelsummary_get = "easystats")
```
```{r data2-wrangle}
#| cache: true
#---Data import
all_data2 <- import(here::here("data", "SWV2012.sav"), setclass = "tibble")
#get all factor levels
all_data2_factor_cols <- read_sav(here::here("data", "SWV2012.sav"))%>%
mutate(across(c("Nation33", "Language", "D_SX", "D_FH", "D_ME", "D_FE"), haven::as_factor)) %>%
select("Nation33", "Language", "D_SX", "D_FH", "D_ME", "D_FE") %>%
tibble()
#attr(all_data2_spss$Nation33, "labels")
#---Globe items rename
#Study1 pattern of Globe names, only pattern is similar the names here don't exactly map on the same variables as in Study 1
s1_GlobeNames <- c("g1igc1", "g1un1", "g1pf1", "g1ic1", "g1igc2", "g1un2", "g1igc3", "g1ge1", "g1po1", "g1un3", "g1igc4", "g2fu", "g2hu", "g2ass1", "g2fu1", "g2hu1", "g2po1", "g2ass2", "g2pf", "g2fu2", "g2hu2", "g2ge", "g2un3", "g2hu3", "g2po2", "g2fu3", "g2hu4", "g2po3", "g2ic2", "g2ge1", "g2pf2", "g2fu4", "g2pf3", "g3ic3", "g3ic4", "g3ic5", "g3ic6", "g3ass", "g3ass2", "g3ge", "g2po5", "g2po6", "g1un4")
#Study2 Globe names, these are replaced with the names above to replicate the coding style from Study1
s2_GlobeNames <- paste0("NO", 1:43)
#tibble(s2_GlobeNames, s1_GlobeNames) %>% View()
all_data2 <- all_data2%>%
rename(!!!setNames(s2_GlobeNames, s1_GlobeNames))
#---Classify items for domains and compute scale scores for each domain
#This ensures that the items are correctly classified within their respective domains
all_vars_d <- tibble(all_vars = names(all_data2)) %>%
mutate(
domain = case_when(
str_starts(all_vars, "ISM") ~ "Isms",
str_detect(all_vars, "^P\\d{1,2}") ~ "Personality6",
str_starts(all_vars, "TL") ~ "Tight_Lose",
str_starts(all_vars, "MT") ~ "Material",
str_starts(all_vars, "AM") ~ "Amoralism",
str_starts(all_vars, "IV") ~ "Fanaticism",
str_detect(all_vars, "^VAL\\d{1,2}$") ~ "Values", #Schwartz
str_detect(all_vars, "^M\\d{1,2}$") ~ "Machiavellianism",
str_detect(all_vars, "^EN") ~ "Ethnonational",
str_starts(all_vars, "IC") ~ "Indv_Collectivism",
str_detect(all_vars, "^MF\\d{2}$") ~ "Moral_Foundations",
str_starts(all_vars, "FV") ~ "Family_Values",
str_detect(all_vars, "^PA\\d{1,2}$") ~ "Prone_Aggr",
str_starts(all_vars, "g1|g2|g3") ~ "Globe",
str_detect(all_vars, "^NO(?:4[4-9])$") ~ "Globe_extra",
str_starts(all_vars, "Ax") ~ "Social_Axioms",
#Duke Religiosity Index
str_detect(all_vars, "^DRI\\d{2}$") ~ "DRI",
TRUE ~ all_vars
))
#get all item numbers for each domain
item_names <- \(domain_name){
all_item_names <- all_vars_d %>%
filter(domain == {{domain_name}}) %>%
pull(all_vars)
return(all_item_names)
}
compute_trait_score <- function(data, pos_items, rev_items = NULL) {
if (is.null(rev_items)) {
# If no reverse items are specified, just sum the positive items
return(rowSums(data %>% select(all_of(pos_items)), na.rm = TRUE))
} else {
# If both positive and reverse items are specified
rowSums(
data %>%
select(all_of(c(pos_items, rev_items))) %>%
#here the items are assumed to be on the scale of 1 to 5;
#this function is only applied after ensuring that all items are
#are measured on or a rescaled to 1to5
mutate(across(all_of(rev_items), ~ 6 - .x)),
na.rm = TRUE
)
}
}
```
```{r data2}
# convert any input into a character string
convert_to_char <- \(...) {
expr <- enquos(...)
char_vector <- map_chr(expr, quo_name)
return(char_vector)
}
#---Compute scale scores
data2 <- bind_cols(
all_data2_factor_cols,
all_data2 %>%
select(D_AGE,
list_c(map(c("Isms","Globe","Globe_extra", "Tight_Lose", "Material","Social_Axioms", "Amoralism", "Fanaticism", "Ethnonational", "Indv_Collectivism", "Prone_Aggr", "Values", "Family_Values","Personality6", "DRI", "Moral_Foundations", "Machiavellianism"), item_names)))
) %>%
rename("country" = "Nation33",
"language_admin" = "Language",
"sex" = "D_SX",
"family_home" = "D_FH",
"mother_edu" = "D_ME",
"father_edu" = "D_FE",
"age" = "D_AGE") %>%
#Retain countries with n >100
group_by(country) %>%
filter(n() >= 100) %>%
ungroup() %>%
# Drop rows with more than 10% missing data
#filter(rowSums(is.na(.)) / ncol(.) <= 0.1)
#Rescale
mutate(across(c(
all_of(item_names("Values")),
all_of(item_names("Isms")),
all_of(item_names("Fanaticism")),
all_of(item_names("Machiavellianism")),
all_of(item_names("Moral_Foundations")),
all_of(item_names("Globe")),
all_of(item_names("Globe_extra")),
all_of(item_names("Prone_Aggr")),
all_of(item_names("Tight_Lose"))),
# rescale all items from 1 to 5 for easy comparison
~scales::rescale(., to = c(1, 5)))) %>%
#-----[CNI ~ Continuous moderators] -------#
# scP; scM = scaled Personality; Mindset scores
mutate(
#compute_trait_score(data, pos_items, rev_items)
#Personality: pty6
scP_Consc.indv = compute_trait_score(.,
paste0("P", c(1,13,25)),
paste0("P", c(7,19,31))),
scP_Hon.indv = compute_trait_score(.,
paste0("P", c(11,17,29)),
paste0("P", c(5,23,35))),
scP_Agree.indv = compute_trait_score(.,
paste0("P", c(8,32)),
paste0("P", c(2,14,20,26))),
scP_Res.indv = compute_trait_score(.,
paste0("P", c(12,36)),
paste0("P", c(6,18,24,30))),
scP_Extra.indv = compute_trait_score(.,
paste0("P", c(3,15,27)),
paste0("P", c(9,21,33))),
scP_Vir.indv = compute_trait_score(.,
paste0("P", c(10,16,22,34)),
paste0("P", c(4,28))),
#Personality: pty2 - scoring here is based on GS' method
#Social self-regulation
scP_Big2_SocialSelfReg.indv = (compute_trait_score(.,
paste0("P", c(1,8,17,25,38)),
paste0("P", c(5,7,14,19,23))))/10,
#Dynamism
scP_Big2_Dynamism.indv = (compute_trait_score(.,
paste0("P", c(10,15,34,36,37,39)),
paste0("P", c(9,28,30,40))))/10,
#Moral Foundations
scM_MFQ_Harm.indv = compute_trait_score(.,
c("MF01","MF07","MF12","MF18")),
scM_MFQ_Fairness.indv = compute_trait_score(.,
c("MF02","MF08","MF13","MF19")),
scM_MFQ_Loyalty.indv = compute_trait_score(.,
c("MF03","MF09","MF14","MF20")),
scM_MFQ_Authority.indv = compute_trait_score(.,
c("MF04","MF10","MF15","MF21")),
scM_MFQ_Purity.indv = compute_trait_score(.,
c("MF05","MF11","MF16","MF22")),
#Mindset: Social Axioms
# Social Cynicism
scM_SocAx_Cynicism.indv = compute_trait_score(.,
c("Ax5", "Ax12", "Ax16", "Ax23")),
# Reward for Application
scM_SocAx_RewardForApp.indv = compute_trait_score(.,
c("Ax4", "Ax9", "Ax13", "Ax21", "Ax28")),
# Social Complexity
scM_SocAx_Complexity.indv = compute_trait_score(.,
c("Ax6", "Ax8", "Ax10", "Ax15", "Ax19", "Ax22", "Ax25",
"Ax27", "Ax30")),
# Fate Control
scM_SocAx_FateControl.indv = compute_trait_score(.,
c("Ax1", "Ax7", "Ax11", "Ax14", "Ax24")),
# Religiosity (importance of religion)
scM_SocAx_Religiosity.indv = compute_trait_score(.,
c("Ax2", "Ax3", "Ax17", "Ax18", "Ax20", "Ax26", "Ax29"
)),
#Mindset: Individualism Collectivism
scM_IndvColl_CollecHori.indv = compute_trait_score(.,
c("IC14", "IC16", "IC4", "IC6")),
scM_IndvColl_CollecVert.indv = compute_trait_score(.,
c("IC10", "IC13", "IC2", "IC7")),
scM_IndvColl_IndivHori.indv = compute_trait_score(.,
c("IC1", "IC11", "IC5", "IC8")),
scM_IndvColl_IndivVert.indv = compute_trait_score(.,
c("IC12", "IC15", "IC3", "IC9")),
#Mindset: Family Values
scM_FamValCohesion.indv = compute_trait_score(., c("FV02","FV04", "FV06", "FV08")),
scM_FamValHierarchy.indv = compute_trait_score(., c("FV01","FV03", "FV05", "FV07" )),
#Mindset: Schwartz Values
#Self-Transcendence: Universalism + Benevolence
scM_Value_SfTran.indv = compute_trait_score(.,c("VAL5", "VAL8")),
#Conservation: Conformity + Tradition + Security
scM_Value_Consrv.indv = compute_trait_score(.,c("VAL4", "VAL1", "VAL7")),
#Self-enhancement: Achievement + Power
scM_Value_SelfEnh.indv = compute_trait_score(.,c("VAL3", "VAL10")),
#Openness to Change: Stimulation + Self-direction
scM_Value_OpnChg.indv = compute_trait_score(.,c("VAL9", "VAL2")),
#Fanaticism
scM_fntc.indv = compute_trait_score(.,
item_names("Fanaticism")[-c(4,6,17)],
item_names("Fanaticism")[c(4,6,17)]),
#dput(convert_to_char())
#Mindset: Isms
#Traditional Religiousness
scM_Isms_TradRelig.indv = compute_trait_score(.,c("ISM12", "ISM17", "ISM22", "ISM35"), c("ISM14", "ISM16", "ISM30", "ISM8")),
# Unmitigated Self-Interest
scM_Isms_UnmiSelfInt.indv = compute_trait_score(.,c("ISM11", "ISM3", "ISM32", "ISM4", "ISM41", "ISM42", "ISM46"), c("ISM21", "ISM26", "ISM7")),
# Communal Rationalism
scM_Isms_CommRation.indv = compute_trait_score(.,c("ISM13", "ISM15", "ISM18", "ISM24", "ISM38", "ISM44"), c("ISM20", "ISM43")),
# Subjective Spirituality
scM_Isms_SubjSpirit.indv = compute_trait_score(.,c("ISM10", "ISM39", "ISM5", "ISM6"), c("ISM2", "ISM23", "ISM25", "ISM27")),
# Inequality Aversion
scM_Isms_IneqAver.indv = compute_trait_score(.,c("ISM19", "ISM31", "ISM33", "ISM34", "ISM45"), c("ISM29", "ISM40")),
#Mindset: Materialism
scM_Material.indv = compute_trait_score(.,
item_names("Material")),
#Mindset: Machiavellianism
scM_Machiaveln.indv = compute_trait_score(., paste0("M", c(1,2,4)), paste0("M", c(3,5))),
#Mindset: Proneness to Aggression - Honor
scM_Prone_Aggr.indv = compute_trait_score(.,item_names("Prone_Aggr")),
#Mindset: Ethnonationalism
scM_Ethnonational.indv = compute_trait_score(., paste0("EN", c(1,3,4,6)), paste0("EN", c(2,5))),
#Mindset: Cultural Tightness
scM_Tightness.indv = compute_trait_score(., paste0("TL0", c(1,2,3,5,6)), "TL04"),
#Mindset: Amorialism
scM_Amoralism.indv = compute_trait_score(., item_names("Amoralism")[-c(1,4,6,8,10,11,13)],
item_names("Amoralism")[c(1,4,6,8,10,11,13)]),
#Mindset: Duke Religion Index (DRI)
# Organizational activity
scM_DRI_OrgRelgAc.indv = compute_trait_score(., "DRI01"),
# Non-organizational activity
scM_DRI_NonOrgRelgAc.indv = compute_trait_score(., "DRI02"),
# Intrinsic religiosity
scM_DRI_IntrnscRelg.indv = compute_trait_score(., paste0("DRI0", c(3,4,5))),
# Civic Multicultural Nationalism
scM_CivicMCNational.indv = compute_trait_score(., c("EN2", "EN5")),
#Mindset: GLOBE
#InGroup Collectivism
scM_globe_igc.indv = compute_trait_score(.,
str_subset(item_names("Globe"), "^g.*igc.*$")),
# Institutional Collectivism
scM_globe_ic.indv = compute_trait_score(.,
str_subset(item_names("Globe"), "^g.*ic.*$")[-c(3,4)],
str_subset(item_names("Globe"), "^g.*ic.*$")[c(3,4)]),
# Performance Orientation/Achievement
scM_globe_pf.indv = compute_trait_score(.,
str_subset(item_names("Globe"), "^g.*pf.*$")[-4],
str_subset(item_names("Globe"), "^g.*pf.*$")[4]),
# Uncertainty Avoidance
scM_globe_un.indv = compute_trait_score(.,
str_subset(item_names("Globe"), "^g.*un.*$")[-5],
str_subset(item_names("Globe"), "^g.*un.*$")[5]),
# Gender Egalitarian
scM_globe_ge.indv = compute_trait_score(.,
str_subset(item_names("Globe"), "^g.*ge.*$")[-c(1,2,4)],
str_subset(item_names("Globe"), "^g.*ge.*$")[c(1,2,4)]),
# Power Distance
scM_globe_po.indv = compute_trait_score(.,
str_subset(item_names("Globe"), "^g.*po.*$")[-c(5,6)],
str_subset(item_names("Globe"), "^g.*po.*$")[c(5,6)]),
# Future Orientation
scM_globe_fu.indv = compute_trait_score(.,
str_subset(item_names("Globe"), "^g.*fu.*$")[-c(3,4)],
str_subset(item_names("Globe"), "^g.*fu.*$")[c(3,4)]),
# Assertiveness
scM_globe_ass.indv = compute_trait_score(.,
str_subset(item_names("Globe"), "^g.*ass.*$")[-c(2,3)],
str_subset(item_names("Globe"), "^g.*ass.*$")[c(2,3)]),
# Humane Orientation
scM_globe_hu.indv = compute_trait_score(.,
str_subset(item_names("Globe"), "^g.*hu.*$"))) %>%
mutate(pID = paste0("p", 1:nrow(.)), .before = everything())%>%
#compute community[here, country] scale scores
with_groups(country,
mutate,
across(
# Select all pty and cog scaled scores
contains(c( "scP", "scM")),
# Compute mean by country
~mean(., na.rm = TRUE),
# Create new name [scP_Hon.indv -> scP_Hon.comm]
.names = "{str_remove(.col, '[.].*')}.comm"
)
) %>%
#-----[CNI ~ Categorical moderators]------#
#Demographic variables
mutate(
mother_edu = case_when(
str_detect(mother_edu, "don't know") ~ NA_character_,
TRUE ~ as.character(mother_edu)
),
father_edu = case_when(
str_detect(father_edu, "don't know") ~ NA_character_,
TRUE ~ as.character(father_edu)
),
age_group = case_when(
age >= 13 & age <= 17 ~ "13-17",
age >= 18 & age <= 25 ~ "18-25",
age >= 26 & age <= 35 ~ "26-35",
age >= 36 & age <= 45 ~ "36-45",
age >= 46 ~ "46+",
TRUE ~ NA_character_
))%>%
relocate(age_group, .before = age)
```
```{r tbl-SampleStatistics, echo=F, results='asis'}
#| tbl-cap: "Sample statistics by Country"
#| label: tbl-SampleStatistics
data2 %>%
group_by(country) %>% #names()
select(age, sex) %>%
summarize(mean_age = mean(age, na.rm=TRUE),
percent.male = mean(sex == "male", na.rm=TRUE)*100,
N = n()) %>%
arrange(desc(N)) %>%
kable(format = "latex",
col.names = c("Country", "Mean Age", "% Male", "n"),
booktabs = TRUE) %>%
kable_styling()
```
```{r tbl-SampleDemographics, echo = F, results='asis'}
#| label: tbl-SampleDemographics
#| tbl-cap: "Sample Demographics"
data2 %>%
select(sex, age_group, family_home, mother_edu, father_edu, language_admin) %>%
mutate(across(everything(), ~as.factor(.))) %>%
pivot_longer(cols = everything(), names_to = "variable", values_to = "category") %>%
group_by(variable, category) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(variable) %>%
mutate(
p = 100 * n / sum(n),
sample = paste0(n, " (", papaja::printnum(p), "\\%)")) %>%
ungroup() %>%
select(-n, -p) %>%
#arrange(variable, category) %>%
mutate(
category = if_else(is.na(category), "missing",
as.character(category)),
variable = factor(variable, levels = c("sex", "age_group", "family_home", "mother_edu", "father_edu", "language_admin"))
) %>%
arrange(variable) %>%
dplyr::select(-variable) %>%
kable(longtable = T,
format = "latex",
align = c("l", "c", "c"), escape = FALSE,
col.names = c("Variable", "$\\textit{n}$ (\\%)"),
booktabs = TRUE) %>%
kable_styling(latex_options = c("repeat_header", "longtable")) %>%
group_rows("Sex", 1, 3) %>%
group_rows("Age", 4, 8) %>%
group_rows("Family Home", 9, 13) %>%
group_rows("Mother's Education", 14, 23) %>%
group_rows("Father's Education", 24, 33) %>%
group_rows("Language of Administration", 34, 51)
```
```{r}
#| include: false
#knitr::opts_chunk$set(eval = FALSE, echo = FALSE)
```
# Analysis
## RQ 1: How to quantify CNI?
Similar to Study 1, CNI is estimated using a Random slope Multilevel model (varying slopes and fixed intercepts) wherein items nested within person \[pID\] (Weston,2024). This model accounts for the nested nature of our data, where individual items are clustered within each participant. As responses are ipsatized (standardized within participant), the resulting slope estimates can be interpreted as correlation estimates.
For study 2, eight CNI models are computed:
::: {#note-CNI-models .callout-note title="CNI Models" icon="false"}
1. Mindset CNI: consisting of all mindset items
2. Isms CNI
3. Globe CNI: Original Globe items by House et al.
4. GlobeEx CNI: House et al. + extra items
5. Moral Foundations CNI
6. Fanaticism CNI
7. Social Axioms CNI
8. Personality 6 CNI: consisting of Big 6 items
:::
The Mindset-CNI and Personality6-CNI encompass distinct domains of psychological measurement, while the remaining CNI types can be viewed as subcomponents or specialized facets within the broader Mindset-CNI framework. Mindset-CNI therefore is the most comprehensive of all the CNI types explored here.
```{r funs.-model_CNI}
# functions for computing CNIs
#item names for every CNI type
get_item_names <- \(CNI_type, data) {
domains <- switch(CNI_type,
"Mindset" = c("Isms", "Fanaticism", "Social_Axioms", "Values", "Tight_Lose", "Material",
"Amoralism", "Ethnonational", "Indv_Collectivism", "Globe", "Globe_extra",
"Moral_Foundations", "Family_Values","Machiavellianism","Prone_Aggr"),
"pty6" = c("Personality6"),
"Fanaticism" =c("Fanaticism"),
"Isms" =c("Isms"),
"Globe" =c("Globe"),
"Globe_extra" =c("Globe", "Globe_extra"),
"Moral_foundations" =c("Moral_Foundations"),
"Social_Axioms" = c("Social_Axioms"),
stop("CNI type not found")
)
item_names <- data %>%
filter(domain %in% domains) %>%
pull(all_vars)
if (length(item_names) == 0) {
warning("No items found for the specified domains")
}
return(item_names)
}
#ipsatize/within person standardize scores
#scores are ipsatized before CNI's are computed
ipsatize <- \(x){
value <- (x-mean(x, na.rm=T))/(sd(x, na.rm=T))
return(value)
}
#compute z_profiles for pID
compute_profiles <- \(var_names){
#compute mean country profiles
country_profiles <- data2 %>%
group_by(country) %>%
summarise(
across(all_of(var_names), ~mean(.x, na.rm = TRUE), .names = "avg_{.col}"),#variables for cogCNI
n_p = n())%>% # number of participants)
ungroup() %>%
#retain participants who have responded to at least 30 items.- all countries have n_p>30
filter(n_p>=30) %>%
pivot_longer(
cols = starts_with("avg_"),
#item names stored in item
names_to = "item",
values_to = "ctry_response"
) %>%
mutate(item = str_remove(item, "avg_")) %>%
ungroup()
self.ctry_profiles <- data2 %>%
select(all_of(var_names), pID, country) %>%
pivot_longer(names_to = "item",
values_to = "response",
cols = all_of(var_names)) %>%
filter(!is.na(response)) %>%
left_join(country_profiles, by = c("country", "item")) %>%
filter(!is.na(ctry_response)) %>%
# ipsatize responses within profile
with_groups(pID, mutate, across(
# mutate self-rating and mean country profiles for each pID at once
contains("response"),
# with the ipsatize function ,
ipsatize,
# create new name
.names = "z_{.col}"))
overallM_profiles <- self.ctry_profiles %>%
select(item, z_response) %>%
filter(!is.na(z_response)) %>%
group_by(item) %>%
summarise(overall.M = mean(z_response)) %>%
ungroup()
return(lst(self.ctry_profiles, overallM_profiles)) #set_names() to the input object names
}
```
```{r CNI-profiles}
#| cache: true
#get all item names for CNI types
mindsetCNI_items <- get_item_names("Mindset", all_vars_d)
pty6CNI_items <- get_item_names("pty6", all_vars_d)
ismCNI_items <- get_item_names("Isms", all_vars_d)
fantcCNI_items <- get_item_names("Fanaticism", all_vars_d)
globeCNI_items <- get_item_names("Globe", all_vars_d)
globeExCNI_items <- get_item_names("Globe_extra", all_vars_d)
mfCNI_items <- get_item_names("Moral_foundations", all_vars_d)
socAxCNI_items <- get_item_names("Social_Axioms", all_vars_d)
#compute pofiles for each CNI type
mindsetCNI_profiles <- compute_profiles(mindsetCNI_items)
pty6CNI_profiles <- compute_profiles(pty6CNI_items)
ismCNI_profiles <- compute_profiles(ismCNI_items)
fantcCNI_profiles <- compute_profiles(fantcCNI_items)
globeCNI_profiles <- compute_profiles(globeCNI_items)
globeExCNI_profiles <- compute_profiles(globeExCNI_items)
mfCNI_profiles <- compute_profiles(mfCNI_items)
socAxCNI_profiles <- compute_profiles(socAxCNI_items)
all_profiles <- lst(
mindsetCNI_profiles,
pty6CNI_profiles,
ismCNI_profiles,
fantcCNI_profiles,
globeCNI_profiles,
globeExCNI_profiles,
mfCNI_profiles,
socAxCNI_profiles
)
```
```{r DONT-RUN-fit-all-cni-mods, eval=FALSE, echo=TRUE}
#formula for computing CNI
add_CNIformula <- \(data){
lmer(z_response ~ z_ctry_response + overall.M + (-1 + z_ctry_response|country/pID), data = data)
}
mindsetCNI_mod <- mindsetCNI_profiles$self.ctry_profiles %>%
full_join(cogCNI_profiles$overallM_profiles, by = "item") %>% add_CNIformula(.)
pty6CNI_mod <- pty6CNI_profiles$self.ctry_profiles %>%
full_join(pty6CNI_profiles$overallM_profiles, by = "item") %>% add_CNIformula(.)
ismCNI_mod <- ismCNI_profiles$self.ctry_profiles %>%
full_join(ismCNI_profiles$overallM_profiles, by = "item") %>% add_CNIformula(.)
fantcCNI_mod <- fantcCNI_profiles$self.ctry_profiles %>%
full_join(fantcCNI_profiles$overallM_profiles, by = "item") %>% add_CNIformula(.)
globeCNI_mod <- globeCNI_profiles$self.ctry_profiles %>%
full_join(globeCNI_profiles$overallM_profiles, by = "item") %>% add_CNIformula(.)
globeExCNI_mod <- globeExCNI_profiles$self.ctry_profiles %>%
full_join(globeExCNI_profiles$overallM_profiles, by = "item") %>% add_CNIformula(.)
mfCNI_mod <- mfCNI_profiles$self.ctry_profiles %>%
full_join(mfCNI_profiles$overallM_profiles, by = "item") %>% add_CNIformula(.)
socAxCNI_mod <- socAxCNI_profiles$self.ctry_profiles %>%
full_join(socAxCNI_profiles$overallM_profiles, by = "item") %>% add_CNIformula(.)
```
```{r save-CNImods, echo = F, eval = F}
save(mindsetCNI_mod, pty6CNI_mod, ismCNI_mod, fantcCNI_mod, globeCNI_mod, globeExCNI_mod, mfCNI_mod, socAxCNI_mod,
file = here("objects/Study2/all_CNI.mods_study2.Rdata"))
```
```{r load-CNImods, echo = F}
load(here::here("objects/Study2/all_CNI.mods_study2.Rdata"))
```
### Model Complexity
The Eight models vary in the number of items and domains of psychological measures. 1. The Parallel Analysis computed here, identify the optimum number of dimensions that the data for every CNI type can be reduced to. The optimum number of dimensions are retained based on the comparison of eigenvalues of the actual data and those of randomly generated data sets. Dimensions are retained when their eigenvalues exceed those of the random data, ensuring that only factors accounting for more variance than would be expected by chance are retained. 2. Principal Components that explain at least 50% of variance are identified.
```{r fun.-heterogeneity}
# Function to calculate heterogeneity
compute_heterogeneity <- function(data) {
# Compute correlation matrix once - memory efficient
cor_matrix <- cor(data, use = "pairwise.complete.obs")
pca_rotation <- function(rotation) {
# Perform PCA using correlation matrix
pca <- principal(cor_matrix, nfactors = ncol(data), rotate = rotation, covar = FALSE)
# Calculate cumulative variance explained
cumulative_var <- cumsum(pca$values) / sum(pca$values)
# Find number of factors needed to explain at least 50% variance
heterogeneity <- which(cumulative_var >= 0.5)[1]
return(heterogeneity)
}
# Perform parallel analysis using correlation matrix
pa_result <- fa.parallel(cor_matrix, n.obs = nrow(data), fa = "fa",
fm = "minres", show.legend = FALSE, plot = FALSE)
tibble(
n_items = ncol(data),
pa = pa_result$nfact,
pca_50 = pca_rotation("none")
)
}
```
```{r DONT-RUN-compute-heterogeneity, eval = F}
heterogeneity_stats <- map_dfr(
list(
Mindset = select(data2,
get_item_names("Mindset",
all_vars_d)),
Globe = select(data2,
get_item_names("Globe",
all_vars_d)),
GlobeEx = select(data2,
get_item_names(c("Globe_extra"),
all_vars_d)),
Isms = select(data2,
get_item_names("Isms",
all_vars_d)),
Fanaticism = select(data2,
get_item_names("Fanaticism",
all_vars_d)),
Moral_foundations = select(data2,
get_item_names("Moral_foundations",
all_vars_d)),
Social_Axioms = select(data2,
get_item_names("Social_Axioms",
all_vars_d)),
Personality6 = select(data2,
get_item_names("pty6",
all_vars_d))
),
compute_heterogeneity, .id = "cni_type"
)
```
```{r save-compute-heterogeneity, echo = F, eval = F}
save(heterogeneity_stats,
file = here("objects/Study2/heterogeneity_stats_study2.Rdata"))
```
```{r load-compute-heterogeneity, echo = F}
load(here("objects/Study2/heterogeneity_stats_study2.Rdata"))
```
```{r tbl-heterogeneity, results='asis', echo=F}
#| label: tbl-heterogeneity
#| tbl-cap: "CNI Model complexity"
#| cache: true
heterogeneity_stats %>%
mutate(cni_type = str_replace_all(cni_type, "_", " ")) %>%
kable(booktabs = T,
escape = F,
format = "latex",
col.names = c("CNI type", "No. items", "Parallel analysis", "PCA")) %>%
kable_styling() %>%
group_rows("Mindset", 1, 7) %>%
group_rows("Personality", 8, 8)
```
```{r plot-heterogeneity, echo=F}
#| fig-height: 6
#| fig-width: 10
#| out-width: "100%"
#| cap-location: top
#| fig-cap: "Model Complexity Measures Across CNI Types"
#| label: plot-heterogeneity
heterogeneity_stats %>%
pivot_longer(cols = c(pa, pca_50), names_to = "metric", values_to = "value") %>%
mutate(cni_type = str_replace_all(cni_type, "_", " ")) %>%
mutate(cni_type = case_when(
cni_type == "Globe" ~ "GLOBE",
cni_type == "GlobeEx" ~ "GLOBE Extra",
cni_type == "Moral foundations" ~ "Moral Foundations",
TRUE ~ cni_type
)) %>%
mutate(x_label = paste0(cni_type, "\n(", n_items, " items)"),
x_label = reorder(x_label, n_items)) %>%
# Create the plot
ggplot(., aes(x = x_label, y = value, color = metric, group = metric)) +
geom_line() +
geom_point(size = 3) +
scale_color_manual(values = c("pa" = "#C7B5E5", "pca_50" = "#91CE83"),
labels = c("pa" = "Parallel Analysis", "pca_50" = "PCA 50%"),
name = "Metric") +
labs(x = "CNI Type",
y = "Number of Factors") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 11.5),
legend.position = "bottom",
text = element_text(family = "serif", size = 12))
```
Heterogeneity figure indicates that as the item count increases, so does the number of identifiable dimensions, indicating higher model complexity. The Fanaticism CNI, despite having only 17 items, has been included here to maintain consistency with Study 1 and provide a basis for comparison. Notably, the Mindset CNI exhibits a markedly higher complexity, with nearly four times as many dimensions as the other CNI models.
### CNI distributions
```{r fun.-plot_CNI, echo=FALSE}
#Functions
#model1_est = tidy(model1, conf.int = T)
plot_CNI <- function(model_obj) {
# Get the name of the object passed to the function
obj_name <- deparse(substitute(model_obj))
# Extract the title (everything before the underscore)
title <- str_extract(obj_name, "^[^_]+")
# Define color mapping
color_map <- tibble(
mod_name = c("mindsetCNI", "ismCNI", "fantcCNI", "globeCNI", "globeExCNI", "mfCNI", "socAxCNI","pty6CNI"),
color = c("#0074B7","#6082B6", "#70d6ff", "#a0c4ff","#189AB4", "#6495ED", "#5D3FD3", "#ff70a6")
)
# Determine the color based on the title
plot_color <- color_map %>%
filter(str_detect(mod_name, title)) %>%
pull(color)
# Calculate slopes
mod_slopes <- slopes(model_obj, variables = "z_ctry_response", by = "pID")
coefficient_plot <- mod_slopes %>%
ggplot(aes(x = reorder(pID, estimate), y = estimate)) +
geom_segment(aes(xend = pID,
y = conf.low,yend = conf.high),
alpha = .08 , color = plot_color) +
geom_point(size = .22, color = plot_color) +
scale_x_discrete(breaks = NULL) +
#coord_cartesian(ylim = c(-0.2, .92)) +
labs(x = "Participant", y = "CNI coefficient", title = title) +
theme_bw()
hist_plot <- mod_slopes %>%
as.data.frame() %>%
ggplot(aes(x = estimate)) +
geom_histogram(binwidth = .03, color = "white", fill = plot_color) +
labs(x = "CNI coefficient", y = "Frequency", title = title) +
scale_x_continuous(limits = c(-0.3, 0.92))+
scale_y_continuous(limits = c(0, 1000))+
theme_bw()
return(lst(coefficient_plot,hist_plot ))
}
```
```{r DONT-RUN-CNIplots, echo = T, eval = FALSE}
#mindsetCNI_mod, pty6CNI_mod, ismCNI_mod, fantcCNI_mod, globeCNI_mod, globeExCNI_mod, mfCNI_mod, socAxCNI_mod
mindsetCNI_plot <- plot_CNI(mindsetCNI_mod)
ismCNI_plot <- plot_CNI(ismCNI_mod)
globeCNI_plot <- plot_CNI(globeCNI_mod)
globeExCNI_plot <- plot_CNI(globeExCNI_mod)
fantcCNI_plot <- plot_CNI(fantcCNI_mod)
mfCNI_plot <- plot_CNI(mfCNI_mod)
socAxCNI_plot <- plot_CNI(socAxCNI_mod)
pty6CNI_plot <- plot_CNI(pty6CNI_mod)
```
```{r save-CNImods_plots, echo = F, eval = F}
save(mindsetCNI_plot, ismCNI_plot, globeCNI_plot, globeExCNI_plot, fantcCNI_plot, mfCNI_plot, socAxCNI_plot, pty6CNI_plot,
file = here("objects/Study2/all_CNImods_plots_Study2.Rdata"))
```
```{r load-CNImods_plots, echo = F}
#| cache: true
#takes time to load
load(here("objects/Study2/all_CNImods_plots_Study2.Rdata"))
```
```{r fun.-patch_CNIplots}
patch_CNIplots <- \(plot_list){
p1 <- plot_list$coefficient_plot &
theme(plot.title = element_blank(),
text = element_text(family = "serif", size = 12))
p2 <- plot_list$hist_plot &
theme(plot.title = element_blank(),
text = element_text(family = "serif", size = 12))
combined_plot <- p1 + p2
return(combined_plot)
}
```
```{r print-patch_CNIplots-details}
#| cache: true
# Create a list of plot objects and their corresponding titles
CNIplot_data <- lst(
list(plot = mindsetCNI_plot, title = "Mindset CNI"),
list(plot = ismCNI_plot, title = "Isms-Mindset CNI"),
list(plot = globeCNI_plot, title = "Globe-Mindset CNI"),
list(plot = globeExCNI_plot, title = "Globe (Extra items)-Mindset CNI"),
list(plot = fantcCNI_plot, title = "Fanaticism-Mindset CNI"),
list(plot = mfCNI_plot, title = "Moral Foundations-Mindset CNI"),
list(plot = socAxCNI_plot, title = "Social Axioms-Mindset CNI"),
list(plot = pty6CNI_plot, title = "Personality6 CNI")
)
CNIplot_data <- lst(
mindsetCNI_plot, ismCNI_plot, globeCNI_plot, globeExCNI_plot,
fantcCNI_plot, mfCNI_plot, socAxCNI_plot, pty6CNI_plot)
CNIplot_data_details <- tibble(
mod_names = names(CNIplot_data),
caption = c("Mindset CNI", "Isms-Mindset CNI", "Globe-Mindset CNI",
"Globe (Extra items)-Mindset CNI", "Fanaticism-Mindset CNI",
"Moral Foundations-Mindset CNI", "Social Axioms-Mindset CNI",
"Personality6 CNI")
)
```
```{r print-patch_CNIplots}
#| cache: true
#| fig-height: 5
#| fig-width: 10
#| out-width: "100%"
#| fig-cap-location: top
#| fig-cap:
#| !expr CNIplot_data_details %>%
#| pull(caption)
# `walk` to apply patch_CNIplots to each plot in the list and print quietly
walk(names(CNIplot_data), \(mod_names) {
cat("\n\n")
print(
pluck(CNIplot_data, mod_names) %>%
patch_CNIplots())
cat("\n\n")
})
```
### Correlations across CNI types
```{r comparing-cni-corr-fig, results='asis', echo=FALSE}
#| fig-cap: "Correlations between CNI types"
#| cache: true
join_with_data <- \(cni_mod){
mod_name <- deparse(substitute(cni_mod))
cni_type <- gsub("_mod$", "", mod_name)
join_data <- coef(cni_mod)$pID %>%
as.data.frame() %>%
rownames_to_column("pID") %>%
select(-`(Intercept)`) %>%
separate(pID, into = c("pID", "country")) %>%
rename_with(~ cni_type, .cols = z_ctry_response) %>%
full_join(data2, by = c("pID", "country")) %>% select(all_of(cni_type), pID, country, contains("sc_swb")) %>%
as_tibble()
return(join_data)
}
#join_with_data(mindsetCNI_mod)
all_cni_mods <- c("mindsetCNI_mod", "pty6CNI_mod", "ismCNI_mod", "fantcCNI_mod",
"globeCNI_mod", "globeExCNI_mod", "mfCNI_mod", "socAxCNI_mod")
all_combined_data <- map(all_cni_mods, function(mod_name) {
mod_obj <- get(mod_name, envir = .GlobalEnv)
#simulate calling the function with the object name
eval(call("join_with_data", as.name(mod_name)))
}) %>%
set_names(names(all_cni_mods)) %>%
#combine all the tibbles in the list [map_dfr will append the tibbles increasing the no. of rows; the reduce approach retains the original set of rows by using full_join]
reduce(full_join, by = c("pID", "country")) %>%
# Ensure all CNI columns are present (in case of any missing data)
mutate(across(ends_with("CNI"), ~ if_else(is.na(.), 0, .)))
all_combined_data %>%
select(contains("CNI")) %>%
rename(
"Mindset" = mindsetCNI,
"Personality6" = pty6CNI,
"Isms" = ismCNI,
"Fanaticism" = fantcCNI,
"GLOBE" = globeCNI,
"GLOBE (Extra)" = globeExCNI,
"Moral Foundations" = mfCNI,
"Social Axioms" = socAxCNI
) %>%
cor() %>%
ggcorrplot(., hc.order = FALSE,
type = "lower",
outline.col = "white",
colors = c("#3633ff", "white", "#c0392b"),
lab = TRUE, # Show correlation coefficients
ggtheme = ggplot2::theme_bw,
legend.title = "Corr") +
theme(text = element_text(family = "serif", size = 10))
```
## Person level slopes
```{r DONT-RUN-person-slopes, eval=FALSE, echo=TRUE}
#computing person level range
all_cni_models <- lst(mindsetCNI_mod, pty6CNI_mod, ismCNI_mod, fantcCNI_mod, globeCNI_mod, globeExCNI_mod, mfCNI_mod, socAxCNI_mod)
compute_slopes <- function(model_obj) {
obj_name <- deparse(substitute(model_obj))
# Calculate slopes
mod_slopes <- slopes(model_obj, variables = "z_ctry_response", by = "pID")
return(mod_slopes)
}
all_slopes <- map_dfr(all_cni_models, compute_slopes, .id = "model_name")
```
```{r save-person-slopes, echo = F, eval = F}
save(all_slopes,
file = here("objects/Study2/person-slopes-cni.Rdata"))
```
```{r load-person-slopes, echo = F}
#takes time to load
load(here("objects/Study2/person-slopes-cni.Rdata"))
```
```{r range-CNI, results='asis', echo=FALSE}
#| tbl-cap: "Range of CNI estimates across CNI-types"
#| cache: true
tibble(all_slopes) %>%
group_by(model_name) %>%
summarize(
low = min(estimate),
high = max(estimate),
range = high - low) %>%
ungroup() %>%
rename(model = model_name) %>%
arrange(desc(range)) %>%
mutate(model = case_when(
model == "mindsetCNI_mod" ~ "Mindset CNI",
model == "ismCNI_mod" ~ "Isms CNI",
model == "globeCNI_mod" ~ "Globe CNI",
model == "globeExCNI_mod" ~ "Globe Extra CNI",
model == "mfCNI_mod" ~ "Moral Foundations CNI",
model == "fantcCNI_mod" ~ "Fanaticism CNI",
model == "pty6CNI_mod" ~ "Personality 6 CNI",
model == "socAxCNI_mod" ~ "Social Axioms CNI",
TRUE ~ model)) %>%
kable(booktabs = T,
escape = F,
format = "latex",
col.names = c("CNI type", "Min value", "Max value", "Range"),
longtable = T) %>%
kable_styling(latex_options = c("repeat_header"))
```
## CNI Models' Summary
```{r DONT-RUN-create-tbl-cni-mods, echo = T, eval = F}
CNImods_table2 <- modelsummary(
list("Mindset" = mindsetCNI_mod,
"Ism" = ismCNI_mod,
"Globe" = globeCNI_mod,
"Globe (Extra)" = globeExCNI_mod,
"Fanaticism" = fantcCNI_mod,
"Moral Foundations" = mfCNI_mod,
"Social Axioms" = socAxCNI_mod,
"Big 6" = pty6CNI_mod),
estimate = "{estimate}{stars}",
statistic = "[{conf.low}, {conf.high}]",
# group = group + term ~ model,
coef_omit = "Intercept|^status",
coef_rename = c(
"z_ctry_response" = "CNI",
"overall.M" = "Overall Mean Profile",
"SD (z_ctry_response pIDcountry)" = "SD of CNI (across people)",
"SD (z_ctry_response country)" = "SD of CNI (across country)",
"SD (Observations)" = ""),
gof_map = c("nobs"),
fmt = "%.2f",
title = "Cultural Normativity Index",
output = "kableExtra") %>%
add_header_above(c(" ", "Mindset" = 7, "Personality" = 1))
```
```{r save-CNItable, echo = F, eval = F}
save(CNImods_table2,
file = here("objects/Study2/CNImods_table2.Rdata"))
```
```{r load-CNItable, echo = F}
load(here("objects/Study2/CNImods_table2.Rdata"))
```
```{r tbl-cni-mods, results='asis', echo=F, eval=T}
#|tbl-cap: "Mindset and Personality Cultural Normativity Index (CNI) Models"
#| cache: true
CNImods_table2
```
## Country trends
```{r DONT-RUN-cni.ctry-trends, eval=F, echo=T}
#CNI country estimates from `slopes` (also can get CI's here) and `coef` (no CI's) are very close, both provide total effects (fixed+random)
#cogCNI_mod_ctry.slopes <- slopes(cogCNI_mod , variables = "z_ctry_response", by = "country")
all_mods_ctry.trends <- map_dfr(lst(
mindsetCNI_mod, ismCNI_mod, fantcCNI_mod, globeCNI_mod, globeExCNI_mod, mfCNI_mod, socAxCNI_mod, pty6CNI_mod),
~slopes(., variables = "z_ctry_response", by = "country"), .id = "model")
```
```{r save-country-trends, echo = F, eval = F}
save(all_mods_ctry.trends,
file = here("objects/Study2/CNI-ctry-trends2.Rdata"))
```
```{r load-country-trends, echo = F}
load(here("objects/Study2/CNI-ctry-trends2.Rdata"))
```
```{r tbl-ctry.trends-allCNImods, results='asis', echo=F, eval=T}
#| tbl-cap: "Country estimates across CNI types"
#| cache: true
group_names <- c("Mindset CNI", "Isms-Mindset CNI", "Fanaticism-Mindset CNI", "Globe-Mindset CNI", "Globe Extra-Mindset CNI", "Moral Foundations-Mindset CNI", "Social Axioms-Mindset CNI", "Personality6 CNI")
all_mods_ctry.trends %>%
tibble() %>%
select(model, country, estimate, conf.low, conf.high, p.value) %>%
mutate(
p.value = map_chr(p.value, papaja::printp),
across(where(is.numeric), ~ round(., 3)),
ci = paste0("[", conf.low, ", ", conf.high,"]"),
estimate = as.character(estimate)) %>%
select(country, estimate, ci, p.value) %>%
kable(booktabs = T,
escape = F,
format = "latex",
longtable = TRUE,
col.names = c("Country", "Estimate", "95\\% CI", "$\\textit{p}$")) %>%
kable_styling(latex_options = c("repeat_header", "longtable")) %>%
{reduce(seq_along(group_names),
.init = .,
~ group_rows(.x, group_names[.y], 27*(.y-1)+1, 27*.y))} %>%
column_spec(1, width = "5cm")
```
```{r tbl-ctry-CNI_ci.range, results='asis', echo=F, eval=T}
#| label: tbl-ctry-CNI_ci.range
#| tbl-cap: "Confidence Interval Range across CNI types"
#| cache: true
all_mods_ctry.trends %>%
tibble() %>%
select(model, country, estimate, conf.low, conf.high, p.value) %>%
mutate(range_CI = round(conf.high - conf.low, 4),
model = str_extract(model, "^[^_]+")) %>%
mutate(model = case_when(
model == "mindsetCNI" ~ "Mindset CNI",
model == "pty6CNI" ~ "Personality 6 CNI",
model == "pty7CNI" ~ "Personality 7 CNI",
model == "ismCNI" ~ "Isms-Mindset CNI",
model == "fantcCNI" ~ "Fanaticism-Mindset CNI",
model == "globeCNI" ~ "Globe-Mindset CNI",
model == "globeExCNI" ~ "Globe-Extra-Mindset CNI",
model == "mfCNI" ~"Moral Foundations-Mindset CNI",
model == "socAxCNI" ~ "Social Axioms-Mindset CNI",
TRUE ~ model
)) %>%
select(model, range_CI) %>%
distinct() %>%
arrange(desc(range_CI)) %>%
kable(booktabs = T,
escape = F,
format = "latex",
longtable = TRUE,
col.names = c("Model", "CI Range")) %>%
kable_styling(latex_options = c("repeat_header", "longtable"))
```
```{r fun.-oneCoef-plot, echo = F, eval = F}
#One coefficient plot at a time!
plotCNI_ctry_trends <- \(slope_data){
obj_name <- deparse(substitute(slope_data))
title <- str_extract(obj_name, "^[^_]+")
plot <- slope_data %>%
as_tibble() %>%
select(country, estimate, conf.low, conf.high) %>%
ggplot(aes(x = reorder(country,estimate), y = estimate)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
color = "grey") +
geom_point() +
scale_y_continuous(limits = c(0, 1)) +
labs(x = NULL,
y = "CNI",
title = paste0(title, " estimates across countries")) +
coord_flip() +
theme_bw() +
theme(plot.title.position = "plot",
text = element_text(family = "serif", size = 12))
return(plot)
}
#plotCNI_ctry_trends(mindsetCNI_mod_ctry.slopes)
```
```{r ctry-trends-2models-plot, echo=F, eval=F}
#| cache: true
#| fig-height: 12
#| fig-width: 10
#| out-width: "100%"
#| fig-cap-location: top
#| fig-cap: "CNI estimates across countries: Two CNIs"
all_mods_ctry.trends %>%
as_tibble() %>%
mutate(model = str_extract(model, "^[^_]+")) %>%
filter(model %in% c("mindsetCNI", "pty6CNI")) %>%
select(model, country, estimate, conf.low, conf.high) %>%
ggplot(aes(x = reorder(country, estimate), y = estimate,
color = model)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
position = position_dodge(width = 0.5),
width = 0.3) +
geom_point(position = position_dodge(width = 0.5),
size = 1.7)+
scale_y_continuous(limits = c(0, 1)) +
scale_color_manual(values = c("mindsetCNI" = "#0074B7",
"pty6CNI" = "#ff70a6")) +
labs(x = NULL,
y = "CNI",
color = "Model") +
coord_flip() +
theme_bw() +
theme(plot.title.position = "plot",
legend.position = "bottom",
text = element_text(family = "serif", size = 12))
```
```{r ctry-trends-8models-plot, echo=F}
#| fig-height: 12
#| fig-width: 10
#| out-width: "100%"
#| fig-cap-location: top
#| fig-cap: "CNI estimates across countries: Eight CNIs"
#| cache: true
all_mods_ctry.trends %>%
as_tibble() %>%
mutate(model = str_extract(model, "^[^_]+")) %>%
select(model, country, estimate, conf.low, conf.high) %>%
mutate(alpha = ifelse(model %in% c("pty6CNI", "mindsetCNI"), 1, 0.5),
errorbar_size = ifelse(model %in% c("pty6CNI", "mindsetCNI"), .8, 0.4)) %>%
ggplot(aes(x = reorder(country, estimate), y = estimate, color = model)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high, alpha = alpha, size = errorbar_size), position = position_dodge(width = 0.5), width = 2.8) +
geom_point(aes(alpha = alpha, shape = model),
position = position_dodge(width = 0.5),
size = 1.7,
na.rm = TRUE) +
scale_shape_manual(values =
c("mindsetCNI" = 16, # solid circle
"pty6CNI" = 17, # solid triangle
"ismCNI" = 1, # hollow circle
"fantcCNI" = 2, # hollow triangle
"globeCNI" = 5, # hollow diamond
"globeExCNI" = 6, # hollow square
"mfCNI" = 0, # hollow square
"socAxCNI" = 4), # hollow dot
labels = c("mindsetCNI" = "Mindset-CNI",
"ismCNI" = "Ism-CNI",
"fantcCNI" = "Fanaticism-CNI",
"globeCNI" = "GLOBE-CNI",
"globeExCNI" = "GLOBE Extra-CNI",
"mfCNI" = "Moral Foundations-CNI",
"socAxCNI" = "Social Axioms-CNI",
"pty6CNI" = "Personality6-CNI")) +
scale_y_continuous(limits = c(0, 1)) +
scale_color_manual(
values = c("mindsetCNI" = "#0070BB",
"ismCNI" = "#1ca9c9",
"fantcCNI" = "#0000CD",
"globeCNI" = "#A3C1AD",
"globeExCNI" = "#008080",
"mfCNI" = "#8000ff",
"socAxCNI" = "#a69cac",
"pty6CNI" = "#DE3163"),
labels = c("mindsetCNI" = "Mindset-CNI",
"ismCNI" = "Ism-CNI",
"fantcCNI" = "Fanaticism-CNI",
"globeCNI" = "GLOBE-CNI",
"globeExCNI" = "GLOBE Extra-CNI",
"mfCNI" = "Moral Foundations-CNI",
"socAxCNI" = "Social Axioms-CNI",
"pty6CNI" = "Personality6-CNI")
) +
scale_alpha_identity() +
scale_size_identity() +
#scale_linetype_manual(values = c("globeCNI" ="dashed", "globeExCNI" ="dashed","ismCNI" ="dashed", "fantcCNI" ="dashed",
#"mfCNI" = "dashed", "socAxCNI" = "dashed",
# "mindsetCNI" = "solid", "pty6CNI" = "solid")) +
labs(x = NULL, y = "CNI",
color = "Model",
shape = "Model") +
coord_flip() +
theme_bw() +
theme(plot.title.position = "plot",
legend.position = "bottom",
text = element_text(family = "serif", size = 13),
axis.text.y = element_text(size = 13),
legend.text = element_text(size = 11.5)) +
guides(linetype = "none",
alpha = "none",
size = "none",
color = guide_legend(override.aes = list(size = 3)),
shape = guide_legend(override.aes = list(size = 3)))
```
# RQ 2: Associations between CNI and personality and mindset variables
```{r fun.-cni-ContMod-effects}
# Moderators here have to continuous, this won't work for categorical moderators
#Explanation for function below
#[which CNI-type]:
#cni_profiles_data <- mindsetCNI_profiles1$self.ctry_profiles
#[which moderator]:
#modtr_pattern <- "^scP.*indv$" (scP = pty; scM = mindset)
#[causal relationship]:
#CNI <- Continuous moderator
#This version of the function is specifically necessary for parallel processing! Here the modtr_patterns can be varied. The cognitive profiles are long and there are multiple moderators that need to be saved. This method works the best for this.
cni.ContMod_effects <- function(cni_profiles_data, modtr_pattern) {
# Determine the level based on the modtr_pattern
level <- if(grepl("indv", modtr_pattern)) "Individual" else if(grepl("comm", modtr_pattern)) "Country" else "Unknown"
#right data format
cni.modtr_data <- data2 %>%
#select moderator scales individual level (pID) scores
#---$select moderators$ - specify regex
select(pID, matches(modtr_pattern)) %>%
#mean centering the moderator score cols
#since this standardization is col-wise, running pty7 is not an issue
mutate_if(is.numeric,
~. - mean(., na.rm = TRUE)) %>%
pivot_longer(-pID,
names_to = "modtr", #moderator name
values_to = "modtr_score") %>% # moderator score
full_join(cni_profiles_data, by = c("pID"), relationship = "many-to-many")
#CNI <- moderator
cni.modtr_models <- cni.modtr_data %>%
tidyr::nest(data = -modtr) %>%
mutate(
#cni ~ mod
mod.linear = map(data, ~lmer(z_response ~ z_ctry_response * modtr_score + (-1 + z_ctry_response | pID), data = .)),
#cni ~ mod + mod^2
mod.quadratic = map(data, ~lmer(z_response ~ z_ctry_response * modtr_score + z_ctry_response * I(modtr_score^2) + (-1 + z_ctry_response | pID), data = .)),
#save results
tidy_mod.linear = map(mod.linear, tidy, conf.int = TRUE),
tidy_mod.quadratic = map(mod.quadratic, tidy, conf.int = TRUE)
)
results_linear <- cni.modtr_models %>%
select(modtr, tidy_mod.linear) %>%
unnest(cols = c(tidy_mod.linear)) %>%
filter(term == "z_ctry_response:modtr_score") %>%
select(modtr, term, estimate, conf.low, conf.high, p.value)
results_quadratic <- cni.modtr_models %>%
select(modtr, tidy_mod.quadratic) %>%
unnest(cols = c(tidy_mod.quadratic)) %>%
filter(term == "z_ctry_response:I(modtr_score^2)") %>%
select(modtr, term, estimate, conf.low, conf.high, p.value)
output <- bind_rows(
!!level := results_linear,
!!level := results_quadratic,
.id = "level"
)
return(output)
}
#cni.ContMod_effects(pty6CNI_profiles$self.ctry_profiles, "scP_Consc.indv")
```
```{r table_contMods, eval = F, echo=FALSE, results='asis'}
df <- data.frame(
`Mindset` = c(
"[45indv + 46comm]",
paste0(c("Mindset", "Isms", "Globe", "Globe (Extra items)", "Fanaticism", "Moral Foundations", "Social Axioms", "Personality 6"), " CNI ", "x Mindset variables")),
`Personality` = c(
"[8indv + 8comm]",
paste0(c("Mindset", "Isms", "Globe", "Globe (Extra items)", "Fanaticism", "Moral Foundations", "Social Axioms", "Personality 6"), " CNI ", "<- Personality variables"))
)
table_contMods <- kable(df,
format = "latex", align = c("l", "l"),
col.names = c("CNI ~ Mindset variables ", "CNI ~ Personality variables")) %>%
kable_styling(full_width = FALSE) %>%
row_spec(0, bold = TRUE,
) %>%
row_spec(1) %>%
add_header_above(c(" " = 2), bold = TRUE)
```
Given the large *n*, especially for mindsetCNI profiles and the total number of moderators, it will be necessary to use parallel processing for this analysis. It significantly reduces computation time for multiple independent tasks. Instead of splitting the data itself, which is common in parallel processing, we split the list of individual difference variables (e.g., those matching the pattern "\^scMC.\*indv\$"). We create chunks of these variables, distributed across the number of cores - 1. The function `cni.ContMod_effects` function is applied to each variable in `var_names` simultaneously across different cores. Each core processes a different moderator variable independently.
The mindsetCNI models were run on the Talapas (120GB RAM and 30 cores). For the rest of teh CNI types the hardware used involved 64GB RAM and 12 core processor, with the longest processing time for a model being 25mins (1529.822 sec elapsed). A motivated reader may choose to try to optimize this even further, in a more constrained environment.
```{r parallel-processing-contMod, echo = F, eval = F}
# Load necessary libraries
library(furrr)
library(purrr)
library(tictoc)
# Set up parallel processing
plan(multisession, workers = availableCores() - 1)
tic()
# Get the list of moderator variables
patterns <- c("^scP.*indv$", "^scP.*comm$")
var_names <- grep(paste(patterns, collapse = "|"), names(data2), value = TRUE)
#checklist
#1. check if the patterns are correct
#["^scM" for mindset moderators; "^scP" for personality moderators]
#2.change the name of the object
#3. change the data_profiles
#4. save the object in Study 2
# Parallel execution using furrr
pty6CNI.pty_effects<- future_map_dfr(var_names, function(var) {
result <- cni.ContMod_effects(pty6CNI_profiles$self.ctry_profiles, var)
result$variable <- var # Add the variable name to the result
return(result)
}, .progress = TRUE)
# Reset to sequential processing
plan(sequential)
toc()
```
```{r save-cni.ContMod-models, echo = F, eval = F}
#checklist
#1. change name of the object
#2. DANGER: change name of the objet being saved
#use the following to save each or collected chunks in the code above
save(pty6CNI.pty_effects,
file = here("objects/Study2/pty6CNI.pty_effects.Rdata"))
```
```{r load-cni.ContMod-models, echo = F, results='hide'}
# Define a vector of file names
cni.contMod_object_names <- paste0(
rep(c("mindset", "ism", "globe", "globeEx", "fantc", "mf", "socAx", "pty6"), each = 2),
"CNI.",
rep(c("mindset", "pty"), 8),
"_effects.Rdata")
# Set file paths and set names
filepaths <- here("objects/Study2", cni.contMod_object_names) %>%
set_names(nm = basename(.) %>% tools::file_path_sans_ext())
# Load all files
all_cni.contMod_effects <- purrr::map(filepaths, ~get(load(.)))
# Assign loaded objects to global environment - important!
purrr::pmap(.l = list(.x = names(all_cni.contMod_effects), .y = all_cni.contMod_effects), .f = ~assign(.x, .y, envir = .GlobalEnv))
# Verify that objects are now in the global environment
#map_lgl(names(files), ~exists(.x, envir = .GlobalEnv))
```
```{r fun.-specify-var-names}
#function to specify the correct variable names for all moderators [Pty + Cog + swb]
specify_var_names <- \(data){
data %>%
mutate(
modtr = case_when(
#Personality
str_detect(modtr, "Consc") ~ "Conscientiousness",
str_detect(modtr, "Hon") ~ "Honesty",
str_detect(modtr, "Agr") ~ "Agreeableness",
str_detect(modtr, "Res") ~ "Resilience",
str_detect(modtr, "Ext") ~ "Extraversion",
str_detect(modtr, "Vir") ~ "Originality/Virtuosity",
str_detect(modtr, "Dis") ~ "Disintegration",
str_detect(modtr, "SocialSelfReg") ~ "Social Self-Regulation",
str_detect(modtr, "Dynamism") ~ "Dynamism",
#Mindset
#Isms
str_detect(modtr, "Isms_TradRelig") ~ "Traditional Religiousness",
str_detect(modtr, "Isms_UnmiSelfInt") ~ "Unmitigated Self-Interest",
str_detect(modtr, "CommRation") ~ "Communal Rationalism",
str_detect(modtr, "SubjSpirit") ~ "Subjective Spirituality",
str_detect(modtr, "IneqAver") ~ "Inequality Aversion",
#GLOBE
str_detect(modtr, "globe_igc") ~ "InGroup Collectivism",
str_detect(modtr, "globe_ic") ~ "Institutional Collectivism",
str_detect(modtr, "globe_pf") ~ "Performance/Achievement Orientation",
str_detect(modtr, "globe_un") ~ "Uncertainty Avoidance",
str_detect(modtr, "globe_ge") ~ "Gender Egalitarian",
str_detect(modtr, "globe_po") ~ "Power Distance",
str_detect(modtr, "globe_fu") ~ "Future Orientation",
str_detect(modtr, "globe_ass") ~ "Assertiveness",
str_detect(modtr, "globe_hu") ~ "Humane Orientation",
#Social Axioms
str_detect(modtr, "SocAx_Cynicism") ~ "Social Cynicism",
str_detect(modtr, "SocAx_RewardForApp") ~ "Reward for Application",
str_detect(modtr, "SocAx_Complexity") ~ "Social Complexity",
str_detect(modtr, "SocAx_FateControl") ~ "Fate Control",
str_detect(modtr, "SocAx_Religiosity") ~ "Religiosity",
#Moral Foundations
str_detect(modtr, "MFQ_Harm") ~ "Harm",
str_detect(modtr, "MFQ_Fairness") ~ "Fairness",
str_detect(modtr, "MFQ_Loyalty") ~ "Loyalty",
str_detect(modtr, "MFQ_Authority") ~ "Authority",
str_detect(modtr, "MFQ_Purity") ~ "Purity",
#Individualism-Collectivism (Gelfland)
str_detect(modtr, "IndvColl_CollecHori") ~ "Horizontal Collectivism",
str_detect(modtr, "IndvColl_CollecVert") ~ "Vertical Collectivism",
str_detect(modtr, "IndvColl_IndivHori") ~ "Horizontal Individualism",
str_detect(modtr, "IndvColl_IndivVert") ~ "Vertical Individualism",
#Family Values
str_detect(modtr, "FamValCohesion") ~ "Cohesion",
str_detect(modtr, "FamValHierarchy") ~ "Hierarchy",
#Schwartz Values
str_detect(modtr, "Value_SfTran") ~ "Self Transcendence",
str_detect(modtr, "Value_Consrv") ~ " Conservation",
str_detect(modtr, "Value_SelfEnh") ~ "Self Enhancement",
str_detect(modtr, "Value_OpnChg") ~ "Openness to Change",
#Fanaticism
str_detect(modtr, "fntc") ~ "Fanaticism",
#Materialism
str_detect(modtr, "Material") ~ "Materialism",
#Machiavellianism
str_detect(modtr, "Machiaveln") ~ "Machiavellianism",
#Prone to Aggression (Honor)
str_detect(modtr, "Prone_Aggr") ~ "Prone to Aggression",
#Ethnonationalism
str_detect(modtr, "Ethnonational") ~ "Ethnonationalism",
#Tightness-Looseness
str_detect(modtr, "Tightness") ~ "Tightness",
#Amoralism
str_detect(modtr, "Amoralism") ~ "Amoralism",
#Duke Religiosity Index
str_detect(modtr, "DRI_OrgRelgAc") ~ "Organizational Religious Activity",
str_detect(modtr, "DRI_NonOrgRelgAc") ~ "Non-organizational Religious Activity",
str_detect(modtr, "DRI_IntrnscRelg") ~ "Intrinsic religiosity",
#Civic Multicultural Nationalism
str_detect(modtr, "CivicMCNational") ~ "Civic Multicultural Nationalism",
TRUE ~ modtr
)
)
}
```
## Descriptives - Effect size estimates
```{r overall-ContEst-trends, results='asis', echo=FALSE}
#| tbl-cap: "Descriptive Statistics for Continuos Moderator Effect size Estimates"
#| cache: true
library(e1071)
list_rbind(all_cni.contMod_effects, names_to = "model") %>%
mutate(level = case_when(
level == "Community" ~ "Country",
TRUE ~ level
),
term = case_when(
term == "z_ctry_response:modtr_score" ~ "Linear",
term == "z_ctry_response:I(modtr_score^2)" ~ "Quadratic",
TRUE ~ term
)) %>%
group_by(term, level) %>%
summarise(
min = min(estimate),
max = max(estimate),
range = max - min,
kurtosis = e1071::kurtosis(estimate),
skew = e1071::skewness(estimate),
n = n(),
mean = mean(estimate),
sd = sd(estimate),
.groups = "drop"
) %>%
mutate(across(where(is.numeric), ~round(., 2))) %>%
select(-term) %>%
kable(booktabs = T,
escape = F,
format = "latex",
col.names = c(" ", "min", "max", "range", "kurtosis", "skew", "n", "mean", "sd")) %>%
group_rows("Linear effects", 1, 2) %>%
group_rows("Quadratic effects", 3, 4)
```
```{r overall-ContEst-trends-plot, eval = T, echo = F, results='asis'}
#| fig-cap: "Distribution of Continuos Moderator Effect size Estimates"
#| fig-width: 8
#| fig-height: 4
#| cache: true
plot_cont_effects <- function(plot_type) {
base_data <- list_rbind(all_cni.contMod_effects, names_to = "model") %>%
mutate(level = case_when(
level == "Community" ~ "Country",
TRUE ~ level
))
plot_data <- switch(plot_type,
"overall" = base_data,
"linear" = base_data %>% filter(term == "z_ctry_response:modtr_score"),
"quadratic" = base_data %>% filter(term == "z_ctry_response:I(modtr_score^2)"),
stop("Invalid plot type. Use 'overall', 'linear', or 'quadratic'.")
)
plot <- ggplot(plot_data, aes(x = estimate, y = after_stat(density))) +
geom_density(color = "grey40", fill = "grey80", alpha = 0.5) +
geom_density(aes(fill = level, color = level), alpha = 0.5) +
scale_fill_manual(values = c("Individual" = "cornflowerblue", "Country" = "pink")) +
scale_color_manual(values = c("Individual" = "cornflowerblue", "Country" = "pink")) +
scale_x_continuous(limits = c(-0.85, 1)) +
scale_y_continuous(limits = c(0, 280))+
theme_bw() +
labs(x = "Effect size Estimates", y = "Density", fill = "Level", color = "Level")
return(plot)
}
plot_cont_effects("overall") +
plot_cont_effects("linear") +
plot_cont_effects("quadratic") +
plot_layout(
axis_titles = "collect",
guides = "collect"
) +
plot_annotation(
tag_levels = list(c("Overall", "Linear", "Quadratic"))
) &
theme(
plot.tag.position = c("top"),
text = element_text(family = "serif", size = 12),legend.position = "bottom")
```
## Tables
```{r cni.contMod-tbl-details}
CNIplot_data_details = tibble(
mod_names = c(
"mindsetCNI_plot",
"ismCNI_plot",
"globeCNI_plot",
"globeExCNI_plot",
"fantcCNI_plot",
"mfCNI_plot",
"socAxCNI_plot",
"pty6CNI_plot"
),
caption = c(
"Mindset CNI",
"Isms-Mindset CNI",
"Globe-Mindset CNI",
"Globe (Extra items)",
"Fanaticism-Mindset CNI",
"Moral Foundations-Mindset CNI",
"Social Axioms-Mindset CNI",
"Personality6 CNI"
)
)
cni.contMod_object_names <- c(
"mindsetCNI.mindset_effects.Rdata",
"mindsetCNI.pty_effects.Rdata",
"ismCNI.mindset_effects.Rdata",
"ismCNI.pty_effects.Rdata",
"globeCNI.mindset_effects.Rdata",
"globeCNI.pty_effects.Rdata",
"globeExCNI.mindset_effects.Rdata",
"globeExCNI.pty_effects.Rdata",
"fantcCNI.mindset_effects.Rdata",
"fantcCNI.pty_effects.Rdata",
"mfCNI.mindset_effects.Rdata",
"mfCNI.pty_effects.Rdata",
"socAxCNI.mindset_effects.Rdata",
"socAxCNI.pty_effects.Rdata",
"pty6CNI.mindset_effects.Rdata",
"pty6CNI.pty_effects.Rdata"
)
cni.contMod_tbl_details <- tibble(
object_names = names(all_cni.contMod_effects)) %>%
mutate(obj_name = str_extract(object_names, "^[^.]+")) %>%
left_join(
CNIplot_data_details %>%
mutate(mod_name = str_remove(mod_names, "_plot$")),
by = c("obj_name" = "mod_name")) %>%
select(object_names, caption) %>%
mutate(
CNI_type = case_when(
str_detect(object_names, ".mindset") ~ "Mindset",
str_detect(object_names, ".pty") ~ "Personality",
TRUE ~ NA_character_ # for any cases that don't match either condition
)) %>%
mutate(caption = paste0("Associations between ",caption, " and ", CNI_type, " variables")) %>%
select(object_names, caption, CNI_type)
```
### Linear and Quadratic effects (CNI)
```{r fun.-cni.ContMod-effects-tbl}
#| cache: true
cni.ContMod_effects_table <- \(data, title){
df2 <- data %>%
#specify the correct variable/moderator name
specify_var_names()%>%
select(-relation) %>%
mutate(
p.value = map_chr(p.value, papaja::printp),
across(where(is.numeric), ~ round(., 3)),
ci = paste0("[", conf.low, ", ", conf.high,"]"),
term = ifelse(str_detect(term, "2"), "Quadratic", "Linear"),
estimate = as.character(estimate)
) %>%
select(level, modtr, term, estimate, ci, p.value) %>%
pivot_longer(cols = c(estimate, ci, p.value), names_to = "stat", values_to = "stat_value") %>%
unite(term, term, stat) %>%
pivot_wider(names_from = term, values_from = stat_value) %>%
arrange(desc(level)) %>%
select(-level) %>%
select(modtr, contains("Linear"), contains("Quadratic"))
# Calculate the number of rows for each level
no.modtrs <- df2 %>%
distinct(modtr) %>%
nrow()
df2 %>%
kable(booktabs = T,
escape = F,
format = "latex",
longtable = TRUE,
col.names = c(" ", "Est", "95\\% CI", "p", "Est", "95\\% CI", "p"),
caption = title ) %>%
kable_styling(latex_options = c("repeat_header", "longtable")) %>%
add_header_above(c(" ", "Linear (b)" = 3, "Quadratic ($b^2$)" = 3),
escape = FALSE,
line = TRUE) %>%
group_rows("Individual level", 1, no.modtrs, bold = TRUE, italic = T) %>%
group_rows("Country level", no.modtrs+1, no.modtrs*2, bold = TRUE, italic = T) %>%
column_spec(1, width = "5.5cm")
}
```
```{r cni.MindsetMod-effects-tbl, results='asis'}
#| cache: true
cni.contMod_tbl_details %>%
filter(CNI_type == "Mindset") %>%
pwalk(function(object_names, caption, ...) {
cat("\n\n")
print(
pluck(all_cni.contMod_effects, object_names) %>%
cni.ContMod_effects_table(., caption)
)
cat("\n\n")
})
```
```{r cni.PersonalityMod-effects-tbl, results='asis'}
#| cache: true
cni.contMod_tbl_details %>%
filter(CNI_type == "Personality") %>%
pwalk(function(object_names, caption, ...) {
cat("\n\n")
print(
pluck(all_cni.contMod_effects, object_names) %>%
cni.ContMod_effects_table(., caption)
)
cat("\n\n")
})
```
### Strongest and Weakest effects (CNI)
CNI focused view across personality and mindset traits as predictors (moderators of CNI)
```{r funk-strong.weak_effects}
#| cache: true
strong.weak_effects_table <- function(data, arrange_for = c("estimate"), title) {
# Validate and get the arrangement variable
arrange_for <- match.arg(arrange_for) # This ensures arrange_for is one of the options
base_data <- data %>%
specify_var_names() %>%
filter(round(p.value, 3) < 0.001) %>%
arrange(desc(abs(!!sym(arrange_for)))) %>%
mutate(
term_type = case_when(
term == "z_ctry_response:modtr_score" ~ "Linear",
term == "z_ctry_response:I(modtr_score^2)" ~ "Quadratic",
TRUE ~ NA_character_
)
)
# Select 15 highest and 15 lowest estimates
highest_15 <- base_data %>% head(15)
lowest_15 <- base_data %>% tail(15) %>%
arrange(abs(!!sym(arrange_for)))
# Combine the data using bind_rows with .id argument
all_data <- bind_rows(list("Highest 15 Estimates" = highest_15,
"Lowest 15 Estimates" = lowest_15),
.id = "group")
# Prepare data for table
df2 <- all_data %>%
select(-relation) %>%
mutate(
p.value = map_chr(p.value, papaja::printp),
across(where(is.numeric), ~ round(., 3)),
ci = paste0("[", conf.low, ", ", conf.high,"]"),
term = ifelse(str_detect(term, "2"), "($b\\textsuperscript{2}$)", "(b)"),
estimate = as.character(estimate),
modtr = ifelse(level == "Country",
paste0(modtr, footnote_marker_symbol(1, "latex")),
paste0(modtr, footnote_marker_symbol(2, "latex")))) %>%
select(group, modtr, term, estimate, ci, p.value)%>%
arrange(group) %>%
select(-group)
# Create and return the table
df2 %>%
kable(booktabs = T,
escape = F,
format = "latex",
longtable = TRUE,
col.names = c("Moderator", "Term", "Est", "95\\% CI", "$\\textit{p}$" ),
caption = title) %>%
group_rows("Strongest 15 Estimates", 1, 15, underline = T, escape = F) %>%
group_rows("Weakest 15 Estimates", 16, 30, escape = F, underline = T) %>%
kable_styling(latex_options = c("repeat_header", "longtable")) %>%
column_spec(1, width = "6.4cm") %>%
footnote(symbol = c("Country level", "Individual level"),
symbol_title = "",
general = "(b) = Linear model, ($b\\textsuperscript{2}$) = Quadratic",
general_title = "Note:",
title_format = "italic",
fixed_small_size = TRUE,
footnote_as_chunk = FALSE,
threeparttable = TRUE,
escape = FALSE)
}
#strong.weak_effects_table(mindsetCNI.mindset_effects, "mindsetCNI.mindset_effects")
cni.contMod_tbl_details <- cni.contMod_tbl_details %>%
mutate(caption2 = paste0("Strongest and Weakest ", caption),
#change arrange_var to arrange according to another variable
arrange_var = "estimate")
```
```{r cni.contMod-strong.weak-effects-tbl, results='asis'}
#| cache: true
cni.contMod_tbl_details %>%
select(-caption) %>%
pwalk(function(object_names, arrange_var, caption2, ...) {
cat("\n\n")
print(
pluck(all_cni.contMod_effects, object_names) %>%
strong.weak_effects_table(., arrange_var, caption2)
)
cat("\n\n")
})
```
### Overall trends
Similar to study 1, three categories of tables are presented here: 1. Strongest moderators: Effects that significantly enhance or diminish the impact of cultural normativity. 2. Weakest moderators: Effects that, while statistically significant, have minimal influence on cultural normativity. 3. Non-significant effects: Factors that show no statistically significant moderation of cultural normativity.
```{r fun.-overall-effects-table}
overall_table <- function(table_type, arrange_for = c("estimate"), title) {
arrange_for <- match.arg(arrange_for) # This ensures arrange_for is one of the options
table <- bind_rows(all_cni.contMod_effects, .id = "model")
if (table_type == "overall highest") {
table <- table %>%
filter(p.value < .001) %>%
arrange(desc(abs(!!sym(arrange_for))))
} else if (table_type == "overall weakest") {
table <- table %>%
filter(p.value < .001) %>%
arrange(abs(!!sym(arrange_for)))
} else if (table_type == "overall n.s") {
table <- table %>%
filter(p.value >= .05) %>%
arrange(abs(!!sym(arrange_for)))
} else {
stop("Invalid table_type. Choose 'overall highest', 'overall lowest', or 'overall n.s'")
}
final_table <- table %>%
head(20) %>%
specify_var_names() %>%
separate(model, c("CNI_type", "Moderator_type"), sep = "\\.") %>%
mutate(
CNI_type = str_remove(CNI_type, "CNI"),
CNI_type = str_to_sentence(CNI_type),
CNI_type = str_replace_all(CNI_type, c(
"cni" = " CNI",
"Ism" = "Isms",
"Fantc" = "Fanaticism",
"Globe" = "GLOBE",
"Pty6" = "Personality 6",
"ex" = "-Extra",
"Mf"= "Moral Foundations",
"Socax" = "Social Axioms")),
Moderator_type = if_else(Moderator_type == "pty_effects", "Personality", "Mindset"),
p.value = map_chr(p.value, papaja::printp),
across(where(is.numeric), ~ round(., 3)),
ci = paste0("[", conf.low, ", ", conf.high,"]"),
term = ifelse(str_detect(term, "2"), "($b\\textsuperscript{2}$)", "(b)"),
estimate = as.character(estimate),
modtr = paste0(modtr, " (", Moderator_type, ")"),
modtr = ifelse(level == "Country",
paste0(modtr, footnote_marker_symbol(1, "latex")),
paste0(modtr, footnote_marker_symbol(2, "latex"))
),
CNI_type = str_replace(CNI_type, "CNI", "")) %>%
select(CNI_type, modtr, term, estimate, ci, p.value)
final_table %>%
kable(booktabs = T,
escape = F,
format = "html",
col.names = c("CNI type", "Moderator", "Term", "Est", "95\\% CI", "$\\textit{p}$"),
caption = title,
longtable = T) %>%
kable_styling(latex_options = c("repeat_header", "longtable")) %>%
footnote(symbol = c("Country", "Individual"),
symbol_title = "",
general = "(b) = Linear model, ($b\\textsuperscript{2}$) = Quadratic model",
general_title = "Note:",
title_format = "italic",
fixed_small_size = TRUE,
footnote_as_chunk = FALSE,
threeparttable = TRUE,
escape = FALSE)
}
overall_table_details <- tibble(
table_type = c("overall highest", "overall weakest", "overall n.s"))%>%
mutate(title = paste0("Variables with ", word(table_type , 2), " effect sizes with each CNI type"),
#change arrange_var to arrange according to another variable
arrange_var = "estimate")
```
```{r overall-effects-table, echo = F, results='asis'}
#| cache: true
overall_table_details %>%
pwalk(function(table_type, arrange_var, title) {
cat("\n\n")
print(overall_table(table_type, arrange_var, title))
cat("\n\n")
})
```
### Results section
#### Overall trends across levels and relationship type
```{r fun.-overall-groupwise-effects-table}
#| cache: true
#| cache-rebuild: true
# resolves the issue of comparing apples vs. oranges
overall_groupwise_table <- function(level_type, arrange_for = c("estimate"), relation_type) {
arrange_for <- match.arg(arrange_for) # This ensures arrange_for is one of the options
table <- bind_rows(all_cni.contMod_effects, .id = "model") %>% #View()
mutate(term = case_when(
term == "z_ctry_response:modtr_score" ~ "Linear",
term == "z_ctry_response:I(modtr_score^2)" ~ "Quadratic",
TRUE ~ term # keep original value if no match
))
#level_type = "Individual";relation_type = "Linear"; arrange_for = "estimate"
final_table <- table %>%
filter(level == level_type & term == relation_type ) %>%
filter(p.value<.001) %>%
#arrange(desc(abs(estimate))) %>%
arrange(desc(abs(!!sym(arrange_for)))) %>%
head(20) %>%
specify_var_names() %>%
separate(model, c("CNI_type", "Moderator_type"), sep = "\\.") %>%
mutate(
CNI_type = str_remove(CNI_type, "CNI"),
CNI_type = str_to_sentence(CNI_type),
CNI_type = str_replace_all(CNI_type, c(
"cni" = " CNI",
"Ism" = "Isms",
"Fantc" = "Fanaticism",
"Globe" = "GLOBE",
"Pty6" = "Personality 6",
"ex" = "-Extra",
"Mf"= "Moral Foundations",
"Socax" = "Social Axioms")),
Moderator_type = if_else(str_detect(Moderator_type, "pty"), "Personality", "Mindset"),
p.value = map_chr(p.value, papaja::printp),
across(where(is.numeric), ~ round(., 3)),
ci = paste0("[", conf.low, ", ", conf.high,"]"),
modtr = ifelse(Moderator_type == "Personality",
paste0(modtr, footnote_marker_symbol(1, "latex")),
paste0(modtr, footnote_marker_symbol(2, "latex"))),
estimate = as.character(estimate),
CNI_type = str_replace(CNI_type, "CNI", "")) %>%
select(CNI_type, modtr, estimate, ci, p.value)
final_table %>%
kable(booktabs = T,
escape = F,
format = "html",
col.names = c("CNI type", "Moderator trait", "Est", "95\\% CI", "$\\textit{p}$"),
caption = paste0("Strongest ", relation_type, " associations with CNI at the ", level_type, " level")) %>%
kable_styling() %>%
footnote(symbol = c("Personality trait", "Mindset trait"),
symbol_title = "",
general_title = "Note:",
title_format = "italic",
fixed_small_size = TRUE,
footnote_as_chunk = FALSE,
threeparttable = TRUE)
}
overall_groupwise_table_details <- tibble(
level_type = rep(c("Individual", "Country"), each = 2),
relation_type = rep(c("Linear", "Quadratic"), 2),
#change arrange_var to arrange according to another variable
arrange_var = "estimate")
#overall_groupwise_table("Individual","estimate", "Quadratic")
```
```{r overall-groupwise-effects-table, echo = F, results='asis'}
#| cache: true
#| cache-rebuild: true
overall_groupwise_table_details %>%
pwalk(function(level_type, arrange_var, relation_type) {
cat("\n\n")
print(overall_groupwise_table(level_type, arrange_var, relation_type))
cat("\n\n")
})
```
#### Eight strongest effect plots
```{r fun.-Eight-highest-quad-effects, eval=TRUE, include=TRUE}
#| cache: true
cni.ContMod_fig <- \(cni_profiles_data, modtr_pattern, x_title, y_title, modtr_type = c("Personality", "Mindset")){
# Match moderator type and set default if not matched
modtr_type <- match.arg(modtr_type)
# Define colors based on moderator type
line_color <- switch(modtr_type,
"Personality" = "#c1000d",
"Mindset" = "darkgreen")
ribbon_color <- switch(modtr_type,
"Personality" = "rosybrown2",
"Mindset" = "darkseagreen2")
# Evaluate cni_profile string as object name
cni_profiles_data <- get(cni_profiles_data)
# Right data format
data <- data2 %>%
select(pID, matches(modtr_pattern)) %>%
mutate_if(is.numeric,
~. - mean(., na.rm = T)) %>%
pivot_longer(-pID,
names_to = "modtr",
values_to = "modtr_score") %>%
full_join(cni_profiles_data[["self.ctry_profiles"]], by = c("pID"), relationship = "many-to-many")
# Fit the model
model <- lmer(z_response ~ z_ctry_response * modtr_score +
z_ctry_response * I(modtr_score^2) +
(-1 + z_ctry_response | pID),
data = data)
# Extract fixed effects and variance-covariance matrix
fe <- fixef(model)
vcov_mat <- vcov(model)
# Create sequence of moderator scores for prediction
# Now using n_points parameter
mod_seq <- seq(min(data$modtr_score, na.rm = TRUE),
max(data$modtr_score, na.rm = TRUE),
length.out = 500)
# Create design matrix for predictions
X_pred <- cbind(1, # Intercept
mod_seq, # Linear term
mod_seq^2) # Quadratic term
# Calculate predicted slopes and standard errors
pred_slopes <- fe["z_ctry_response"] +
fe["z_ctry_response:modtr_score"] * mod_seq +
fe["z_ctry_response:I(modtr_score^2)"] * mod_seq^2
# Calculate standard errors for the slope predictions
slope_vcov <- vcov_mat[c("z_ctry_response",
"z_ctry_response:modtr_score",
"z_ctry_response:I(modtr_score^2)"),
c("z_ctry_response",
"z_ctry_response:modtr_score",
"z_ctry_response:I(modtr_score^2)")]
se_slopes <- sqrt(diag(X_pred %*% slope_vcov %*% t(X_pred)))
# Create prediction data frame with confidence intervals
pred_data <- data.frame(
modtr_score = mod_seq,
predicted_slope = pred_slopes,
ci_lower = pred_slopes - 1.96 * se_slopes,
ci_upper = pred_slopes + 1.96 * se_slopes
)
# Extract random effects and calculate total slopes for individual points
ranef_values <- ranef(model)$pID[,"z_ctry_response"]
slope_data <- data %>%
distinct(pID, modtr_score) %>%
mutate(
total_slope = fe["z_ctry_response"] +
fe["z_ctry_response:modtr_score"] * modtr_score +
fe["z_ctry_response:I(modtr_score^2)"] * modtr_score^2 +
ranef_values[as.character(pID)]
)
# Add timing measurement
start_time <- Sys.time()
# Create the plot
p <- ggplot() +
geom_ribbon(data = pred_data,
aes(x = modtr_score,
ymin = ci_lower,
ymax = ci_upper),
fill = ribbon_color,
alpha = 0.3) +
geom_point(data = slope_data,
aes(x = modtr_score, y = total_slope),
alpha = 0.3,
color = "gray50") +
geom_line(data = pred_data,
aes(x = modtr_score, y = predicted_slope),
color = line_color,
size = 1) +
labs(x = x_title,
y = y_title) +
# scale_x_continuous(limits = c(-.7, .3)) +
#scale_y_continuous(limits = c(.4, .75)) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5),
panel.grid.minor = element_blank(),
text = element_text(family = "serif", size = 10),
axis.text.x = element_text(size = 7),
axis.text.y = element_text(size = 7)
)
# Remove n_points from timing message since it's fixed at 500
end_time <- Sys.time()
cat(sprintf("Plot generation took %.2f seconds\n",
as.numeric(end_time - start_time)))
return(p)
}
cni.ContMod_fig_emm <- \(cni_profiles_data, modtr_pattern, x_title, y_title, modtr_type = c("Personality", "Mindset")){
# Match moderator type and set default if not matched
modtr_type <- match.arg(modtr_type)
# Define colors based on moderator type
line_color <- switch(modtr_type,
"Personality" = "#c1000d",
"Mindset" = "darkgreen")
ribbon_color <- switch(modtr_type,
"Personality" = "rosybrown2",
"Mindset" = "darkseagreen2")
# Evaluate cni_profile string as object name
cni_profiles_data <- get(cni_profiles_data)
# Right data format
data <- data2 %>%
select(pID, matches(modtr_pattern)) %>%
mutate_if(is.numeric,
~. - mean(., na.rm = T)) %>%
pivot_longer(-pID,
names_to = "modtr",
values_to = "modtr_score") %>%
full_join(cni_profiles_data[["self.ctry_profiles"]],
by = c("pID"),
relationship = "many-to-many")
# Fit the model
model <- lmer(z_response ~ z_ctry_response * modtr_score +
z_ctry_response * I(modtr_score^2) +
(-1 + z_ctry_response | pID),
data = data)
# Get sequence for moderator scores
mod_seq <- seq(min(data$modtr_score, na.rm = TRUE),
max(data$modtr_score, na.rm = TRUE),
length.out = 500)
# Get slopes and confidence intervals using emtrends
slopes <- emmeans::emtrends(
model,
specs = "modtr_score",
var = "z_ctry_response",
at = list(modtr_score = mod_seq)
) %>%
as.data.frame()
# Create prediction dataframe
pred_data <- data.frame(
modtr_score = mod_seq,
predicted_slope = slopes$z_ctry_response.trend,
ci_lower = slopes$z_ctry_response.trend - 1.96 * slopes$SE,
ci_upper = slopes$z_ctry_response.trend + 1.96 * slopes$SE
)
# Extract random effects for individual points
ranef_values <- ranef(model)$pID[,"z_ctry_response"]
slope_data <- data %>%
distinct(pID, modtr_score) %>%
mutate(
total_slope = pred_data$predicted_slope[findInterval(modtr_score, mod_seq)] +
ranef_values[as.character(pID)]
)
# Create the plot
p <- ggplot() +
geom_ribbon(data = pred_data,
aes(x = modtr_score,
ymin = ci_lower,
ymax = ci_upper),
fill = ribbon_color,
alpha = 0.3) +
geom_point(data = slope_data,
aes(x = modtr_score, y = total_slope),
alpha = 0.3,
color = "gray50") +
geom_line(data = pred_data,
aes(x = modtr_score, y = predicted_slope),
color = line_color,
size = 1) +
labs(x = x_title,
y = y_title) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5),
panel.grid.minor = element_blank(),
text = element_text(family = "serif", size = 10),
axis.text.x = element_text(size = 7),
axis.text.y = element_text(size = 7))
return(p)
}
highest_eight_effects_details <- bind_rows(all_cni.contMod_effects, .id = "model") %>%
mutate(term = case_when(
term == "z_ctry_response:modtr_score" ~ "Linear",
term == "z_ctry_response:I(modtr_score^2)" ~ "Quadratic",
TRUE ~ term # keep original value if no match
)) %>%
#filter(level == level_type & term == relation_type ) %>%
filter(level == "Country" & term == "Quadratic" ) %>%
filter(p.value < .001) %>%
arrange(desc(abs(estimate))) %>%
filter(abs(estimate) > .1) %>%
head(8) %>%
mutate(modtr_name = modtr) %>%
specify_var_names() %>%
separate(model, c("CNI_type", "Moderator_type"), sep = "\\.") %>%
mutate(cni_profile = paste0(CNI_type, "_profiles")) %>%
dplyr::select(cni_profile,modtr_name, CNI_type, Moderator_type, level, modtr, term) %>%
mutate(
# CNI_type = str_remove(CNI_type, "CNI"),
CNI_type = str_to_sentence(CNI_type),
CNI_type = str_replace_all(CNI_type, c(
"cni" = "-CNI",
"Ism" = "Isms",
"Fantc" = "Fanaticism",
"Globe" = "GLOBE",
"Pty6" = "Personality6",
"Mf" = "Moral Foundations",
"ex" = " Extra")),
Moderator_type = if_else(Moderator_type == "pty_effects", "Personality", "Mindset"))
plot_list <- highest_eight_effects_details %>% #pull(modtr_name)
pmap(function(cni_profile, modtr_name, CNI_type, Moderator_type, level, modtr, term) {
cni.ContMod_fig(
cni_profiles_data = cni_profile,
modtr_pattern = modtr_name,
x_title = modtr,
y_title = CNI_type,
modtr_type = Moderator_type
)
})
combined_plot <- wrap_plots(plot_list, ncol = 2) +
plot_annotation(
caption = "Variable domain: Personality and Mindset predictors"
) &
theme(
plot.caption = element_markdown(family = "serif",hjust = 0, size = 10)
)
```
```{r Eight-highest-quad-effects, eval=TRUE, include=TRUE, results='asis'}
#| cache: true
#| fig-height: 10
#| fig-width: 7
#| out-width: "100%"
#| fig-cap-location: top
#| fig-cap: "Eight Strongest Country-level Nonlinear Effects with CNI"
combined_plot
```
#### Twenty Highest Quadratic Effects
```{r fun.-quad-slopes, eval=FALSE, include=TRUE, warning=FALSE}
#| cache: true
emm_options(lmerTest.limit = 176090)
cni.ContMod_quad_slopes <- \(cni_profiles_data, modtr_pattern, CNI_type, modtr_type = c("Personality", "Mindset")){
# Evaluate cni_profile string as object name
cni_profiles_data <- get(cni_profiles_data)
# Right data format
data <- data2 %>%
select(pID, matches(modtr_pattern)) %>%
mutate_if(is.numeric,
~. - mean(., na.rm = T)) %>%
pivot_longer(-pID,
names_to = "modtr",
values_to = "modtr_score") %>%
full_join(cni_profiles_data[["self.ctry_profiles"]],
by = c("pID"),
relationship = "many-to-many")
# Fit the model
model <- lmer(z_response ~ z_ctry_response * modtr_score +
z_ctry_response * I(modtr_score^2) +
(-1 + z_ctry_response | pID),
data = data)
# Get sequence for moderator scores
mod_seq <- seq(min(data$modtr_score, na.rm = TRUE),
max(data$modtr_score, na.rm = TRUE),
length.out = 500)
# Get slopes and confidence intervals using emtrends
slopes <- emmeans::emtrends(
model,
specs = "modtr_score",
var = "z_ctry_response",
at = list(modtr_score = mod_seq)
) %>%
as.data.frame()
#slopes$z_ctry_response.trend
# Create prediction dataframe
pred_data <- tibble(
CNI_type ={{CNI_type}},
modtr = {{modtr_pattern}},
modtr_type = modtr_type,
modtr_score = mod_seq, # x-axis
predicted_slope = slopes$z_ctry_response.trend ) #y-axis; predicted CNI
return(pred_data)
}
highest_20_effects_details <- bind_rows(all_cni.contMod_effects, .id = "model") %>%
mutate(term = case_when(
term == "z_ctry_response:modtr_score" ~ "Linear",
term == "z_ctry_response:I(modtr_score^2)" ~ "Quadratic",
TRUE ~ term )) %>%
filter(level == "Country" & term == "Quadratic" ) %>%
filter(p.value < .001) %>%
arrange(desc(abs(estimate))) %>%
head(20) %>%
mutate(modtr_name = modtr) %>%
specify_var_names() %>%
separate(model, c("CNI_type", "Moderator_type"), sep = "\\.") %>%
mutate(cni_profile = paste0(CNI_type, "_profiles")) %>%
dplyr::select(cni_profile,modtr_name, CNI_type, Moderator_type, level, modtr, term) %>%
mutate(
CNI_type = str_remove(CNI_type, "CNI"),
CNI_type = str_to_sentence(CNI_type),
CNI_type = str_replace_all(CNI_type, c(
"cni" = "-CNI",
"Ism" = "Isms",
"Fantc" = "Fanaticism",
"Globe" = "GLOBE",
"Pty6" = "Personality6",
"Mf" = "MoralFoundation",
"GLOBEex" = "GLOBE(Extra)" )),
Moderator_type = if_else(Moderator_type == "pty_effects", "Personality", "Mindset"))
quad_slopes_data <- highest_20_effects_details %>%
select(cni_profile, modtr_name, CNI_type, Moderator_type) %>%
pmap(function(cni_profile, modtr_name, CNI_type, Moderator_type) {
cni.ContMod_quad_slopes(
cni_profiles_data = cni_profile,
modtr_pattern = modtr_name,
CNI_type = CNI_type,
modtr_type = Moderator_type
)
})
```
```{r save-slopes, echo = F, eval = F}
save(quad_slopes_data,
file = here("objects/Study2/quad_slopes_data.Rdata"))
```
```{r echo = F}
load("~/Desktop/CNI_Dissertation/objects/Study2/quad_slopes_data.Rdata")
```
```{r quad-slopes, eval=TRUE, echo=FALSE, results='asis'}
#| cache: true
#| fig-height: 10
#| fig-width: 7
#| out-width: "100%"
#| fig-cap-location: top
#| fig-cap: "Twenty Strongest Country-level Nonlinear Effects with CNI"
plot_data <- quad_slopes_data %>%
list_c() %>%
specify_var_names()%>%
mutate(
# Add "CNI" to CNI_type for y-axis label
plot_label = paste(CNI_type, "CNI"),
# Create shortened trait names
trait_short = case_when(
modtr == "Dynamism" ~ "D",
modtr == "Social Self-Regulation" ~ "SSR",
modtr == "Openness to Change" ~ "OpnChange",
modtr == "Self Transcendence" ~ "SelfTran",
TRUE ~ modtr
),
CNI_type = str_replace_all(CNI_type, c( "Socax" = "SocialAxioms" )) #,
# modtr_type = factor(modtr_type, levels = c("Personality", "Mindset"))
)
# Get end points for text labels
data_ends <- plot_data %>%
group_by(CNI_type, modtr, modtr_type, trait_short) %>%
slice(n()) %>%
ungroup()
library(ggrepel)
ggplot(plot_data, aes(x = modtr_score,
y = predicted_slope)) +
geom_line(aes(color = modtr_type, group = modtr), show.legend = TRUE) +
geom_text_repel(
data = data_ends,
aes(label = trait_short, color = modtr_type),
segment.color = "grey50",
segment.size = 0.3,
nudge_x = 0.2,
direction = "y",
hjust = 0,
show.legend = FALSE # Hide legend for text labels
) +
# Use facet_wrap with custom y-axis labels
facet_wrap(~CNI_type,
scales = "free_y",
ncol = 2,
strip.position = "left") +
scale_color_manual(
values = c(
"Personality" = "#c1000d",
"Mindset" = "darkgreen"),
drop= FALSE) +
scale_x_continuous(limits = c(-1.9, 1.5)) +
scale_y_continuous(limits = c(.3, .75)) +
theme_bw() +
theme(
text = element_text(family = "serif", size = 17.5),
legend.position = "bottom",
panel.grid.minor = element_blank(),
strip.background = element_blank(),
strip.placement = "outside",
strip.text = element_text(angle = 0, hjust = 0.5),
axis.title.y = element_blank()
) +
labs(
x = "Attribute Score",
color = "Attribute Type")
```
#### R-squared
Adding R-squared for each models. These were computed to make comparisons across models across linear and quadratic variations for every moderator variable. Camparing b and b^2 across linear and quadratic models would mean comparing apples and oranges, computing r^2 addresses this issue.
Marginal R²: This represents the variance explained by the fixed effects only (z_ctry_response, modtr_score, and their interaction)
Conditional R²: This represents the variance explained by both fixed and random effects (including random effect of z_ctry_response varying by participant/pID)
```{r compute-R2-funk}
# Build a function similar to the one built for computing continuous effects
r2_cni.ContMod_effects <- \(cni_profiles_data, modtr_pattern){
level <- if(grepl("indv", modtr_pattern)) "Individual" else if(grepl("comm", modtr_pattern)) "Country" else "Unknown"
cni.modtr_data <- data2 %>%
#select moderator scales individual level (pID) scores
#---$select moderators$ - specify regex
select(pID, matches(modtr_pattern)) %>%
#mean centering the moderator score cols
#since this standardization is col-wise, running pty7 is not an issue
mutate_if(is.numeric,
~. - mean(., na.rm = TRUE)) %>%
pivot_longer(-pID,
names_to = "modtr", #moderator name
values_to = "modtr_score") %>% # moderator score
full_join(cni_profiles_data, by = c("pID"), relationship = "many-to-many")
mod.linear = lmer(z_response ~ z_ctry_response * modtr_score + (-1 + z_ctry_response | pID), data = cni.modtr_data)
library(performance)
r2_linear <- r2(mod.linear) %>%
unlist() %>%
as_tibble_row() %>%
rename_with(~str_remove(., "\\..*")) %>%
mutate(relation = "linear")
#cni ~ mod + mod^2
mod.quadratic = lmer(z_response ~ z_ctry_response * modtr_score + z_ctry_response * I(modtr_score^2) + (-1 + z_ctry_response | pID), data = cni.modtr_data)
r2_quadratic <- r2(mod.quadratic) %>%
unlist() %>%
as_tibble_row() %>%
rename_with(~str_remove(., "\\..*")) %>%
mutate(relation = "quadratic")
output <- bind_rows(
!!level := r2_linear ,
!!level := r2_quadratic ,
.id = "level"
) %>%
mutate(modtr = modtr_pattern)
return(output)
}
#Parallel analysis for all models - even if this takes time it won't take up a lot of RAM space
#Add the results to the CNI-wise objects
```
```{r parallel-processing-R2-contMod, echo = F, eval = F}
# Load necessary libraries
library(furrr)
library(purrr)
library(tictoc)
# Set up parallel processing
plan(multisession, workers = availableCores() - 1)
tic()
# Get the list of moderator variables
patterns <- c("^scM.*indv$", "^scM.*comm$")
var_names <- grep(paste(patterns, collapse = "|"), names(data2), value = TRUE)
#checklist
#1. check if the patterns are correct
#["^scM" for mindset moderators; "^scP" for personality moderators]
#2.change the name of the object
#3. change the data_profiles
#4. save the object in Study 2
# Parallel execution using furrr
socAxCNI.mindset_R2effects<- future_map_dfr(var_names, function(var) {
result <- r2_cni.ContMod_effects(socAxCNI_profiles$self.ctry_profiles, var)
result$variable <- var # Add the variable name to the result
return(result)
}, .progress = TRUE)
# Reset to sequential processing
plan(sequential)
toc()
```
```{r echo=F, eval=F}
#-----LOAD ALL PREVIOUS CONTINUOUS EFFECT .RDATA
# Define a vector of file names
cni.contMod_object_names <- paste0(
rep(c("mindset", "ism", "globe", "globeEx", "fantc", "mf", "socAx", "pty6"), each = 2),
"CNI.",
rep(c("mindset", "pty"), 8),
"_effects.Rdata")
# Set file paths and set names
filepaths <- here("objects/Study2", cni.contMod_object_names) %>%
set_names(nm = basename(.) %>% tools::file_path_sans_ext())
# Load all files
all_cni.contMod_effects <- purrr::map(filepaths, ~get(load(.)))
# Assign loaded objects to global environment - important!
purrr::pmap(.l = list(.x = names(all_cni.contMod_effects), .y = all_cni.contMod_effects), .f = ~assign(.x, .y, envir = .GlobalEnv))
# Verify that objects are now in the global environment
#map_lgl(names(files), ~exists(.x, envir = .GlobalEnv))
```
```{r addingR2-cont.effects, echo=F, eval=F}
#Add data to the models
#---cni.contMod_object_names:
#| 1 mindsetCNI.mindset_effects.Rdata
#| 2 mindsetCNI.pty_effects.Rdata
#| 3 ismCNI.mindset_effects.Rdata
#| 4 ismCNI.pty_effects.Rdata
#| 5 globeCNI.mindset_effects.Rdata
#| 6 globeCNI.pty_effects.Rdata
#| 7 globeExCNI.mindset_effects.Rdata
#| 8 globeExCNI.pty_effects.Rdata
#| 9 fantcCNI.mindset_effects.Rdata
#| 10 fantcCNI.pty_effects.Rdata
#| 11 mfCNI.mindset_effects.Rdata
#| 12 mfCNI.pty_effects.Rdata
#| 13 socAxCNI.mindset_effects.Rdata
#| 14 socAxCNI.pty_effects.Rdata
#| 15 pty6CNI.mindset_effects.Rdata
#| 16 pty6CNI.pty_effects.Rdata
#---------
new_R2_df <- \(original_contMod_df, R2_contMod_df){
original_contMod_df <-
original_contMod_df %>%
left_join(
R2_contMod_df %>%
select(level, modtr, relation, R2_conditional, R2_marginal),
by = c("level", "modtr", "relation"))
return(original_contMod_df)
}
socAxCNI.mindset_effects <-
new_R2_df(socAxCNI.mindset_effects,
socAxCNI.mindset_R2effects)
```
```{r save-cni.ContMod-models, echo = F, eval = F}
#checklist
#1. change name of the object
#2. DANGER: change name of the object being saved
#use the following to save each or collected chunks in the code above
save(ismCNI.pty_effects,
file = here("objects/Study2/
ismCNI.pty_effects.Rdata"))
```
## Plots
### Strongest and Weakest effects (CNI)
Given that there are `r (length(str_subset(names(data2), "\\.(indv|comm)$"))-1)*2` (individual level scale scores = `r (length(str_subset(names(data2), "\\.(indv)$")))-1`, country level scale scores = `r (length(str_subset(names(data2), "\\.(comm)$")))-1` and their linear and quadratic associations with CNI types) effect sizes to evaluate in total for each CNI type, it was necessary to identify a starting point. In order to assess the continuous moderator effects for each of the CNI types, the highest 15 and lowest 15 associations (across linear and quadratic) were considered.
```{r fun.plot-high-low-effects}
plot_high.low_mindset_effects <- function(data) {
base_data <- data %>%
specify_var_names() %>%
filter(round(p.value, 3) < 0.05) %>% #distinct(level) %>%
arrange(desc(abs(estimate))) %>%
mutate(
color = case_when(
str_detect(variable, "MFQ") ~ "Moral Foundations",
str_detect(variable, "SocAx") ~ "Social Axioms",
str_detect(variable, "globe") ~ "Globe",
str_detect(variable, "Isms|CommRation|SubjSpirit|IneqAver") ~ "Isms",
str_detect(variable, "IndvColl") ~ "Individualism-Collectivism",
str_detect(variable, "FamVal") ~ "Family Values",
str_detect(variable, "Value") ~ "Schwartz Values",
str_detect(variable, "fntc") ~ "Fanaticism",
str_detect(variable, "Material") ~ "Materialism",
str_detect(variable, "Machiaveln") ~ "Machiavellianism",
str_detect(variable, "Prone_Aggr") ~ "Prone to Aggression",
str_detect(variable, "Ethnonational") ~ "Ethnonationalism",
str_detect(variable, "Tightness") ~ "Tightness-Looseness",
str_detect(variable, "Amoralism") ~ "Amoralism",
str_detect(variable, "DRI") ~ "Duke Religiosity Index",
str_detect(variable, "CivicMCNational") ~ "Civic Multicultural Nationalism",
TRUE ~ "OTHER"),
term_type = case_when(
term == "z_ctry_response:modtr_score" ~ "Linear",
term == "z_ctry_response:I(modtr_score^2)" ~ "Quadratic",
TRUE ~ NA_character_,
.default = NA_character_)
)
# Select 15 highest and 15 lowest estimates and tag them as 'Highest' or 'Lowest'
highest_15 <- base_data %>% head(15) %>% mutate(group = "Highest 15 Estimates")
lowest_15 <- base_data %>% tail(15) %>% mutate(group = "Lowest 15 Estimates")
# Combine the data
all_data <- bind_rows(highest_15, lowest_15)
# Define the overall range of estimates for consistent axis scaling
estimate_range <- range(c(all_data$conf.high, all_data$conf.low))
# Define the plot
combined_plot <- ggplot(all_data,
aes(x = reorder(modtr, estimate), y = estimate, color = color, linetype = term_type, shape = level, group = interaction(modtr, term_type, level))) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, alpha = 0.3, size = 3, position = position_dodge(width = 0.4), linetype = "solid") +
geom_point(size = 2, position = position_dodge(width = 0.4)) +
geom_linerange(aes(ymin = 0, ymax = estimate), position = position_dodge(width = 0.4)) +
scale_shape_manual(values = c("Country" = 15, "Individual" = 17), drop = FALSE) +
scale_color_manual(values = c(
"Isms" = "#7F0000",
"Globe" = "#377EB8",
"Social Axioms" = "#4DAF4A",
"Moral Foundations" = "#984EA3",
"Individualism-Collectivism" = "#FF7F00",
"Family Values" = "#DAA520",
"Schwartz Values" = "#5A3CB4",
"Fanaticism" = "#196569",
"Materialism" = "#F781BF",
"Machiavellianism" = "#8f8d2f",
"Prone to Aggression" = "#E41A1C",
"Ethnonationalism" = "#C67B53",
"Tightness-Looseness" = "#8B4513",
"Amoralism" = "#00FFFF",
"Duke Religiosity Index" = "#FF1493",
"Civic Multicultural Nationalism" = "#1da87e"
), drop = FALSE) +
scale_linetype_manual(values = c("Linear" = "dotdash", "Quadratic" = "solid"), drop = FALSE) +
scale_y_continuous(limits = estimate_range) +
coord_flip() +
facet_wrap(~group, ncol = 1, scales = "free_y")+
theme_minimal(base_family = "serif", base_size = 12) +
labs(
x = NULL,
y = "Estimate",
color = "Domain",
shape = "Level",
linetype = "Association") +
theme_bw() +
theme(text = element_text(family = "serif", size = 12), legend.position = "bottom",
legend.text = element_text(size = 8.1),
legend.title = element_text(size = 8.5))+
guides(
color = guide_legend(ncol = 4, title.position = "top"), # 2 columns for Domain
shape = guide_legend(ncol = 1, title.position = "top"), # 1 column for Level
linetype = guide_legend(ncol = 1, title.position = "top") # 1 column for Association
)
return(combined_plot)
}
plot_high.low_pty_effects <- function(data) {
base_data <- data %>%
# cni.ContMod_tables$ismCNI.pty %>%
specify_var_names() %>%
filter(round(p.value, 3) < 0.05) %>%
arrange(desc(abs(estimate))) %>%
mutate(
color = case_when(
modtr == "Conscientiousness" ~ "Big 6",
modtr == "Honesty" ~ "Big 6",
modtr == "Agreeableness" ~ "Big 6",
modtr == "Resilience" ~ "Big 6",
modtr == "Extraversion" ~ "Big 6",
modtr == "Originality/Virtuosity" ~ "Big 6",
modtr == "Social Self-Regulation"~ "Big 2",
modtr == "Dynamism"~ "Big 2",
TRUE ~ "Other"),
term_type = case_when(
term == "z_ctry_response:modtr_score" ~ "Linear",
term == "z_ctry_response:I(modtr_score^2)" ~ "Quadratic",
TRUE ~ NA_character_,
.default = NA_character_)
)
# Select 15 highest and 15 lowest estimates and tag them as 'Highest' or 'Lowest'
highest_15 <- base_data %>% head(15) %>% mutate(group = "Highest 15 Estimates")
lowest_15 <- base_data %>% tail(15) %>% mutate(group = "Lowest 15 Estimates")
# Combine the data
all_data <- bind_rows(highest_15, lowest_15)
# Define the overall range of estimates for consistent axis scaling
estimate_range <- range(c(all_data$conf.low, all_data$conf.high))
# Define the plot
combined_plot <- ggplot(all_data,
aes(x = reorder(modtr, estimate), y = estimate, color = color, linetype = term_type, shape = level, group = interaction(modtr, term_type, level))) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, alpha = 0.3, size = 3, position = position_dodge(width = 0.4), stat = "identity", linetype = "solid") +
geom_point(aes(shape = level), size = 2, position = position_dodge(width = 0.4)) +
geom_linerange(aes(ymin = 0, ymax = estimate, linetype = term_type), position = position_dodge(width = 0.4)) +
scale_shape_manual(values = c("Country" = 15, "Individual" = 17), drop = FALSE) +
scale_color_manual(values = c(
"Big 2" = "#377EB8",
"Big 6" = "#008f7a")) +
scale_linetype_manual(values = c("Linear" = "dotdash", "Quadratic" = "solid"), drop = FALSE) +
scale_y_continuous(limits = estimate_range) + # Ensure consistent y-axis limits
coord_flip() +
facet_wrap(~group, ncol = 1, scales = "free_y") + # Facet the plot into Highest and Lowest
theme_minimal(base_family = "serif", base_size = 12) +
labs(
x = NULL,
y = "Estimate",
color = "Domain",
fill = "Domain",
shape = "Level",
linetype = "Association") +
theme_bw() +
theme(
text = element_text(family = "serif", size = 12), legend.position = "bottom") +
guides(
color = guide_legend(ncol = 2, title.position = "top"), # 2 columns for Domain
fill = guide_legend(ncol = 2, title.position = "top"), # 2 columns for Domain (for the ribbon)
shape = guide_legend(ncol = 1, title.position = "top"), # 1 column for Level
linetype = guide_legend(ncol = 1, title.position = "top") # 1 column for Association
)
return(combined_plot)
}
```
```{r echo=F}
#| cache: true
#| fig-height: 10
#| fig-width: 7
#| out-width: "100%"
#| fig-cap-location: top
#| fig-cap: !expr |-
#| cni.contMod_tbl_details %>%
#| filter(CNI_type == "Mindset") %>%
#| mutate(caption = paste0("Overall trends: ",
#| str_extract(caption, "Associations.*"))) %>%
#| pull(caption)
cni.contMod_tbl_details %>%
filter(CNI_type == "Mindset") %>%
pull(object_names) %>%
walk(~ {
cat("\n\n")
# Extract the table and generate the plot
table_data <- pluck(all_cni.contMod_effects, .x )
plot <-
plot_high.low_mindset_effects(table_data)
print(plot)
cat("\n\n")
})
```
```{r echo=F}
#| cache: true
#| fig-height: 8
#| fig-width: 7
#| out-width: "100%"
#| fig-cap-location: top
#| fig-cap: !expr |
#| cni.contMod_tbl_details %>%
#| filter(CNI_type == "Personality") %>%
#| mutate(caption = paste0("Overall trends: ", str_extract(caption, "Associations.*"))) %>%
#| pull(caption)
cni.contMod_tbl_details %>%
filter(CNI_type == "Personality") %>%
pull(object_names) %>%
walk(~ {
cat("\n\n")
# Extract the table and generate the plot
table_data <- pluck(all_cni.contMod_effects, .x )
plot <-
plot_high.low_pty_effects(table_data)
# Print the plot
print(plot)
cat("\n\n")
})
```
### Linear and Quadratic effects
```{r fun.-cni.ContMod-effects-plot}
cni.pty_effects_plot <-\(data, association_type){
if (association_type == "linear") {
data_ass.type <- data %>% #cni.ContMod_tables$globeCNI.pty7
#linear effects
filter(term == "z_ctry_response:modtr_score")
} else if (association_type == "quadratic") {
data_ass.type <- data %>%
#quadratic effects
filter(term == "z_ctry_response:I(modtr_score^2)")
} else {
stop("Invalid association type. Choose 'linear' or 'quadratic'")
}
processed_data <- data_ass.type %>%
specify_var_names() %>%
with_groups(modtr, ~ mutate(.x, avg_abs_estimate = mean(abs(estimate)))) %>%
mutate(
modtr = fct_reorder(modtr, avg_abs_estimate, .desc = FALSE),
domain = case_when(
modtr == "Conscientiousness" ~ "Big 6",
modtr == "Honesty" ~ "Big 6",
modtr == "Agreeableness" ~ "Big 6",
modtr == "Resilience" ~ "Big 6",
modtr == "Extraversion" ~ "Big 6",
modtr == "Originality/Virtuosity" ~ "Big 6",
modtr == "Social Self-Regulation"~ "Big 2",
modtr == "Dynamism"~ "Big 2",
TRUE ~ "Other"),
sig = case_when(
p.value < .05 & estimate > 0 ~ "Pos",
p.value < .05 & estimate < 0 ~ "Neg",
TRUE ~ "NOT"),
colors = case_when(
domain == "Big 2" ~ "#377EB8",
domain == "Big 6" ~ "#008f7a")
)
base_plot <- ggplot(processed_data, aes(x = modtr, y = estimate, fill = domain, alpha = level)) +
geom_col(position = position_dodge(width = 0.8)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
position = position_dodge(width = 0.8),
width = 0.25) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
scale_fill_manual(values = setNames(processed_data$colors, processed_data$domain)) +
scale_alpha_manual(values = c("Individual" = 1, "Country" = 0.4)) +
labs(x = "Moderator", y = "Estimate",
fill = "Trait", alpha = "Level")
get_color_mapping <- function(data) {
data %>%
arrange(desc(avg_abs_estimate)) %>%
distinct(modtr, .keep_all = TRUE) %>%
select(modtr, colors) %>%
deframe()
}
# Get the color mapping
color_mapping <- get_color_mapping(processed_data)
final_plot <- base_plot +
# domain independent order
theme_bw() +
coord_flip()+
theme(text = element_text(family="serif", size = 11),
legend.position = "bottom",
axis.text.y = element_markdown(
color = color_mapping[levels(processed_data$modtr)],
size = 11))
return(final_plot)
}
cni.mindset_effects_plot <- function(data, association_type, plot_type = c("domain dependent", "domain independent")) {
if (association_type == "linear") {
data_ass.type <- data %>%
#linear effects
filter(term == "z_ctry_response:modtr_score")
} else if (association_type == "quadratic") {
data_ass.type <- data %>%
#quadratic effects
filter(term == "z_ctry_response:I(modtr_score^2)")
} else {
stop("Invalid association type. Choose 'linear' or 'quadratic'")
}
processed_data <- data_ass.type %>%
specify_var_names() %>%
with_groups(modtr, ~ mutate(.x, avg_abs_estimate = mean(abs(estimate)))) %>%
mutate(
modtr = fct_reorder(modtr, avg_abs_estimate, .desc = FALSE),
domain = case_when(
str_detect(variable, "MFQ") ~ "Moral Foundations",
str_detect(variable, "SocAx") ~ "Social Axioms",
str_detect(variable, "globe") ~ "Globe",
str_detect(variable, "Isms|CommRation|SubjSpirit|IneqAver") ~ "Isms",
str_detect(variable, "IndvColl") ~ "Individualism-Collectivism",
str_detect(variable, "FamVal") ~ "Family Values",
str_detect(variable, "Value") ~ "Schwartz Values",
str_detect(variable, "fntc") ~ "Fanaticism",
str_detect(variable, "Material") ~ "Materialism",
str_detect(variable, "Machiaveln") ~ "Machiavellianism",
str_detect(variable, "Prone_Aggr") ~ "Prone to Aggression",
str_detect(variable, "Ethnonational") ~ "Ethnonationalism",
str_detect(variable, "Tightness") ~ "Tightness-Looseness",
str_detect(variable, "Amoralism") ~ "Amoralism",
str_detect(variable, "DRI") ~ "Duke Religiosity Index",
str_detect(variable, "CivicMCNational") ~ "Civic Multicultural Nationalism",
TRUE ~ "OTHER"),
sig = case_when(
p.value < .05 & estimate > 0 ~ "Pos",
p.value < .05 & estimate < 0 ~ "Neg",
TRUE ~ "NOT")
) %>%
mutate(
colors = case_when(
str_detect(domain, "Isms") ~ "#7F0000",
str_detect(domain, "Globe") ~ "#377EB8",
str_detect(domain, "Social Axioms") ~ "#4DAF4A",
str_detect(domain, "Moral Foundations") ~ "#984EA3",
str_detect(domain, "Individualism-Collectivism") ~ "#FF7F00",
str_detect(domain, "Family Values") ~ "#DAA520",
str_detect(domain, "Schwartz Values") ~ "#5A3CB4",
str_detect(domain, "Fanaticism") ~ "#196569",
str_detect(domain, "Materialism") ~ "#F781BF",
str_detect(domain, "Machiavellianism") ~ "#8f8d2f",
str_detect(domain, "Prone to Aggression") ~ "#E41A1C",
str_detect(domain, "Ethnonationalism") ~ "#C67B53",
str_detect(domain, "Tightness-Looseness") ~ "#8B4513",
str_detect(domain, "Amoralism") ~ "#00FFFF",
str_detect(domain, "Duke Religiosity Index") ~ "#FF1493",
str_detect(domain, "Civic Multicultural Nationalism") ~ "#1da87e",
TRUE ~ "#000000"
)
)
base_plot <- ggplot(processed_data, aes(x = modtr, y = estimate, fill = domain, alpha = level)) +
geom_col(position = position_dodge(width = 0.8)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(width = 0.8), width = 0.25) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
scale_fill_manual(
values = setNames(processed_data$colors, processed_data$domain),
guide = guide_legend(
position = "bottom",
ncol = 4,
direction = "vertical")) +
scale_alpha_manual(
values = c("Individual" = 1, "Country" = 0.4),
guide = guide_legend(
position = "bottom",
ncol = 1,
direction = "vertical")) +
labs(x = NULL, y = "Estimate", fill = "Trait", alpha = "Level")
#add here different plot versions
if(plot_type=="domain dependent"){
final_plot <- base_plot +
labs(caption = paste(tools::toTitleCase(plot_type), "view"))+
#domain dependent order
scale_x_discrete(
limits = unique(processed_data$modtr),
labels = function(x) {
sapply(x, function(label) {
color <- processed_data$colors[processed_data$modtr == label][1]
sprintf("%s", color, label)
})
}
) +
theme_bw() +
theme(axis.text.y = element_markdown(size = 11),
text = element_text(family="serif", size = 10)) +
coord_flip()
}else if(plot_type=="domain independent"){
get_color_mapping <- function(data) {
data %>%
arrange(desc(avg_abs_estimate)) %>%
distinct(modtr, .keep_all = TRUE) %>%
select(modtr, colors) %>%
deframe()
}
# Get the color mapping
color_mapping <- get_color_mapping(processed_data)
final_plot <- base_plot +
labs(caption = paste(toTitleCase(plot_type), "view"))+
# domain independent order
theme_bw() +
coord_flip()+
theme(text = element_text(family="serif", size = 10),
axis.text.y = element_markdown(size = 11,
color = color_mapping[levels(processed_data$modtr)]
))
}else{
stop("Invalid plot type. Choose either 'domain dependent' or 'domain independent'.")
}
return(final_plot)
}
cni.contMod_plot_details <-
cni.contMod_tbl_details %>%
split(.$CNI_type) %>%
map_dfr(~ {
if (first(.x$CNI_type) == "Mindset") {
crossing(.x,
view = c("Domain Independent", "Domain Dependent"),
association_type = c("Linear", "Quadratic"))
} else if (first(.x$CNI_type) == "Personality") {
crossing(.x,
association_type = c("Linear", "Quadratic")) %>%
mutate(view = NA_character_)
} else {
# For any other CNI_type, keep rows as is
.x %>% mutate(view = NA_character_, association_type = NA_character_)
}
})
```
```{r walk-cni.mindsetContMod-effects-plot, results='asis', echo=FALSE}
#| cache: true
#| fig-height: 15
#| fig-width: 10
#| out-width: "100%"
#| cap-location: top
#| fig-cap: !expr |
#| cni.contMod_plot_details %>%
#| filter(CNI_type == "Mindset") %>%
#| mutate(caption = paste0(view, " view for ", association_type, " effects: ", caption)) %>%
#| pull(caption)
cni.contMod_plot_details %>%
filter(CNI_type == "Mindset") %>%
mutate(across(c(view, association_type), str_to_lower)) %>%
select(object_names, association_type, view) %>%
pwalk(function(object_names, association_type, view) {
cat("\n\n")
print(
pluck(all_cni.contMod_effects, object_names) %>%
cni.mindset_effects_plot(association_type, view)
)
cat("\n\n")
})
```
```{r walk-cni.ptyContMod-effects-plot, results='asis', echo=FALSE}
#| cache: true
#| fig-height: 8
#| fig-width: 10
#| out-width: "100%"
#| cap-location: top
#| fig-cap: !expr |
#| cni.contMod_plot_details %>%
#| filter(CNI_type == "Personality") %>%
#| mutate(caption = paste0(association_type, " effects: ", caption)) %>%
#| pull(caption)
#retain models that have cognitive moderators
cni.contMod_plot_details %>%
filter(CNI_type == "Personality") %>%
mutate(across(c(association_type), str_to_lower)) %>%
select(object_names, association_type) %>%
pwalk(function(object_names, association_type) {
cat("\n\n")
print(
pluck(all_cni.contMod_effects, object_names) %>%
cni.pty_effects_plot(association_type)
)
cat("\n\n")
})
```
### Five Strongest moderators
```{r five-strongest-moderators}
cni_types <- unique(sub("CNI.*", "CNI", names(all_cni.contMod_effects)))
strongest5_contMods <- function(cni_type, data_list) {
df_names <- grep(paste0("^", cni_type),
names(data_list), value = TRUE)
bind_rows(data_list[df_names], .id = "modtr_type") %>%
arrange(desc(abs(estimate))) %>%
filter(p.value<.05) %>%
slice_head(n = 5)
}
# Apply the function to each CNI type
strongest5_contMods_CNI <-
map_dfr(cni_types, ~strongest5_contMods(.x, all_cni.contMod_effects)) %>%
specify_var_names() %>%
with_groups(modtr, ~ mutate(.x, avg_abs_estimate = mean(abs(estimate)))) %>%
mutate(modtr = fct_reorder(modtr, avg_abs_estimate, .desc = FALSE)) %>%
select(-term, -variable) %>%
mutate(cni_type = str_extract(modtr_type, "^[^.]+")) %>%
mutate(label_color = case_when(
str_detect(modtr_type, "\\.pty_effects$") ~ "#c1000d",
str_detect(modtr_type, "\\.mindset_effects$") ~ "darkgreen",
TRUE ~ "black"
))
modtr_order <- strongest5_contMods_CNI %>%
group_by(modtr) %>%
summarize(mean_estimate = mean(estimate)) %>%
arrange(mean_estimate) %>%
pull(modtr)
strongest5_contMods_CNI <- strongest5_contMods_CNI %>%
mutate(modtr = factor(modtr, levels = modtr_order))
color_mapping <- strongest5_contMods_CNI %>%
select(modtr, label_color) %>%
distinct() %>%
deframe()
```
```{r tbl-five-strongest-moderators, results='asis'}
#| tbl-cap: "Five Strongest Estimates across CNI types"
#| cache: true
final_table <-
map_dfr(cni_types, ~strongest5_contMods(.x, all_cni.contMod_effects)) %>%
specify_var_names() %>%
with_groups(modtr, ~ mutate(.x, avg_abs_estimate = mean(abs(estimate)))) %>%
mutate(modtr = fct_reorder(modtr, avg_abs_estimate, .desc = FALSE)) %>%
specify_var_names() %>%
separate(modtr_type, c("CNI_type", "Moderator_type"), sep = "\\.") %>%
mutate(
CNI_type = str_remove(CNI_type, "CNI"),
CNI_type = str_to_sentence(CNI_type),
CNI_type = str_replace_all(CNI_type, c(
"cni" = " CNI",
"Ism" = "Isms",
"Fantc" = "Fanaticism",
"Globe" = "GLOBE",
"Pty6" = "Personality 6",
"Pty7" = "Personality 7")),
Moderator_type = if_else(Moderator_type == "pty7", "Personality", "Mindset"),
p.value = map_chr(p.value, papaja::printp),
across(where(is.numeric), ~ round(., 3)),
ci = paste0("[", conf.low, ", ", conf.high,"]"),
term = ifelse(str_detect(term, "2"), "(b2)", "(b)"),
estimate = as.character(estimate),
modtr = paste0(modtr, " (", Moderator_type, ")"),
modtr = ifelse(level == "Country",
paste0(modtr, footnote_marker_symbol(1, "latex")),
paste0(modtr, footnote_marker_symbol(2, "latex"))),
CNI_type = str_replace(CNI_type, "CNI", "")) %>%
select(CNI_type, modtr, term, estimate, ci, p.value)
final_table %>%
kable(booktabs = T,
escape = F,
format = "latex",
col.names = c("CNI type", "Moderator", "Term", "Est", "95\\% CI", "$\\textit{p}$"),
longtable = T,
align = c('l', 'r', 'c', 'c', 'c', 'c')) %>%
kable_styling(latex_options = c("repeat_header")) %>%
collapse_rows(columns = 1, latex_hline = "linespace") %>%
footnote(symbol = c("Country", "Individual"),
symbol_title = "",
general = "(b) = Linear model, (b2) = Quadratic",
general_title = "Note:",
title_format = "italic",
fixed_small_size = TRUE,
footnote_as_chunk = FALSE,
threeparttable = TRUE)
```
```{r fig-five-strongest-moderators, results='asis'}
#| fig-height: 7
#| fig-width: 10
#| out-width: "100%"
#| fig-cap-location: top
#| fig-cap: "Five Strongest Estimates across CNI types"
#| cache: true
strongest5_contMods_CNI %>%
mutate(relation = str_to_sentence(relation)) %>%
ggplot(aes(x = modtr, y = estimate, color = cni_type, group = interaction(modtr))) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, alpha = 0.3, size = 3.8, position = position_dodge(width = 0.4), linetype = "solid")+
geom_point(size= 3) +
facet_wrap(~relation)+
scale_shape_manual(
values = c("FALSE" = 4, "TRUE" = 16),
labels = c("FALSE" = expression(paste("Linear model (", italic(beta)["1t"], italic(mu)[11], ")")),
"TRUE" = expression(paste("Quadratic (", italic(beta)["1t"],italic(mu)[12]^2, ")"))),
name = "Term Type"
)+
scale_color_manual(
values = c("mindsetCNI" = "#0070BB",
"ismCNI" = "#1ca9c9",
"fantcCNI" = "#0000CD",
"globeCNI" = "#A3C1AD",
"globeExCNI" = "#008080",
"mfCNI" = "#8000ff",
"socAxCNI" = "#a69cac",
"pty6CNI" = "#DE3163"),
labels = c("mindsetCNI" = "Mindset-CNI",
"ismCNI" = "Ism-CNI",
"fantcCNI" = "Fanaticism-CNI",
"globeCNI" = "GLOBE-CNI",
"globeExCNI" = "GLOBE Extra-CNI",
"mfCNI" = "Moral Foundations-CNI",
"socAxCNI" = "Social Axioms-CNI",
"pty6CNI" = "Personality6-CNI")
)+
labs(x = NULL,
y = "CNI estimate",
color = "CNI Types",
caption = "Variable domain: Personality and Mindset moderators on Y axis") +
coord_flip() +
theme_bw() +
theme(plot.title.position = "plot",
legend.position = "bottom",
text = element_text(family = "serif", size = 13),
legend.box = "vertical",
legend.title = element_text(size = 11),
legend.spacing.y = unit(0.2, "cm"),
legend.margin = margin(t = 0.2, r = 0, b = 0, l = 0),
legend.key.size = unit(0.8, "lines"),
axis.ticks.y = element_blank(),
axis.text.y = element_markdown(face = "italic",
color = color_mapping[levels(strongest5_contMods_CNI$modtr)],
size = 14),
plot.caption = element_markdown(hjust =.5, size = 11))+
guides(
color = guide_legend(order = 1, label.hjust = 0)
)
```
```{r}
#| echo: false
#| include: false
#knitr::knit_exit()
```
# RQ 3: To what extent is CNI related to person characteristics?
The primary goal is to assess whether and how specific categorical moderators (modtr_cat) influence the alignment between individual psychological responses (z_response) and country-level psychological profiles (z_ctry_response). The focus is on examining whether the relationship between `z_response` and `z_ctry_response` (i.e., the slope of `z_ctry_response`) changes depending on the level of `modtr_score`. By including the interaction term `z_ctry_response * modtr_score` in the model, the effect of `z_ctry_response` on `z_response` is allowed to vary across different levels of `modtr_score`.
By modeling this relationship, the aim is to:
- Determine if the degree of cultural normativity varies across levels of the moderator.
- Understand the interaction between individual and cultural profiles in the context of the moderator.
- Identify significant factors that contribute to variations in CNI.
*Model fit* is assessed using `car::Anova` wherein the significance of the interaction and the overall model is of interest. Type III sum of squares (focus is on the extent to which the moderator interaction explains variability in the model, above and beyond other main and random effects; Type III is invariant to order effects unlike Type I and accounts for interaction effects unlike Type II). The significant contrasts for each model are retained.
*Marginal effects* or *trends* (slopes) were estimated, adjusted for other variables in the model. This is especially useful when dealing with interactions where direct interpretation of coefficients can be challenging. The `emtrends()` function calculates the estimated slopes of `z_ctry_response` at each level of `modtr_score`, adjusting for other variables in the model (here,`overall.M`) and accounting for the random effects. By specifying `pairwise ~ modtr_score`, pairwise comparisons of these slopes between different levels of `modtr_score` are computed., The `adjust = "holm"` argument ensures that the p-values for these pairwise comparisons are adjusted to reduce the chance of Type I errors due to multiple testing (multiple comparisons adjustment).
## Tables - All differences
```{r table_catMods, echo=FALSE, eval=FALSE}
df <- data.frame(
`models` = c(
paste0(c("Mindset", "Isms", "Globe", "Globe (Extra items)", "Fanaticism", "Moral Foundations", "Social Axioms", "Personality 6"), " CNI <- ", c("Sex", "Family Home", "Mother's education", "Father's education", "Age group"))
))
table_catMods <- kable(df, format = "latex", align = c("l", "l")) %>%
kable_styling(full_width = FALSE) %>%
row_spec(0, bold = TRUE) %>%
row_spec(1, italic = TRUE) %>%
add_header_above(c(" " = 2), bold = TRUE)
```
Mindset CNI being the most complex model, it required a different approach for running parallel processing. Here the models are computed first followed by assesing the contrasts between different groups.
```{r mindsetCNI.catMod-models}
cni.CatMod_effects_models <- function(cni_profiles_data, modtr_cat) {
#consists of self,ctry and overall.M profiles
joined_profiles <- pluck(cni_profiles_data, "self.ctry_profiles") %>%
left_join(pluck(cni_profiles_data, "overallM_profiles"), by = "item")
#right data format
cni.modtr_data <- data2 %>%
#select moderator scales individual level (pID) scores
select(pID, modtr_cat) %>%
pivot_longer(-pID,
names_to = "modtr", #moderator name
values_to = "modtr_score") %>% # moderator score
full_join(joined_profiles, by = c("pID"), relationship = "many-to-many")
#CNI <- moderator
cni.modtr_model <- cni.modtr_data %>%
#cni ~ mod [linear]
lmer(z_response ~ z_ctry_response * modtr_score + overall.M + (-1 + z_ctry_response | pID), data = .)
return(cni.modtr_model)
}
```
```{r parallel-processing-mindsetCNI.catMod, echo = F, eval = F}
# Hardware used: (12-1) cores, 62GB RAM
library(furrr);library(purrr);library(tictoc)
library(car)
emm_options(pbkrtest.limit = 521388) #emtrends setting
# Set up cores
plan(multisession, workers = availableCores() - 1)
tic()
# All categorical moderator variables
var_names <- c("sex", "age_group",
"mother_edu","father_edu",
"family_home")
# Parallel execution using furrr
mindsetCNI.CatMod_effects_models<- future_map(var_names, function(var) {
result <- cni.CatMod_effects_models(mindsetCNI_profiles, var)
return(result)
}, .progress = TRUE) %>%
set_names(var_names)
mindsetCNI.CatMod_effects_models$father_edu %>% #object.size()
tidy() %>%
object.size()
# Reset to sequential processing
plan(sequential)
toc()
#object diagnostics:
#library(pryr)
#rm(globeExCNI.cat_effects)
#object_size(fantcCNI.cat_effects)
#object_size(globeExCNI.cat_effects$father_edu)
```
```{r save-mindset.cni.CatMod-models, echo = F, eval = F}
#checklist
#1. change name of the object
#2. DANGER: change name of the object being saved
#use the following to save each or collected chunks in the code above
save(mindsetCNI.CatMod_effects_models,
file = here("objects/Study2/mindsetCNI.CatMod_effects_models.Rdata"))
```
```{r echo = F}
load("~/Desktop/CNI_Dissertation/objects/Study2/mindsetCNI.CatMod_effects_models.Rdata")
```
The focus here is to understand how the relationship between `z_ctry_response` and `z_response` changes across different levels of `modtr_score`. Slopes of z_response with respect to z_ctry_response at different values of modtr_score are computed. `emtrends()` from the `emmeans` package is used to compute estimated marginal means (EMMs) of trends. The p-values for multiple comparisons are adjusted using Holm's method. The limit for the parametric bootstrap Kenward-Roger method is increased here.
```{r mindsetCNI.catMod-trends, message=FALSE, warning=FALSE, eval=FALSE, echo=TRUE}
#emm_options(pbkrtest.limit = 1623160)
cni.CatMod_effects_trends <- \(cni.modtr_model) {
#2. contrasts
trends <- emtrends(cni.modtr_model, pairwise ~ modtr_score,
var = "z_ctry_response",
adjust = "holm")
return(trends)
}
mindsetCNI.CatMod_effects_trends <- map(mindsetCNI.CatMod_effects_models, cni.CatMod_effects_trends)
```
```{r save-mindsetCNI.catMod-trends, eval=FALSE, echo=FALSE}
save(mindsetCNI.CatMod_effects_trends, file = here("objects/Study2/mindsetCNI.CatMod_effects_trends.Rdata"))
```
```{r load-mindsetCNI.catMod-trends, echo = F}
load("~/Desktop/CNI_Dissertation/objects/Study2/mindsetCNI.CatMod_effects_trends.Rdata")
```
```{r fun.-sig_groups1-contrasts, echo = FALSE}
sig_groups1 <- \(df, var, cni_name, sig.only = FALSE) {
# Extract the part of cni_name before the first period
cni_name_cap <- str_extract(cni_name, "^[^.]+")
#age_group -> "Age group"
var_name_cap <- var %>%
str_replace_all("_", " ") %>%
str_to_title()
trends <- df %>%
pluck(var, "contrasts") %>%
as.data.frame()
if (sig.only) {
trends <- filter(trends, p.value < .05)
cap <- paste0("Significant differences in ", cni_name_cap, " by ", var_name_cap, ".")
} else {
cap <- paste0("Differences in ", cni_name_cap , " by ", var_name_cap, ".")
}
if (nrow(trends) > 0) {
trends %>%
mutate(
p.value = papaja::printp(p.value),
across(where(is.numeric), printnum)
) %>%
select(-df) %>%
kable(booktabs = T,
caption = cap,
longtable = T,
format = "latex",
col.names = c("Contrast", "Est", "SE", "Z", "p")) %>%
kable_styling(latex_options = c("repeat_header")) %>%
column_spec(1, width = "7cm")
}
}
```
```{r mindsetCNI.catMod-trends-tables, results='asis'}
#| cache: true
cat_mod_names <- c("family_home", "mother_edu", "father_edu", "age_group") #"sex",
# Print all the tables
walk(cat_mod_names, function(var) {
cat("\n\n")
print(sig_groups1(mindsetCNI.CatMod_effects_trends, var, "Mindset CNI"))
cat("\n\n")
})
```
For degrees of freedom, Kenward-Roger with Satterthwaite's approximation for df was initially preferred as it provides adjusted degrees of freedom and adjusts the covariance matrix of the fixed effects, aiming for more accurate Type I error rates, especially in complex models. When dealing with highly unbalanced data, heteroscedasticity, or models with unequal variances, and multiple random effects, Kenward-Roger provides better variance estimation, improving the accuracy of hypothesis tests.
```{r fun.-cni-CatMod-effects, echo= T, eval=F}
#Explanation for function below
#[which CNI-type]:
#cni_profiles_data <- mindsetCNI_profiles1 [self.ctry_profiles,overallM_profiles]
#[which moderator]:
# "sex","family_home","mother_edu","father_edu","age_group"
#Model: CNI <- Categorical moderator
#levels(factor(data2$father_edu)) #how many levels?
#choose(9, 2) #check for number of contrasts
#This version of the function is necessary for parallel processing! Moderators are varied.
#cni.CatMod_effects(mindsetCNI_profiles, "family_home")
cni.CatMod_effects <- function(cni_profiles_data, modtr_cat) {
#consists of self,ctry and overall.M profiles
joined_profiles <- pluck(cni_profiles_data, "self.ctry_profiles") %>%
left_join(pluck(cni_profiles_data, "overallM_profiles"), by = "item")
#right data format
cni.modtr_data <- data2 %>%
#select moderator scales individual level (pID) scores
select(pID, all_of(modtr_cat)) %>%
pivot_longer(-pID,
names_to = "modtr", #moderator name
values_to = "modtr_score") %>% # moderator score
full_join(joined_profiles, by = c("pID"), relationship = "many-to-many")
#1. CNI <- moderator
cni.modtr_model <- cni.modtr_data %>%
#cni ~ mod [linear]
lmer(z_response ~ z_ctry_response * modtr_score + overall.M + (-1 + z_ctry_response | pID), data = .)
#save results
#2. model fit
#Deafult = Kenward-Roger "F" tests with Satterthwaite degrees of freedom
model_fit <- Anova(cni.modtr_model, type = "III", test.statistic = "F") %>% tidy()
#3. contrasts
trends <- emtrends(cni.modtr_model, pairwise ~ modtr_score,
var = "z_ctry_response",
adjust = "holm")
return(lst(model_fit,trends))
}
```
```{r parallel-processing-catMods, echo = F, eval = F}
#This code chunk is representative of the
# Hardware specifications: 12 cores, 62 GB RAM
# Load necessary libraries
library(furrr);library(purrr);library(tictoc)
library(car)
emm_options(pbkrtest.limit = 521388) #emtrends setting
# Set up cores
plan(multisession, workers = availableCores() - 1)
tic()
# All categorical moderator variables
var_names <- c("sex", "age_group" , "mother_edu","father_edu",
"family_home")
#checklist
#1. change the name of the object according to the CNI being used
#2. change the data_profiles object name
#3. save the object in Study 2
#pty6
# Parallel execution using furrr
pty6CNI.cat_effects<- future_map(var_names, function(var) {
result <- cni.CatMod_effects(pty6CNI_profiles, var)
return(result)
}, .progress = TRUE) %>%
set_names(var_names)
# Reset to sequential processing
plan(sequential)
toc()
#object diagnostics:
#library(pryr)
#rm(globeExCNI.cat_effects)
#object_size(fantcCNI.cat_effects)
#object_size(globeExCNI.cat_effects$father_edu)
```
```{r save-cni.CatMod-models, echo = F, eval = F}
#checklist
#1. change name of the object
#2. DANGER: change name of the object being saved
#use the following to save each or collected chunks in the code above
save(pty6CNI.cat_effects,
file = here("objects/Study2/pty6CNI.cat_effects.Rdata"))
```
```{r load-cni.CatMod-models, echo=F, results='hide'}
#| cache: true
options(scipen = 999)
# Define a vector of file names
cni.catMod_object_names <- paste0(c(
"ismCNI",
"globeCNI",
"globeExCNI",
"fantcCNI",
"mfCNI",
"socAxCNI",
"pty6CNI"), ".cat_effects.Rdata")
# Set file paths and set names
filepaths <- here("objects/Study2", cni.catMod_object_names) %>%
set_names(nm = basename(.) %>% tools::file_path_sans_ext())
# Load all files
all_cni.catMod_effects <- purrr::map(filepaths, ~get(load(.)))
# Assign loaded objects to global environment - important!
purrr::pmap(.l = list(.x = names(all_cni.catMod_effects), .y = all_cni.catMod_effects), .f = ~assign(.x, .y, envir = .GlobalEnv))
# Verify that objects are now in the global environment
#map_lgl(names(all_cni.catMod_effects), ~exists(.x, envir = .GlobalEnv))
all_cni.catMod_effects$ismCNI.cat_effects$age_group$trends$emtrends
```
### Model fit
```{r anova-cni.catMod_effects}
cat_var_names <- c("age_group" , "mother_edu","father_edu",
"family_home")#"sex",
cat_fit_tbl <- \(cat_var){
title_cap <- function(s){
s %>%
str_replace_all("_", " ") %>%
str_to_sentence() %>%
str_replace("edu", "Education")}
map_dfr(all_cni.catMod_effects,
pluck, cat_var, "model_fit",
.id = "CNI_type") %>%
filter(!str_detect(tolower(term), "intercept")) %>% #distinct(term)
mutate(term = case_when(
term == "z_ctry_response" ~ "CNI",
term == "modtr_score" ~ "Predictor",
term == "overall.M" ~ "Overall Mean",
grepl(":", term, fixed = TRUE) ~ "Interaction",
TRUE ~ term
),
stars = case_when(
p.value < .001 ~ "***",
p.value < .01 ~ "**",
p.value < .05 ~ "*",
TRUE ~ ""),
statistic = sprintf("F(%d, %.0f) = %.2f%s",
df, Df.res, statistic, stars),
CNI_type = case_when(
CNI_type == "ismCNI.cat_effects" ~ "Isms CNI",
CNI_type == "globeCNI.cat_effects" ~ "Globe CNI",
CNI_type == "globeExCNI.cat_effects" ~ "Globe Extra CNI",
CNI_type == "mfCNI.cat_effects" ~ "Moral Foundations CNI",
CNI_type == "fantcCNI.cat_effects" ~ "Fanaticism CNI",
CNI_type == "pty6CNI.cat_effects" ~ "Personality 6 CNI",
CNI_type == "socAxCNI.cat_effects" ~ "Social Axioms CNI",
TRUE ~ CNI_type)) %>%
select(CNI_type, term, statistic) %>%
pivot_wider(names_from = term, values_from = statistic) %>%
kable(booktabs = T,
longtable = T,
format = "html",
caption = paste0("ANOVA summary of ", title_cap({{cat_var}}), " individual characteristic models")) %>%
kable_styling() %>%
column_spec(2, width = "3cm") %>%
column_spec(3, width = "3cm") %>%
column_spec(4, width = "3cm") %>%
column_spec(5, width = "3cm")
}
```
```{r walk-anova-cni.catMod_effects, results = 'asis', eval=TRUE, echo=FALSE}
#Define the categorical moderators
cat_mod_names <- c("family_home", "mother_edu", "father_edu", "age_group")
# Use walk to iterate over the categorical moderators
walk(cat_mod_names, function(var) {
result <- cat_fit_tbl(var)
print(result)
cat("\n\n") # Add some space between tables
})
```
For assessing model fit.
Issues with running model fit for Mindset-CNI were encountered. The data was too large for `car::Anova` with `statistic = "F"` (tried with 23 cores and 120 GB RAM on Talapas); code kept crashing. As an alternative, chi square test comparing the H0(null model; no moderator) and H1(alternative model; with moderator) was computed for every moderator. `car::Anova` with `statistic = "Chisq"` is used to facilitate model comparison across all CNI models.
```{r cni.CatMod_AnovaChisq}
cni.CatMod_AnovaChisq <- function(cni_profiles_data, modtr_cat) {
#consists of self,ctry and overall.M profiles
joined_profiles <-
pluck(all_profiles, cni_profiles_data, "self.ctry_profiles") %>%
left_join(pluck(all_profiles , cni_profiles_data, "overallM_profiles"), by = "item")
#right data format
cni.modtr_data <- data2 %>%
#select moderator scales individual level (pID) scores
select(pID, all_of(modtr_cat)) %>%
pivot_longer(-pID,
names_to = "modtr", #moderator name
values_to = "modtr_score") %>% # moderator score
full_join(joined_profiles, by = c("pID"), relationship = "many-to-many")
#1. CNI <- moderator
cni.modtr_model <- cni.modtr_data %>%
#cni ~ mod [linear]
lmer(z_response ~ z_ctry_response*modtr_score + overall.M + (-1 + z_ctry_response | pID), data = .)
#save results
#2. model fit
#Deafult = Kenward-Roger with Satterwhiaite df Approximation
model_fit <- Anova(cni.modtr_model, type = "III", test.statistic = "Chisq") %>% broom::tidy()
return(model_fit)
}
#cni.CatMod_AnovaChisq("pty6CNI_profiles", "sex")
cni.CatMod_AnovaChisq_details <- expand_grid(
cni = names(all_profiles),
modtr = c("sex", "age_group" , "mother_edu","father_edu",
"family_home"))
```
```{r DONT-RUN-cni.CatMod_AnovaChisq, echo = F, eval=FALSE}
#Hardware specifications:
# 16GB RAM, 8 cores
# Load libraries
library(furrr);library(purrr);library(tictoc);library(car); library(broom)
# Set up cores
plan(multisession, workers = availableCores() - 1)
# Run the analysis
tic()
cni.cat_AnovaChisq_fit <- future_pmap(
cni.CatMod_AnovaChisq_details,
function(cni, modtr) {
result <- cni.CatMod_AnovaChisq(cni, modtr)
return(result)
},
.progress = TRUE
) %>%
set_names(with(cni.CatMod_AnovaChisq_details, paste(cni, modtr, sep = "_")))
# Reset to sequential processing
plan(sequential)
toc()
cni.cat_AnovaChisq_fit
```
```{r save-cni.CatMod_AnovaChisq, echo = F, eval = F}
save(cni.cat_AnovaChisq_fit,
file = here("objects/Study2/cni.cat_AnovaChisq_fit.Rdata"))
```
```{r load-cni.CatMod_AnovaChisq, echo = F, eval = T}
load("~/Desktop/CNI_Dissertation/objects/Study2/cni.cat_AnovaChisq_fit.Rdata")
```
```{r cat_fit.chisq_tbl}
cat_fit.chisq_tbl <- function(cat_var) {
title_cap <- function(s) {
s %>%
str_replace_all("_", " ") %>%
str_to_sentence() %>%
str_replace("edu", "Education")
}
bind_rows(cni.cat_AnovaChisq_fit, .id = "CNI_type") %>%
separate(CNI_type, into = c("CNI_type", "profiles", "moderator"), sep = "_", extra = "merge") %>%
filter(moderator == cat_var) %>%
select(-profiles, -moderator)%>%
filter(!str_detect(tolower(term), "intercept")) %>%
mutate(
term = case_when(
term == "z_ctry_response" ~ "CNI",
term == "modtr_score" ~ "Predictor",
term == "overall.M" ~ "Overall Mean",
grepl(":", term, fixed = TRUE) ~ "Interaction",
TRUE ~ term),
stars = case_when(
p.value < .001 ~ "***",
p.value < .01 ~ "**",
p.value < .05 ~ "*",
TRUE ~ ""
),
statistic = sprintf("χ²(%d) = %.2f%s",
df, statistic, stars),
CNI_type = case_when(
CNI_type == "mindsetCNI" ~ "Mindset CNI",
CNI_type == "ismCNI" ~ "Isms-Mindset CNI",
CNI_type == "fantcCNI" ~ "Fanaticism-Mindset CNI",
CNI_type == "globeCNI" ~ "Globe-Mindset CNI",
CNI_type == "globeExCNI" ~ "Globe Extra-Mindset CNI",
CNI_type == "mfCNI" ~ "Moral Foundations-Mindset CNI",
CNI_type == "socAxCNI" ~ "Social Axioms-Mindset CNI",
CNI_type == "pty6CNI" ~ "Personality 6 CNI",
TRUE ~ CNI_type
)) %>%
select(CNI_type, term, statistic) %>%
pivot_wider(names_from = term, values_from = statistic) %>%
kable(booktabs = T,
longtable = T,
format = "html",
caption = paste0("ANOVA summary of ", title_cap({{cat_var}}), " individual characteristic models")) %>%
kable_styling() %>%
column_spec(2, width = "3cm") %>%
column_spec(3, width = "3cm") %>%
column_spec(4, width = "3cm") %>%
column_spec(5, width = "3cm")}
```
```{r walk-cat_fit.chisq_tbl, results = 'asis', eval=TRUE, echo=FALSE}
cat_var_names <- c("age_group" , "mother_edu","father_edu",
"family_home")#"sex",
# Use walk to iterate over the categorical moderators
walk(cat_var_names, function(var) {
result <- cat_fit.chisq_tbl(var)
print(result)
cat("\n\n") # Add some space between tables
})
```
### Marginal effects
#### Estimated slopes (trends)
```{r emtrends-cni.catMod_effects}
cat_emtrends_tbl <- \(cat_var){
title_cap <- function(s){
s %>%
str_replace_all("_", " ") %>%
str_to_sentence() %>%
str_replace("edu", "Education")}
group_names <- c("Mindset CNI", "Isms-Mindset CNI", "Fanaticism-Mindset CNI", "Globe-Mindset CNI", "Globe Extra-Mindset CNI", "Moral Foundations-Mindset CNI", "Social Axioms-Mindset CNI", "Personality6 CNI")
full_df <-
rbind(
mindsetCNI.CatMod_effects_trends %>%
pluck(cat_var, "emtrends") %>%
as_tibble() %>%
mutate(CNI_type = "mindsetCNI.cat_effects",
.before = modtr_score),
map_df(all_cni.catMod_effects,
~as_tibble(summary(pluck(.x, cat_var, "trends", "emtrends"))),
.id = "CNI_type"))
full_df %>%
mutate(CNI_type = case_when(
CNI_type == "mindsetCNI.cat_effects" ~ "Mindset CNI",
CNI_type == "ismCNI.cat_effects" ~ "Isms-Mindset CNI",
CNI_type == "globeCNI.cat_effects" ~ "Globe-Mindset CNI",
CNI_type == "globeExCNI.cat_effects" ~ "Globe Extra-Mindset CNI",
CNI_type == "fantcCNI.cat_effects" ~ "Fanaticism-Mindset CNI",
CNI_type == "mfCNI.cat_effects" ~ "Moral Foundations-Mindset CNI",
CNI_type == "socAxCNI.cat_effects" ~ "Social Axioms-Mindset CNI",
CNI_type == "pty6CNI.cat_effects" ~ "Personality6 CNI",
TRUE ~ CNI_type)) %>%
select(-CNI_type, -df) %>%
mutate(across(where(is.numeric), printnum)) %>%
kable(booktabs = T,
longtable = T,
format = "latex",
col.names = c("Group", "CNI", "SE", "Lower", "Upper"),
caption = paste0("Expected CNI by ", title_cap({{cat_var}}), ".")) %>% #
kable_styling(latex_options = c("repeat_header", "longtable")) %>%
{reduce(seq_along(group_names),
.init = .,
~ group_rows(.x, group_names[.y], n_distinct(full_df$modtr_score)*(.y-1)+1, n_distinct(full_df$modtr_score)*.y))} %>%
column_spec(1, width = "7.2cm") %>%
add_header_above(c(" " = 3, "95% CI" = 2))
}
```
```{r walk-emtrends-cni.catMod_effects, results = 'asis', echo=FALSE, eval=TRUE}
#| cache: true
#Define the categorical moderators
cat_mod_names <- c("family_home", "mother_edu", "father_edu", "age_group")
# Use walk to iterate over the categorical moderators
walk(cat_mod_names, function(var) {
result <- cat_emtrends_tbl(var)
print(result)
cat("\n\n") # Add some space between tables
})
```
```{r emtrends-cni.catMod-emtrends_plot}
#| cache: true
cat_emtrends_plot <- \(cat_var){
xtitle <- function(s){
s %>%
str_replace_all("_", " ") %>%
str_to_sentence() %>%
str_replace("edu", "Education")}
plot_data <- rbind(
mindsetCNI.CatMod_effects_trends %>%
pluck(cat_var, "emtrends") %>%
as_tibble() %>%
mutate(CNI_type = "mindsetCNI.cat_effects",
.before = modtr_score),
map_df(all_cni.catMod_effects,
~as_tibble(summary(pluck(.x, cat_var, "trends", "emtrends"))),
.id = "CNI_type")) %>%
mutate(CNI_type = case_when(
CNI_type == "mindsetCNI.cat_effects" ~ "Mindset CNI",
CNI_type == "ismCNI.cat_effects" ~ "Isms-Mindset CNI",
CNI_type == "globeCNI.cat_effects" ~ "Globe-Mindset CNI",
CNI_type == "globeExCNI.cat_effects" ~ "Globe Extra-Mindset CNI",
CNI_type == "fantcCNI.cat_effects" ~ "Fanaticism-Mindset CNI",
CNI_type == "mfCNI.cat_effects" ~ "Moral Foundations-Mindset CNI",
CNI_type == "socAxCNI.cat_effects" ~ "Social Axioms-Mindset CNI",
CNI_type == "pty6CNI.cat_effects" ~ "Personality6 CNI",
TRUE ~ CNI_type # keep original if no match
))
# Create the plot
ggplot(plot_data, aes(x = modtr_score, y = z_ctry_response.trend, fill = modtr_score)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL),
width = 0.2, position = position_dodge(0.9)) +
scale_y_continuous(limits = c(0, 1), expand = c(0, 0)) +
coord_flip()+
facet_wrap(~ CNI_type, ncol = 2) +
labs(x = xtitle({{cat_var}}),
y = "Estimated Marginal Trend",
fill = "Levels") +
theme_bw() +
theme(text = element_text(family = "serif"),
legend.position = "bottom",
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
# # Control legend text size
# legend.text = element_text(size = 6.4),
# legend.title = element_text(size = 7.6),
# # Control legend key size
# legend.key.size = unit(0.6, "lines"),
# # Control spacing
# legend.spacing.x = unit(0.1, 'cm'),
# legend.spacing.y = unit(0.1, 'cm'),
# # Optional: reduce margin around legend
# legend.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt")
) +
scale_fill_brewer(palette = "Set3")#,
#guide = guide_legend(ncol = 2))
}
#cat_emtrends_plot("mother_edu")
cat_emtrends_plot_details <- tibble(
cat_mod = c("family_home", "mother_edu", "father_edu", "age_group")) %>%
mutate(cat_mod = map_chr(cat_mod, function(x) {
x %>%
str_replace_all("_", " ") %>%
str_to_sentence() %>%
str_replace("edu", "Education")}),
caption = paste0("Expected CNI by ", cat_mod))
```
```{r walk-emtrends-cni.catMod-emtrends_plot, results = 'asis', echo=FALSE, eval=TRUE}
#| cache: true
#| fig-height: 12
#| fig-width: 10
#| out-width: "100%"
#| fig-cap-location: top
#| fig-cap: !expr |
#| cat_emtrends_plot_details %>%
#| pull(caption)
cat_mod_names <- c("family_home", "mother_edu", "father_edu", "age_group")
# Use walk to iterate over the categorical moderators
walk(cat_mod_names, function(var) {
result <- cat_emtrends_plot(var)
print(result)
cat("\n\n") # Add some space between tables
})
```
#### Contrasts
```{r fun.-sig_groups2}
#| cache: true
sig_groups2 <- function(CNI.cat_effects, var, cni_name, sig.only = FALSE) {
# Extract the part of cni_name before the first period
cni_name_cap <- str_extract(cni_name, "^[^.]+")
#age_group -> "Age group"
var_name_cap <- var %>%
str_replace_all("_", " ") %>%
str_to_title()
trends <- CNI.cat_effects %>%
pluck(cni_name,var, "trends", "contrasts") %>%
as.data.frame()
if (sig.only) {
trends <- filter(trends, p.value < .05)
cap <- paste0("Significant differences in ", cni_name_cap, " by ", var_name_cap, ".")
} else {
cap <- paste0("Differences in ", cni_name_cap, " by ", var_name_cap, ".")
}
if (nrow(trends) > 0) {
trends %>%
mutate(
p.value = papaja::printp(p.value),
across(where(is.numeric), printnum)) %>%
select(-df) %>%
kable(booktabs = T,
caption = cap,
longtable = T,
format = "latex",
col.names = c("Contrast", "Est", "SE", "Z", "p")) %>%
kable_styling(latex_options = c("repeat_header")) %>%
column_spec(1, width = "7.2cm")
}
}
```
```{r cni.catMod_effects_tables, results = 'asis', echo=FALSE, eval=TRUE}
#| cache: true
# categorical moderators
cat_mod_names <- c("family_home", "mother_edu", "father_edu", "age_group")
# CNI names (adjust this list according to your needs)
cni_names <- names(all_cni.catMod_effects)
# Iterate over both cni_names and cat_mod_names
walk(cni_names, function(cni) {
walk(cat_mod_names , function(var) {
result <- sig_groups2(all_cni.catMod_effects, var, cni)
if (!is.null(result)) {
print(result)
cat("\n\n") # Add some space between tables
}
})
})
```
## Tables - Contrasts( _p_ < .05)
```{r sig-only.catMod-trends-plots, results='asis', echo=FALSE, eval=TRUE}
#| cache: true
# Moderator names
cat_mod_names <- c("family_home", "mother_edu", "father_edu", "age_group") #"sex",
#--Mindset CNI--#
# Print all the tables
walk(cat_mod_names, function(var) {
cat("\n\n")
print(sig_groups1(mindsetCNI.CatMod_effects_trends, var, "Mindset CNI", sig.only = TRUE))
cat("\n\n")
})
#--others CNI's--#
# CNI names (adjust this list according to your needs)
cni_names <- names(all_cni.catMod_effects)
# Iterate over both cni_names and cat_mod_names
walk(cni_names, function(cni) {
walk(cat_mod_names , function(var) {
result <- sig_groups2(all_cni.catMod_effects, var, cni, sig.only = TRUE)
if (!is.null(result)) {
print(result)
cat("\n\n") # Add some space between tables
}
})
})
```
# Discussion section
## MFQ responding across countries
If global south countries emphasize moral virtues more, compelling more desirable responding in reference to them, (higher SSR) but also get more confused with survey methods and many unfamiliar questions and variables (low CNI), this could be explained perhaps.
Discussion: In contrast SSR demonstrated a more complex pattern. At the country level, it often exhibited negative quadratic relationships with CNI, while at the individual level, the associations were positive but smaller in magnitude.
```{r eval=FALSE, include=TRUE}
global_north <- c(
"Canada", "England", "Germany", "Greece", "Japan",
"Poland", "Spain", "USA", "Singapore", "Taiwan")
data2 %>%
select(country, matches("^scM_MFQ.*\\.comm$")) %>%
group_by(country) %>%
distinct(across(contains(".comm")), .keep_all = TRUE) %>%
ungroup() %>%
rowwise() %>%
mutate(avg_score = mean(c_across(matches("^scM_MFQ.*\\.comm$"))),
region = if_else(country %in% global_north, "Global North", "Global South")) %>%
arrange(desc(avg_score)) %>%
select(country, region, everything()) #%>%
#print(n = Inf) %>% View()
```
## Religiosity responding across countries
```{r eval=FALSE, include=TRUE}
data2 %>%
select(country, scM_DRI_OrgRelgAc.comm, scM_DRI_NonOrgRelgAc.comm) %>%
group_by(country) %>%
distinct(across(contains(".comm")), .keep_all = TRUE) %>%
ungroup() %>%
rowwise() %>%
mutate(avg_score = mean(c_across(matches("^scM_DRI_.*\\.comm$"))),
region = if_else(country %in% global_north, "Global North", "Global South")) %>%
arrange(desc(avg_score)) %>%
select(country, region, everything())
cni.ContMod_effects(pty6CNI_profiles$self.ctry_profiles, "scP_Consc.indv")
cni.ContMod_effects <- function(cni_profiles_data, modtr_pattern) {
# Determine the level based on the modtr_pattern
level <- if(grepl("indv", modtr_pattern)) "Individual" else if(grepl("comm", modtr_pattern)) "Country" else "Unknown"
#right data format
cni.modtr_data <- data2 %>%
#select moderator scales individual level (pID) scores
#---$select moderators$ - specify regex
select(pID, matches(modtr_pattern)) %>%
#mean centering the moderator score cols
#since this standardization is col-wise, running pty7 is not an issue
mutate_if(is.numeric,
~. - mean(., na.rm = TRUE)) %>%
pivot_longer(-pID,
names_to = "modtr", #moderator name
values_to = "modtr_score") %>% # moderator score
full_join(cni_profiles_data, by = c("pID"), relationship = "many-to-many")
#CNI <- moderator
cni.modtr_models <- cni.modtr_data %>%
tidyr::nest(data = -modtr) %>%
mutate(
#cni ~ mod
mod.linear = map(data, ~lmer(z_response ~ z_ctry_response * modtr_score + (-1 + z_ctry_response | pID), data = .)),
#cni ~ mod + mod^2
mod.quadratic = map(data, ~lmer(z_response ~ z_ctry_response * modtr_score + z_ctry_response * I(modtr_score^2) + (-1 + z_ctry_response | pID), data = .)),
#save results
tidy_mod.linear = map(mod.linear, tidy, conf.int = TRUE),
tidy_mod.quadratic = map(mod.quadratic, tidy, conf.int = TRUE)
)
results_linear <- cni.modtr_models %>%
select(modtr, tidy_mod.linear) %>%
unnest(cols = c(tidy_mod.linear)) %>%
filter(term == "z_ctry_response:modtr_score") %>%
select(modtr, term, estimate, conf.low, conf.high, p.value)
results_quadratic <- cni.modtr_models %>%
select(modtr, tidy_mod.quadratic) %>%
unnest(cols = c(tidy_mod.quadratic)) %>%
filter(term == "z_ctry_response:I(modtr_score^2)") %>%
select(modtr, term, estimate, conf.low, conf.high, p.value)
output <- bind_rows(
!!level := results_linear,
!!level := results_quadratic,
.id = "level"
)
return(output)
}
```
```{r DONT-RUN-country-slopes, eval=FALSE, echo=FALSE}
## Post hoc
#computing person level range
cni.ContMod_effects <- function(cni_profiles_data, modtr_pattern) {
# Determine the level based on the modtr_pattern
level <- if(grepl("indv", modtr_pattern)) "Individual" else if(grepl("comm", modtr_pattern)) "Country" else "Unknown"
#right data format
cni.modtr_data <- data2 %>%
#select moderator scales individual level (pID) scores
#---$select moderators$ - specify regex
select(pID, matches(modtr_pattern)) %>%
#mean centering the moderator score cols
#since this standardization is col-wise, running pty7 is not an issue
mutate_if(is.numeric,
~. - mean(., na.rm = TRUE)) %>%
pivot_longer(-pID,
names_to = "modtr", #moderator name
values_to = "modtr_score") %>% # moderator score
full_join(cni_profiles_data, by = c("pID"), relationship = "many-to-many")
#CNI <- moderator
cni.modtr_models <- cni.modtr_data %>%
tidyr::nest(data = -modtr) %>%
mutate(
#cni ~ mod
mod.linear = map(data, ~lmer(z_response ~ z_ctry_response * modtr_score + (-1 + z_ctry_response | pID), data = .)),
#cni ~ mod + mod^2
mod.quadratic = map(data, ~lmer(z_response ~ z_ctry_response * modtr_score + z_ctry_response * I(modtr_score^2) + (-1 + z_ctry_response | pID), data = .)),
#save results
tidy_mod.linear = map(mod.linear, tidy, conf.int = TRUE),
tidy_mod.quadratic = map(mod.quadratic, tidy, conf.int = TRUE)
)
results_linear <- cni.modtr_models %>%
select(modtr, tidy_mod.linear) %>%
unnest(cols = c(tidy_mod.linear)) %>%
filter(term == "z_ctry_response:modtr_score") %>%
select(modtr, term, estimate, conf.low, conf.high, p.value)
results_quadratic <- cni.modtr_models %>%
select(modtr, tidy_mod.quadratic) %>%
unnest(cols = c(tidy_mod.quadratic)) %>%
filter(term == "z_ctry_response:I(modtr_score^2)") %>%
select(modtr, term, estimate, conf.low, conf.high, p.value)
output <- bind_rows(
!!level := results_linear,
!!level := results_quadratic,
.id = "level"
)
return(output)
}
cni.modtr_data <- data2 %>% #names()
#select moderator scales individual level (pID) scores
#---$select moderators$ - specify regex
select(pID, #"scP_Big2_SocialSelfReg.indv",
"scP_Big2_SocialSelfReg.comm"#, #"scP_Big2_Dynamism.indv",
#"scP_Big2_Dynamism.comm" ) %>%
#mean centering the moderator score cols
#since this standardization is col-wise, running pty7 is not an issue
# mutate_if(is.numeric,
# ~. - mean(., na.rm = TRUE)) %>%
pivot_longer(-pID,
names_to = "modtr", #moderator name
values_to = "modtr_score") %>% # moderator score
full_join(pty6CNI_profiles$self.ctry_profiles, by = c("pID"), relationship = "many-to-many")
pty6_SSR_ctry <- lmer(z_response ~ z_ctry_response * modtr_score + z_ctry_response * I(modtr_score^2) +
(-1 + z_ctry_response:I(modtr_score^2) | country)+
(-1 + z_ctry_response | pID),
data = cni.modtr_data)
pty6_Dynamism_ctry <- lmer(z_response ~ z_ctry_response * modtr_score + z_ctry_response * I(modtr_score^2) +
(-1 + z_ctry_response:I(modtr_score^2) | country),
data = cni.modtr_data)
coef(pty6_SSR_ctry)$country %>%
as.data.frame() %>%
rownames_to_column(var = "country") %>%
#tibble() %>%
select(country , `z_ctry_response:I(modtr_score^2)`) %>%
mutate(
category = case_when(
country %in% c(
"Tanzania", "Kenya", "Ethiopia", "Morocco", "Nepal",
"India", "Bangladesh", "Malaysia", "Philippines", "Thailand",
"China", "Turkey", "Peru", "Argentina", "Brazil", "Mexico"
) ~ "Global South",
country %in% c(
"USA", "Canada", "England", "Germany", "Spain", "Japan",
"Singapore", "Taiwan", "Greece", "Poland", "Ukraine"
) ~ "Global North",
TRUE ~ "Uncategorized"
),
.before = everything()
) %>%
arrange(desc(abs(`z_ctry_response:I(modtr_score^2)`)))
all_cni_models <- lst(mindsetCNI_mod, pty6CNI_mod, ismCNI_mod, fantcCNI_mod, globeCNI_mod, globeExCNI_mod, mfCNI_mod, socAxCNI_mod)
compute_ctry_slopes <- function(model_obj) {
obj_name <- deparse(substitute(model_obj))
# Calculate slopes
mod_slopes <- slopes(model_obj, variables = "z_ctry_response", by = "pID")
return(mod_slopes)
}
all_slopes <- map_dfr(all_cni_models, compute_slopes, .id = "model_name")
```