How to interpret the association values in R stan_jm fitted models for clinical use

Modeling
1

Hi everyone! I am new to this forum so hopefully I am following the appropriate etiquette.

I have a dataset of N = 134 subjects. 21 of these will have an event by the end of the longitudinal follow-up period. The time to event is measured in months from the baseline visit, coded by the variable months_from_baseline.

The prognosis is usually computed at the baseline visit using a score (integer, 0 to 6, higher score = worse prognosis), here coded in the variable prognostic_baseline_score. It makes sense to predict the event in the event submodel through an univariable Cox regression using prognostic_baseline_score as sole predictor. The median prognostic_baseline_score is 1 (IQR: 0-6).

In this study, a biomarker score (coded by the variable biomarker_score is measured at the baseline visit and every 3 months after the baseline visit. biomarker_score is a continuous variable scaled in such way that a 1-unit increase in biomarker_score corresponds to a clinically meaningful change in the underlying biomarker. The mean biomarker_score across longitudinal visit is 1.32 in those who will have an event and 0.52 in those who will not (SD = 1.2).

Here a dataset summary at the baseline visit:

Screenshot 2024-04-06 at 17.48.03
Screenshot 2024-04-06 at 17.48.031248×536 36.6 KB

The hypothesis is that a high biomarker score worsen the prognosis of the subject independent of the prognostic_baseline_score.

I used the following joint model with a current association structure:

mod1 <- stan_jm(
    #Longitudinal submodel: biomarker score depends on 2-knot spline subject-specific intercept & slope + global intercept & slope
    formulaLong = biomarker_score ~ 1 + months_from_baseline_visit + (1 + ns(months_from_baseline_visit,2) | subject_id),
    dataLong = df.toy,
    #Event submodel: a 3-B splines approximation of the survival function depends on the baseline_prognostic_score
    formulaEvent = survival::Surv(months_to_event, event) ~ baseline_prognostic_score, 
    dataEvent = df.toy.baseline,
    family = gaussian,
    basehaz = "bs",
    basehaz_ops = list(df = 3),
    #Current association structure
    assoc = c("etavalue"),
    id_var = "subject_id",
    time_var = "months_from_baseline_visit",
    chains = 3,
    cores = 3, 
    refresh = 200, 
    iter = 5000,
    seed = seed)

The summary ouput of this model is:

summary(mod1, probs = c(.025,.975))

Screenshot 2024-04-06 at 17.50.50
Screenshot 2024-04-06 at 17.50.501920×1764 231 KB

My question is - is it correct to interpret that each 1-unit increase in a biomarker score observed during follow-up is associated with a HR = exp(0.477) = 1.61 independent of the baseline_prognostic_score and with a 95% credible interval between exp(0.096) = 1.10 and exp(0.828) = 2.29?

I am using R version 4.3.2 (2023-10-31) with rstanarm_2.32.1
Thank you so much for your help!

2

So sorry that nobody responded to this question. I don’t know rstanarm, but I can help a bit. @bgoodri is the one to answer both about the stats and about rstanarm.

Not quite. You are only seeing the reports of the marginal posteriors here. The issue is that there can be correlation among the parameters in the posteriors, which proper predictive inference will take into account. For example, in a simple linear regression, I might have a univariate slope and an intercept parameter. If all the covariates are positive, you’re going to have strong (negative) correlations between the intercept and slope. If you just treated the intercept and slope as independent, you’d get the wrong posterior inferences—instead, the negative correlations will make sure you’re not getting inconsistent values of the slope and intercept.

Also, whenever you interpret a regression, it’s always relative to all the other covariates being held constant. If you were to use a different set of interactions or other covariates, the effect of baseline_prognostic_score will be different.