12
$\begingroup$

In a recent paper Norton et al. (2018)$^{[1]}$ state that

Different odds ratios from the same study cannot be compared when the statistical models that result in odds ratio estimates have different explanatory variables because each model has a different arbitrary scaling factor. Nor can the magnitude of the odds ratio from one study be compared with the magnitude of the odds ratio from another study, because different samples and different model specifications will have different arbitrary scaling factors. A further implication is that the magnitudes of odds ratios of a given association in multiple studies cannot be synthesized in a meta-analysis.

A small simulation illustrates this (R code is at the bottom of the question). Suppose the true model is: $$ \mathrm{logit}(y_{i})=1 + \log(2)x_{1i} + \log(2.5)x_{2i} + \log(3)x_{3i} + 0x_{4i} $$ Imagine further that the same data generated by the above model is analyzed by four different researchers using a logistic regression. Researcher 1 only includes $x_{1}$ as a covariate, researcher 2 includes both $x_{1}$ and $x_{2}$ and so forth. The average simulated estimates of the odds ratio for $x_{1}$ of the four researchers were:

res_1    res_2    res_3    res_4 
1.679768 1.776200 2.002157 2.004077

It's apparent that only researchers 3 and 4 get the correct odds ratio of about $2$ whereas researchers 1 and 2 do not. This does not happen in linear regression, which can be easily shown by a similar simulation (not shown here). I must confess that this result was quite surprising to me, although this problem seems to be well known$^{[2]}$. Hernán et al. (2011)$^{[3]}$ call this a "mathematical oddity" instead of a bias.

My questions:

  1. If odds ratios are basically uncomparable across studies and models, how can we combine the results of different studies for binary outcomes?
  2. What can be said about the countless meta-analyses that did combine the odds ratios from different studies where each study possibly adjusted for a different set of covariates? Are they essentially useless?

References

[1]: Norton EC, Dowd BE, Maciejewski ML (2018): Odds Ratios - Current Best Practice and Use. JAMA 320(1): 84-85.

[2]: Norton EC, Dowd BE (2017): Log Odds and the Interpretation of Logit Models. Health Serv Res. 53(2): 859-878.

[3]: Hernán MA, Clayton D, Keiding N (2011): The Simpson's paradox unraveled. Int J Epidemiol 40: 780-785.

Disclosure

The question (including the R code) is a modified version of a question posed by the user timdisher on datamethods.

R code

set.seed(142857)

n_sims <- 1000 # number of simulations

out <- data.frame(
  treat_1 = rep(NA, n_sims)
  , treat_2 = rep(NA, n_sims)
  , treat_3 = rep(NA, n_sims)
  , treat_4 = rep(NA, n_sims)
)

n <- 1000 # number of observations in each simulation

coef_sim <- "x1" # Coefficient of interest

# Coefficients (log-odds)

b0 <- 1
b1 <- log(2)
b2 <- log(2.5)
b3 <- log(3)
b4 <- 0

for(i in 1:n_sims){

  x1 <- rbinom(n, 1, 0.5)
  x2 <- rnorm(n)
  x3 <- rnorm(n) 
  x4 <- rnorm(n) 

  z <-  b0 + b1*x1 + b2*x2 + b3*x3 + b4*x4

  pr <- 1/(1 + exp(-z))  

  y <-  rbinom(n, 1, pr)

  df <-  data.frame(y = y, x1 = x1, x2 = x2, x3 = x3, x4 = x4)
  model1 <- glm(y ~ x1, data = df, family = "binomial")
  model2 <- glm(y ~ x1 + x2, data = df, family = "binomial")
  model3 <- glm(y ~ x1 + x2 + x3, data = df, family = "binomial")
  model4 <- glm(y ~ x1 + x2 + x3 + x4, data = df, family = "binomial")


  out$treat_1[i] <- model1$coefficients[coef_sim]
  out$treat_2[i] <- model2$coefficients[coef_sim]
  out$treat_3[i] <- model3$coefficients[coef_sim]
  out$treat_4[i] <- model4$coefficients[coef_sim]

}

# Coefficients

colMeans(out)
exp(colMeans(out)) # Odds ratios
$\endgroup$
11
+50
$\begingroup$

There are a number of alternative effects one can derive from the logistic regression model that do not suffer this same problem. One of the easiest is the average marginal effect of the variable. Assume the following logistic regression model:

\begin{equation} \ln\Bigg[\frac{p}{1-p}\Bigg]=X\beta + \gamma d \end{equation}

where $X$ is an $n$ (cases) by $k$ (covariates) matrix, $\beta$ are the regression weights for the $k$ covariates, $d$ is the treatment variable of interest and $\gamma$ is its effect.

