data extraction review

library(tidyverse)

## Warning: 패키지 'tidyverse'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'ggplot2'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'tibble'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'readr'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'purrr'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'dplyr'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'lubridate'는 R 버전 4.2.3에서 작성되었습니다

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── 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

library(readxl)

## Warning: 패키지 'readxl'는 R 버전 4.2.3에서 작성되었습니다

library(tidymodels)

## Warning: 패키지 'tidymodels'는 R 버전 4.2.3에서 작성되었습니다

## ── Attaching packages ────────────────────────────────────── tidymodels 1.1.1 ──
## ✔ broom        1.0.5     ✔ rsample      1.2.0
## ✔ dials        1.2.0     ✔ tune         1.1.2
## ✔ infer        1.0.4     ✔ workflows    1.1.3
## ✔ modeldata    1.2.0     ✔ workflowsets 1.0.1
## ✔ parsnip      1.1.1     ✔ yardstick    1.2.0
## ✔ recipes      1.0.8

## Warning: 패키지 'broom'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'dials'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'modeldata'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'parsnip'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'recipes'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'rsample'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'tune'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'workflows'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'workflowsets'는 R 버전 4.2.3에서 작성되었습니다

## Warning: 패키지 'yardstick'는 R 버전 4.2.3에서 작성되었습니다

## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ scales::discard() masks purrr::discard()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ recipes::fixed()  masks stringr::fixed()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ yardstick::spec() masks readr::spec()
## ✖ recipes::step()   masks stats::step()
## • Search for functions across packages at https://www.tidymodels.org/find/

theme_set(theme_light())

rm(list = ls())

QEX

folder_path <- "d:/r/data extraction process/"

data <- 
  tibble(filename = list.files(path = folder_path, full.names = TRUE, pattern = "\\.xlsx$")) %>%
  mutate(data = map(filename, ~read_excel(.x, sheet = "DATA EXTRACTION SHEET", skip = 2))) %>%
  unnest(data) %>%
  select(-c(n_T_revised:'ROB Category')) %>%
  filter(!is.na(`Coder name`)) %>%
  filter(!`Coder name` == "Record your name") %>%
  janitor::clean_names() %>% 
  mutate(batch = str_extract(filename, "(?<=Group )\\d+")) %>%
  select(-filename) %>%
  group_by(batch) %>%
  arrange(coder_name) %>%
  mutate(coders = paste0("coder_", dense_rank(coder_name))) %>%
  ungroup() %>%
  select(-coder_name)

data_long <- data %>%
  pivot_longer(
    cols = -c(batch, coders, estimate_id), 
    names_to = "variable", 
    values_to = "value"
  ) %>%
  mutate(value = na_if(value, "NA"),
         value = na_if(value, "N/A"))

data_pivoted <- data_long %>%
  pivot_wider(
    names_from = coders,
    values_from = value
  ) %>%
    mutate(matched = case_when(
      is.na(coder_1) & is.na(coder_2) ~ TRUE,
      is.na(coder_1) & !is.na(coder_2) ~ FALSE,
      !is.na(coder_1) & is.na(coder_2) ~ FALSE,
      coder_1 == coder_2 ~ TRUE,
      TRUE ~ FALSE))

each question level

data_pivoted %>%
  count(batch, estimate_id, variable, matched) %>%
  mutate(matched = ifelse(matched, "count_true", "count_false")) %>%
  pivot_wider(
    names_from = matched,
    values_from = n,
    values_fill = 0 # Fill in 0 if there are missing values
  ) %>%
  count(variable, count_true) %>%
  filter(count_true == 1) %>%
  mutate(pct = n / max(n)) %>%
  filter(pct != 1) %>%
  mutate(variable = fct_reorder(variable, pct)) %>%
  ggplot(aes(pct, variable)) +
  geom_col() +
  scale_x_continuous(label = scales::percent_format()) +
  labs(y = "",
       x = "between-coder reliability")

each estimate_id level

each_estimate_id <-
  data_pivoted %>%
  group_by(estimate_id, batch) %>%
  summarize(n = n(),
            agreement_rate = sum(matched == 1)/n,
            .groups = "drop")

each_estimate_id

