Multi-project assessment of sample quality in cross-national surveys

Replication materials

library(tidyverse) # for manipulating data
library(labelled) # for using labels
library(flextable) # for formatting output tables
library(readr) # for reading csv files
library(broom) # for dealing with annoying output chunks
library(scales) # for percentage labelling
library(countrycode) # for working with country names
library(ggrepel) # for labels
library(ggthemes) # for plot formatting
library(ggpubr) # for plot arrangements

SUPLEMENTARY MATERIALS

As accompanying material to “Multi-project assessments of sample quality in cross-national surveys: The role of weights in applying external and internal measures of sample bias”, this document presents details of data transformations, statistical procedures, and exact results of analyses performed in the paper’s main body. Code chunks are hidden by default but can be accessed by clicking the relevant buttons. All calculations rely on the publicly available integrated .csv data file.

ESM2e02 <- read_delim("ESM2e02f.csv", ";",
                      escape_double = FALSE, 
                      na = c("Not applicable", "Not available"), 
                      trim_ws = TRUE, 
                      col_types = cols(
  .default = col_character(),
  T_SURVEY_YEAR = col_double(),
  Q_PROP_FEM_POPULATION_18_74 = col_double(),
  Q_EFFECT_SIZE_PROP_FEMALE_EXTERNAL_18_74_psweights = col_double(),
  Q_EFFECT_SIZE_PROP_FEMALE_EXTERNAL_18_74_dweights_ESS = col_double(),
  Q_EFFECT_SIZE_PROP_FEMALE_EXTERNAL_18_74_no_weights = col_double(),
  Q_PROP_FEM_INTERNAL_KOHLER_psweights_strict = col_double(),
  Q_PROP_FEM_INTERNAL_KOHLER_no_weights_strict = col_double(),
  Q_SUBSAMPLE_SIZE_INTERNAL_KOHLER_strict = col_double(),
  Q_EFFECT_SIZE_PROP_FEMALE_INTERNAL_KOHLER_psweights_strict = col_double(),
  Q_EFFECT_SIZE_PROP_FEMALE_INTERNAL_KOHLER_no_weights_strict = col_double(),
  Q_ABS_BIAS_INTERNAL_KOHLER_psweights_strict = col_double(),
  Q_ABS_BIAS_INTERNAL_KOHLER_no_weights_strict = col_double(),
  Q_BIAS_INTERNAL_KOHLER_psweights_strict = col_double(),
  Q_BIAS_INTERNAL_KOHLER_no_weights_strict = col_double(),
  Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_dweights_ESS = col_double(),
  Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_no_weights = col_double(),
  Q_ABS_BIAS_INTERNAL_KOHLER_dweights_ESS = col_double(),
  Q_ABS_BIAS_INTERNAL_KOHLER_no_weights = col_double(),
  Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_psweights = col_double(),
  Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_no_weights = col_double(),
  Q_ABS_BIAS_INTERNAL_KOHLER_psweights = col_double(),
  Q_ABS_BIAS_INTERNAL_KOHLER_no_weights = col_double()
))

DISTRIBUTION OF THE TYPES OF WEIGHTS AVAILABLE ACROSS THE FOUR PROJECTS

Figure 1 - with embedded R code

