# Nonlinear effects of time-varying covariates from marginal rates model

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(>|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?

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.

• Thank you for the clarification @EdM. I was following the code by Amorim & Cai (2015). In their supplementary data, they specified the Marginal means and rates model as model.2 = coxph(Surv(tstart,tstop,status) ~ var1 + var2+ … + vark + cluster(id), method=”breslow”, data = example1). I also wonder if their code was erroneous as the only difference in the code between AG and marginal means models is the addition of cluster(id). So do you know the correct way to specify marginal means model?
– cliu
Commented Nov 2, 2022 at 1:15
• @cliu the distinction here is between marginal and conditional models, as Amorim & Cai discuss in the text. The"frailty model" they distinguish from a "marginal means/rates model" is a type of conditional model, as are mixed models with random effects. Amorim & Cai have it right for a marginal model; including the cluster term is a way to adjust the coefficient covariance matrix in a marginal model for the dependence of observations within individuals.
– EdM
Commented Nov 2, 2022 at 13:10
• @cliu see this page for an introduction to the differences in what's estimated by marginal versus conditional models. The generalized estimating equation (GEE) method mentioned there is another way to handle within-subject correlations in a marginal model.
– EdM
Commented Nov 2, 2022 at 13:16
• Thanks @EdM. I read some other references after I posted my last comment. Indeed that was the code for marginal means/rates model.
– cliu
Commented Nov 2, 2022 at 14:50