56
$\begingroup$

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?

Improve this question
$\endgroup$
1
  • 2
    $\begingroup$ Multivariate = multiple dependent variables. What you are describing is multivarible. $\endgroup$
    – Frank Harrell
    Oct 2, 2021 at 12:08

7 Answers 7

Reset to default
47
$\begingroup$

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...

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ +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 $\endgroup$
    – lockedoff
    Mar 23, 2011 at 15:19
  • $\begingroup$ Does anyone know which one is implemented in Stata? $\endgroup$
    – altroware
    Jul 20, 2020 at 15:30
  • $\begingroup$ 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? $\endgroup$
    – Mark
    Apr 13 at 18:59
46
$\begingroup$

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 ***
Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ 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? $\endgroup$
    – SabreWolfy
    Mar 23, 2011 at 13:28
  • 3
    $\begingroup$ @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. $\endgroup$
    – caracal
    Mar 23, 2011 at 17:04
  • 5
    $\begingroup$ 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. $\endgroup$
    – rawr
    Jan 9, 2015 at 17:55
  • $\begingroup$ 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. $\endgroup$
    – Nicolas Molano
    Jul 24, 2020 at 15:51
23
$\begingroup$

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

Improve this answer
$\endgroup$
2
7
$\begingroup$

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
Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ Is there any way to combine logistic display with a latex wrapper like outreg or xtable? $\endgroup$
    – user226134
    Nov 12, 2019 at 14:56
1
$\begingroup$

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
```
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – altroware
    Jul 20, 2020 at 15:23
  • $\begingroup$ It seems exp(cbind(coef(glmFit), confint.default(glmFit))) is the same as epiDisplay::logistic.display(glmFit, decimal=6) $\endgroup$
    – altroware
    Jul 20, 2020 at 15:50
  • $\begingroup$ @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 $\endgroup$
    – Richard DiSalvo
    Nov 22, 2020 at 17:07
1
$\begingroup$

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)
dd <- datadist(mydata); options(datadist='dd')
# 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.

Improve this answer
$\endgroup$
0
0
$\begingroup$

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

Current choosen answer:

require(MASS)
exp(cbind(coef(x), confint(x)))

In one line:

summ(my_model, exp = T)

(I could not make a comment since I'm still a newbi and don't have enough reputation points on the website).

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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