This R Markdown-notebook contains the final annotated analysis code for the manuscript: “Overall bias and sample sizes were unchanged in ICU trials over time: a meta-epidemiological study” by Carl Thomas Anthon, Anders Granholm, Jon Henrik Laake, Anders Perner and Morten Hylander Møller.

The annotated analysis code was written by Anders Granholm according to the pre-published protocol and statistical analysis plain (doi: 10.1111/aas.13258). For further details, please see the full article and its supplementary files, which also includes the CSV-file containing the final dataset and a data dictionary.

For correspondence related to the manuscript, contact the corresponding author Carl Thomas Anthon. Questions directly related to the analysis code can also be sent to Anders Granholm.

Figures displayed in this document are scaled differently than in the final paper and supplement.

1. Setup

Load the packages/functions required for the analyses.

library(tidyverse) # Load the Tidyverse packages, used for data processing and visualisation.
plot_grid <- cowplot::plot_grid # Load this function from the cowplot package; only this function will be used from this package, and the entire library will not be loaded as this is unnecessary and cowplot overrides the default ggplot2 theme, which then has to be reset.

Import the file containing data for all the included trials.

Note: This code block uses RStudio’s API to allow manual selection of the final dataset. This requires the script to be run in a recent version of RStudio; otherwise this has to be manually replaced with the appropriate file path.

included_trials <- read_csv(rstudioapi::selectFile(filter = "CSV files (*.csv)"), col_types = cols(trial_year = col_integer(), sample_size = col_integer()))

Define the output directory where all results will be saved.

Note: This block uses RStudio’s API to allow manual selection of the directory. This requires the script to be run in a recent version of RStudio; otherwise this has to be manually replaced with the appropriate path.

output_path <- rstudioapi::selectDirectory()

Define plot theme and image settings for saving images, to allow easy changing later.

# File settings
img_dpi <- 600 # DPI - not used for vectorized formats, e.g. PDFs, but argument is kept if this is changed.
img_format <- "pdf" # File format
# Widths and heights (in millimeters) for single run chart plots
img_width <- 150
img_height <- 80
# Widths and heights (in millimeters) for the triple run chart plots
img_width_triple <- 150
img_height_triple <- 150
# Width and height (in millimeters) for the risk of bias overview plot
img_width_bias <- 200
img_height_bias <- 80

# Set ggplot2 theme
theme_set(theme_linedraw())

2. Data inspection and preparation

2.1 Helper function

The following helper function filters a dataset so only years with a certain number of trials are included, and optionally prints the excluded years. Arguments:

  • .dta: the dataset to filter.
  • minimum: the minimum number of required trials per year, default is 5.
  • verbose: if TRUE (default), the years with too few trials are printed.
filter_years <- function(.dta, minimum = 5, verbose = TRUE){
  years <- .dta %>% 
    group_by(trial_year) %>% 
    summarise(n = n()) %>% 
    filter(n < minimum) %>% 
    pull(trial_year)
  
  if (verbose){
    print(str_glue("The following years contain less than {minimum} trials and have beeen filtered from this dataset: {str_c(years, collapse = ', ')}."))
  }
  
  .dta %>% 
    filter(!(trial_year %in% years))
}

2.2 Risk of bias

Create overall_risk_of_bias variable and convert all risk of bias domains to factors.

included_trials <- included_trials %>% 
  mutate(overall_risk_of_bias = case_when(is.na(random_sequence_generation) | is.na(allocation_concealment) | is.na(blinding_of_participants_and_personnel) | is.na(blinding_of_outcome_assessment) | is.na(incomplete_outcome_data) | is.na(selective_reporting) | is.na(other_sources_of_bias) ~ NA_character_, # If any risk of bias domain is missing, the overall risk of bias will be missing
                                          random_sequence_generation == "Low" & allocation_concealment == "Low" & blinding_of_participants_and_personnel == "Low" & blinding_of_outcome_assessment == "Low" & incomplete_outcome_data == "Low" & selective_reporting == "Low" & other_sources_of_bias == "Low" ~ "Low", # If all domains are low risk, the overall risk of bias will be low as well
                                          TRUE ~ "High" # If none of the above criteria apply, the overall risk of bias will be high
         )) %>% 
  mutate_at(vars(random_sequence_generation:overall_risk_of_bias), ~factor(., levels = c("Low", "Unclear", "High")))

Assess how many trials have missing data for the overall risk of bias domain (which will be missing if any of the other domains are missing) and create dataset that only includes trials with available risk of bias ratings.

included_trials %>% 
  filter(is.na(overall_risk_of_bias)) %>% 
  str_glue_data("In total, {nrow(.)} trials have missing data for one or more risk of bias domains and will be excluded from these analyses.")
