03_CRI_Spec_Analysis.Rmd

---
title: "Tech 3. Specification Curves"
output: html_document
---

```{r setup, include=FALSE}
rm(list = ls())
library("pacman")

pacman::p_load("devtools","ggplot2","tidyverse","readr","ExPanDaR","plotscale","lattice","tidyr","mlogit","knitr","grid","zoo","ggpubr","ragg","factoextra","cluster","kableExtra","skimr","WeightedCluster")
```

## Load Data

```{r load, warning=FALSE,message=FALSE}
cri <- read.csv(file = here::here("data/cri.csv"), header = TRUE)
cri <- subset(cri, u_teamid != 0)
cri_str <- read.csv(file = here::here("data/cri_str.csv"), header = TRUE)
cri_str <- subset(cri_str, u_teamid!=0)
cri_team <- read.csv(file = here::here("data/cri_team.csv"), header = TRUE)
cri_team <- subset(cri_team, u_teamid!=0)
```

## Specification Curves

### Data Prep

```{r data_prep, warning=FALSE, message=FALSE, include=FALSE}

# prep data
crispectest <- dplyr::select(cri_str, id, u_teamid, dv_type, iv_type, software, indepv, mator, dichotomize, twowayfe, cluster_any, Jobs:Scale,  w1985:w2016, p, AME_Z, lower_Z, upper_Z, id, eeurope, allavailable, mlm_any, Hresult, orig13, BELIEF_HYPOTHESIS, PRO_IMMIGRANT, TOPIC_KNOWLEDGE, STATISTICS_SKILL, MODEL_SCORE, main_IV_measurement, emplrate_ivC:gdp_ivC, inv_weight)

# create background degree variable (1 = other, 2 = methods focus, 3 = soc, 4 = polsci, 5 = psych)
crimrg <- cri %>%
  mutate(backgr_degree1 = ifelse(is.na(backgr_degree1), 6, backgr_degree1), # replace NA
         backgr_degree2 = ifelse(is.na(backgr_degree2), 6, backgr_degree2),
         backgr_degree3 = ifelse(is.na(backgr_degree3), 6, backgr_degree3),
         degree = car::recode(backgr_degree1, "1 = 1; 2 = 1; 6 = 1; 7 = 2; 3 = 3; 4 = 4; 5 = 5"),
         degree = ifelse(backgr_degree2 == 7, 2, ifelse(backgr_degree3 == 7, 2, degree)),
         # if any methods, then make methods for team
         )

crimrg <- dplyr::select(crimrg, degree, id)

crispectest <- left_join(crispectest, crimrg, by = "id")

# create sample variables 
crispectest <- crispectest %>%
  mutate(w1985 = ifelse(is.na(w1985), "Other", w1985),
         w1990 = ifelse(is.na(w1990), "Other", w1990),
         w1996 = ifelse(is.na(w1996), "Other", w1996),
         w2006 = ifelse(is.na(w2006), "Other", w2006),
         w2016 = ifelse(is.na(w2016), "Other", w2016),
         orig13 = ifelse(is.na(orig13), "Other", orig13),
         eeurope = ifelse(is.na(eeurope), "Other", eeurope),
         allavailable = ifelse(is.na(allavailable), "Other", allavailable))

crispectest <- crispectest %>%
  mutate(SURVEY_WAVES = ifelse(w1985 == "Other" & w1990 == "Other" & w1996 == "w1996" & w2006 == "w2006" & w2016 == "Other", "1996 & 2006", ifelse(w1985 == "Other" & w1990 == "Other" & w1996 == "w1996" & w2006 == "w2006" & w2016 == "w2016", "1996, 2006 & 2016", ifelse(w1985 == "Other" & w1990 == "Other" & w1996 == "Other" & w2006 == "w2006" & w2016 == "w2016", "2006 & 2016", "Other"))),
         COUNTRIES = ifelse(orig13 == "orig13", "13 Rich Democracies Only", ifelse(allavailable == "allavailable", "All Available (~26)", ifelse(eeurope == "eeurope", "Includes Eastern/CE Europe", "Other"))))

# Remove NAs
crispectest <- subset(crispectest, !is.na(crispectest$AME_Z))
crispectest <- subset(crispectest, !is.na(crispectest$lower_Z))
  
  # some upper and lower bounds are exactly zero, adjust fix this
crispectest$lower_Z <- ifelse(crispectest$lower_Z > -0.00000001 & crispectest$lower_Z < 0.00000001, -0.0001, crispectest$lower_Z)
crispectest$upper_Z <- ifelse(crispectest$upper_Z > -0.00000001 & crispectest$upper_Z < 0.00000001, 0.0001, crispectest$upper_Z)

  # trim to have better plot range
  crispectest$lower_Z <- ifelse(crispectest$lower_Z < -0.75, -0.75, crispectest$lower_Z)
  crispectest$upper_Z <- ifelse(crispectest$upper_Z > 0.75, 0.75, crispectest$upper_Z)
  crispectest$upper_Z <- ifelse(crispectest$upper_Z < -0.75, -0.745, crispectest$upper_Z)
  crispectest$AME_Z <- ifelse(crispectest$AME_Z < -0.75, -0.747, crispectest$AME_Z)
  crispectest$AME_Z <- ifelse(crispectest$AME_Z > 0.75, 0.747, crispectest$AME_Z)
  
  # data prep
crispectest$est <- round(crispectest$AME_Z, 6)
crispectest$lb <- round(crispectest$lower_Z, 6)
crispectest$ub <- round(crispectest$upper_Z, 6)

# set up subjective conclusions at team level as 1 support, 0 not testable and -1 reject
cri_team <- cri_team %>%
  mutate(est = ifelse(Hsup == 1, 0.1, ifelse(Hrej == 1, -0.1, ifelse(Hno == 1, 0, NA))),
         ub = est+0.01,
         lb = est-0.01)

# import researcher aspects to crispectest
cri_merge <- dplyr::select(cri, id, belief_ipred, pro_immigrant, topic_ipred, stats_ipred, total_score)

crispectest <- left_join(crispectest, cri_merge, by = "id")

# remove original study
crispectest <- subset(crispectest, u_teamid != 0)
cri_team <- subset(cri_team, u_teamid !=0)

# order rows to make counter
crispectest <- crispectest[order(crispectest$est),]
crispectest$count <- 1:nrow(crispectest)
```

