Logistic Regression in R: glm(), Odds Ratios, and Diagnostics

To fit a logistic regression in R, use glm(outcome ~ predictors, family = binomial, data = df), then convert the log-odds coefficients to odds ratios with exp(coef(model)) and attach profile-likelihood confidence intervals with exp(confint(model)). The model object also carries everything you need to check it: deviance for likelihood-ratio tests, predict(type = "response") for fitted probabilities, and convergence warnings that flag separation.

This page is the procedure: the exact calls, the actual console output, and the interpretation of every number in it. If you want the underlying method — what the logit link is, why odds ratios are not risk ratios, how the likelihood is maximised — read logistic regression and the logit model first and come back here to run it. The two pages are deliberately split: that one explains the model, this one explains R.

Every code block below was executed in R 4.6.1 and every block of output is the real console result, not a stylised approximation. Where R prints a warning or a message, it is reproduced verbatim, because those strings are what you will actually search for when something goes wrong.

The example data: MASS::birthwt

The examples use birthwt, which ships with the MASS package. MASS is a recommended package, meaning it is installed with every standard R distribution — nothing to download. The data frame has 189 rows and 10 columns, collected at Baystate Medical Center, Springfield, Massachusetts during 1986, and originates with Hosmer and Lemeshow’s Applied Logistic Regression (Wiley, 1989). The outcome low indicates a birth weight below 2.5 kg.

Recode the categorical columns as factors before fitting. R will happily treat a 1/2/3 race code as a continuous number and give you a meaningless single slope if you do not.

library(MASS)

bw <- birthwt
bw$low   <- factor(bw$low,   levels = c(0, 1), labels = c("normal", "low"))
bw$race  <- factor(bw$race,  levels = 1:3,     labels = c("white", "black", "other"))
bw$smoke <- factor(bw$smoke, levels = c(0, 1), labels = c("no", "yes"))
bw$ht    <- factor(bw$ht,    levels = c(0, 1), labels = c("no", "yes"))
bw$ui    <- factor(bw$ui,    levels = c(0, 1), labels = c("no", "yes"))
bw$ptd   <- factor(bw$ptl > 0, levels = c(FALSE, TRUE), labels = c("no", "yes"))

table(bw$low)

normal    low
   130     59

59 events in 189 observations. Hold on to that number — it constrains how many predictors the model can support, and we return to it below.

Fitting the model with glm()

fit <- glm(low ~ age + lwt + race + smoke + ptd + ht + ui,
           family = binomial(link = "logit"), data = bw)

Three things about the family argument that trip people up:

• family = binomial, family = binomial() and family = binomial(link = "logit") are all the same model. The signature is binomial(link = "logit"), so the logit link is the default and the long form is only worth writing when you want it to be explicit in a script someone else will read.
• family = "binomial" (a quoted string) also works, because glm() resolves a character family name. It is not an error, despite what some tutorials imply.
• Omitting family entirely gives you gaussian — ordinary linear regression on a 0/1 outcome, silently. There is no warning. This is the single most common way a “logistic regression” script turns out not to be one.

The outcome-coding trap: R models the second factor level

When the response is a factor, glm() models the probability of the second level, treating the first as the reference. Nothing in the output states which one that is. Get the level order backwards and every odds ratio on the page is its own reciprocal.

contrasts(bw$low)

       low
normal   0
low      1

Here normal is 0 and low is 1, so the model estimates the odds of a low birth weight. Reverse the level order and watch the effect of smoking invert:

fit_fwd <- glm(low ~ smoke, family = binomial, data = bw)
exp(coef(fit_fwd))["smokeyes"]          # outcome level = "low"
#> 2.022

bw$low_rev <- factor(as.character(bw$low), levels = c("low", "normal"))
fit_rev <- glm(low_rev ~ smoke, family = binomial, data = bw)
exp(coef(fit_rev))["smokeyes"]          # outcome level = "normal"
#> 0.495

