Odds ratios vs hazard ratios

library("epitools")

How do we compare/combine the results of several studies, some of which measure odds ratios and some of which measure hazard ratios?

Odds are defined as $p/(1-p)$, where $p$ is the probability of success; we are often interested in log-odds, which are additive rather than multiplicative and can range over all real values: $\log(p/(1-p))$. The inverse transformation is logistic. In R the log-odds (or logit) transformation is qlogis(), the inverse-logit/logistic transformation is plogis(). A probability of 0.5 corresponds to even or “1:1” odds, or log-odds of $\log(1)=0$. (I’m using the natural logarithm throughout.) Odds are symmetric: if you want to measure or model probabilities of survival rather than mortality, you just take the reciprocal (or change the sign if you’re working with log-odds). Odds ratios are multiplicative differences between the odds in different treatments, or per unit of a continuous predictor variable.

par(las=1,bty="l")
curve(qlogis,from=0.01,to=0.99,
      xlab="probability",ylab="log-odds")
abline(h=0,lty=2)

In contrast, hazards are the probabilities of mortality (failure) per unit time. For a single time step, you can go from hazard $\beta$ to probability $p$ via $p=1-\exp(-\beta)$. The inverse-transform is the complementary log $\beta=-\log(1-p)$. Like going from odds to log-odds, it is often more convenient to work with log-hazards, which can be positive or negative (hazards must be $\geq 0$) and which combine additively; then we have the complementary log-log transformation, $\log(-\log(1-p))$. (Unlike the logit, the complementary log-log is asymmetric: $p=1-\exp(-1) \approx 0.63$ gives a hazard of 1 or a log-hazard of 0.) Hazard ratios, analogous to odds ratios, are the ratios of hazards under different conditions. Hazards have units of (1/time), but hazard ratios are unitless.

par(las=1,bty="l")
curve(log(-log(1-x)),from=0.01,to=0.9999,
      xlab="probability",ylab="log-hazard")
abline(h=0,lty=2)
abline(v=1-exp(-1),lty=2)

Odds and hazards represent similar transformations, but they are not quite the same. If you want to translate an effect (e.g. the difference in log-odds between a control and a treatment group), you have to translate the odds or hazards back to the original probability scale, then translate

One way to get the functions we need is to extract them from GLM family() objects, which contain link functions (translating from the (0,1) interval to the real line, e.g. the logit or log-odds functions) in their $linkfun slots and inverse link functions (translating the other way, e.g. the logistic or inverse-logit function) in their $linkinv slots.

ff.log <- binomial(link="log")
ff.logit <- binomial(link="logit")
ff.cloglog <- binomial(link="cloglog")

We could also define these ourselves, for more clarity (the internal functions above are written in C for efficiency, and may have some tweaks to make them more robust, e.g. not returning 0 or 1 for extreme values of the input).

## same as plogis() or ff.logit$linkfun:
logit <- function(x) log(x/(1-x))
## same as qlogis() or ff.logit$linkinv:
logistic <- function(x) 1/(1+exp(-x))
all.equal(logit(logistic(1.375)),1.375)

## [1] TRUE

## same as ff.cloglog$linkfun:
cloglog <- function(x) log(-log(1-x))
## same as ff.cloglog$linkinv
gompertz <- function(x) 1-exp(-exp(x))
all.equal(cloglog(gompertz(1.375)),1.375)

## [1] TRUE

So in order to convert from log-odds to log-hazards, you have to (1) convert the log-odds for both the baseline and treatment groups back to the probability scale (using the logistic/inverse-logit function) (2) convert both probabilities to the log-hazard scale (3) compute the difference between the log-hazards. Conversion of estimates this way is not too hard, but conversions of standard errors are a little bit trickier (I’ll try to do this example below).

convhaz <- function(base,trt.eff) {
    ## translate base and trt values from log-odds to probabilities
    L1 <- ff.logit$linkinv(c(base,base+trt.eff))
    L2 <- ff.cloglog$linkfun(L1)
    return(diff(L2))
}

Suppose the base log-odds is 1 and the difference in log-odds between control and treatment is 2:

convhaz(1,2)

## [1] 0.8421644

If the base log-odds is only 0 (the difference is still 2):

convhaz(0,2)

## [1] 1.121192

(To do: demonstration of confidence intervals, standard errors on log-hazard vs. hazard scale. Bottom line: if someone gives you 95% CI on hazard ratios, you should probably guess that they computed them as $\mu \pm 1.96 \sigma$ on the log-hazard scale, and invert (solve for $\sigma$) as follows: $\sigma = \left(\log(\textrm{upper})-\log(\textrm{lower})\right)/(2 \times 1.96)$.)

Standard errors

This is a little bit trickier. In principle we can use a version of the delta method (actually the Wikipedia link is to “Taylor expansions for the moments of functions of random variables”, but the general idea of linearization is the same) to translate the standard errors, but we have to use the full multivariate version. Suppose we have a multivariate, nonlinear function $(x_1',x_2') = F(x_1,x_2)$ (two-dimensional for notational simplicity, and because that’s all we need, but the stuff below is more general). Call the gradient matrix \[ G=\partial F(x)_j/\partial x_i. \] Then the variance-covariance matrix of the transformed vector $x'$ is