ESM2e02  %>%
  mutate(T_Weighting_factor2 = factor(ESM2e02$T_Weighting_factor2, 
          levels = c("No weights", "Post-stratification weights", "Post-stratification weights with design weights"), 
          labels = c("No weights", "Post-stratification", "Post-stratification and design"))) %>%
  mutate(T_SURVEY_NAME = factor(ESM2e02$T_SURVEY_NAME, 
                                levels = c("EB", "ESS", "EQLS", "ISSP"), 
                                labels = c("EB", "ESS", "EQLS", "ISSP"))) %>%
  mutate(cntry_order = fct_infreq(as.factor(T_COUNTRY), ordered = FALSE)) %>%
  mutate(cntry_factor = dense_rank(cntry_order)) %>%
  ggplot() +
  geom_point(aes(x = ifelse(T_SURVEY_NAME == "EB", T_SURVEY_YEAR - 0.3,
                            ifelse(T_SURVEY_NAME == "ESS", T_SURVEY_YEAR - 0.1,
                                   ifelse(T_SURVEY_NAME == "EQLS", T_SURVEY_YEAR + 0.1, T_SURVEY_YEAR + 0.3))), 
                 y = cntry_factor, 
                 fill = T_Weighting_factor2,
                 color = T_Weighting_factor2,
                 shape = T_Weighting_factor2),
             show.legend = TRUE,
             size = 2.5,
             alpha = 0.8) +
  geom_text(aes(x = ifelse(T_SURVEY_NAME == "EB", T_SURVEY_YEAR - 0.3,
                           ifelse(T_SURVEY_NAME == "ESS", T_SURVEY_YEAR - 0.1,
                                  ifelse(T_SURVEY_NAME == "EQLS", T_SURVEY_YEAR + 0.1, T_SURVEY_YEAR + 0.3))),
                y = -0.9,
                label = ifelse(T_SURVEY_NAME == "EB", "EB",
                               ifelse(T_SURVEY_NAME == "ESS", "ESS",
                                      ifelse(T_SURVEY_NAME == "EQLS", "EQLS", "ISSP")))),
            size = 3.2,
            angle = 90) +
  scale_fill_manual(values = c("#F2CB3D", "#4F9DDB", "#5A1816")) +
  scale_color_manual(values = c("#000000", "#4F9DDB", "#5A1816")) +
  scale_x_continuous(breaks = c(2002:2016), 
                     labels = c("'02", "'03", "'04", "'05", "'06", "'07", "'08", "'09", "'10", "'11","'12", "'13", "'14", "'15","'16"), 
                     limits = c(2001.7, 2016.1), minor_breaks = seq(2001.5, 2016.5, 1)) +
  scale_y_continuous(breaks = c(1:40),
                     labels = c("Germany",  "Finland",  "France",   "United Kingdom",   "Sweden",   "Slovenia",     "Czechia",  "Denmark",  "Netherlands",  "Belgium",  "Spain",    "Hungary",  "Poland",   "Portugal",     "Slovakia",     "Austria",  "Ireland",  "Latvia",   "Bulgaria",     "Lithuania",    "Cyprus",   "Estonia",  "Croatia",  "Turkey",   "Italy",    "Greece",   "Norway",   "Switzerland",  "Luxembourg",   "Romania",  "Israel",   "Malta",    "Russia",   "Iceland",  "North Macedonia",  "Montenegro",   "Serbia",   "Ukraine",  "Albania",  "Georgia"), sec.axis = sec_axis(~.*1, breaks = c(1:40), labels = c("DE", "FI", "FR", "GB", "SE", "Sl", "CZ", "DK", "NL", "BE", "ES", "HU", "PL", "PT", "SK", "AT", "IE", "LV", "BG", "LT", "CY", "EE", "HR", "TR", "IT", "GR", "NO", "CH", "LU", "RO", "IL", "MT", "RU", "IS", "MK", "ME", "RS", "UA", "AL", "GE"))) +
  guides(shape = guide_legend("Available weights"), fill = guide_legend("Available weights"), color = guide_legend("Available weights")) +
  scale_shape_manual(values = c(21, 22, 23)) +
  coord_cartesian(xlim =  c(2002, 2016), ylim = c(-1, 39)) +
  theme_tufte() + 
  theme(text = element_text(family = "serif"), legend.position = "bottom", legend.text = element_text(size = 12, face = "bold", color = "black"),legend.title = element_text(size = 12, face = "bold", color = "black"), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.y = element_text(size = 12, face = "bold", margin = margin(0,0,0,0, "pt"), color = "black"), axis.text.x = element_text(size = 12, face = "bold", color = "black"), panel.grid.major.y = element_line(size = 0.1), panel.grid.minor.x = element_line(size = 0.3, linetype = "dashed"), plot.margin = unit(c(1,1,1,1),"mm"), panel.border = element_blank())

Raw counts and row percentages of data represented in the paper on Figure 1.

ESM2e02 %>%
  count(T_SURVEY_NAME, T_Weighting_factor2) %>%
  pivot_wider(names_from = T_SURVEY_NAME, values_from = n, values_fill = 0) %>%
  rename("Type of weight" = T_Weighting_factor2) %>%
  qflextable() %>%
  align(j = 2:5, align = "center") %>%
  align_nottext_col(align = "center")

