Chris,

Thanks a lot for the help and your investigation by simulations.

> 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)

This formulation works but doesn't provide the beta estimate of Tx on the

second event type (we need to sum the two beta). An equivalent formulation

that provides beta estimates on the two types of events is to create

intermediate variables as proposed page 177 of Therneau and Grambsch and

attached here.

<

https://r.789695.n4.nabble.com/file/t388337/example.png>

> 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").

To clarify the context of this analysis, we are interested to estimate

(&test) the effect of a genetic variant on two time-to-complication traits

(non-competitive based on expert's knowledge), while accounting for

potential unexplained dependency between the traits. Multi-state models can

effectively be an interesting

alternative model formulation that I need to think a bit more in this

context.

> 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).

I modified your simulation script to plot the beta estimates from r=100 data

replicates obtained from the following model comparisons:

modf : Cox PH frailty model (model of interest)

modc : Cox PH including the simulated gamma variable (for comparison with

modf )

modcf : Cox PH frailty model + gamma (for comparisons with modf & modc; the

introduction of the simulated gamma variable should reduced the variance of

the frailty term; )

modsep: separate Cox PH model for each time-to-event trait (no frailty, no

gamma)

It Overall, I got some expected results:

-modf reduces the bias of Tx effect on event 1 and event2 compared to modsep

(see boxplots)

-modsep show an increased bias, while not for modf when I increase the

variance of the simulated gamma term (see boxplots)

-when the gamma simulated covariate is introduced as a covariate in the

model (modc), the variance of the frailty term is reduced (modcf versus modf

results; see for the last replicates).

So, it validates some conclusions I had before, but :

1) modf returns an estimate of the variance of the frailty term that tends

to be lower than the one specified for the data simulation

2) To be sure the specified model fits exactly what I want, I need to check

the likelihood function especially because they use an equivalence with a

penalized likelihood for the gamma frailty (see help of frailty()).

Regarding the model I specified for the two time-to-event traits, and the

conditional assumption I assume (the two time-to-event traits are assumed

independent given the frailty term),

the contribution of each subject i to the likelihood of the model of

interest should have the form:

Li = Li_ev1(beta1|ui) * Li_ev2(beta2|ui) * p(ui|zz),

with ui = subject's frailty term and ui~gamma(shape,scale) & p the gamma

distribution

Li_ev1 and Li_ev2 are the contribution of i to the likelihood for event1 &

event2

There is a description of the frailty model in Therneau's et al paper:

https://www.mayo.edu/research/documents/frailtypdf/doc-10027273 starting

page 15, detailed page 38.

The formulation is described in a different context, with N individuals

(1<=i<=N) in q independent clusters j (1<=j<=q) with a frailty term that

accounts for the dependence between individuals within each cluster.

In my case, each cluster corresponds to an individual i and dependencies

between the two observations for each subject are accounted by the frailty

term. Using their notations, j= individual and i= type of event.

I still need to clarify if their likelihood covers the likelihood of my

model formulation.

Any comments, suggestions or experiences are welcome. :)

If some have other suggestions of packages to fit this type of model,

please let me know (other potential R packages: parfm, frailtySurv, cph from

rms, mets, MST,..?).

Denise,

I attached below the extension of Chris's code.

#######################

# Modified Chris's code

#######################

library(flexsurv)

library(survival)

set.seed(20190729)

res <- NULL

for(r in 1:100) {

# 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 <- 2^2

#tau <- 0.4^2

tau <- 0.8^2

#tau <- 1^2

#tau <- 0.6^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)

coef(mod0)

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

coef(mod1)

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

coef(mod2)

# Transform long data & create contrasts terms

# create tx_ev1 & tx_ev2 to estimate effects of TX on event 1 & on event 2

long2 <- long

long2$event <- long$event-1

long2$tx_ev1 <- long2$tx*(1-long2$event)

long2$tx_ev2 <- long2$tx*(long2$event)

# mod1 corresponds to separate analysis of each time-to-event trait

modsep <- coxph(Surv(times)~tx_ev1 + tx_ev2 + strata(event), data=long2)

coef(modsep) # should match exactly Chris's results for mod1 & mod2

coef(mod1) # for verification

coef(mod2) # for verification

# Frailty model for time-to-event1 & time-to-event2 (frailty term captures

unexplained dependency)

