# Logistic Regression in R (Odds Ratio)

I'm trying to undertake a logistic regression analysis in R. I have attended courses covering this material using STATA. I am finding it very difficult to replicate functionality in R. Is it mature in this area? There seems to be little documentation or guidance available. Producing odds ratio output seems to require installing epicalc and/or epitools and/or others, none of which I can get to work, are outdated or lack documentation. I've used glm to do the logistic regression. Any suggestions would be welcome.

I'd better make this a real question. How do I run a logistic regression and produce odds rations in R?

Here's what I've done for a univariate analysis:

x = glm(Outcome ~ Age, family=binomial(link="logit"))

And for multivariate:

y = glm(Outcome ~ Age + B + C, family=binomial(link="logit"))

I've then looked at x, y, summary(x) and summary(y).

Is x$coefficients of any value? • Multivariate = multiple dependent variables. What you are describing is multivarible. Oct 2, 2021 at 12:08 ## 7 Answers if you want to interpret the estimated effects as relative odds ratios, just do exp(coef(x)) (gives you$e^\beta$, the multiplicative change in the odds ratio for$y=1$if the covariate associated with$\beta$increases by 1). For profile likelihood intervals for this quantity, you can do require(MASS) exp(cbind(coef(x), confint(x)))  EDIT: @caracal was quicker... • +1 for @fabian's suggestion. An inferior way to do this that usually yields similar intervals is to compute the interval on the logit scale and then transform to the odds scale: cbind( exp(coef(x)), exp(summary(x)$coefficients[,1] - 1.96*summary(x)$coefficients[,2]), exp(summary(x)$coefficients[,1] + 1.96*summary(x)$coefficients[,2]) ). There's also the delta method: ats.ucla.edu/stat/r/faq/deltamethod.htm Mar 23, 2011 at 15:19 • Does anyone know which one is implemented in Stata? Jul 20, 2020 at 15:30 • I've been using this: cbind( coef(summary(model))[, 4], exp(coef(model)), exp(confint(model))), because I also want a p-value, but it takes a long time, "Waiting for profiling to be done..." any suggestions as to how to fix this? – Mark Apr 13, 2023 at 18:59 You are right that R's output usually contains only essential information, and more needs to be calculated separately. N <- 100 # generate some data X1 <- rnorm(N, 175, 7) X2 <- rnorm(N, 30, 8) X3 <- abs(rnorm(N, 60, 30)) Y <- 0.5*X1 - 0.3*X2 - 0.4*X3 + 10 + rnorm(N, 0, 12) # dichotomize Y and do logistic regression Yfac <- cut(Y, breaks=c(-Inf, median(Y), Inf), labels=c("lo", "hi")) glmFit <- glm(Yfac ~ X1 + X2 + X3, family=binomial(link="logit"))  coefficients() gives you the estimated regression parameters$b_{j}$. It's easier to interpret$exp(b_{j})$though (except for the intercept). > exp(coefficients(glmFit)) (Intercept) X1 X2 X3 5.811655e-06 1.098665e+00 9.511785e-01 9.528930e-01  To get the odds ratio, we need the classification cross-table of the original dichotomous DV and the predicted classification according to some probability threshold that needs to be chosen first. You can also see function ClassLog() in package QuantPsyc (as chl mentioned in a related question). # predicted probabilities or: predict(glmFit, type="response") > Yhat <- fitted(glmFit) > thresh <- 0.5 # threshold for dichotomizing according to predicted probability > YhatFac <- cut(Yhat, breaks=c(-Inf, thresh, Inf), labels=c("lo", "hi")) > cTab <- table(Yfac, YhatFac) # contingency table > addmargins(cTab) # marginal sums YhatFac Yfac lo hi Sum lo 41 9 50 hi 14 36 50 Sum 55 45 100 > sum(diag(cTab)) / sum(cTab) # percentage correct for training data [1] 0.77  For the odds ratio, you can either use package vcd or do the calculation manually. > library(vcd) # for oddsratio() > (OR <- oddsratio(cTab, log=FALSE)) # odds ratio [1] 11.71429 > (cTab[1, 1] / cTab[1, 2]) / (cTab[2, 1] / cTab[2, 2]) [1] 11.71429 > summary(glmFit) # test for regression parameters ... # test for the full model against the 0-model > glm0 <- glm(Yfac ~ 1, family=binomial(link="logit")) > anova(glm0, glmFit, test="Chisq") Analysis of Deviance Table Model 1: Yfac ~ 1 Model 2: Yfac ~ X1 + X2 + X3 Resid. Df Resid. Dev Df Deviance P(>|Chi|) 1 99 138.63 2 96 110.58 3 28.045 3.554e-06 ***  • Thanks -- I'll need to look through your answer carefully. In STATA one can just run logit and logistic and get odds ratios and confidence intervals easily. I am somewhat frustrated that this appears to be so complicated and non-standard in R. Can I just use exp(cbind(coef(x), confint(x))) from fabians' answer below to get the OD and CI? I'm not clear what your answer is providing? Mar 23, 2011 at 13:28 • @SabreWolfy I wasn't sure what OR you are referring to: originally, I thought you meant the OR from the classification table that compares actual category membership with predicted membership (the cTab part in my answer). But now I see you probably just mean the$exp(b_{j})$: like fabians explained,$exp(b_{j})$equals the factor by which the predicted odds change when$X_{j}\$ increases by 1 unit. This is the ratio of the odds "after the 1-unit-increase" to "before the 1-unit-increase". So yes, you can just use fabians answer. Mar 23, 2011 at 17:04
• actually @SabreWolfy I find it frustrating that people can click a single button in stata/sas/spss etc, and obtain odds ratios (insert fit statistics, type III SS, whatever you like here) without having a clue as to what it means/how to calculate it/whether it is meaningful in a particular situation/and (perhaps more importantly) without having a working knowledge of the language itself.
– rawr
Jan 9, 2015 at 17:55
• The thing is that for a given data set there exists several different ways to compute OR’s. Consider a dichotomous variable (0,1 values) and a categorical variable with 3 categories. Depending on the reference coding several OR’s can be computed. Thus, care must be taken in setting appropriate reference levels for both variables in order to obtain the desired OR’s from the glm function. This can be achieved if the user knows how glm works for the case of the binomial family and the meaning of the coefficients for the (dummy, reference encoded) categorical variable used as covariate. Jul 24, 2020 at 15:51

