## Load surival package
library(survival)
## Load Ovarian Cancer Survival Data
data(ovarian)
## Fit Cox regression model
fit <- coxph(Surv(futime, fustat) ~ age + ecog.ps,
data = ovarian)
## Test the Proportional Hazards Assumption of a Cox Regression
fitZph <- cox.zph(fit)
## Print test results
fitZph
rho chisq p
age -0.243 0.856 0.355
ecog.ps 0.520 2.545 0.111
GLOBAL NA 3.195 0.202
## Show plots
plot(fitZph)
See if the hazard ratio changes before and after 200 days.
## split follow-up interval at 200 days
ovarian.split <- survSplit(data = ovarian,
cut = 200,
end = "futime",
event = "fustat",
start = "start",
id = "id",
zero = 0,
episode = NULL)
## Create interval indicator and survival object
ovarian.split <- within(ovarian.split, {
## Interval indicator (or can just use factor(start))
interval.indicator <- factor(start, c(0,200), c("<=200",">200"))
## Surv(start, stop, event) format
SurvObj <- Surv(start, futime, fustat)
})
## Fit model with interaction. Use cluster(id) to account for clustering by subject
fit2 <- coxph(SurvObj ~ age + ecog.ps + (age + ecog.ps):interval.indicator + cluster(id),
data = ovarian.split)
## Show results: Both interaction terms are not significant, i.e., the effect is stable over time.
summary(fit2)
Call:
coxph(formula = SurvObj ~ age + ecog.ps + (age + ecog.ps):interval.indicator +
cluster(id), data = ovarian.split)
n= 49, number of events= 12
coef exp(coef) se(coef) robust se z Pr(>|z|)
age 0.2573 1.2934 0.1154 0.0507 5.07 0.0000004 ***
ecog.ps -1.3896 0.2492 1.4418 1.0538 -1.32 0.19
age:interval.indicator>200 -0.1364 0.8725 0.1301 0.0934 -1.46 0.14
ecog.ps:interval.indicator>200 1.8198 6.1707 1.6108 1.2528 1.45 0.15
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.293 0.773 1.1710 1.43
ecog.ps 0.249 4.013 0.0316 1.97
age:interval.indicator>200 0.872 1.146 0.7266 1.05
ecog.ps:interval.indicator>200 6.171 0.162 0.5296 71.89
Concordance= 0.803 (se = 0.091 )
Rsquare= 0.289 (max possible= 0.76 )
Likelihood ratio test= 16.7 on 4 df, p=0.00218
Wald test = 33.3 on 4 df, p=0.00000103
Score (logrank) test = 15.1 on 4 df, p=0.00455, Robust = 6.97 p=0.138
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).
## Fit model without the main effects
fit3 <- coxph(SurvObj ~ (age + ecog.ps):interval.indicator + cluster(id),
data = ovarian.split)
## Show results: Now the interaction terms represent interval-specific effect (not effect difference).
summary(fit3)
Call:
coxph(formula = SurvObj ~ (age + ecog.ps):interval.indicator +
cluster(id), data = ovarian.split)
n= 49, number of events= 12
coef exp(coef) se(coef) robust se z Pr(>|z|)
age:interval.indicator<=200 0.2573 1.2934 0.1154 0.0507 5.07 0.0000004 ***
age:interval.indicator>200 0.1209 1.1285 0.0600 0.0695 1.74 0.082 .
interval.indicator<=200:ecog.ps -1.3896 0.2492 1.4418 1.0538 -1.32 0.187
interval.indicator>200:ecog.ps 0.4302 1.5376 0.7182 0.6501 0.66 0.508
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age:interval.indicator<=200 1.293 0.773 1.1710 1.43
age:interval.indicator>200 1.128 0.886 0.9847 1.29
interval.indicator<=200:ecog.ps 0.249 4.013 0.0316 1.97
interval.indicator>200:ecog.ps 1.538 0.650 0.4300 5.50
Concordance= 0.803 (se = 0.091 )
Rsquare= 0.289 (max possible= 0.76 )
Likelihood ratio test= 16.7 on 4 df, p=0.00218
Wald test = 33.3 on 4 df, p=0.00000103
Score (logrank) test = 15.1 on 4 df, p=0.00455, Robust = 6.97 p=0.138
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).