Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

In R, is there a predefined function that will give me the log hazard ratio and its standard error for a black male (as shown in the example below) given the output of coxph regression?


#Kidney transplant data from Klein and Moeshberger. Massage data to make
#results look like those in book
data2 <- kidtran
data2$Gender <- "male"
    data2[data2$gender==2,7] <- "female"
data2$Race <- "white"
    data2[data2$race==2,8] <- "black"
data2$Gender <- as.factor(data2$Gender)
data2$Race <- as.factor(data2$Race)
data2$Race <- relevel(data2$Race,ref="white")

fit2 <- coxph(Surv(time,delta) ~ Gender * Race, data=data2)

#Relative log risk for a black male (reference white female) from
#page 252 in Klein and Moeshberger
(coef(fit2)[3] + coef(fit2)[2] + coef(fit2)[1])
sqrt(sum(diag(fit2$var)) + 2*fit2$var[2,1] + 2*fit2$var[3,1] + 2*fit2$var[3,2])

#Let's use predict
black.male <- data.frame(

white.female <- data.frame(

bm <- predict(fit2,newdata=black.male,se.fit=TRUE)

#bm in terms of original coefficients
coef(fit2)[3]*(1-fit2$means[3]) + coef(fit2)[2]*(1-fit2$means[2]) + coef(fit2)[1]*(1-fit2$means[1])

wf <- predict(fit2,newdata=white.female,se.fit=TRUE)

#Relative log risk and se for a black male (reference white female) 
bm$fit - wf$fit
sqrt(bm$se.fit^2 + wf$se.fit^2)
share|improve this question
add comment

1 Answer

up vote 2 down vote accepted

I do not think there exists a function like logHazardRatio("black males", "white females") which would return an estimate and a standard error. However, you can easily recover them as below.

Answer prior to your edit

Your Cox model is as follows, $$h_j(t) = h_0(t) \exp(\beta_1 \textrm{gender}_j + \beta_2 \textrm{race}_j + \beta_3 \textrm{gender}_j * \textrm{race}_j),$$ with $$\textrm{gender} = \left\{ \begin{array}{l} 0 \textrm{ if female} \\ 1 \textrm{ if male,} \\ \end{array} \right. $$ $$ \textrm{race} = \left\{ \begin{array}{l} 0 \textrm{ if white} \\ 1 \textrm{ if black,} \\ \end{array} \right.$$ and we have $$\frac{\textrm{hazard in black males}}{\textrm{hazard in white females}} = \frac{h_0(t) \exp(\beta_1 + \beta_2 + \beta_3)}{h_0(t)}$$ so that $$\log \left( \frac{\textrm{hazard in black males}}{\textrm{hazard in white females}} \right) = \beta_1 + \beta_2 + \beta_3.$$

Your estimate coef(fit2)[3] + coef(fit2)[2] + coef(fit2)[1] is thus correct :-)

To get a standard error, note that $$ \begin{align*} \textrm{Var}(\hat{\beta}_1 + \hat{\beta}_2 + \hat{\beta}_3) & = \textrm{Var}(\hat{\beta}_1) + \textrm{Var}(\hat{\beta}_2) + \textrm{Var}(\hat{\beta}_3) \\ & + 2 \textrm{Cov}(\hat{\beta}_1, \hat{\beta}_2) + 2 \textrm{Cov}(\hat{\beta}_1, \hat{\beta}_3) + 2 \textrm{Cov}(\hat{\beta}_2, \hat{\beta}_3). \end{align*} $$ Estimates for all of these terms can be found in the estimated covariance matrix of $(\hat{\beta}_1, \hat{\beta}_2, \hat{\beta}_3)'$:

> fit2$var
            [,1]        [,2]        [,3]
[1,]  0.03941978  0.02576211 -0.03945987
[2,]  0.02576211  0.09731841 -0.09729282
[3,] -0.03945987 -0.09729282  0.18242564

For example, $\widehat{\textrm{Cov}}(\hat{\beta}_1, \hat{\beta}_3) = -0.03945987.$

share|improve this answer
Yes, but is there a function in R that does all of this for you? I thought that predict might do it, but I have been unable to get it to give me answers like the ones you provided. –  Tom Nov 23 '12 at 20:36
I do not think there exists a function like 'logHazardRatio("black males", "white females")' which would return an estimate and a standard error. But it is not too tedious to go through calculations like in my answer... –  ocram Nov 23 '12 at 20:38
I could not figure out how to use predict to do this calculation. Apparently predict gives the risk relative to a person who has mean values for all covariates, so perhaps I can just calculate the log hazard ratio for black male and subtract out the log hazard ratio for white female. This seems to work, but the SE is different than what you calculated above. –  Tom Nov 24 '12 at 18:26
@Tom: 1) Can you say what 'bm' (in your edited question) is supposed to represent in terms of beta coefficients? 2) You computed Var('bm' - 'wf') assuming 'bm' and 'wf' being uncorrelated... which is probably false because your result differs from what we know to be true... –  ocram Nov 24 '12 at 19:36
I added what bm is in terms of the original coefficients. I don't understand completely why R centers on the mean values of the factors for predicted values. Yes, it is unfortunate that the SE is wrong because otherwise this approach is giving me what I wanted. –  Tom Nov 24 '12 at 21:42
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.