45
$\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?

$\endgroup$
37
$\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...

| cite | improve this answer | |
$\endgroup$
  • 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 '11 at 15:19
  • $\begingroup$ Does anyone know which one is implemented in Stata? $\endgroup$ – altroware Jul 20 at 15:30
44
$\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 ***
| cite | improve this answer | |
$\endgroup$
  • 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 '11 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 '11 at 17:04
  • 4
    $\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 '15 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 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.

| cite | improve this answer | |
$\endgroup$
5
$\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
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Is there any way to combine logistic display with a latex wrapper like outreg or xtable? $\endgroup$ – user226134 Nov 12 '19 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
```
| cite | improve this answer | |
$\endgroup$
  • $\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 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 at 15:50

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.