The UCLA stats page has a nice walk-through of performing logistic regression in R. It includes a brief section on calculating odds ratios.

The epiDisplay package does this very easily.

library(epiDisplay)
data(Wells, package="carData")
glm1 <- glm(switch~arsenic+distance+education+association,
family=binomial, data=Wells)
logistic.display(glm1)
Logistic regression predicting switch : yes vs no

crude OR(95%CI)         adj. OR(95%CI)         P(Wald's test) P(LR-test)
arsenic (cont. var.)   1.461 (1.355,1.576)     1.595 (1.47,1.731)     < 0.001        < 0.001

distance (cont. var.)  0.9938 (0.9919,0.9957)  0.9911 (0.989,0.9931)  < 0.001        < 0.001

education (cont. var.) 1.04 (1.021,1.059)      1.043 (1.024,1.063)    < 0.001        < 0.001

association: yes vs no 0.863 (0.746,0.999)     0.883 (0.759,1.027)    0.1063         0.1064

Log-likelihood = -1953.91299
No. of observations = 3020
AIC value = 3917.82598

• Is there any way to combine logistic display with a latex wrapper like outreg or xtable?
– user226134
Nov 12, 2019 at 14:56

R has been mature with regard to odds ratio calculations more more than two decades. It's best to think about this in general terms. For example what if x1 and x2 are continuous and have nonlinear effects and interact with each other? Here is example code where the inter-quartile-range effect of x1 is computed, adjusted to x2=1.5.

require(rms)
# restricted cubic spline in x1, quadratic in x2, interactions
f <- lrm(y ~ rcs(x1, 4) * pol(x2, 2), data=mydata)
summary(f, x2=1.5)  # gives IQR odds ratio for x1
summary(f, x1=c(1, 3), x2=1.5)  # OR for x1=1 vs. x1=3, still x2=1.5


You can see that dealing with individual coefficients is not the general solution.

I tried @fabians's answer. It gave different results compared to @lockedoff's and @Edward answer when using a binary predictor. Please be careful when choosing the method.

For my own model, using @fabian's method, it gave Odds ratio 4.01 with confidence interval [1.183976, 25.038871] while @lockedoff's answer gave odds ratio 4.01 with confidence interval [0.94,17.05].

My model summary is as the following:

Call:
glm(formula = slicc_leuk ~ binf, family = binomial, data = kk)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.6823   0.3482   0.4625   0.4625   0.4625

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   2.1812     0.1783  12.232   <2e-16 ***
binf1         0.5914     0.6213   0.952   0.3412
binf2         1.3883     0.7388   1.879   0.0602 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 273.43  on 468  degrees of freedom
Residual deviance: 267.65  on 466  degrees of freedom
AIC: 273.65

Number of Fisher Scoring iterations: 6

> exp(cbind(coef(fit1), confint(fit1)))
Waiting for profiling to be done...
2.5 %    97.5 %
(Intercept) 8.857143 6.340297 12.782462
binf1       1.806452 0.618633  7.697776
binf2       4.008065 1.183976 25.038871

> logistic.display(fit1)

Logistic regression predicting slicc_leuk : 1 vs 0

OR(95%CI)          P(Wald's test) P(LR-test)
binf: ref.=0                                   0.056
1         1.81 (0.53,6.1)    0.341
2         4.01 (0.94,17.05)  0.06

Log-likelihood = -133.8248
No. of observations = 469
AIC value = 273.6496
$$$$

• I second this. It seems there are different methods (approximations) to get the confidence intervals. I guess the question is: Which method does Stata implement? Jul 20, 2020 at 15:23
• It seems exp(cbind(coef(glmFit), confint.default(glmFit)))  is the same as epiDisplay::logistic.display(glmFit, decimal=6) Jul 20, 2020 at 15:50
• @altroware based on the stata manual stata.com/manuals/rlogit.pdf, page 3, "Standard errors and confidence intervals are similarly transformed," which suggests to me that the method there is also to apply exp(.) to everything Nov 22, 2020 at 17:07

Similar to the choosen answer, but there is a direct command to get the exp(coefficients) and the intervals in one line.

require(MASS)

summ(my_model, exp = T)
`