In total, 2 trials have missing data for one or more risk of bias domains and will be excluded from these analyses.
trials_bias_available <- included_trials %>% 
  filter(!is.na(overall_risk_of_bias))

Create dataset that only contains trials from years with at least 5 trials with available risk of bias ratings (for the primary analyses).

trials_bias_available_5plus <- filter_years(trials_bias_available)
The following years contain less than 5 trials and have beeen filtered from this dataset: 1977, 1978, 1981, 1982, 1983, 1984, 1986, 2017, 2018.

For all trials with available risk of bias, calculate overall distributions for each domain.

trials_bias_available %>% 
  gather(key = "bias_domain", value = "rating", random_sequence_generation:overall_risk_of_bias, factor_key = TRUE) %>% 
  select(bias_domain, rating) %>% 
  group_by(bias_domain) %>% 
  summarise(low = sum(rating == "Low") / n() * 100,
            unclear = sum(rating == "Unclear") / n() * 100,
            high = sum(rating == "High") / n() * 100) %>% 
  mutate_if(is.numeric, ~paste0(round(., 1), "%"))

Create plot visualising distributions of risk of bias for each domain (only include trials where risk of bias ratings are available).

trials_bias_available %>% 
  gather(key = "bias_domain", value = "rating", random_sequence_generation:overall_risk_of_bias, factor_key = TRUE) %>% 
  select(bias_domain, rating) %>%
  mutate(rating = factor(rating, levels = c("High", "Unclear", "Low")),
         bias_domain = fct_rev(bias_domain)) %>%  # Reverse the order to ensure proper printing in the table
  ggplot(aes(x = bias_domain, fill = rating)) + # Create the plot
  geom_bar(position = "fill", width = 0.95) +
  coord_flip() +
  labs(x = NULL, y = NULL) +
  scale_fill_manual(values = c("#C00000", "#E0E000", "#00C000"), name = "Risk of bias rating:", guide = guide_legend(reverse = TRUE)) +
  scale_y_continuous(expand = c(0, 0), breaks = 0.1 * 0:10, minor_breaks = 0.05 * 0:20, labels = paste0(10 * 0:10, "%")) +
  scale_x_discrete(expand = c(0, 0), labels = rev(c( # Use rev() to ensure that order matches the domains above
    "Random sequence generation",
    "Allocation concealment",
    "Blinding of participants and personnel",
    "Blinding of outcome assessment",
    "Incomplete outcome data",
    "Selective reporting",
    "Other sources of bias",
    "Overall risk of bias"
  ))) +
  theme(plot.margin = margin (t = 5, l = 10, r = 15, b = 5, unit = "pt"), legend.position = "bottom") # Fix plot margins to ensure all labels are printed and change legend

ggsave(paste0(output_path, "/Risk of bias overview.", img_format), width = img_width_bias, height = img_height_bias, dpi = img_dpi, units = "mm") # Save file

2.3 Sample sizes

Assess how many trials have missing data for sample size and create dataset that only includes trials with available sample sizes.

included_trials %>% 
  filter(is.na(sample_size)) %>% 
  str_glue_data("In total, {nrow(.)} trials have missing data for sample sizes and will be excluded from these analyses.")
In total, 2 trials have missing data for sample sizes and will be excluded from these analyses.
trials_size_available <- included_trials %>% 
  filter(!is.na(sample_size))

Create dataset that only contains trials from years with at least 5 trials with available sample sizes (for the primary analyses).

trials_size_available_5plus <- filter_years(trials_size_available)
The following years contain less than 5 trials and have beeen filtered from this dataset: 1977, 1978, 1981, 1982, 1983, 1984, 1986, 2017, 2018.

For all trials with available sample sizes, calculate the overall distribution.