The formula for the average marginal effect of $d$ would be:

\begin{equation} \frac{1}{n}\sum_{i=1}^n\Bigg[{\Big(1+e^{-(X\beta + \gamma)}\Big)^{-1} - \Big(1+e^{-X\beta}\Big)^{-1}}\Bigg] \end{equation}

This effect would be the average probability difference in the outcome between the treatment and control group for those who have the same values on other predictors (see Gelman & Hill, 2007, p. 101).

The corresponding R syntax given OP's example would be:

dydx_bin <- function(fit, coef) {
  mod.mat <- model.matrix(fit) # Obtain model matrix
  coefs <- coef(fit)
  oth_coefs <- coefs[!(names(coefs) == coef)] # Coefs bar focal predictor
  # Get model matrix excluding focal predictor
  X_nb <- as.matrix(mod.mat[, names(oth_coefs)])
  # Predictions for all data ignoring focal predictor
  Xb_nb <- X_nb %*% oth_coefs
  mean(plogis(Xb_nb + coefs[coef]) - plogis(Xb_nb))
}

I modified OP's syntax to demonstrate that it is not affected by which variables are in the model, as long as the predictor variable of interest is unrelated to other predictors.

I modified the results data frame thus:

out <- data.frame(
  treat_1 = rep(NA, n_sims), treat_2 = rep(NA, n_sims),
  treat_3 = rep(NA, n_sims), treat_4 = rep(NA, n_sims),
  treat_11 = rep(NA, n_sims), treat_21 = rep(NA, n_sims),
  treat_31 = rep(NA, n_sims), treat_41 = rep(NA, n_sims)
)

And within the simulation, I saved the computed average probability difference:

out$treat_11[i] <- dydx_bin(model1, coef_sim)
out$treat_21[i] <- dydx_bin(model2, coef_sim)
out$treat_31[i] <- dydx_bin(model3, coef_sim)
out$treat_41[i] <- dydx_bin(model4, coef_sim)

And the new results:

colMeans(out)[5:8]
 treat_11  treat_21  treat_31  treat_41 
0.1019574 0.1018248 0.1018544 0.1018642 

The estimated effect was consistent regardless of model specification. And adding covariates improved efficiency as with the linear regression model:

apply(out[, 5:8], 2, sd)
  treat_11   treat_21   treat_31   treat_41 
0.02896480 0.02722519 0.02492078 0.02493236 

There are additional effects that OP can compute like the average probability ratio between the two groups. The average probability difference computed above is available from the margins package in R and margins command in Stata. The average probability ratio is only available in Stata.

Onto the other question about trusting meta-analysis results. For one, the direction of the effect should not be useless. The problem with odds ratios does not affect the sign of the coefficients. So if a bulk of studies have an odds ratio above one, there is no reason to doubt this effect because of this particular problem.

As for the exact estimate, there is no reason to believe it. The nice thing is that if constituent studies are randomized controlled trials, then the odds ratios are conservative estimates and the actual results are even larger. This is because the effect OP demonstrated shrinks the odds ratios towards one. So if the bulk of studies have an odds ratio above 1 and the meta-analysis is pointing in this direction, then the actual OR once all relevant covariates are adjusted for is even larger. So these meta-analyses are not entirely useless.

But I would rather other effect estimates be used in meta-analysis. The average probability difference is one approach, and there are others.


Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press.

$\endgroup$
  • 1
    $\begingroup$ @COOLSerdash Thanks. There's so much more to explore here. It gets even more interesting when the odds ratio comes from a continuous variable that was dichotomized, especially if there was heteroskedasticity in the original relationships. See Achim Zeileis answer to this question - stats.stackexchange.com/questions/370876/… $\endgroup$ – Heteroskedastic Jim Oct 17 '18 at 14:01
  • $\begingroup$ Thanks for the link. I must confess that the derivation of the logistic model using an underlying latent continuous variable is new to me. I'm coming from biostatistics and the seminal sources in this field do not seem to mention these problems (e.g. the book by Lemeshow & Hosmer "Applied logistic regression"). I will award you the bounty as soon as I can (tomorrow). $\endgroup$ – COOLSerdash Oct 17 '18 at 14:14
  • $\begingroup$ I think derivation under underlying continuous variable is strange if you assume logistic errors. If you assume normal errors, it's more justifiable thanks to CLT. So for probit regression used a lot in econometrics, it is a common derivation. But if you will be dichotomizing a continuous variable, then derivation under errors is very helpful. Moreover, this derivation allows one to better explore the model generally and discover certain quirks. And thanks for the retrospective bounty. $\endgroup$ – Heteroskedastic Jim Oct 17 '18 at 14:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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