[R] CoxPH multivariate frailty model with coxph (survival)
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) - 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) - Status=c(Event_outcome1, Event_outcome2) - Event_type=c(rep(0,length(Time_outcome1)),rep(1,length(Time_outcome2)) - 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). 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 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.
Re: [R] CoxPH multivariate frailty model with coxph (survival)
,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 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: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 + strat
Re: [R] Score test for a subset of parameters estimated with coxph
>
> Dear R help list,
>
> I would like to perform a score test for a subset of the parameters
> estimated with coxph using the frailty() option.
> As illustrated in the following reproducible example, I am able to perform
> the score test with the standard coxph() but not in presence of frailty()
> argument. See Examples 1 and 2
>
> library(survival)
>
> # Example 1: score test using standard coxph()
> # score test for ph.karno & pat.karno coefficients
> fit1<-coxph(Surv(time,status)~sex+ ph.karno + pat.karno , data=lung)
> fit0<-coxph(Surv(time,status)~pat.karno , data=lung) # fit under null
> scorestat <- coxph(Surv(time,status)~sex+ ph.karno + pat.karno,
> init=c(0,0,coefficients(fit0)) ,iter=0, data=lung)
> scorestat$score
> [1] 10.50409
>
> # Example 2: coxph with frailty term
> # score test for ph.karno & pat.karno coefficients
> frailtyfit1<-coxph(Surv(time,status)~sex+ ph.karno + pat.karno
> +frailty(inst), data=lung)
> frailtyfit0<-coxph(Surv(time,status)~sex+ frailty(inst), data=lung)
> scorefrailtyfit1<-coxph(Surv(time,status)~sex+ ph.karno + pat.karno
> +frailty(inst), data=lung,
> init=c(0,0,coefficients(frailtyfit0)) ,iter=0)
>
> I get the following error message:
> Error in coxph(Surv(time, status) ~ sex + ph.karno + pat.karno +
> frailty(inst), :
> wrong length for init argument
>
>
> When I checked in the code I found
> if (length(init) != ncol(X))
> stop("wrong length for init argument")
> Which seems that the length of the init argument doesnt match the number of
> predictors in the model.
>
> It seems that the code is expecting that I specify another parameter for
> frailty term.
> So, I tried to specify 2 as an inital value for frailty but got another
> error message later:
>
> scorefrailtyfit1<-coxph(Surv(time,status)~sex+ ph.karno + pat.karno
> +frailty(inst), data=lung,
> + init=c(0,0,coefficients(frailtyfit0),4.5),iter=0)
> Error in coxpenal.fit(X, Y, strats, offset, init = init, control, weights =
> weights, :
> Wrong length for inital values
>
> I have also checked coxme() but it seems that there is no option
> implemented
> for the score test.
> Do you have an idea of what could be the problem?
>
> Thanks,
> 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.
Re: [R] Score test for a subset of parameters estimated with coxph
An update regarding my question to compute the score test for a subset of
parameters with coxph,
when the option frailty is used:
After further explorations, based on my previous examples, I realized that
in presence of frailty term,
coxph returns only the results from the likelihood ratio test, whereas when
there is no frailty term in the model specified, it returns the likelihood
ratio, wald and score test as illustrated in the following example.
It could be due to the penalized-likelihood approach used by coxph for
frailty model in coxph.
Is there any other potential explanation?
# R code
library(survival)
fit1<-coxph(Surv(time,status)~sex+ ph.karno + pat.karno , data=lung)
summary(fit1)
Call:
coxph(formula = Surv(time, status) ~ sex + ph.karno + pat.karno,
data = lung)
n= 224, number of events= 161
(4 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
sex -0.511878 0.599369 0.169275 -3.024 0.00249 **
ph.karno -0.006155 0.993864 0.006822 -0.902 0.36695
pat.karno -0.017020 0.983124 0.006535 -2.604 0.00920 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
sex 0.5994 1.6680.43010.8352
ph.karno 0.9939 1.0060.98071.0072
pat.karno0.9831 1.0170.97060.9958
Concordance= 0.646 (se = 0.026 )
Rsquare= 0.097 (max possible= 0.999 )
Likelihood ratio test= 22.91 on 3 df, p=4.208e-05
Wald test= 22.64 on 3 df, p=4.789e-05
Score (logrank) test = 23.2 on 3 df, p=3.673e-05
frailtyfit1<-coxph(Surv(time,status)~sex+ ph.karno + pat.karno
+frailty(inst), data=lung)
summary(frailtyfit1)
Call:
coxph(formula = Surv(time, status) ~ sex + ph.karno + pat.karno +
frailty(inst), data = lung)
n= 223, number of events= 160
(5 observations deleted due to missingness)
coef se(coef) se2 Chisq DF p
sex -0.50984 0.169534 0.169498 9.04 1.00 0.0026
ph.karno -0.00620 0.006850 0.006842 0.82 1.00 0.3700
pat.karno -0.01703 0.006544 0.006541 6.77 1.00 0.0093
frailty(inst)0.25 0.21 0.3600
exp(coef) exp(-coef) lower .95 upper .95
sex 0.6006 1.6650.43080.8373
ph.karno 0.9938 1.0060.98061.0073
pat.karno0.9831 1.0170.97060.9958
Iterations: 5 outer, 21 Newton-Raphson
Variance of random effect= 0.001494078 I-likelihood = -711.9
Degrees of freedom for terms= 1.0 1.0 1.0 0.2
Concordance= 0.649 (se = 0.026 )
Likelihood ratio test= 23.22 on 3.21 df, p=4.703e-05
Le jeu. 15 août 2019 à 18:11, Denise b a écrit :
> Dear R help list,
>
> I would like to perform a score test for a subset of the parameters
> estimated with coxph using the frailty() option.
> As illustrated in the following reproducible example, I am able to perform
> the score test with the standard coxph() but not in presence of frailty()
> argument. See Examples 1 and 2
>
> library(survival)
>
> # Example 1: score test using standard coxph()
> # score test for ph.karno & pat.karno coefficients
> fit1<-coxph(Surv(time,status)~sex+ ph.karno + pat.karno , data=lung)
> fit0<-coxph(Surv(time,status)~pat.karno , data=lung) # fit under null
> scorestat <- coxph(Surv(time,status)~sex+ ph.karno + pat.karno,
> init=c(0,0,coefficients(fit0)) ,iter=0, data=lung)
> scorestat$score
> [1] 10.50409
>
> # Example 2: coxph with frailty term
> # score test for ph.karno & pat.karno coefficients
> frailtyfit1<-coxph(Surv(time,status)~sex+ ph.karno + pat.karno
> +frailty(inst), data=lung)
> frailtyfit0<-coxph(Surv(time,status)~sex+ frailty(inst), data=lung)
> scorefrailtyfit1<-coxph(Surv(time,status)~sex+ ph.karno + pat.karno
> +frailty(inst), data=lung,
> init=c(0,0,coefficients(frailtyfit0)) ,iter=0)
>
> I get the following error message:
> Error in coxph(Surv(time, status) ~ sex + ph.karno + pat.karno +
> frailty(inst), :
> wrong length for init argument
>
>
> When I checked in the code I found
> if (length(init) != ncol(X))
> stop("wrong length for init argument")
> Which seems that the length of the init argument doesnt match the number of
> predictors in the model.
>
> It seems that the code is expecting that I specify another parameter for
> frailty term.
> So, I tried to specify 2 as an inital value for frailty but got another
> error message later:
>
> scorefrailtyfit1<-coxph(Surv(time,status)~sex+ ph.karno + pat.karno
> +frailty(inst), data=lung,
> + init=c(0,0,coefficients(frailtyfit0),4.5),iter=0)
> Error in coxpenal.fit(X, Y, strats, offset, init = init, control, weights =
> weights, :
> Wrong length for inital values
>
> I have also checked coxme() but it seems that there is no option
>