## Model-Level Effect Specifications

### Main Plot

Setup

```{r fig1_main_prep}
# create identifier of neg, ns, pos
crispectest <- crispectest %>%
  mutate(sig_group = as.factor(ifelse(est < 0 & ub > 0, NA, ifelse(est > 0 & lb < 0, NA, ifelse(est < 0 & ub < 0, 1, ifelse(est > 0 & lb > 0, 3, 2))))),
         est_sig = ifelse(!is.na(sig_group), est, NA),
         est_ns = ifelse(is.na(sig_group), est, NA),
         sig_group2 = as.factor(ifelse(is.na(sig_group), 2, sig_group)))

#compress values for better visualization
# 0-0.05 = 0-1
# 0.05-0.505 = 1-2

crispectest <- crispectest %>%
  mutate(est_ns_scl = abs(est_ns),
         est_ns_scl = ifelse(est_ns_scl > 0 & est_ns_scl <= 0.05, est_ns_scl/0.05, ifelse(est_ns_scl > 0.05 & est_ns_scl <=0.505, ((est_ns_scl - 0.05)/0.505)+1, est_ns_scl+1.505)),
         est_ns_scl = ifelse(est_ns > 0, est_ns_scl, -1*est_ns_scl),
         est_sig_scl = abs(est_sig),
         est_sig_scl = ifelse(est_sig_scl > 0 & est_sig_scl <= 0.05, est_sig_scl/0.05, ifelse(est_sig_scl > 0.05 & est_sig_scl <=0.505, ((est_sig_scl - 0.05)/0.505)+1, est_sig_scl+1.505)),
         est_sig_scl = ifelse(est_sig > 0, est_sig_scl, -1*est_sig_scl))
```

Generate weighted percentage of each effect

```{r fig1_pcts}
sig_neg <- 100*round(sum(crispectest$inv_weight[crispectest$sig_group2 == 1])/sum(crispectest$inv_weight),3)

sig_pos <- 100*round(sum(crispectest$inv_weight[crispectest$sig_group2 == 3])/sum(crispectest$inv_weight),3)

ns <- 100*round(sum(crispectest$inv_weight[crispectest$sig_group2 == 2])/sum(crispectest$inv_weight),3)
```

#### Save Original PNG version

```{r fig1_main}
# for plotting a circle
dat = data.frame(x=runif(1), y=runif(1))

agg_png(filename = here::here("results/Fig1.png"), width = 800, height = 500, res = 144)
ggplot(crispectest) +
  #geom_errorbar(aes(ymin = lb, ymax = ub), color = "thistle1") +
  geom_point(aes(x = count, y = est_ns_scl), color = "grey55", shape = "|", size = 2, show.legend =F) + 
  geom_point(aes(x = count, y = est_sig_scl, color = sig_group), shape = "|", size = 3) +
  scale_color_manual(values = c("#E6AB02","NA", "#7570B3"), labels = c("Negative","Not sig.","Positive"," ")) +
  #geom_vline(xintercept = 1214, color = "grey", linetype = "dashed") +
  #geom_vline(xintercept = 52, color = "grey", linetype = "dashed") +
  #ylim(-1.8, 1.8) +
  annotate(geom = "text", x = 340, y = 1.45, label = "NEGATIVE (95% CI)", color = "#E6AB02", fontface = "bold", size = 3) +
  annotate(geom = "text", x = 940, y = 1.45, label = "POSITIVE (95% CI)", color = "#7570B3", fontface = "bold", size = 3) +
  annotate(geom = "text", x = 640, y = 1.45, label = "NOT STAT.\nSIGNIFICANT", fontface = "bold", color = "grey55", size = 3) +
  #annotate(geom = "text", x = 680, y = -1.5, label = "grey indicates not sig.\ndifferent from zero", color = "grey55", size = 3) +
  geom_vline(xintercept = -30) +
  annotate(geom = "segment", x = -35, xend = -25, y = -1.1, yend = -1.1, color = "white", size = 2) +
  annotate(geom = "segment", x = -35, xend = -25, y = 1.1, yend = 1.1, color = "white", size = 2) +
  annotate(geom = "segment", x = -38, xend = -22, y = 1.12, yend = 1.18, color = "black", size = 0.5) +
  annotate(geom = "segment", x = -38, xend = -22, y = -1.18, yend = -1.12, color = "black", size = 0.5) +
  labs(color = "Effect at 95% CI", x = "Models Ordered by AME", y = "Average Marginal Effect (AME)") +
  geom_point(aes(x=340, y=-2), data=dat, size=18, shape=1, color="#E6AB02") +
  geom_point(aes(x=640, y=-2), data=dat, size=26, shape=1, color="grey55") +
  geom_point(aes(x=940, y=-2), data=dat, size=14, shape=1, color="#7570B3") +
  annotate("text", x=340, y=-2, label = paste(sig_neg), color="#E6AB02", fontface = "bold") +
  annotate("text", x=640, y=-2, label = paste(ns), color="grey55", fontface = "bold") +
  annotate("text", x=940, y=-2, label = paste(sig_pos), color="#7570B3", fontface = "bold") +
  annotate("text", x=640, y=-1.1, label = "Weighted % of Distribution:", color="grey20", size = 2.5, fontface = 3) +
  annotate("text", x=640, y=0.5, label = "Point Estimates:", color="grey20", size = 2.5, fontface = 3) +
  scale_y_continuous(labels = c("-0.50","-0.05","0.00","0.05","0.50"), breaks = c(-1.75, -1, 0, 1, 1.75), limits = c(-2.6,1.8)) +
  coord_cartesian(xlim = c(-15,1255), clip = "off") +
  scale_x_continuous(expand = c(0,15)) +
  theme_classic() +
  guides(color = guide_legend(override.aes = list(size=7, color=c("#E6AB02","grey55", "#7570B3","NA")))) +
  theme(
    legend.key.size = unit(2,"line"),
    legend.position = "none",
    #axis.title.x = element_blank()
    axis.line.y = element_blank(),
    )
dev.off()

knitr::include_graphics(here::here("results/Fig1.png"))
```