2.022 and 0.495 are reciprocals of each other. Both models are correct; they answer opposite questions. Always print levels(df$outcome) before you interpret a coefficient, and state in your write-up which level was modelled. A plain 0/1 integer outcome avoids the ambiguity entirely and produces byte-identical coefficients — verified above with all.equal().

Reading summary()

summary(fit)

Call:
glm(formula = low ~ age + lwt + race + smoke + ptd + ht + ui,
    family = binomial(link = "logit"), data = bw)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.63691    1.23028    0.52   0.6047
age         -0.03775    0.03781   -1.00   0.3181
lwt         -0.01491    0.00704   -2.12   0.0342 *
raceblack    1.21274    0.53248    2.28   0.0228 *
raceother    0.80412    0.44843    1.79   0.0729 .
smokeyes     0.84640    0.40806    2.07   0.0381 *
ptdyes       1.22175    0.46301    2.64   0.0083 **
htyes        1.83869    0.70324    2.61   0.0089 **
uiyes        0.71113    0.46311    1.54   0.1247
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 234.67  on 188  degrees of freedom
Residual deviance: 196.83  on 180  degrees of freedom
AIC: 214.8

Number of Fisher Scoring iterations: 4

What each part actually means:

• Estimate is on the log-odds scale, not the probability scale and not the odds scale. smokeyes = 0.84640 is a log-odds difference; it becomes an odds ratio only after exponentiation.
• Factor level names are concatenated: raceblack is the black level of race, contrasted against white, the first level. raceother is likewise versus white. There is no coefficient for the reference level, by construction.
• z value, not t value. Logistic regression uses the Wald z statistic because the binomial dispersion is fixed at 1 — which is exactly what the “dispersion parameter taken to be 1” line is telling you. If that line ever says something other than 1, you are not looking at a binomial fit.
• Null vs residual deviance: 234.67 on 188 df for the intercept-only model, 196.83 on 180 df for this one. The difference, 37.84 on 8 df, is the global likelihood-ratio statistic. R does not test it for you in summary() — you have to ask (below).
• Fisher Scoring iterations: 4. A small number is healthy. Something like 25 is a red flag for separation.

Odds ratios: exponentiating the coefficients

round(exp(coef(fit)), 3)

(Intercept)         age         lwt   raceblack   raceother    smokeyes      ptdyes
      1.891       0.963       0.985       3.363       2.235       2.331       3.393
      htyes       uiyes
      6.288       2.036

Read these as multiplicative effects on the odds of the modelled outcome, holding the other predictors fixed. Smoking multiplies the odds of low birth weight by 2.33. A history of hypertension multiplies them by 6.29. lwt = 0.985 means each additional pound of maternal weight multiplies the odds by 0.985 — a 1.5% reduction per pound.

Rescale continuous predictors before you report them

An odds ratio of 0.985 per pound is technically correct and practically useless. Because the model is linear in the log-odds, you can rescale by multiplying the coefficient before exponentiating — you do not need to refit:

exp(10 * coef(fit)["lwt"])                      # OR per 10 lb
#> 0.8615

exp(10 * confint.default(fit)["lwt", ])         # its Wald 95% CI
#> 0.7504 0.9889

“Each 10 lb of maternal weight is associated with a 14% reduction in the odds of low birth weight (OR 0.86, 95% CI 0.75 to 0.99)” is a reportable sentence. Note that you must rescale the coefficient and the interval bounds on the log scale, then exponentiate — raising the odds ratio itself to a power is the same arithmetic but easier to get wrong.

Confidence intervals: confint() is not confint.default()

This is the distinction most tutorials skip, and it is the one that changes numbers.

round(exp(confint(fit)), 3)            # profile likelihood
#> Waiting for profiling to be done...

            2.5 % 97.5 %
(Intercept) 0.175 22.277
age         0.892  1.036
lwt         0.971  0.998
raceblack   1.186  9.734
raceother   0.938  5.501
smokeyes    1.058  5.293
ptdyes      1.382  8.602
htyes       1.640 27.258
uiyes       0.812  5.058

round(exp(confint.default(fit)), 3)    # Wald

            2.5 % 97.5 %
