library(tidyverse) # Data wrangling and plotting
library(janitor) # To change quickly change names of columns
library(metafor) # To create priors from meta-analysis
library(ggdist) # To plot distributions
library(patchwork) # To arrange plots
library(kableExtra) # To create nice tables
BaSICS was kind of a big deal earlier this month. You may have read about it on Twitter when all the excitement was being live-streamed…
The choice of which fluid to use and when is an argument that intensivists love, and combining multiple studies seemed to suggest that using balanced crystalloids (fluids which more closely mimic the composition of blood plasma) have better outcomes in critically unwell adults than the “old faithful” of 0.9% sodium chloride (“abnormal” saline), which has a high chloride load (and hence is acidotic) and lacks any potassium, calcium, or other electrolytes.
Then along came BaSICS. Aiming to randomise 11,000 critically ill patients to medical brine vs a balanced solution (in this case, Plasma-Lyte 148, beloved of neurointensivists), it reported recently via alifestream on the amazing Critical Care Reviews website that there was no difference between mortality rates, or rates of acute kidney injury requiring renal replacement therapy, between the two (p = 0.47). Cue the “balanced crystalloids are pointless” editorials…
Hmmm. A yes/no answer supported by an arbitrary p-value? Sounds like an excuse for yet more Bayesian statistics chat.
A couple of months ago Arthur Albuquerque, a medical student and Bayesian stats wizard, published a fantastic Bayesian renalaysis of the RECOVERY trial’s toilizumab results, so we’ll borrow his code and do the same thing for BaSICS. For simplicity, we’ll simply look at mortality in all patients, but AKI and/or subgroup analysis would be performed the same way.
Get some data
We’ll start with importing packages we need:
We’ll borrow the data from Hammond et al’s excellent pre-BaSICS meta-analysis, focusing just on the randomised controlled trials and just the mortality outcomes, and add the outcomes from BaSICS to it.
d = read_csv("extracted-data.csv")
d = clean_names(d) | trial | events_balanced | total_balanced | events_control | total_control |
|---|---|---|---|---|
| Annane 2013 | 22 | 72 | 275 | 1035 |
| Semler 2017 | 72 | 520 | 68 | 454 |
| Semler 2018 | 928 | 7942 | 975 | 7860 |
| Verma 2016 | 5 | 33 | 2 | 34 |
| Young 2014 | 3 | 32 | 4 | 33 |
| Young 2015 | 87 | 1152 | 95 | 1110 |
| Zampieri 2021 | 1381 | 5230 | 1439 | 5290 |
Calculate distributions
First step in Bayesian-ing this up is working out the posterior distributions (what do BaSICS tell us, and all the other trials that Hammond analyses, tell us in isolation) and the prior distributions (what did we think before BaSICS).
Lets start with working out the odds ratio for both BaSICS, and everything else combined:
# Log(odds ratio) for BASICS
logOR_basics_mortality =
rma(
# Outcome
measure = "OR",
# Balanced crystalloid group
ai = events_balanced,
n1i = total_balanced,
# Control group
ci = events_control,
n2i = total_control,
# Mechanics
data = d %>% filter(outcome == "mortality",
trial == "Zampieri 2021"),
method = "REML",
slab = paste(trial)
)
# Log(odds ratio) for all other trials
logOR_prior_mortality =
rma(
# Outcome
measure = "OR",
# Balanced crystalloid group
ai = events_balanced,
n1i = total_balanced,
# Control group
ci = events_control,
n2i = total_control,
# Mechanics
data = d %>% filter(outcome == "mortality",
trial != "Zampieri 2021"),
method = "REML",
slab = paste(trial)
)So now we build a range of prior distributions - what we knew (or thought) before BaSICS was on the scene. As in the Albuquerque paper, we’ll use a number of priors to determine the strength of the outcomes after BaSICS. These include:
- The “evidence based prior” (i.e. what we think after Hammond’s summary of the evidence)
- An “optimistic” prior (what the authors were hoping for, which was a 10% reduction in the odds ratio1)
- A “pessimistic” prior (so the opposite of the optimistic, that balanced crystalloids cause harm)
- A skeptical (no difference, but same variance as the evidence-based prior) and non-informative (no information at all) prior.
To do this we’ll re-use Albuquerque’s scripts, available from Github, modified for our ORs:
# Normal conjugate posterior probabilities
source("conjugate_normal_posterior.R")
# Generate posterior probabilities with an evidence-based prior
source("normal_approximation_logOR.R")
# Generate posterior probabilities with multiple priors
source("logOR_normal_approximation_multiple_priors.R")
# Compile all posteriors into a tibble
source("tibble_all_posteriors_logOR.R") Now we can build our prior distributions, and calculate the posterior distributions from applying the results of BaSICS to each of these:
# Build evidence-based priors
posterior_mortality = normal_approximation_logOR(
logOR_basics_mortality, # Data object with BASICS data
logOR_prior_mortality, # Output from the rma() in the previous code chunk
posterior_mortality
)
# Build other priors
posteriors_mortality = logOR_normal_approximation_multiple_priors(
posterior_mortality, # Data object with evidence-based posterior from normal_approximation()
"mortality", # Outcome (within quotes)
posteriors_mortality # Name for the final output with data+prior+posterior
)
# Calculate OR for posterior distributions based on the priors
samples_logOR_posteriors_mortality =
tibble_all_posteriors_logOR(posteriors_mortality)
samples_posteriors_mortality =
samples_logOR_posteriors_mortality %>%
pivot_longer(
cols = c('Non-informative':'Pessimistic'),
names_to = "Underlying_Prior", values_to = "log(OR)"
) %>%
mutate(OR = exp(`log(OR)`))| type | data.mean | data.sd | prior.mean | prior.sd | post.mean | post.sd |
|---|---|---|---|---|---|---|
| non-informative | -0.04062107 | 0.0440286 | 0.0000000 | 100.0000000 | -0.0406211 | 0.0440286 |
| evidence-based | -0.04062107 | 0.0440286 | -0.0653904 | 0.0444135 | -0.0528979 | 0.0312681 |
| skeptical | -0.04062107 | 0.0440286 | 0.0000000 | 0.0444135 | -0.0204873 | 0.0312681 |
| optimistic | -0.04062107 | 0.0440286 | -0.1053605 | 0.0444135 | -0.0727091 | 0.0312681 |
| pessimistic | -0.04062107 | 0.0440286 | 0.1053605 | 0.0444135 | 0.0317344 | 0.0312681 |
Visualise some results
So what does all this mean? Remember that a Bayesian analysis is taking your prior beliefs (in this case, either the evidence based belief from the Hammond paper, or the pre-specified other priors), and applying the likelihood from the new results (BaSICS) to update your posterior beliefs. Let’s visualise this.
Firstly, we’ll extract the priors…
# Extract priors
priors_mortality =
d %>%
filter(
trial != "Zampieri 2021",
outcome == "mortality"
)
basics_mortality =
d %>%
filter(
trial == "Zampieri 2021",
outcome == "mortality"
)
# Set seed for reproducibility
set.seed = 69 #dudes… and combine these into a dataframe of distributions…
### Create distributions dataframe ----
distributions =
posteriors_mortality %>%
filter(type == "evidence-based") %>%
summarise(
BASICS = list(rnorm(10e4,
mean = data.mean,
sd = data.sd
)),
Prior = list(rnorm(10e4,
mean = prior.mean,
sd = prior.sd
))
) %>%
unnest(BASICS:Prior)
distributions = bind_cols(
distributions,
# Bind corresponding posterior distribution data
samples_posteriors_mortality %>%
filter(`Underlying_Prior` == "Evidence-based") %>%
select(`log(OR)`) %>%
rename("Posterior" = `log(OR)`)
)… which we can then plot:
# Prepare arrows (After https://github.com/rdboyes/forester/blob/master/R/forester.R)
arrow_labels = c("Balanced solution beneficial", "Balanced solution harmful")
xlim = c(0, 20)
xlab_df = data.frame(text = arrow_labels,
x = c(2,18),
y = c(0, 0),
hjust = c(0, 1))
a_small_amount = abs(xlim[1] - xlim[2])/35
null_line_at = 10
arrow_df = data.frame(id = c(1,2),
xstart = c(null_line_at - a_small_amount,
null_line_at + a_small_amount),
xend = c(xlim[1] + a_small_amount, xlim[2] - a_small_amount),
y = c(1, 1))
arrows_plot = ggplot() +
geom_segment(data = arrow_df,
aes(x = .data$xstart,
xend = .data$xend,
y = .data$y,
yend = .data$y),
arrow = arrow(angle = 15, type = "closed", length = grid::unit(0.1, "in"))) +
geom_text(data = xlab_df,
aes(x = .data$x,
y = .data$y,
label = .data$text,
hjust = .data$hjust), size = 3.5) +
scale_y_continuous(expand = c(0,0), limits = c(-0.5, 1.75)) +
scale_x_continuous(expand = c(0,0), limits = xlim) +
theme(panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent", color = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.background = element_rect(fill = "transparent"),
legend.box.background = element_rect(fill = "transparent"),
panel.border = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank())
# Construct plot
plot1 <- distributions %>%
# make it tidy for ggplot
pivot_longer(
cols = c(BASICS:Posterior),
names_to = "dist", values_to = "samples"
) %>%
mutate(
# Set order of distributions
dist = factor(dist,
levels = c(
"Posterior",
"BASICS",
"Prior"
)
)
) %>%
# exp() to convert log-odds ratio into odds ratio
ggplot(aes(
x = exp(samples), y = dist,
fill = dist
)) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(position = "dodge",
point_interval = median_hdi, # Highest density interval
.width = c(0.95)) +
geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
scale_fill_manual(' ',
breaks=c("Prior","BASICS", "Posterior"),
values = c(
"#E0C6B6", # Prior
"#9FA464", # BASICS
"#9B5446" # Posterior
)) +
labs(
x = "\nOdds Ratio",
y = " "
) +
scale_x_continuous(breaks = seq(from = 0.25, to = 1.75, 0.25)) +
coord_cartesian(xlim = c(0.25, 1.75),
clip = 'off') +
scale_y_discrete(expand = c(0, 0.5)) +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 11),
axis.text.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = 'right',
#legend.text.align = 0.5,
legend.margin = margin(0, 0, 0, 10),
legend.text = element_text(size=12),
legend.key= element_blank(),
plot.margin = margin(20, 0, 60, 25)
) +
inset_element(arrows_plot,
ignore_tag = TRUE,
align_to = "full",
left = unit(0, 'cm'),
bottom = unit(0.4, 'cm'),
right = unit(23, 'cm'),
top = unit(2.3, 'cm'))plot1
What does it all mean?
So how do we interpret this?
The probability distribution (probability of the real value, given the data available) is shown as a bell curve. The 95% credible interval (kind of like the 95% confidence interval, but in this case actually giving you the probability of the real value, rather than the range of values given the null hypothesis) is shown as the dark line.
The prior, pink, shows the range of potential odds ratios given the pre-existing evidence base. The likelihood of the true odds ratio from the BaSICS data is in green, and the posterior distribution (the prior, updated by the new data) is in red, seeming to show that BaSICS in fact confirmed our pre-existing suspicions that balanced crystalloids are better.
But what is we think the evidence base on which BaSICS was designed is flawed? Well, we can compare the effect of different priors:
or_mortality = samples_posteriors_mortality %>%
select(`Underlying_Prior`, OR)
plot2 <- or_mortality %>%
# Define orders
mutate(`Underlying_Prior` =
factor(`Underlying_Prior`,
levels = c("Pessimistic",
"Optimistic",
"Skeptical",
"Non-informative",
"Evidence-based"
)
)
) %>%
# Plot!
ggplot(aes(
x = OR, y = `Underlying_Prior`,
# https://github.com/mjskay/ggdist/issues/71
fill_ramp = stat(x < 1))
) +
geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(position = "dodge",
point_interval = median_hdi, # Highest density interval
.width = 0.95,
# https://github.com/mjskay/ggdist/issues/70
interval_size_range = c(0.4,1.2)
) +
scale_fill_ramp_discrete(from = "gray85", range = c(0,1)) +
scale_fill_manual(values = c("#516F88", "#8DBEE9")) + # Using
labs(x = "\nOdds Ratio",
y = "Underlying Prior\n \n"
) +
theme(
strip.background = element_rect(fill = "#E4E6E7"),
strip.text.x = element_text(size = 14),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 16),
axis.text.y = element_text(size = 15),
axis.text.x = element_text(size = 12),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = "none",
plot.title.position = 'plot',
plot.margin = margin(20, 25, 60, 20)
) +
scale_x_continuous(breaks = seq(from = 0.25, to = 1.75, 0.25)) +
coord_cartesian(xlim = c(0.25, 1.75)) +
scale_y_discrete(expand = c(0, 0.3)) +
inset_element(arrows_plot,
ignore_tag = TRUE,
align_to = "full",
left = unit(5, 'cm'),
bottom = unit(0, 'cm'),
right = unit(26, 'cm'),
top = unit(2.3, 'cm'))plot2
Holding the results of BaSICS in isolation (non-informative and skeptical priors) seems to show no conclusive benefit to using balanced crystalloids, although the probability does lean towards balanced solutions being beneficial - more data needed. However, using the evidence-based prior (based on what we knew before BaSICS), BaSICS actually narrows down the credible interval, providing more evidence that balanced crystalloids confer a survival benefit, and the same finding is seen given an optimistic prior (that the BaSCIS investigators powered their study around). So based on these, BaSICS again supports the use of balanced crystalloids.
In light of this then, what BaSICS has done is helped narrow down what we already knew - balanced crystalloids are (probably) better for the critically unwell. This then feeds into the other criticisms of the paper (high number of elective - and hence “well critically unwell” patients, low numbers of septic patients, etc) to suggest that actually that although a great paper, especially given the number of patients collected in the middle of a COVID-19 pandemic, this isn’t the fluids debate paper to end all fluids debates paper that Twitter seems to think it is.
One last thing…
More worryingly, BaSICS also seemed to find a signal for harm using balanced solutions in patients with traumatic brain injuries, so my next step will be looking into this in a bit more detail. This will require a bit more digging into the literature, and the oversight of a proper statistician, so get in touch if you want to help out with this…
Footnotes
Technically the hazard ratio, but it’s close to 1 so we can assume HR ~ OR.↩︎