## # A tibble: 267 × 4
##    estimate_id batch     n agreement_rate
##    <chr>       <chr> <int>          <dbl>
##  1 73628889_1  48       90          0.789
##  2 73628889_2  48       90          0.789
##  3 73795904_1  47       90          0.756
##  4 73795904_10 47       90          0.744
##  5 73795904_11 47       90          0.711
##  6 73795904_12 47       90          0.711
##  7 73795904_13 47       90          0.711
##  8 73795904_14 47       90          0.711
##  9 73795904_15 47       90          0.7  
## 10 73795904_16 47       90          0.7  
## # ℹ 257 more rows

study level

study_level <-
  each_estimate_id %>%
  mutate(id = str_extract(estimate_id, "\\d+(?=_)")) %>%
  group_by(id, batch) %>%
  summarize(mean = mean(agreement_rate),
            .groups = "drop")

study_level

## # A tibble: 44 × 3
##    id       batch   mean
##    <chr>    <chr>  <dbl>
##  1 73628889 48    0.789 
##  2 73795904 47    0.723 
##  3 73796232 47    0.8   
##  4 73797297 47    0.789 
##  5 73797942 51    0.8   
##  6 73798224 48    0.0111
##  7 73798778 51    0.8   
##  8 73801283 48    0.756 
##  9 73801284 48    0.5   
## 10 73801962 52    0.833 
## # ℹ 34 more rows

• These study level data can be merged with the study feature such as evaluation design, the number of outcomes and others

the number of outcomes in each study and between-coders reliability

each_estimate_id %>%
  mutate(id = str_extract(estimate_id, "\\d+(?=_)")) %>%
  group_by(id, batch) %>%
  summarize(n = n(),
            mean = mean(agreement_rate),
            .groups = "drop") %>%
  ggplot(aes(n, mean)) +
  geom_point(size = 2, alpha = .7) +
  scale_x_log10() +
  geom_smooth(se = F, size = 2) +
  labs(title = "he number of outcomes in each study and between-coders reliability")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.30103

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.30103

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1.8832e-16

- it does not seem related

# study design and QEX between-coders reliability

data %>%
  select(estimate_id, evaluation_method) %>%
  filter(!str_detect(evaluation_method, regex("linked", ignore_case = T))) %>%
  distinct(estimate_id, .keep_all = T) %>%
  inner_join(each_estimate_id, by = "estimate_id") %>%
  mutate(evaluation_method = fct_reorder(evaluation_method, agreement_rate, .fun = mean)) %>%
  ggplot(aes(agreement_rate, evaluation_method, color = evaluation_method)) +
  geom_boxplot(size = 1, show.legend = F) +
  scale_x_continuous(label = percent_format()) +
  guides(color = guide_legend(reverse = TRUE)) +
  labs(y = "",
        x = "between-coders reliability",
       color = "") +
  theme(axis.text.y = element_text(hjust = .5)) +
  labs(title = "study design and QEX between-coders reliability")

batch level

batch_level <-
  study_level %>%
  group_by(batch) %>%
  summarize(mean = mean(mean),
            .groups = "drop")

batch_level

## # A tibble: 4 × 2
##   batch  mean
##   <chr> <dbl>
## 1 47    0.637
## 2 48    0.613
## 3 51    0.779
## 4 52    0.698

• These batch level data can be merged with the time spent on each batch and the number of studies and outcomes to quantify the difficulty of each studies

ROB

read_rob_sheet <- function(filename) {
  all_sheets <- excel_sheets(filename)
  
  # Find sheets that contain "QED" but not "group"
  qed_sheet <- all_sheets[
    str_detect(all_sheets, regex("QED", ignore_case = TRUE)) &
    !str_detect(all_sheets, regex("group", ignore_case = TRUE))
  ]
  
  # If matching sheet(s) found, read the first one
  if (length(qed_sheet) >= 1) {
    return(read_excel(filename, sheet = qed_sheet[1]))
  } else {
    return(NULL)
  }
}

folder_path_rob <- "d:/r/data extraction_ROB/"