#### Save PNAS TIF version

```{r fig1_main_tif}
# for plotting a circle
dat = data.frame(x=runif(1), y=runif(1))

agg_tiff(filename = here::here("results/Fig1.tiff"), width = 800, height = 500, res = 144)
ggplot(crispectest) +
  #geom_errorbar(aes(ymin = lb, ymax = ub), color = "thistle1") +
  geom_point(aes(x = count, y = est_ns_scl), color = "grey55", shape = "|", size = 2, show.legend =F) + 
  geom_point(aes(x = count, y = est_sig_scl, color = sig_group), shape = "|", size = 3) +
  scale_color_manual(values = c("#E6AB02","NA", "#7570B3"), labels = c("Negative","Not sig.","Positive"," ")) +
  #geom_vline(xintercept = 1214, color = "grey", linetype = "dashed") +
  #geom_vline(xintercept = 52, color = "grey", linetype = "dashed") +
  #ylim(-1.8, 1.8) +
  annotate(geom = "text", x = 320, y = 1.44, label = " \nNEGATIVE\n ", color = "#E6AB02", fontface = "bold", size = 3.5) +
  annotate(geom = "text", x = 960, y = 1.44, label = " \nPOSITIVE\n ", color = "#7570B3", fontface = "bold", size = 3.5) +
  annotate(geom = "text", x = 640, y = 1.44, label = " \nINCLUDES\nZERO", fontface = "bold", color = "grey55", size = 3.5) +
  geom_vline(xintercept = -30) +
  annotate(geom = "segment", x = -35, xend = -25, y = -1.1, yend = -1.1, color = "white", size = 2) +
  annotate(geom = "segment", x = -35, xend = -25, y = 1.1, yend = 1.1, color = "white", size = 2) +
  annotate(geom = "segment", x = -38, xend = -22, y = 1.12, yend = 1.18, color = "black", size = 0.5) +
  annotate(geom = "segment", x = -38, xend = -22, y = -1.18, yend = -1.12, color = "black", size = 0.5) +
  labs(color = "Effect at 95% CI", x = "Models Ordered by AME", y = "           Average Marginal Effect (AME)") +
  geom_point(aes(x=320, y=-2), data=dat, size=18, shape=1, color="#E6AB02") +
  geom_point(aes(x=640, y=-2), data=dat, size=26, shape=1, color="grey55") +
  geom_point(aes(x=960, y=-2), data=dat, size=14, shape=1, color="#7570B3") +
  annotate("text", x=320, y=-2, label = paste(sig_neg), color="#E6AB02", fontface = "bold") +
  annotate("text", x=640, y=-2, label = paste(ns), color="grey55", fontface = "bold") +
  annotate("text", x=960, y=-2, label = paste(sig_pos), color="#7570B3", fontface = "bold") +
  annotate("text", x=640, y=-1.1, label = "Weighted % of Distribution", color="grey20", size = 2.5, fontface = 3) +
  annotate("text", x=640, y=0.42, label = "Point Estimates", color="grey20", size = 2.5, fontface = 3) +
  annotate("text", x=640, y=1.8, label = "95% Confidence Interval", color="grey20", size = 2.5, fontface = 3) +
  scale_y_continuous(labels = c("-0.50","-0.05","0.00","0.05","0.50"), breaks = c(-1.75, -1, 0, 1, 1.75), limits = c(-2.6,1.8)) +
  coord_cartesian(xlim = c(-15,1255), clip = "off") +
  scale_x_continuous(expand = c(0,15)) +
  theme_classic() +
  guides(color = guide_legend(override.aes = list(size=7, color=c("#E6AB02","grey55", "#7570B3","NA")))) +
  theme(
    legend.key.size = unit(2,"line"),
    legend.position = "none",
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10),
    axis.line.y = element_blank(),
    )
dev.off()

```
### Supplementary Plots

#### Fig 1_supplementary. 

This figure combines many spec curves. Two big ones with researcher characteristics, and several small ones with the sub-domains of specifications.

First two are ordered by effect size, the other by confidence interval minimum/maximum relative to the hypothesis.

#### Panel A. By Effect Size

```{r S6panelA}
g1 <- ggplot(crispectest) +
  #geom_errorbar(aes(ymin = lb, ymax = ub), color = "thistle1") +
  geom_point(aes(x = count, y = est_ns), color = "grey55", shape = "|", size = 2, show.legend =F) + 
  geom_point(aes(x = count, y = est_sig, color = sig_group), shape = "|", size = 3) +
  scale_color_manual(values = c("#E6AB02","grey55", "#7570B3"), labels = c("Negative","Not sig.","Positive"," ")) +
  ylim(-0.21, 0.21) +
  labs(color = "Effect at 95%CI", x = "Model", y = " \n \nStandardized\nEffect Size") +
  theme_classic() +
  guides(color = guide_legend(override.aes = list(size=7, color=c("#E6AB02","grey55", "#7570B3","NA")))) +
  theme(
    legend.key.size = unit(2,"line"),
    legend.position = "none",
    axis.title.x = element_blank())

# make an adjusted version for the smaller plots
g1a <- ggplot(crispectest) +
  #geom_errorbar(aes(ymin = lb, ymax = ub), color = "thistle1") +
  geom_point(aes(x = count, y = est_ns), color = "grey55", shape = "|", size = 3, show.legend =F) + 
  geom_point(aes(x = count, y = est_sig, color = sig_group), shape = "|", size = 4) +
  scale_color_manual(values = c("#E6AB02","grey55", "#7570B3"), labels = c("Negative","Not sig.","Positive"," ")) +
  ylim(-0.21, 0.21) +
  labs(color = "Effect at 95%CI", x = "by Standardized Effect Size", y = " \n \n \n ") +
  theme_classic() +
  guides(color = guide_legend(override.aes = list(size=7, color=c("#E6AB02","grey55", "#7570B3","NA")))) +
  theme(
    legend.position = "none",
    axis.text.x = element_blank())
```

