cox proportional hazard model
1
4
Entering edit mode
4.4 years ago
liu4gre ▴ 210

Hi All,

I am wondering how to derive HR, CI and p values for each factor from cox model like follows. Using coxph only gives these values for groups such as BRCA status, TUmor stage ...

THanks.

survival cox • 7.2k views
8
Entering edit mode
4.3 years ago

Update 13th March, 2019:

I posted a related tutorial: Survival analysis with gene expression

## ------------------------------

Update 24th September, 2018:

Note that the actual plot does not match the data, and neither do the stat values. Everything here is purely for display purposes only.

## -------------------------------------

This would have been performed in the realm of survival analysis, looking at overall survival (OS) and progression-free survival (PFS), as you can probably see.

The starting point for the Cox Proportional Hazards Regression (Cox) is data in this format:

head(df)
OS Event  Group
1 1065      0 group1
2    0      0 group2
3  883      0 group1
4   33      1 group2
5  790      0 group1
6 2517      1 group2


The columns are

• OS: overall survival (days, weeks, months, years - just needs to be consistent)
• Event: e.g. death, diagnosis, or some other event
• Group: the categories of interest - can be anything such as ER status, IHC scores for CD20, race, or something else

Cox is run with coxph in R, and it needs to be performed on a survival object, e..g, produced by Surv

As per the table (above), there is a reference level for the category of interest, e.g., BRCA wild-type. Thus, we must also choose a reference category against which all other categories will be compared (here group1 is the reference):