round(summary(trials_size_available$sample_size), 1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    3.0    42.0    82.5   224.0   181.2  7000.0 

3. Define run chart functions

Define the function that will produce the run charts and associated statistics.

This contains two helper functions, run_chart_stat and run_chart_caption which should only be called from the main function, run_chart, and not directly.

The run_chart function has the following arguments:

The function returns an annotated plot.

An additional, simplified helper function is present, triple_bias_run_chart, with the following arguments:

# Function that calculates the run chart statistics using the "val" column in the aggregated dataset from the main function that has already been sorted by year.
# This function should not be called directly
run_chart_stat <- function(val, fixed_median = NULL){
  n_total <- length(val)
  if (n_total < 1){
    stop("ERROR: run_chart_stat: number of observations (years) < 1.")
  }
  if (!(is.null(fixed_median) | is.numeric(fixed_median))) {
    stop("ERROR: run_chart_stat: fixed_median must be either NULL or a numeric value.")
  }
  med <- ifelse(is.null(fixed_median), median(val), fixed_median)
  n_useful <- sum(!near(val, med))
  
  run_longest <- 0
  run_current <- 0
  pos <- NA
  cross_n <- NA
  
  for (i in 1:n_total) {
    cur <- val[i]
    if (!near(cur, med)) {
      
      if (is.na(pos)){ # Will be true for the first value that is not equal to the median (and only once)
        pos <- !(cur > med) # Sets current position as the opposite of what it is, to start new run
        cross_n <- -1 # The first run will not count as a crossing
      }
      
      if ((cur > med ) == pos) { # If the current value is on the same side as the last one
        run_current <- run_current + 1
      } else { # If the current value is on the opposite side as the last one
        if (run_current > run_longest) {
          run_longest <- run_current
        }
        
        run_current <- 1
        cross_n <- cross_n +1
        pos <- cur > med
      }
      
    }
  }
  if (run_current > run_longest) {
    run_longest <- run_current
  }
  
  # Return results (upper and lower prediction limits calculated according to Anhøj's rules)
  list(n_total = n_total, n_useful = n_useful, median = med,
       run_UPL = round(log2(n_useful) + 3), run_longest = run_longest,
       cross_LPL = qbinom(0.05, n_useful - 1, 0.5), cross_n = cross_n, fixed_median = !is.null(fixed_median))
}

# Function that formats the caption used in each run chart
# This function should not be called directly
run_chart_caption <- function(res, percentage = NULL, dec = 1){
  res$base <- paste0(res$n_total, " observations (", res$n_useful, " useful). ", ifelse(res$fixed_median, "Fixed median: ", "Median: "), format(res$median, digits = dec, nsmall = dec), percentage, ". Longest run:")
  res$long <- paste0("[max. ", res$run_UPL, "].")
  res$cross <- paste0("[min. ", res$cross_LPL, "].")
  
  if (res$run_longest > res$run_UPL & res$cross_n < res$cross_LPL){
    # Both results are significant (non-random variation)
    res <- map(res, as.character)
    cap <- substitute(base~bold(run_longest)~long~"Crossings:"~bold(cross_n)~cross, res)
  } else if (res$run_longest > res$run_UPL) {
    # Only the longest run is significant (non-random variation)
    res <- map(res, as.character)
    cap <- substitute(base~bold(run_longest)~long~"Crossings:"~cross_n~cross, res)
  } else if (res$cross_n < res$cross_LPL) {
    # Only the number of crossings is significant (non-random variation)
    res <- map(res, as.character)
    cap <- substitute(base~run_longest~long~"Crossings:"~bold(cross_n)~cross, res)
  } else {
    # No significance (only random variation)
    res <- map(res, as.character)
    cap <- substitute(base~run_longest~long~"Crossings:"~cross_n~cross, res)
  }
  # Return caption
  cap
}

# The main function - see description above
run_chart <- function(.dta, variable, category = NULL, title = NULL, delete_titles = TRUE, fixed_median = NULL, num_agg_mean = FALSE){
  # Quote variable name and get the type of it
  var <- enquo(variable)
  var_pull <- pull(.dta, !!var)
  type <- case_when(is.factor(var_pull) ~ "factor",
                    is.integer(var_pull) ~ "integer",
                    TRUE ~ "other")
  
  # Input checks
  if (type == "factor") {
    if (is.null(category)) {
      stop("ERROR: run_chart: 'variable' is a factor (risk of bias domain?) and 'category' is NULL. 'category' must be either 'Low', 'Unclear' or 'High'.", call. = FALSE)
    } else if (!(category %in% c("Low", "Unclear", "High"))){
      stop(paste0("ERROR: run_chart: 'variable' is a factor (risk of bias domain?), 'category' is '", category, "', but must be either 'Low', 'Unclear' or 'High'."), call. = FALSE)
    }
  } else if (type == "integer") {
    if (!is.null(category)){
      warning(paste0("WARNING: run_chart: 'variable' is an integer (sample size?). 'category' is '", category, "' and will NOT be used."), call. = FALSE)
    }
  } else {
    stop(paste0("ERROR: run_chart: 'variable' must be either a factor (risk of bias category) or an integer (sample size), '", quo_name(var), "' is neither."))
  }
  
  # Remove missing data
  n_miss <- .dta %>% 
    filter(is.na(trial_year) | is.na(!!var)) %>% 
    nrow()
  if (n_miss > 0) {
    warning(paste0("WARNING: run_chart: ", n_miss, " rows with missing data in either trial_year or ", quo_name(var), " removed from the dataset prior to analyses."))
    .dta <- .dta %>% 
      filter(!is.na(trial_year) & !is.na(!!var))
  }
  
  # If delete_titles == TRUE, then set title to null
  if (delete_titles){
    title <- NULL
  }
  
  # Aggregate data
  if (type == "factor") {
    agg <- .dta %>% 
      group_by(trial_year) %>% 
      summarise(val = sum(!!var == category) / n() * 100) %>% 
      arrange(trial_year)
  } else if (type == "integer") {
    if (!num_agg_mean){
      agg <- .dta %>% 
        group_by(trial_year) %>% 
        summarise(val = median(!!var),
                  y_high = quantile(!!var, 0.75),
                  y_low = quantile(!!var, 0.25)) %>%
        arrange(trial_year)      
    } else if (num_agg_mean) {
      agg <- .dta %>% 
        group_by(trial_year) %>% 
        summarise(val = mean(!!var),
                  se_val = sd(!!var) / sqrt(n()),
                  y_high = val + 1.96 * se_val,
                  y_low_tmp = val - 1.96 * se_val,
                  y_low = ifelse(y_low_tmp < 0, 0, y_low_tmp)) %>%  # 95% confidence intervals (CIs) below 0 are impossible
        arrange(trial_year)
    }

  }
  
  # Calculate stats
  stats <- agg %>%
    pull(val) %>% 
    run_chart_stat(fixed_median = fixed_median)
  
  # General plot aspects
  gg <- ggplot(data = agg, aes(x = trial_year, y = val)) +
    labs(caption = run_chart_caption(stats, if_else(type == "factor", "%", "")), x = NULL, y = NULL, subtitle = title) + # The title argument provided is used for a ggplot2 subtitle, due to the default smaller font size
    theme(plot.caption = element_text(hjust = 0.5)) + # Align caption in the middle
    scale_x_continuous(limits = c(1975, 2020), breaks = 1975 + 5 * 0:8, minor_breaks = 1975:2020, expand = c(0, 0))
  
  # Add different aspects for risk of bias domains and sample sizes
  if (type == "factor") {
    # Risk of bias domains - different y-axis for overall risk of bias plots and subdomains
    if (quo_name(var) == "overall_risk_of_bias") {
      gg <- gg +
        scale_y_continuous(limits = c(0, 50), breaks = 10 * 0:5, minor_breaks = 5 * 0:10, labels = paste0(10 * 0:5, "%"), expand = c(0.015, 0))
    } else {
      gg <- gg +
        scale_y_continuous(limits = c(0, 100), breaks = 20 * 0:5, minor_breaks = 5 * 0:20, labels = paste0(20 * 0:5, "%"), expand = c(0.015, 0))
    }
  } else if (type == "integer") {
    # Sample sizes - different y-axis according to aggregation function
    gg <- gg +
      scale_y_continuous(breaks = 100 * 0:15, minor_breaks = 25 * 0:60, expand = c(0, 0)) + # Axis long enough for both types of plots - different limits for median/mean values are set in next line
      coord_cartesian(ylim = c(0, ifelse(!num_agg_mean, 600, 1400))) + # Limits are set here to avoid outliers for 95% CIs of yearly means to be removed
      geom_ribbon(aes(ymin = y_low, ymax = y_high), fill = "#B5B5B5", alpha = 0.80) # Add grey shaded area
  }
  
  # Add points and lines
  gg <- gg + geom_point(size = 2) +
    geom_line() +
    geom_hline(yintercept = stats$median, size = 1.25)
  
  # Return plot
  gg
}

# Simple helper function that creates a combined plot with all 3 risk of bias ratings
triple_bias_run_chart <- function(.dta, variable, title_add = NULL){
  variable <- enquo(variable)
  title_base <- str_replace_all(str_to_title(quo_name(variable)), "_", " ")
  plot_grid(
    run_chart(.dta, !!variable, category = "Low", title = paste0(title_base, title_add)) + ylab("Low"),
    run_chart(.dta, !!variable, category = "Unclear") + ylab("Unclear"),
    run_chart(.dta, !!variable, category = "High") + ylab("High"),
    nrow = 3, ncol = 1)
}

4. Primary analysis

For all primary analyses, only years with 5 or more trials are used.

Analysis 1.1: Percentage of overall low risk of bias trials per year

trials_bias_available_5plus %>% 
  run_chart(overall_risk_of_bias, "Low", title = "Overall low risk of bias")

ggsave(paste0(output_path, "/1.1 Primary analysis - overall low risk of bias.", img_format), width = img_width, height = img_height, dpi = img_dpi, units = "mm")

Analysis 1.2: Annual median sample sizes

trials_size_available_5plus %>% 
  run_chart(sample_size, title = "Median (IQR) sample sizes")

ggsave(paste0(output_path, "/1.2 Primary analysis - sample sizes.", img_format), width = img_width, height = img_height, dpi = img_dpi, units = "mm")