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:

**.dta**: the dataset used.
**variable**: the variable (e.g. a risk of bias domain or the sample sizes variable), which will be assessed.
**category**: if the variable is a risk of bias domain, the rating (e.g. “Low”) that will be assessed.
**title**: the plot title, if any is wanted. See argument below.
**delete_title**: this argument was added after the majority of the code was written. If TRUE (current default), no plot titles will be used regardless of the above argument.
**fixed_median**: defaults to NULL, and then the median value will be calculated using the data provided. If a numeric value is provided, the median will be fixed to this value. This argument is currently not used in any analyses.
**num_agg_mean**: defaults to FALSE. If FALSE, for numeric values (i.e. yearly samply sizes) the median value will be used in the run chart and the interquartile range will be displayed. If TRUE, mean values will be used and the 95% confidence intervals for the mean will be displayed instead. Ignored for risk of bias domains.

The function returns an annotated plot.

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

**.dta**: the dataset used.
**variable**: the risk of bias domain that will be assessed.
**title_add**: extra text added as part of the title (added directly after the bias domain), if NULL (default), nothing will be added.

```
# 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")
```