(Intercept) 0.170 21.078
age         0.894  1.037
lwt         0.972  0.999
raceblack   1.184  9.548
raceother   0.928  5.382
smokeyes    1.048  5.187
ptdyes      1.369  8.408
htyes       1.585 24.954
uiyes       0.822  5.047

confint() on a glm object gives profile-likelihood intervals. confint.default() gives Wald intervals — the familiar estimate ± 1.96 × SE, computed on the log-odds scale and then exponentiated. They are different methods answering the same question, and R gives you the better one by default.

The Wald interval assumes the log-likelihood is quadratic around the estimate. The profile interval does not: it inverts the likelihood-ratio test, finding the parameter values at which the deviance rises by the critical χ² amount. When the likelihood is skewed — sparse cells, strong effects, small samples — the two diverge, and the profile interval is the one to trust.

Look at htyes: the profile upper bound is 27.26, the Wald upper bound 24.95, roughly 9% narrower. That is not rounding. Only 12 of 189 mothers had a history of hypertension, so that coefficient’s likelihood is visibly asymmetric. Compare age, which is well-populated and near-null: 0.892–1.036 profile against 0.894–1.037 Wald, agreeing to three decimals. The rule: the sparser the cell and the larger the effect, the more the two disagree, and the more the Wald interval understates the upper bound.

In an extreme case the gap becomes decisive. Fitting glm(am ~ wt + hp, family = binomial, data = mtcars) — 32 observations, near-separated — gives a profile interval for wt of −17.21 to −3.77 against a Wald interval of −14.10 to −2.07 on the log-odds scale. Exponentiate those and they are not the same finding.

The library(MASS) instruction is out of date

Almost every published tutorial tells you to load MASS before calling confint() on a GLM. That has not been necessary since R 4.4.0. The R 4.4.0 release notes state that the confint() methods for "glm" and "nls" objects “have been copied to the stats package. Previously, they were stubs which called versions in package MASS. The MASS namespace is no longer loaded if you invoke (say) confint(glmfit)“, along with the profile() method for "glm". Confirmed in R 4.6.1:

environmentName(environment(getS3method("confint", "glm")))
#> "stats"

The practical consequence is not that the old advice breaks — loading MASS is harmless — but that a script which fails with “could not find function” is failing for some other reason, and that confint() now works in a bare session with no packages attached at all.

A publication-ready odds-ratio table in four lines

or_tab <- data.frame(
  OR    = exp(coef(fit)),
  lower = exp(confint.default(fit))[, 1],
  upper = exp(confint.default(fit))[, 2],
  p     = summary(fit)$coefficients[, 4]
)
round(or_tab, 3)

               OR lower  upper     p
(Intercept) 1.891 0.170 21.078 0.605
age         0.963 0.894  1.037 0.318
lwt         0.985 0.972  0.999 0.034
raceblack   3.363 1.184  9.548 0.023
raceother   2.235 0.928  5.382 0.073
smokeyes    2.331 1.048  5.187 0.038
ptdyes      3.393 1.369  8.408 0.008
htyes       6.288 1.585 24.954 0.009
uiyes       2.036 0.822  5.047 0.125

Swap in confint(fit) for the profile version. Drop the intercept row before publication: an exponentiated intercept is the odds when every continuous predictor equals zero, which here would be a mother of age 0 weighing 0 lb. It is not a meaningful quantity and reporting it invites a reviewer’s question you cannot answer.

Model fit: deviance, likelihood-ratio tests, and pseudo-R²

The global test

summary() gives you the two deviances but not the test. Ask for it explicitly:

anova(update(fit, . ~ 1), fit, test = "LRT")

Analysis of Deviance Table

Model 1: low ~ 1
Model 2: low ~ age + lwt + race + smoke + ptd + ht + ui
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       188        235
2       180        197  8     37.8  8.1e-06 ***

χ²(8) = 37.8, p = 8.1 × 10⁻⁶. This is the logistic analogue of the overall F test in linear regression and belongs in your results.

Testing individual terms — use drop1(), not anova()