df$Group <- factor(df$Group, levels=c("group1","group2","group3","group4"))
df$Group [1] group1 group2 group1 group2 group1 group2 group3 group2 group3 group2 [11] group1 group4 group4 group3 group4 group4 group2 group3 group1 group3 [21] group4 group4 group4 group1 group3 group3 group2 group1 group3 group4 [31] group1 group1 group4 group2 group3 group3 group4 group3 group2 group4 [41] group4 group3 group3 group4 group4 group4 group3 group2 group2 group1 *et cetera* Levels: group1 group2 group3 group4  Now we can actually generate hazard ratios (including CIs) and P values: coxmodel <- coxph(Surv(time = OS, event = Event) ~ Group, data=df) summary(coxmodel) Call: coxph(formula = Surv(time = OS, event = Event) ~ Group, data = df) n= 106, number of events= 106 coef exp(coef) se(coef) z Pr(>|z|) Groupgroup2 0.15929 1.17267 0.29957 0.532 0.595 Groupgroup3 0.03724 1.03794 0.27747 0.134 0.893 Groupgroup4 -0.14772 0.86267 0.28570 -0.517 0.605 exp(coef) exp(-coef) lower .95 upper .95 Groupgroup2 1.1727 0.8528 0.6519 2.109 Groupgroup3 1.0379 0.9634 0.6025 1.788 Groupgroup4 0.8627 1.1592 0.4928 1.510 Concordance= 0.515 (se = 0.032 ) Rsquare= 0.011 (max possible= 0.999 ) Likelihood ratio test= 1.19 on 3 df, p=0.7566 Wald test = 1.18 on 3 df, p=0.7575 Score (logrank) test = 1.19 on 3 df, p=0.7563  The P values for each category are given by Pr(>|z|). The HRs are given by exp(coef). and you can probably guess the CIs. Just to be sure, here are the HRs with 2.5% and 97.5% CIs: exp(confint(coxmodel)) 2.5 % 97.5 % Groupgroup2 0.6519067 2.109444 Groupgroup3 0.6025433 1.787944 Groupgroup4 0.4927892 1.510188  ## ---------------------- Finally, you can then actually plot the Kaplan-Meier survival curve for this using a wrapper, km.coxph.plot: km.coxph.plot(formula.s=Surv(time=OS, event = Event) ~ Group, data.s=df, mark.time=TRUE, x.label="Time (days)", y.label="Overall survival", main.title="", leg.text=c("Group1","Group2","Group3", "Group4"), leg.pos="topright", leg.bty="n", leg.inset=0, .col=c("limegreen","royalblue","purple","red1"), o.text="", .lty=c(1,1,1,1), .lwd=c(1.75,1.75,1.75,1.75), show.n.risk=TRUE, n.risk.step=500, n.risk.cex=0.8, verbose=FALSE) mtext(side=3, line=-1, adj=-0.25, "Cox PH survival", cex=3) mtext(side=3, line=-13, adj=0.95, "HR=2.95 (0.52, 16.62), p=0.2", cex=0.8, col="red")  ADD COMMENT 0 Entering edit mode Thanks, Kevin. It looks cool. But it is for univariable test? How about multivariable, e.g. BRCA status, Tumor stage and Residual tumor? ADD REPLY 1 Entering edit mode If you want to further subdivide your cohort, then you could just create new categories, like, for example: BRCA1mutation.StageI BRCA1mutation.StageII BRCA1mutation.StageIII BRCA1mutation.StageIV BRCA2mutation.StageI BRCA2mutation.StageII BRCA2mutation.StageIII BRCA2mutation.StageIV  You can also adjust for other factors / covariates in the Cox model, or do interactions, which is an alternative hypothesis to the above but still useful: coxph(Surv(OS, Event) ~ Group + TumourStage, data=df) coxph(Surv(OS, Event) ~ Group:TumourStage, data=df)  I recall having a conversation in this regard last year amongst a group of statisticians (i.e., interaction terms in a Cox model - from what I recall, it's a somewhat unexplored area). ADD REPLY 0 Entering edit mode Thanks. But it looks the reference for each category can not be defined? E.g. Stage I as reference in Tumor stage and BRCA wildtype in BRCA status? Or in above table, the reference was defined individually, which means three testes were performed to generate the table of PFS or OS? ADD REPLY 1 Entering edit mode In that table that you posted, the reference levels are 'BRCA wild-type', 'Tumour Stage II', and 'Residual Tumour 0'. They did not have data for Stage I or they just did not consider it for the study. They would have used 6 tests (3 for Os; 3 for PFS) to generate the results in that table. ADD REPLY 1 Entering edit mode These would have been the models: coxph(Surv(OS, Event) ~ BRCAstatus, data=df.OS) coxph(Surv(OS, Event) ~ TumourStage, data=df.OS) coxph(Surv(OS, Event) ~ ResidualTumour, data=df.Os) coxph(Surv(PFS, Event) ~ BRCAstatus, data=df.PFS) coxph(Surv(PFS, Event) ~ TumourStage, data=df.PFS) coxph(Surv(PFS, Event) ~ ResidualTumour, data=df.PFS)  ADD REPLY 0 Entering edit mode Great, thanks. Sorry for one more question. df.OS you have put into the models are: BRCA1mutation BRCA2mutation BRCA1methylation, or: BRCA1mutation.StageII BRCA1mutation.StageIII and IV BRCA2mutation.StageII BRCA2mutation.StageIII and IV ADD REPLY 1 Entering edit mode For the first OS model, the data would be: head(df.OS) OS Event BRCAstatus 1065 0 BRCA1mutation 0 0 BRCA1mutation 883 0 BRCA1methylation 33 1 BRCAwildtype ... levels(df.OS$BRCAstatus)
BRCAwildtype, BRCA1mutation, BRCA2mutation, BRCA1methylation


For PFS, it would be:

head(df.PFS)
PFS  Event  BRCAstatus
501  0       BRCA1mutation
0    0       BRCA1mutation
38   0       BRCA1methylation
10   1       BRCAwildtype
...


Then, there are 4 more different models:

1. OS with Tumour Stage
2. OS with Residual Tumour
3. PFS with Tumour Stage
4. PFS with Residual Tumour

Hope that this helps.

1
Entering edit mode

Thanks. Nice tutorial.