ESM2e02 %>%
  count(T_SURVEY_NAME, T_Weighting_factor2) %>%
  group_by(T_SURVEY_NAME) %>%
  mutate(Percent = percent(n/sum(n), accuracy = 0.1)) %>% 
  select(T_SURVEY_NAME, T_Weighting_factor2, Percent) %>% 
  ungroup() %>%
  pivot_wider(names_from = T_SURVEY_NAME, values_from = Percent, values_fill = "0%") %>% 
  rename("Type of weight" = T_Weighting_factor2) %>%
  qflextable() %>% 
  theme_booktabs() %>%
  align(j = 1, align = "left")

STRICT VS LENIENT APPROACH TO INTERNAL CRITERIA: SEE SECTION 4.3.

Pearson correlation coefficients between biases resulting internal criteria without weights according to the strict approach vs the lenient approach. Calculations performed for surveys within the two projects where data for implementing the strict approach was available.

ESM2e02 %>% filter(T_SURVEY_NAME %in% c("ESS", "EQLS")) %>%
group_by(T_SURVEY_NAME) %>% 
summarize(`Perason correlation` = (cor(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, Q_ABS_BIAS_INTERNAL_KOHLER_no_weights_strict, method = 'pearson', use = "pairwise.complete.obs")), .groups = 'drop') %>% 
  rename("Project name" = T_SURVEY_NAME) %>%
qflextable() %>% 
theme_booktabs() %>% 
align_nottext_col(align = "center") %>% 
colformat_num(digits = 2, na_str = "N/A")

DIFFERENCES IN THE VALUES OF BIAS ACCORDING TO INTERNAL AND EXTERNAL CRITERIA: NO WEIGHT VS DESIGN WEIGHT

Figure 2 - with embedded R code

mean_vals <- ESM2e02 %>%
  filter(T_SURVEY_NAME == "ESS") %>%
  filter(!T_SURVEY_EDITION_COUNTRY %in% c("ESS4_SK", "ESS5_SK", "ESS6_SK", "ESS3_LV", "ESS3_RO", "ESS4_LT")) %>%
  filter(T_DEFFp_ESS == "DEFF>1") %>%
  group_by(T_SURVEY_EDITION) %>%
  summarise(int_no = mean(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, na.rm = TRUE), int_ds = mean(Q_ABS_BIAS_INTERNAL_KOHLER_dweights_ESS, na.rm = TRUE), ext_no = mean(Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_no_weights, na.rm = TRUE), ext_ds = mean(Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_dweights_ESS,na.rm = TRUE),  .groups = 'drop')