anova(fit, test = "LRT") adds terms sequentially, in the order they appear in the formula, so its p-values depend on that order. drop1() tests each term against the full model, which is almost always what you want — and it is the only correct way to test a multi-level factor such as race, whose two coefficients must be tested jointly.

drop1(fit, test = "LRT")

Single term deletions

Model:
low ~ age + lwt + race + smoke + ptd + ht + ui
       Df Deviance AIC  LRT Pr(>Chi)
<none>         197 215
age     1      198 214 1.02   0.3130
lwt     1      202 218 5.00   0.0254 *
race    2      203 217 6.41   0.0406 *
smoke   1      201 217 4.41   0.0357 *
ptd     1      204 220 7.11   0.0076 **
ht      1      204 220 7.18   0.0074 **
ui      1      199 215 2.32   0.1279

Note the race row: 2 df, p = 0.041 for the factor as a whole. Neither individual Wald p-value in summary() tells you that, and reporting only raceblack would misstate the finding. The sequential anova() on the same model returns p = 0.107 for race, because there it is entered third rather than tested last — a live demonstration of why the two functions are not interchangeable.

Pseudo-R², computed from first principles

There is no R² in logistic regression, and R does not print one. Several pseudo-R² measures exist; all of them are computable from two log-likelihoods and the sample size, with no extra package:

LL_full <- as.numeric(logLik(fit))
LL_null <- as.numeric(logLik(update(fit, . ~ 1)))
n       <- nobs(fit)

mcfadden   <- 1 - LL_full / LL_null
coxsnell   <- 1 - exp((2 / n) * (LL_null - LL_full))
nagelkerke <- coxsnell / (1 - exp((2 / n) * LL_null))

c(LL_full = LL_full, LL_null = LL_null, n = n,
  McFadden = mcfadden, CoxSnell = coxsnell, Nagelkerke = nagelkerke)

   LL_full    LL_null          n   McFadden   CoxSnell Nagelkerke
  -98.4170  -117.3360   189.0000     0.1612     0.1814     0.2551

Three numbers, one model, and they range from 0.16 to 0.26. That spread is the whole argument against quoting a pseudo-R² without naming which one it is. McFadden’s is the most commonly reported and the most conservative; the adjusted form, 1 - (LL_full - k) / LL_null where k is the number of coefficients, gives 0.0845 here. None of them is comparable to a linear-regression R², to each other, or across models fitted to different outcomes.

Discrimination: ROC and AUC with pROC

pROC is the standard R package for this (Robin et al., BMC Bioinformatics 2011;12:77, doi:10.1186/1471-2105-12-77). Install with install.packages("pROC").

library(pROC)

p_hat <- predict(fit, type = "response")
roc_obj <- roc(bw$low, p_hat, levels = c("normal", "low"), direction = "<")

roc_obj
ci.auc(roc_obj)

Call:
roc.default(response = bw$low, predictor = p_hat, levels = c("normal", "low"), direction = "<")

Data: p_hat in 130 controls (bw$low normal) < 59 cases (bw$low low).
Area under the curve: 0.758

95% CI: 0.685-0.83 (DeLong)

Set levels and direction explicitly. Omit them and pROC guesses, printing two messages to the console:

Setting levels: control = normal, case = low
Setting direction: controls < cases

Those are messages, not warnings, so they vanish in a knitted report or a script run with --quiet. pROC infers direction from the data, which means a model that discriminates backwards can silently be flipped and reported as AUC 0.76 when it should be 0.24. Specifying both arguments makes the orientation a stated assumption rather than an inference.

AUC 0.758 (95% CI 0.685–0.830) is honest discrimination for a handful of routinely-collected risk factors. It is also an apparent AUC, computed on the same 189 rows used to fit the model, and therefore optimistic. Report it as such, or correct it by bootstrap or cross-validation — see bootstrapping in statistics for the resampling machinery.

The 0.5 threshold is a modelling choice, not a default

table(predicted = ifelse(p_hat > 0.5, "low", "normal"), observed = bw$low)

         observed
predicted normal low
   low        14  25
   normal    116  34

