ordinal PARAMEDICS2

ordinal_PARAMEDICS2.md

This code was inspired by this post: https://discourse.datamethods.org/t/a-randomized-trial-of-epinephrine-in-out-of-hospital-cardiac-arrest

Install/Load packages

pacman::p_load(brms,
               tidybayes,
               ggplot2,
               dplyr,
               tidyr,
               rms,
               parameters,
               gt)

Load data from https://www.nejm.org/doi/full/10.1056/NEJMoa1806842

a <- c(rep(0,15), rep(1,10), rep(2,29), rep(3,20), rep(4,8), rep(5,8), rep(6,3904))
b <- c(rep(0,12), rep(1,17), rep(2,23), rep(3,35), rep(4,12), rep(5,27), rep(6,3881))
x <- c(rep('placebo', length(a)), rep('epinephrine', length(b)))
y <- c(a, b)

d = data.frame(x, y)

d$y <-factor(d$y, ordered = T)
d$x<-factor(d$x, levels = c("placebo", "epinephrine"))

Fit model with deafult vague priors

m <- brm(
  family = cumulative("logit"),
  formula = y ~ 1 + x,
  data = d,
  backend = "cmdstanr",
  cores = 4
)
#> Start sampling

Bayesian Cumulative Odds Ratio (mean and 95% CrI):

model_parameters(m, centrality = "mean") |> 
  tail(1) |>
  data.frame() |> 
  select("Mean", "CI_low", "CI_high") |> 
  mutate(across(1:3, ~exp(.)))
#> Possible multicollinearity between b_Intercept[3] and b_Intercept[2] (r = 0.81), b_Intercept[4] and b_Intercept[2] (r = 0.72), b_Intercept[5] and b_Intercept[2] (r = 0.71), b_Intercept[4] and b_Intercept[3] (r = 0.9), b_Intercept[5] and b_Intercept[3] (r = 0.88), b_Intercept[6] and b_Intercept[3] (r = 0.85), b_xepinephrine and b_Intercept[4] (r = 0.73), b_xepinephrine and b_Intercept[5] (r = 0.75), b_xepinephrine and b_Intercept[6] (r = 0.77). This might lead to inappropriate results. See 'Details' in '?rope'.
#>        Mean    CI_low   CI_high
#> 7 0.7156337 0.5378477 0.9450243

which is very similar to the frequentist one:

lrm(y ~ x, data = d) |> 
  model_parameters() |>
  tail(1) |> 
  data.frame() |> 
  select("Coefficient", "CI_low", "CI_high") |> 
  mutate(across(1:3, ~exp(.)))
#>   Coefficient    CI_low  CI_high
#> 7   0.7140849 0.5428909 0.939263

Back to Bayesian results, here is the code to estimate 95% CrIs and posterior probabilities on the risk difference at each mRS score

cri = 
  m |> 
  add_epred_draws(newdata = data.frame(x = c("placebo", "epinephrine"))) |> 
  pivot_wider(values_from = c(.row, .epred), names_from = x) |> 
  mutate(estimate = 100*(.epred_epinephrine - .epred_placebo)) |> 
  group_by(.category) |> 
  median_hdi(estimate) |> 
  select(1:4) |> 
  mutate(across(2:4, ~round(., 2)))

prob = 
  m |> 
  add_epred_draws(newdata = data.frame(x = c("placebo", "epinephrine"))) |> 
  pivot_wider(values_from = c(.row, .epred), names_from = x) |> 
  mutate(estimate = 100*(.epred_epinephrine - .epred_placebo)) |> 
  group_by(.category) |> 
  summarise("Pr(> 0)" = mean(estimate > 0),
            "Pr(> 0.1%)" = mean(estimate > 0.1),
            "Pr(< -0.2%)" = mean(estimate < -0.2)) |> 
  mutate(across(2:4, ~round(., 2)))

left_join(cri, prob) |> 
  rename("mRS" = .category,
         "Median" = estimate,
         "2.5th" = .lower,
         "97.5th" = .upper) |> 
  gt()
#> Joining, by = ".category"

Model diagnostics

Trace plots

plot(m)

Posterior Predictive Check

pp_check(m, group = "x", type = "bars_grouped", freq = F) 
#> Using 10 posterior draws for ppc type 'bars_grouped' by default.

^{Created on 2022-08-08 by the reprex package (v2.0.1)}