data_rob <- 
  tibble(filename = list.files(path = folder_path_rob, full.names = TRUE, pattern = "\\.xlsx$")) %>%
  mutate(data = purrr::map(filename, read_rob_sheet)) %>%
  mutate(sheet = "QED") %>%
  unnest(data) %>%
  mutate(batch = str_extract(filename, "(?<=Group )\\d+")) %>% 
  group_by(batch) %>%
  arrange(Coder) %>%
  filter(!is.na(Coder)) %>%
  mutate(coders = paste0("coder_", dense_rank(Coder))) %>%
  ungroup() %>%
  janitor::clean_names() %>%
  select(-coder, -general_3, -general_5, -general_7, -general_6, -code, -x8_external_validity, -x9_limitations) %>%
  rename(estimate_id = "general_4",
         study_design = "general_8",
         eval_design = "general_9"
         ) %>%
  mutate(estimate_id = str_remove_all(estimate_id, "\r\n"))

## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## • `General` -> `General...3`
## • `General` -> `General...4`
## • `General` -> `General...5`
## • `General` -> `General...6`
## • `General` -> `General...7`
## • `General` -> `General...8`
## • `General` -> `General...9`
## • `1: Selection bias - Justification` -> `1: Selection bias -
##   Justification...11`
## • `1: Selection bias - Justification` -> `1: Selection bias -
##   Justification...12`
## • `1: Selection bias - Justification` -> `1: Selection bias -
##   Justification...13`
## • `2: Confounding - Justification` -> `2: Confounding - Justification...15`
## • `2: Confounding - Justification` -> `2: Confounding - Justification...16`
## • `2: Confounding - Justification` -> `2: Confounding - Justification...17`
## • `2: Confounding - Justification` -> `2: Confounding - Justification...18`
## • `2: Confounding - Justification` -> `2: Confounding - Justification...19`

data_rob_long <- data_rob %>%
  select(-contains("justification")) %>%
  pivot_longer(
    cols = -c(batch, coders, estimate_id, sheet), 
    names_to = "variable", 
    values_to = "value"
  ) %>%
  mutate(value = na_if(value, "NA"),
         value = na_if(value, "N/A")) %>%
  mutate(value = case_when(
    str_detect(variable, "x") & value %in% c("1", "2") ~ "yes",
    str_detect(variable, "x") & value %in% c("3", "4") ~ "no",
    str_detect(variable, "x") & value == "8" ~ "unclear",
    TRUE ~ value
  ))

extract_answers <- function(text) {
  answers <- list()
  for(letter in c('a', 'b', 'c', 'd')) {
    pattern <- paste0(letter, "\\)\\s*([\\w\\s]+)")
    match <- str_match(text, pattern)
    answer <- ifelse(is.na(match[1, 2]), NA, match[1, 2])
    if (!is.na(answer)) {
      answer <- str_remove_all(answer, "\r\n")
      answer <- str_replace_all(answer, regex("N/A", ignore_case = T), "na")
      answer <- str_replace_all(answer, regex("probably not",ignore_case = T), "no")  
      answer <- str_replace_all(answer, regex("probably no",ignore_case = T), "no")  
      answer <- str_replace_all(answer, regex("probably yes",ignore_case = T), "yes") 
      answer <- str_trim(answer)
      answer <- str_to_lower(answer)
      answer <- str_replace_all(answer, "na\\s+[a-d]", "na")  
      answer <- str_trim(answer)
      answer <- str_replace_all(answer,"^n$", "na")
      answer <- str_replace_all(answer, "unsure", "unclear")
    }
    answers[[letter]] <- answer
  }
  return(answers)
}


selection_sub <-
  data_rob %>%
  pivot_longer(
    cols = -c(batch, coders, estimate_id, sheet), 
    names_to = "variable", 
    values_to = "value"
  ) %>%
  mutate(value = na_if(value, "NA"),
         value = na_if(value, "N/A")) %>%
  mutate(value = case_when(
    str_detect(variable, "x") & value %in% c("1", "2") ~ "yes",
    str_detect(variable, "x") & value %in% c("3", "4") ~ "no",
    str_detect(variable, "x") & value == "8" ~ "unclear",
    TRUE ~ value
  )) %>%
  filter(str_detect(variable, "x1")) %>%
  filter(!str_detect(value, regex("linked", ignore_case = T))) %>%
  filter(!is.na(value)) %>%
  mutate(var_num = str_extract(variable, "\\d+$")) %>%
  mutate(variable = str_replace(variable, "justification_\\d+$", "justification")) %>%
  pivot_wider(
    id_cols = c(estimate_id, sheet, batch, coders),
    names_from = variable,
    values_from = c(value, var_num),
    names_glue = "{.value}_{variable}"
  ) %>%
  select(-var_num_x1_selection_bias_assessment) %>%
  add_count(estimate_id, var_num_x1_selection_bias_justification) %>%
  rename(justification = "value_x1_selection_bias_justification") %>%
  filter(n == 2) %>%
  select(-n) %>%
  arrange(estimate_id) %>%
  filter(var_num_x1_selection_bias_justification == 12) %>%
  mutate(answers = purrr::map(justification, extract_answers)) %>%
  unnest_wider(answers) %>%
  mutate_if(is.character, as.factor)

