# loading packages library(haven) library(tidyverse) library(emmeans) # for estimated marginal means library(jtools) # for summary model with confidence intervals library(semTools) # for calculating Cronbach alpha library(ggpubr) # # for stat_pvalue_manual library(rstatix) # for pairwise comparisons library(ggh4x) # for nested facets # I. Loading data --------------------------------------------------------- d <- read_sav("data/democratic_inoculation_data_cleaned.sav") # II. Descriptive statistics ------------------------------------------------ attr(d$Gender, "labels") round((prop.table(table(d$Gender)) * 100), 2) # male: 48.77%, female: 51.23% sum(!is.na(d$Gender)) # 1175 valid responses psych::describe(d$Age) # M = 43.12, SD = 12.66, Median = 43, Min = 18, Max = 69 attr(d$Education, "labels") round((prop.table(table(d$Education)) * 100), 2) # elementary/vocational: 18.30%, secondary: 48.85%, tertiary: 32.85% sum(!is.na(d$Education)) # 1175 valid responses attr(d$Settlement, "labels") round((prop.table(table(d$Settlement)) * 100), 2) # capital: 24.00%, county seat: 26.30%, city: 31.66%, village: 18.04%, sum(!is.na(d$Settlement)) # 1175 valid responses # party preference attr(d$Partisanship, "labels") table(d$Partisanship) round((prop.table(table(d$Partisanship)) * 100), 2) # kormánypárti: 447 (38.04%) # ellenzéki: 728 (61.96%) d %>% filter(cond != "other") %>% group_by(cond) %>% summarise(n = n()) %>% mutate(prop = n / sum(n), perc = sprintf("%.2f", prop * 100)) # control: 388 (33.02%), # treatment cracks ("reped"): 342 (29.11%), # treatment build-up ("kiep"): 445 (37.87%) # III. Manipulation checks ------------------------------------------------- # gender table(d$cond, d$Gender, useNA = "always") model <- glm(Gender ~ cond, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p = 0.725 # cond2 p = 0.187 means <- emmeans(model, specs = "cond") means pairs(means, reverse = TRUE, adjust = "none") # ps > 0.18 # age model <- lm(Age ~ cond, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p = 0.714 # cond2 p = 0.047 summary(d$Age) means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 43.9 0.642 1172 42.6 45.1 # 1 43.5 0.684 1172 42.2 44.9 # 2 42.1 0.559 1172 41.0 43.3 pairs(means, reverse = TRUE, adjust = "none") # contrast estimate SE df t.ratio p.value # cond1 - cond0 -0.344 0.938 1172 -0.367 0.7138 # cond2 - cond0 -1.744 0.878 1172 -1.985 0.0474 # cond2 - cond1 -1.400 0.909 1172 -1.539 0.1241 # education table(d$cond, d$Education, useNA = "always") model <- glm(Education ~ cond, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p = 0.285 # cond2 p = 0.002 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 2.22 0.0354 1172 2.15 2.29 # 1 2.16 0.0377 1172 2.09 2.24 # 2 2.07 0.0331 1172 2.00 2.13 pairs(means, reverse = TRUE, adjust = "none") # contrast estimate SE df t.ratio p.value # cond1 - cond0 -0.0553 0.0518 1172 -1.069 0.2855 # cond2 - cond0 -0.1517 0.0485 1172 -3.128 0.0018 # cond2 - cond1 -0.0963 0.0502 1172 -1.919 0.0553 # settlement table(d$cond, d$Settlement, useNA = "always") model <- glm(Settlement ~ cond, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p = 0.157 # cond2 p = 0.899 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 2.40 0.0529 1172 2.30 2.51 # 1 2.51 0.0564 1172 2.40 2.62 # 2 2.41 0.0494 1172 2.31 2.51 pairs(means, reverse = TRUE, adjust = "none") # ps > 0.15 # partisanship table(d$Partisanship) prop.table(table(d$Partisanship)) table(d$cond, d$Partisanship, useNA = "always") model <- glm(Partisanship ~ cond, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p = 0.265 # cond2 p = 0.708 means <- emmeans(model, specs = "cond") means pairs(means, reverse = TRUE, adjust = "none") # ps > 0.26 # IV. Analyses ------------------------------------------------------------- ## IV/1. Creating necessary variables ---- ### 1.1. Authoritarian Cue Recognition Index ---- # agreement (“To what extent do you agree with this statement?”), # perceived authoritarian tone (“To what extent does this statement reflect authoritarian patterns?”), # concern (“How concerning would it be if a country’s political leaders regularly used such rhetoric?”), # moral outrage (“To what extent does this statement evoke moral outrage in you?”), and # perceived threat (“When a leading politician says something like this, how threatening do you find it to the presence of an authoritarian political system?”). # Responses are given on 6-point Likert-type scales. d <- d %>% mutate(agreement = rowMeans(dplyr::select(., matches("^out\\d+a$")), na.rm = TRUE), authoritarian_tone = rowMeans(dplyr::select(., matches("^out\\d+c$")), na.rm = TRUE), concern = rowMeans(dplyr::select(., matches("^out\\d+e$")), na.rm = TRUE), moral_outrage = rowMeans(dplyr::select(., matches("^out\\d+b$")), na.rm = TRUE), perceived_threat = rowMeans(dplyr::select(., matches("^out\\d+d$")), na.rm = TRUE)) psych::describe(d$agreement) psych::describe(d$authoritarian_tone) psych::describe(d$concern) psych::describe(d$moral_outrage) psych::describe(d$perceived_threat) # Internal consistency (alpha) calculation agree <- ' agreem =~ out2a + out6a + out8a + out9a + out11a' fit_agree <- cfa(agree, data = d, estimator = 'MLR') agree_rel <- as.data.frame(compRelSEM(fit_agree)) round(agree_rel[1,], 2) # 0.78 auth_tone <- ' authtone =~ out2c + out6c + out8c + out9c + out11c' fit_auth_tone <- cfa(auth_tone, data = d, estimator = 'MLR') auth_tone_rel <- as.data.frame(compRelSEM(fit_auth_tone)) round(auth_tone_rel[1,], 2) # 0.86 concern <- ' concerned =~ out2e + out6e + out8e + out9e + out11e' fit_concern <- cfa(concern, data = d, estimator = 'MLR') concern_rel <- as.data.frame(compRelSEM(fit_concern)) round(concern_rel[1,], 2) # 0.86 moral_outrage <- ' moralout =~ out2b + out6b + out8b + out9b + out11b' fit_moral_outrage <- cfa(moral_outrage, data = d, estimator = 'MLR') moral_outrage_rel <- as.data.frame(compRelSEM(fit_moral_outrage)) round(moral_outrage_rel[1,], 2) # 0.80 perceived_threat <- ' percthreat =~ out2d + out6d + out8d + out9d + out11d' fit_perceived_threat <- cfa(perceived_threat, data = d, estimator = 'MLR') perceived_threat_rel <- as.data.frame(compRelSEM(fit_perceived_threat)) round(perceived_threat_rel[1,], 2) # 0.85 ### 1.2. Cracks Recognition Index (primary) ---- # Ratings of regime instability indicators (e.g., elite defection, protests, # propaganda losing credibility, inconsistent government communication). # Items will be averaged into a single score. d <- d %>% mutate(regime_instability = rowMeans(dplyr::select(., matches("^crack[1-7]$")), na.rm = TRUE), regime_distractors = rowMeans(dplyr::select(., matches("^crack([8-9]|1[0-9]|20)$")), na.rm = TRUE)) psych::describe(d$regime_instability) psych::describe(d$regime_distractors) # Internal consistency (alpha) calculation regins <- ' reginsta =~ crack1 + crack2 + crack3 + crack4 + crack5 + crack6 + crack7' fit_regins <- cfa(regins, data = d, estimator = 'MLR') regins_rel <- as.data.frame(compRelSEM(fit_regins)) round(regins_rel[1,], 2) # 0.88 names(d) ### 1.3. Government Negativity Index ---- # Average of anger, disgust, contempt, shame, and indignation toward government actors. d$lawview <- as.numeric(d$lawview) d <- d %>% mutate(government_negativity = rowMeans(across(c(govern_anger, govern_disgust, govern_contempt, govern_shame, govern_indignation)), na.rm = TRUE)) psych::describe(d$government_negativity) # Internal consistency (alpha) calculation govnega <- ' govnegat =~ govern_anger + govern_disgust + govern_contempt + govern_shame + govern_indignation' fit_govnega <- cfa(govnega, data = d, estimator = 'MLR') govnega_rel <- as.data.frame(compRelSEM(fit_govnega)) round(govnega_rel[1,], 2) # 0.96 ### 1.4. Civil Society/Media Positivity Index ---- # Average of respect, pride, and empathy toward civil society and media actors. # Items will be averaged and checked for internal consistency. Individual emotion scores will also be reported in exploratory analyses. d <- d %>% mutate(civil_positivity = rowMeans(across(c(target_respect, target_pride, target_empathy)), na.rm = TRUE)) psych::describe(d$civil_positivity) # Internal consistency (alpha) calculation civposi<- ' civilposit =~ target_respect + target_pride + target_empathy' fit_civposi <- cfa(civposi, data = d, estimator = 'MLR') civposi_rel <- as.data.frame(compRelSEM(fit_civposi)) round(civposi_rel[1,], 2) # 0.91 # Competition Attitude Indices (control-domain check). # - Self-developmental competition: average of items highlighting growth-oriented aspects of competition. # # - Hypercompetitive orientation: average of items highlighting a “winning at any cost” mindset (reverse-coded where needed). # # Engagement Index (secondary). Combines (a) the number of words in the open-ended response (log-transformed if needed) and (b) the willingness-to-share item. Each will be standardized and then averaged. # # Treatment-on-the-treated compliance proxy (for sensitivity). Based on time spent on intervention pages, completion of reflection prompt, and attention check performance. These indicators will be standardized and combined to identify an “engaged” subsample for sensitivity analyses. ## IV/2. Analysis based on the Analysis Plan ------------------------------------------------ ### 2.1. Primary outcomes and models ---- #### Authoritarian Cue Recognition Index ---- # OLS with the two condition dummies (no covariates). # agreement model <- lm(agreement ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p < 0.001 # cond2 p = 0.140 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 2.83 0.0514 1170 2.72 2.93 # 1 2.53 0.0546 1170 2.42 2.64 # 2 2.72 0.0480 1170 2.63 2.82 pairs(means, reverse = TRUE, adjust = "bonferroni") # contrast estimate SE df t.ratio p.value # cond1 - cond0 -0.295 0.0749 1170 -3.935 0.0003 # cond2 - cond0 -0.104 0.0705 1170 -1.479 0.4185 # cond2 - cond1 0.190 0.0728 1170 2.617 0.0269 model <- lm(scale(agreement) ~ cond + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") pairs(means, reverse = TRUE, adjust = "bonferroni") # authoritarian_tone model <- lm(authoritarian_tone ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p < 0.001 # cond2 p = 0.081 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 3.89 0.0603 1170 3.77 4.00 # 1 4.21 0.0641 1170 4.08 4.33 # 2 4.03 0.0564 1170 3.92 4.14 pairs(means, reverse = TRUE, adjust = "bonferroni") # contrast estimate SE df t.ratio p.value # cond1 - cond0 0.320 0.0880 1170 3.639 0.0009 # cond2 - cond0 0.145 0.0828 1170 1.746 0.2431 # cond2 - cond1 -0.175 0.0855 1170 -2.054 0.1207 model <- lm(scale(authoritarian_tone) ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") pairs(means, reverse = TRUE, adjust = "bonferroni") # concern model <- lm(concern ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p < 0.001 # cond2 p = 0.187 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 4.08 0.0592 1170 3.96 4.20 # 1 4.42 0.0629 1170 4.29 4.54 # 2 4.19 0.0553 1170 4.08 4.30 pairs(means, reverse = TRUE, adjust = "bonferroni") # contrast estimate SE df t.ratio p.value # cond1 - cond0 0.337 0.0863 1170 3.901 0.0003 # cond2 - cond0 0.112 0.0812 1170 1.380 0.5034 # cond2 - cond1 -0.224 0.0838 1170 -2.678 0.0225 model <- lm(scale(concern) ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # SMD=0.286 means <- emmeans(model, specs = "cond") pairs(means, reverse = TRUE, adjust = "bonferroni") # moral_outrage model <- lm(moral_outrage ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p < 0.001 # cond2 p = 0.017 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 3.36 0.0604 1170 3.24 3.48 # 1 3.75 0.0642 1170 3.62 3.87 # 2 3.56 0.0565 1170 3.45 3.67 pairs(means, reverse = TRUE, adjust = "bonferroni") # contrast estimate SE df t.ratio p.value # cond1 - cond0 0.385 0.0881 1170 4.376 <.0001 # cond2 - cond0 0.199 0.0829 1170 2.400 0.0497 # cond2 - cond1 -0.186 0.0856 1170 -2.179 0.0886 model <- lm(scale(moral_outrage) ~ cond + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") pairs(means, reverse = TRUE, adjust = "bonferroni") # SMDcracks=0.321 # SMDbuiltup=0.166 # perceived_threat model <- lm(perceived_threat ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p < 0.001 # cond2 p = 0.064 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 3.97 0.0570 1170 3.86 4.09 # 1 4.32 0.0605 1170 4.20 4.43 # 2 4.12 0.0533 1170 4.01 4.22 pairs(means, reverse = TRUE, adjust = "bonferroni") # contrast estimate SE df t.ratio p.value # cond1 - cond0 0.341 0.0831 1170 4.106 0.0001 # cond2 - cond0 0.145 0.0782 1170 1.857 0.1908 # cond2 - cond1 -0.196 0.0807 1170 -2.426 0.0462 model <- lm(scale(perceived_threat) ~ cond + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") pairs(means, reverse = TRUE, adjust = "bonferroni") #SMD=0.300 #### Cracks Recognition Index ---- # same OLS specification model <- lm(regime_instability ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p < 0.001 # cond2 p = 0.706 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 4.43 0.0469 1170 4.34 4.52 # 1 4.70 0.0498 1170 4.60 4.80 # 2 4.46 0.0438 1170 4.37 4.54 pairs(means, reverse = TRUE, adjust = "bonferroni") # contrast estimate SE df t.ratio p.value # cond1 - cond0 0.2649 0.0684 1170 3.874 0.0003 # cond2 - cond0 0.0243 0.0644 1170 0.378 1.0000 # cond2 - cond1 -0.2405 0.0664 1170 -3.621 0.0009 model <- lm(scale(regime_instability) ~ cond + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") pairs(means, reverse = TRUE, adjust = "bonferroni") # SMD=0.284 ### 2.2. Secondary outcomes ---- #### Agreement with the government's initiative ---- model <- lm(lawview ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p = 0.080 # cond2 p = 0.685 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 3.85 0.140 1170 3.58 4.13 # 1 3.49 0.149 1170 3.20 3.78 # 2 3.77 0.131 1170 3.52 4.03 pairs(means, reverse = TRUE, adjust = "bonferroni") # contrast estimate SE df t.ratio p.value # cond1 - cond0 -0.358 0.204 1170 -1.752 0.2404 # cond2 - cond0 -0.078 0.193 1170 -0.405 1.0000 # cond2 - cond1 0.280 0.199 1170 1.410 0.4763 model <- lm(scale(lawview) ~ cond + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") pairs(means, reverse = TRUE, adjust = "bonferroni") #### Policy Attitude Index and Affective Response indices ---- # identical OLS specification (two condition dummies vs Control) ##### Government Negativity Index ---- model <- lm(government_negativity ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p < 0.001 # cond2 p = 0.011 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 3.80 0.1050 1170 3.59 4.01 # 1 4.35 0.1120 1170 4.13 4.57 # 2 4.17 0.0985 1170 3.98 4.36 pairs(means, reverse = TRUE, adjust = "bonferroni") # contrast estimate SE df t.ratio p.value # cond1 - cond0 0.554 0.154 1170 3.604 0.0010 # cond2 - cond0 0.369 0.145 1170 2.551 0.0326 # cond2 - cond1 -0.185 0.149 1170 -1.237 0.6487 model <- lm(scale(government_negativity) ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") pairs(means, reverse = TRUE, adjust = "bonferroni") # SMDc=0.265 # SMDb=0.176 d$party <- as.factor(d$Partisanship) # # model <- lm(scale(government_negativity) ~ cond + Age + Education + party, data = d) # # summary model with confidence intervals # summ(model, confint = TRUE, ci.width = .95, digits = 3) # # SMDc=0.210 (p<0.001) # # SMDb=0.152 (p=0.006) ##### Civil Society/Media Positivity Index ---- # identical OLS specification (two condition dummies vs Control). model <- lm(civil_positivity ~ cond + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1 p = 0.040 # cond2 p = 0.055 means <- emmeans(model, specs = "cond") means # cond emmean SE df lower.CL upper.CL # 0 3.19 0.0968 1170 3.00 3.38 # 1 3.48 0.1030 1170 3.28 3.68 # 2 3.44 0.0904 1170 3.27 3.62 pairs(means, reverse = TRUE, adjust = "bonferroni") # contrast estimate SE df t.ratio p.value # cond1 - cond0 0.2897 0.141 1170 2.054 0.1207 # cond2 - cond0 0.2551 0.133 1170 1.920 0.1652 # cond2 - cond1 -0.0346 0.137 1170 -0.253 1.0000 model <- lm(scale(civil_positivity) ~ cond + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") pairs(means, reverse = TRUE, adjust = "bonferroni") # SMDc=0.151 (p=0.040) # SMDb=0.133 (p=0.055) # model <- lm(scale(civil_positivity) ~ cond + Age + Education + party, data = d) # summ(model, confint = TRUE, ci.width = .95, digits = 3) # # SMDc=0.115 (p=0.090) # # SMDb=0.117 (p=0.066) ##### Specific emotions toward government actors ---- model <- lm(govern_anger ~ cond*party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(govern_anger ~ cond*scale(populism) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(govern_anger ~ cond*scale(crt_total) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(govern_anger ~ party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "party") means model <- lm(scale(govern_contempt) ~ cond*party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(govern_contempt ~ cond*scale(populism) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(govern_contempt ~ cond*scale(crt_total) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(govern_contempt ~ party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "party") means model <- lm(scale(govern_pride) ~ cond*party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(govern_pride) ~ cond*scale(populism) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(govern_pride) ~ cond*scale(crt_total) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(govern_pride ~ party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "party") means model <- lm(scale(govern_respect) ~ cond*party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(govern_respect) ~ cond*scale(populism) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(govern_respect) ~ cond*scale(crt_total) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(govern_respect ~ party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "party") means model <- lm(scale(govern_shame) ~ cond*party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(govern_shame) ~ cond*scale(populism) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(govern_shame) ~ cond*scale(crt_total) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(govern_shame ~ party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "party") means ##### Specific emotions toward civil society/media actors ---- model <- lm(target_anger ~ cond*party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(target_anger ~ cond*scale(populism) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(target_anger ~ cond*scale(crt_total) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(target_anger ~ party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "party") means model <- lm(scale(target_contempt) ~ cond*party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means #sig a cond2 es a party kozt model <- lm(target_contempt ~ cond*scale(populism) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(target_contempt ~ cond*scale(crt_total) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(target_contempt ~ party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "party") means model <- lm(scale(target_pride) ~ cond*party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(target_pride) ~ cond*scale(populism) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(target_pride) ~ cond*scale(crt_total) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(target_pride ~ party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "party") means model <- lm(scale(target_respect) ~ cond*party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(target_respect) ~ cond*scale(populism) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(target_respect) ~ cond*scale(crt_total) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(target_respect ~ party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "party") means model <- lm(scale(target_shame) ~ cond*party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means #itt is van egy sig interakcio a cond 2 vel model <- lm(scale(target_shame) ~ cond*scale(populism) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(scale(target_shame) ~ cond*scale(crt_total) + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "cond") means model <- lm(target_shame ~ party + Age + Education, data = d) summ(model, confint = TRUE, ci.width = .95, digits = 3) means <- emmeans(model, specs = "party") means # V. Exploratory analyses -------------------------------------------------- # - Political orientation: Self-placement on a 7-point liberal–conservative scale. # # - Partisanship: Voting preference (in terms of pro vs. anti-governmental status) and party closeness (specific party and level of sympathy). # # - Further covariates: Analytical thinking, fact relativism, authoritarian illiberalism, populism, and conspiracy belief. (As exploratory analyses, we will examine whether they will function as moderator variables). # Moderation. # We test political orientation and populism (continuous, mean-centered) as moderators on the two primary outcomes by adding the moderator main effect and its interactions with each condition dummy. We probe significant interactions with simple-slope estimates at conventional values of the moderator (e.g., −1 SD, mean, +1 SD) and by plotting predicted values. Effect sizes are reported as standardized simple-slope differences (and as Δd across moderator values). ## 1. Partisanship ---- # Partisanship: Voting preference (in terms of pro vs. anti-governmental status) table(d$Partisanship) attr(d$Partisanship, "labels") # agreement model <- lm(agreement ~ cond*Partisanship + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:Partisanship p = 0.426 # cond2:Partisanship p = 0.074 means <- emmeans(model, specs = pairwise ~ cond | Partisanship) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # authoritarian_tone model <- lm(authoritarian_tone ~ cond*Partisanship + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:Partisanship p = 0.720 # cond2:Partisanship p = 0.297 means <- emmeans(model, specs = pairwise ~ cond | Partisanship) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # concern model <- lm(concern ~ cond*Partisanship + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:Partisanship p = 0.881 # cond2:Partisanship p = 0.378 means <- emmeans(model, specs = pairwise ~ cond | Partisanship) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # moral_outrage model <- lm(moral_outrage ~ cond*Partisanship + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:Partisanship p = 0.871 # cond2:Partisanship p = 0.662 means <- emmeans(model, specs = pairwise ~ cond | Partisanship) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # perceived_threat model <- lm(perceived_threat ~ cond*Partisanship + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:Partisanship p = 0.577 # cond2:Partisanship p = 0.312 means <- emmeans(model, specs = pairwise ~ cond | Partisanship) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # - Cracks Recognition Index: same OLS specification. model <- lm(regime_instability ~ cond*Partisanship + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:Partisanship p = 0.205 # cond2:Partisanship p = 0.051 means <- emmeans(model, specs = pairwise ~ cond | Partisanship) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # lawview - agreement with the government's initiative model <- lm(lawview ~ cond*Partisanship + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:Partisanship p = 0.896 # cond2:Partisanship p = 0.027 means <- emmeans(model, specs = pairwise ~ cond | Partisanship) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # government_negativity model <- lm(government_negativity ~ cond*Partisanship + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:Partisanship p = 0.405 # cond2:Partisanship p = 0.326 means <- emmeans(model, specs = pairwise ~ cond | Partisanship) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # civil_positivity model <- lm(civil_positivity ~ cond*Partisanship, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:Partisanship p = 0.705 # cond2:Partisanship p = 0.194 means <- emmeans(model, specs = pairwise ~ cond | Partisanship) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") ## 2. Populism ---- # be careful we overwrite the original variables! d <- d %>% # reverse the selected items (1–6 → 6–1) mutate(across(c(pop_1_1, pop_1_3, pop_1_5, pop_1_7), ~ 7 - as.numeric(.))) %>% # compute mean of pop_1_1 – pop_1_12 mutate(populism = rowMeans(dplyr::select(., matches("^pop_1_([1-9]|1[0-2])$")),na.rm = TRUE)) # d <- d %>% # mutate(populism = rowMeans(dplyr::select(., matches("^pop_1_([1-9]|1[0-2])$")), na.rm = TRUE)) psych::describe(d$populism) # Internal consistency (alpha) calculation popul <- ' populi =~ pop_1_1 + pop_1_2 + pop_1_3 + pop_1_4 + pop_1_5 + pop_1_6 + pop_1_7 + pop_1_8 + pop_1_9 + pop_1_10 + pop_1_11 + pop_1_12' fit_popul <- cfa(popul, data = d, estimator = 'MLR') popul_rel <- as.data.frame(compRelSEM(fit_popul)) round(popul_rel[1,], 2) # 0.70 # agreement model <- lm(agreement ~ cond*scale(populism) + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:populism p = 0.042 # cond2:populism p = 0.053 means <- emmeans(model, specs = pairwise ~ cond | populism) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # authoritarian_tone model <- lm(authoritarian_tone ~ cond*populism + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:populism p = 0.577 # cond2:populism p = 0.678 means <- emmeans(model, specs = pairwise ~ cond | populism) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # concern model <- lm(concern ~ cond*populism + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:populism p = 0.925 # cond2:populism p = 0.989 means <- emmeans(model, specs = pairwise ~ cond | populism) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # moral_outrage model <- lm(moral_outrage ~ cond*populism + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:populism p = 0.528 # cond2:populism p = 0.407 means <- emmeans(model, specs = pairwise ~ cond | populism) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # perceived_threat model <- lm(perceived_threat ~ cond*populism + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:populism p = 0.632 # cond2:populism p = 0.817 means <- emmeans(model, specs = pairwise ~ cond | populism) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # - Cracks Recognition Index: same OLS specification. model <- lm(regime_instability ~ cond*populism + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:populism p = 0.270 # cond2:populism p = 0.863 means <- emmeans(model, specs = pairwise ~ cond | populism) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # lawview - agreement with the government's initiative model <- lm(lawview ~ cond*populism + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:populism p = 0.207 # cond2:populism p = 0.126 means <- emmeans(model, specs = pairwise ~ cond | populism) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # government_negativity model <- lm(government_negativity ~ cond*populism + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:populism p = 0.290 # cond2:populism p = 0.875 means <- emmeans(model, specs = pairwise ~ cond | populism) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # civil_positivity model <- lm(civil_positivity ~ cond*populism + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:populism p = 0.212 # cond2:populism p = 0.524 means <- emmeans(model, specs = pairwise ~ cond | populism) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") ## 3. CRT ---- table(d$crt1) # 5000 table(d$crt2) # 5 table(d$crt3) # 47 d <- d %>% mutate(crt1_out = if_else(crt1 == 5000, 1, 0), crt2_out = if_else(crt2 == 5, 1, 0), crt3_out = if_else(crt3 == 47, 1, 0)) d <- d %>% mutate(crt_total = crt1_out + crt2_out + crt3_out) table(d$crt_total) # agreement model <- lm(agreement ~ cond*scale(crt_total) + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:crt_total p = 0.210 # cond2:crt_total p = 0.881 means <- emmeans(model, specs = pairwise ~ cond | crt_total) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # authoritarian_tone model <- lm(authoritarian_tone ~ cond*crt_total + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:crt_total p = 0.748 # cond2:crt_total p = 0.192 means <- emmeans(model, specs = pairwise ~ cond | crt_total) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # concern model <- lm(concern ~ cond*crt_total + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:crt_total p = 0.515 # cond2:crt_total p = 0.708 means <- emmeans(model, specs = pairwise ~ cond | crt_total) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # moral_outrage model <- lm(moral_outrage ~ cond*crt_total + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:crt_total p = 0.798 # cond2:crt_total p = 0.645 means <- emmeans(model, specs = pairwise ~ cond | crt_total) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # perceived_threat model <- lm(perceived_threat ~ cond*crt_total + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:crt_total p = 0.245 # cond2:crt_total p = 0.868 means <- emmeans(model, specs = pairwise ~ cond | crt_total) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # - Cracks Recognition Index: same OLS specification. model <- lm(regime_instability ~ cond*crt_total + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:crt_total p = 0.051 # cond2:crt_total p = 0.735 means <- emmeans(model, specs = pairwise ~ cond | crt_total) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # lawview - agreement with the government's initiative model <- lm(lawview ~ cond*crt_total + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:crt_total p = 0.299 # cond2:crt_total p = 0.461 means <- emmeans(model, specs = pairwise ~ cond | crt_total) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # government_negativity model <- lm(government_negativity ~ cond*crt_total + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:crt_total p = 0.143 # cond2:crt_total p = 0.105 means <- emmeans(model, specs = pairwise ~ cond | crt_total) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # civil_positivity model <- lm(civil_positivity ~ cond*crt_total + Age + Education, data = d) # summary model with confidence intervals summ(model, confint = TRUE, ci.width = .95, digits = 3) # cond1:crt_total p = 0.468 # cond2:crt_total p = 0.382 means <- emmeans(model, specs = pairwise ~ cond | crt_total) means$emmeans pairs(means$emmeans, reverse = TRUE, adjust = "none") # VI/1. Visualization of all outcome variables using zoomed y-axis scales ---- # Recode condition variable for plotting d <- d %>% mutate(cond = factor(cond, levels = c(0, 1, 2), labels = c("Control", "Cracks", "Gradual Build-Up"))) # Define outcomes and their scale ranges psych::describe(d$agreement) psych::describe(d$authoritarian_tone) psych::describe(d$concern) psych::describe(d$moral_outrage) psych::describe(d$perceived_threat) psych::describe(d$regime_instability) psych::describe(d$lawview) psych::describe(d$government_negativity) psych::describe(d$civil_positivity) outcomes <- tibble::tribble( ~var, ~outcome, ~ymin, ~ymax, "agreement", "Agreement", 1, 6, "authoritarian_tone", "Perceived authoritarian tone", 1, 6, "concern", "Concern", 1, 6, "moral_outrage", "Moral outrage", 1, 6, "perceived_threat", "Perceived threat", 1, 6, "regime_instability", "Regime instability", 1, 6, "lawview", "Agreement with proposed law", 1, 10, "government_negativity","Government negativity", 1, 7, "civil_positivity", "Civil-society positivity", 1, 7 ) # Calculate means, SDs, and ±2 SD ranges for each outcome # y_limits <- outcomes %>% # rowwise() %>% # mutate( # mean = mean(d[[var]], na.rm = TRUE), # sd = sd(d[[var]], na.rm = TRUE), # y_lower = max(ymin, mean - 2 * sd), # y_upper = min(ymax, mean + 2 * sd) # ) %>% # ungroup() # # y_limits %>% # select(var, ymin, ymax, mean, sd, y_lower, y_upper) # # # pure_sd_ranges <- bind_rows(lapply(outcomes$var, function(v) { # x <- d[[v]] # mu <- mean(x, na.rm = TRUE) # sdx <- sd(x, na.rm = TRUE) # # tibble( # outcome = v, # n = sum(!is.na(x)), # mean = mu, # sd = sdx, # lower_2sd = mu - 2 * sdx, # upper_2sd = mu + 2 * sdx # # lower_2_5sd = mu - 2.5 * sdx, # # upper_2_5sd = mu + 2.5 * sdx # ) # })) # # pure_sd_ranges # Calculate means, SDs, and ±1 SD ranges for each outcome y_limits <- outcomes %>% rowwise() %>% mutate( n = sum(!is.na(d[[var]])), mean = mean(d[[var]], na.rm = TRUE), sd = sd(d[[var]], na.rm = TRUE), # Pure ±1 SD (no scale constraints) pure_lower_1sd = mean - 1 * sd, pure_upper_1sd = mean + 1 * sd, # Bounded ±1 SD (for plotting within scale limits) bounded_lower_1sd = max(ymin, mean - 1 * sd), bounded_upper_1sd = min(ymax, mean + 1 * sd) ) %>% ungroup() y_limits %>% select( var, ymin, ymax, n, mean, sd, pure_lower_1sd, pure_upper_1sd, bounded_lower_1sd, bounded_upper_1sd ) # One-outcome function (model -> EMMs -> Bonferroni pairs -> bracket y positions) analyze_one <- function(df, var, outcome_name) { # Use complete cases for THIS outcome + covariates + condition df2 <- df %>% select(cond, Age, Education, all_of(var)) %>% filter(!is.na(cond), !is.na(Age), !is.na(Education), !is.na(.data[[var]])) %>% mutate(cond = droplevels(cond)) # If an outcome is missing in one condition, skip it (prevents lm contrasts error) if (nlevels(df2$cond) < 2) { message("Skipping ", outcome_name, " (<2 condition levels after dropping NAs).") return(NULL) } # Fit model fml <- as.formula(paste0(var, " ~ cond + Age + Education")) model <- lm(fml, data = df2) # EMMs and 95% CI means <- emmeans(model, specs = "cond") sum_df <- as.data.frame(means) %>% transmute( outcome = outcome_name, cond, mean = emmean, ci_low = lower.CL, ci_high = upper.CL ) # Pairwise comparisons with Bonferroni pairs_df <- as.data.frame(pairs(means, adjust = "bonferroni")) %>% mutate(outcome = outcome_name) %>% tidyr::separate(contrast, into = c("group1", "group2"), sep = " - ") %>% mutate( group1 = gsub("^\\(|\\)$", "", trimws(group1)), group2 = gsub("^\\(|\\)$", "", trimws(group2)), xmin = group2, xmax = group1 ) # Bracket positions: above the highest CI, with specific spacing top_y <- max(sum_df$ci_high, na.rm = TRUE) span <- max(sum_df$ci_high, na.rm = TRUE) - min(sum_df$ci_low, na.rm = TRUE) step <- ifelse(is.finite(span) && span > 0, 0.3 * span, 0.4) pairs_df <- pairs_df %>% arrange(p.value) %>% mutate( y.position = top_y + row_number() * step, p.adj.signif = case_when( p.value < 0.001 ~ "***", p.value < 0.01 ~ "**", p.value < 0.05 ~ "*", TRUE ~ "ns" ), y.position = pmin(y.position, 5.5) ) list(sum_df = sum_df, pairs_df = pairs_df) } # Run function across outcomes res <- pmap(outcomes[, c("var", "outcome")], ~ analyze_one(d, ..1, ..2)) %>% compact() sum_all <- bind_rows(map(res, "sum_df")) pairs_all <- bind_rows(map(res, "pairs_df")) # Keep facet order identical to outcomes table sum_all <- sum_all %>% mutate(outcome = factor(outcome, levels = outcomes$outcome)) pairs_all <- pairs_all %>% mutate(outcome = factor(outcome, levels = outcomes$outcome)) # Build per-facet y-scales (1–6, 1–7, 1–10) using ggh4x # Create a list of y scales in the same order as facet levels # y_scales <- lapply(seq_len(nrow(outcomes)), function(i) { # scale_y_continuous( # limits = c(outcomes$ymin[i], outcomes$ymax[i]), # breaks = outcomes$ymin[i]:outcomes$ymax[i], # expand = c(0, 0) # ) # }) y_scales <- lapply(seq_len(nrow(outcomes)), function(i) { scale_y_continuous( limits = c(2, 5.6), breaks = 2:5, expand = c(0, 0) ) }) # Plotting all outcomes with faceted y-scales plot <- ggplot(sum_all, aes(x = cond, y = mean)) + geom_col(aes(fill = cond), width = 0.6, color = NA) + geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.15, linewidth = 0.6) + stat_pvalue_manual( pairs_all, label = "p.adj.signif", xmin = "xmin", xmax = "xmax", y.position = "y.position", tip.length = 0.02, size = 3.6 ) + facet_wrap(~ outcome, ncol = 3, scales = "free_y") + ggh4x::facetted_pos_scales(y = y_scales) + labs( x = "", y = "Estimated marginal mean (95% CI)" ) + scale_fill_manual( values = c( "Control" = "grey80", "Cracks" = "grey60", "Gradual Build-Up" = "grey40" ), guide = "none" # hides legend (since x-axis already shows condition) ) + theme_minimal(base_size = 13) + theme( text = element_text(family = "Times New Roman"), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 12, angle = 0, hjust = 0.5), axis.title.x = element_text(size = 13), axis.title.y = element_text(size = 13), strip.text = element_text(size = 14, face = "bold"), panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.minor.y = element_blank(), panel.spacing.y = unit(1.5, "lines"), panel.spacing.x = unit(0.8, "lines"), axis.line.x = element_line(linewidth = 0.6, colour = "black"), axis.line.y = element_line(linewidth = 0.6, colour = "black") ) plot ggsave("all_outcomes_faceted_bonferroni_zoomed.png", plot, width = 11, height = 9.5, dpi = 300) # VI/2. Visualization of all outcome variables using their original y-axis scales ---- # Recode condition variable for plotting d <- d %>% mutate(cond = factor(cond, levels = c(0, 1, 2), labels = c("Control", "Cracks", "Gradual Build-Up"))) # Define outcomes and their scale ranges psych::describe(d$agreement) psych::describe(d$authoritarian_tone) psych::describe(d$concern) psych::describe(d$moral_outrage) psych::describe(d$perceived_threat) psych::describe(d$regime_instability) psych::describe(d$lawview) psych::describe(d$government_negativity) psych::describe(d$civil_positivity) outcomes <- tibble::tribble( ~var, ~outcome, ~ymin, ~ymax, "agreement", "Agreement", 1, 6, "authoritarian_tone", "Perceived authoritarian tone", 1, 6, "concern", "Concern", 1, 6, "moral_outrage", "Moral outrage", 1, 6, "perceived_threat", "Perceived threat", 1, 6, "regime_instability", "Regime instability", 1, 6, "lawview", "Agreement with proposed law", 1, 10, "government_negativity","Government negativity", 1, 7, "civil_positivity", "Civil-society positivity", 1, 7 ) # Calculate means, SDs, and ±2 SD ranges for each outcome # y_limits <- outcomes %>% # rowwise() %>% # mutate( # mean = mean(d[[var]], na.rm = TRUE), # sd = sd(d[[var]], na.rm = TRUE), # y_lower = max(ymin, mean - 2 * sd), # y_upper = min(ymax, mean + 2 * sd) # ) %>% # ungroup() # # y_limits %>% # select(var, ymin, ymax, mean, sd, y_lower, y_upper) # # # pure_sd_ranges <- bind_rows(lapply(outcomes$var, function(v) { # x <- d[[v]] # mu <- mean(x, na.rm = TRUE) # sdx <- sd(x, na.rm = TRUE) # # tibble( # outcome = v, # n = sum(!is.na(x)), # mean = mu, # sd = sdx, # lower_2sd = mu - 2 * sdx, # upper_2sd = mu + 2 * sdx # # lower_2_5sd = mu - 2.5 * sdx, # # upper_2_5sd = mu + 2.5 * sdx # ) # })) # # pure_sd_ranges # Calculate means, SDs, and ±1 SD ranges for each outcome y_limits <- outcomes %>% rowwise() %>% mutate( n = sum(!is.na(d[[var]])), mean = mean(d[[var]], na.rm = TRUE), sd = sd(d[[var]], na.rm = TRUE), # Pure ±1 SD (no scale constraints) pure_lower_1sd = mean - 1 * sd, pure_upper_1sd = mean + 1 * sd, # Bounded ±1 SD (for plotting within scale limits) bounded_lower_1sd = max(ymin, mean - 1 * sd), bounded_upper_1sd = min(ymax, mean + 1 * sd) ) %>% ungroup() y_limits %>% select( var, ymin, ymax, n, mean, sd, pure_lower_1sd, pure_upper_1sd, bounded_lower_1sd, bounded_upper_1sd ) # One-outcome function (model -> EMMs -> Bonferroni pairs -> bracket y positions) analyze_one <- function(df, var, outcome_name) { # Use complete cases for THIS outcome + covariates + condition df2 <- df %>% select(cond, Age, Education, all_of(var)) %>% filter(!is.na(cond), !is.na(Age), !is.na(Education), !is.na(.data[[var]])) %>% mutate(cond = droplevels(cond)) # If an outcome is missing in one condition, skip it (prevents lm contrasts error) if (nlevels(df2$cond) < 2) { message("Skipping ", outcome_name, " (<2 condition levels after dropping NAs).") return(NULL) } # Fit model fml <- as.formula(paste0(var, " ~ cond + Age + Education")) model <- lm(fml, data = df2) # EMMs and 95% CI means <- emmeans(model, specs = "cond") sum_df <- as.data.frame(means) %>% transmute( outcome = outcome_name, cond, mean = emmean, ci_low = lower.CL, ci_high = upper.CL ) # Pairwise comparisons with Bonferroni pairs_df <- as.data.frame(pairs(means, adjust = "bonferroni")) %>% mutate(outcome = outcome_name) %>% tidyr::separate(contrast, into = c("group1", "group2"), sep = " - ") %>% mutate( group1 = gsub("^\\(|\\)$", "", trimws(group1)), group2 = gsub("^\\(|\\)$", "", trimws(group2)), xmin = group2, xmax = group1 ) # Bracket positions: above the highest CI, with specific spacing top_y <- max(sum_df$ci_high, na.rm = TRUE) span <- max(sum_df$ci_high, na.rm = TRUE) - min(sum_df$ci_low, na.rm = TRUE) step <- ifelse(is.finite(span) && span > 0, 0.5 * span, 0.6) pairs_df <- pairs_df %>% arrange(p.value) %>% mutate( y.position = top_y + row_number() * step, p.adj.signif = case_when( p.value < 0.001 ~ "***", p.value < 0.01 ~ "**", p.value < 0.05 ~ "*", TRUE ~ "ns" ) ) list(sum_df = sum_df, pairs_df = pairs_df) } # Run function across outcomes res <- pmap(outcomes[, c("var", "outcome")], ~ analyze_one(d, ..1, ..2)) %>% compact() sum_all <- bind_rows(map(res, "sum_df")) pairs_all <- bind_rows(map(res, "pairs_df")) # Keep facet order identical to outcomes table sum_all <- sum_all %>% mutate(outcome = factor(outcome, levels = outcomes$outcome)) pairs_all <- pairs_all %>% mutate(outcome = factor(outcome, levels = outcomes$outcome)) # Build per-facet y-scales (1–6, 1–7, 1–10) using ggh4x # Create a list of y scales in the same order as facet levels y_scales <- lapply(seq_len(nrow(outcomes)), function(i) { scale_y_continuous( limits = c(outcomes$ymin[i], outcomes$ymax[i]), breaks = outcomes$ymin[i]:outcomes$ymax[i], expand = c(0, 0) ) }) # Plotting all outcomes with faceted y-scales plot <- ggplot(sum_all, aes(x = cond, y = mean)) + geom_col(aes(fill = cond), width = 0.6, color = NA) + geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.15, linewidth = 0.6) + stat_pvalue_manual( pairs_all, label = "p.adj.signif", xmin = "xmin", xmax = "xmax", y.position = "y.position", tip.length = 0.02, size = 3.6 ) + facet_wrap(~ outcome, ncol = 3, scales = "free_y") + ggh4x::facetted_pos_scales(y = y_scales) + labs( x = "", y = "Estimated marginal mean (95% CI)" ) + scale_fill_manual( values = c( "Control" = "grey80", "Cracks" = "grey60", "Gradual Build-Up" = "grey40" ), guide = "none" # hides legend (since x-axis already shows condition) ) + theme_minimal(base_size = 13) + theme( text = element_text(family = "Times New Roman"), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 12, angle = 0, hjust = 0.5), axis.title.x = element_text(size = 13), axis.title.y = element_text(size = 13), strip.text = element_text(size = 14, face = "bold"), panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.minor.y = element_blank(), panel.spacing.y = unit(1.5, "lines"), panel.spacing.x = unit(0.8, "lines"), axis.line.x = element_line(linewidth = 0.6, colour = "black"), axis.line.y = element_line(linewidth = 0.6, colour = "black") ) plot ggsave("all_outcomes_faceted_bonferroni.png", plot, width = 11, height = 9.5, dpi = 300)