I want to fit a Bayesian Constrained Partial Proportional Odds Model with rmsb. In the example below, my dependent variable is a 7-level ordinal score (y) from 0 to 6. Y represents the observed outcome, and X the exposure group (tx = placebo or epinephrine) and another covariate “age”.
I am interested in these estimands:
- P(Y = y \mid X = tx, age)
- P(Y \geq y \mid X = tx, age)
For example, I want to have access to all posterior draws of P(Y = 6) and P(Y \geq 2) (one vector for each estimand) given tx = placebo and age = 60.
I have thoroughly read all rmsb-related vignettes kindly designed by @f2harrell, but I haven’t figured out how to get access to these estimands, as explained above.
Can anyone help me to write the code corresponding to this task? See below the R code I have generated thus far.
Thank you in advance!
library(rmsb)
library(dplyr)
a <- c(rep(0,15), rep(1,10), rep(2,29), rep(3,20), rep(4,8), rep(5,8), rep(6,3903))
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"))
# Data
a <- c(rep(0,15), rep(1,10), rep(2,29), rep(3,20), rep(4,8), rep(5,8), rep(6,3903))
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)
# Simulated covariate
set.seed(123)
age_p = rnorm(4000, 69.7, 16.6)
age_e = rnorm(4000, 69.8, 16.4)
age = c(age_p, age_e)
d = data.frame(x, y, age)
# Reduce sample size
dd = sample_frac(d, size = 0.1)
model <- blrm(y ~ 1 + x + age,
ppo = ~ x,
cppo=function(y) y == 6,
priorsd = 1.5,
priorsdppo = 0.246,
seed = 123,
iter = 2000,
chains = 4,
data = dd)