##### Panel B. By Confidence Interval

```{r S6_panelB}
# create identifier of neg, ns, pos
crispectest <- crispectest %>%
  mutate(ci = ifelse(sig_group2 == 1, ub, ifelse(sig_group2 == 3, lb, ifelse(est <= 0 & sig_group2 == 2, ub, ifelse(est > 0 & sig_group2 == 2, lb, NA)))),
         ub_order = ifelse(est <= 0 & sig_group2 == 1, 0, ifelse(est <= 0 & sig_group2==2, 1, ifelse(est > 0 & sig_group2==2, 2, ifelse(est > 0 & sig_group2 == 3, 3, NA)))),
         )

# order rows to make counter

crispectest <- crispectest[order(crispectest$ub_order, crispectest$ci),]
crispectest$count2 <- 1:nrow(crispectest)

# make an estimate where error bars show only for significant at 95% results
crispectest <- crispectest %>%
  mutate(est2 = ifelse(sig_group2 == 2, 0, est))

g2 <- 
  ggplot(crispectest) +
  geom_errorbar(aes(x = count2, ymin = lb, ymax = ub), color = "grey55") +
  geom_point(aes(x = count2, y = est, color = sig_group2), size = 1) +
  scale_color_manual(values = c("#E6AB02","grey55", "#7570B3"), labels = c("Negative","Not sig.","Positive")) +
  coord_cartesian(ylim = c(-0.21,0.21)) +
  labs(color = "Effect\nat 95%CI", x = "Model", y = "Effect Range", size = 2) +
  annotate("text", x = 450, y = 0.13,label = "Effects < 0\nordered by\nupper bound", size = 2.5, color = "grey10") +
  annotate("text", x = 900, y = -0.13,label = "Effects > 0\nordered by\nlower bound", size = 2.5, color = "grey10") +
  geom_segment(aes(x = 600, y = 0.08, xend = 350, yend = 0.08), arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "grey10") +
  geom_vline(xintercept = 688) +
  geom_segment(aes(x = 776, y = -0.08, xend = 1000, yend = -0.08), arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "grey10") +
  geom_vline(xintercept = 688) +
  theme_classic() +
  theme(
    axis.title.x = element_blank()
  )

# make an adjusted version for the smaller plots
g2a <- ggplot(crispectest) +
  geom_errorbar(aes(x = count2, ymin = lb, ymax = ub), color = "grey55") +
  geom_point(aes(x = count2, y = est, color = sig_group2), size = 1) +
  scale_color_manual(values = c("#E6AB02","grey55", "#7570B3"), labels = c("Negative","Not sig.","Positive")) +
  coord_cartesian(ylim = c(-0.21,0.21)) +
  labs(color = "Effect\nat 95%CI", x = "by Confidence Interval Minimum", size = 2) +
  geom_vline(xintercept = 688) +
  theme_classic() +
  theme(
    legend.position = "none",
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank()
  )

```

#### Panels B & D. Researcher Characteristics

Heatmap Concept

We create a 'rolling average' of researcher characteristics across rows ordered by effect size. This is visually pleasing and helps identify any patterns (one tick per model for so many models is really hard visually with a continuous color scale).