Sensitivity is 25/59 = 42.4%; specificity is 116/130 = 89.2%. Overall accuracy is 74.6%, and a model that is “75% accurate” while missing more than half the cases it exists to detect is not useful at that cut-point. Because only 31% of the sample are cases, a 0.5 threshold is severe by default. Choose the threshold from the clinical or operational cost of each error type — coords(roc_obj, "best") in pROC will find one by Youden’s index — and see sensitivity vs. specificity for the trade-off in full.

Calibration: are the predicted probabilities right?

AUC measures ranking only. A model that ranks perfectly but predicts 0.9 where the truth is 0.4 has an AUC of 1.0 and is unusable. Calibration is the separate question, and a decile table answers it with no package at all:

grp <- cut(p_hat, breaks = quantile(p_hat, probs = seq(0, 1, 0.1)),
           include.lowest = TRUE, labels = FALSE)

data.frame(decile    = 1:10,
           n         = as.vector(table(grp)),
           mean_pred = round(tapply(p_hat, grp, mean), 3),
           observed  = round(tapply(as.integer(bw$low == "low"), grp, mean), 3))

 decile  n mean_pred observed
      1 19     0.057    0.000
      2 19     0.102    0.105
      3 19     0.154    0.263
      4 19     0.207    0.158
      5 19     0.243    0.211
      6 18     0.283    0.389
      7 19     0.341    0.316
      8 19     0.450    0.421
      9 19     0.551    0.579
     10 19     0.731    0.684

Predicted and observed track each other across the range — the deviations in deciles 3 and 6 are two or three subjects in a cell of 19, which is noise at this sample size. This is the same grouping the Hosmer-Lemeshow test formalises; the table is more informative than the test’s single p-value, which is sensitive to the arbitrary choice of group count and has been criticised for exactly that reason.

Separation: what R tells you, and what to do

Complete or quasi-complete separation is the failure mode specific to logistic regression: a predictor (or combination) perfectly predicts the outcome, the maximum-likelihood estimate runs off to infinity, and the algorithm stops when it hits a tolerance rather than a maximum. R warns, but the warning is easy to miss because the model object still prints.

sep <- data.frame(y = c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1),
                  x = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
sep_fit <- glm(y ~ x, family = binomial, data = sep)

Warning messages:
1: glm.fit: algorithm did not converge
2: glm.fit: fitted probabilities numerically 0 or 1 occurred

summary(sep_fit)$coefficients

            Estimate Std. Error     z value Pr(>|z|)
(Intercept)  -245.85     337834 -0.00072772  0.99942
x              44.70      61172  0.00073072  0.99942

Note the signature: an enormous coefficient with an enormous standard error and a p-value of essentially 1. The exponentiated coefficient would be astronomical and the Wald interval would span everything. That combination — not the warning text — is the reliable diagnostic, because quasi-separation on one level of one factor can produce it without triggering the convergence warning at all.

Three practical checks and remedies:

• Cross-tabulate every categorical predictor against the outcome before fitting. A zero cell is separation. In birthwt, table(bw$ht, bw$low) returns 125/52 for no hypertension and 5/7 for yes — only 12 subjects and no zero cell, which is why htyes has a wide but finite interval rather than a broken one.
• Watch the Fisher scoring iteration count. Four is normal. Twenty-five, the default maxit, means it never converged.
• Use penalised likelihood rather than dropping the variable. Firth’s bias-reduction method (Firth, Biometrika 1993) gives finite estimates under separation and is implemented in the logistf package; Heinze and Schemper’s 2002 Statistics in Medicine paper is the standard reference for applying it to this problem. Deleting the offending predictor discards a variable that is, by definition, strongly associated with the outcome.

Predicted probabilities, and the interval mistake almost everyone makes

predict() has two scales and they behave differently:

nd <- data.frame(age = 25, lwt = 130,
                 race  = factor("white", levels = levels(bw$race)),
                 smoke = factor("yes",   levels = levels(bw$smoke)),
                 ptd   = factor("no",    levels = levels(bw$ptd)),
                 ht    = factor("no",    levels = levels(bw$ht)),
                 ui    = factor("no",    levels = levels(bw$ui)))

