Motivated reasoning and intelligence: renalysis
26 March, 2026
Init
library(kirkegaard)## Loading required package: tidyverse## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.1.9000 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.2.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: magrittr
##
##
## Attaching package: 'magrittr'
##
##
## The following object is masked from 'package:purrr':
##
## set_names
##
##
## The following object is masked from 'package:tidyr':
##
## extract
##
##
## Loading required package: weights
##
## Loading required package: assertthat
##
##
## Attaching package: 'assertthat'
##
##
## The following object is masked from 'package:tibble':
##
## has_name
##
##
## Loading required package: psych
##
##
## Attaching package: 'psych'
##
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
##
##
## Loading required package: robustbase
##
##
## Attaching package: 'kirkegaard'
##
##
## The following object is masked from 'package:psych':
##
## rescale
##
##
## The following object is masked from 'package:assertthat':
##
## are_equal
##
##
## The following object is masked from 'package:purrr':
##
## is_logical
##
##
## The following object is masked from 'package:base':
##
## +load_packages(
"tidyverse",
"metafor",
"patchwork",
"ggrepel"
)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Loading required package: metadat
## Loading required package: numDeriv
##
## Loading the 'metafor' package (version 4.8-0). For an
## introduction to the package please type: help(metafor)theme_set(theme_bw())
options(
digits = 3
)
#multithreading
#library(future)
#plan(multisession(workers = 8))Functions
Data
Persson et al. (2021) replication
d_repl <- read_csv("data/orig/mn-data.csv", show_col_types = FALSE)
#filter to attention check passers, as authors did
d_repl <- d_repl %>% filter(attention == 1)
#condition: 1=rash decrease, 2=rash increase, 3=crime increase, 4=crime decrease
#correct: 1=correct, 0=incorrect
#numeracy: 0-9 correct responses
#political: z-scored composite of ideology + party (positive = conservative/republican)
d_repl %>% str()## spc_tbl_ [3,073 × 32] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Progress : num [1:3073] 100 100 100 100 100 100 100 100 100 100 ...
## $ Duration : num [1:3073] 842 722 776 311 389 ...
## $ Finished : logi [1:3073] TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ consent : num [1:3073] 1 1 1 1 1 1 1 1 1 1 ...
## $ condition : num [1:3073] 4 3 2 4 2 2 3 4 2 3 ...
## $ rashDecrease : num [1:3073] NA NA NA NA NA NA NA NA NA NA ...
## $ rashIncrease : num [1:3073] NA NA 1 NA 1 2 NA NA 1 NA ...
## $ crimeIncrease: num [1:3073] NA 2 NA NA NA NA 1 NA NA 2 ...
## $ crimeDecrease: num [1:3073] 1 NA NA 2 NA NA NA 2 NA NA ...
## $ correct : num [1:3073] 1 1 0 0 0 1 0 0 0 1 ...
## $ num1 : num [1:3073] 500 500 60 500 166 500 3 500 500 500 ...
## $ num2 : num [1:3073] 10 10 0.1 10 100 10 10 10 10 100 ...
## $ num3 : num [1:3073] 1e-01 1e-01 1e-01 1e-02 1e+02 1e-01 1e-01 1e-01 1e-01 1e-01 ...
## $ num4 : num [1:3073] 20 20 2 0.5 20 20 20 20 20 20 ...
## $ num5 : num [1:3073] 100 100 100 10 100 100 100 100 100 100 ...
## $ num6 : num [1:3073] 0.5 0.9 0.19 0.1 0.1 ...
## $ num7 : num [1:3073] 4 5.5 18 9 9 9 9 3 4 5 ...
## $ num8 : num [1:3073] 0.25 0.25 0.5 0.5 0.5 0.25 0.5 0.25 0.5 0.5 ...
## $ num9 : num [1:3073] 27 27 56 14 14 27 14 27 27 27 ...
## $ numeracy : num [1:3073] 8 7 2 2 2 8 4 8 7 5 ...
## $ ideology : chr [1:3073] "Moderate" "Very Liberal" "Very Liberal" "Moderate" ...
## $ party : chr [1:3073] "Independent" "Independent" "Republican" "Independent Lean Democrat" ...
## $ political : num [1:3073] 0.146 -0.693 -0.134 -0.134 -0.413 ...
## $ political_d : num [1:3073] 1 0 0 0 0 0 0 0 0 0 ...
## $ sex : num [1:3073] 0 1 0 1 1 0 1 0 1 0 ...
## $ age : num [1:3073] 49 69 36 36 48 29 34 39 44 33 ...
## $ income : chr [1:3073] "$50,000 - $59,999" "$20,000 - $29,999" "$10,000 - $19,999" "$80,000 - $99,999" ...
## $ education : chr [1:3073] "2-year college" "2-year college" "4-year college" "Some college" ...
## $ attention : num [1:3073] 1 1 1 1 1 1 1 1 1 1 ...
## $ duplicate_id : num [1:3073] 0 0 1 0 0 0 0 0 0 0 ...
## $ IP_block : num [1:3073] 0 NA NA 0 NA 0 2 0 0 0 ...
## $ IP_country : chr [1:3073] "United States" NA NA "United States" ...
## - attr(*, "spec")=
## .. cols(
## .. Progress = col_double(),
## .. Duration = col_double(),
## .. Finished = col_logical(),
## .. consent = col_double(),
## .. condition = col_double(),
## .. rashDecrease = col_double(),
## .. rashIncrease = col_double(),
## .. crimeIncrease = col_double(),
## .. crimeDecrease = col_double(),
## .. correct = col_double(),
## .. num1 = col_double(),
## .. num2 = col_double(),
## .. num3 = col_double(),
## .. num4 = col_double(),
## .. num5 = col_double(),
## .. num6 = col_double(),
## .. num7 = col_double(),
## .. num8 = col_double(),
## .. num9 = col_double(),
## .. numeracy = col_double(),
## .. ideology = col_character(),
## .. party = col_character(),
## .. political = col_double(),
## .. political_d = col_double(),
## .. sex = col_double(),
## .. age = col_double(),
## .. income = col_character(),
## .. education = col_character(),
## .. attention = col_double(),
## .. duplicate_id = col_double(),
## .. IP_block = col_double(),
## .. IP_country = col_character()
## .. )
## - attr(*, "problems")=<externalptr>d_repl %>% count(condition)Analysis
Persson et al. replication: interaction models
Table 3 from their paper: logistic regression correct ~ numeracy * political per condition.
cond_labels <- c(
"3" = "Crime increase (cons-congruent)",
"4" = "Crime decrease (lib-congruent)",
"2" = "Rash got worse (control)",
"1" = "Rash got better (control)"
)
models_repl <- map(names(cond_labels), function(cond) {
glm(correct ~ numeracy * political,
data = d_repl %>% filter(condition == as.integer(cond)),
family = binomial)
}) %>% set_names(names(cond_labels))
#extract interaction coefficients (log-odds scale) for meta-analysis
repl_int <- map_dfr(names(cond_labels), function(cond) {
m <- models_repl[[cond]]
s <- summary(m)$coefficients
tibble(
study = "Persson et al. (2021)",
condition = cond_labels[cond],
cond_num = as.integer(cond),
b = s["numeracy:political", "Estimate"],
se = s["numeracy:political", "Std. Error"],
z = s["numeracy:political", "z value"],
p = s["numeracy:political", "Pr(>|z|)"],
n = nrow(d_repl %>% filter(condition == as.integer(cond)))
)
})
repl_int#verify: ORs match Table 3
repl_int %>% mutate(OR = exp(b))Full model summaries:
models_repl %>% map(summary)## $`3`
##
## Call:
## glm(formula = correct ~ numeracy * political, family = binomial,
## data = d_repl %>% filter(condition == as.integer(cond)))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0050 0.2087 -4.81 1.5e-06 ***
## numeracy 0.1443 0.0338 4.26 2.0e-05 ***
## political 0.4330 0.1929 2.25 0.025 *
## numeracy:political -0.0222 0.0321 -0.69 0.490
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1041.9 on 755 degrees of freedom
## Residual deviance: 1013.0 on 752 degrees of freedom
## AIC: 1021
##
## Number of Fisher Scoring iterations: 4
##
##
## $`4`
##
## Call:
## glm(formula = correct ~ numeracy * political, family = binomial,
## data = d_repl %>% filter(condition == as.integer(cond)))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8863 0.2361 -7.99 1.3e-15 ***
## numeracy 0.2056 0.0379 5.43 5.6e-08 ***
## political 0.9584 0.2187 4.38 1.2e-05 ***
## numeracy:political -0.1724 0.0356 -4.84 1.3e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 992.23 on 767 degrees of freedom
## Residual deviance: 938.29 on 764 degrees of freedom
## AIC: 946.3
##
## Number of Fisher Scoring iterations: 4
##
##
## $`2`
##
## Call:
## glm(formula = correct ~ numeracy * political, family = binomial,
## data = d_repl %>% filter(condition == as.integer(cond)))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1809 0.2263 -9.64 <2e-16 ***
## numeracy 0.3447 0.0364 9.46 <2e-16 ***
## political 0.0749 0.2238 0.33 0.74
## numeracy:political -0.0229 0.0361 -0.63 0.53
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1054.79 on 768 degrees of freedom
## Residual deviance: 940.13 on 765 degrees of freedom
## AIC: 948.1
##
## Number of Fisher Scoring iterations: 4
##
##
## $`1`
##
## Call:
## glm(formula = correct ~ numeracy * political, family = binomial,
## data = d_repl %>% filter(condition == as.integer(cond)))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.9437 0.2044 -4.62 3.9e-06 ***
## numeracy 0.1384 0.0338 4.10 4.2e-05 ***
## political 0.5851 0.1969 2.97 0.0030 **
## numeracy:political -0.0886 0.0330 -2.69 0.0072 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1077.6 on 779 degrees of freedom
## Residual deviance: 1054.1 on 776 degrees of freedom
## AIC: 1062
##
## Number of Fisher Scoring iterations: 4Stagnaro et al. (2023) PNAS replication
N=2,355 from AmeriSpeak probability panel. Five political topics (gun, health, immigration, police, taxes). Same paradigm: 2x2 contingency table, binary correct/incorrect DV.
d_stag <- read_csv("data/Stagnaro/Unmotivated_numeracy_selfGov_data_OSF.csv", show_col_types = FALSE)
#filter: exclude quality-flagged participants
d_stag <- d_stag %>%
filter(is.na(qual)) %>%
filter(!is.na(Correct), !is.na(RepCon), !is.na(Numericy))
#Numericy is 0-1 proportion; RepCon is raw (2-12)
d_stag %>% count(VinVersoin)#fit per-condition models for all topics
stag_models <- d_stag %>%
filter(!is.na(VinVersoin), !is.na(OutcomeIsConserv)) %>%
group_by(VinVersoin, OutcomeIsConserv) %>%
group_map(function(dd, keys) {
m <- glm(Correct ~ Numericy * RepCon, data = dd, family = binomial)
s <- summary(m)$coefficients
tibble(
study = paste0("Stagnaro (2023) - ", keys$VinVersoin),
condition = ifelse(keys$OutcomeIsConserv == 1, "Cons-congruent outcome", "Lib-congruent outcome"),
b = s["Numericy:RepCon", "Estimate"],
se = s["Numericy:RepCon", "Std. Error"],
n = nrow(dd)
)
}) %>% bind_rows()
stag_modelsStagnaro: replication of their key results (standardized betas)
They report standardized coefficients. We replicate exactly: concordance β=0.093, numeracy β=0.111, interaction β=0.002 (p=.94).
d_stag <- d_stag %>%
mutate(id_conc = ifelse(ID_Concord == "ID_Con", 1, 0))
#main model (standardized betas)
m_stag_main <- lm(scale(Correct) ~ scale(id_conc) * scale(Numericy),
data = d_stag %>% filter(!is.na(id_conc)))
summary(m_stag_main)$coefficients %>% round(4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0000 0.0214 0.0001 1.000
## scale(id_conc) 0.0969 0.0214 4.5276 0.000
## scale(Numericy) 0.1017 0.0214 4.7488 0.000
## scale(id_conc):scale(Numericy) -0.0002 0.0214 -0.0113 0.991#three-way with political ideology
m_stag_3way <- lm(Correct ~ id_conc * scale(Numericy) * scale(RepCon),
data = d_stag %>% filter(!is.na(id_conc), !is.na(RepCon)))
cat("Three-way: concordance × numeracy × ideology\n")## Three-way: concordance × numeracy × ideologysummary(m_stag_3way)$coefficients %>% round(4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4350 0.0151 28.799 0.0000
## id_conc 0.0959 0.0214 4.492 0.0000
## scale(Numericy) 0.0473 0.0150 3.144 0.0017
## scale(RepCon) -0.0106 0.0150 -0.705 0.4808
## id_conc:scale(Numericy) 0.0028 0.0214 0.129 0.8976
## id_conc:scale(RepCon) -0.0178 0.0214 -0.831 0.4060
## scale(Numericy):scale(RepCon) -0.0425 0.0150 -2.826 0.0048
## id_conc:scale(Numericy):scale(RepCon) 0.0205 0.0217 0.948 0.3433#split by ideology (median split on RepCon)
med_rep <- median(d_stag$RepCon, na.rm = TRUE)
cat("=== Liberals (RepCon below median) ===\n")## === Liberals (RepCon below median) ===lm(Correct ~ id_conc * scale(Numericy),
data = d_stag %>% filter(!is.na(id_conc), RepCon < med_rep)) %>%
summary() %>% .$coefficients %>% round(4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4388 0.0210 20.9338 0.000
## id_conc 0.1298 0.0300 4.3306 0.000
## scale(Numericy) 0.0886 0.0209 4.2480 0.000
## id_conc:scale(Numericy) 0.0008 0.0300 0.0261 0.979cat("\n=== Conservatives (RepCon above median) ===\n")##
## === Conservatives (RepCon above median) ===lm(Correct ~ id_conc * scale(Numericy),
data = d_stag %>% filter(!is.na(id_conc), RepCon > med_rep)) %>%
summary() %>% .$coefficients %>% round(4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4044 0.0249 16.253 0.0000
## id_conc 0.0909 0.0350 2.599 0.0095
## scale(Numericy) 0.0213 0.0247 0.863 0.3883
## id_conc:scale(Numericy) -0.0042 0.0350 -0.119 0.9049Stagnaro: replication of Figure 1A
d_fig <- d_stag %>%
filter(!is.na(id_conc), !is.na(Numericy)) %>%
mutate(concordance = ifelse(id_conc == 1, "Concordant", "Discordant"))
obs_means <- d_fig %>%
group_by(Numericy, concordance) %>%
summarise(prop_correct = mean(Correct), n = n(), .groups = "drop")
ggplot(obs_means, aes(x = Numericy, y = prop_correct, color = concordance, fill = concordance)) +
geom_point(size = 2.5) +
geom_smooth(data = d_fig, aes(y = Correct), method = "lm", se = TRUE, alpha = 0.15, linewidth = 1) +
geom_smooth(data = d_fig, aes(y = Correct), method = "loess", se = FALSE, linetype = "dotted", linewidth = 0.8) +
scale_color_manual(values = c("Concordant" = "forestgreen", "Discordant" = "darkorange")) +
scale_fill_manual(values = c("Concordant" = "forestgreen", "Discordant" = "darkorange")) +
scale_y_continuous(labels = scales::percent) +
coord_cartesian(ylim = c(0.2, 0.7)) +
labs(x = "Numeracy", y = "Proportion Correct",
title = "Stagnaro et al. (2023): replication of Figure 1A",
subtitle = "P(correct) by numeracy and identity concordance. Solid = linear, dotted = LOESS.",
color = NULL, fill = NULL) +
theme(legend.position = c(0.2, 0.85), legend.background = element_blank())## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Numeracy × ideology (ignoring concordance)
Does numeracy benefit liberals more than conservatives overall, regardless of condition?
#Stagnaro (standardized)
cat("=== Stagnaro: Correct ~ Numeracy * Ideology (standardized) ===\n")## === Stagnaro: Correct ~ Numeracy * Ideology (standardized) ===lm(scale(Correct) ~ scale(Numericy) * scale(RepCon),
data = d_stag %>% filter(!is.na(RepCon))) %>%
summary() %>% .$coefficients %>% round(4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0026 0.0214 0.119 0.9054
## scale(Numericy) 0.0987 0.0215 4.587 0.0000
## scale(RepCon) -0.0362 0.0215 -1.686 0.0919
## scale(Numericy):scale(RepCon) -0.0670 0.0217 -3.086 0.0021#Persson (standardized)
cat("\n=== Persson: Correct ~ numeracy * political (standardized) ===\n")##
## === Persson: Correct ~ numeracy * political (standardized) ===lm(correct ~ scale(numeracy) * scale(political), data = d_repl) %>%
summary() %>% .$coefficients %>% round(4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4189 0.0089 47.13 0.00
## scale(numeracy) 0.1120 0.0090 12.51 0.00
## scale(political) 0.0207 0.0089 2.33 0.02
## scale(numeracy):scale(political) -0.0380 0.0087 -4.36 0.00#Stagnaro
d_stag_plot <- d_stag %>%
filter(!is.na(RepCon), !is.na(Numericy), !is.na(Correct)) %>%
mutate(ideology = ifelse(RepCon <= median(RepCon, na.rm = TRUE), "Liberal", "Conservative"))
obs_s <- d_stag_plot %>%
group_by(Numericy, ideology) %>%
summarise(prop_correct = mean(Correct), n = n(), .groups = "drop")
p_s <- ggplot(obs_s, aes(x = Numericy, y = prop_correct, color = ideology, fill = ideology)) +
geom_point(aes(size = n), alpha = 0.7) +
geom_smooth(data = d_stag_plot, aes(y = Correct),
method = "glm", method.args = list(family = "binomial"),
se = TRUE, alpha = 0.15, linewidth = 1) +
scale_color_manual(values = c("Conservative" = "#E41A1C", "Liberal" = "#377EB8")) +
scale_fill_manual(values = c("Conservative" = "#E41A1C", "Liberal" = "#377EB8")) +
scale_size_continuous(range = c(1, 4), guide = "none") +
scale_y_continuous(labels = scales::percent) +
labs(x = "Numeracy (0-1)", y = "Proportion Correct", color = NULL, fill = NULL,
title = "Stagnaro (2023)", subtitle = "b = -0.065, p = .003") +
theme(legend.position = c(0.2, 0.85), legend.background = element_blank())
#Persson
d_repl_plot <- d_repl %>%
mutate(ideology = ifelse(political > 0, "Conservative", "Liberal"))
obs_p <- d_repl_plot %>%
group_by(numeracy, ideology) %>%
summarise(prop_correct = mean(correct), n = n(), .groups = "drop")
p_p <- ggplot(obs_p, aes(x = numeracy, y = prop_correct, color = ideology, fill = ideology)) +
geom_point(aes(size = n), alpha = 0.7) +
geom_smooth(data = d_repl_plot, aes(y = correct),
method = "glm", method.args = list(family = "binomial"),
se = TRUE, alpha = 0.15, linewidth = 1) +
scale_color_manual(values = c("Conservative" = "#E41A1C", "Liberal" = "#377EB8")) +
scale_fill_manual(values = c("Conservative" = "#E41A1C", "Liberal" = "#377EB8")) +
scale_size_continuous(range = c(1, 4), guide = "none") +
scale_y_continuous(labels = scales::percent) +
labs(x = "Numeracy (0-9)", y = "Proportion Correct", color = NULL, fill = NULL,
title = "Persson (2021)", subtitle = "b = -0.038, p < .001") +
theme(legend.position = c(0.2, 0.85), legend.background = element_blank())
p_p + p_s + plot_annotation(
title = "Numeracy x ideology interaction (all conditions pooled)",
subtitle = "Liberals benefit more from numeracy than conservatives. Dot size = n. Logistic fit.",
theme = theme(plot.title = element_text(size = 13, face = "bold"),
plot.subtitle = element_text(size = 10))
)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Concordance × numeracy × ideology: per topic and combined
The three-way interaction tests whether the motivated numeracy effect (concordance × numeracy) differs by political orientation.
d_s <- d_stag %>%
filter(!is.na(id_conc), !is.na(Correct), !is.na(Numericy), !is.na(RepCon), !is.na(VinVersoin)) %>%
mutate(z_num = as.numeric(scale(Numericy)),
z_pol = as.numeric(scale(RepCon)))
#per-topic
cat("=== Stagnaro: per-topic three-way interactions ===\n")## === Stagnaro: per-topic three-way interactions ===map_dfr(unique(d_s$VinVersoin), function(t) {
m <- lm(Correct ~ id_conc * z_num * z_pol, data = d_s %>% filter(VinVersoin == t))
s <- summary(m)$coefficients
tibble(
topic = t,
conc_num_b = s["id_conc:z_num", "Estimate"],
conc_num_p = s["id_conc:z_num", "Pr(>|t|)"],
threeway_b = s["id_conc:z_num:z_pol", "Estimate"],
threeway_p = s["id_conc:z_num:z_pol", "Pr(>|t|)"],
n = sum(d_s$VinVersoin == t)
)
})#combined with topic as covariate
cat("\n=== Stagnaro: combined model ===\n")##
## === Stagnaro: combined model ===m_sc <- lm(Correct ~ id_conc * z_num * z_pol + VinVersoin, data = d_s)
round(summary(m_sc)$coefficients[c("id_conc", "z_num", "z_pol", "id_conc:z_num",
"id_conc:z_pol", "z_num:z_pol", "id_conc:z_num:z_pol"), ], 4)## Estimate Std. Error t value Pr(>|t|)
## id_conc 0.0956 0.0214 4.472 0.0000
## z_num 0.0473 0.0150 3.147 0.0017
## z_pol -0.0113 0.0150 -0.749 0.4539
## id_conc:z_num 0.0022 0.0215 0.103 0.9183
## id_conc:z_pol -0.0170 0.0214 -0.792 0.4283
## z_num:z_pol -0.0423 0.0151 -2.808 0.0050
## id_conc:z_num:z_pol 0.0211 0.0217 0.975 0.3297#Persson: crime conditions only (rash has no concordance meaning)
d_p <- d_repl %>%
filter(condition %in% c(3, 4)) %>%
mutate(
id_conc = case_when(
condition == 3 & political > 0 ~ 1L,
condition == 3 & political <= 0 ~ 0L,
condition == 4 & political > 0 ~ 0L,
condition == 4 & political <= 0 ~ 1L
),
z_num = as.numeric(scale(numeracy)),
z_pol = as.numeric(scale(political))
)
cat("=== Persson: crime conditions combined ===\n")## === Persson: crime conditions combined ===m_pc <- lm(correct ~ id_conc * z_num * z_pol, data = d_p)
round(summary(m_pc)$coefficients[c("id_conc", "z_num", "z_pol", "id_conc:z_num",
"id_conc:z_pol", "z_num:z_pol", "id_conc:z_num:z_pol"), ], 4)## Estimate Std. Error t value Pr(>|t|)
## id_conc 0.0779 0.0250 3.1162 0.0019
## z_num 0.0612 0.0182 3.3670 0.0008
## z_pol -0.0151 0.0183 -0.8235 0.4104
## id_conc:z_num 0.0587 0.0252 2.3257 0.0202
## id_conc:z_pol 0.0953 0.0251 3.7925 0.0002
## z_num:z_pol -0.0456 0.0177 -2.5711 0.0102
## id_conc:z_num:z_pol -0.0019 0.0242 -0.0773 0.9384#Stagnaro: concordance × numeracy by ideology and topic
d_plot_s <- d_s %>%
mutate(ideology = ifelse(RepCon <= median(RepCon, na.rm = TRUE), "Liberal", "Conservative"),
concordance = ifelse(id_conc == 1, "Concordant", "Discordant"))
obs_s <- d_plot_s %>%
group_by(VinVersoin, Numericy, ideology, concordance) %>%
summarise(prop_correct = mean(Correct), n = n(), .groups = "drop")
ggplot(obs_s, aes(x = Numericy, y = prop_correct, color = concordance, fill = concordance)) +
geom_point(aes(size = n), alpha = 0.6) +
geom_smooth(data = d_plot_s, aes(y = Correct),
method = "glm", method.args = list(family = "binomial"),
se = TRUE, alpha = 0.1, linewidth = 0.8) +
facet_grid(ideology ~ VinVersoin) +
scale_color_manual(values = c("Concordant" = "forestgreen", "Discordant" = "darkorange")) +
scale_fill_manual(values = c("Concordant" = "forestgreen", "Discordant" = "darkorange")) +
scale_size_continuous(range = c(0.5, 3), guide = "none") +
scale_y_continuous(labels = scales::percent) +
labs(x = "Numeracy", y = "Proportion Correct", color = NULL, fill = NULL,
title = "Stagnaro (2023): concordance × numeracy by ideology and topic",
subtitle = "Motivated numeracy predicts diverging lines. Parallel lines = no interaction.") +
theme(legend.position = "bottom", strip.text = element_text(size = 9))## `geom_smooth()` using formula = 'y ~ x'
#Persson: concordance × numeracy by ideology and condition
d_plot_p <- d_repl %>%
filter(condition %in% c(3, 4)) %>%
mutate(
ideology = ifelse(political > 0, "Conservative", "Liberal"),
concordance = case_when(
condition == 3 & political > 0 ~ "Concordant",
condition == 3 & political <= 0 ~ "Discordant",
condition == 4 & political > 0 ~ "Discordant",
condition == 4 & political <= 0 ~ "Concordant"
),
cond_label = ifelse(condition == 3, "Crime increase", "Crime decrease")
)
obs_p <- d_plot_p %>%
group_by(cond_label, numeracy, ideology, concordance) %>%
summarise(prop_correct = mean(correct), n = n(), .groups = "drop")
ggplot(obs_p, aes(x = numeracy, y = prop_correct, color = concordance, fill = concordance)) +
geom_point(aes(size = n), alpha = 0.6) +
geom_smooth(data = d_plot_p, aes(y = correct),
method = "glm", method.args = list(family = "binomial"),
se = TRUE, alpha = 0.1, linewidth = 0.8) +
facet_grid(ideology ~ cond_label) +
scale_color_manual(values = c("Concordant" = "forestgreen", "Discordant" = "darkorange")) +
scale_fill_manual(values = c("Concordant" = "forestgreen", "Discordant" = "darkorange")) +
scale_size_continuous(range = c(0.5, 3), guide = "none") +
scale_y_continuous(labels = scales::percent) +
labs(x = "Numeracy (0-9)", y = "Proportion Correct", color = NULL, fill = NULL,
title = "Persson (2021): concordance × numeracy by ideology and condition",
subtitle = "Concordant = correct answer matches participant's politics.") +
theme(legend.position = "bottom")## `geom_smooth()` using formula = 'y ~ x'
Numeracy × ideology: does it survive a quadratic numeracy term?
The numeracy-accuracy relationship is visibly nonlinear (J-shaped) in Persson. A linear model may falsely attribute this curvature to an ideology interaction.
#Persson: linear vs quadratic
m_lin_p <- lm(correct ~ numeracy * political, data = d_repl)
m_quad_p <- lm(correct ~ numeracy + I(numeracy^2) + political +
numeracy:political + I(numeracy^2):political, data = d_repl)
cat("=== Persson: linear ===\n")## === Persson: linear ===round(summary(m_lin_p)$coefficients, 4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1578 0.0229 6.89 0
## numeracy 0.0476 0.0038 12.49 0
## political 0.1096 0.0221 4.97 0
## numeracy:political -0.0162 0.0037 -4.36 0cat("\n=== Persson: with quadratic ===\n")##
## === Persson: with quadratic ===round(summary(m_quad_p)$coefficients, 4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4630 0.0404 11.457 0.000
## numeracy -0.1013 0.0164 -6.165 0.000
## I(numeracy^2) 0.0145 0.0015 9.380 0.000
## political 0.0364 0.0376 0.967 0.334
## numeracy:political 0.0023 0.0157 0.146 0.884
## I(numeracy^2):political -0.0009 0.0015 -0.588 0.556cat("\nModel comparison:\n")##
## Model comparison:anova(m_lin_p, m_quad_p)#Stagnaro: linear vs quadratic
d_s <- d_stag %>% filter(!is.na(RepCon), !is.na(Correct))
m_lin_s <- lm(Correct ~ Numericy * RepCon, data = d_s)
m_quad_s <- lm(Correct ~ Numericy + I(Numericy^2) + RepCon +
Numericy:RepCon + I(Numericy^2):RepCon, data = d_s)
cat("=== Stagnaro: linear ===\n")## === Stagnaro: linear ===round(summary(m_lin_s)$coefficients, 4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2808 0.0582 4.83 0.0000
## Numericy 0.4827 0.1009 4.78 0.0000
## RepCon 0.0160 0.0083 1.93 0.0541
## Numericy:RepCon -0.0442 0.0143 -3.09 0.0021cat("\n=== Stagnaro: with quadratic ===\n")##
## === Stagnaro: with quadratic ===round(summary(m_quad_s)$coefficients, 4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2672 0.0922 2.8978 0.0038
## Numericy 0.4879 0.3866 1.2619 0.2071
## I(Numericy^2) 0.0127 0.3581 0.0355 0.9717
## RepCon 0.0344 0.0137 2.5109 0.0121
## Numericy:RepCon -0.1276 0.0563 -2.2658 0.0236
## I(Numericy^2):RepCon 0.0772 0.0517 1.4916 0.1359cat("\nModel comparison:\n")##
## Model comparison:anova(m_lin_s, m_quad_s)#visualize with quadratic fit
d_repl_plot <- d_repl %>%
filter(!is.na(political)) %>%
mutate(ideology = ifelse(political > 0, "Conservative", "Liberal"))
obs_p <- d_repl_plot %>%
group_by(numeracy, ideology) %>%
summarise(prop_correct = mean(correct), n = n(), .groups = "drop")
d_stag_plot <- d_stag %>%
filter(!is.na(RepCon), !is.na(Numericy), !is.na(Correct)) %>%
mutate(ideology = ifelse(RepCon <= median(RepCon, na.rm = TRUE), "Liberal", "Conservative"))
obs_s <- d_stag_plot %>%
group_by(Numericy, ideology) %>%
summarise(prop_correct = mean(Correct), n = n(), .groups = "drop")
p_p <- ggplot(obs_p, aes(x = numeracy, y = prop_correct, color = ideology, fill = ideology)) +
geom_point(aes(size = n), alpha = 0.7) +
geom_smooth(data = d_repl_plot, aes(y = correct),
method = "lm", formula = y ~ x + I(x^2),
se = TRUE, alpha = 0.15, linewidth = 1) +
scale_color_manual(values = c("Conservative" = "#E41A1C", "Liberal" = "#377EB8")) +
scale_fill_manual(values = c("Conservative" = "#E41A1C", "Liberal" = "#377EB8")) +
scale_size_continuous(range = c(1, 4), guide = "none") +
scale_y_continuous(labels = scales::percent) +
labs(x = "Numeracy (0-9)", y = "Proportion Correct", color = NULL, fill = NULL,
title = "Persson (2021)", subtitle = "Quadratic fit: interaction p = .88") +
theme(legend.position = c(0.2, 0.85), legend.background = element_blank())
p_s <- ggplot(obs_s, aes(x = Numericy, y = prop_correct, color = ideology, fill = ideology)) +
geom_point(aes(size = n), alpha = 0.7) +
geom_smooth(data = d_stag_plot, aes(y = Correct),
method = "lm", formula = y ~ x + I(x^2),
se = TRUE, alpha = 0.15, linewidth = 1) +
scale_color_manual(values = c("Conservative" = "#E41A1C", "Liberal" = "#377EB8")) +
scale_fill_manual(values = c("Conservative" = "#E41A1C", "Liberal" = "#377EB8")) +
scale_size_continuous(range = c(1, 4), guide = "none") +
scale_y_continuous(labels = scales::percent) +
labs(x = "Numeracy (0-1)", y = "Proportion Correct", color = NULL, fill = NULL,
title = "Stagnaro (2023)", subtitle = "Quadratic fit: interaction p = .063") +
theme(legend.position = c(0.2, 0.85), legend.background = element_blank())
p_p + p_s + plot_annotation(
title = "Numeracy x ideology with quadratic numeracy term",
subtitle = "The linear interaction vanishes once curvature is accounted for.",
theme = theme(plot.title = element_text(size = 13, face = "bold"),
plot.subtitle = element_text(size = 10))
)
Kahan et al. (2017) original
No raw data available. Coefficients extracted from Table A1 (Stage 3, combined logistic regression model). Reference condition = crime decrease. N=1111.
numeracy × Conserv_Repub: b = -0.33, z = -1.89numeracy × Conserv_Repub × crime_increases: b = 0.54, z = 2.17numeracy × Conserv_Repub × rash_increases: b = 0.26, z = 1.08numeracy × Conserv_Repub × rash_decreases: b = 0.33, z = 1.40
Per-condition interaction (numeracy on raw 0-9 scale, Conserv_Repub z-scored):
se_base <- 0.33 / 1.89
kahan_orig <- tibble(
study = "Kahan et al. (2017)",
condition = c("Crime increase (cons-congruent)",
"Crime decrease (lib-congruent)",
"Rash got worse (control)",
"Rash got better (control)"),
b = c(-0.33 + 0.54, -0.33, -0.33 + 0.26, -0.33 + 0.33),
#reference condition SE is exact; others are upper bounds (ignore covariance)
se = c(
sqrt(se_base^2 + (0.54/2.17)^2),
se_base,
sqrt(se_base^2 + (0.26/1.08)^2),
sqrt(se_base^2 + (0.33/1.40)^2)
),
n = rep(NA_integer_, 4)
)
kahan_origKahan & Peters (2017 SSRN) replication
Their own replication (N=780). Coefficients from Table 1 (combined logistic model, gun ban conditions). left_right is z-scored political orientation. Numeracy centered at 0 (raw 0-9 scale). Reference = crime decrease.
#gun ban model coefficients from SSRN Table 1
kp_ssrn_coefs <- tibble(
term = c("numeracy", "LR", "LR_x_num", "crime_increase",
"LR_x_crime_increase", "crime_increase_x_num",
"LR_x_crime_increase_x_num", "constant"),
b = c(0.20, -0.19, -0.16, 0.37, 0.43, 0.17, 0.23, 0.31),
z = c(1.88, -3.05, -2.43, 2.26, 4.76, 1.02, 2.43, 2.78)
) %>% mutate(se = abs(b / z))
kp_ssrn_coefsVerify: predicted probabilities vs their Figure 2
#reconstruct predictions from the model
b_kp <- setNames(kp_ssrn_coefs$b, kp_ssrn_coefs$term)
num_mean_kp <- 3.7 #approximate, same lab as original
predict_kp <- function(num_raw, LR_val, crime_inc) {
nc <- num_raw - num_mean_kp
logit <- b_kp["constant"] +
b_kp["numeracy"] * nc +
b_kp["LR"] * LR_val +
b_kp["LR_x_num"] * LR_val * nc +
b_kp["crime_increase"] * crime_inc +
b_kp["LR_x_crime_increase"] * LR_val * crime_inc +
b_kp["crime_increase_x_num"] * crime_inc * nc +
b_kp["LR_x_crime_increase_x_num"] * LR_val * crime_inc * nc
plogis(logit)
}
pred_kp <- expand.grid(numeracy = 0:9, LR = c(-0.95, 0.95), crime_inc = c(0, 1)) %>%
mutate(
prob = pmap_dbl(list(numeracy, LR, crime_inc), predict_kp),
group = ifelse(LR < 0, "Liberal Democrat", "Conservative Republican"),
condition = ifelse(crime_inc == 1,
"Crime Increase\n(identity affirmed for cons)",
"Crime Decrease\n(identity affirmed for libs)")
)
ggplot(pred_kp, aes(x = numeracy, y = prob, color = group, linetype = group)) +
geom_line(linewidth = 1) + geom_point(size = 1.5) +
facet_wrap(~condition) +
scale_color_manual(values = c("Conservative Republican" = "#E41A1C", "Liberal Democrat" = "#377EB8")) +
scale_y_continuous(limits = c(0, 1), labels = scales::percent) +
labs(x = "Numeracy score", y = "P(correct)", color = NULL, linetype = NULL,
title = "Predicted probabilities from Kahan & Peters (2017 SSRN) Table 1",
subtitle = "Gun ban model. Political orientation at ±0.95 SD.") +
theme(legend.position = "bottom")
Identity-affirmed vs identity-threatened (collapsed across party)
pred_affirm <- expand.grid(numeracy = 0:9, scenario = c("Identity affirmed", "Identity threatened")) %>%
rowwise() %>%
mutate(prob = if (scenario == "Identity affirmed") {
mean(c(predict_kp(numeracy, 0.95, 1), predict_kp(numeracy, -0.95, 0)))
} else {
mean(c(predict_kp(numeracy, 0.95, 0), predict_kp(numeracy, -0.95, 1)))
}) %>%
ungroup()
pred_wide_kp <- pred_affirm %>%
pivot_wider(names_from = scenario, values_from = prob) %>%
mutate(gap = `Identity affirmed` - `Identity threatened`)
p_aff <- ggplot(pred_affirm, aes(x = numeracy, y = prob, color = scenario, linetype = scenario)) +
geom_line(linewidth = 1.2) + geom_point(size = 2) +
scale_color_manual(values = c("Identity affirmed" = "forestgreen", "Identity threatened" = "darkorange")) +
scale_y_continuous(limits = c(0, 1), labels = scales::percent) +
labs(x = "Numeracy score", y = "P(correct)", color = NULL, linetype = NULL,
title = "Predicted from SSRN Table 1", subtitle = "Collapsed across party") +
theme(legend.position = "bottom")
p_gap_kp <- ggplot(pred_wide_kp, aes(x = numeracy, y = gap)) +
geom_line(linewidth = 1.2) + geom_point(size = 2) +
geom_hline(yintercept = 0, lty = 2, color = "grey50") +
scale_y_continuous(labels = scales::percent) +
labs(x = "Numeracy score", y = "Gap (affirmed − threatened)",
title = "Motivated numeracy gap", subtitle = "Affirmed minus threatened")
p_aff + p_gap_kp +
plot_annotation(title = "Kahan & Peters (2017 SSRN) replication: gun ban conditions",
theme = theme(plot.title = element_text(size = 14, face = "bold")))
Meta-analysis
We use two approaches to meta-analyze the interaction effects:
- Logistic interaction coefficients (z_numeracy × z_political): from per-condition logistic regressions with z-scored predictors. Kahan coefficients approximated from Table A1 combined model.
- Probability-scale difference-in-differences: change in partisan gap (percentage points) from low (3) to high (7) numeracy, evaluated at ±1 SD political orientation. Uses Kahan Figure 8 values (as read by Persson et al.).
For both approaches, sign is flipped for lib-congruent conditions so that positive = motivated numeracy (polarization amplified by numeracy).
Approach 1: Logistic interaction coefficients
#Persson: z-score both predictors
d_repl <- d_repl %>%
mutate(z_num = as.numeric(scale(numeracy)),
z_pol = as.numeric(scale(political)))
repl_z <- map_dfr(names(cond_labels), function(cond) {
m <- glm(correct ~ z_num * z_pol,
data = d_repl %>% filter(condition == as.integer(cond)),
family = binomial)
s <- summary(m)$coefficients
tibble(
study = "Persson et al. (2021)",
condition = cond_labels[cond],
b = s["z_num:z_pol", "Estimate"],
se = s["z_num:z_pol", "Std. Error"],
n = nrow(d_repl %>% filter(condition == as.integer(cond)))
)
})
#Stagnaro: z-score both predictors
d_stag <- d_stag %>%
mutate(z_num = as.numeric(scale(Numericy)),
z_pol = as.numeric(scale(RepCon)))
stag_z <- d_stag %>%
filter(!is.na(VinVersoin), !is.na(OutcomeIsConserv)) %>%
nest(data = -c(VinVersoin, OutcomeIsConserv)) %>%
mutate(
model = map(data, ~glm(Correct ~ z_num * z_pol, data = .x, family = binomial)),
tidy = map(model, ~{
s <- summary(.x)$coefficients
tibble(b = s["z_num:z_pol", "Estimate"], se = s["z_num:z_pol", "Std. Error"])
}),
n = map_int(data, nrow)
) %>%
unnest(tidy) %>%
mutate(
study = paste0("Stagnaro (2023) - ", VinVersoin),
condition = ifelse(OutcomeIsConserv == 1, "Cons-congruent", "Lib-congruent")
) %>%
select(study, condition, b, se, n)
#Kahan: convert from raw numeracy to z-scored scale (multiply by SD_numeracy=2.1)
se_base <- 0.33 / 1.89
kahan_z <- tibble(
study = "Kahan et al. (2017)",
condition = c("Crime inc (cons-congruent)", "Crime dec (lib-congruent)",
"Rash inc (control)", "Rash dec (control)"),
b = c(0.21, -0.33, -0.07, 0.00) * 2.1,
se = c(sqrt(se_base^2 + (0.54/2.17)^2), se_base,
sqrt(se_base^2 + (0.26/1.08)^2), sqrt(se_base^2 + (0.33/1.40)^2)) * 2.1,
n = rep(NA_integer_, 4)
)
#Kahan & Peters SSRN replication: from Table 1
#Crime decrease (reference): LR_x_num = -0.16, z = -2.43
#Crime increase: LR_x_num + LR_x_crime_increase_x_num = -0.16 + 0.23 = 0.07
se_lr_num <- 0.16 / 2.43
se_3way <- 0.23 / 2.43
kp_z <- tibble(
study = "Kahan & Peters (2017 SSRN)",
condition = c("Crime inc (cons-congruent)", "Crime dec (lib-congruent)"),
b = c(0.07, -0.16) * 2.1, #convert to z_numeracy scale
se = c(sqrt(se_lr_num^2 + se_3way^2), se_lr_num) * 2.1,
n = rep(NA_integer_, 2)
)
meta_logistic <- bind_rows(repl_z, stag_z, kahan_z, kp_z) %>%
mutate(
type = case_when(
grepl("cons-congruent|Cons-congruent", condition, ignore.case = TRUE) ~ "cons-congruent",
grepl("lib-congruent|Lib-congruent", condition, ignore.case = TRUE) ~ "lib-congruent",
TRUE ~ "control"
),
yi = b,
vi = se^2,
#flip lib-congruent so positive = motivated numeracy
yi_flip = ifelse(type == "lib-congruent", -b, b)
)
meta_logistic %>% select(study, condition, type, b, se, yi_flip)#political conditions only (positive = motivated numeracy)
pol_log <- meta_logistic %>% filter(type != "control")
cat("=== All political conditions ===\n")## === All political conditions ===res_pol <- rma(yi = yi_flip, vi = vi, data = pol_log, method = "REML")
res_pol##
## Random-Effects Model (k = 16; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.0348 (SE = 0.0212)
## tau (square root of estimated tau^2 value): 0.1864
## I^2 (total heterogeneity / total variability): 64.18%
## H^2 (total variability / sampling variability): 2.79
##
## Test for Heterogeneity:
## Q(df = 15) = 44.9537, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0825 0.0617 1.3382 0.1808 -0.0384 0.2034
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n=== Cons-congruent only ===\n")##
## === Cons-congruent only ===res_cc <- rma(yi = yi_flip, vi = vi, data = pol_log %>% filter(type == "cons-congruent"), method = "REML")
res_cc##
## Random-Effects Model (k = 8; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0089)
## tau (square root of estimated tau^2 value): 0
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 1.00
##
## Test for Heterogeneity:
## Q(df = 7) = 6.2072, p-val = 0.5158
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0939 0.0477 -1.9685 0.0490 -0.1874 -0.0004 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n=== Lib-congruent only ===\n")##
## === Lib-congruent only ===res_lc <- rma(yi = yi_flip, vi = vi, data = pol_log %>% filter(type == "lib-congruent"), method = "REML")
res_lc##
## Random-Effects Model (k = 8; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.0095 (SE = 0.0159)
## tau (square root of estimated tau^2 value): 0.0977
## I^2 (total heterogeneity / total variability): 31.94%
## H^2 (total variability / sampling variability): 1.47
##
## Test for Heterogeneity:
## Q(df = 7) = 10.2156, p-val = 0.1767
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2521 0.0624 4.0391 <.0001 0.1298 0.3744 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n=== Controls only ===\n")##
## === Controls only ===res_ctrl <- rma(yi = yi_flip, vi = vi, data = meta_logistic %>% filter(type == "control"), method = "REML")
res_ctrl##
## Random-Effects Model (k = 4; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.0044 (SE = 0.0155)
## tau (square root of estimated tau^2 value): 0.0666
## I^2 (total heterogeneity / total variability): 19.13%
## H^2 (total variability / sampling variability): 1.24
##
## Test for Heterogeneity:
## Q(df = 3) = 1.8571, p-val = 0.6026
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.1331 0.0730 -1.8230 0.0683 -0.2761 0.0100 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest(res_pol, slab = pol_log$study %>% paste("-", pol_log$condition),
xlab = "Interaction (positive = motivated numeracy)", cex = 0.65,
main = "Logistic interaction: political conditions")
Approach 2: Probability-scale difference-in-differences
Change in partisan gap (cons - lib, percentage points) from numeracy=3 to numeracy=7. Values from Kahan Figure 8 (read off by Persson et al.) and computed from fitted models.
#Kahan Figure 8 values
kahan_fig <- tibble(
condition = rep(c("Rash dec", "Rash inc", "Crime dec", "Crime inc"), 2),
numeracy_level = c(rep("low", 4), rep("high", 4)),
gap = c(-4, 0, -23, 27, -4, -5, -50, 40),
ci_lo = c(-18, -14, -37, 13, -26, -30, -70, 20),
ci_hi = c(10, 14, -9, 41, 18, 20, -30, 60)
) %>% mutate(se = (ci_hi - ci_lo) / (2 * 1.96))
#Persson Figure 4 values (from mn-code-3-figures.R)
persson_fig <- tibble(
condition = rep(c("Rash dec", "Rash inc", "Crime dec", "Crime inc"), 2),
numeracy_level = c(rep("low", 4), rep("high", 4)),
gap = c(14.8, 0.2, 15.1, 16.8, -1.8, -4.2, -11.8, 13.8),
ci_lo = c(5.0, -9.0, 7.2, 7.3, -10.5, -12.9, -20.5, 5.0),
ci_hi = c(24.6, 9.4, 23.0, 26.2, 7.0, 4.5, -3.0, 22.5)
) %>% mutate(se = (ci_hi - ci_lo) / (2 * 1.96))
#compute difference-in-differences
compute_did <- function(df, study_name) {
df %>%
select(condition, numeracy_level, gap, se) %>%
pivot_wider(names_from = numeracy_level, values_from = c(gap, se)) %>%
mutate(
study = study_name,
did = gap_high - gap_low,
did_se = sqrt(se_low^2 + se_high^2) #upper-bound SE (assumes independence)
) %>%
select(study, condition, gap_low, gap_high, did, did_se)
}
kahan_did <- compute_did(kahan_fig, "Kahan et al. (2017)")
persson_did <- compute_did(persson_fig, "Persson et al. (2021)")
#Stagnaro: compute from fitted models
rep_sd <- sd(d_stag$RepCon, na.rm = TRUE)
rep_mean <- mean(d_stag$RepCon, na.rm = TRUE)
stag_did_list <- list()
for (topic in c("gun", "health", "immigration", "police", "tax")) {
for (oc in c(0, 1)) {
dd <- d_stag %>% filter(VinVersoin == topic, OutcomeIsConserv == oc)
if (nrow(dd) < 20) next
m <- glm(Correct ~ Numericy * RepCon, data = dd, family = binomial)
newdata <- expand.grid(Numericy = c(3/9, 7/9), RepCon = rep_mean + c(-1, 1) * rep_sd)
preds <- predict(m, newdata, type = "response", se.fit = TRUE)
newdata$pred <- preds$fit
newdata$se_pred <- preds$se.fit
gap_low <- newdata$pred[newdata$Numericy == 3/9 & newdata$RepCon > rep_mean] -
newdata$pred[newdata$Numericy == 3/9 & newdata$RepCon < rep_mean]
gap_high <- newdata$pred[newdata$Numericy == 7/9 & newdata$RepCon > rep_mean] -
newdata$pred[newdata$Numericy == 7/9 & newdata$RepCon < rep_mean]
se_low <- sqrt(sum(newdata$se_pred[newdata$Numericy == 3/9]^2))
se_high <- sqrt(sum(newdata$se_pred[newdata$Numericy == 7/9]^2))
did <- gap_high - gap_low
did_se <- sqrt(se_low^2 + se_high^2)
stag_did_list[[length(stag_did_list) + 1]] <- tibble(
study = paste0("Stagnaro (2023) - ", topic),
condition = ifelse(oc == 1, "Cons-congruent", "Lib-congruent"),
gap_low = gap_low * 100, gap_high = gap_high * 100,
did = did * 100, did_se = did_se * 100
)
}
}
stag_did <- bind_rows(stag_did_list)
all_did <- bind_rows(kahan_did, persson_did, stag_did) %>%
mutate(
type = case_when(
grepl("Crime inc|Cons-congruent", condition) ~ "cons-congruent",
grepl("Crime dec|Lib-congruent", condition) ~ "lib-congruent",
TRUE ~ "control"
),
yi = did, vi = did_se^2,
#flip lib-congruent so positive = motivated numeracy
yi_flip = ifelse(type == "lib-congruent", -did, did),
label = paste(study, "-", condition)
)
all_did %>% select(label, type, did, did_se, yi_flip)pol_did <- all_did %>% filter(type != "control")
cat("=== All political conditions (pp scale) ===\n")## === All political conditions (pp scale) ===res_pol_pp <- rma(yi = yi_flip, vi = vi, data = pol_did, method = "REML")
res_pol_pp##
## Random-Effects Model (k = 14; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 169.1855 (SE = 118.3592)
## tau (square root of estimated tau^2 value): 13.0071
## I^2 (total heterogeneity / total variability): 59.08%
## H^2 (total variability / sampling variability): 2.44
##
## Test for Heterogeneity:
## Q(df = 13) = 34.1215, p-val = 0.0012
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 5.3175 4.6731 1.1379 0.2552 -3.8416 14.4766
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n=== Cons-congruent only ===\n")##
## === Cons-congruent only ===rma(yi = yi_flip, vi = vi, data = pol_did %>% filter(type == "cons-congruent"), method = "REML")##
## Random-Effects Model (k = 7; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.0000 (SE = 64.9494)
## tau (square root of estimated tau^2 value): 0.0017
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 1.00
##
## Test for Heterogeneity:
## Q(df = 6) = 6.9318, p-val = 0.3272
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -5.9854 4.0630 -1.4731 0.1407 -13.9487 1.9780
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n=== Lib-congruent only ===\n")##
## === Lib-congruent only ===rma(yi = yi_flip, vi = vi, data = pol_did %>% filter(type == "lib-congruent"), method = "REML")##
## Random-Effects Model (k = 7; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 18.9049 (SE = 76.9279)
## tau (square root of estimated tau^2 value): 4.3480
## I^2 (total heterogeneity / total variability): 13.39%
## H^2 (total variability / sampling variability): 1.15
##
## Test for Heterogeneity:
## Q(df = 6) = 5.9781, p-val = 0.4257
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 19.1772 4.4232 4.3356 <.0001 10.5079 27.8465 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n=== Controls only ===\n")##
## === Controls only ===rma(yi = yi_flip, vi = vi, data = all_did %>% filter(type == "control"), method = "REML")##
## Random-Effects Model (k = 4; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 7.8977 (SE = 65.3119)
## tau (square root of estimated tau^2 value): 2.8103
## I^2 (total heterogeneity / total variability): 8.81%
## H^2 (total variability / sampling variability): 1.10
##
## Test for Heterogeneity:
## Q(df = 3) = 2.3230, p-val = 0.5081
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -8.6829 4.5176 -1.9220 0.0546 -17.5373 0.1715 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest(res_pol_pp, slab = pol_did$label,
xlab = "Change in partisan gap, pp (positive = motivated numeracy)", cex = 0.65,
main = "Probability-scale: political conditions")
Plots
Forest plot: logistic interaction (all political conditions)
pol_all <- meta_logistic %>%
filter(type != "control") %>%
arrange(type, study) %>%
mutate(short_label = paste0(
str_replace(study, "Stagnaro \\(2023\\) - ", "Stagnaro: "), " [", condition, "]"
))
res_cc <- rma(yi = yi_flip, vi = vi, data = pol_all %>% filter(type == "cons-congruent"), method = "REML")
res_lc <- rma(yi = yi_flip, vi = vi, data = pol_all %>% filter(type == "lib-congruent"), method = "REML")
pooled_df <- tibble(
short_label = c("Pooled: Cons-congruent", "Pooled: Lib-congruent"),
type = c("cons-congruent", "lib-congruent"),
yi_flip = c(as.numeric(res_cc$beta), as.numeric(res_lc$beta)),
se = c(res_cc$se, res_lc$se),
ci_lo = c(res_cc$ci.lb, res_lc$ci.lb),
ci_hi = c(res_cc$ci.ub, res_lc$ci.ub),
is_pooled = TRUE
)
forest_df <- pol_all %>%
mutate(ci_lo = yi_flip - 1.96 * se, ci_hi = yi_flip + 1.96 * se, is_pooled = FALSE) %>%
select(short_label, type, yi_flip, se, ci_lo, ci_hi, is_pooled) %>%
bind_rows(pooled_df) %>%
mutate(
type_label = ifelse(type == "cons-congruent", "Cons-congruent conditions", "Lib-congruent conditions"),
short_label = factor(short_label, levels = rev(short_label))
)
ggplot(forest_df, aes(x = yi_flip, y = short_label, xmin = ci_lo, xmax = ci_hi)) +
geom_vline(xintercept = 0, lty = 2, color = "grey50") +
geom_pointrange(aes(shape = is_pooled, size = is_pooled)) +
scale_shape_manual(values = c("FALSE" = 16, "TRUE" = 18), guide = "none") +
scale_size_manual(values = c("FALSE" = 0.4, "TRUE" = 0.7), guide = "none") +
facet_wrap(~type_label, scales = "free_y", ncol = 1) +
labs(x = "Interaction coefficient (positive = motivated numeracy)", y = NULL,
title = "Meta-analysis: numeracy × political orientation interaction",
subtitle = "Logistic regression, z-scored predictors. 4 studies (Kahan originals approx. from tables).") +
theme(strip.text = element_text(face = "bold"))
Summary plot: pooled meta-analytic effects
res_all_pol <- rma(yi = yi_flip, vi = vi, data = pol_all, method = "REML")
res_ctrl <- rma(yi = yi_flip, vi = vi, data = meta_logistic %>% filter(type == "control"), method = "REML")
summary_df <- tibble(
subset = c(
paste0("All political\n(k=", nrow(pol_all), ")"),
paste0("Cons-congruent\n(k=", sum(pol_all$type == "cons-congruent"), ")"),
paste0("Lib-congruent\n(k=", sum(pol_all$type == "lib-congruent"), ")"),
paste0("Controls\n(k=", sum(meta_logistic$type == "control"), ")")
),
estimate = c(as.numeric(res_all_pol$beta), as.numeric(res_cc$beta),
as.numeric(res_lc$beta), as.numeric(res_ctrl$beta)),
ci_lo = c(res_all_pol$ci.lb, res_cc$ci.lb, res_lc$ci.lb, res_ctrl$ci.lb),
ci_hi = c(res_all_pol$ci.ub, res_cc$ci.ub, res_lc$ci.ub, res_ctrl$ci.ub),
p = c(res_all_pol$pval, res_cc$pval, res_lc$pval, res_ctrl$pval),
I2 = c(res_all_pol$I2, res_cc$I2, res_lc$I2, res_ctrl$I2)
) %>%
mutate(
subset = factor(subset, levels = rev(subset)),
sig = ifelse(p < .05, "p < .05", "n.s."),
p_label = ifelse(p < .001, "p < .001", paste0("p = ", formatC(p, format = "f", digits = 3))),
i2_label = paste0("I² = ", round(I2, 0), "%")
)
ggplot(summary_df, aes(x = estimate, y = subset, xmin = ci_lo, xmax = ci_hi, color = sig)) +
geom_vline(xintercept = 0, lty = 2, color = "grey50") +
geom_pointrange(size = 0.8, linewidth = 1) +
geom_text(aes(label = paste0(p_label, ", ", i2_label), x = ci_hi + 0.02),
hjust = 0, size = 3.2, color = "grey30") +
scale_color_manual(values = c("p < .05" = "black", "n.s." = "grey50"), guide = "none") +
scale_x_continuous(limits = c(-0.3, 0.65)) +
labs(x = "Pooled interaction (positive = motivated numeracy)", y = NULL,
title = "Meta-analytic evidence for motivated numeracy",
subtitle = "4 studies: Kahan (2017), Kahan & Peters (2017), Persson (2021), Stagnaro (2023)")
Numeracy slopes by ideology and condition
#compute numeracy slopes at ±1 SD political for each study × condition
num_sd_persson <- sd(d_repl$numeracy)
num_sd_stagnaro <- sd(d_stag$Numericy, na.rm = TRUE)
#Persson
slopes_persson <- map_dfr(names(cond_labels), function(cond) {
m <- models_repl[[cond]]
pol_sd <- sd(d_repl$political)
b_num <- coef(m)["numeracy"]
b_int <- coef(m)["numeracy:political"]
tibble(
study = "Persson (2021)",
condition = cond_labels[cond],
type = ifelse(grepl("Crime inc", cond_labels[cond]), "cons-congruent",
ifelse(grepl("Crime dec", cond_labels[cond]), "lib-congruent", "control")),
group = c("Conservative", "Liberal"),
slope = c(b_num + b_int * pol_sd, b_num + b_int * (-pol_sd))
)
})
#Stagnaro
rep_sd <- sd(d_stag$RepCon, na.rm = TRUE)
slopes_stagnaro <- d_stag %>%
filter(!is.na(VinVersoin), !is.na(OutcomeIsConserv)) %>%
nest(data = -c(VinVersoin, OutcomeIsConserv)) %>%
mutate(slopes = map(data, function(dd) {
m <- glm(Correct ~ Numericy * RepCon, data = dd, family = binomial)
b_num <- coef(m)["Numericy"]
b_int <- coef(m)["Numericy:RepCon"]
tibble(group = c("Conservative", "Liberal"),
slope = c(b_num + b_int * rep_sd, b_num + b_int * (-rep_sd)))
})) %>%
unnest(slopes) %>%
mutate(
study = paste0("Stagnaro: ", VinVersoin),
condition = ifelse(OutcomeIsConserv == 1, "Cons-congruent", "Lib-congruent"),
type = ifelse(OutcomeIsConserv == 1, "cons-congruent", "lib-congruent")
) %>%
select(study, condition, type, group, slope)
all_slopes <- bind_rows(slopes_persson, slopes_stagnaro) %>%
filter(type != "control") %>%
mutate(slope_std = case_when(
grepl("Persson", study) ~ slope * num_sd_persson,
TRUE ~ slope * num_sd_stagnaro
))
ggplot(all_slopes, aes(x = group, y = slope_std, color = group)) +
geom_line(aes(group = study), color = "grey60", linewidth = 0.5) +
geom_point(size = 3) +
geom_text_repel(data = all_slopes %>% filter(group == "Liberal"),
aes(label = str_replace(study, "Stagnaro: ", "") %>%
str_replace("Persson \\(2021\\)", "crime (Persson)")),
nudge_x = 0.3, size = 3, direction = "y", hjust = 0,
segment.color = "grey80", color = "grey30") +
scale_color_manual(values = c("Conservative" = "#E41A1C", "Liberal" = "#377EB8"), guide = "none") +
facet_wrap(~type, labeller = labeller(type = c(
"cons-congruent" = "Cons-congruent conditions",
"lib-congruent" = "Lib-congruent conditions"))) +
geom_hline(yintercept = 0, lty = 2, color = "grey50") +
scale_x_discrete(expand = expansion(mult = c(0.1, 0.5))) +
labs(x = NULL, y = "Numeracy slope (per 1 SD, log-odds)",
title = "Effect of numeracy on correct responses by political orientation",
subtitle = "Each line connects conservative and liberal slopes within one study × condition")
Session info
#versions
write_sessioninfo()## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Linux Mint 21.1
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0 LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_DK.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_DK.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_DK.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Brussels
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggrepel_0.9.6 patchwork_1.3.2.9000 metafor_4.8-0
## [4] numDeriv_2016.8-1.1 metadat_1.4-0 Matrix_1.7-4
## [7] kirkegaard_2025-11-19 robustbase_0.99-6 psych_2.5.6
## [10] assertthat_0.2.1 weights_1.1.2 magrittr_2.0.4
## [13] lubridate_1.9.4 forcats_1.0.1 stringr_1.6.0
## [16] dplyr_1.1.4 purrr_1.2.0 readr_2.1.5
## [19] tidyr_1.3.1 tibble_3.3.0 ggplot2_4.0.1.9000
## [22] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] Rdpack_2.6.4 mnormt_2.1.1 gridExtra_2.3 rlang_1.1.6.9000
## [5] compiler_4.5.2 mgcv_1.9-3 gdata_3.0.1 vctrs_0.6.5
## [9] pkgconfig_2.0.3 shape_1.4.6.1 crayon_1.5.3 fastmap_1.2.0
## [13] backports_1.5.0 labeling_0.4.3 rmarkdown_2.30 tzdb_0.5.0
## [17] nloptr_2.2.1 bit_4.6.0 xfun_0.57 glmnet_4.1-10
## [21] jomo_2.7-6 cachem_1.1.0 jsonlite_2.0.0 pan_1.9
## [25] broom_1.0.11 parallel_4.5.2 cluster_2.1.8.2 R6_2.6.1
## [29] bslib_0.9.0 stringi_1.8.7 RColorBrewer_1.1-3 boot_1.3-32
## [33] rpart_4.1.24 jquerylib_0.1.4 Rcpp_1.1.0 iterators_1.0.14
## [37] knitr_1.50 base64enc_0.1-3 splines_4.5.2 nnet_7.3-20
## [41] timechange_0.3.0 tidyselect_1.2.1 rstudioapi_0.17.1 yaml_2.3.10
## [45] codetools_0.2-20 lattice_0.22-7 withr_3.0.2 S7_0.2.1
## [49] evaluate_1.0.5 foreign_0.8-91 archive_1.1.12 survival_3.8-6
## [53] pillar_1.11.1 mice_3.18.0 checkmate_2.3.3 foreach_1.5.2
## [57] reformulas_0.4.1 generics_0.1.4 vroom_1.6.6 mathjaxr_2.0-0
## [61] hms_1.1.3 scales_1.4.0 minqa_1.2.8 gtools_3.9.5
## [65] glue_1.8.0 Hmisc_5.2-4 tools_4.5.2 data.table_1.17.8
## [69] lme4_1.1-37 grid_4.5.2 rbibutils_2.3 colorspace_2.1-2
## [73] nlme_3.1-168 htmlTable_2.4.3 Formula_1.2-5 cli_3.6.5
## [77] gtable_0.3.6 DEoptimR_1.1-4 sass_0.4.10 digest_0.6.39
## [81] htmlwidgets_1.6.4 farver_2.1.2 htmltools_0.5.8.1 lifecycle_1.0.4
## [85] mitml_0.4-5 bit64_4.6.0-1 MASS_7.3-65#write data to file for reuse
list(d_repl = d_repl, d_stag = d_stag, meta_logistic = meta_logistic, all_did = all_did) %>%
write_rds("data/data_for_reuse.rds")
#OSF
if (F) {
library(osfr)
#login
osf_auth(readr::read_lines("~/.config/osf_token"))
#the project we will use
osf_proj = osf_retrieve_node("https://osf.io/XXX/")
#upload all files in project
#overwrite existing (versioning)
osf_upload(
osf_proj,
path = c("data", "figures", "papers", "notebook.Rmd", "notebook.html", "sessions_info.txt"),
conflicts = "overwrite"
)
}