Question: cox proportional hazard model
gravatar for liu4gre
20 months ago by
United States
liu4gre200 wrote:

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


enter image description here

survival cox • 2.3k views
ADD COMMENTlink modified 20 months ago by Kevin Blighe50k • written 20 months ago by liu4gre200
gravatar for Kevin Blighe
20 months ago by
Kevin Blighe50k
Kevin Blighe50k wrote:

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:

    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"))
  [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)
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:

                  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,
  .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 COMMENTlink modified 7 months ago • written 20 months ago by Kevin Blighe50k

Thanks, Kevin. It looks cool. But it is for univariable test? How about multivariable, e.g. BRCA status, Tumor stage and Residual tumor?

ADD REPLYlink written 20 months ago by liu4gre200

If you want to further subdivide your cohort, then you could just create new categories, like, for example:


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 REPLYlink modified 7 months ago • written 20 months ago by Kevin Blighe50k

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 REPLYlink written 20 months ago by liu4gre200

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 REPLYlink written 20 months ago by Kevin Blighe50k

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 REPLYlink modified 7 months ago • written 20 months ago by Kevin Blighe50k

Great, thanks.

Sorry for one more question. df.OS you have put into the models are: BRCA1mutation BRCA2mutation BRCA1methylation,


BRCA1mutation.StageII BRCA1mutation.StageIII and IV BRCA2mutation.StageII BRCA2mutation.StageIII and IV

ADD REPLYlink written 20 months ago by liu4gre200

For the first OS model, the data would be:

OS    Event  BRCAstatus
1065  0       BRCA1mutation
0     0       BRCA1mutation
883   0       BRCA1methylation
33    1       BRCAwildtype

BRCAwildtype, BRCA1mutation, BRCA2mutation, BRCA1methylation

For PFS, it would be:

    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.

ADD REPLYlink modified 7 months ago • written 20 months ago by Kevin Blighe50k

Thanks. Nice tutorial.

ADD REPLYlink written 20 months ago by liu4gre200
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2639 users visited in the last hour