rmsb Random effect model

gistfile1.txt

pacman::p_load(rmsb,
               dplyr,
               magrittr,
               ggdist,
               tidyr,
               brms)

# Data
a <- c(rep(0,0), rep(1,3), rep(2,5), rep(3,5), rep(4,25), rep(5,40), rep(6,24))
b <- c(rep(0,2), rep(1,3), rep(2,9), rep(3,17), rep(4,33), rep(5,18), rep(6,18))
x <- c(rep('medical', length(a)), rep('endovascular', length(b)))
y <- c(a, b)
d1 = data.frame(x, y, study = "RESCUE-Japan LIMIT")


a <- c(rep(0,0), rep(1,3), rep(2,9), rep(3,20), rep(4,36), rep(5,32), rep(6,71))
b <- c(rep(0,2), rep(1,9), rep(2,25), rep(3,31), rep(4,27), rep(5,15), rep(6,68))
x <- c(rep('medical', length(a)), rep('endovascular', length(b)))
y <- c(a, b)
d2 = data.frame(x, y, study = "SELECT2")

a <- c(rep(0,0), rep(1,8), rep(2,18), rep(3,49), rep(4,60), rep(5,45), rep(6,45))
b <- c(rep(0,9), rep(1,19), rep(2,41), rep(3,39), rep(4,45), rep(5,27), rep(6,50))
x <- c(rep('medical', length(a)), rep('endovascular', length(b)))
y <- c(a, b)
d3 = data.frame(x, y, study = "ANGEL-ASPECT")

d = rbind(d1, d2, d3)

d$y <-factor(d$y, ordered = T)
d$x<-factor(d$x, levels = c("endovascular", "medical"))
dd <- datadist(d); options(datadist='dd')

## Proportional Odds Model----

# brms
PO_brms <- brms::brm(
  family = cumulative("logit"),
  formula = y ~ x + (1|study) ,
  prior = prior(normal(0, 1.5), class = "b"),
  data = d,
  backend = "cmdstanr",
  cores = 4,
  seed = 123
)


PO_brms

# rmsb

PO_rmsbb <- rmsb::blrm(y ~ x + cluster(study), 
                      priorsd = 1.5, 
                      seed = 123,
                      data=d)

PO_rmsb