[R] CoxPH multivariate frailty model with coxph (survival)

[Andrews, Chris]

Denise, 

Interacting a fixed effect with a strata variable is not uncommon.  No main effect of the strata variable is possible but the interaction term allows the fixed effect variable to have different values in different strata.  I've more often seen it coded as a direct interaction:

coxph( Surv(Time, Status) ~ treatment*strata(Event_type) + frailty(ID), data=example)

(e.g., page 47 of Therneau and Grambsch)

I don't see anything inherently wrong with your interpretation of the model.  While it seems that the assumption of no effect of one event on the other is very strong, I don't know the context of your analysis. I visualized it as a 4-state model.  Everyone starts in 0 ("neither event").  There are hazards of moving from state 0 to state 1 (lambda01(t)) and from state 0 to state 2 (lambda02(t)).  The last state is "both events".  There are hazards of moving from state 1 to state 3 (lambda13(t), which you are assuming is identical to lambda02(t)) and from state 2 to state 3 (lambda23(t), which you are assuming is identical to lambda01(t)).

I did not observe much gain from the frailty term (unmeasured covariate) with only 2 events in the short simulation I tried (except when the effect was very strong and I got convergence warnings).

Chris


library("flexsurv")

set.seed(20190729)

# Multiple non-competing outcomes, connected only by frailty (unmeasured covariate)

nn <- 1000
kk <- 2

# frailty, 1 per individual
# variance of random effect
#tau <- 0.01^2 
tau <- 0.1^2
#tau <- 0.4^2
#tau <- 1^2
gamma <- rgamma(nn, shape=1/tau, scale=tau) # mean is 1
hist(gamma)
sd(gamma) # ~ tau

# covariate
tx <- sample(0:1, nn, replace = TRUE)
# covariate effect on log hazard
beta <- 1
# might want to allow different treatment effects for different events
beta <- seq(kk)/kk

short <- data.frame(id=seq(nn), tx=tx, gamma=gamma)

# survival times, kk per individual
# tt <- rweibullPH(kk*nn, shape=2, scale=exp(beta*tx)*gamma)
# hist(tt)
# might want to allow different shapes or scales for different events
for (k in seq(kk)) {
	short[[paste0("time", k)]] <- rweibullPH(nn, shape=2, scale=exp(beta[k]*tx)*gamma)
}
# might want to allow censoring

long <- reshape(short, direction="long", varying = paste0("time", seq(kk)), v.names="times", timevar="event", idvar="id", times=seq(kk))
long <- long[order(long$id, long$event),]

mod <- coxph(Surv(times) ~ tx * strata(event) + frailty(id), data=long)
#summary(mod)

mod0 <- coxph(Surv(times) ~ tx * strata(event), data=long)
#summary(mod0)

mod1 <- coxph(Surv(times) ~ tx, data=long, subset=event==1)
#summary(mod1)

mod2 <- coxph(Surv(times) ~ tx, data=long, subset=event==2)
#summary(mod2)

coef(mod)
coef(mod0)
coef(mod1)
coef(mod2) - coef(mod1)

coef(summary(mod))
coef(summary(mod0))
coef(summary(mod1))



-----Original Message-----
From: David Winsemius [mailto:dwinsemius using comcast.net] 
Sent: Sunday, July 28, 2019 2:03 AM
To: r-help using r-project.org
Subject: Re: [R] CoxPH multivariate frailty model with coxph (survival)


On 7/19/19 10:19 AM, Denise b wrote:
>   Dear R users,
>
> I am interested in estimating the effects of a treatment on two
> time-to-event traits (on simulated data), accounting for the dependency
> between the two time-to-event outcomes.
>
> I precise that the events are NOT recurrent, NOT competitive, NOT ordered.
> The individuals are NOT related and can have 0, 1 or the 2 events (in any
> ordered).
>
> So, I specified a time-to-event model with one Cox PH hazard function for
> each outcome.
> The two hazard functions are linked by a common subject-specific frailty
> term (gamma distribution) to account for the dependency. The likelihood
> function of that model would include two risks sets (1 for eachoutcome)
> connected via the shared frailty term.
>
> To fit that model, I used the coxph function (survival R package, T.
> Thernaud):
> coxph( Surv(Time, Status) ~ treatment*Event_type + strata(Event_type) +
> frailty(ID), data=example)
If you have Event_type as a strata variable, it seems problematic that 
it be also interacting as a fixed effect.
> - where example is my dataset, with 2*N individuals (2 rows for each
> individual, corresponding to each time-to-event outcome)
> - Time = c(Time_outcome1, Time_outcome2)
How exactly do you expect the Surv-function to handle such a vector?
> - Status=c(Event_outcome1, Event_outcome2)

The help page for Surv() says:  "For multiple endpoint data the event 
variable will be a factor, whose first level is treated as censoring."

So it, too, would not be "expecting" a two-item vector. Seems that the 
multi-state formulation would make more sense.

> - Event_type=c(rep(0,length(Time_outcome1)),rep(1,length(Time_outcome2))


I'm reasonably sure 0 would be interpreted as a censored observation.

> - ID=c(ID,ID)
>
> So, the model implies different baseline hazard function for each outcome
> (argument strata) and estimate a treatment effect for each event type(
> treatment*Event_type).


I doubt it. I'm pretty sure there would be a conflation of effects and 
stratification.

>
> Although the model returns some estimates close to what I expected
> (simulation study), the problem is that most of the examples I found with
> this type of formulation are for competitive time-to-event models (as
> presented in Therneau's and Grambush's book "Modeling survival
> data"  (pages 240-260 and 175-185) and also in some info I found on


I cannot lay my hands on my copy of that text, but I'm fairly sure it 
never mixed strata(.) and fixed effects for a variable in the same model.


I suppose this is where I should admit  that I'm trained as an 
epidemiologist and not as a statistician, but I've done a fair amount of 
work with these models. Seems to me that your code are sufficiently out 
of the mainstream practice that you would at least want to create a 
simulated data-set and see if this results in agreement with assumptions.


-- 

David.

> internet).
>
> So, I am wondering if some of you have the same or different interpretation
> of mime.
> Otherwise do you have other functions to recommend me to fit the desired
> model, such as coxme...? I also checked the package " mets" which also
> describes several examples for competitive time-to-event traits, but not
> for NON competitive events.
>
> Thanks in advance for your help,
> Denise
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


**********************************************************
Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues