### ======================================================================== ### ### ### ### Bringing Social Identity Theory Back In: ### ### A New Measure of Issue-Based Affective Polarization ### ### ### ### REPLICATION SCRIPT - FINAL (submission version) ### ### ### ### Authors: [anonymized for review] ### ### ### ### Data: [Anonymized for review] (2024). Swiss population survey on ### ### opinion-based polarization. [Institution anonymized.] ### ### Fielded November 14 - December 15, 2023 (N = 4,001). ### ### ### ### Requirements: R >= 4.2 (developed under 4.5.1); packages listed in ### ### Section 0.1 are installed automatically if missing. ### ### ### ### Layout expected by this script: ### ### ./data/Data_polarization_aufbereitet_2.dta ### ### ./output/ (created automatically) ### ### ### ### ---------------------------------------------------------------------- ### ### OUTPUT -> MANUSCRIPT TABLE/FIGURE MAP ### ### ### ### MAIN TEXT ### ### Table1_CFA_fit.csv .............. Table 1 (CFA model fit) ### ### Table2_reliability.csv .......... Table 2 (scale reliability) ### ### Table3_APSIT_descriptives.csv ... Table 3 (APSIT descriptives) ### ### Figure1_validation_combined ..... Figure 1 (criterion validity) ### ### Figure2_incremental_validity .... Figure 2 (incremental validity) ### ### Figure3_group_effects ........... Figure 3 (camp-specific slopes) ### ### ### ### APPENDIX ### ### TableB1_loadings.csv ............ Table B1 ### ### TableB2_model_comparison.csv .... Table B2 ### ### TableC1_reliability_extended.csv Table C1 ### ### FigureD1_APSIT_distribution ..... Figure D1 ### ### FigureD2_aggregation_* .......... Figure D2 ### ### FigureD3_aggregation_* .......... Figure D3 ### ### TableG1_sample_descriptives.csv . Table G1 ### ### TableH1_H2_convergent_validity .. Tables H1 and H2 ### ### TableH3_criterion_validity.csv .. Table H3 ### ### TablesI1_I4_separate_models.csv . Tables I1-I4 ### ### TablesJ1_J4_full_models.csv ..... Tables J1-J4 ### ### TableJ5_slope_differences.csv ... Table J5 ### ### FigK_parallel_*.png ............. Figures K1-K4 ### ### ### ### Structure: ### ### 0. Setup & data preparation ### ### 1. Part I: Scale construction (CFA, reliability, APSIT index) ### ### 2. Part II: Validation (convergent, criterion) ### ### 3. Part III: Democratic attitudes & camp-specific effects ### ### 4. Figures ### ### 5. Appendix tables & figures ### ### 6. Robustness checks ### ### ### ### ======================================================================== ### # ══════════════════════════════════════════════════════════════════════════════ # 0. SETUP & DATA PREPARATION # ══════════════════════════════════════════════════════════════════════════════ # --- 0.1 Packages ----------------------------------------------------------- required_packages <- c( "haven", "tidyr", "tibble", "purrr", "ggplot2", "patchwork", "lavaan", "semTools", "psych", "sandwich", "lmtest", "MASS", "lme4", "dplyr" ) new_pkgs <- setdiff(required_packages, rownames(installed.packages())) if (length(new_pkgs)) install.packages(new_pkgs, dependencies = TRUE) # Load dplyr LAST to avoid select() being masked invisible(lapply(required_packages, library, character.only = TRUE)) select <- dplyr::select filter <- dplyr::filter set.seed(12345) # --- 0.2 Data path ----------------------------------------------------------- # Replication package layout (relative to this script's folder): # ./data/Data_polarization_aufbereitet_2.dta (survey data, see README) # ./output/ (all tables and figures) data_path <- file.path("data", "Data_polarization_aufbereitet_2.dta") # --- 0.3 Output directory ---------------------------------------------------- output_dir <- "output" dir.create(output_dir, showWarnings = FALSE, recursive = TRUE) out <- function(filename) { file.path(output_dir, filename) } cat("Outputs will be saved to:\n", output_dir, "\n") # Console log: everything printed below is mirrored to output/console_log.txt # (split = TRUE keeps the live console output). If the script aborts midway, # run `while (sink.number() > 0) sink()` once before re-running. log_con <- file(out("console_log.txt"), open = "wt") sink(log_con, split = TRUE) cat("Run started:", format(Sys.time()), "\n") # --- 0.4 Load & clean data -------------------------------------------------- dat <- haven::read_dta(data_path) |> as.data.frame() # Recode missing values # NOTE on missing codes (verified against the codebook): the prepared .dta # encodes item non-response as -77/-88/-99. In addition, 99 = "don't know" in # ALL issue batteries: 1-5 item batteries (ident*/resent*/deprivation*/ # threat*/char*/politik*/hochzeit*), the 0-10 feeling thermometers (thermo*), # and the 0-10 self-placements (position*). The earlier assumption that # thermometers run 0-100 (making 99 a legitimate value) was incorrect; the # survey thermometers are 0-10 scales. 99 is therefore recoded to NA for all # issue-battery variables in Section 0.5b below (with a per-family count so # any impact on previously reported numbers is visible), followed by hard # range assertions. 99 is NOT recoded globally (e.g., alter can be 99). missing_codes <- c(-77, -88, -99) for (v in names(dat)) { if (is.numeric(dat[[v]])) dat[[v]][dat[[v]] %in% missing_codes] <- NA } # Gender: haven_labelled -> numeric factor levels, recode to 0/1 # Original coding: 1 = male, 2 = female, 3/4 = other/missing # Gender: ROBUST recode to 1 = male, 0 = female. # FIX (v3.1): the previous recode assumed raw survey coding 1 = male, # 2 = female (per codebook core_geschlecht). The prepared .dta evidently uses # a different coding (the v3 run produced Female_pct = 0 and halved every # regression N), so one category fell through to NA. This block detects the # actual coding from the data/labels, prints a diagnosis, and hard-stops if # the recode would lose more than 10% of cases. g_raw <- dat$geschlecht g_lab <- attr(g_raw, "labels", exact = TRUE) g_num <- suppressWarnings(as.numeric(haven::zap_labels(g_raw))) cat("0.4 geschlecht: raw value table:\n"); print(table(g_num, useNA = "ifany")) if (!is.null(g_lab)) { cat("0.4 geschlecht: value labels:\n"); print(g_lab) } find_code <- function(patterns) { if (is.null(g_lab)) return(NA_real_) nm <- tolower(names(g_lab)) hit <- Reduce(`|`, lapply(patterns, function(p) grepl(p, nm))) if (any(hit)) unname(g_lab[which(hit)[1]]) else NA_real_ } # Order matters: match female first so "male" cannot swallow "female". code_female <- find_code(c("weib", "frau", "female", "femme", "donna")) code_male <- find_code(c("m\u00e4nn", "maenn", "mann\\b", "\\bmale\\b", "homme", "uomo")) vals <- sort(unique(stats::na.omit(g_num))) if (!is.na(code_male) && !is.na(code_female)) { dat$geschlecht <- ifelse(g_num == code_male, 1, ifelse(g_num == code_female, 0, NA)) cat("0.4 geschlecht: recoded via value labels (male =", code_male, ", female =", code_female, ")\n") } else if (identical(vals, c(1, 2))) { dat$geschlecht <- ifelse(g_num == 1, 1, ifelse(g_num == 2, 0, NA)) cat("0.4 geschlecht: recoded assuming raw codebook coding 1=male, 2=female\n") } else if (identical(vals, c(0, 1))) { dat$geschlecht <- ifelse(g_num %in% c(0, 1), g_num, NA) cat("0.4 geschlecht: 0/1 coding WITHOUT labels detected; kept as is,\n", " ASSUMING 1 = male, 0 = female. VERIFY against the preparation\n", " syntax; if reversed, only the control's sign flips (substantive\n", " coefficients are unaffected).\n") } else { stop("Unexpected geschlecht coding: values ", paste(vals, collapse = ", "), " - inspect the .dta and adapt this block.") } stopifnot("geschlecht recode lost >10% of all cases (coding mismatch?)" = sum(!is.na(dat$geschlecht)) >= 0.9 * length(g_num)) cat("0.4 geschlecht: n male =", sum(dat$geschlecht == 1, na.rm = TRUE), "| n female =", sum(dat$geschlecht == 0, na.rm = TRUE), "| NA =", sum(is.na(dat$geschlecht)), "\n") # Now: 1 = male, 0 = female, NA = other/missing # Age, education: ensure numeric dat$alter <- suppressWarnings(as.numeric(dat$alter)) dat$bildung <- suppressWarnings(as.numeric(dat$bildung)) # Income x <- dat$einkommenssituation if (inherits(x, "haven_labelled")) x <- haven::zap_labels(x) dat$einkommenssituation <- suppressWarnings(as.numeric(x)) # Language region (reference = German) dat$sprachregion <- factor(dat$sprache, levels = 1:3, labels = c("German", "French", "Italian")) dat$sprachregion <- relevel(dat$sprachregion, ref = "German") # Left-right self-placement (z-standardized) dat$rile <- suppressWarnings(as.numeric(haven::zap_labels(dat$rile))) dat$rile_z <- as.numeric(scale(dat$rile)) # --- 0.5 Issue domain definitions -------------------------------------------- # Uses PRE-AGGREGATED variables from the dataset where available issues <- list( Migration = list( ident = c("identmigration1", "identmigration3", "identmigration4", "identmigration5"), comp = c("resentmigration4", "deprivationmigration1", "threatmigration2", "threatmigration3"), ident_m = "identmigrationm", comp_m = "resentmigrationm", pos = "positionmigration", opp = "oppmigration", procon = "proconmigration", tol_items = paste0("poltoleranzmigration", 1:5), tol_m = "poltoleranzmigrationm", thermo_b = "thermomigrationbefuerworter", thermo_g = "thermomigrationgegner", ap1 = "ap1_migration", char_b = paste0("charmigrationbefuerworter", 1:4), char_g = paste0("charmigrationgegner", 1:4), pol_b = "politikmigrationbefuerworter", pol_g = "politikmigrationgegner", priv_b = "hochzeitmigrationbefuerworter", priv_g = "hochzeitmigrationgegner" ), Climate = list( ident = c("identclimate1", "identclimate3", "identclimate4", "identclimate5"), comp = c("resentclimate4", "deprivationclimate1", "threatclimate2", "threatclimate3"), ident_m = "identclimatem", comp_m = "resentclimatem", pos = "positionclimate", opp = "oppclimate", procon = "proconclimate", tol_items = paste0("poltoleranzclimate", 1:5), tol_m = "poltoleranzclimatem", thermo_b = "thermoclimatebefuerworter", thermo_g = "thermoclimategegner", ap1 = "ap1_climate", char_b = paste0("charclimatebefuerworter", 1:4), char_g = paste0("charclimategegner", 1:4), pol_b = "politikclimatebefuerworter", pol_g = "politikclimategegner", priv_b = "hochzeitclimatebefuerworter", priv_g = "hochzeitclimategegner" ), Vaccination = list( ident = c("identimpfung1", "identimpfung3", "identimpfung4", "identimpfung5"), comp = c("resentimpfung4", "deprivationimpfung1", "threatimpfung2", "threatimpfung3"), ident_m = "identimpfungm", comp_m = "resentimpfungm", pos = "positionimpfung", opp = "oppimpfung", procon = "proconimpfung", tol_items = paste0("poltoleranzimpfung", 1:5), tol_m = "poltoleranzimpfungm", thermo_b = "thermoimpfungbefuerworter", thermo_g = "thermoimpfunggegner", ap1 = "ap1_impfung", char_b = paste0("charimpfungbefuerworter", 1:4), char_g = paste0("charimpfunggegner", 1:4), pol_b = "politikimpfungbefuerworter", pol_g = "politikimpfunggegner", priv_b = "hochzeitimpfungbefuerworter", priv_g = "hochzeitimpfunggegner" ), Gender = list( ident = c("identgender1", "identgender3", "identgender4", "identgender5"), comp = c("resentgender4", "deprivationgender1", "threatgender2", "threatgender3"), ident_m = "identgenderm", comp_m = "resentgenderm", pos = "positiongender", opp = "oppgender", procon = "procongender", tol_items = paste0("poltoleranzgender", 1:5), tol_m = "poltoleranzgenderm", thermo_b = "thermogenderbefuerworter", thermo_g = "thermogendergegner", ap1 = "ap1_gender", char_b = paste0("chargenderbefuerworter", 1:4), char_g = paste0("chargendergegner", 1:4), pol_b = "politikgenderbefuerworter", pol_g = "politikgendergegner", priv_b = "hochzeitgenderbefuerworter", priv_g = "hochzeitgendergegner" ) ) # --- 0.5b Recode 99 ("don't know") in issue batteries; range assertions ------ # Codebook: 99 = don't know for all issue-battery variables. Expected ranges: # 1-5 : ident*, resent*/deprivation*/threat* (comp), char*, politik*, hochzeit* # 0-10 : thermo* (feeling thermometers), position* (self-placements) issue_vars_15 <- unique(unlist(lapply(issues, function(cfg) c(cfg$ident, cfg$comp, cfg$char_b, cfg$char_g, cfg$pol_b, cfg$pol_g, cfg$priv_b, cfg$priv_g)))) issue_vars_010 <- unique(unlist(lapply(issues, function(cfg) c(cfg$thermo_b, cfg$thermo_g, cfg$pos)))) recode99 <- function(vars, lo, hi, label) { vars <- vars[vars %in% names(dat)] n99 <- 0L; n_oob <- 0L for (v in vars) { x <- suppressWarnings(as.numeric(haven::zap_labels(dat[[v]]))) n99 <- n99 + sum(x == 99, na.rm = TRUE) x[x == 99] <- NA n_oob <- n_oob + sum(x < lo | x > hi, na.rm = TRUE) dat[[v]] <<- x } cat(sprintf("0.5b %-28s: %d values of 99 recoded to NA; %d out-of-range [%g,%g]\n", label, n99, n_oob, lo, hi)) invisible(n_oob) } oob1 <- recode99(issue_vars_15, 1, 5, "1-5 item batteries") oob2 <- recode99(issue_vars_010, 0, 10, "0-10 thermometers/positions") # Hard stop if anything remains outside the documented ranges: stopifnot("Out-of-range values remain in issue batteries" = (oob1 + oob2) == 0) # NOTE: if the 99-counts above are > 0, previously reported estimates that # used these variables were contaminated and WILL shift; compare old vs. new # outputs cell by cell. # --- 0.6 Helper functions ---------------------------------------------------- row_mean_min <- function(df, items, min_valid = 3) { m <- as.matrix(df[, items, drop = FALSE]) m[] <- suppressWarnings(apply(m, 2, as.numeric)) nn <- rowSums(!is.na(m)) out <- rowMeans(m, na.rm = TRUE) out[nn < min_valid] <- NA_real_ out } alpha_ord <- function(df, items) { X <- df[, items, drop = FALSE] X[] <- lapply(X, function(x) suppressWarnings(as.numeric(x))) R <- psych::polychoric(X)$rho psych::alpha(R, n.obs = nrow(na.omit(X)))$total$raw_alpha } extract_hc3 <- function(model, var) { ct <- lmtest::coeftest(model, vcov. = sandwich::vcovHC(model, type = "HC3")) if (!var %in% rownames(ct)) return(list(b = NA_real_, se = NA_real_, p = NA_real_)) list(b = unname(ct[var, 1]), se = unname(ct[var, 2]), p = unname(ct[var, 4])) } controls_rhs <- "geschlecht + alter + bildung + einkommenssituation + rile_z + sprachregion" # --- 0.7 Construct dependent variables (shared across issues) ---------------- # Democratic support (reversed items) items_support <- c("demsupport1", "demsupport2", "demsupport3") for (v in items_support) dat[[v]] <- suppressWarnings(as.numeric(dat[[v]])) dat$demsupport2_rev <- ifelse(is.na(dat$demsupport2), NA, 6 - dat$demsupport2) dat$demsupport3_rev <- ifelse(is.na(dat$demsupport3), NA, 6 - dat$demsupport3) dat$demo_support_mean <- row_mean_min( dat, c("demsupport1", "demsupport2_rev", "demsupport3_rev"), min_valid = 1 ) # Democracy satisfaction (evaluation - conception, using pre-aggregated vars) dat$democracy_satisfaction <- as.numeric(dat$democracyevaluationm) - as.numeric(dat$democracyconceptionm) # Institutional trust trust_items <- c("poltrustbundesrat", "poltrustparlament", "poltrustjustiz") trust_vars <- trust_items[trust_items %in% names(dat)] X_trust <- dat[, trust_vars, drop = FALSE] X_trust[] <- lapply(X_trust, function(x) { suppressWarnings(as.numeric( if (inherits(x, "haven_labelled")) haven::zap_labels(x) else x )) }) dat$trust_mean <- rowMeans(X_trust, na.rm = TRUE) # Political tolerance (per issue): 6 - mean of items (higher = more tolerant) for (issue_name in names(issues)) { cfg <- issues[[issue_name]] sfx <- tolower(issue_name) for (v in cfg$tol_items) { dat[[v]] <- suppressWarnings(as.numeric(dat[[v]])) dat[[v]] <- dplyr::na_if(dat[[v]], 99) } dat[[paste0("poltolerance_", sfx)]] <- 6 - rowMeans( dat[, cfg$tol_items, drop = FALSE], na.rm = TRUE ) } # --- 0.8 Construct issue-specific predictors --------------------------------- for (issue_name in names(issues)) { cfg <- issues[[issue_name]] sfx <- tolower(issue_name) # Ensure raw items are numeric for (v in c(cfg$ident, cfg$comp)) { dat[[v]] <- suppressWarnings(as.numeric(dat[[v]])) } # Build subscale means from theoretically-justified items # (per manuscript Tables A1/A2; item selection documented in the # measurement appendix [reference anonymized for review]) # # Identity items (4): ident1, _3, _4, _5 # (Item 2 dropped per EFA; see manuscript Appendix A) # Comparison items (4): resent4, deprivation1, # threat2, threat3 # (mixed-family, theoretically derived from SIT) # # Per-person mean if at least 3 valid responses (manuscript Table 3 note) min_valid <- 3 id_mat <- as.matrix(dat[, cfg$ident, drop = FALSE]) co_mat <- as.matrix(dat[, cfg$comp, drop = FALSE]) n_valid_id <- rowSums(!is.na(id_mat)) n_valid_co <- rowSums(!is.na(co_mat)) id_mean <- ifelse(n_valid_id >= min_valid, rowMeans(id_mat, na.rm = TRUE), NA_real_) co_mean <- ifelse(n_valid_co >= min_valid, rowMeans(co_mat, na.rm = TRUE), NA_real_) # Store raw means for downstream use (replaces cfg$ident_m / cfg$comp_m) dat[[paste0("identity_", sfx, "_mean")]] <- id_mean dat[[paste0("comparison_", sfx, "_mean")]] <- co_mean # Z-standardized subscales (pool-wide) dat[[paste0("identity_", sfx, "_z")]] <- as.numeric(scale(id_mean)) dat[[paste0("comparison_", sfx, "_z")]] <- as.numeric(scale(co_mean)) # APSIT index: mean of z-standardized subscales, then z-standardized id_z <- dat[[paste0("identity_", sfx, "_z")]] co_z <- dat[[paste0("comparison_", sfx, "_z")]] ap <- rep(NA_real_, nrow(dat)) ok <- complete.cases(id_z, co_z) ap[ok] <- rowMeans(cbind(id_z[ok], co_z[ok])) dat[[paste0("APSIT_", sfx)]] <- as.numeric(scale(ap)) # Group variable from opp (0 = supporter, 1 = opponent) opp_var <- cfg$opp if (opp_var %in% names(dat)) { v <- suppressWarnings(as.numeric(haven::zap_labels(dat[[opp_var]]))) dat[[paste0("group_", sfx)]] <- factor( ifelse(is.na(v), NA_character_, ifelse(v == 0, "Supporter", ifelse(v == 1, "Opponent", NA_character_))), levels = c("Supporter", "Opponent") ) } grp <- dat[[paste0("group_", sfx)]] id_z_val <- dat[[paste0("identity_", sfx, "_z")]] co_z_val <- dat[[paste0("comparison_", sfx, "_z")]] # Group-specific predictors (for Figure 3) dat[[paste0("id_sup_", sfx)]] <- ifelse(grp == "Supporter", id_z_val, 0) dat[[paste0("id_opp_", sfx)]] <- ifelse(grp == "Opponent", id_z_val, 0) dat[[paste0("co_sup_", sfx)]] <- ifelse(grp == "Supporter", co_z_val, 0) dat[[paste0("co_opp_", sfx)]] <- ifelse(grp == "Opponent", co_z_val, 0) # Camp-consistent Feeling Thermometer (using procon variable) tb <- cfg$thermo_b; tg <- cfg$thermo_g procon_var <- cfg$procon if (all(c(tb, tg, procon_var) %in% names(dat))) { camp <- suppressWarnings(as.numeric(haven::zap_labels(dat[[procon_var]]))) tb_num <- suppressWarnings(as.numeric(dat[[tb]])) tg_num <- suppressWarnings(as.numeric(dat[[tg]])) # procon: 1 = befuerworter camp, 0 = gegner camp dat[[paste0("Tin_", sfx)]] <- ifelse(camp == 1, tb_num, ifelse(camp == 0, tg_num, NA_real_)) dat[[paste0("Tout_", sfx)]] <- ifelse(camp == 1, tg_num, ifelse(camp == 0, tb_num, NA_real_)) dat[[paste0("FT_", sfx)]] <- dat[[paste0("Tin_", sfx)]] - dat[[paste0("Tout_", sfx)]] } # Position extremity (symmetric folding around 5) if (cfg$pos %in% names(dat)) { pos <- suppressWarnings(as.numeric(dat[[cfg$pos]])) dat[[paste0("pos_extreme_", sfx)]] <- abs(pos - 5) } # Construct social distance DVs if (all(c(cfg$char_b, cfg$char_g) %in% names(dat))) { char_diffs <- sapply(1:4, function(i) { abs(as.numeric(dat[[cfg$char_b[i]]]) - as.numeric(dat[[cfg$char_g[i]]])) }) dat[[paste0("char_dist_", sfx)]] <- rowMeans(char_diffs, na.rm = TRUE) } if (all(c(cfg$pol_b, cfg$pol_g) %in% names(dat))) { dat[[paste0("pol_dist_", sfx)]] <- abs( as.numeric(dat[[cfg$pol_b]]) - as.numeric(dat[[cfg$pol_g]]) ) } if (all(c(cfg$priv_b, cfg$priv_g) %in% names(dat))) { dat[[paste0("priv_dist_", sfx)]] <- abs( as.numeric(dat[[cfg$priv_b]]) - as.numeric(dat[[cfg$priv_g]]) ) } } cat("=== Data preparation complete. N =", nrow(dat), "===\n") # ══════════════════════════════════════════════════════════════════════════════ # 1. PART I: SCALE CONSTRUCTION # ══════════════════════════════════════════════════════════════════════════════ # --- 1.1 Confirmatory Factor Analysis (Table 1, Table B1, Table B2) ---------- cfa_results <- list() for (issue_name in names(issues)) { cfg <- issues[[issue_name]] all_items <- c(cfg$ident, cfg$comp) dat_cfa <- dat[, all_items, drop = FALSE] dat_cfa[] <- lapply(dat_cfa, function(x) { factor(suppressWarnings(as.numeric(x)), ordered = TRUE) }) # Two-factor model (correlated) — primary measurement model mod_2f <- paste0( "Identity =~ ", paste(cfg$ident, collapse = " + "), "\n", "Comparison =~ ", paste(cfg$comp, collapse = " + "), "\n", "Identity ~~ Comparison" ) # One-factor model (for comparison) mod_1f <- paste0("G =~ ", paste(all_items, collapse = " + ")) fit_2f <- cfa(mod_2f, data = dat_cfa, ordered = all_items, estimator = "WLSMV") fit_1f <- cfa(mod_1f, data = dat_cfa, ordered = all_items, estimator = "WLSMV") get_fit <- function(fit) { # NOTE: Table 1 in the manuscript reports the mean- and variance-adjusted # (scaled) chi-square together with the conventional CFI/TLI/RMSEA/SRMR of # the WLSMV solution. Keep this extraction unchanged so the published # table values remain exactly reproducible; switching to the ".scaled"/ # robust index variants would change the reported numbers. m <- fitMeasures(fit, c("chisq.scaled", "df", "cfi", "tli", "rmsea", "srmr")) tibble(chi2 = m["chisq.scaled"], df = m["df"], CFI = m["cfi"], TLI = m["tli"], RMSEA = m["rmsea"], SRMR = m["srmr"]) } # Standardized loadings from the 2F model loads <- parameterEstimates(fit_2f, standardized = TRUE) |> subset(op == "=~", select = c("lhs", "rhs", "std.all")) # Factor correlation from 2F model fcor <- parameterEstimates(fit_2f, standardized = TRUE) |> subset(op == "~~" & lhs == "Identity" & rhs == "Comparison") # Composite reliability (rho_C) from the 2F model rhoC_all <- semTools::compRelSEM(fit_2f, return.total = FALSE) rhoC_id <- as.numeric(rhoC_all["Identity"]) rhoC_co <- as.numeric(rhoC_all["Comparison"]) cfa_results[[issue_name]] <- list( fit_2f = fit_2f, fit_1f = fit_1f, fit_2f_idx = get_fit(fit_2f), fit_1f_idx = get_fit(fit_1f), loadings = loads, fcor = fcor$std.all, rhoC_identity = rhoC_id, rhoC_comparison = rhoC_co ) cat("\n===", issue_name, "CFA ===\n") cat("Two-factor: CFI =", round(get_fit(fit_2f)$CFI, 3), "| TLI =", round(get_fit(fit_2f)$TLI, 3), "| RMSEA =", round(get_fit(fit_2f)$RMSEA, 3), "| SRMR =", round(get_fit(fit_2f)$SRMR, 3), "\n") cat("One-factor: CFI =", round(get_fit(fit_1f)$CFI, 3), "| RMSEA =", round(get_fit(fit_1f)$RMSEA, 3), "\n") cat("Factor correlation =", round(fcor$std.all, 3), "\n") cat("rhoC: Identity =", round(rhoC_id, 3), "| Comparison =", round(rhoC_co, 3), "\n") } # --- TABLE 1: CFA Model Fit Comparison (main text) --------------------------- table1 <- bind_rows(lapply(names(cfa_results), function(nm) { r <- cfa_results[[nm]] bind_rows( r$fit_2f_idx |> mutate(Domain = nm, Model = "Two-factor", .before = 1), r$fit_1f_idx |> mutate(Domain = nm, Model = "One-factor", .before = 1) ) })) cat("\n=== Table 1: CFA Model Fit ===\n") print(table1 |> mutate(across(where(is.numeric), ~round(.x, 3)))) write.csv(table1, out("Table1_CFA_fit.csv"), row.names = FALSE) # --- TABLE B2: Two- vs. One-Factor Model Comparison (Appendix B) -------------- # REPLACES the previous hand-written Table B2 in the manuscript, which # contained placeholder values. Fit indices come from the SAME get_fit() # extraction as Table 1, so Table 1 and Table B2 are identical by construction # on the shared columns. The chi-square difference test uses the scaled # difference method appropriate for WLSMV (lavTestLRT, satorra.2000); # arithmetic differences of scaled chi-squares are NOT a valid test. table_b2 <- bind_rows(lapply(names(cfa_results), function(nm) { r <- cfa_results[[nm]] lrt <- tryCatch( lavTestLRT(r$fit_1f, r$fit_2f, method = "satorra.2000"), error = function(e) lavTestLRT(r$fit_1f, r$fit_2f) ) # Version-safe column access (names are "Chisq diff", "Df diff", # "Pr(>Chisq)" in current lavaan; matched by pattern as a safeguard) cn <- colnames(lrt) col_dchi <- cn[grepl("Chisq diff", cn, fixed = TRUE)][1] col_ddf <- cn[grepl("Df diff", cn, fixed = TRUE)][1] col_p <- cn[grepl("Pr(>", cn, fixed = TRUE)][1] diff_row <- lrt[2, , drop = FALSE] bind_rows( r$fit_2f_idx |> mutate(Domain = nm, Model = "Two-factor", .before = 1) |> mutate(Dchi2_scaled = NA_real_, Ddf = NA_real_, p_diff = NA_real_), r$fit_1f_idx |> mutate(Domain = nm, Model = "One-factor", .before = 1) |> mutate(Dchi2_scaled = if (!is.na(col_dchi)) as.numeric(diff_row[[col_dchi]]) else NA_real_, Ddf = if (!is.na(col_ddf)) as.numeric(diff_row[[col_ddf]]) else NA_real_, p_diff = if (!is.na(col_p)) as.numeric(diff_row[[col_p]]) else NA_real_) ) })) cat("\n=== Table B2: Two- vs. One-Factor Comparison (scaled diff test) ===\n") print(table_b2 |> mutate(across(where(is.numeric), ~round(.x, 3)))) write.csv(table_b2, out("TableB2_model_comparison.csv"), row.names = FALSE) # --- TABLE B1: Standardized Factor Loadings (Appendix) ----------------------- table_b1 <- bind_rows(lapply(names(cfa_results), function(nm) { cfa_results[[nm]]$loadings |> mutate(Domain = nm, .before = 1) })) cat("\n=== Table B1: Standardized Loadings ===\n") print(table_b1 |> mutate(std.all = round(std.all, 2))) # --- 1.2 Scale Reliability (Table 2, Table C1) -------------------------------- table2 <- bind_rows(lapply(names(issues), function(nm) { cfg <- issues[[nm]] r <- cfa_results[[nm]] tibble( Domain = nm, `Identity alpha_ord` = alpha_ord(dat, cfg$ident), `Identity rho_C` = r$rhoC_identity, `Comparison alpha_ord` = alpha_ord(dat, cfg$comp), `Comparison rho_C` = r$rhoC_comparison ) })) cat("\n=== Table 2: Scale Reliability ===\n") print(table2 |> mutate(across(where(is.numeric), ~round(.x, 3)))) # Compute omega (Appendix Table C1) # NOTE: semTools::reliability() is deprecated (still functional as of # semTools 0.5-7). If a future semTools release removes it, replace with the # corresponding compRelSEM() call and verify Table C1 values stay identical. for (nm in names(cfa_results)) { fit <- cfa_results[[nm]]$fit_2f rel <- suppressMessages(suppressWarnings(semTools::reliability(fit))) cfa_results[[nm]]$omega_identity <- as.numeric(rel["omega", "Identity"]) cfa_results[[nm]]$omega_comparison <- as.numeric(rel["omega", "Comparison"]) } # --- TABLE C1: Extended Reliability with omega, Mean, SD (Appendix) ----------- table_c1 <- bind_rows(lapply(names(issues), function(nm) { cfg <- issues[[nm]] r <- cfa_results[[nm]] sfx <- tolower(nm) # Use the manuscript-consistent subscale means (min_valid = 3, Section 0.8) id_vals <- dat[[paste0("identity_", sfx, "_mean")]] co_vals <- dat[[paste0("comparison_", sfx, "_mean")]] tibble( Domain = c(nm, ""), Dimension = c("Identity", "Comparison"), Items = c(4, 4), alpha_ord = c(alpha_ord(dat, cfg$ident), alpha_ord(dat, cfg$comp)), rho_C = c(r$rhoC_identity, r$rhoC_comparison), omega = c(r$omega_identity, r$omega_comparison), Mean = c(mean(id_vals, na.rm = TRUE), mean(co_vals, na.rm = TRUE)), SD = c(sd(id_vals, na.rm = TRUE), sd(co_vals, na.rm = TRUE)) ) })) cat("\n=== Table C1: Extended Scale Reliability ===\n") print(table_c1 |> mutate(across(where(is.numeric), ~round(.x, 3)))) write.csv(table_c1, out("TableC1_reliability_extended.csv"), row.names = FALSE) # --- TABLE 3: APSIT Descriptive Statistics by Issue Domain (Main Text) -------- table3_apsit <- bind_rows(lapply(names(issues), function(nm) { sfx <- tolower(nm) apsit_var <- paste0("APSIT_", sfx) vals <- dat[[apsit_var]] vals <- vals[!is.na(vals)] tibble( Domain = nm, N = length(vals), Mean = mean(vals), SD = sd(vals), Min = min(vals), Max = max(vals) ) })) cat("\n=== Table 3: APSIT Descriptive Statistics by Issue Domain ===\n") print(table3_apsit |> mutate(across(where(is.numeric), ~round(.x, 2)))) write.csv(table3_apsit, out("Table3_APSIT_descriptives.csv"), row.names = FALSE) # ══════════════════════════════════════════════════════════════════════════════ # 2. PART II: SCALE VALIDATION # ══════════════════════════════════════════════════════════════════════════════ # --- 2.1 Convergent Validity ------------------------------------------------- conv_coefs <- list() for (issue_name in names(issues)) { sfx <- tolower(issue_name) cfg <- issues[[issue_name]] ap1_var <- cfg$ap1 pos_var <- paste0("pos_extreme_", sfx) if (!ap1_var %in% names(dat)) next dat[[ap1_var]] <- suppressWarnings(as.numeric(dat[[ap1_var]])) # --- AP1 ORIENTATION CHECK (read this console block before final wording) -- # The pre-aggregated ap1_* index has an undocumented orientation. This block # settles it empirically against the self-built thermometer variables from # Section 0.8 (Tin_ = ingroup warmth, Tout_ = outgroup warmth, FT_ = Tin-Tout): # r(ap1, FT) close to +1 -> ap1 = ingroup-minus-outgroup difference # r(ap1, FT) close to -1 -> ap1 = outgroup-minus-ingroup difference # r(ap1, Tout) ~ +1 -> ap1 = raw OUTGROUP warmth (higher = warmer) # Align the Figure 1 note, the Table F1/F2 titles, and the manuscript text # with whichever definition this check confirms. No sign flipping is done # anywhere in this script; coefficients are reported as estimated. ftv <- paste0("FT_", sfx); tinv <- paste0("Tin_", sfx); toutv <- paste0("Tout_", sfx) if (all(c(ftv, tinv, toutv) %in% names(dat))) { a <- dat[[ap1_var]] cat("\n--- AP1 ORIENTATION CHECK:", issue_name, "---\n") cat(" range(ap1) :", paste(range(a, na.rm = TRUE), collapse = " to "), "\n") cat(" r(ap1, FT=in-out):", round(cor(a, dat[[ftv]], use = "complete.obs"), 3), "\n") cat(" r(ap1, Tin) :", round(cor(a, dat[[tinv]], use = "complete.obs"), 3), "\n") cat(" r(ap1, Tout) :", round(cor(a, dat[[toutv]], use = "complete.obs"), 3), "\n") grp_tmp <- dat[[paste0("group_", sfx)]] cat(" mean(ap1) Supporter:", round(mean(a[grp_tmp == "Supporter"], na.rm = TRUE), 3), "| Opponent:", round(mean(a[grp_tmp == "Opponent"], na.rm = TRUE), 3), "\n") } # Build item-level subscale means # DECISION IMPLEMENTED (option b): min_valid = 3 throughout, harmonized with # Section 0.8 and the Table 3 note ("at least three valid items per # subscale"). This resolves the N inconsistency between Table 3 (761/772/ # 817/764) and the previous Tables F1/F2 (892/876/884/869): Part II now uses # the same inclusion rule as Part III. Figure 1 and Tables F1-F3 must be # re-set from this run's outputs. id_items_mean <- paste0("id_items_mean_", sfx) co_items_mean <- paste0("co_items_mean_", sfx) dat[[id_items_mean]] <- row_mean_min(dat, cfg$ident, min_valid = 3) dat[[co_items_mean]] <- row_mean_min(dat, cfg$comp, min_valid = 3) # Build item-level APSIT index ap_items_var <- paste0("APSIT_items_", sfx) id_z_tmp <- as.numeric(scale(dat[[id_items_mean]])) co_z_tmp <- as.numeric(scale(dat[[co_items_mean]])) ap_tmp <- rep(NA_real_, nrow(dat)) ok_tmp <- complete.cases(id_z_tmp, co_z_tmp) ap_tmp[ok_tmp] <- rowMeans(cbind(id_z_tmp[ok_tmp], co_z_tmp[ok_tmp])) dat[[ap_items_var]] <- as.numeric(scale(ap_tmp)) # Model A: APSIT -> FT index f_ap <- as.formula(paste(ap1_var, "~", ap_items_var, "+", controls_rhs)) m_ap <- lm(f_ap, data = dat) # Model B: Identity + Comparison -> FT index f_ic <- as.formula(paste(ap1_var, "~", id_items_mean, "+", co_items_mean, "+", controls_rhs)) m_ic <- lm(f_ic, data = dat) # Model C: APSIT -> Position extremity f_pos_ap <- as.formula(paste(pos_var, "~", ap_items_var, "+", controls_rhs)) m_pos_ap <- lm(f_pos_ap, data = dat) # Model D: Identity + Comparison -> Position extremity f_pos_ic <- as.formula(paste(pos_var, "~", id_items_mean, "+", co_items_mean, "+", controls_rhs)) m_pos_ic <- lm(f_pos_ic, data = dat) cc_row <- function(outcome, predictor, model, var) { e <- extract_hc3(model, var) tibble(issue = issue_name, outcome = outcome, predictor = predictor, b = e$b, se = e$se, p = e$p, N = nobs(model), adj_R2 = summary(model)$adj.r.squared) } conv_coefs[[issue_name]] <- bind_rows( cc_row("Feeling Thermometer", "APSIT", m_ap, ap_items_var), cc_row("Feeling Thermometer", "Identity", m_ic, id_items_mean), cc_row("Feeling Thermometer", "Comparison", m_ic, co_items_mean), cc_row("Self-placement", "APSIT", m_pos_ap, ap_items_var), cc_row("Self-placement", "Identity", m_pos_ic, id_items_mean), cc_row("Self-placement", "Comparison", m_pos_ic, co_items_mean) ) cat("\n=== Convergent Validity:", issue_name, "===\n") cat("APSIT -> FT index: b =", round(extract_hc3(m_ap, ap_items_var)$b, 3), "\n") cat("Identity -> FT index: b =", round(extract_hc3(m_ic, id_items_mean)$b, 3), "\n") cat("Comparison -> FT index: b =", round(extract_hc3(m_ic, co_items_mean)$b, 3), "\n") } conv_coefs_df <- bind_rows(conv_coefs) # --- 2.2 Criterion Validity: Social Distance --------------------------------- crit_coefs <- list() for (issue_name in names(issues)) { sfx <- tolower(issue_name) cfg <- issues[[issue_name]] ap_var <- paste0("APSIT_items_", sfx) id_var <- paste0("id_items_mean_", sfx) co_var <- paste0("co_items_mean_", sfx) dist_dvs <- c( `Character traits` = paste0("char_dist_", sfx), `Political distance` = paste0("pol_dist_", sfx), `Social distance` = paste0("priv_dist_", sfx) ) for (dv_name in names(dist_dvs)) { dv_var <- dist_dvs[dv_name] if (!dv_var %in% names(dat)) next f_ap <- as.formula(paste(dv_var, "~", ap_var, "+", controls_rhs)) m_ap <- lm(f_ap, data = dat) f_ic <- as.formula(paste(dv_var, "~", id_var, "+", co_var, "+", controls_rhs)) m_ic <- lm(f_ic, data = dat) cr_row <- function(predictor, model, var) { e <- extract_hc3(model, var) tibble(issue = issue_name, outcome = dv_name, predictor = predictor, b = e$b, se = e$se, p = e$p, N = nobs(model), adj_R2 = summary(model)$adj.r.squared) } crit_coefs[[paste(issue_name, dv_name)]] <- bind_rows( cr_row("APSIT", m_ap, ap_var), cr_row("Identity", m_ic, id_var), cr_row("Comparison", m_ic, co_var) ) } } crit_coefs_df <- bind_rows(crit_coefs) # ══════════════════════════════════════════════════════════════════════════════ # 3. PART III: DEMOCRATIC ATTITUDES # ══════════════════════════════════════════════════════════════════════════════ # --- 3.1 Incremental Validity (Tables G, H; Figure 2) ------------------------ outcomes_shared <- list( "Democratic support" = "demo_support_mean", "Democracy evaluation" = "democracyevaluationm", "Democracy satisfaction" = "democracy_satisfaction", "Institutional trust" = "trust_mean" ) inc_results <- list() fig2_coefs <- list() full_model_rows <- list() # full coefficient tables for Appendix Tables H1-H4 # Camp-control variant: identical models plus camp membership (Supporter vs. # Opponent) as a covariate. The Figure 2/3 notes in the current manuscript # claim this control, the base models do not include it. Both variants are # estimated and exported so the choice can be made from the numbers; note that # adding camp drops neutral respondents (group_ is NA for them), so Ns differ. inc_results_camp <- list() full_model_rows_camp <- list() sep_model_rows <- list() # Appendix Tables I1-I4 (separate models) tidy_hc3 <- function(model, model_label, issue, dv) { ct <- lmtest::coeftest(model, vcov. = sandwich::vcovHC(model, type = "HC3")) tibble( Domain = issue, Outcome = dv, Model = model_label, term = rownames(ct), b = ct[, 1], se = ct[, 2], p = ct[, 4], N = nobs(model), adj_R2 = summary(model)$adj.r.squared ) } for (issue_name in names(issues)) { sfx <- tolower(issue_name) cfg <- issues[[issue_name]] # Use the manuscript-consistent subscale means built in section 0.8 id_raw <- paste0("identity_", sfx, "_mean") co_raw <- paste0("comparison_", sfx, "_mean") ft_raw <- paste0("FT_", sfx) tin_raw <- paste0("Tin_", sfx) tout_raw <- paste0("Tout_", sfx) ap_raw <- paste0("APSIT_", sfx) outcomes_all <- c( outcomes_shared, list("Political tolerance" = paste0("poltolerance_", sfx)) ) for (dv_name in names(outcomes_all)) { dv_var <- outcomes_all[[dv_name]] if (!dv_var %in% names(dat)) { cat("WARNING: Missing DV", dv_var, "\n"); next } needed <- c(ap_raw, id_raw, co_raw, ft_raw, tin_raw, tout_raw, dv_var) if (any(!needed %in% names(dat))) { cat("WARNING: Missing vars for", issue_name, dv_name, "\n"); next } dv_term <- paste0("scale(", dv_var, ")") d <- dat |> mutate( AP_z = as.numeric(scale(.data[[ap_raw]])), I_z = as.numeric(scale(.data[[id_raw]])), C_z = as.numeric(scale(.data[[co_raw]])), FT_z = as.numeric(scale(.data[[ft_raw]])), Tin_z = as.numeric(scale(.data[[tin_raw]])), Tout_z = as.numeric(scale(.data[[tout_raw]])), camp = .data[[paste0("group_", sfx)]] ) f4a <- as.formula(paste(dv_term, "~ AP_z +", controls_rhs)) f4b <- as.formula(paste(dv_term, "~ FT_z +", controls_rhs)) f1 <- as.formula(paste(dv_term, "~ AP_z + FT_z +", controls_rhs)) f2 <- as.formula(paste(dv_term, "~ I_z + C_z + FT_z +", controls_rhs)) f3 <- as.formula(paste(dv_term, "~ I_z + C_z + Tin_z + Tout_z +", controls_rhs)) m4a <- lm(f4a, data = d); m4b <- lm(f4b, data = d) m1 <- lm(f1, data = d); m2 <- lm(f2, data = d); m3 <- lm(f3, data = d) # Camp-control variant (adds camp; drops neutrals via listwise deletion) f1c <- update(f1, . ~ . + camp) f2c <- update(f2, . ~ . + camp) f3c <- update(f3, . ~ . + camp) m1c <- lm(f1c, data = d); m2c <- lm(f2c, data = d); m3c <- lm(f3c, data = d) r1_ap <- extract_hc3(m1, "AP_z"); r1_ft <- extract_hc3(m1, "FT_z") r2_i <- extract_hc3(m2, "I_z"); r2_c <- extract_hc3(m2, "C_z") r2_ft <- extract_hc3(m2, "FT_z") r3_i <- extract_hc3(m3, "I_z"); r3_c <- extract_hc3(m3, "C_z") r3_tin <- extract_hc3(m3, "Tin_z"); r3_tout <- extract_hc3(m3, "Tout_z") r4a_ap <- extract_hc3(m4a, "AP_z"); r4b_ft <- extract_hc3(m4b, "FT_z") # Full coefficient tables (Appendix Tables J1-J4) full_model_rows[[paste(issue_name, dv_name)]] <- bind_rows( tidy_hc3(m1, "Model 1", issue_name, dv_name), tidy_hc3(m2, "Model 2", issue_name, dv_name), tidy_hc3(m3, "Model 3", issue_name, dv_name) ) full_model_rows_camp[[paste(issue_name, dv_name)]] <- bind_rows( tidy_hc3(m1c, "Model 1 + camp", issue_name, dv_name), tidy_hc3(m2c, "Model 2 + camp", issue_name, dv_name), tidy_hc3(m3c, "Model 3 + camp", issue_name, dv_name) ) # Appendix Tables I1-I4 (separate models: APSIT only vs. FT only). # The manuscript prints beta AND adjusted R-squared per model, so both # single-predictor models are stored with their own fit statistics. sep_model_rows[[paste(issue_name, dv_name)]] <- tibble( Domain = issue_name, Outcome = dv_name, M4a_APSIT_b = r4a_ap$b, M4a_APSIT_se = r4a_ap$se, M4a_APSIT_p = r4a_ap$p, M4a_AdjR2 = round(summary(m4a)$adj.r.squared, 3), M4a_N = nobs(m4a), M4b_FT_b = r4b_ft$b, M4b_FT_se = r4b_ft$se, M4b_FT_p = r4b_ft$p, M4b_AdjR2 = round(summary(m4b)$adj.r.squared, 3), M4b_N = nobs(m4b) ) inc_results[[paste(issue_name, dv_name)]] <- tibble( Domain = issue_name, Outcome = dv_name, M4a_APSIT_b = r4a_ap$b, M4b_FT_b = r4b_ft$b, M1_APSIT_b = r1_ap$b, M1_FT_b = r1_ft$b, M2_I_b = r2_i$b, M2_C_b = r2_c$b, M2_FT_b = r2_ft$b, N = nobs(m3), AdjR2_M1 = round(summary(m1)$adj.r.squared, 3), AdjR2_M2 = round(summary(m2)$adj.r.squared, 3), AdjR2_M3 = round(summary(m3)$adj.r.squared, 3), DeltaR2 = round(summary(m2)$adj.r.squared - summary(m1)$adj.r.squared, 3), DeltaR2_M3M2 = round(summary(m3)$adj.r.squared - summary(m2)$adj.r.squared, 3) ) inc_results_camp[[paste(issue_name, dv_name)]] <- tibble( Domain = issue_name, Outcome = dv_name, M2c_I_b = extract_hc3(m2c, "I_z")$b, M2c_C_b = extract_hc3(m2c, "C_z")$b, M2c_FT_b = extract_hc3(m2c, "FT_z")$b, M2c_camp_b = extract_hc3(m2c, "campOpponent")$b, N_camp = nobs(m2c), AdjR2_M1c = round(summary(m1c)$adj.r.squared, 3), AdjR2_M2c = round(summary(m2c)$adj.r.squared, 3), AdjR2_M3c = round(summary(m3c)$adj.r.squared, 3) ) fig2_coefs[[paste(issue_name, dv_name)]] <- bind_rows( tibble(issue=issue_name, outcome=dv_name, model="Model 1", predictor="APSIT", coef=r1_ap$b, se=r1_ap$se), tibble(issue=issue_name, outcome=dv_name, model="Model 1", predictor="Feeling Thermometer", coef=r1_ft$b, se=r1_ft$se), tibble(issue=issue_name, outcome=dv_name, model="Model 2", predictor="Identity", coef=r2_i$b, se=r2_i$se), tibble(issue=issue_name, outcome=dv_name, model="Model 2", predictor="Comparison", coef=r2_c$b, se=r2_c$se), tibble(issue=issue_name, outcome=dv_name, model="Model 2", predictor="Feeling Thermometer", coef=r2_ft$b, se=r2_ft$se), tibble(issue=issue_name, outcome=dv_name, model="Model 3", predictor="Identity", coef=r3_i$b, se=r3_i$se), tibble(issue=issue_name, outcome=dv_name, model="Model 3", predictor="Comparison", coef=r3_c$b, se=r3_c$se), tibble(issue=issue_name, outcome=dv_name, model="Model 3", predictor="Ingroup Temp.", coef=r3_tin$b, se=r3_tin$se), tibble(issue=issue_name, outcome=dv_name, model="Model 3", predictor="Outgroup Temp.", coef=r3_tout$b, se=r3_tout$se) ) } } inc_table <- bind_rows(inc_results) inc_table_camp <- bind_rows(inc_results_camp) fig2_df <- bind_rows(fig2_coefs) cat("\n=== Tables J1-J4: Incremental Validity (Model 2 key coefficients) ===\n") print(inc_table |> select(Domain, Outcome, M1_APSIT_b, M1_FT_b, M2_I_b, M2_C_b, M2_FT_b, N, AdjR2_M1, AdjR2_M2, AdjR2_M3, DeltaR2) |> mutate(across(where(is.numeric), ~round(.x, 3)))) sep_table <- bind_rows(sep_model_rows) cat("\n=== Tables I1-I4: Separate Models (APSIT only vs. FT only) ===\n") print(sep_table |> select(Domain, Outcome, M4a_APSIT_b, M4a_AdjR2, M4b_FT_b, M4b_AdjR2, M4a_N) |> mutate(across(where(is.numeric), ~round(.x, 3))) |> as.data.frame()) write.csv(sep_table, out("TablesI1_I4_separate_models.csv"), row.names = FALSE) cat("\n=== Tables J1-J4 (camp-control variant): Model 2 + camp ===\n") print(inc_table_camp |> mutate(across(where(is.numeric), ~round(.x, 3)))) write.csv(inc_table_camp, out("TablesJ1_J4_summary_campcontrol.csv"), row.names = FALSE) write.csv(bind_rows(full_model_rows_camp), out("TablesJ1_J4_full_models_campcontrol.csv"), row.names = FALSE) # --- 3.2 Differential Effects by Group (Figure 3) ---------------------------- fig3_coefs <- list() fig3_coefs_camp <- list() # variant with a camp main effect (separate intercepts) fig3_slopediff <- list() # Brambor-style interaction parametrization: # camp x dimension coefficients = formal tests of # the supporter-opponent slope difference for (issue_name in names(issues)) { sfx <- tolower(issue_name) preds <- paste0(c("id_sup_", "id_opp_", "co_sup_", "co_opp_"), sfx) if (any(!preds %in% names(dat))) next outcomes_fig3 <- c( outcomes_shared, list("Political tolerance" = paste0("poltolerance_", sfx)) ) for (dv_name in names(outcomes_fig3)) { dv_var <- outcomes_fig3[[dv_name]] if (!dv_var %in% names(dat)) next d3 <- dat d3$camp <- dat[[paste0("group_", sfx)]] f <- as.formula(paste(dv_var, "~", paste(preds, collapse = " + "), "+", controls_rhs)) m <- lm(f, data = d3) # Camp main effect added: without it, both camps share one intercept, # which the current split-slope design silently imposes. f_c <- update(f, . ~ . + camp) m_c <- lm(f_c, data = d3) # Interaction parametrization (Brambor, Clark, and Golder 2006): # DV ~ camp + I_z + C_z + camp:I_z + camp:C_z + controls. This spans the # SAME column space as the split-slope + camp model above (identical fit), # but the camp:dimension coefficients directly estimate the opponent- # minus-supporter slope DIFFERENCE with an HC3 standard error - the formal # test behind the manuscript claim that comparison effects are "more # consistent and stronger for opponents". d3$I_zf <- dat[[paste0("identity_", sfx, "_z")]] d3$C_zf <- dat[[paste0("comparison_", sfx, "_z")]] f_i <- as.formula(paste(dv_var, "~ camp*I_zf + camp*C_zf +", controls_rhs)) m_i <- lm(f_i, data = d3) for (tm in c("campOpponent:I_zf", "campOpponent:C_zf")) { e <- extract_hc3(m_i, tm) fig3_slopediff[[paste(issue_name, dv_name, tm)]] <- tibble( issue = issue_name, outcome = dv_name, dimension = ifelse(grepl("I_zf", tm), "Identity", "Comparison"), diff_opp_minus_sup = e$b, se = e$se, p = e$p, N = nobs(m_i)) } for (pred in preds) { r <- extract_hc3(m, pred) rc <- extract_hc3(m_c, pred) label <- dplyr::case_when( grepl("id_sup", pred) ~ "Identity x Supporter", grepl("id_opp", pred) ~ "Identity x Opponent", grepl("co_sup", pred) ~ "Comparison x Supporter", grepl("co_opp", pred) ~ "Comparison x Opponent" ) fig3_coefs[[paste(issue_name, dv_name, pred)]] <- tibble( issue = issue_name, outcome = dv_name, predictor = label, b = r$b, se = r$se ) fig3_coefs_camp[[paste(issue_name, dv_name, pred)]] <- tibble( issue = issue_name, outcome = dv_name, predictor = label, b = rc$b, se = rc$se, camp_main_b = extract_hc3(m_c, "campOpponent")$b ) } } } fig3_df <- bind_rows(fig3_coefs) fig3_df_camp <- bind_rows(fig3_coefs_camp) write.csv(fig3_df_camp, out("Figure3_coefficients_campmain.csv"), row.names = FALSE) fig3_slopediff_df <- bind_rows(fig3_slopediff) cat("\n=== Figure 3: opponent-minus-supporter slope differences (HC3) ===\n") print(fig3_slopediff_df |> mutate(across(where(is.numeric), ~round(.x, 3))) |> as.data.frame()) write.csv(fig3_slopediff_df, out("TableJ5_slope_differences.csv"), row.names = FALSE) # ══════════════════════════════════════════════════════════════════════════════ # 4. FIGURES # ══════════════════════════════════════════════════════════════════════════════ theme_campos <- function(base_size = 10) { theme_minimal(base_size = base_size, base_family = "sans") + theme( panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(color = "grey85", linewidth = 0.3), panel.background = element_rect(fill = "grey97", color = NA), panel.spacing = unit(0.8, "lines"), strip.background = element_rect(fill = "grey85", color = NA), strip.text = element_text(face = "bold", size = 8, color = "grey20"), axis.text.y = element_text(size = 7, color = "grey30"), axis.text.x = element_text(size = 6, color = "grey30"), axis.title.x = element_text(size = 7, color = "grey20", margin = margin(t = 5)), axis.title.y = element_blank(), legend.position = "bottom", legend.title = element_blank(), legend.text = element_text(size = 7), legend.key.size = unit(0.7, "lines"), legend.margin = margin(t = 5, b = 0), plot.title = element_text(face = "bold", size = 9, hjust = 0.5), plot.subtitle = element_text(size = 7, hjust = 0.5, color = "grey40"), plot.margin = margin(8, 8, 4, 8) ) } issue_colors <- c("Migration" = "#D55E00", "Climate" = "#009E73", "Vaccination" = "#E69F00", "Gender" = "#0072B2") # --- FIGURE 1: Combined Convergent & Criterion Validity ---------------------- fig1_combined_data <- bind_rows( conv_coefs_df, crit_coefs_df ) |> mutate( # NO sign flipping: coefficients are shown as estimated. The orientation # of the pre-aggregated ap1_* outcome is established empirically by the # "AP1 ORIENTATION CHECK" console block in Part II; align the Figure 1 # note and the Table F1/F2 titles in the manuscript with that result. ci_low = b - 1.96 * se, ci_high = b + 1.96 * se, issue = factor(issue, levels = c("Migration", "Climate", "Vaccination", "Gender")), outcome = factor(outcome, levels = c( "Self-placement", "Feeling Thermometer", "Character traits", "Political distance", "Social distance" )), predictor = factor(predictor, levels = c("Comparison", "Identity", "APSIT")) ) fig1 <- ggplot( fig1_combined_data, aes(x = b, y = predictor, xmin = ci_low, xmax = ci_high, colour = issue) ) + geom_vline(xintercept = 0, linetype = "dashed", colour = "grey50", linewidth = 0.4) + geom_errorbarh( aes(xmin = ci_low, xmax = ci_high), height = 0, linewidth = 0.35, position = position_dodge(width = 0.65) ) + geom_point( size = 1.2, position = position_dodge(width = 0.65) ) + facet_wrap(~ outcome, nrow = 1, scales = "free_x") + scale_colour_manual(values = issue_colors) + scale_x_continuous(breaks = scales::pretty_breaks(n = 3)) + labs( x = "Regression Coefficient (b, 95% CI)", y = NULL ) + theme_campos() + theme( legend.position = "bottom", legend.direction = "horizontal", strip.text = element_text(face = "bold", size = 8), panel.spacing.x = unit(0.6, "lines"), axis.text.y = element_text(size = 8) ) + guides(colour = guide_legend(nrow = 1)) ggsave(out("Figure1_validation_combined.pdf"), fig1, width = 11, height = 3.5, dpi = 300) ggsave(out("Figure1_validation_combined.png"), fig1, width = 11, height = 3.5, dpi = 300, bg = "white") cat("\n=== Figure 1 (Combined Validation) saved ===\n") # --- FIGURE 2: Incremental Validity ------------------------------------------ outcome_levels <- c("Democratic Support", "Democracy Evaluation", "Democracy Satisfaction", "Institutional Trust", "Political Tolerance") fig2_plot <- fig2_df |> mutate(ci_low = coef - 1.96*se, ci_high = coef + 1.96*se, issue = factor(issue, levels = names(issues)), outcome = dplyr::recode(outcome, "Democratic support" = "Democratic Support", "Democracy evaluation" = "Democracy Evaluation", "Democracy satisfaction" = "Democracy Satisfaction", "Institutional trust" = "Institutional Trust", "Political tolerance" = "Political Tolerance"), outcome = factor(outcome, levels = outcome_levels)) create_fig2_panel <- function(data, model_name, subtitle, pred_levels, show_legend = FALSE) { d <- data |> filter(model == model_name) |> mutate(predictor = factor(predictor, levels = rev(pred_levels))) p <- ggplot(d, aes(x = coef, y = predictor, color = issue)) + geom_vline(xintercept = 0, linetype = "dashed", color = "grey50", linewidth = 0.4) + geom_errorbarh(aes(xmin = ci_low, xmax = ci_high), height = 0, linewidth = 0.35, position = position_dodge(width = 0.65)) + geom_point(size = 1.1, position = position_dodge(width = 0.65)) + facet_wrap(~ outcome, ncol = 5) + scale_color_manual(values = issue_colors) + scale_x_continuous(breaks = seq(-0.6, 0.4, 0.2), labels = function(x) sprintf("%.1f", x)) + coord_cartesian(xlim = c(-0.75, 0.5)) + # zoom without dropping points labs(title = model_name, subtitle = subtitle, x = "Standardized Coefficient") + theme_campos() if (!show_legend) p <- p + theme(legend.position = "none") p } p2_m1 <- create_fig2_panel(fig2_plot, "Model 1", "APSIT Index + Feeling Thermometer Index", c("APSIT", "Feeling Thermometer")) p2_m2 <- create_fig2_panel(fig2_plot, "Model 2", "Identity + Comparison + Feeling Thermometer Index", c("Identity", "Comparison", "Feeling Thermometer")) p2_m3 <- create_fig2_panel(fig2_plot, "Model 3", "Identity + Comparison + Ingroup Temp. + Outgroup Temp.", c("Identity", "Comparison", "Ingroup Temp.", "Outgroup Temp."), show_legend = TRUE) # No embedded figure title: APSR sets figure numbers/captions from the # manuscript (previous TODO resolved). fig2 <- p2_m1 / p2_m2 / p2_m3 + plot_layout(heights = c(1, 1.2, 1.5)) ggsave(out("Figure2_incremental_validity.pdf"), fig2, width = 13, height = 10) ggsave(out("Figure2_incremental_validity.png"), fig2, width = 13, height = 10, dpi = 300, bg = "white") cat("\n=== Figure 2 (Incremental Validity) saved ===\n") # --- FIGURE 3: Differential Effects by Group --------------------------------- # IMPORTANT: Figure 3 in the manuscript is plotted from the camp-main-effect # specification (Brambor, Clark, and Golder 2006), NOT from the split-slope # model without a camp intercept. fig3_df_camp is therefore the plot source; # fig3_df (no camp main effect) is exported for comparison only. fig3_plot <- fig3_df_camp |> mutate( ci_low = b - 1.96 * se, ci_high = b + 1.96 * se, issue = factor(issue, levels = c("Migration", "Climate", "Vaccination", "Gender")), outcome = dplyr::recode(outcome, "Democratic support" = "Democratic Support", "Democracy evaluation" = "Democracy Evaluation", "Democracy satisfaction" = "Democracy Satisfaction", "Institutional trust" = "Institutional Trust", "Political tolerance" = "Political Tolerance"), outcome = factor(outcome, levels = c( "Democratic Support", "Democracy Evaluation", "Democracy Satisfaction", "Institutional Trust", "Political Tolerance")), predictor = factor(predictor, levels = rev(c( "Identity x Supporter", "Identity x Opponent", "Comparison x Supporter", "Comparison x Opponent"))) ) fig3 <- ggplot( fig3_plot, aes(x = b, y = predictor, xmin = ci_low, xmax = ci_high, colour = issue) ) + geom_vline(xintercept = 0, linetype = "dashed", colour = "grey50", linewidth = 0.4) + geom_errorbarh( aes(xmin = ci_low, xmax = ci_high), height = 0, linewidth = 0.35, position = position_dodge(width = 0.65) ) + geom_point( size = 1.1, position = position_dodge(width = 0.65) ) + facet_wrap(~ outcome, nrow = 1, scales = "free_x") + scale_colour_manual(values = issue_colors) + scale_x_continuous(breaks = scales::pretty_breaks(n = 3)) + labs( x = "Regression Coefficient (b, 95% CI)", y = NULL ) + theme_campos() + theme( legend.position = "bottom", legend.direction = "horizontal", strip.text = element_text(face = "bold", size = 8), panel.spacing.x = unit(0.6, "lines"), axis.text.y = element_text(size = 8) ) + guides(colour = guide_legend(nrow = 1)) ggsave(out("Figure3_group_effects.pdf"), fig3, width = 13, height = 4, dpi = 300) ggsave(out("Figure3_group_effects.png"), fig3, width = 13, height = 4, dpi = 300, bg = "white") cat("\n=== Figure 3 (Group Effects) saved ===\n") # --- FIGURE D1 (Appendix): APSIT Distribution -------------------------------- apsit_dist <- bind_rows(lapply(names(issues), function(issue_name) { sfx <- tolower(issue_name) tibble( issue = issue_name, APSIT = dat[[paste0("APSIT_", sfx)]] ) })) figD1 <- ggplot(apsit_dist, aes(x = APSIT)) + geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "grey70", color = "white", alpha = 0.6) + geom_density(color = "black", linewidth = 0.8) + facet_wrap(~issue, ncol = 2) + labs( x = "APSIT Index (standardized)", y = "Density" ) + theme_campos() ggsave(out("FigureD1_APSIT_distribution.pdf"), figD1, width = 8, height = 6) ggsave(out("FigureD1_APSIT_distribution.png"), figD1, width = 8, height = 6, dpi = 300, bg = "white") cat("\n=== Figure D1 (APSIT Distribution) saved ===\n") cat("\nAll figures saved.\n") # ══════════════════════════════════════════════════════════════════════════════ # 5. APPENDIX # ══════════════════════════════════════════════════════════════════════════════ # --- Table G1: Sample Descriptives ------------------------------- # Previous version was wrong in three ways: Female/Age were computed on the # FULL sample for every domain, neutrals were excluded from N, and the # supporter/opponent percentages had no neutral base. New version: # - Domain membership: respondent answered the domain's self-placement # (between-subjects design; expected Ns ~ 1,012/998/1,003/988). # - Camps from self-placement per Appendix A3: 0-4 Supporter, 5 Neutral, # 6-10 Opponent (0 = strong SUPPORTER; the main text sentence claiming # 0 = rejection is the error, not A3). # - Descriptives computed WITHIN each domain; language shares added. sample_desc <- bind_rows(lapply(names(issues), function(nm) { cfg <- issues[[nm]] pos <- suppressWarnings(as.numeric(dat[[cfg$pos]])) inD <- !is.na(pos) cat3 <- cut(pos, breaks = c(-Inf, 4.5, 5.5, Inf), labels = c("Supporter", "Neutral", "Opponent")) # Cross-check the A3-based camps against the opp-based group_ variable: grp <- dat[[paste0("group_", tolower(nm))]] mismatch <- sum(inD & !is.na(grp) & ((cat3 == "Supporter" & grp != "Supporter") | (cat3 == "Opponent" & grp != "Opponent")), na.rm = TRUE) tibble( Domain = nm, N = sum(inD), Female_pct = round(100 * mean(dat$geschlecht[inD] == 0, na.rm = TRUE), 1), Age_mean = round(mean(dat$alter[inD], na.rm = TRUE), 1), Age_sd = round(sd(dat$alter[inD], na.rm = TRUE), 1), Education_mean = round(mean(dat$bildung[inD], na.rm = TRUE), 2), German_pct = round(100 * mean(dat$sprachregion[inD] == "German", na.rm = TRUE), 1), French_pct = round(100 * mean(dat$sprachregion[inD] == "French", na.rm = TRUE), 1), Italian_pct = round(100 * mean(dat$sprachregion[inD] == "Italian", na.rm = TRUE), 1), Supporters_pct = round(100 * mean(cat3[inD] == "Supporter", na.rm = TRUE), 1), Neutral_pct = round(100 * mean(cat3[inD] == "Neutral", na.rm = TRUE), 1), Opponents_pct = round(100 * mean(cat3[inD] == "Opponent", na.rm = TRUE), 1), Camp_mismatch_vs_opp = mismatch ) })) cat("\n=== Table G1: Sample Descriptives (per domain, incl. neutrals) ===\n") print(as.data.frame(sample_desc)) cat("Note: Camp_mismatch_vs_opp > 0 would indicate that the opp-based group_\n", "variable and the A3 self-placement rule disagree; investigate if so.\n") # --- Parallel Analysis Plots (Appendix Figures K1-K4) ------------------------ for (issue_name in names(issues)) { cfg <- issues[[issue_name]] sfx <- tolower(issue_name) for (dim_label in c("Identity", "Comparison", "Combined")) { items <- if (dim_label == "Identity") cfg$ident else if (dim_label == "Comparison") cfg$comp else c(cfg$ident, cfg$comp) X <- dat[, items, drop = FALSE] X[] <- lapply(X, function(x) suppressWarnings(as.numeric(x))) R <- psych::polychoric(X)$rho nobs <- nrow(na.omit(X)) fname <- out(paste0("FigK_parallel_", sfx, "_", tolower(dim_label), ".png")) png(fname, width = 600, height = 400) psych::fa.parallel(R, fm = "pa", fa = "fa", n.obs = nobs, main = paste("Parallel Analysis -", dim_label, "-", issue_name)) dev.off() } } # ── Index builder: three aggregation rules from raw subscale means ─────────── agg_indices <- function(id_mean, co_mean) { ok <- complete.cases(id_mean, co_mean) # Published additive index: arithmetic mean of the z-standardized subscales add <- rep(NA_real_, length(id_mean)) add[ok] <- rowMeans(cbind(scale(id_mean)[ok], scale(co_mean)[ok])) list( additive = as.numeric(scale(add)), geometric = as.numeric(scale(sqrt(id_mean * co_mean))), minimum = as.numeric(scale(pmin(id_mean, co_mean))) ) } rules <- c("additive", "geometric", "minimum") # ── Re-estimate Figure 1 and Figure 2 models under each rule ───────────────── robust_long <- list() cor_long <- list() for (issue_name in names(issues)) { cfg <- issues[[issue_name]] sfx <- tolower(issue_name) # Subscale means: min_valid = 3 throughout (harmonized; see Part II). id1 <- row_mean_min(dat, cfg$ident, min_valid = 3) co1 <- row_mean_min(dat, cfg$comp, min_valid = 3) id3 <- dat[[paste0("identity_", sfx, "_mean")]] co3 <- dat[[paste0("comparison_", sfx, "_mean")]] idx1 <- agg_indices(id1, co1) # for Figure 1 idx3 <- agg_indices(id3, co3) # for Figure 2 # Correlations among the index variants (reported in the appendix note) idxmat <- as.data.frame(idx3) cor_long[[issue_name]] <- tibble( issue = issue_name, r_add_geom = cor(idxmat$additive, idxmat$geometric, use = "complete.obs"), r_add_min = cor(idxmat$additive, idxmat$minimum, use = "complete.obs"), r_geom_min = cor(idxmat$geometric, idxmat$minimum, use = "complete.obs") ) ft_z <- as.numeric(scale(dat[[paste0("FT_", sfx)]])) # FT control for Figure 2 # ----- Figure 1: convergent + criterion (raw outcome, standardized APSIT) -- f1_dvs <- c( `Feeling Thermometer` = paste0("ap1_", sfx), `Self-placement` = paste0("pos_extreme_", sfx), `Character traits` = paste0("char_dist_", sfx), `Political distance` = paste0("pol_dist_", sfx), `Social distance` = paste0("priv_dist_", sfx) ) for (dv_name in names(f1_dvs)) { dv <- f1_dvs[[dv_name]] if (!dv %in% names(dat)) next for (rule in rules) { d <- dat; d$IDX <- idx1[[rule]] f <- as.formula(paste(dv, "~ IDX +", controls_rhs)) r <- extract_hc3(lm(f, data = d), "IDX") robust_long[[paste("F1", issue_name, dv_name, rule)]] <- tibble( figure = "Figure 1", issue = issue_name, outcome = dv_name, rule = rule, b = r$b, se = r$se, p = r$p) } } # ----- Figure 2: incremental, Model 1 (z outcome, standardized APSIT + FT_z) f2_dvs <- c( `Democratic support` = "demo_support_mean", `Democracy evaluation` = "democracyevaluationm", `Democracy satisfaction` = "democracy_satisfaction", `Institutional trust` = "trust_mean", `Political tolerance` = paste0("poltolerance_", sfx) ) for (dv_name in names(f2_dvs)) { dv <- f2_dvs[[dv_name]] if (!dv %in% names(dat)) next for (rule in rules) { d <- dat; d$IDX <- idx3[[rule]]; d$FT_z <- ft_z f <- as.formula(paste0("scale(", dv, ") ~ IDX + FT_z + ", controls_rhs)) r <- extract_hc3(lm(f, data = d), "IDX") robust_long[[paste("F2", issue_name, dv_name, rule)]] <- tibble( figure = "Figure 2", issue = issue_name, outcome = dv_name, rule = rule, b = r$b, se = r$se, p = r$p) } } } robust_df <- bind_rows(robust_long) cor_df <- bind_rows(cor_long) # ── Wide table for the appendix (b with significance stars) ────────────────── star <- function(p) ifelse(p < .001, "***", ifelse(p < .01, "**", ifelse(p < .05, "*", ""))) robust_wide <- robust_df |> mutate(cell = sprintf("%+.3f%s", b, star(p))) |> select(figure, issue, outcome, rule, cell) |> tidyr::pivot_wider(names_from = rule, values_from = cell) |> arrange(figure, issue, outcome) cat("\n=== Aggregation robustness: index-variant correlations (min_valid = 3) ===\n") print(cor_df |> mutate(across(where(is.numeric), ~round(.x, 3)))) cat("\n=== Aggregation robustness: APSIT coefficients by rule ===\n") print(robust_wide, n = nrow(robust_wide)) # ── Agreement summary (geometric / minimum vs. additive baseline) ──────────── agree <- robust_df |> select(figure, issue, outcome, rule, b, p) |> tidyr::pivot_wider(names_from = rule, values_from = c(b, p)) summary_tbl <- agree |> group_by(figure) |> summarise( n_coef = dplyr::n(), sign_geom = mean(sign(b_geometric) == sign(b_additive)), sign_min = mean(sign(b_minimum) == sign(b_additive)), sig_geom = mean((p_geometric < .05) == (p_additive < .05)), sig_min = mean((p_minimum < .05) == (p_additive < .05)), .groups = "drop" ) cat("\n=== Sign / significance agreement vs. additive baseline ===\n") print(summary_tbl |> mutate(across(where(is.numeric), ~round(.x, 3)))) write.csv(robust_df, out("TableX_aggregation_robustness_long.csv"), row.names = FALSE) write.csv(robust_wide, out("TableX_aggregation_robustness.csv"), row.names = FALSE) write.csv(cor_df, out("TableX_aggregation_index_correlations.csv"), row.names = FALSE) # ── Appendix figures: APSIT coefficient under the three rules ──────────────── rule_colors <- c(additive = "#000000", geometric = "#0072B2", minimum = "#D55E00") plot_df <- robust_df |> mutate( ci_low = b - 1.96 * se, ci_high = b + 1.96 * se, rule = factor(rule, levels = c("additive", "geometric", "minimum")), issue = factor(issue, levels = c("Migration", "Climate", "Vaccination", "Gender")) ) make_robust_fig <- function(fig_label, dv_levels, title) { d <- plot_df |> filter(figure == fig_label) |> mutate(outcome = factor(outcome, levels = dv_levels)) ggplot(d, aes(x = b, y = issue, colour = rule)) + geom_vline(xintercept = 0, linetype = "dashed", colour = "grey50", linewidth = 0.4) + geom_errorbarh(aes(xmin = ci_low, xmax = ci_high), height = 0, linewidth = 0.35, position = position_dodge(width = 0.6)) + geom_point(size = 1.2, position = position_dodge(width = 0.6)) + facet_wrap(~ outcome, nrow = 1, scales = "free_x") + scale_colour_manual(values = rule_colors) + scale_x_continuous(breaks = scales::pretty_breaks(n = 3)) + labs(x = "APSIT coefficient (b, 95% CI)", y = NULL, title = title) + theme_campos() + theme(legend.position = "bottom", strip.text = element_text(face = "bold", size = 8)) } # Appendix labels D2/D3 (previous "Figure X1/X2" TODO resolved); no embedded # titles, per APSR figure policy. figD2 <- make_robust_fig( "Figure 1", c("Self-placement", "Feeling Thermometer", "Character traits", "Political distance", "Social distance"), NULL) figD3 <- make_robust_fig( "Figure 2", c("Democratic support", "Democracy evaluation", "Democracy satisfaction", "Institutional trust", "Political tolerance"), NULL) ggsave(out("FigureD2_aggregation_robustness_validity.pdf"), figD2, width = 13, height = 4) ggsave(out("FigureD2_aggregation_robustness_validity.png"), figD2, width = 13, height = 4, dpi = 300, bg = "white") ggsave(out("FigureD3_aggregation_robustness_incremental.pdf"), figD3, width = 13, height = 4) ggsave(out("FigureD3_aggregation_robustness_incremental.png"), figD3, width = 13, height = 4, dpi = 300, bg = "white") cat("\n=== Aggregation-robustness tables and figures saved to:", output_dir, "===\n") ### End of aggregation-robustness block ### # ══════════════════════════════════════════════════════════════════════════════ # 5b. APPENDIX J: ROBUSTNESS CHECKS (previously claimed but not coded) # ══════════════════════════════════════════════════════════════════════════════ # Manuscript Appendix J claims robustness to (i) ordered logistic regression, # (ii) multilevel models clustering by language region, (iii) exclusion of # neutral respondents (self-placement = 5), (iv) separate estimation by # language region. This section produces those numbers for the core Model 2 # specification (I_z + C_z + FT_z + controls). If any check fails to support # the claim, soften or drop the corresponding sentence in Appendix J. controls_rhs_noreg <- "geschlecht + alter + bildung + einkommenssituation + rile_z" build_d_j <- function(cfg, sfx) { dat |> mutate( I_z = as.numeric(scale(.data[[paste0("identity_", sfx, "_mean")]])), C_z = as.numeric(scale(.data[[paste0("comparison_", sfx, "_mean")]])), FT_z = as.numeric(scale(.data[[paste0("FT_", sfx)]])), pos_raw = suppressWarnings(as.numeric(.data[[cfg$pos]])) ) } rob_rows <- list() add_rob <- function(key, issue, dv, spec, term, b, se, p, note = "") { rob_rows[[key]] <<- tibble(issue = issue, outcome = dv, spec = spec, term = term, b = b, se = se, p = p, note = note) } for (issue_name in names(issues)) { cfg <- issues[[issue_name]] sfx <- tolower(issue_name) d <- build_d_j(cfg, sfx) outcomes_j <- c(outcomes_shared, list("Political tolerance" = paste0("poltolerance_", sfx))) for (dv_name in names(outcomes_j)) { dv_var <- outcomes_j[[dv_name]] if (!dv_var %in% names(dat)) next dv_term <- paste0("scale(", dv_var, ")") # (0) base: identical to Tables H Model 2 (reference row) m_base <- lm(as.formula(paste(dv_term, "~ I_z + C_z + FT_z +", controls_rhs)), data = d) for (tm in c("I_z", "C_z")) { e <- extract_hc3(m_base, tm) add_rob(paste("base", issue_name, dv_name, tm), issue_name, dv_name, "base_M2", ifelse(tm == "I_z", "Identity", "Comparison"), e$b, e$se, e$p) } # (iii) exclusion of neutral respondents (self-placement = 5) m_non <- lm(as.formula(paste(dv_term, "~ I_z + C_z + FT_z +", controls_rhs)), data = d |> filter(!is.na(pos_raw) & pos_raw != 5)) for (tm in c("I_z", "C_z")) { e <- extract_hc3(m_non, tm) add_rob(paste("non", issue_name, dv_name, tm), issue_name, dv_name, "no_neutrals", ifelse(tm == "I_z", "Identity", "Comparison"), e$b, e$se, e$p) } # (i) ordered logit on the raw outcome (rounded to half scale points) d_ol <- d d_ol$.y_ord <- ordered(round(suppressWarnings(as.numeric(d_ol[[dv_var]])) * 2) / 2) m_ol <- tryCatch( MASS::polr(as.formula(paste(".y_ord ~ I_z + C_z + FT_z +", controls_rhs)), data = d_ol, Hess = TRUE), error = function(e) NULL) for (tm in c("I_z", "C_z")) { if (is.null(m_ol) || !tm %in% rownames(summary(m_ol)$coefficients)) { add_rob(paste("ol", issue_name, dv_name, tm), issue_name, dv_name, "ordered_logit", ifelse(tm == "I_z", "Identity", "Comparison"), NA_real_, NA_real_, NA_real_, "polr failed / term missing") } else { co <- summary(m_ol)$coefficients b <- co[tm, "Value"]; s <- co[tm, "Std. Error"] p <- 2 * pnorm(-abs(b / s)) add_rob(paste("ol", issue_name, dv_name, tm), issue_name, dv_name, "ordered_logit", ifelse(tm == "I_z", "Identity", "Comparison"), b, s, p, "log-odds; normal-approx p") } } # (ii) multilevel: random intercept for language region m_ml <- tryCatch( lme4::lmer(as.formula(paste(dv_term, "~ I_z + C_z + FT_z +", controls_rhs_noreg, "+ (1 | sprachregion)")), data = d, REML = TRUE), error = function(e) NULL) for (tm in c("I_z", "C_z")) { if (is.null(m_ml)) { add_rob(paste("ml", issue_name, dv_name, tm), issue_name, dv_name, "multilevel_region", ifelse(tm == "I_z", "Identity", "Comparison"), NA_real_, NA_real_, NA_real_, "lmer failed") } else { co <- summary(m_ml)$coefficients b <- co[tm, "Estimate"]; s <- co[tm, "Std. Error"] p <- 2 * pnorm(-abs(b / s)) add_rob(paste("ml", issue_name, dv_name, tm), issue_name, dv_name, "multilevel_region", ifelse(tm == "I_z", "Identity", "Comparison"), b, s, p, "normal-approx p") } } # (iv) separate estimation by language region for (reg in levels(d$sprachregion)) { d_r <- d |> filter(sprachregion == reg) m_r <- tryCatch( lm(as.formula(paste(dv_term, "~ I_z + C_z + FT_z +", controls_rhs_noreg)), data = d_r), error = function(e) NULL) for (tm in c("I_z", "C_z")) { if (is.null(m_r)) { add_rob(paste("reg", reg, issue_name, dv_name, tm), issue_name, dv_name, paste0("region_", reg), ifelse(tm == "I_z", "Identity", "Comparison"), NA_real_, NA_real_, NA_real_, "lm failed") } else { e <- extract_hc3(m_r, tm) add_rob(paste("reg", reg, issue_name, dv_name, tm), issue_name, dv_name, paste0("region_", reg), ifelse(tm == "I_z", "Identity", "Comparison"), e$b, e$se, e$p, paste0("N=", nobs(m_r))) } } } } } robJ_df <- bind_rows(rob_rows) write.csv(robJ_df, out("TableJ_robustness.csv"), row.names = FALSE) # Console summary: sign/significance agreement with the base Model 2 robJ_base <- robJ_df |> filter(spec == "base_M2") |> select(issue, outcome, term, b_base = b, p_base = p) robJ_cmp <- robJ_df |> filter(!spec %in% "base_M2", !is.na(b)) |> left_join(robJ_base, by = c("issue", "outcome", "term")) |> group_by(spec) |> summarise( n_coefs = dplyr::n(), sign_match = mean(sign(b) == sign(b_base)), sig_match = mean((p < .05) == (p_base < .05), na.rm = TRUE), .groups = "drop" ) cat("\n=== Appendix J robustness: agreement with base Model 2 ===\n") print(robJ_cmp |> mutate(across(where(is.numeric), ~round(.x, 3)))) cat("Note: ordered-logit coefficients are log-odds; only sign/significance\n", "are comparable to the OLS base, not magnitudes.\n") # --- Identity x Comparison interactions with BH correction (Appendix D) ------ int_rows <- list() for (issue_name in names(issues)) { cfg <- issues[[issue_name]] sfx <- tolower(issue_name) d <- build_d_j(cfg, sfx) outcomes_j <- c(outcomes_shared, list("Political tolerance" = paste0("poltolerance_", sfx))) for (dv_name in names(outcomes_j)) { dv_var <- outcomes_j[[dv_name]] if (!dv_var %in% names(dat)) next dv_term <- paste0("scale(", dv_var, ")") m_int <- lm(as.formula(paste(dv_term, "~ I_z * C_z + FT_z +", controls_rhs)), data = d) e <- extract_hc3(m_int, "I_z:C_z") int_rows[[paste(issue_name, dv_name)]] <- tibble( issue = issue_name, outcome = dv_name, b_interaction = e$b, se = e$se, p_raw = e$p) } } int_df <- bind_rows(int_rows) |> mutate(p_BH = p.adjust(p_raw, method = "BH")) cat("\n=== Identity x Comparison interactions (HC3), BH-corrected across", nrow(int_df), "tests ===\n") print(int_df |> mutate(across(where(is.numeric), ~round(.x, 4)))) cat("Significant after BH (p_BH < .05):", sum(int_df$p_BH < .05), "of", nrow(int_df), "| min p_BH =", round(min(int_df$p_BH), 3), "\n") write.csv(int_df, out("TableD_interactions_BH.csv"), row.names = FALSE) # --- Export tables as CSV ---------------------------------------------------- write.csv(sample_desc, out("TableG1_sample_descriptives.csv"), row.names = FALSE) # (Table 1 and Table B2 are exported in Section 1.1, same run, same extraction) write.csv(table2, out("Table2_reliability.csv"), row.names = FALSE) write.csv(table_b1, out("TableB1_loadings.csv"), row.names = FALSE) write.csv(fig1_combined_data, out("Figure1_coefficients.csv"), row.names = FALSE) write.csv(conv_coefs_df, out("TableH1_H2_convergent_validity.csv"), row.names = FALSE) write.csv(crit_coefs_df, out("TableH3_criterion_validity.csv"), row.names = FALSE) write.csv(inc_table, out("TablesJ1_J4_summary_incremental.csv"), row.names = FALSE) write.csv(bind_rows(full_model_rows), out("TablesJ1_J4_full_models.csv"), row.names = FALSE) write.csv(fig3_df, out("Figure3_coefficients.csv"), row.names = FALSE) # --- Known limitations of this replication package ---------------------------- # (1) Off-diagonal profile analysis (share of high-identity/low-comparison # respondents) is reported in the text but not reproduced here: the # cutoff rule is not documented in the analysis files. Port the exact # rule from the original scripts before including it in the package. # (2) Table D1 (democratic-attitude operationalization) is hand-written in # the manuscript, not generated from data. # Quick DV range check for the A4-vs-D1 scale conflict: cat("\n=== DV scale ranges (settles manuscript A4 vs. Table D1 conflict) ===\n") for (v in c("demo_support_mean", "democracyevaluationm", "democracyconceptionm", "democracy_satisfaction", "trust_mean")) { if (v %in% names(dat)) { x <- suppressWarnings(as.numeric(dat[[v]])) cat(sprintf(" %-24s: %s to %s\n", v, round(min(x, na.rm = TRUE), 2), round(max(x, na.rm = TRUE), 2))) } } # Reproducibility: record package versions writeLines(capture.output(sessionInfo()), out("sessionInfo.txt")) cat("\n=== VERIFICATION CHECKLIST (compare against the manuscript) ===\n") cat(" Table 1 <- Table1_CFA_fit.csv\n") cat(" Table 2 <- Table2_reliability.csv\n") cat(" Table 3 <- Table3_APSIT_descriptives.csv (N should read 761/772/817/764)\n") cat(" Table G1 <- TableG1_sample_descriptives.csv (Female% near 50, NOT 0)\n") cat(" Tables H1/H2 <- TableH1_H2_convergent_validity.csv\n") cat(" NOTE: APSIT -> feeling-thermometer coefficients are POSITIVE here.\n") cat(" The thermometer outcome is the ingroup-minus-outgroup difference\n") cat(" index, so higher values mean greater relative outgroup coldness.\n") cat(" Tables I1-I4 <- TablesI1_I4_separate_models.csv (beta and Adj. R2)\n") cat(" Tables J1-J4 <- TablesJ1_J4_full_models.csv\n") cat(" Table J5 <- TableJ5_slope_differences.csv\n") cat(" (exactly one significant comparison difference expected)\n") cat(" Figure 3 <- plotted from the camp-main specification\n") cat("\n=== All tables and figures saved to:", output_dir, "===\n") cat("Files created:\n") print(list.files(output_dir)) cat("Run finished:", format(Sys.time()), "\n") cat("=== Replication script complete. ===\n") # Close the console log while (sink.number() > 0) sink() close(log_con) ### ======================================================================== ### ### End of Script ### ### ======================================================================== ###