# Logistic Regression: [Your study question or course] > A method card for binary logistic regression: when to use it, the logit > link, odds and odds ratios, how to read each coefficient, Wald and > likelihood-ratio tests, profile and Wald confidence intervals, predicted > probabilities and average marginal effects, classification metrics with > their threshold caveats, calibration, separation and events per variable. > It is for students preparing for a course or exam, thesis writers with a > yes/no outcome and analysts who build or check simple risk models. It > speeds up learning by fixing the order of decisions, giving a fully worked > example in R and Python, then fading the support in a half-solved example > and an exercise, and testing you before you see the answers. The worked > examples use fictional data; replace them with your own. > Based on generalised linear model theory (Nelder and Wedderburn 1972; Cox > 1958), the NIST/SEMATECH e-Handbook of Statistical Methods (public domain) > and the method papers listed at the end. ## Document control | Field | Value | |---|---| | Course or project | [STAT 310 Generalised linear models / thesis chapter 5] | | Learner or analyst | [Your name] | | Data set | [*Example:* Lakemoor University pass study, n = 30, fictional] | | Software and versions | [*Example:* R 4.5.1 (stats); Python 3.14 with numpy 2.5.4, pandas 3.1.0, scipy 1.18.1, statsmodels 0.15.0] | | Date | YYYY-MM-DD | | Version | [0.1] | ## When to use it **One-sentence rule:** use logistic regression when the outcome has exactly two categories (pass or fail, defect or no defect) and you want to model how the probability of the event changes with one or more predictors. | Use when | Do not use when | |---|---| | The outcome is binary and coded 0/1 (1 = the event) | The outcome is continuous: use linear regression | | You want odds ratios adjusted for other predictors | The outcome is a count of events or a rate: use Poisson or negative binomial regression | | You want predicted probabilities for new cases | The outcome has three or more categories: use multinomial or ordinal logistic regression | | Observations are independent | Rows are clustered (students within classes, repeated tests per person): use a mixed or GEE logistic model | | You have enough events for the number of predictors | One predictor perfectly separates the two outcomes: use Firth's penalised logistic regression or exact methods | **Alternatives:** a chi-square or Fisher test for one categorical predictor without adjustment; the log-binomial or modified Poisson model when you want risk ratios rather than odds ratios (common outcomes, where odds ratios exaggerate); probit regression gives nearly the same fitted probabilities with a different link; classification trees or penalised models for pure prediction with many predictors. ## The question and the data **Example question:** at Lakemoor University (fictional), how do weekly study hours and attending the optional pre-exam review session relate to the probability of passing the introductory statistics exam? | Variable | Role | Measurement level | Coding | |---|---|---|---| | `passed` | Outcome | Binary | 1 = passed, 0 = failed | | `hours` | Predictor | Continuous (ratio) | Mean weekly study hours | | `review` | Predictor | Binary | 1 = attended the review session, 0 = did not (reference) | **Design:** cross-sectional observational study, sample size n = 30 students, 19 passed and 11 failed. The smaller of the two counts (11) is the number that limits model size: with 2 predictors that is 11 / 2 = 5.5 events per variable (EPV), below the traditional rule of thumb of 10 (Peduzzi et al. 1996). The example is fine for learning, but in a real study it would be a reason for caution about the precision of the estimates. | Student | hours | review | passed | Student | hours | review | passed | |---|---|---|---|---|---|---|---| | 1 | 1 | 0 | 0 | 16 | 7 | 1 | 0 | | 2 | 2 | 0 | 0 | 17 | 8 | 0 | 1 | | 3 | 2 | 1 | 0 | 18 | 8 | 1 | 1 | | 4 | 3 | 0 | 0 | 19 | 8 | 1 | 1 | | 5 | 3 | 1 | 1 | 20 | 9 | 0 | 0 | | 6 | 4 | 0 | 0 | 21 | 9 | 1 | 1 | | 7 | 4 | 0 | 0 | 22 | 10 | 0 | 1 | | 8 | 5 | 1 | 0 | 23 | 10 | 1 | 1 | | 9 | 5 | 0 | 1 | 24 | 11 | 0 | 1 | | 10 | 5 | 1 | 1 | 25 | 11 | 1 | 1 | | 11 | 6 | 0 | 0 | 26 | 12 | 1 | 1 | | 12 | 6 | 1 | 1 | 27 | 12 | 0 | 1 | | 13 | 6 | 0 | 0 | 28 | 13 | 1 | 1 | | 14 | 7 | 1 | 1 | 29 | 14 | 0 | 1 | | 15 | 7 | 0 | 1 | 30 | 15 | 1 | 1 | The chart shows the outcome (0 or 1, slightly jittered by colour) against study hours, with the fitted probability of passing from the one-predictor model `passed ~ hours` (Example B). The S-shaped curve is the point of the logit link: probabilities stay between 0 and 1, and the effect of one more hour on the probability is largest near p = 0.5. ```vega-lite { "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "description": "Pass or fail against study hours with the fitted logistic curve", "width": 380, "height": 240, "layer": [ { "data": {"values": [ {"hours": 1, "passed": 0}, {"hours": 2, "passed": 0}, {"hours": 2, "passed": 0}, {"hours": 3, "passed": 0}, {"hours": 3, "passed": 1}, {"hours": 4, "passed": 0}, {"hours": 4, "passed": 0}, {"hours": 5, "passed": 0}, {"hours": 5, "passed": 1}, {"hours": 5, "passed": 1}, {"hours": 6, "passed": 0}, {"hours": 6, "passed": 1}, {"hours": 6, "passed": 0}, {"hours": 7, "passed": 1}, {"hours": 7, "passed": 1}, {"hours": 7, "passed": 0}, {"hours": 8, "passed": 1}, {"hours": 8, "passed": 1}, {"hours": 8, "passed": 1}, {"hours": 9, "passed": 0}, {"hours": 9, "passed": 1}, {"hours": 10, "passed": 1}, {"hours": 10, "passed": 1}, {"hours": 11, "passed": 1}, {"hours": 11, "passed": 1}, {"hours": 12, "passed": 1}, {"hours": 12, "passed": 1}, {"hours": 13, "passed": 1}, {"hours": 14, "passed": 1}, {"hours": 15, "passed": 1} ]}, "mark": {"type": "point", "filled": true, "size": 60, "opacity": 0.6}, "encoding": { "x": {"field": "hours", "type": "quantitative", "title": "Weekly study hours"}, "y": {"field": "passed", "type": "quantitative", "title": "Passed (1) or probability of passing"} } }, { "data": {"values": [ {"hours": 0, "p": 0.037}, {"hours": 1, "p": 0.064}, {"hours": 2, "p": 0.109}, {"hours": 3, "p": 0.179}, {"hours": 4, "p": 0.280}, {"hours": 5, "p": 0.410}, {"hours": 6, "p": 0.554}, {"hours": 7, "p": 0.689}, {"hours": 8, "p": 0.798}, {"hours": 9, "p": 0.876}, {"hours": 10, "p": 0.927}, {"hours": 11, "p": 0.958}, {"hours": 12, "p": 0.976}, {"hours": 13, "p": 0.986}, {"hours": 14, "p": 0.992}, {"hours": 15, "p": 0.996}, {"hours": 16, "p": 0.998} ]}, "mark": {"type": "line", "color": "#555555", "interpolate": "monotone"}, "encoding": { "x": {"field": "hours", "type": "quantitative"}, "y": {"field": "p", "type": "quantitative"} } } ] } ``` ## Assumptions and how to check each | Assumption | How to check | R | Python | If violated | |---|---|---|---|---| | Binary outcome, one row per independent unit | Study design; `table(y)` | `table(d$passed)` | `d["passed"].value_counts()` | Clustered data: mixed logistic (`lme4::glmer`) or GEE | | Linearity of the log-odds in each continuous predictor | Add a squared term or spline and compare by likelihood ratio; binned residual plot | `anova(fit, update(fit, . ~ . + I(hours^2)), test = "LRT")` | Compare `llf` of the two fits; `2 * (llf1 - llf0)` | Spline or polynomial term, or a transformation (log of a skewed predictor) | | Enough events per variable | Count the smaller outcome group and divide by the number of coefficients (not counting the intercept) | `min(table(d$passed)) / 2` | `d["passed"].value_counts().min() / 2` | Fewer predictors, penalised (Firth) regression, collect more data; for prediction models use formal sample size methods | | No (quasi-)complete separation | Huge coefficients and SEs, warning "fitted probabilities numerically 0 or 1" | Read warnings; `summary(fit)` | `PerfectSeparationWarning` or huge `bse` | Firth regression (`logistf` package), exact logistic, or merge sparse categories | | No severe multicollinearity | VIF on the same predictors | `car::vif(fit)` | `variance_inflation_factor` on the design matrix | As for linear regression | | No overly influential points | Cook's distance and dfbeta for GLMs | `plot(fit, which = 5)`; `influence.measures(fit)` | `fit_glm.get_influence().summary_frame()` | Check for errors, report with and without | Notice what is **not** in the table: residuals are not expected to be normal, and their variance is not constant (it is $p(1-p)$ by design). Normality tests on logistic residuals are a common wasted step. ## How it works **In plain language.** A straight line for a probability would give values below 0 and above 1. Logistic regression instead models the **log-odds** of the event as a straight line, and then turns the log-odds back into a probability with an S-shaped curve. Each coefficient is the change in log-odds per unit of the predictor, holding the others constant. Exponentiating it gives an **odds ratio**: the factor by which the odds of the event are multiplied for a one-unit increase. The coefficients are found by maximum likelihood: the values that make the observed passes and fails most probable. **Three scales.** For probability $p$: $$ \text{odds} = \frac{p}{1 - p}, \qquad \operatorname{logit}(p) = \ln\!\left(\frac{p}{1 - p}\right), \qquad p = \frac{1}{1 + e^{-\eta}} $$ **Model** for observation $i$, with linear predictor $\eta_i$: $$ \ln\!\left(\frac{p_i}{1 - p_i}\right) = \eta_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik}, \qquad Y_i \sim \text{Bernoulli}(p_i) $$ **Odds ratio** for a one-unit increase in $x_j$, and for a $c$-unit increase: $$ \text{OR}_j = e^{\beta_j}, \qquad \text{OR}_j(c) = e^{c\,\beta_j} $$ **Log-likelihood** maximised by the software (iteratively reweighted least squares): $$ \ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left[ y_i \ln p_i + (1 - y_i)\ln(1 - p_i) \right] $$ **Tests.** Wald test of one coefficient, and likelihood-ratio (LR) test of a full model against a reduced one with $q$ fewer coefficients, where $D$ is the deviance $-2\ell$: $$ z = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)}, \qquad G^2 = D_{\text{reduced}} - D_{\text{full}} = 2(\ell_{\text{full}} - \ell_{\text{reduced}}) \sim \chi^2_q $$ **Average marginal effect (AME)** of a continuous predictor on the probability scale, and of a binary predictor (discrete change): $$ \text{AME}_j = \frac{1}{n}\sum_{i=1}^{n} \hat{\beta}_j\, \hat{p}_i (1 - \hat{p}_i), \qquad \text{AME}_{\text{binary}} = \frac{1}{n}\sum_{i=1}^{n} \left[\hat{p}_i(x = 1) - \hat{p}_i(x = 0)\right] $$ | Symbol | Meaning | |---|---| | $p_i$ | Probability of the event (Y = 1) for observation $i$ | | $\eta_i$ | Linear predictor, the log-odds of the event | | $\beta_j$ | Change in log-odds per unit of $x_j$, others held constant | | $e^{\beta_j}$ | Odds ratio per unit of $x_j$ | | $\ell$, $D$ | Log-likelihood and deviance ($D = -2\ell$ for binary data) | | $G^2$ | Likelihood-ratio statistic, compared with $\chi^2$ on $q$ df | | AME | Average change in predicted probability, in probability units | | EPV | Events per variable: smaller outcome count divided by the number of coefficients | **Profile vs Wald intervals.** The Wald interval is $\hat{\beta} \pm 1.96\,SE$, symmetric on the log-odds scale. The profile likelihood interval follows the actual shape of the likelihood and is usually more accurate in small samples or near separation. In R, `confint()` on a `glm` gives **profile** intervals (since R 4.4.0 the method lives in `stats`; earlier versions borrowed it from `MASS`), while `confint.default()` gives Wald intervals. In statsmodels, `conf_int()` gives **Wald** intervals. So R and Python will print different CIs for the same model, and both are "correct" for their method. **Classification.** Turning probabilities into yes/no predictions needs a threshold $t$: predict 1 when $\hat{p} \ge t$. Sensitivity (true positive rate) and specificity (true negative rate) depend on $t$; the ROC curve plots them for every threshold, and the area under it (AUC, or the c-statistic) is the probability that a randomly chosen event has a higher predicted probability than a randomly chosen non-event. **Calibration** asks a different question: among cases given about 70%, do about 70% have the event? ```mermaid flowchart TD A["Binary outcome coded 0 and 1"] --> B["Count events and compute EPV"] B --> C["Fit glm or Logit"] C --> D{"Warnings or huge SEs?"} D -->|Yes| E["Separation: Firth or exact methods"] D -->|No| F["Check log-odds linearity and influence"] F --> G{"Purpose?"} G -->|Explain| H["Odds ratios with profile CIs and LR tests"] G -->|Predict| I["AUC, calibration and validation"] H --> J["Translate to probabilities and AMEs"] I --> J J --> K["Report"] ``` ## Worked example A (full): passing the exam **Hypotheses.** For each predictor, $H_0: \beta_j = 0$ (odds ratio 1, no association with passing once the other predictor is held constant) against $H_1: \beta_j \neq 0$, two-sided, α = 0.05. Overall: $H_0$: neither predictor is related to passing. **Checks before reading the results.** Events: 19 passes and 11 fails, so EPV = 5.5 (low; noted as a limitation). No separation warning appears, and the standard errors are of normal size. Adding `I(hours^2)` lowers the deviance by only 0.06 (LR test *p* = .803), so log-odds linearity in hours is acceptable. **R code.** ```r d <- data.frame( hours = c(1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 14, 15), review = c(0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1), passed = c(0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)) fit <- glm(passed ~ hours + review, family = binomial, data = d) summary(fit) confint(fit) # profile likelihood CIs (stats::confint, R >= 4.4) confint.default(fit) # Wald CIs: estimate +/- 1.96 SE round(exp(cbind(OR = coef(fit), confint(fit))), 3) drop1(fit, test = "LRT") # likelihood-ratio test for each term # Predicted probabilities and average marginal effects (AME) p <- fitted(fit) predict(fit, data.frame(hours = c(4, 8), review = 1), type = "response") mean(p * (1 - p)) * coef(fit)["hours"] # AME of one more hour d1 <- transform(d, review = 1); d0 <- transform(d, review = 0) mean(predict(fit, d1, type = "response") - predict(fit, d0, type = "response")) # Classification at threshold 0.5, and AUC via the rank (Mann-Whitney) formula table(predicted = as.integer(p >= 0.5), actual = d$passed) r <- rank(p); n1 <- sum(d$passed); n0 <- sum(1 - d$passed) (sum(r[d$passed == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0) # Calibration: observed vs mean predicted in three risk groups grp <- cut(p, c(0, 1/3, 2/3, 1)) aggregate(cbind(observed = d$passed, predicted = p) ~ grp, FUN = mean) ``` **R output** (R 4.5.1; output abridged: some calls are shown with a `>` prompt and a short comment so you can see which line printed what, and R's "Waiting for profiling to be done..." messages are omitted). ```text Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.1784 1.6771 -2.491 0.01272 * hours 0.5961 0.2273 2.623 0.00871 ** review 1.7397 1.1060 1.573 0.11572 Null deviance: 39.429 on 29 degrees of freedom Residual deviance: 22.248 on 27 degrees of freedom AIC: 28.248 > confint(fit) # profile 2.5 % 97.5 % (Intercept) -8.2697881 -1.441756 hours 0.2326828 1.159566 review -0.2981269 4.211865 > confint.default(fit) # Wald 2.5 % 97.5 % (Intercept) -7.4654175 -0.8913465 hours 0.1507225 1.0415589 review -0.4280001 3.9073975 OR 2.5 % 97.5 % (Intercept) 0.015 0.000 0.237 hours 1.815 1.262 3.189 review 5.696 0.742 67.482 Df Deviance AIC LRT Pr(>Chi) 22.248 28.248 hours 1 35.740 39.740 13.4918 0.0002396 *** review 1 25.021 29.021 2.7735 0.0958393 . > predict(fit, data.frame(hours = c(4, 8), review = 1), type = "response") 1 2 0.4864732 0.9113671 > mean(p * (1 - p)) * coef(fit)["hours"] hours 0.06996238 > mean(predict(fit, d1, ...) - predict(fit, d0, ...)) [1] 0.2263289 actual predicted 0 1 0 8 3 1 3 16 [1] 0.9043062 grp observed predicted 1 (0,0.333] 0.1428571 0.1284975 2 (0.333,0.667] 0.5714286 0.4939522 3 (0.667,1] 0.8750000 0.9151782 ``` **Python code.** ```python import numpy as np import pandas as pd import statsmodels.api as sm import statsmodels.formula.api as smf from scipy.stats import mannwhitneyu d = pd.DataFrame({ "hours": [1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 14, 15], "review": [0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1], "passed": [0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]}) fit = smf.logit("passed ~ hours + review", data=d).fit(disp=0) print(fit.summary().tables[1]) print("LR chi2 =", round(fit.llr, 3), " p =", round(fit.llr_pvalue, 7), " McFadden R2 =", round(fit.prsquared, 4)) print("OR:", np.exp(fit.params).round(3).to_dict()) print("Wald 95% CI for OR:\n", np.exp(fit.conf_int()).round(3)) # Average marginal effects; dummy=True gives the discrete 0 -> 1 change for review print(fit.get_margeff(at="overall", dummy=True).summary().tables[1]) # Same model as a GLM (identical estimates; deviance and AIC as in R) glm = smf.glm("passed ~ hours + review", data=d, family=sm.families.Binomial()).fit() print("GLM deviance =", round(glm.deviance, 3), " AIC =", round(glm.aic, 3)) p = fit.predict(d) print(pd.crosstab((p >= 0.5).astype(int), d["passed"], rownames=["predicted"])) u = mannwhitneyu(p[d.passed == 1], p[d.passed == 0]).statistic print("AUC =", round(u / (19 * 11), 4)) ``` **Python output** (Python 3.14, statsmodels 0.15.0, scipy 1.18.1; output abridged: the `====` table borders are omitted). ```text coef std err z P>|z| [0.025 0.975] ------------------------------------------------------------------------------ Intercept -4.1784 1.677 -2.491 0.013 -7.465 -0.891 hours 0.5961 0.227 2.623 0.009 0.151 1.042 review 1.7397 1.106 1.573 0.116 -0.428 3.907 LR chi2 = 17.181 p = 0.0001858 McFadden R2 = 0.4358 OR: {'Intercept': 0.015, 'hours': 1.815, 'review': 5.696} Wald 95% CI for OR: 0 1 Intercept 0.001 0.410 hours 1.163 2.834 review 0.652 49.769 dy/dx std err z P>|z| [0.025 0.975] ------------------------------------------------------------------------------ hours 0.0700 0.011 6.176 0.000 0.048 0.092 review 0.2263 0.136 1.662 0.096 -0.041 0.493 GLM deviance = 22.248 AIC = 28.248 passed 0 1 predicted 0 8 3 1 3 16 AUC = 0.9043 ``` **Where R and Python differ, and why.** The estimates, SEs, Wald z tests, deviance and AIC agree exactly. The confidence intervals do not: R's `confint()` gives profile intervals (OR for review [0.742, 67.482]) while statsmodels gives Wald intervals (OR [0.652, 49.769]); the Wald intervals in R (`confint.default`) match Python. The marginal effect of `review` also depends on a setting: with `dummy=True` statsmodels computes the discrete change from 0 to 1 (0.2263, the same as the R calculation); without it, statsmodels treats `review` as continuous and reports a derivative of 0.2042. For a 0/1 predictor the discrete change is the one to report. **Calculation by hand for the key statistics.** - Odds ratio for hours: $e^{0.5961} = 1.815$. Each extra weekly hour multiplies the odds of passing by 1.82, an 82% increase in the odds. - For 3 extra hours: $e^{3 \times 0.5961} = e^{1.788} = 5.98$. - Predicted probability for a review attender who studies 4 hours: $\eta = -4.1784 + 0.5961 \times 4 + 1.7397 = -0.0543$, so $p = 1/(1 + e^{0.0543}) = 0.486$. - LR test for the whole model: $G^2 = 39.429 - 22.248 = 17.181$ on 2 df, $p < .001$. - From the confusion matrix (threshold 0.5): sensitivity $= 16/19 = 0.842$, specificity $= 8/11 = 0.727$, accuracy $= 24/30 = 0.80$. **Wald vs likelihood-ratio for review.** The Wald test gives *p* = .116, the LR test *p* = .096. They disagree because the Wald test relies on a symmetric normal approximation that is poor when a coefficient is large relative to the information in the data (here only 30 students and a wide, skewed likelihood). In small samples prefer the LR test and the profile CI; both say the review effect is imprecisely estimated (OR 5.70, profile 95% CI [0.74, 67.48]). **Effect size with CI.** Odds ratios with profile CIs: hours OR = 1.82, 95% CI [1.26, 3.19]; review OR = 5.70, 95% CI [0.74, 67.48]. On the probability scale (AMEs), one more weekly hour raises the probability of passing by about 7.0 percentage points on average (95% CI [4.8, 9.2] from the delta method in statsmodels), and attending the review session by about 22.6 points (95% CI [-4.1, 49.3]). Discrimination: AUC = 0.90. **Interpretation.** Holding review attendance constant, students who study more are much more likely to pass: the odds of passing rise by a factor of 1.82 for each extra weekly hour. Review attendance is associated with higher odds of passing, but the interval is very wide and includes 1, so these data cannot tell a large benefit from none. The model separates passers from failers well (AUC 0.90) and is roughly calibrated in three risk groups (observed 14%, 57% and 88% against predicted 13%, 49% and 92%), but with EPV = 5.5 and no external validation, it should not be used to predict outcomes for other cohorts. **APA-style report.** A logistic regression predicted passing the exam from weekly study hours and review session attendance, χ²(2) = 17.18, *p* < .001, McFadden *R*² = .44. Each additional study hour was associated with higher odds of passing, OR = 1.82, 95% CI [1.26, 3.19] (profile), LR χ²(1) = 13.49, *p* < .001; the average marginal effect was 7.0 percentage points per hour. Review attendance was not clearly associated with passing, OR = 5.70, 95% CI [0.74, 67.48], LR χ²(1) = 2.77, *p* = .096. The model discriminated well (AUC = .90); with 11 failures and 2 predictors, estimates are imprecise. ## Worked example B (half solved): one predictor Model `passed ~ hours` on the same 30 students. R output (abridged: `coef(summary(fit1))`, and the last line summarises `exp(confint(fit1))` in words): ```text Estimate Std. Error z value Pr(>|z|) (Intercept) -3.2632042 1.3955360 -2.338316 0.019370860 hours 0.5798946 0.2155692 2.690062 0.007143865 profile 95% CI for the odds ratio of hours: [1.264, 3.039] ``` 1. Odds ratio for hours: $e^{0.5799} =$ [ ]. 2. In words: each extra weekly hour multiplies the odds of passing by [ ], that is, the odds increase by [ ]%. 3. Log-odds for a student with 5 hours: $\eta = -3.2632 + 0.5799 \times 5 =$ [ ]. 4. Probability: $p = 1/(1 + e^{-\eta}) =$ [ ]. 5. The hours at which the predicted probability is exactly 0.5 (where $\eta = 0$): $3.2632 / 0.5799 =$ [ ] hours. 6. Report sentence: OR = [ ], 95% CI [ ], *z* = [ ], *p* = [ ]. ## Exercise C (unsolved) Brightwell Bakery (fictional) records, for 20 batches of rye loaves, how far the oven temperature was from the set point (°C, negative = too cold) and whether the batch failed the quality check (`defect` = 1). The baker suspects that deviations in **either** direction cause defects. | Batch | temp_dev | defect | Batch | temp_dev | defect | |---|---|---|---|---|---| | 1 | -12 | 1 | 11 | 0 | 0 | | 2 | -10 | 1 | 12 | 1 | 0 | | 3 | -9 | 0 | 13 | 2 | 1 | | 4 | -7 | 1 | 14 | 3 | 0 | | 5 | -6 | 0 | 15 | 4 | 0 | | 6 | -5 | 0 | 16 | 5 | 0 | | 7 | -4 | 1 | 17 | 6 | 1 | | 8 | -3 | 0 | 18 | 8 | 1 | | 9 | -2 | 0 | 19 | 10 | 1 | | 10 | -1 | 0 | 20 | 12 | 1 | Tasks: (a) explain why `defect ~ temp_dev` would be the wrong model for the baker's idea, and create `abs_dev = abs(temp_dev)`; (b) fit `defect ~ abs_dev`; report the odds ratio per °C with its profile 95% CI and the LR test; (c) give the predicted probability of a defect at 0 °C and at 10 °C of deviation; (d) compute EPV and comment. ## Reading software output | Output line (R `summary(glm)`) | statsmodels equivalent | Meaning | Report it? | |---|---|---|---| | `Estimate` | `coef` | Change in log-odds per unit, others held constant | Usually as OR = exp(estimate) | | `Std. Error`, `z value`, `Pr(>\|z\|)` | `std err`, `z`, `P>\|z\|` | Wald test of $\beta_j = 0$ | Prefer the LR test in small samples | | `confint(fit)` | not built in (Wald only) | Profile likelihood CI for $\beta_j$; exponentiate for OR | Yes, say "profile" | | `confint.default(fit)` | `conf_int()` | Wald CI | Say "Wald" if you use it | | `Null deviance`, `Residual deviance` | `-2 * llnull`, `-2 * llf` (or GLM `deviance`) | Model fit; their difference is the LR statistic | Yes, as χ²(df) | | `AIC` | `aic` | Comparison of non-nested models on the same data | When comparing models | | `drop1(fit, test = "LRT")` | compare `llf` of nested fits | LR test for each term | Yes | | no direct line | `prsquared` | McFadden pseudo-R², not comparable with OLS R² | Optional, name the type | | `Number of Fisher Scoring iterations` | `converged` | Convergence; very many iterations can signal separation | Check, do not report | ## Reporting template **Fill-in sentence.** A logistic regression modelled the probability of [event] from [predictors] (n = [n], [events] events), χ²([df]) = [G²], *p* [=/<] [p]. Holding [other predictors] constant, each [unit] increase in [predictor] was associated with [higher/lower] odds of [event], OR = [OR], 95% CI [lower, upper] ([profile/Wald]), LR χ²(1) = [value], *p* [=/<] [p]; the average marginal effect was [AME] percentage points. [For prediction: AUC = [value]; calibration [result]; threshold [t] gave sensitivity [ ] and specificity [ ].] | Predictor | *b* | SE | OR | 95% CI for OR | LR χ²(1) | *p* | |---|---|---|---|---|---|---| | *Example:* Study hours (per hour) | 0.596 | 0.227 | 1.82 | [1.26, 3.19] | 13.49 | < .001 | | *Example:* Review attended (vs not) | 1.740 | 1.106 | 5.70 | [0.74, 67.48] | 2.77 | .096 | | [Predictor] | [ ] | [ ] | [ ] | [ ] | [ ] | [ ] | **Results checklist.** - [ ] Event coded 1 and reference levels stated - [ ] Number of events and non-events, and EPV, reported - [ ] Separation checked (no warnings, no huge SEs) - [ ] Odds ratios with CIs, CI method named (profile or Wald) - [ ] LR tests used, or Wald tests justified - [ ] At least one result translated to probabilities (predicted probabilities or AMEs) - [ ] For prediction aims: AUC, calibration, threshold and validation reported - [ ] Odds ratios not described as risk ratios - [ ] Software and versions named ## Common mistakes and misconceptions - **"OR = 1.82 means passing is 1.82 times as likely."** That is a risk ratio. An odds ratio multiplies the **odds**. When the event is common (here 63% passed), the odds ratio is further from 1 than the risk ratio. - **Reading a coefficient as a change in probability.** 0.596 is a change in log-odds. The change in probability depends on where you start: it is largest near p = 0.5 and small near 0 or 1. Report AMEs or predicted probabilities for that. - **Comparing odds ratios or coefficients across models with different predictors** as if only the predictor of interest changed. In logistic regression adding a predictor that is unrelated to the others (but related to the outcome) still changes the other coefficients (non-collapsibility). - **Trusting Wald p-values when SEs are large.** Use LR tests and profile intervals, especially in small samples. - **Ignoring separation warnings.** In a small perfectly separated example (8 students, all with 5 or more hours passed), R reports an estimate of 45.8 with SE 80,644 and *p* = 1.00. That is not "no effect": the maximum likelihood estimate does not exist. Use Firth regression. - **Reporting accuracy alone.** With 63% passers, predicting "pass" for everyone already gives 63% accuracy. Report sensitivity and specificity, and say which threshold you used and why. - **Treating the 0.5 threshold as part of the model.** The threshold is a decision about costs, not a statistical result. A tutor who wants to catch every student at risk will choose a different threshold from one who has few tutoring places. - **Using a high AUC as evidence of good calibration.** AUC measures ranking only. A model can rank perfectly and still predict 90% for a group where 60% pass. - **Too many predictors for the events.** With 11 failures, a model with 8 predictors is badly overfitted, whatever its apparent AUC. ## Self-test questions 1. Convert a probability of 0.80 to odds and to log-odds. 2. The coefficient of `hours` is 0.5961. State its meaning on the log-odds scale and on the odds scale in one sentence each. 3. Why do R's `confint(fit)` and statsmodels' `conf_int()` give different intervals for the same logistic model, and which would you report for n = 30? 4. With a threshold of 0.7 instead of 0.5, the confusion matrix becomes: predicted 0 / actual 0 = 9, predicted 0 / actual 1 = 5, predicted 1 / actual 0 = 2, predicted 1 / actual 1 = 14. Compute sensitivity and specificity, and say what changed. 5. You fit a model and see the warning "fitted probabilities numerically 0 or 1 occurred" and a coefficient of 38 with SE 9,000. What is happening and what do you do? 6. A study has 400 participants and 24 events. How many predictors can you reasonably include by the EPV = 10 guide, and what else should you consider? ## My notes and next steps | Item | Entry | |---|---| | Confidence with this method (1-5) | [3] | | What I can now do | [Fit a logistic model and explain an odds ratio] | | Still unclear | [When to use risk ratios instead] | | Mistakes to log in academic/stats-error-log | [Called an OR a risk ratio] | | Next review date (spacing) | YYYY-MM-DD | | Next template | [academic/stats-count-regression] | ## Sources and further reading Structure aligned with the method card used across this statistics library; formulas are standard generalised linear model results in our own notation. - Cox, D. R. (1958). The regression analysis of binary sequences. *Journal of the Royal Statistical Society, Series B*, 20(2), 215-232. - Nelder, J. A., and Wedderburn, R. W. M. (1972). Generalized linear models. *Journal of the Royal Statistical Society, Series A*, 135(3), 370-384. - Albert, A., and Anderson, J. A. (1984). On the existence of maximum likelihood estimates in logistic regression models. *Biometrika*, 71(1), 1-10. - Firth, D. (1993). Bias reduction of maximum likelihood estimates. *Biometrika*, 80(1), 27-38. - Hanley, J. A., and McNeil, B. J. (1982). The meaning and use of the area under a receiver operating characteristic (ROC) curve. *Radiology*, 143(1), 29-36. - Peduzzi, P., Concato, J., Kemper, E., Holford, T. R., and Feinstein, A. R. (1996). A simulation study of the number of events per variable in logistic regression analysis. *Journal of Clinical Epidemiology*, 49(12), 1373-1379. - Van Calster, B., McLernon, D. J., van Smeden, M., Wynants, L., and Steyerberg, E. W. (2019). Calibration: the Achilles heel of predictive analytics. *BMC Medicine*, 17, 230. - Riley, R. D., Ensor, J., Snell, K. I. E., and colleagues (2020). Calculating the sample size required for developing a clinical prediction model. *BMJ*, 368, m441. - Hosmer, D. W., Lemeshow, S., and Sturdivant, R. X. (2013). *Applied Logistic Regression* (3rd ed.). Wiley. - NIST/SEMATECH e-Handbook of Statistical Methods, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain. - Further reading by name: de Smith, M. J., *Statistical Analysis Handbook* (StatsRef, 2024 edition); OpenIntro Statistics; APA *Publication Manual* (7th ed.) and APA JARS-Quant; Wasserstein and Lazar (2016), The ASA statement on p-values. ## Answer key **Example B.** 1. $e^{0.5799} = 1.786$. 2. By 1.79; the odds increase by about 79%. 3. $-3.2632 + 2.8995 = -0.3637$. 4. $1/(1 + e^{0.3637}) = 0.410$ (R `predict(..., type = "response")`: 0.4100566). 5. $3.2632 / 0.5799 = 5.63$ hours. 6. OR = 1.79, 95% CI [1.26, 3.04] (profile), *z* = 2.69, *p* = .007. **Exercise C** (R 4.5.1; output abridged and relabelled: the values come from `coef(summary())`, `exp(cbind(coef(), confint()))`, `anova(fit, test = "LRT")` and `predict(type = "response")`, condensed to one line each): ```text Estimate Std. Error z value Pr(>|z|) (Intercept) -2.7269241 1.2355324 -2.207084 0.02730817 abs_dev 0.4542173 0.2005087 2.265325 0.02349277 odds ratio and profile 95% CI: abs_dev 1.57494013 [1.134215917, 2.593523] LR test: Deviance 8.2491 on 1 df, Resid. Dev 19.276, Pr(>Chi) 0.004077 predicted probability at abs_dev 0: 0.0614032 at abs_dev 10: 0.8599950 ``` (a) `temp_dev` assumes the log-odds rise steadily from too cold to too hot, so cold and hot deviations would cancel out; the absolute deviation (or a squared term) matches a U-shaped risk. (b) OR = 1.57 per °C of deviation, 95% CI [1.13, 2.59], LR χ²(1) = 8.25, *p* = .004. (c) 6.1% at the set point, 86.0% at 10 °C off. (d) 9 defects and 11 good batches, so the smaller group is 9 and EPV = 9 / 1 = 9, close to the guide of 10; one predictor is the most these data support. **Self-test.** 1. Odds $= 0.80 / 0.20 = 4$; log-odds $= \ln 4 = 1.386$. 2. Log-odds: each extra weekly hour raises the log-odds of passing by 0.596, review attendance held constant. Odds: each extra hour multiplies the odds of passing by $e^{0.596} = 1.82$. 3. R's `confint` on a glm gives profile likelihood intervals; statsmodels gives Wald intervals. For n = 30, report the profile interval (more accurate when the likelihood is not symmetric) and name the method. 4. Sensitivity $= 14/19 = 0.737$, specificity $= 9/11 = 0.818$. Raising the threshold makes the model predict "pass" less often, so it misses more passers (sensitivity falls from 0.842) and wrongly predicts fewer failers as passers (specificity rises from 0.727). 5. Complete or quasi-complete separation: some combination of predictors predicts the outcome perfectly, so the maximum likelihood estimate goes to infinity. Look for the predictor or category that separates, then use Firth penalised regression (`logistf` in R) or exact logistic regression, or merge sparse categories if that is substantively justified. 6. 24 / 10 = about 2 predictors (counting dummies and interaction terms as separate coefficients). EPV is only a rough guide: for prediction models, use a formal sample size method (Riley and colleagues, 2020), and consider penalised estimation. ## Revision history | Version | Date | Author | Change | |---|---|---|---| | 0.1 | YYYY-MM-DD | [Name] | First draft from the template | | [0.2] | YYYY-MM-DD | [Name] | [Switched to profile CIs and LR tests] |