plot2a <- ESM2e02 %>%
  filter(T_SURVEY_NAME == "ESS") %>%
  filter(!T_SURVEY_EDITION_COUNTRY %in% c("ESS4_SK", "ESS5_SK", "ESS6_SK", "ESS3_LV", "ESS3_RO", "ESS4_LT")) %>%
  filter(T_DEFFp_ESS == "DEFF>1") %>%
  filter(T_External_and_internal_criteria != "Only external applicable") %>%
  ggplot() +
  ggtitle("Sample bias according to internal criteria") +
  ylab("") +
  scale_y_continuous(breaks = c(0, 1.96, 3, 6, 9, 12), 
                     labels = c(0, 1.96, 3, 6, 9, 12), 
                     limits = c(0, 12)) +
  geom_point(aes(x = 1,  y = Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, fill = as.factor(1), shape = as.factor(1)), 
             alpha = 0.5, show.legend = TRUE, size = 2.5) +
  geom_violin(aes(group = T_SURVEY_YEAR, x = 1, y = Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, fill = as.factor(1)), 
              alpha = 0.01, trim = TRUE, size = 0.8, na.rm = TRUE, draw_quantiles = 0.5, show.legend = TRUE) + 
  geom_point(aes(x = 2,  y = Q_ABS_BIAS_INTERNAL_KOHLER_dweights_ESS, fill = as.factor(2), shape = as.factor(2)), 
             alpha =  0.5, size = 2.5, show.legend = TRUE) + 
  geom_violin(aes(group = T_SURVEY_YEAR, x = 2, y = Q_ABS_BIAS_INTERNAL_KOHLER_dweights_ESS, fill = as.factor(2)),
              alpha = 0.01, trim = TRUE, size = 0.8, draw_quantiles = 0.5, show.legend = TRUE) +
  geom_hline(linetype = "twodash", yintercept = 1.96, alpha = 0.6) +
  scale_color_manual(values = c("#F2CB3D" , "#A6080B"), labels = c("Without Weights", "Design Weights")) +
  scale_fill_manual(values = c("#F2CB3D" , "#A6080B"), labels = c("Without Weights", "Design Weights")) +
  scale_shape_manual(values = c(21 , 24), labels = c("Without Weights", "Design Weights")) +
  geom_point(data = mean_vals, aes(x = 2, y = int_ds), shape = "X", size = 4, show.legend = FALSE) +
  geom_point(data = mean_vals, aes(x = 1, y = int_no), shape = "X", size = 4, show.legend = FALSE) +
  guides(shape = guide_legend("Median values indicated by horizontal line and the mean by X "), 
         colour = guide_legend("Median values indicated by horizontal line and the mean by X "), 
         fill = guide_legend("Median values indicated by horizontal line and the mean by X ")) +
  coord_cartesian(ylim = c(0.1, 9)) +
  theme_bw() + theme(axis.text.x = element_blank(), text = element_text(family = "serif"), 
                     strip.background = element_rect(fill = "white"), 
                     axis.title.x = element_blank(),
                     axis.ticks.x = element_blank(),
                     title = element_text(size=12, face = "bold", color = "black"), 
                     strip.text = element_text(size=12, face = "bold", color = "black"),
                     axis.text.y = element_text(size=12, face = "bold", color = "black"),
                     legend.position = "bottom",
                     legend.box = "vertical",
                     legend.justification = "center",
                     legend.text = element_text(size=12, face = "bold", color = "black"),
                     plot.margin = unit(c(1,1,1,1),"mm",),
                     legend.key.size = unit(1.2,"line")) +
  facet_wrap(~ T_SURVEY_EDITION, nrow = 1)

plot2b <- ESM2e02 %>%
  filter(T_SURVEY_NAME == "ESS") %>%
  filter(!T_SURVEY_EDITION_COUNTRY %in% c("ESS4_SK", "ESS5_SK", "ESS6_SK", "ESS3_LV", "ESS3_RO", "ESS4_LT")) %>%
  filter(T_DEFFp_ESS == "DEFF>1") %>%
  ggplot() +
  ggtitle("Sample bias according to external criteria") +
  ylab("") +
  scale_y_continuous(breaks = c(0, 1.96, 3, 6, 9), 
                     labels = c(0, 1.96, 3, 6, 9), 
                     limits = c(0, 9)) +
  geom_point(aes(x = 1,  y = Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_no_weights, fill = as.factor(1), shape = as.factor(1)), 
             alpha = 0.5, size = 2.5, show.legend = FALSE) +
  geom_violin(aes(group = T_SURVEY_YEAR, x = 1, y = Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_no_weights, fill = as.factor(1)), 
              alpha = 0.01, trim = TRUE, size = 0.8, draw_quantiles = 0.5, show.legend = FALSE) +
  geom_point(aes(x = 2,  y = Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_dweights_ESS, fill = as.factor(2), shape = as.factor(2)), 
             alpha = 0.5, size = 2.5, show.legend = FALSE) +
  geom_violin(aes(group = T_SURVEY_YEAR, x = 2, y = Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_dweights_ESS, fill = as.factor(2)), 
              alpha = 0.01, trim = TRUE, size = 0.8, draw_quantiles = 0.5, show.legend = FALSE) +
  geom_hline(linetype = "twodash", yintercept = 1.96, alpha = 0.6) +
  scale_color_manual(values = c("#F2CB3D" , "#A6080B"), labels = c("Without Weights", "Design Weights")) +
  scale_fill_manual(values = c("#F2CB3D" , "#A6080B"), labels = c("Without Weights", "Design Weights")) +
  scale_shape_manual(values = c(21 , 24), labels = c("Without Weights", "Design Weights")) +
  geom_point(data = mean_vals, aes(x = 2, y = ext_ds), shape = "X", size = 4, show.legend = FALSE) +
  geom_point(data = mean_vals, aes(x = 1, y = ext_no), shape = "X", size = 4, show.legend = FALSE) +
  coord_cartesian(ylim = c(0.1, 9)) +
  guides(colour = FALSE, shape = FALSE, fill = FALSE) +
  theme_bw() + theme(axis.text.x = element_blank(),
                     text = element_text(family = "serif"), 
                     strip.background = element_rect(fill = "white"),
                     axis.title.x = element_blank(),
                     axis.ticks.x = element_blank(),
                     title = element_text(size = 12, face = "bold", color = "black"), 
                     strip.text = element_text(size = 12, face = "bold", color = "black"),
                     axis.text.y = element_text(size = 12, face = "bold", color = "black"),
                     legend.position = "bottom",
                     legend.box = "vertical",
                     legend.justification = "left",
                     legend.text = element_text(size = 12, face = "bold", color = "black"),
                     plot.margin = unit(c(1,1,1,1),"mm")) + 
  facet_wrap(~ T_SURVEY_EDITION, nrow = 1)