library(patchwork)
library(viridis)

## 필요한 패키지를 로딩중입니다: viridisLite

## Warning: 패키지 'viridisLite'는 R 버전 4.2.3에서 작성되었습니다

## 
## 다음의 패키지를 부착합니다: 'viridis'

## The following object is masked from 'package:scales':
## 
##     viridis_pal

library(ggrepel)

custom_colors <- viridis_pal()(4)

selection_a <-
  selection_sub %>% 
  add_count(value_x1_selection_bias_assessment, name = "total") %>%
  count(value_x1_selection_bias_assessment, a, total) %>%
  mutate(pct = n/ total) %>%
  mutate(a = fct_relevel(a, "yes", "no", "na", "unclear")) |> 
  ggplot(aes(pct, value_x1_selection_bias_assessment, fill = a)) +
  geom_col(width = .5) +
  scale_x_continuous(label = percent_format(),
                     expand = c(0,0)) +
  scale_fill_manual(values = custom_colors) +
  labs(y = "Sub-criteria a)",
       fill = "",
       x = "") +
  theme(legend.position = "top",
        axis.text.x = element_blank(),
        axis.ticks = element_blank(),
        axis.text.y = element_text(size = 10, hjust =.5))
  
  
selection_b <- 
  selection_sub %>% 
  add_count(value_x1_selection_bias_assessment, name = "total") %>%
  count(value_x1_selection_bias_assessment, b, total) %>%
  mutate(pct = n/ total) %>%
  mutate(b = fct_relevel(b, "yes", "no")) |> 
  ggplot(aes(pct, value_x1_selection_bias_assessment, fill = b)) +
  geom_col(width = .5, show.legend = F) +
  scale_x_continuous(label = percent_format(),
                     expand = c(0,0)) +
    scale_fill_manual(values = custom_colors) +
  labs(y = "Sub-criteria b)",
       fill = NULL,
       x = "")  +
  theme(legend.position = "top",
        axis.text.x = element_blank(),
        axis.ticks = element_blank(),
        axis.text.y = element_text(size = 10, hjust =.5))


selection_c <-
  selection_sub %>% 
  add_count(value_x1_selection_bias_assessment, name = "total") %>%
  mutate(either_yes = ifelse(a == "yes" | b == "yes", "yes", "no")) %>%
  count(value_x1_selection_bias_assessment, either_yes, total) %>%
  mutate(pct = n/ total) %>%
  mutate(either_yes = fct_relevel(either_yes, "yes", "no")) |> 
  ggplot(aes(pct, value_x1_selection_bias_assessment, fill = either_yes)) +
  geom_col(width = .5, show.legend = F) +
  scale_x_continuous(label = percent_format(),
                     expand = c(0,0)) +
    scale_fill_manual(values = custom_colors) +
  labs(y = "Either yes in a) or b)",
       fill = NULL,
       x = "")  +
  theme(legend.position = "top",
        axis.ticks = element_blank(),
        axis.text.y = element_text(size = 10, hjust =.5))

selection_a + selection_b + selection_c + plot_layout(nrow = 3)

Focusing on the selection bias assessment in PSM (N = 150 for two coders) - the current decision rules indicate either yes in a) or b) make “yes” for the overall analysis. The decision rules would be great to revisit the relevance between sub criteria and the final assessment outcome. - For instance, There is no sub-criteria for the self-selection bias. the sub criteria needs to revised to make another sub-criteria for the self-selection bias

