Class 8 · Multiple Logistic Regression and Lab 7.6

STAT 517: Advanced Statistical Models · Fall 2026

This interactive companion follows the Class 8 notes and lecture script: multiple logistic regression, adjusted odds ratios, inference, model selection, and Lab 7.6: What changes when we adjust?

Each code box runs base R directly in your browser. Press Run Code, inspect the output, edit the code, and run it again. Run the lecture cells in order because later cells reuse earlier objects. The two data files are included; no installation or upload is needed. Profile-likelihood intervals and stepwise selection can take a little longer to run.

For the group activity, you can go directly to Section 8. Its embedded data make the lab independent of the lecture cells. Set group to your assigned number (1–6), then run the lab cells in order.

Goals

• Define the modeled event and reference categories.
• Interpret adjusted coefficients and odds ratios, holding other predictors fixed.
• Distinguish the overall likelihood-ratio test from a test of one coefficient.
• Construct confidence intervals for odds ratios and event probabilities.
• Compare models on the same observations using AIC.
• Explain what changes when correlated predictors enter a model.

1. Framingham: prepare the data and fit the model

The response is the ten-year coronary heart disease outcome. With a two-level factor response, binomial() models the probability of the second level: here, AtRisk. The five binary predictors use 0 as their reference level. A numerical 0/1 predictor gives the same slope as a two-level factor with reference 0; factor coding makes the categories explicit.

Use the same 3,800 complete observations for every Framingham model below. Otherwise differences between fits could also reflect differences in samples. Complete-case analysis can introduce bias if the missingness mechanism is incompatible with that analysis.

For predictors \(x_1,\ldots,x_p\), the model is

\[ \log\left\{\frac{p(\mathbf{x})}{1-p(\mathbf{x})}\right\} =\beta_0+\beta_1x_1+\cdots+\beta_px_p. \]

2. Fitted probabilities and diagnostic displays

type = "link" returns fitted log-odds; type = "response" returns fitted event probabilities. plogis() transforms log-odds to probabilities stably.

Residuals and influence

Binary responses naturally produce two bands in residual displays. Normal errors and constant error variance are not logistic-regression assumptions, so the usual normal Q–Q interpretation from linear regression does not apply. Cook’s distance can still help identify observations worth investigating.

Explore fitted log-odds against continuous predictors

The lecture overlays smooth curves on fitted log-odds. These are descriptive plots: the fitted log-odds already come from the assumed model, and correlations with other predictors affect each marginal curve. They cannot establish that the true conditional logit is linear. A more direct assessment would compare prespecified nonlinear terms against the linear specification.

3. Inference and adjusted odds ratios

Test all slopes together

The overall likelihood-ratio test compares the fitted model with an intercept-only model: \(H_0:\beta_1=\cdots=\beta_p=0\). The statistic is the null deviance minus the residual deviance. Its degrees of freedom count slope coefficients, including each factor contrast. Rejection indicates improvement over the intercept-only model; it does not establish adequate absolute fit or predictive performance.

Test individual coefficients

Each Wald test concerns one coefficient conditional on the other predictors in the model. A large p-value does not prove absence of association.

Marginal and conditional associations

Compare the BMI and heart-rate slopes from single-predictor models with their slopes in the full model. The comparison uses identical observations.

Predictors can contain overlapping information. The following summaries explore that overlap; they do not by themselves establish why a particular coefficient has a large p-value. In logistic regression, marginal and conditional odds ratios can also differ because odds ratios are noncollapsible.

Odds ratios and confidence intervals

For a one-unit increase in \(x_j\), \(\exp(\hat\beta_j)\) is the adjusted odds multiplier. For a factor contrast, it compares that category with its reference. The percentage change in odds is \(100\{\exp(\hat\beta_j)-1\}\); this is not a percentage change in probability or risk.

Interpret the male and age odds ratios while holding the other covariates fixed. Profile-likelihood intervals and Wald tests use different approximations, so their conclusions need not agree exactly near a cutoff.

4. Model selection with AIC

AIC is \(-2\ell(\hat\beta)+2k\), with \(k\) including the intercept. Smaller AIC indicates a better balance of fitted likelihood and parameter count among models for the same response and observations. It does not mean that every retained coefficient must pass a 5% Wald test.

Stepwise selection is an exploratory procedure. Ordinary p-values and intervals computed after selection do not account for the search over models.

5. Horseshoe crabs: multiple predictors

The event is having at least one satellite male, satell > 0. Width is in centimeters and the original weight is in grams; dividing by 1,000 gives kilograms. Color and spine condition are factors, with LightMedium and BothGood as the reference categories.

The formula lists predictors explicitly so the satellite count used to define the response is not inadvertently included as a predictor.

6. Crab model selection and interpretation

Interpret the width odds ratio and its confidence interval after adjusting for color. Individual color contrasts compare one level with LightMedium; none of those individual tests is a joint test of the entire color factor. The following likelihood-ratio comparison removes all color contrasts together from the selected model. Treat this post-selection comparison as exploratory.

7. Worked Examples 7.2–7.5

These hypothetical examples isolate the calculations used in the notes.

Example 7.2

Example 7.3

Example 7.4

Example 7.5

8. Lab 7.6: What changes when we adjust?

Set group to your assigned group (1–6). Run this section in order. Every group uses the same 173 crabs, the same event, and the same prediction profile; only the adjustment variables change. The data are embedded so this lab section can run independently of the lecture examples. Only the assigned model is fitted by default.

Group	Predictors
1	Width
2	Width + color
3	Width + weight
4	Width + spine
5	Width + weight + color
6	Width + weight + color + spine