```{r rolling_average}
# create rolling average for effect ordering
crispectest <- crispectest[order(crispectest$count),]

# mean replacement, then rolling average, replace ends of rolling average w/ reg scores
crispectest <- crispectest %>%
  mutate(belief_ipredm = scale(ifelse(is.na(belief_ipred), mean(belief_ipred, na.rm=T), belief_ipred)),
         pro_immigrantm = scale(ifelse(is.na(pro_immigrant), mean(pro_immigrant, na.rm=T), pro_immigrant)),
         topic_ipredm = scale(ifelse(is.na(topic_ipred), mean(topic_ipred, na.rm=T), topic_ipred)),
         stats_ipredm = scale(ifelse(is.na(stats_ipred), mean(stats_ipred, na.rm=T), stats_ipred)),
         total_scorem = scale(ifelse(is.na(total_score), mean(total_score, na.rm=T), total_score)),
         belief_roll = zoo::rollmean(belief_ipredm, k = 30, fill = NA),
         pro_roll = zoo::rollmean(pro_immigrantm, k = 30, fill = NA),
         topic_roll = zoo::rollmean(topic_ipredm, k = 30, fill = NA),
         stats_roll = zoo::rollmean(stats_ipredm, k = 30, fill = NA),
         total_roll = zoo::rollmean(total_scorem, k = 30, fill = NA),
         belief_roll = ifelse(is.na(belief_roll), mean(belief_ipredm[count>1230,]), ifelse(is.na(belief_roll) & count < 30, mean(belief_ipredm[count<40,]), belief_roll)),
         pro_roll = ifelse(is.na(pro_roll) & count > 1230, mean(pro_immigrantm[count>1230,]), ifelse(is.na(pro_roll) & count < 30, mean(pro_immigrantm[count < 40,]), pro_roll)),
         topic_roll = ifelse(is.na(topic_roll) & count > 1230, mean(topic_ipredm[count>1230,]), ifelse(is.na(topic_roll) & count < 30, mean(topic_ipredm[count<30,]), topic_roll)),
         stats_roll = ifelse(is.na(stats_roll) & count > 1230, mean(stats_ipredm[count>1215,]), ifelse(is.na(stats_roll) & count < 30, mean(stats_ipredm[count<30,]), stats_roll)),
         total_roll = ifelse(is.na(total_roll) & count > 1230, mean(total_scorem[count>1230,]), ifelse(is.na(total_roll) & count < 30, mean(total_scorem[count<30,]), total_roll))
                              )

# create rolling average for effect CI order
crispectest <- crispectest[order(crispectest$count2),]

# mean replacement, then rolling average then replace tails w/ mean
crispectest <- crispectest %>%
  mutate(belief_roll2 = zoo::rollmean(belief_ipredm, k = 30, fill = NA),
         pro_roll2 = zoo::rollmean(pro_immigrantm, k = 30, fill = NA),
         topic_roll2 = zoo::rollmean(topic_ipredm, k = 30, fill = NA),
         stats_roll2 = zoo::rollmean(stats_ipredm, k = 30, fill = NA),
         total_roll2 = zoo::rollmean(total_scorem, k = 30, fill = NA),         
         belief_roll2 = ifelse(is.na(belief_roll2), mean(belief_ipredm[count2>1230,]), ifelse(is.na(belief_roll2) & count2 < 30, mean(belief_ipredm[count2<30,]), belief_roll2)),
         pro_roll2 = ifelse(is.na(pro_roll2) & count2 > 1230, mean(pro_immigrantm[count2>1230,]), ifelse(is.na(pro_roll2) & count2 < 30, mean(pro_immigrantm[count2 < 30,]), pro_roll2)),
         topic_roll2 = ifelse(is.na(topic_roll2) & count2 > 1230, mean(topic_ipredm[count2>1230,]), ifelse(is.na(topic_roll2) & count2 < 30, mean(topic_ipredm[count2<30,]), topic_roll2)),
         stats_roll2 = ifelse(is.na(stats_roll2) & count2 > 1230, mean(stats_ipredm[count2>1215,]), ifelse(is.na(stats_roll2) & count2 < 30, mean(stats_ipredm[count2<30,]), stats_roll2)),
         total_roll2 = ifelse(is.na(total_roll2) & count2 > 1230, mean(total_scorem[count2>1230,]), ifelse(is.na(total_roll2) & count2 < 30, mean(total_scorem[count2<30,]), total_roll2))
                              )
labs_y <- c("Model Score","Stats-Skill","Topic\nKnowledge","Pro Immigrant","Belief Hyp. True")

g3 <- ggplot(crispectest) +
  geom_tile(aes(x = count, y = 5, fill = belief_roll), height = 0.5, width = 2.00) +
  geom_tile(aes(x = count, y = 4, fill = pro_roll), height = 0.5, width = 2.00) +
  geom_tile(aes(x = count, y = 3, fill = topic_roll), height = 0.5, width = 2.00) +
  geom_tile(aes(x = count, y = 2, fill = stats_roll), height = 0.5, width = 2.00) +
  geom_tile(aes(x = count, y = 1, fill = total_roll), height = 0.5, width = 2.00) +
  scale_y_continuous(breaks = c(1,2,3,4,5), labels=labs_y) +
    theme_classic() +
  xlab("Models ordered by Effect Size") +
  theme(
    legend.position = "none",
    axis.title.y = element_blank()
  )
labs_y2 <- c("Model","l","e","t","e")

g4 <- ggplot(crispectest) +
  geom_tile(aes(x = count2, y = 5, fill = belief_roll2), height = 0.5, width = 2.00) +
  geom_tile(aes(x = count2, y = 4, fill = pro_roll2), height = 0.5, width = 2.00) +
  geom_tile(aes(x = count2, y = 3, fill = topic_roll2), height = 0.5, width = 2.00) +
  geom_tile(aes(x = count2, y = 2, fill = stats_roll2), height = 0.5, width = 2.00) +
  geom_tile(aes(x = count2, y = 1, fill = total_roll2), height = 0.5, width = 2.00) +
  labs(fill = "Score") +
  scale_y_continuous(breaks = c(1,2,3,4,5), labels=labs_y2) +
  xlab("Models ordered by Minimum Confidence Boundary (95%CI)") +
    theme_classic() +
  theme(
    plot.margin=unit(c(0.1,0.9,0.1,0.3),"cm"),
    axis.text.y = element_text(color = "NA"),
    axis.title.y = element_blank(),
    legend.text = element_text(size = 6),
    legend.title = element_text(size = 8)
  )
```

#### Fig S6. 

```{r figS6, warning = FALSE}
agg_png(filename = here::here("results/FigS6.png"), res = 144, height = 900, width = 1700)
ggarrange(g1,g2,g3,g4, labels = c("A ","B ","C ","D "))
dev.off()

knitr::include_graphics(here::here("results/FigS6.png"))
```

#### Fig 1. Sub-plots

We include two smaller versions of the spec-curves with four sub-domains of specifications below. This time without a rolling average.

First we recode our variables of interest into 1/0s for plotting:

```{r recode_order}
crispectest <- crispectest %>%
  mutate(DV_Jobs = ifelse(is.na(Jobs), 0, 1), # Measurement
         DV_Unemp = ifelse(is.na(Unemp), 0, 1),
         DV_IncDiff = ifelse(is.na(IncDiff), 0, 1),
         DV_OldAge = ifelse(is.na(OldAge), 0, 1),
         DV_House = ifelse(is.na(House), 0, 1),
         DV_Health = ifelse(is.na(Health), 0, 1),
         DV_Scale = ifelse(is.na(Scale), 0, 1), 
         Test_Stock = ifelse(iv_type == "Stock", 1, 0),
         Test_Flow = ifelse(iv_type == "Flow", 1, 0),
         Test_ChangeFlow = ifelse(iv_type == "ChangeFlow", 1, 0),
         M_Gross = ifelse(main_IV_measurement == "Immigrant, foreign-born" | main_IV_measurement == "Refugee" | main_IV_measurement == "Non-Western Immigrant", 1, 0),
         M_Net = ifelse(main_IV_measurement == "Net Migration", 1, 0),
         W_Other = ifelse(SURVEY_WAVES == "Other", 1, 0), # Sample 
         W_96_06 = ifelse(SURVEY_WAVES == "1996 & 2006", 1, 0),
         W_96_06_16 = ifelse(SURVEY_WAVES == "1996, 2006 & 2016", 1, 0),
         W_06_16 = ifelse(SURVEY_WAVES == "2006 & 2016", 1, 0),
         C_13 = ifelse(COUNTRIES == "13 Rich Democracies Only", 1, 0),
         C_Other = ifelse(COUNTRIES == "Other", 1, 0),
         C_All = ifelse(COUNTRIES == "All Available (~26)", 1, 0),
         C_CEE = ifelse(COUNTRIES == "Includes Eastern/CE Europe", 1, 0),
         IV_Emp = ifelse(emplrate_ivC == "emplrate_ivC" | unemprate_ivC == "unemprate_ivC", 1, 0),
         IV_Emp = ifelse(is.na(IV_Emp), 0, IV_Emp),
         IV_GDP = ifelse(is.na(gdp_ivC), 0, 1),
         IV_Socx = ifelse(is.na(socx_ivC), 0, 1),
         IV_None = ifelse((IV_Emp + IV_GDP + IV_Socx) == 0, 1, 0),
         Est_ols = ifelse(mator == "ols", 1, 0),
         Est_logit = ifelse(mator == "logit", 1, 0),
         Est_mlglm = ifelse(mator == "ml_glm", 1, 0),
         Est_ologit = ifelse(mator == "ologit", 1, 0),
         Est_bayes = ifelse(mator == "bayes", 1, 0))
```