predict(fit, newdata = nd, type = "response")
#> 0.19801

A 25-year-old smoker weighing 130 lb with no other risk factors has a predicted 19.8% probability of low birth weight. Every factor column in newdata must carry the same levels as the fitted data — that is what the levels = levels(...) boilerplate is for. A level that is not in the fitted data produces factor smoke has new level ... and no prediction at all.

For an interval, compute on the link scale and transform, never the reverse:

pl <- predict(fit, newdata = nd, type = "link", se.fit = TRUE)
plogis(pl$fit + c(-1.96, 1.96) * pl$se.fit)
#> 0.1124 0.3249

Doing it the other way round — type = "response" with se.fit = TRUE, then ±1.96 SE — returns 0.0921 to 0.3040 for the same subject. Both look plausible, which is what makes this mistake durable. The response-scale version is a symmetric interval on a bounded quantity: wrong by a little in the middle of the range, and wrong by a lot near the edges, where it will happily hand back a probability below zero. plogis() is R’s inverse-logit; qlogis() goes the other way.

Average marginal effects: from odds ratios to something a reader can use

An odds ratio of 2.33 for smoking does not mean smokers are 2.33 times as likely. To get a probability-scale answer, predict the whole sample twice — once as if nobody smoked, once as if everybody did — and average:

d0 <- bw; d0$smoke <- factor("no",  levels = levels(bw$smoke))
d1 <- bw; d1$smoke <- factor("yes", levels = levels(bw$smoke))

mean(predict(fit, newdata = d0, type = "response"))   #> 0.2520
mean(predict(fit, newdata = d1, type = "response"))   #> 0.4039

An average marginal effect of +15.2 percentage points, and an implied risk ratio of 1.60 against an odds ratio of 2.33. This is the gap between odds and risk made concrete: with an outcome this common, the odds ratio overstates the relative risk substantially. The mechanism is covered in the logit model guide and in absolute risk reduction, relative risk and NNT.

Two silent problems worth checking every time

Missing data are dropped without comment

bwna <- bw
bwna$lwt[1:20] <- NA
fna <- glm(low ~ lwt + smoke, family = binomial, data = bwna)

nobs(fna)                 #> 169
length(fna$na.action)     #> 20

R’s default na.action = na.omit removes any row with a missing value in any model variable, prints nothing, and reports the reduced n only if you ask. Two models fitted to “the same data” with different predictor sets can therefore be fitted to different rows — which makes anova() and AIC() comparisons between them invalid. Always print nobs() for every model you compare, and see tidy data rules for structuring research data for getting the frame right before it reaches glm().

Events per variable

sum(bw$low == "low") / (length(coef(fit)) - 1)
#> 7.375

59 events across 8 estimated coefficients is 7.4 events per variable, below the conventional minimum of 10 (Peduzzi et al., Journal of Clinical Epidemiology 1996). The consequence here is visible in the output: wide intervals, and coefficients such as htyes resting on 12 subjects. The 10-EPV heuristic has been superseded for prediction models by the sample-size approach of Riley and colleagues, but as a quick sanity check on an explanatory model it still earns its place. Influence is worth a glance too — cooks.distance(fit) peaks at 0.048 here, well under any conventional cut-off, so no single mother is driving the fit.

What to report

For an explanatory model, a journal will expect: the number of observations analysed and the number of events; which outcome level was modelled; the full list of predictors with no post-hoc selection left undeclared; odds ratios with confidence intervals and the interval method named; and the global likelihood-ratio test. Discrimination and calibration statistics belong with any model presented for prediction.

STROBE item 12 requires the statistical methods, including those used to control for confounding, and item 16 requires unadjusted and confounder-adjusted estimates with their precision. For a prediction model, TRIPOD+AI is the applicable checklist and asks explicitly for how predictors were handled, how the model was developed, and both discrimination and calibration. Interval method and pseudo-R² variant are exactly the details these checklists exist to stop authors from leaving out — and see how to report confidence intervals in a manuscript and how to report p values for the formatting conventions.

