Survival models with time-varying specification: the first covariate will be modeled as time-varying variable

Summary:

Time-varying covariate specification seems to always apply to the first covariate in the model.

Description:

model:

mod = survival object ~ p2 + p3 + tve(p1)

p2 is modelled as a time-varying covariate.

However, changing covariate order will fix the problem.

mod2 = survival object ~ tve(p1) + p2 + p3

Thus, it seems that tve()-specified covariates always have to the first in the model.

Reproducible Steps:

#data, its size reduzing and standardization of 2 variables
library(tidyverse)
df = read_delim("https://www.dropbox.com/s/ufpmejd6do8saq6/shared_data.csv?dl=1", ";", escape_double = FALSE, trim_ws = TRUE) %>% mutate(p3 = scale(p3), p4 = scale(p4)) %>%  group_by(p1) %>% sample_n(50)

This model specification seems to work as supposed to.

library(rstanarm)
options(mc.cores = parallel::detectCores())
#model 1
m1 = stan_surv(Surv(time, status) ~ tve(p1) + p3 + p4, df)
m1

stan_surv
 baseline hazard: M-splines on hazard scale
 formula:         Surv(time, status) ~ tve(p1) + p3 + p4
 observations:    100
 events:          45 (45%)
 right censored:  55 (55%)
 delayed entry:   no
------

                 Median MAD_SD exp(Median)
(Intercept)      -0.6    0.2     NA       
**p1B               0.2    0.4    1.2**       
p3                0.6    0.2    1.8       
p4                0.3    0.1    1.4       
**p1B:tve-bs-coef1 -0.3    0.7    0.7       
p1B:tve-bs-coef2 -0.5    0.7    0.6       
p1B:tve-bs-coef3 -0.6    0.8    0.6**       
smooth_sd[p1B]    0.5    0.5    1.7       
m-splines-coef1   0.1    0.0     NA       
m-splines-coef2   0.2    0.1     NA       
m-splines-coef3   0.2    0.1     NA       
m-splines-coef4   0.2    0.1     NA       
m-splines-coef5   0.1    0.1     NA       
m-splines-coef6   0.1    0.1     NA       

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

This models the first covariate p4 as time-varying.

#model 2
m2 = stan_surv(Surv(time, status) ~ p4 + p3 + tve(p1), df)
m2

There were 5 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Examine the pairs() plot to diagnose sampling problems
stan_surv
 baseline hazard: M-splines on hazard scale
 formula:         Surv(time, status) ~ p4 + p3 + tve(p1)
 observations:    100
 events:          45 (45%)
 right censored:  55 (55%)
 delayed entry:   no
------
                Median MAD_SD exp(Median)
(Intercept)     -0.7    0.2     NA       
**p4               0.3    0.2    1.4**       
p3               0.6    0.2    1.8       
p1B              0.0    0.3    1.0       
**p4:tve-bs-coef1 -0.2    0.5    0.8       
p4:tve-bs-coef2 -0.1    0.5    0.9       
p4:tve-bs-coef3  0.3    0.5    1.3**       
smooth_sd[p4]    0.5    0.5    1.7       
m-splines-coef1  0.1    0.0     NA       
m-splines-coef2  0.2    0.1     NA       
m-splines-coef3  0.2    0.1     NA       
m-splines-coef4  0.2    0.1     NA       
m-splines-coef5  0.1    0.1     NA       
m-splines-coef6  0.1    0.1     NA       

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

tve plot for the first covariate p4

m2 %>% plot(, plotfun = "tve", pars = "p4")

RStanARM Version:

‘2.21.2’

R Version:

‘4.0.4’

Operating System:

Windows 10 Enterprise, 21H1, 19043.985, Windows Feature Experience Pack 120.2212.2020.0

Start by reproducing the two stan_surv formulas from the report with the supplied data and compare which covariate receives the time-varying terms. Trace the stan_surv formula handling and tve processing, then add or run a regression check showing that tve(p1) remains attached to p1 regardless of covariate order.