#### Measurement Specs

For each set of explanatory variables we plot two versions. One ordered by effect size (point estimate) and a second by confidence interval minimum. Four columns in total.

Test Format

```{r figS7d, warning = FALSE}
testlabs <- c("       Stock", "Flow", "Change\nFlow")
g01 <- 
  ggplot(crispectest) +
  geom_tile(aes(x = count, y = 0.75, fill = factor(Test_Stock)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 0.5, fill = factor(Test_Flow)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 0.25, fill = factor(Test_ChangeFlow)), height = 0.133, width = 2.00) +
  scale_y_continuous(breaks = c(0.75,0.5,0.25), labels=testlabs) +
  scale_fill_manual(values = c("white","blue")) +
  theme_classic() +
  ylab("Immigration Test\nOperationalization") +
  theme(
    axis.title.x = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
  )

g01ci <- 
  ggplot(crispectest) +
  geom_tile(aes(x = count2, y = 0.75, fill = factor(Test_Stock)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 0.5, fill = factor(Test_Flow)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 0.25, fill = factor(Test_ChangeFlow)), height = 0.133, width = 2.00) +
  scale_y_continuous(breaks = c(0.75,0.5,0.25), labels=testlabs) +
  scale_fill_manual(values = c("white","blue")) +
  theme_classic() +
  theme(
    axis.title.x = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank()
  )
```

Measurement of Immigration

```{r figS7c, warning = FALSE}
mlabs <- c("  Immigra-\nnts (%)", "Net\nMigration")
g02 <- 
  ggplot(crispectest) +
  geom_tile(aes(x = count, y = 0.5, fill = factor(M_Gross)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 0.25, fill = factor(M_Net)), height = 0.133, width = 2.00) +
  scale_y_continuous(breaks = c(0.5,0.25), labels=mlabs) +
  scale_fill_manual(values = c("white","darkblue")) +
  theme_classic() +
  ylab("Immigration Test\nMeasurement") +
  theme(
    axis.title.x = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
  )

g02ci <- 
  ggplot(crispectest) +
  geom_tile(aes(x = count2, y = 0.5, fill = factor(M_Gross)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 0.25, fill = factor(M_Net)), height = 0.133, width = 2.00) +
  scale_y_continuous(breaks = c(0.5,0.25), labels=mlabs) +
  scale_fill_manual(values = c("white","darkblue")) +
  theme_classic() +
  theme(
    axis.title.x = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank()
  )
```

Dependent Variables

```{r figS7b, warning = FALSE}
# set labels
dvlabs = c("             Jobs","Unemp","IncDiff","OldAge","House","Health","Scale")

g03 <- 
  ggplot(crispectest) +
  geom_tile(aes(x = count, y = 1.75, fill = factor(DV_Jobs)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 1.5, fill = factor(DV_Unemp)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 1.25, fill = factor(DV_IncDiff)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 1, fill = factor(DV_OldAge)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 0.75, fill = factor(DV_House)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 0.5, fill = factor(DV_Health)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 0.25, fill = factor(DV_Scale)), height = 0.133, width = 2.00) +
  scale_y_continuous(breaks = c(1.75,1.5,1.25,1,0.75,0.5,0.25), labels=dvlabs) +
  scale_fill_manual(values = c("white","blue")) +
    theme_classic() +
    ylab("Dependent Variable") +
  theme(
    axis.title.x = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
  )
  
g03ci <- 
  ggplot(crispectest) +
  geom_tile(aes(x = count2, y = 1.75, fill = factor(DV_Jobs)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 1.5, fill = factor(DV_Unemp)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 1.25, fill = factor(DV_IncDiff)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 1, fill = factor(DV_OldAge)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 0.75, fill = factor(DV_House)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 0.5, fill = factor(DV_Health)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 0.25, fill = factor(DV_Scale)), height = 0.133, width = 2.00) +
  scale_y_continuous(breaks = c(1.75,1.5,1.25,1,0.75,0.5,0.25), labels=dvlabs) +
  scale_fill_manual(values = c("white","blue")) +
    theme_classic() +
  theme(
    axis.title.x = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank()
  )
```

Independent Variables

```{r figS7a, warning = FALSE}
ivlabs <- c(" (Un)Emp\nRate","GDP\nper cap", "Social\nSpending","None\n/other")

g04 <- ggplot(crispectest) +
  geom_tile(aes(x = count, y = 1, fill = factor(IV_Emp)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 0.75, fill = factor(IV_GDP)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 0.5, fill = factor(IV_Socx)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count, y = 0.25, fill = factor(IV_None)), height = 0.133, width = 2.00) +
  scale_y_continuous(breaks = c(1,0.75,0.5,0.25), labels=ivlabs) +
  scale_fill_manual(values = c("white","darkblue")) +
    theme_classic() +
  ylab("Other Independent\nVariables in Model") +
  xlab("Models ordered by Effect Size") +
  theme(
    legend.position = "none")

# version for right-side column 
g04ci <- ggplot(crispectest) +
  geom_tile(aes(x = count2, y = 1, fill = factor(IV_Emp)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 0.75, fill = factor(IV_GDP)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 0.5, fill = factor(IV_Socx)), height = 0.133, width = 2.00) +
  geom_tile(aes(x = count2, y = 0.25, fill = factor(IV_None)), height = 0.133, width = 2.00) +
  scale_y_continuous(breaks = c(1,0.75,0.5,0.25), labels=ivlabs) +
  scale_fill_manual(values = c("white","darkblue")) +
    theme_classic() +
  xlab("Models ordered by 95%CI Minimum") +
  theme(
    legend.position = "none",
    axis.title.y = element_blank(),
    axis.text.y = element_blank())
```

