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?
 A: 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

A: 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 ***

A: 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...
A: The UCLA stats page has a nice walk-through of performing logistic regression in R. It includes a brief section on calculating odds ratios.
A: 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.
A: 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).