8.0 Lab setup: hold the data and event fixed

Original 2024 crab data, embedded so this file needs no download or package. Source: Sources_2024/R code for lectures/crabs.txt Raw columns: color, spine, width (cm), satell (count), weight (grams). The multiline text below is the entire dataset. read.table(…, header=TRUE) uses its first line as column names and reads the remaining lines as rows. Do not filter different rows for different groups: a comparison of models must not quietly become a comparison of different samples as well.

The comparison satell>=1 is TRUE or FALSE for each crab. as.integer turns these into 1 and 0, so glm models P(at least one satellite male | covariates). Unlike Lab 6.2, EVERY group now uses this same event definition.

Dividing the original grams by 1000 changes the unit to kilograms. Keep width in centimeters. Units determine what a one-unit OR means.

factor() makes separate indicators rather than treating the labels as a numerical trend. The first level is the reference: LightMedium for color, BothGood for spine. A four-level factor contributes three slopes; a three-level factor contributes two. Keep these references the same for every group so that coefficient comparisons use the same definitions.

Each model has an intercept and additive predictors, without interactions. A formula lists the response to the left of ~ and the included predictors to the right. The + signs add predictors; they do not add interaction terms. The list stores six FORMULAS. It does not fit all six models. [[group]] below selects exactly one formula. Group 1 supplies the unadjusted width baseline.

Every group predicts for this same hypothetical crab. A model uses only the columns named in its formula; supplying the other profile values does not add predictors to that model. Factor labels and levels must agree with the training data, so we copy levels(crabs\(color) and levels(crabs\)spine).

8.1 Lab Part 1: fit and interpret width (7 minutes)

stopifnot checks that the group is a single integer between 1 and 6. binomial() uses the logit link. glm estimates the coefficients by maximum likelihood; the convergence check confirms that its iteration terminated.

Width: change in log odds for one additional cm, holding the other covariates in this group’s model fixed. Group 1 is unadjusted. coef(…)[“width”] selects the coefficient by NAME so the code still works when different groups have different numbers and orders of predictors. The diagonal of vcov() contains estimated coefficient variances. Taking the square root gives the standard error of the width coefficient.

For a one-cm increase, the adjusted odds multiplier is exp(width_beta). First form beta_hat +/- 1.96*SE on the log-odds scale, then exponentiate both endpoints. This preserves their order because exp() is increasing. Under H0: beta_width=0, beta_hat/SE is approximately standard normal. pnorm(-abs(…)) gives one tail; multiplying by 2 gives the two-sided test.

WRITE: “Holding [your model’s other covariates] fixed, an additional cm of width is associated with [OR] times the odds of at least one satellite.” State whether the CI includes 1 and explain the conditional conclusion. A multiplier of 1 would mean equal odds. An OR is not a multiplier for the probability itself, and a nonsignificant test does not prove no association.

8.2 Lab Part 2: test the model as a whole (4 minutes)

Overall LR test: all slope coefficients are zero versus at least one nonzero slope. Factor levels contribute separate slope coefficients. The null deviance comes from an intercept-only fit to these observations; the residual deviance comes from your assigned model. Their difference measures the improvement in fit from including that model’s slopes. rank counts all fitted coefficients, including the intercept. Subtract 1 to get the number of slopes. Do not count a multi-level factor as one slope. lower.tail=FALSE requests the upper chi-squared tail directly. It keeps very small p-values accurate instead of subtracting a number near 1 from 1.

WRITE: your overall null hypothesis, LR statistic, df, p-value, and decision at 0.05. Compare this ALL-slopes null with the ONE-width-coefficient null from Part 1. Rejection of one hypothesis does not force rejection of the other.

8.3 Lab Part 3: predict a mean probability (6 minutes)

Common profile: width 26.3 cm, weight 2.5 kg, Medium color, BothGood spine. Form an interval on the link scale and apply the inverse logit to both endpoints. This is a CI for a mean probability, not a binary outcome. type=“link” returns eta_hat = x_new’ beta_hat rather than a probability. se.fit=TRUE also returns the standard error of that fitted linear predictor. This standard error uses the whole coefficient covariance matrix; it is not the width coefficient’s standard error from Part 1. plogis(eta) = exp(eta)/(1+exp(eta)). Applying it to eta_hat and to the two logit-interval endpoints produces probabilities strictly between 0 and 1.

WRITE: the estimated probability and its 95% confidence interval. The interval concerns the mean event probability at the covariate values your model includes; it does not mean a crab’s observed outcome lies between those endpoints. That individual outcome must be either 0 or 1.

8.4 Report your results and compare groups (5 minutes)

These three tables supply the numbers for your group’s GitHub comment. Use four decimals for estimates/CIs; keep scientific notation for tiny p-values so a value such as 2e-8 is not misreported as exactly zero.

DISCUSS with the class:

1. Compare with Group 1 (Group 1 compares with an adjusted model). What changes in the width estimate and its SE after weight is included? How might overlap between the two body-size measurements explain this?
2. Can the overall test reject while the conditional width test does not? Use the two different null hypotheses, not just the two p-values.
3. Does a larger width OR guarantee a larger probability for this profile? The intercept, other coefficients, and covariate values also enter eta.

Report fitted associations; these observational comparisons do not by themselves establish what would happen if width were changed by intervention.

8.5 Optional extension: AIC

Optional after the group comparison: lower AIC balances fit and model size. This fits only a width baseline in addition to your assigned model. Both models use the same event and observations. AIC=-2logLik+2k, with k including the intercept. Adding predictors can improve likelihood yet increase AIC if the improvement is smaller than the added parameter penalty.