ggarrange(plot2b, plot2a, nrow = 2)

Number of national surveys evaluated per edition for data presented on Figure 2.

ESM2e02 %>%
  filter(T_SURVEY_NAME == "ESS") %>%
  filter(!T_SURVEY_EDITION_COUNTRY %in% c("ESS4_SK", "ESS5_SK", "ESS6_SK", "ESS3_LV", "ESS3_RO", "ESS4_LT")) %>%
  filter(T_DEFFp_ESS == "DEFF>1") %>%
  count(T_SURVEY_EDITION) %>%
rename("Project Edition" = T_SURVEY_EDITION,
       "Number of national surveys \n under evaluation" = n) %>%
  qflextable() %>% 
  theme_booktabs() %>% 
  align_nottext_col(align = "center") 

Mean and median values for data represented on Figure 2.

ESM2e02 %>% filter(T_SURVEY_NAME == "ESS", T_DEFFp_ESS == "DEFF>1") %>%
  mutate(delta1 = Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_dweights_ESS  - Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_no_weights,
         delta2 = Q_ABS_BIAS_INTERNAL_KOHLER_dweights_ESS - Q_ABS_BIAS_INTERNAL_KOHLER_no_weights) %>%
  group_by(T_SURVEY_EDITION) %>%
  summarise(meanXdifference_external = mean(delta1),
            medianXdifference_external = median(delta1),
            meanXdifference_internal = mean(delta2),
            medianXdifference_internal = round(median(delta2), digits = 2), .groups = 'drop') %>%
  rename("ESS wave" = T_SURVEY_EDITION, 
         "Mean difference \n external" = meanXdifference_external,
         "Median difference \n external"  = medianXdifference_external,
         "Mean difference \n internal" = meanXdifference_internal,
         "Median difference \n internal" = medianXdifference_internal) %>%
  qflextable() %>% 
  theme_booktabs() %>%
  align_nottext_col(align = "center", header = T) %>%
  colformat_num(j = 2:5, digits = 1, na_str = "N/A")

Mean and median differences as referenced in section 5.1.

ESM2e02 %>% filter(T_SURVEY_NAME == "ESS", T_DEFFp_ESS == "DEFF>1") %>%
  mutate(delta1 = Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_dweights_ESS  - Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_no_weights,
         delta2 = Q_ABS_BIAS_INTERNAL_KOHLER_dweights_ESS - Q_ABS_BIAS_INTERNAL_KOHLER_no_weights) %>%
  summarise(meanXdifference_external = mean(delta1),
            medianXdifference_external = median(delta1),
            meanXdifference_internal = mean(delta2),
            medianXdifference_internal = round(median(delta2), digits = 2), .groups = 'drop') %>%
  pivot_longer(1:4, names_to = c("fun", "Difference"), values_to = "values", names_sep = "X") %>%
  pivot_wider(names_from = fun, values_from = values) %>%
  qflextable() %>% 
  theme_booktabs() %>%
  colformat_num(j = 2:3, digits = 1, na_str = "N/A")

WITHIN-PROJECT DIFFERENCES BETWEEN BIAS WITH AND WITHOUT POST-STRATIFICATION WEIGHTS