Record your sessionInfo() alongside the results. Package versions change defaults, and the R 4.4.0 confint() migration above is a live example of a change that alters which library a reproducible script depends on.

Error and warning quick reference

What R says	What it means
glm.fit: algorithm did not converge	Almost always separation. Check for zero cells; consider Firth penalisation.
glm.fit: fitted probabilities numerically 0 or 1 occurred	Same cause. Can also appear on its own with quasi-separation.
factor smoke has new level maybe	A newdata factor level absent from the fitted data. Build newdata factors with levels = levels(original).
y values must be 0 <= y <= 1	A numeric outcome outside [0, 1] passed to binomial. Recode to 0/1 or to a factor.
non-integer #successes in a binomial glm!	Proportions passed without weights. Supply the denominators as weights.
Waiting for profiling to be done...	Not an error. confint() is computing profile-likelihood intervals.
Coefficient shown as NA	Perfect collinearity; R aliased the term out. See multicollinearity and VIF.

Frequently asked questions

How do I get odds ratios in R?

exp(coef(model)) for the point estimates and exp(confint(model)) for profile-likelihood confidence intervals. Exponentiate the log-odds coefficients and their interval bounds — never the standard errors, which have no interpretation on the odds scale.

What does family = binomial do in glm()?

It specifies a Bernoulli/binomial error distribution and, by default, a logit link, which is what makes the fit a logistic regression. The default is binomial(link = "logit"); "probit", "cloglog" and others are available and change the coefficient interpretation — only the logit gives you odds ratios.

Why does confint() give different numbers from confint.default()?

confint() inverts the likelihood-ratio test to produce profile-likelihood intervals; confint.default() produces Wald intervals from estimate ± 1.96 × SE. They agree when the log-likelihood is close to quadratic and diverge with sparse data or strong effects. Prefer the profile interval and say which you used.

How do I calculate pseudo-R² in R?

From logLik() on the fitted and null models: McFadden’s is 1 - LL_full / LL_null; Cox-Snell and Nagelkerke follow from the same two values plus nobs(). Packages such as DescTools and pscl wrap this, but the arithmetic is three lines and reporting which measure you used matters more than which function produced it.

How do I plot a ROC curve for a glm in R?

plot(roc(outcome, predict(model, type = "response"))) using pROC, with levels and direction set explicitly. Add ci.auc() for a DeLong confidence interval on the AUC. Remember that an AUC computed on the training data is optimistic.

Does the outcome need to be a factor?

No. A 0/1 integer works and gives identical coefficients. A factor is safer for labelling but introduces the level-order trap: R models the probability of the second level, so check levels() before interpreting anything.

How many predictors can my model support?

A long-standing rule of thumb is at least 10 events per estimated coefficient, counting the smaller of the two outcome classes. The example above sits at 7.4 and shows the symptom — wide intervals on sparse predictors. Modern prediction-model guidance replaces the rule with a formal sample-size calculation, but as a pre-fit sanity check it remains useful.

What is the R equivalent of SPSS’s logistic regression output?

summary() covers the coefficient table and model fit; exp(coef()) plus exp(confint()) replaces the Exp(B) column with its interval; drop1(test = "LRT") replaces the omnibus tests; and Nagelkerke’s pseudo-R² must be computed from logLik(). See SPSS vs. R for the wider comparison, or what SPSS is for the other side.

Related CASRAI resources

• Logistic regression and the logit model — the method behind the code: link functions, odds ratios versus risk ratios, multinomial and ordinal variants.
• Regression analysis: assumptions, interpretation and reporting
• Multicollinearity and VIF in regression
• Sensitivity vs. specificity, PPV/NPV and ROC curves
• Case-control study design and odds ratios
• Poisson models for count outcomes — the other common glm() family
• Mixed-effects models — when observations are clustered and glm() is not enough
• R vs. Stata and SPSS vs. R
• Research tools — the full cluster, including statistical and qualitative software.