Nonlinear effects of time-varying covariates from marginal rates model

Question

This is a followup analysis on this post, inspired by the comments from @EdM. I fitted a marginal rates model (Lin, Wei, Yang, & Ying, 2002) on recurrent event data but don't know how to include and then intepret the nonlinear effects of time-varying covariates. Here is a reproducible example:

library(survival)
library(survsim) #package to simulate survival data
N=100 #number of patients
set.seed(123)
df.tf<-simple.surv.sim(#baseline time fixed
  n=N, foltime=500,
  dist.ev=c('llogistic'),
  anc.ev=c(0.68), beta0.ev=c(5.8),
  anc.cens=1.2,
  beta0.cens=7.4,
  z=list(c("unif", 0.8, 1.2)),
  beta=list(c(-0.4),c(0)),
  x=list(c("bern", 0.5),
         c("normal", 70, 13)))
names(df.tf)[c(1,6,7)]<-c("id","grp","age")
nft<-sample(1:10, N,replace=TRUE)#number of follow up time points
crp<-round(abs(rnorm(sum(nft)+N,
                     mean=100,sd=40)),1)
time<-NA
id<-NA
i=0
for(n in nft){
  i=i+1
  time.n<-sample(1:500,n)
  time.n<-c(0,sort(time.n))
  time<-c(time,time.n)
  id.n<-rep(i,n+1)
  id<-c(id,id.n)
}
df.td <- cbind(data.frame(id,time)[-1,],crp) #time-varying covariate
df<-tmerge(df.tf,df.tf,id=id,
           endpt=event(stop,status)) 
df <- tmerge(df,df.td,id=id,
             crp=tdc(time,crp))
df <-df[,c(1,6:11)] 
#fit marginal rates model:
model.fit<-coxph(Surv(tstart, tstop, endpt) ~ grp + age + crp + cluster(id), method="breslow", data = df)
summary(model.fit)
Call:
coxph(formula = Surv(tstart, tstop, endpt) ~ grp + age + crp, 
    data = df, method = "breslow", cluster = id)
n= 378, number of events= 67
     coef exp(coef)  se(coef) robust se     z Pr(&gt;|z|)  

grp 0.5011417 1.6506047 0.2525867 0.2553780 1.962   0.0497 *
age 0.0004950 1.0004951 0.0081309 0.0072486 0.068   0.9456

crp 0.0009027 1.0009031 0.0027431 0.0023715 0.381   0.7035

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95

grp     1.651     0.6058    1.0006     2.723
age     1.000     0.9995    0.9864     1.015
crp     1.001     0.9991    0.9963     1.006
Concordance= 0.553  (se = 0.04 )
Likelihood ratio test= 4.22  on 3 df,   p=0.2
Wald test            = 4.38  on 3 df,   p=0.2
Score (logrank) test = 4.19  on 3 df,   p=0.2,   Robust = 4.58  p=0.2
(Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).

Question 1): If I suspect the effect of grp is nonlinear, how do I test that and maybe plot it? Do I need to worry about the violation of proportional hazards assumption? Question 2): How do I interpret the coefficients? Are those mean rates, not hazard rates?

Answer 1

This seems to be an Andersen-Gill repeated events model, with a coefficient covariance matrix based on clusters to take into account repeated measures for inference. See Section 3.2 of the R survival vignette.

Coefficient interpretation

The coefficients themselves aren't hazards or rates. Each is the estimated log of the hazard ratio for an event, given a unit change in the covariate. Those hazard ratios work around a baseline cumulative hazard that can be estimated after the model is fit. In this analysis there is a single baseline cumulative hazard assumed to apply to all events and individuals.

Nonlinearity and proportional hazards

You can flexibly model continuous covariates with the pspline() function in the survival package. There's an example in Section 3.1 of the above-referenced vignette, showing how to plot predictions and test for non-linearity. You also could use regression splines.

If you are doing proportional hazards modeling, you do have to examine proportional hazards (PH). Remember that, at all event times, you are modeling the current hazard associated with the current value of a covariate, regardless of its prior values. Thus having time-varying covariate values doesn't provide a dispensation from dealing with the PH assumption. You are assuming that the association of a covariate with event hazard is the same for all times evaluated.

That said, with a large model it's possible to get a statistically significant violation of PH that is of minor practical significance. You have to apply your understanding of the subject matter to evaluate.