#### FigS7

```{r figS7, warning = FALSE}
agg_png(filename = here::here("results/FigS7.png"), height = 1600, width = 1600, res = 144)
ggarrange(g1a, g2a, g01, g01ci, g02, g02ci, g03, g03ci, g04, g04ci, nrow = 5, ncol = 2, heights = c(0.6,.75,0.5,1.75,1.2), widths = c(1,0.9))
dev.off()

knitr::include_graphics(here::here("results/FigS7.png"))
```

```{r figS8c, warning = FALSE}
samplabs <- c("1996 &\n2006","2006 &\n2016","1996, 2006\n& 2016","Other")

 g9 <-
ggplot(crispectest) +
  geom_tile(aes(x = count, y = 4, fill = factor(W_96_06)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 3, fill = factor(W_96_06_16)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 2, fill = factor(W_06_16)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 1, fill = factor(W_Other)), height = 0.5, width = 2) +
  scale_y_continuous(breaks = c(4,3,2,1), labels=samplabs) +
  scale_fill_manual(values = c("NA","blue")) +
    theme_classic() +
  ylab("Survey Waves in Sample") +
  theme(
    axis.title.x = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
  )
 
 g9ci <-
ggplot(crispectest) +
  geom_tile(aes(x = count2, y = 4, fill = factor(W_96_06)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 3, fill = factor(W_96_06_16)), height = 0.5, width = 2) +
  geom_tile(aes(x = count2, y = 2, fill = factor(W_06_16)), height = 0.5, width = 2) +
  geom_tile(aes(x = count2, y = 1, fill = factor(W_Other)), height = 0.5, width = 2) +
  scale_y_continuous(breaks = c(4,3,2,1), labels=samplabs) +
  scale_fill_manual(values = c("NA","blue")) +
    theme_classic() +
  ylab("Survey Waves in Sample") +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
  )
```

```{r figS8b, warning = FALSE}
clabs <- c("Rich 13"," CEE/East\nEurope","All\nAvailable","Other")

g10 <-
ggplot(crispectest) +
  geom_tile(aes(x = count, y = 4, fill = factor(C_13)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 3, fill = factor(C_CEE)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 2, fill = factor(C_All)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 1, fill = factor(C_Other)), height = 0.5, width = 2) +
  scale_y_continuous(breaks = c(4,3,2,1), labels=clabs) +
  scale_fill_manual(values = c("NA","darkblue")) +
    theme_classic() +
  ylab("Countries in Sample") +
  theme(
    axis.title.x = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
  )

g10ci <-
ggplot(crispectest) +
  geom_tile(aes(x = count2, y = 4, fill = factor(C_13)), height = 0.5, width = 2) +
  geom_tile(aes(x = count2, y = 3, fill = factor(C_CEE)), height = 0.5, width = 2) +
  geom_tile(aes(x = count2, y = 2, fill = factor(C_All)), height = 0.5, width = 2) +
  geom_tile(aes(x = count2, y = 1, fill = factor(C_Other)), height = 0.5, width = 2) +
  scale_y_continuous(breaks = c(4,3,2,1), labels=clabs) +
  scale_fill_manual(values = c("NA","darkblue")) +
    theme_classic() +
  ylab("Countries in Sample") +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
  )
```

```{r figS8a, warning = FALSE}
estlabs <- c("OLS","Logit"," ML/GLM","Ordered\nlogit","Bayes")

g11 <-
ggplot(crispectest) +
  geom_tile(aes(x = count, y = 5, fill = factor(Est_ols)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 4, fill = factor(Est_logit)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 3, fill = factor(Est_mlglm)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 2, fill = factor(Est_ologit)), height = 0.5, width = 2) +
  geom_tile(aes(x = count, y = 1, fill = factor(Est_bayes)), height = 0.5, width = 2) +
  scale_y_continuous(breaks = c(5,4,3,2,1), labels=estlabs) +
  scale_fill_manual(values = c("NA","blue")) +
    theme_classic() +
  ylab("Estimation Method") +
  xlab("Models ordered by Effect Size") +
  theme(
    #axis.title.x = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
  )

g11ci <-
ggplot(crispectest) +
  geom_tile(aes(x = count2, y = 5, fill = factor(Est_ols)), height = 0.5, width = 2) +
  geom_tile(aes(x = count2, y = 4, fill = factor(Est_logit)), height = 0.5, width = 2) +
  geom_tile(aes(x = count2, y = 3, fill = factor(Est_mlglm)), height = 0.5, width = 2) +
  geom_tile(aes(x = count2, y = 2, fill = factor(Est_ologit)), height = 0.5, width = 2) +
  geom_tile(aes(x = count2, y = 1, fill = factor(Est_bayes)), height = 0.5, width = 2) +
  scale_y_continuous(breaks = c(5,4,3,2,1), labels=estlabs) +
  scale_fill_manual(values = c("NA","blue")) +
    theme_classic() +
  ylab("Estimation Method") +
  xlab("Models ordered by 95%CI Minimum") +
  theme(
    axis.title.y = element_blank(),
    legend.position = "none",
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
  )
```

#### FigS8

