Selection of patients for the PBC datasets

[wlandau]

From ?pbcSurv:

Longitudinal biomarker and time-to-event survival data for 40 patients with primary biliary cirrhosis who participated in a randomised placebo controlled trial of D-penicillamine conducted at the Mayo Clinic between 1974 and 1984.

The actual study includes many more patients. In fact, survival::pbc has 312 randomized patients and 106 non-randomized patients.

How were the 40 patients selected for pbcSurv and pbcLong?

[jgabry]

@sambrilleman do you know?

[sambrilleman]

Gosh I really can't remember. And I searched the repo and don't think the selection happens in there - I probably just committed a static subset of the full data. It was purely to run the code examples, not to make any reasonable interpretations about. If it matters, is there an ID or metric column you can use to triangulate against the original data? But again, you wouldn't want to interpret the results from fitting anything to that tiny dataset, it is just to use in the examples showing how to run the code!

[wlandau]

Thanks for clarifying. Indeed, I see at https://mc-stan.org/rstanarm/articles/jm.html#dataset-used-in-the-examples:

We use the Mayo Clinic’s primary biliary cirrhosis (PBC) dataset in the examples below. The dataset contains 312 individuals with primary biliary cirrhosis who participated in a randomised placebo controlled trial of D-penicillamine conducted at the Mayo Clinic between 1974 and 1984 [19]. However, to ensure the examples run quickly, we use a small random subset of just 40 patients from the full data.

Searching further, I see full versions of the same datasets in the survival package. According to the documentation of both packages, rstanarm::pbcSurv should be a subset of survival::pbc, and rstanarm::pbcLong should be a subset of survival::pbcseq.

I’ll just note that an analysis of the subset in rstanarm supports the opposite conclusion from the original Mayo Clinic study. In the original study by Dickson et al. (1985), penicillamine was found to be ineffective against primary biliary cholangitis, but the hazard ratio from modeling the subset suggests a strong effect. The hazard ratio centers at 1 if I analyze the full datasets from survival in the same way.

library(dplyr)
library(posterior)
library(rstanarm)
library(survival)

fit <- stan_jm(
  formulaLong = logBili ~ sex + trt + year + (year | id),
  dataLong = pbcLong,
  formulaEvent = Surv(futimeYears, event) ~ sex + trt,
  dataEvent = pbcSurv |>
    # futimeYears is time until death or transplant (or censoring)
    mutate(event = status > 0),
  time_var = "year",
  chains = 4,
  refresh = 4000,
  cores = 4,
  seed = 12345
)
fit |>
  as_draws_df() |>
  pull(`Event|trt`) |>
  exp() |>
  hist(
    main = "Hazard ratio from data subset in rstanarm",
    xlab = "hazard ratio",
    ylab = "Marginal posterior density",
    freq = FALSE
  )

References

Dickson, E. R., T. R. Fleming, R. H. Wiesner, W. P. Baldus, C. R. Fleming, J. Ludwig, and J. T. McCall. 1985. “Trial of Penicillamine in Advanced Primary Biliary Cirrhosis.” New England Journal of Medicine 312 (16): 1011–15. https://doi.org/10.1056/NEJM198504183121602.