Estimands from Bayesian Constrained Partial Proportional Odds Model in rmsb

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)

1 Like

I have found a way to extract all posterior draws for any P(Y = y \mid X = tx, age) of interest (code below), thanks to the rms::contrast() function.

Still need help with P(Y \geq y \mid X = tx, age) though.

Anyone?

# Install/load packages
pacman::p_load(ggdist,
               dplyr,
               tidyr,
               rmsb)

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

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)

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

dd = sample_frac(d, size = 0.1, weight = x)

# Fit Bayesian Constrained Partial Proportional Odds Model
b_adj <- blrm(y ~ 1 + x + age,
              ppo = ~ x,
              cppo=function(y) y == 6,
              priorsd = 1.5,
              priorsdppo = 0.246,
              seed = 123, 
              iter = 1000,
              chains = 4,
              data = dd)


nd = data.frame(x = c("placebo", "epinephrine"),
               age = 60)

# y = {0, 1, 2, 3, 4, 5, 6}

# P( Y = y | x = placebo, age = 60)
placebo_epred = 
  rms::contrast(b_adj,
                subset(nd, x == "placebo"),
                ycut=0:6)

# P( Y = y | x = epinephrine, age = 60)
epinephrine_epred = 
  rms::contrast(b_adj,
                subset(nd, x == "epinephrine"),
                ycut=0:6)

# Extract posterior draws (1 column per y)

placebo_epred_draws = as.data.frame(placebo_epred$cdraws)
epinephrine_epred_draws = as.data.frame(epinephrine_epred$cdraws)

# Rename columns according to y
outcome = 
  c("No symptoms",                   # y = 0
  "No significant disability",       # y = 1
  "Slight disability",               # y = 2
  "Moderate disability",             # y = 3
  "Moderately severe disability",    # y = 4
  "Severe disability",               # y = 5
  "Dead")                            # y = 6

colnames(placebo_epred_draws) = outcome
colnames(epinephrine_epred_draws) = outcome

# Transform from wide to long format
placebo_epred_draws = placebo_epred_draws |> tidyr::pivot_longer(dplyr::everything(), names_to = "y")
epinephrine_epred_draws = epinephrine_epred_draws |> tidyr::pivot_longer(dplyr::everything(), names_to = "y")

# Add column regarding "x" and join both data frame...
placebo_epred_draws$x = "placebo"
epinephrine_epred_draws$x = "epinephrine"
epred_draws = rbind(placebo_epred_draws, epinephrine_epred_draws)

# Arrange x and y column
epred_draws$x = factor(epred_draws$x, levels = c("placebo", "epinephrine"))
epred_draws$y = factor(epred_draws$y, levels = outcome)

# Show mean and 95% CI of both estimands (P( Y = y | x = {placebo, epinephrine}, age = 60))
epred_draws |> 
  dplyr::group_by(x, y) |>
  ggdist::mean_hdi(value) |>
  dplyr::select(-c(.width:.interval))
#> # A tibble: 14 × 5
#>    x           y                              value .lower .upper
#>    <fct>       <fct>                          <dbl>  <dbl>  <dbl>
#>  1 placebo     No symptoms                   0.0508 -0.346  0.397
#>  2 placebo     No significant disability     0.0508 -0.346  0.397
#>  3 placebo     Slight disability             0.0508 -0.346  0.397
#>  4 placebo     Moderate disability           0.0508 -0.346  0.397
#>  5 placebo     Moderately severe disability  0.0508 -0.346  0.397
#>  6 placebo     Severe disability             0.0508 -0.346  0.397
#>  7 placebo     Dead                          0.0508 -0.346  0.397
#>  8 epinephrine No symptoms                  -0.0341 -0.561  0.387
#>  9 epinephrine No significant disability    -0.0341 -0.561  0.387
#> 10 epinephrine Slight disability            -0.0341 -0.561  0.387
#> 11 epinephrine Moderate disability          -0.0341 -0.561  0.387
#> 12 epinephrine Moderately severe disability -0.0341 -0.561  0.387
#> 13 epinephrine Severe disability            -0.0341 -0.561  0.387
#> 14 epinephrine Dead                          0.0228 -0.397  0.444


# It matches summaries generated by rms::contrast()
placebo_epred
#>               x age   Contrast     S.E.      Lower     Upper Pr(Contrast>0)
#> 1   y=0 placebo  60 0.05082771 0.190282 -0.3459689 0.3970608         0.6175
#> 2*  y=1 placebo  60 0.05082771 0.190282 -0.3459689 0.3970608         0.6175
#> 3*  y=2 placebo  60 0.05082771 0.190282 -0.3459689 0.3970608         0.6175
#> 4*  y=3 placebo  60 0.05082771 0.190282 -0.3459689 0.3970608         0.6175
#> 5*  y=4 placebo  60 0.05082771 0.190282 -0.3459689 0.3970608         0.6175
#> 6*  y=5 placebo  60 0.05082771 0.190282 -0.3459689 0.3970608         0.6175
#> 7*  y=6 placebo  60 0.05082771 0.190282 -0.3459689 0.3970608         0.6175
#> 
#> Redundant contrasts are denoted by *
#> 
#> Intervals are 0.95 highest posterior density intervals
#> Contrast is the posterior mean
epinephrine_epred
#>                   x age    Contrast      S.E.      Lower     Upper
#> 1   y=0 epinephrine  60 -0.03409503 0.2443753 -0.5607075 0.3865915
#> 2*  y=1 epinephrine  60 -0.03409503 0.2443753 -0.5607075 0.3865915
#> 3*  y=2 epinephrine  60 -0.03409503 0.2443753 -0.5607075 0.3865915
#> 4*  y=3 epinephrine  60 -0.03409503 0.2443753 -0.5607075 0.3865915
#> 5*  y=4 epinephrine  60 -0.03409503 0.2443753 -0.5607075 0.3865915
#> 6*  y=5 epinephrine  60 -0.03409503 0.2443753 -0.5607075 0.3865915
#> 7   y=6 epinephrine  60  0.02280599 0.2207783 -0.3966771 0.4443049
#>         Pr(Contrast>0)
#> 1   y=0          0.454
#> 2*  y=1          0.454
#> 3*  y=2          0.454
#> 4*  y=3          0.454
#> 5*  y=4          0.454
#> 6*  y=5          0.454
#> 7   y=6          0.533
#> 
#> Redundant contrasts are denoted by *
#> 
#> Intervals are 0.95 highest posterior density intervals
#> Contrast is the posterior mean

Created on 2022-10-04 with reprex v2.0.2