\[ V' = G^\top V G. \]

If we want to split the transformation $F$ up into several steps $F_1 F_2 \ldots$ then the chain rule says we can just compute $G$ as $G_1 G_2 \ldots$. Here we have four transformations (individually pretty simple): (1) add baseline and treatment diff; (2) inverse-logit transform; (3) cloglog-transform; (4) take the difference.

Suppose we know the standard errors of the log-odds for the baseline level and the treatment effect and the correlations between them (in general we don’t know this; we’ll probably have to pretend they’re independent, but it’s worth working out the general case first). Conveniently, the derivatives for the inverse-link functions are stored in a family object ff as ff$mu.eta, and the derivatives for the link functions are the reciprocals.

I’m not sure of the effects of assuming independence; in general I’m guessing it will inflate the standard errors slightly (assuming the intercept and treatment effect are generally negatively correlated?)

Relative risk/crude rates

What if we are given, not log-odds, but crude rates (probabilities) or relative risks (ratios of probabilities)? This is basically a subset of the previous problem — one way or the other we have to get back to the original effects on the probability scale (and their standard errors), and then go forward to the log-hazard scale. Can I adapt the function above to make it more general?

convhaz3 <- function(b1,b2,b1.se,b2.se,
                     b12.cor=NULL,
                     diffvals=TRUE,
                     link.in="logit",
                     link.out="cloglog",debug=FALSE) {
    ## step 1: {base, diff} to {ttt1, ttt2}
    ## step 2: log-odds (or whatever) to probs
    ## step 3: probs to cloglog
    ## step 4: {ttt1, ttt2} to {base, diff}
    ## for log-odds (b1=base, b2=diff) we need all 4
    ## if we start from probs and SE of probs we only need (3,4)
    ## translate base and trt values from log-odds to probabilities
    if (diffvals) {
        x1 <- c(b1,b1+b2)
        G1 <- matrix(c(1,0,1,1),nrow=2)
    } else {
        x1 <- c(b1,b2)
        G1 <- diag(2)
    }
    trans <- binomial(link=link.in)
    x2 <- trans$linkinv(x1)
    out <- binomial(link=link.out)
    x3 <- out$linkfun(x2)
    res <- diff(x3)
    if (debug) {
        print(rbind(input=c(b1,b2),
                    tttlev=x1,
                    trans=x2,
                    cloglog=x3,
                    diff=res))
    }

    V <- diag(c(b1.se^2,b2.se^2)) ## ASSUME independence
    if (!is.null(b12.cor)) {
        V[1,2] <- V[2,1] <- b12.cor*b1.se*b2.se
    }
    if (debug) print(V)
    ## convert standard errors
    ## gradient of (base,trt=base+trt diff) wrt (base,trt diff)
    ## change notation from above (opposite numbering order)
    ## gradient of (prob(base),prob(trt))
    ## mu.eta is the derivative of the inverse-link function
    ## 
    G2 <- diag(trans$mu.eta(x1))
    G3 <- diag(1/out$mu.eta(x3)) 
    G4 <- matrix(c(-1,1),ncol=1) ## gradient of diff
    G.all <- G1 %*% G2 %*% G3 %*% G4
    if (debug) { G.trans <- G2 %*% G3
                 print(t(G.trans) %*% G.trans)
                 G.diff <- G1 %*% G4
                 t(G.diff) %*% G.diff
                 print(G.all)
                }
    res.var <- t(G.all) %*% V %*% G.all
    return(c(val=res,SE=sqrt(res.var)))
}

Testing

Test on some data:

dd <- data.frame(dead=c(100,200),alive=c(300,400),ttt=c("control","treat"))
g0 <- glm(cbind(dead,alive)~ttt,data=dd,family=binomial(link="cloglog"))
g1 <- glm(cbind(dead,alive)~ttt,data=dd,family=binomial(link="logit"))
tmpf <- function(link="logit",...) {
    g <- update(g0,family=binomial(link=link))
    c1 <- coef(summary(g))[,c("Estimate","Std. Error")]
    c1 <- setNames(as.data.frame(c1),c("est","se"))
    convhaz3(c1[1,"est"],c1[2,"est"],c1[1,"se"],c1[2,"se"],
             link.in=link,...)
}
## direct log-hazards estimate
c1 <- setNames(coef(summary(g0))["ttttreat",c("Estimate","Std. Error")],
               c("est","se"))

Testing:

cc <- cov2cor(vcov(g1))[2,1]  ## logit cor
cc2 <- cov2cor(vcov(g0))[2,1]  ## cloglog cor
(resMat <- rbind(orig=c1,
                logit=tmpf(),
                logit.cor=tmpf(b12=cc),
                cloglog=tmpf("cloglog"),
                cloglog.cor=tmpf("cloglog",b12.cor=cc2),
                log=tmpf()))

##                   est        se
## orig        0.3431789 0.1230366
## logit       0.3431789 0.1187836
## logit.cor   0.3431789 0.1230367
## cloglog     0.3431789 0.1230366
## cloglog.cor 0.3431789 0.1230366
## log         0.3431789 0.1187836

This looks pretty good now. We get the estimates right; the SEs are close (and exact if we specify the correlations). The only puzzle here is why the cloglog case is right even without putting in the correlation … ? Perhaps because the transformation matrix is essentially just picking out the second component?