confounding_sub <-
  data_rob %>%
  pivot_longer(
    cols = -c(batch, coders, estimate_id, sheet), 
    names_to = "variable", 
    values_to = "value"
  ) %>%
  mutate(value = na_if(value, "NA"),
         value = na_if(value, "N/A")) %>%
  mutate(value = case_when(
    str_detect(variable, "x") & value %in% c("1", "2") ~ "yes",
    str_detect(variable, "x") & value %in% c("3", "4") ~ "no",
    str_detect(variable, "x") & value == "8" ~ "unclear",
    TRUE ~ value
  )) %>%
  filter(str_detect(variable, "x2")) %>%
  filter(!str_detect(value, regex("linked", ignore_case = T))) %>%
  filter(!is.na(value)) %>%
  mutate(var_num = str_extract(variable, "\\d+$")) %>%
  mutate(variable = str_replace(variable, "justification_\\d+$", "justification"))%>%
  pivot_wider(
    id_cols = c(estimate_id, sheet, batch, coders),
    names_from = variable,
    values_from = c(value, var_num),
    names_glue = "{.value}_{variable}",
    values_fill = list(0),
    values_fn = list(value = list, var_num = list)) %>%
  unnest() %>%
  select(-var_num_x2_confounding_assessment) %>%
  add_count(estimate_id, var_num_x2_confounding_justification) %>%
  rename(justification = "value_x2_confounding_justification") %>%
  filter(n == 2) %>%
  select(-n) %>%
  arrange(estimate_id) %>%
  filter(var_num_x2_confounding_justification == 19) %>%
  mutate(answers = purrr::map(justification, extract_answers)) %>%
  unnest_wider(answers) %>%
  mutate_if(is.character, as.factor)

## Warning: `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(value_x2_confounding_assessment,
##   value_x2_confounding_justification, var_num_x2_confounding_assessment,
##   var_num_x2_confounding_justification)`.

data_rob_pivoted <- data_rob_long %>%
  add_count(batch,estimate_id, coders, name = "checker") %>%
  filter(checker != max(checker)) %>%
  select(-checker) %>%
  pivot_wider(
    names_from = coders,
    values_from = value
  ) %>%
    mutate(matched = case_when(
      is.na(coder_1) & is.na(coder_2) ~ NA,
      is.na(coder_1) & !is.na(coder_2) ~ FALSE,
      !is.na(coder_1) & is.na(coder_2) ~ FALSE,
      coder_1 == coder_2 ~ TRUE,
      TRUE ~ FALSE)) %>%
   filter(!is.na(matched))

data_rob_pivoted %>%
  count(batch, estimate_id, variable, matched) %>%
  mutate(matched = ifelse(matched, "count_true", "count_false")) %>%
  pivot_wider(
    names_from = matched,
    values_from = n,
    values_fill = 0 # Fill in 0 if there are missing values
  ) %>%
  count(variable, count_true) %>%
  filter(count_true == 1) %>%
  mutate(pct = n / max(n)) %>%
  filter(pct != 1) %>%
  mutate(variable = fct_reorder(variable, pct)) %>%
  ggplot(aes(pct, variable)) +
  geom_col() +
  scale_x_continuous(label = scales::percent_format()) +
  labs(y = "",
       x = "between-coder reliability",
       title = "ROB only QED")

estimate id and evaluation design

reference_table <- tibble(
  eval_design = c(1, 2, 3, 4, 5, 6, 7),
  Variable_Name = c(
    "Regression discontinuity",
    "Instrumental Variable",
    "Statistical Matching",
    "Difference-in-difference",
    "Interrupted time series",
    "Natural experiments",
    "Endogenous treatment-effect model"
  )
) %>%
  mutate(eval_design = as.character(eval_design))


data_rob_pivoted %>%
  inner_join(data_rob %>%
               distinct(estimate_id, .keep_all = T) %>%
               select(estimate_id, eval_design) %>% 
               filter(!str_detect(eval_design, regex("linked", ignore_case = T))), by = "estimate_id") %>%
  inner_join(reference_table, count_true, by = "eval_design") %>%
  add_count(variable, Variable_Name, name = "total") %>%
  count(batch, estimate_id, variable, Variable_Name, matched, total) %>%
  filter(variable != "filename") %>%
  mutate(matched = ifelse(matched, "count_true", "count_false")) %>%
  count(Variable_Name, variable, matched, total) %>%
  arrange(Variable_Name) %>%
  filter(matched == "count_true") %>%
  mutate(pct = n / total) %>%
  mutate(variable = fct_reorder(variable, pct)) %>%
  ggplot(aes(pct, Variable_Name, fill = Variable_Name)) +
  geom_col() +
  facet_wrap(~variable) + 
  scale_x_continuous(label = scales::percent_format(),
                     expand = c(0.01,0)) +
  labs(y = "",
       x = "between-coder reliability",
       title = "ROB only QED",
       fill = "") +
  theme(legend.position = "bottom",
        axis.text.y = element_text(hjust = .5),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 15, hjust = .5))