Figure 3 - with embedded R code
Note: In the paper, the internal and external criteria visualisations have been arranged on parallel panels in order to minimize overall figure size.

 ESM2e02 %>% 
  filter(!T_SURVEY_EDITION_COUNTRY %in% c("ISSP2007_NL", "ISSP2015_DK")) %>%
  filter(T_Weighting_factor == "Total weights present in dataset") %>%
  ggplot(aes(x = T_SURVEY_YEAR)) + 
  ggtitle("Sample bias according to external criteria") +
  ylab("") +
  xlab("") +
  scale_x_continuous(breaks = c(2002:2016), 
                     labels = c("'02", "'03", "'04", "'05", "'06", "'07", "'08", "'09", "'10", "'11","'12", "'13", "'14", "'15","'16"), 
                     limits = c(2001.3, 2016.7)) +
  scale_y_continuous(breaks = c(0, 1.96, 3, 6, 9, 12), 
                     labels = c(0, 1.96, 3, 6, 9, 12), 
                     limits = c(0, 13))  +
  geom_violin(aes(group = T_SURVEY_YEAR, x = T_SURVEY_YEAR - 0.2,  y = Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_no_weights, fill = as.factor(1)),
              alpha = 0.4,
              trim = TRUE,
              na.rm = T,
              size = 0.8,
              draw_quantiles = 0.5) +
  geom_violin(aes(group = T_SURVEY_YEAR, x = T_SURVEY_YEAR + 0.2, y = Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_psweights,
              fill = as.factor(2)),
              alpha = 0.8,
              trim = TRUE,
              na.rm = T,
              size = 0.8,
              draw_quantiles = 0.5) +
  scale_fill_manual(values = c("#F2CB3D" , "#162AF2"), labels = c("Without weights", "Post-stratification weights")) +
  geom_hline(linetype = "twodash", yintercept = 1.96, alpha = 0.6) +
  guides(shape = FALSE, fill = guide_legend(title = "Absolute bias per wave",
                                            override.aes = list(alpha = c(0.2, 8))), color = FALSE, alpha = FALSE) +
  theme_bw() + theme(text = element_text(family = "serif"), strip.background = element_rect(fill = "white"),
                     axis.text.x = element_text(size = 11, face = "bold", color = "black"),
                     axis.title.x = element_blank(),
                     axis.ticks.x = element_blank(),
                     title = element_text(size = 11, face = "bold", color = "black"), 
                     strip.text = element_text(size = 11, face = "bold", color = "black"),
                     axis.text.y = element_text(size = 11, face = "bold", color = "black"),
                     legend.position = "bottom",
                     legend.box = "vertical",
                     legend.justification = "left",
                     legend.text = element_text(size = 11, face = "bold", color = "black"),
                     plot.margin = unit(c(1,1,1,1),"mm")) +
  facet_wrap(~ T_SURVEY_NAME, nrow = 4)

 ESM2e02 %>%
  filter(T_External_and_internal_criteria %in% c("Both applicable", "Only external applicable")) %>%
  filter(T_Weighting_factor == "Total weights present in dataset") %>%
  ggplot(aes(x = T_SURVEY_YEAR)) + 
  ggtitle("Sample bias according to internal criteria") +
  ylab("") +
  xlab("") +
  scale_x_continuous(breaks = c(2002:2016), 
                     labels = c("'02", "'03", "'04", "'05", "'06", "'07", "'08", "'09", "'10", "'11","'12", "'13", "'14", "'15","'16"), 
                     limits = c(2001.3, 2016.7)) + 
  scale_y_continuous(breaks = c(0, 1.96, 3, 6, 9, 12), 
                     labels = c(0, 1.96, 3, 6, 9, 12), 
                     limits = c(0, 13))  +
  geom_violin(aes(group = T_SURVEY_YEAR, x = T_SURVEY_YEAR - 0.2,
                  y = Q_ABS_BIAS_INTERNAL_KOHLER_no_weights,  
                  fill = as.factor(1)),
              na.rm = T,
              alpha = 0.4, 
              trim = TRUE, 
              size = 0.8,
              draw_quantiles = 0.5) +
  geom_violin(aes(group = T_SURVEY_YEAR,
                  x = T_SURVEY_YEAR + 0.2,
                  y = Q_ABS_BIAS_INTERNAL_KOHLER_psweights,
                  fill = as.factor(2)),
              na.rm = T,
              alpha = 0.9, 
              trim = TRUE, 
              size = 0.8,
              draw_quantiles = 0.5) +
  scale_fill_manual(values = c("#F2CB3D" , "#162AF2"), labels = c("Without weights", "Post-stratification weights")) +
  geom_hline(linetype = "twodash", yintercept = 1.96, alpha = 0.6) +
  guides(shape = FALSE, fill = guide_legend(title = "Absolute bias per wave",
                                            override.aes = list(alpha = c(0.2, 8))), color = FALSE, alpha = FALSE) +
  theme_bw() + theme(text = element_text(family = "serif"), strip.background = element_rect(fill = "white"),
                     axis.text.x = element_text(size = 11, face = "bold", color = "black"),
                     axis.title.x = element_blank(),
                     axis.ticks.x = element_blank(),
                     title = element_text(size = 11, face = "bold", color = "black"), 
                     strip.text = element_text(size = 11, face = "bold", color = "black"),
                     axis.text.y = element_text(size = 11, face = "bold", color = "black"),
                     legend.position = "bottom",
                     legend.box = "vertical",
                     legend.justification = "left",
                     legend.text = element_text(size = 11, face = "bold", color = "black"),
                     plot.margin = unit(c(1,1,1,1),"mm")) +
  facet_wrap(~ T_SURVEY_NAME, nrow = 4)

