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?

up vote 34 down vote accepted

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
    +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 – lockedoff Mar 23 '11 at 15:19

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 ***
  • 1
    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? – SabreWolfy Mar 23 '11 at 13:28
  • 3
    @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. – caracal Mar 23 '11 at 17:04
  • 3
    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 '15 at 17:55

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

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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