modf <- coxph(Surv(times)~ tx_ev1 + tx_ev2 + strata(event) + frailty(id,

distribution="gamma"), data=long2)

coef(modf) # should match results of coef(mod)

#### for comparisons

# Corresponds to full model (includes simulated gamma as covariate, for

comparisons)

modc <- coxph(Surv(times)~tx_ev1 + tx_ev2 + gamma + strata(event) ,

data=long2) # corresponds to "true" model

coef(modc)

# Includes frailty term + gamma as a covariate (for comparisons)

modcf <- coxph(Surv(times)~tx_ev1 + tx_ev2 + gamma + strata(event) +

frailty(id, distribution="gamma") , data=long2) # corresponds to "true"

model

coef(modcf)

res <- rbind(res,list(coef(modc), coef(modf),coef(modcf),coef(modsep)))

}

pdf(paste("Boxplot.beta.event1.event2.tau",tau,".pdf",sep=""))

par(mfrow=c(2,4))

#boxplot for beta of Tx on event1

beta.tx.ev1.TM <- unlist(lapply(res[,1],function(x) {x[1] }))

boxplot(beta.tx.ev1.TM, main="Complete Model",ylab="beta of TX on event

1",col="grey")

abline(h=beta[1],col="red")

mean(beta.tx.ev1.TM)- beta[1] #bias

beta.tx.ev1.FM <-unlist(lapply(res[,2],function(x) {x[1] }))

boxplot(beta.tx.ev1.FM, main="Frailty Model",ylab="beta of TX on event

1",col="grey")

abline(h=beta[1],col="red")

mean(beta.tx.ev1.FM)- beta[1] #bias

beta.tx.ev1.FGM <-unlist(lapply(res[,3],function(x) {x[1] }))

boxplot(beta.tx.ev1.FGM, main="Frailty Model with sim. gamma

covariate",ylab="beta of TX on event 1",col="grey")

abline(h=beta[1],col="red")

mean(beta.tx.ev1.FGM)- beta[1] #bias

beta.tx.ev1.SM <- unlist(lapply(res[,4],function(x) {x[1] }))

boxplot(beta.tx.ev1.SM,main="Separate Models",ylab="beta of TX on event

1",col="grey")

abline(h=beta[1],col="red")

mean(beta.tx.ev1.SM)- beta[1] #bias

#beta of TX on event 2

beta.tx.ev2.TM <- unlist(lapply(res[,1],function(x) {x[2] }))

boxplot(beta.tx.ev2.TM,main="Complete Model",ylab="beta of TX on event

2",col="grey")

abline(h=beta[2],col="red")

mean(beta.tx.ev2.TM)- beta[2] #bias

beta.tx.ev2.FM <-unlist(lapply(res[,2],function(x) {x[2] }))

boxplot(beta.tx.ev2.FM,main="Frailty Model",ylab="beta of TX on event

2",col="grey")

abline(h=beta[2],col="red")

mean(beta.tx.ev2.FM)- beta[2] #bias

beta.tx.ev2.FGM <-unlist(lapply(res[,3],function(x) {x[2] }))

boxplot(beta.tx.ev2.FGM,main="Frailty Model with sim. gamma

covariate",ylab="beta of TX on event 2",col="grey")

abline(h=beta[2],col="red")

mean(beta.tx.ev2.FGM)- beta[2] #bias

beta.tx.ev2.SM <- unlist(lapply(res[,3],function(x) {x[2] }))

boxplot(beta.tx.ev2.SM,main="Separate Models",ylab="beta of TX on event

2",col="grey")

abline(h=beta[2],col="red")

mean(beta.tx.ev2.SM)- beta[2] #bias

dev.off()

# Print summary for the last data replicate

summary(modc)

summary(modf)

summary(modcf)

summary(modsep)

Le lun. 29 juil. 2019 à 12:00, Andrews, Chris <

[hidden email]> a

écrit :

> 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:

[hidden email]]

> Sent: Sunday, July 28, 2019 2:03 AM

> To:

[hidden email]
> 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]]

> >

> > ______________________________________________

> >

[hidden email] 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

> ______________________________________________

>

[hidden email] 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.

>

[[alternative HTML version deleted]]

______________________________________________

[hidden email] mailing list -- To UNSUBSCRIBE and more, see

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.