Number of national surveys evaluated per edition for data presented on Figure 3.

ESM2e02 %>%
  filter(T_External_and_internal_criteria %in% c("Both applicable", "Only external applicable", "Only external applicable")) %>%
  filter(T_Weighting_factor == "Total weights present in dataset") %>%
  group_by(T_SURVEY_YEAR) %>%
count(T_SURVEY_NAME) %>%
  pivot_wider(names_from = T_SURVEY_NAME, values_from = n) %>%
  mutate(Year = as.character(T_SURVEY_YEAR)) %>%
  ungroup() %>%
  select(Year, EB, EQLS, ESS, ISSP) %>%
  qflextable() %>% 
  theme_booktabs() %>% 
  align_nottext_col(align = "center") 

See Table 2 in the paper.

ESM2e02 %>% filter(T_Weighting_factor == "Total weights present in dataset") %>% 
  mutate(delta1 = Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_psweights - Q_ABS_BIAS_PROP_FEMALE_EXTERNAL_no_weights,
         delta2 = Q_ABS_BIAS_INTERNAL_KOHLER_psweights - Q_ABS_BIAS_INTERNAL_KOHLER_no_weights) %>%
  group_by(T_SURVEY_NAME) %>%
  summarise(mean_external = mean(delta1),
            median_external = median(delta1),
            mean_internal = mean(delta2, na.rm = T),
            median_internal = median(delta2, na.rm = T), .groups = 'drop') %>%
  rename("Project name" = T_SURVEY_NAME,
         "Mean difference \n external" = mean_external,
         "Median difference \n external" = median_external,
         "Mean difference \n internal" = mean_internal,
         "Median difference \n internal" = median_internal) %>%
  qflextable() %>% 
  theme_booktabs() %>% 
  align_nottext_col(align = "center") %>%
  colformat_num(j = 2:3, big.mark = ",", digits = 3, na_str = "N/A")

INTERNAL CRITERIA WITH NO WEIGHTS: THE OUTLIERS

Descriptive statistics of data referenced in section 6.

ESM2e02 %>%
  group_by(T_SURVEY_NAME) %>%
  summarise(Q1 = quantile(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, probs = 0.25, na.rm = T),
            median = quantile(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, probs = 0.25, na.rm = T),
            Q3 = quantile(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, probs = 0.75, na.rm = T),
            IQR = IQR(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, na.rm = T),
            Outlier_threshold = quantile(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, probs = 0.75, na.rm = T) + 1.5 * IQR(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, na.rm = T), .groups = 'drop') %>%
  rename("Project name" = T_SURVEY_NAME, "Outlier \n threshold" = Outlier_threshold) %>%
  qflextable() %>% 
  theme_booktabs() %>% 
  align_nottext_col(align = "center") %>%
  colformat_num(j = 2:3, digits = 2, na_str = "N/A") %>%
  colformat_num(j = 4:6, digits = 3, na_str = "N/A")