- selection bias (IV, DiD) - confounding (statistical matching, IV) - reporting bias (IV)

No sub-criteria is likely to make the between-coders reliability greater. For instance, IV seems the lowest, whereas the PSM seems the highest.

each estimate_id level

each_estimate_id_rob <-
  data_rob_pivoted %>%
  group_by(estimate_id, batch) %>%
  summarize(n = n(),
            agreement_rate = sum(matched == 1)/n,
            .groups = "drop")

each_estimate_id_rob

## # A tibble: 322 × 4
##    estimate_id batch     n agreement_rate
##    <chr>       <chr> <int>          <dbl>
##  1 73628889_1  48       10            0.6
##  2 73628889_2  48       10            0.6
##  3 73795470_1  46       10            0  
##  4 73795470_2  46       10            0  
##  5 73795470_3  46       10            0  
##  6 73795470_4  46       10            0  
##  7 73795470_5  46       10            0  
##  8 73795470_6  46       10            0  
##  9 73795904_1  47       10            0.4
## 10 73795904_10 47       10            0.4
## # ℹ 312 more rows

scatter plot between estimate_id and eval_design

each_estimate_id_rob %>%
  inner_join(data_rob %>%
               distinct(estimate_id, .keep_all = T) %>%
               select(estimate_id, eval_design) %>% 
               filter(!str_detect(eval_design, regex("linked", ignore_case = T))), by = "estimate_id") %>%
  inner_join(reference_table , by = "eval_design") %>%
  ggplot(aes(agreement_rate, fill = Variable_Name)) +
  geom_density(alpha = .7, linewidth = 0) +
  labs(fill = "",
       x = "Between-coders reliability",
       title = "Density plot depending on the evaluation method")

- IV seems the most difficult evaluation method to assess the risk of bias the mean is smaller than others. Endogenous treatment-effect model is the highest between-coders reliability

study level

study_level_rob <-
  each_estimate_id_rob %>%
  mutate(id = str_extract(estimate_id, "\\d+(?=_)")) %>%
  group_by(id, batch) %>%
  summarize(mean = mean(agreement_rate),
            .groups = "drop")

study_level_rob

## # A tibble: 49 × 3
##    id       batch  mean
##    <chr>    <chr> <dbl>
##  1 73628889 48    0.6  
##  2 73795470 46    0    
##  3 73795904 47    0.4  
##  4 73796232 47    0.4  
##  5 73797668 46    0.425
##  6 73797942 51    0.7  
##  7 73798036 46    0.5  
##  8 73798224 48    0    
##  9 73798237 46    0    
## 10 73798251 46    0    
## # ℹ 39 more rows

• These study level data can be merged with the study feature such as evaluation design, the number of outcomes and others

batch level

batch_level_rob <-
  study_level_rob %>%
  group_by(batch) %>%
  summarize(mean = mean(mean),
            .groups = "drop")

batch_level_rob

## # A tibble: 5 × 2
##   batch  mean
##   <chr> <dbl>
## 1 46    0.222
## 2 47    0.617
## 3 48    0.627
## 4 51    0.617
## 5 52    0.670

• These batch level data can be merged with the time spent on each batch and the number of studies and outcomes to quantify the difficulty of each studies

QEX + ROB

batch_level_rob |> 
  rename(rob = "mean") |> 
  inner_join(batch_level |> 
               rename(qex = "mean"), 
             by = "batch") %>%
  mutate(mean = (rob + qex) / 2) 

## # A tibble: 4 × 4
##   batch   rob   qex  mean
##   <chr> <dbl> <dbl> <dbl>
## 1 47    0.617 0.637 0.627
## 2 48    0.627 0.613 0.620
## 3 51    0.617 0.779 0.698
## 4 52    0.670 0.698 0.684