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:dwinsem...@comcast.net]
Sent: Sunday, July 28, 2019 2:03 AM
To: r-help@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@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
______________________________________________
R-help@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.