Figure 4 - with embedded R code

ESM2e02 %>% 
  filter(T_External_and_internal_criteria %in% c("Both applicable", "Only internal applicable")) %>%
  ggplot(aes(x = T_SURVEY_NAME)) + 
  ylab("") +
  scale_y_continuous(breaks = c(0, 1.96, 3, 6, 9, 12), 
                     labels = c(0, 1.96,  3, 6, 9, 12), 
                     limits = c(0, 12)) +
  geom_boxplot(aes(y = Q_ABS_BIAS_INTERNAL_KOHLER_no_weights), size = 1, width = 0.5, notch = FALSE, outlier.alpha = 0) +
  geom_jitter(data = filter(ESM2e02, Outliers != "outliers"), aes(y = Q_ABS_BIAS_INTERNAL_KOHLER_no_weights), width = 0.1, alpha = 0.2, na.rm = T) +
  geom_label_repel(data = filter(ESM2e02, Outliers == "outliers"),
                   aes(y =  Q_ABS_BIAS_INTERNAL_KOHLER_no_weights,
                       label = T_SURVEY_LABEL,
                       fill = Outlier_fill),
                   label.size = unit(0.1, "mm"), 
                   label.padding = unit(0.5, "mm"),
                   box.padding = unit(0.01, "mm"),
                   label.r = unit(1, "mm"), 
                   color = "#A6080B",
                   direction = "x",
                   nudge_y = 0,
                   segment.alpha = 0, 
                   force = 0.6,
                   max.overlaps = 100) +
  scale_fill_manual(values =  c("white", "gray80")) +
  geom_hline(linetype = "twodash", yintercept = 1.96, alpha = 0.6) +
  coord_cartesian(ylim = c(0.4,11.6)) +
  theme_bw() + theme(text = element_text(family = "serif"),
                     axis.text.x = element_text(size = 12, face = "bold", color = "black"),
                     axis.title.x = element_blank(),
                     title = element_text(size = 12, face = "bold", color = "black"), 
                     strip.text = element_text(size = 12, face = "bold", color = "black"),
                     axis.text.y = element_text(size = 12, face = "bold", color = "black"),
                     legend.position = "none",
                     plot.margin = unit(c(1,1,1,1),"mm"))

Number of national surveys evaluated per edition for data presented on Figure 4.

ESM2e02 %>% 
  filter(T_External_and_internal_criteria %in% c("Both applicable", "Only internal applicable")) %>%
count(T_SURVEY_NAME) %>%
    rename("Project name" = T_SURVEY_NAME,
         "Number of national surveys \n under evaluation" = n) %>%
  qflextable() %>% 
  theme_booktabs() %>% 
  align_nottext_col(align = "center")

A list of outliers as represented on Figure 4.

ESM2e02 %>% 
  group_by(T_SURVEY_NAME) %>%
  filter(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights > quantile(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, probs = 0.75, na.rm = T) + 1.5 * IQR(Q_ABS_BIAS_INTERNAL_KOHLER_no_weights, na.rm = T)) %>%
  select(T_SURVEY_NAME, T_SURVEY_YEAR, T_COUNTRY, Q_ABS_BIAS_INTERNAL_KOHLER_no_weights) %>%
  rename("Project name" = T_SURVEY_NAME, "Year" = T_SURVEY_YEAR, "Country code" = T_COUNTRY,  "Internal bias" = Q_ABS_BIAS_INTERNAL_KOHLER_no_weights) %>%
  ungroup() %>%
  mutate(`Country name` = countrycode(`Country code`, origin = "iso2c", destination = "country.name")) %>%
  relocate(`Project name`, `Year`, `Country code`, `Country name`, `Internal bias`) %>%
  arrange(`Project name`, `Country name`) %>%
  qflextable() %>% 
  theme_booktabs() %>%
  align(j = 2:4, align = "center") %>%
  align_nottext_col(align = "center", header = T) %>%
  colformat_num(digits = 3, big.mark = "", na_str = "N/A")