```{r figS8, warning = FALSE}
agg_png(filename = here::here("results/FigS8.png"), height = 1200, width = 1600, res = 144)
ggarrange(g1a, g2a, g9, g9ci, g10, g10ci, g11, g11ci, nrow = 4, ncol = 2, heights = c(0.6,1,1,1.2), widths = c(1,0.9))
dev.off()

knitr::include_graphics(here::here("results/FigS8.png"))
```

### Fig S9. Team-Level Supplementary Plot

Order results by subjective conclusion and then percent ratio of effect sizes at p05 neg/pos.

```{r figS9a, warning = FALSE}
cri_team <- cri_team %>%
  mutate(subj_concl = factor(Hresult, levels = c("Reject","No test","Support")),
         pct_negpos = ifelse(subj_concl == "Reject", AME_neg_p05, 
                             ifelse(subj_concl == "No test", AME_ns_p05, AME_sup_p05)))
  
# order rows to make counter (by negpos)
cri_team <- cri_team[order(cri_team$subj_concl, cri_team$pct_negpos),]
cri_team$count <- 1:nrow(cri_team)

# alternative counter (by average effect size within conclusions)
cri_team <- cri_team[order(cri_team$subj_concl, cri_team$AME_Z),]
cri_team$count2 <- 1:nrow(cri_team)

# create trimmed AME for plotting
cri_team <- cri_team %>%
  mutate(AME_Zt = ifelse(AME_Z < -0.03, -0.03, ifelse(AME_Z > 0.02, 0.02, AME_Z)))

cri_team <- cri_team %>%
  mutate(belief_ipredm = scale(ifelse(is.na(belief_ipred), mean(belief_ipred, na.rm=T), belief_ipred)),
         pro_immigrantm = scale(ifelse(is.na(pro_immigrant), mean(pro_immigrant, na.rm=T), pro_immigrant)),
         topic_ipredm = scale(ifelse(is.na(topic_ipred), mean(topic_ipred, na.rm=T), topic_ipred)),
         stats_ipredm = scale(ifelse(is.na(stats_ipred), mean(stats_ipred, na.rm=T), stats_ipred)),
         total_scorem = scale(ifelse(is.na(total_score), mean(total_score, na.rm=T), total_score)),
         percentpos = scale(pos_test_pct_p05),
         belief_ipredm = ifelse(belief_ipredm > 2, 2, ifelse(belief_ipredm < -2, -2, belief_ipredm)),
         pro_immigrantm = ifelse(pro_immigrantm > 2, 2, ifelse(pro_immigrantm < -2, -2, pro_immigrantm)),
         topic_ipredm = ifelse(topic_ipredm > 2, 2, ifelse(topic_ipredm < -2, -2, topic_ipredm)),
         stats_ipredm = ifelse(stats_ipredm > 2, 2, ifelse(stats_ipredm < -2, -2, belief_ipredm)),
         total_scorem = ifelse(total_scorem > 2, 2, ifelse(total_scorem < -2, -2, total_scorem)),
         percentpos = ifelse(percentpos > 2, 2, ifelse(percentpos < -2, -2, percentpos)))

labs_y <- c("Model Score","Stats-Skill","Topic\nKnowledge","Pro Immigrant","Belief Hyp. True") 

labs_y_new <- c("Model Score","Stats-Skill","Topic\nKnowledge","Pro Immigrant","Belief Hyp. True", "% Positive\nTest Results")
```

#### Figure S9

```{r figS9, warning = FALSE}
g5a_new <- 
  ggplot(cri_team) +
  geom_point(aes(x = count2, y = AME_Zt, color = subj_concl), shape = 18, size = 4) +   scale_color_manual(values = c("#7570B3", "grey55", "#E6AB02"), labels = c("Rejected","Not testable","Supported")) +
  labs(color = "Conclusion\nHypothesis is:", x = "Team Conclusions", y = "Average Effect Size\nof Team's Test Models") +
  theme_classic() +
  theme( 
    plot.margin=unit(c(0.2,0.2,0.2,1.05),"cm"),
    axis.title.y = element_text(size = 12),
    axis.title.x = element_text(size = 12),
    axis.text.x = element_blank(),
  )

g5b_new <- ggplot(cri_team) +
  geom_tile(aes(x = count2, y = 6, fill = percentpos), height = 0.5) +
  geom_tile(aes(x = count2, y = 5, fill = belief_ipredm), height = 0.5) +
  geom_tile(aes(x = count2, y = 4, fill = pro_immigrantm), height = 0.5) +
  geom_tile(aes(x = count2, y = 3, fill = topic_ipredm), height = 0.5) +
  geom_tile(aes(x = count2, y = 2, fill = stats_ipredm), height = 0.5) +
  geom_tile(aes(x = count2, y = 1, fill = total_scorem), height = 0.5) +
  scale_y_continuous(breaks = c(1,2,3,4,5,6), labels=labs_y_new) +
    theme_classic() +
  labs(fill = "Standardized\nScore") +
  theme(
    #legend.position = "none",
    axis.title.y = element_blank(),
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    plot.margin=unit(c(0.2,0.3,0.2,0.5),"cm")
  )
agg_png(filename = here::here("results/FigS9.png"), res = 144, height = 900, width = 1200)
ggarrange(g5a_new,g5b_new,ncol = 1, heights = c(1,1.2))
dev.off()

knitr::include_graphics(here::here("results/FigS9.png"))
```

#### Note for FigS9

We find that `r round(length(cri_team$pos_test_pct_p05[cri_team$pos_test_pct_p05>0.6])/88,3)` percent of teams had a majority of positive test results (a negative effect of immigration at p<0.05) and `r round(length(cri_team$neg_test_pct_p05[cri_team$neg_test_pct_p05>0.6])/88,3)` percent had a majority negative.

```{r saveout}
write.csv(crispectest, file = here::here("results/crispectest.csv"))
```


## Colophon

This file is part of [https://github.com/nbreznau/CRI](https://github.com/nbreznau/CRI), the reproduction materials for [*Observing Many Researchers using the Same Data and Hypothesis Reveals a Hidden Universe of Uncertainty*](https://doi.org/10.31222/osf.io/cd5j9).

```{r colophon, echo=FALSE}
sessionInfo()
```