--- title: "Categorical predictors" description: > How categorical predictors enter a regression table and how to work with them: dummy coding and the reference level, joint (multiple-degree-of-freedom) tests of a factor, R's polynomial contrasts for ordered factors, the scores-versus-dummies decision for ordinal predictors, successive-difference coding, and why continuous predictors should not be categorized. output: rmarkdown::html_vignette: toc: true toc_depth: 2 vignette: > %\VignetteIndexEntry{Categorical predictors} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) source("_pkgdown-helpers.R") ``` ```{r setup} library(spicy) ``` Every regression vignette in this collection is organised by the type of *outcome* — continuous, binary, count, ordinal, multinomial, survival. This one is organised by the type of *predictor*. How a categorical predictor is coded, tested, simplified, and presented is the same decision whether the model is an `lm()`, a logistic regression, or a Cox model, so it has one home rather than a copy in every family vignette. The examples use the bundled `sochealth` survey data (`?sochealth`) with `lm()` and `glm()` fits; everything shown carries over to the other supported classes. The running order follows the analysis itself: how a factor enters the model, how to choose what it is compared against, how to test whether it matters at all, what changes when its levels are ordered, and — the reverse question — why a continuous predictor should not be turned into a factor. ## How a factor enters the model R expands a factor with `k` levels into `k − 1` **indicator (dummy) variables**, one per non-reference level; the remaining level becomes the **reference** absorbed by the intercept (Gelman, Hill & Vehtari 2020, §10.4). `table_regression()` groups the indicator rows under the parent variable and renders the reference level explicitly, with `(ref.)` and an en dash in the statistic columns: ```{r basic} fit <- lm(wellbeing_score ~ age + employment_status, data = sochealth) table_regression(fit) ``` Each coefficient is a comparison **with the reference level, holding the other predictors constant** — not an overall effect of the variable. Reading the rows against that reference: unemployed respondents score on average 3.6 points lower on the 0–100 wellbeing index than *employed* respondents of the same age (B = -3.60, 95% CI [-6.17, -1.03], p = .006); for students and inactive respondents the data show no clear difference *from the employed* (p = .441 and .854). Nothing in this table compares students with the unemployed — that contrast is the difference between the two coefficients, and its standard error is not printed here. If one specific non-reference contrast is the research question, the cleanest route is to make one of its ends the reference. ## Choosing and changing the reference level The reference level is a **presentation choice, not a modeling one**: refitting with another reference is the same model reparameterized — identical fit, likelihood, AIC, and fitted values; only the comparisons displayed change (Gelman, Hill & Vehtari 2020, §10.4). ```{r relevel} d2 <- sochealth d2$employment_status <- relevel(d2$employment_status, ref = "Unemployed") fit_u <- lm(wellbeing_score ~ age + employment_status, data = d2) c(AIC(fit), AIC(fit_u)) table_regression(fit_u, keep = "employment_status") ``` Every group now sits 3.6 to 4.7 points *above* the unemployed — the same fit read from the other end. Two practical rules: * **Do not let alphabetical order decide.** `factor()` sorts levels alphabetically by default, so the reference is whichever category happens to sort first — rarely the one you would have chosen, and a silent source of confusing signs when data come from text files. * **Choose a reference that makes every displayed contrast readable**: the largest category (stable comparisons), an unexposed / control category (epidemiology), or the substantively natural baseline. A rare category makes a poor reference — every contrast inherits its noise. ## Does the factor matter at all? The joint test The per-level p-values above cannot answer "does employment status matter?": each one tests a single level against the reference, and the set of them changes with the reference choice. The question about the factor *as a whole* is a **joint test of all `k − 1` coefficients at once** — a multiple-degree-of-freedom test, invariant to the reference (Harrell 2015, §2.6). Harrell's rule is to run it *before* interpreting any individual contrast: when levels straddle significance, as here, the joint test is the honest headline. The change statistics of a nested comparison give exactly that test — adding the factor to a model that lacks it. One prerequisite: change statistics are only valid when both models are fit to the same rows, and R's listwise deletion silently drops *different* rows when the added variable carries its own missing values — so the code below first restricts the data to cases complete on every variable involved. `table_regression()` enforces this: models with different `nobs` are refused, with instructions to refit on the common subset. ```{r joint} d_cc <- sochealth[complete.cases( sochealth[, c("wellbeing_score", "age", "employment_status")] ), ] m0 <- lm(wellbeing_score ~ age, data = d_cc) m1 <- lm(wellbeing_score ~ age + employment_status, data = d_cc) table_regression(list(m0, m1), nested = TRUE, show_columns = c("b", "p")) ``` The `F-change` row is the joint 3-df test: F = 3.13, p = .025 — employment status improves the model, driven by the unemployed contrast. For `glm` fits, the `partial_chi2` column reports the same joint question per term as a likelihood-ratio chi-square with its df, without needing a second model (see the main vignette's *Term-level partial chi-square*). The rule cuts the other way too: ```{r joint-region} d_r <- sochealth[complete.cases( sochealth[, c("wellbeing_score", "age", "region")] ), ] anova(lm(wellbeing_score ~ age, data = d_r), lm(wellbeing_score ~ age + region, data = d_r)) ``` For `region`, the joint test gives F(5, 1193) = 0.89, p = .48 — no evidence that region matters at all (`anova()` on two nested fits is the same F-change test the table reports). Fitting the factor anyway and scanning its five per-level p-values for a "significant" contrast to report alone would be a multiplicity artifact — "if the overall test is not significant, it can be dangerous to rely on individual pairwise comparisons" (Harrell 2015, §2.6). ## Ordered factors: R's polynomial contrasts An **ordinal** predictor — levels with a natural order, like `self_rated_health` (Poor < Fair < Good < Very good) — raises a coding decision that R quietly makes for you. For an `ordered` factor, R does *not* use dummy coding: it uses **orthogonal polynomial contrasts**, and the table shows trend components instead of level comparisons. spicy flags this the first time it happens and the footer names the coding: ```{r ordered} fit_srh <- lm(wellbeing_score ~ age + self_rated_health, data = sochealth) table_regression(fit_srh, show_columns = c("b", "se", "p")) ``` These rows answer trend questions across the ordered levels: `.L` (linear) is the dominant component here — wellbeing rises steeply across health ratings — while `.Q` (quadratic, B = -3.93, p < .001) says the rise is **concave**: the steps get smaller toward the top of the scale. `.C` (cubic, p = .235) adds nothing. This decomposition is a legitimate reading — trend tests are exactly what polynomial contrasts are for — but it is *not* the level-versus-reference reading most regression tables intend. For that, convert to a plain factor (`factor(x, ordered = FALSE)`), as the message suggests and as the other vignettes in this collection do. ## Scores or dummies? An ordinal predictor poses a substantive decision: should it enter as **numeric scores** (1, 2, 3, … — one slope, read as "per step up the scale") or as **dummies** (k − 1 coefficients, no order imposed)? The scores model is nested in the dummy model, so a likelihood-ratio test compares them: fit both codings and let the data say whether the level coefficients depart from a linear trend in the scores (Agresti 2007, §4.4.3). `nested = TRUE` runs the comparison. The two examples below land on opposite sides of the decision. **Case 1 — the scores suffice.** Education (3 levels) predicting smoking: ```{r scores-education} sh <- sochealth sh$education <- factor(sh$education, ordered = FALSE) sh$educ_score <- as.numeric(sh$education) sh <- sh[complete.cases(sh[, c("smoking", "sex", "age", "education")]), ] g_lin <- glm(smoking ~ sex + age + educ_score, data = sh, family = binomial) g_fac <- glm(smoking ~ sex + age + education, data = sh, family = binomial) table_regression(list(g_lin, g_fac), nested = TRUE, show_columns = c("b", "p")) ``` Freeing education from the linear trend buys a chi-square of 0.03 on 1 df, p = .852 — no evidence the dummies improve on the scores — and AIC agrees (1198.9 vs 1200.9). The scores model is simpler to report (one coefficient: B = -0.45 per education step, on the log-odds scale) and more powerful against a genuine trend, because it spends one parameter where the dummies spend two. With 3 levels the test has a familiar counterpart: it is the same 1-df question as the Wald test on the `.Q` contrast of the ordered coding (z-test p = .852 here; in a GLM the Wald and likelihood-ratio answers agree asymptotically, and exactly in a linear model) — the quadratic component *is* the departure from linearity. **Case 2 — linearity is rejected.** Self-rated health (4 levels) predicting wellbeing: ```{r scores-srh} d3 <- sochealth[complete.cases( sochealth[, c("wellbeing_score", "age", "self_rated_health")] ), ] d3$srh_score <- as.numeric(d3$self_rated_health) d3$srh_factor <- factor(d3$self_rated_health, ordered = FALSE) s_lin <- lm(wellbeing_score ~ age + srh_score, data = d3) s_fac <- lm(wellbeing_score ~ age + srh_factor, data = d3) table_regression(list(s_lin, s_fac), nested = TRUE, show_columns = c("b", "p")) ``` Here the test rejects the single slope: F-change = 12.96 on 2 df, p < .001. The dummy coefficients show why: the steps from Poor are +14.97 (Fair), +28.01 (Good), +35.13 (Very good) — increments of about 15, 13, then 7 points. The effect flattens at the top of the scale, exactly the concavity the `.Q` contrast detected above (the 2-df test here is identical to the joint test of `.Q` and `.C`). Forcing one slope of 10.96 per step would overstate the last step and understate the first. Keep the dummies. Three caveats keep the decision honest (Harrell 2015, §§2.6 and 2.7.3; Grambsch & O'Brien 1991): * **A non-significant test does not prove linearity** — in modest samples several codings fit adequately, and integer scores assume the steps are *equally spaced*, a substantive claim about the scale, not a statistical default. * **Test-then-simplify inflates error.** Reporting the scores model as if it had been chosen in advance makes the subsequent test's distribution wrong. Grambsch and O'Brien showed, for the quadratic-versus-linear case, that the pre-specified 2-df joint test remains close to optimal *even when the truth is linear*; the same logic carries over to the factor's full-df test — so when the factor is a control variable rather than the focus, the dummies cost little and assume nothing. * **The stakes rise with the number of levels.** With a handful of levels — three to five, in Harrell's account — the variable is usually modeled as a polytomous factor: the dummies spend a parameter or two and assume nothing. The advantage of scores grows when `k` is larger or the model is dense relative to the sample. ## A third coding: successive differences Dummy coding answers "how far is each level from the reference"; scores answer "how much per step, forced equal". A third question — "how big is **each** step?" — has its own coding: MASS's successive-difference contrasts, whose coefficients are the differences between *adjacent* levels (Venables & Ripley 2002). Harrell (2015, §2.7.3) reaches the same reading with threshold indicators `[X ≥ level]`; the two parameterizations give identical fits. ```{r sdif} d3$srh_sdif <- d3$self_rated_health contrasts(d3$srh_sdif) <- MASS::contr.sdif(4) table_regression( lm(wellbeing_score ~ age + srh_sdif, data = d3), keep = "srh_sdif" ) ``` Each row is one step of the scale: Fair sits 14.97 points above Poor, Good 13.04 above Fair, Very good 7.12 above Good — the flattening gradient read directly, with a standard error per step. This is the natural display when the steps themselves are the finding. (Two notes: the rows group under the factor but carry no `(ref.)` row — no reference level exists under this coding; and the joint test of the factor is unchanged — same model, different questions asked of it.) ## The reverse temptation: categorizing a continuous predictor Turning a *categorical* variable into a number is sometimes justified (the scores above); turning a **continuous** predictor into categories almost never is (Harrell 2015, §2.4.1). `sochealth` ships both `bmi` and a conventional `bmi_category` (normal weight / overweight / obesity), so the cost can be shown directly: ```{r categorize} d4 <- sochealth[complete.cases( sochealth[, c("wellbeing_score", "bmi", "bmi_category")] ), ] d4$bmi_category <- factor(d4$bmi_category, ordered = FALSE) # per-level dummies, as above b_cont <- lm(wellbeing_score ~ bmi, data = d4) b_cat <- lm(wellbeing_score ~ bmi_category, data = d4) table_regression(list(b_cont, b_cat), show_columns = c("b", "p"), show_fit_stats = c("nobs", "r2", "aic")) ``` The continuous fit explains twice the variance with half the parameters (R² = .020 with 1 df against .010 with 2 df; AIC 9884.2 vs 9898.1). That is Harrell's argument in one table: categorization discards the variation *within* intervals — every BMI from 25.0 to 29.9 is treated as the same exposure — and assumes wellbeing jumps at the cut-points and is flat between them, a far stronger assumption than linearity. The estimates also stop being portable: the "Overweight vs Normal weight" coefficient (B = -1.69, p = .082) depends on where the sample sits inside the 25.0–29.9 interval, where the per-unit slope (B = -0.60 points per BMI unit, p < .001) does not. And if a *nonlinear* effect is the worry, the remedy is a flexible continuous fit — a polynomial or a spline — not bins: ```{r quadratic} AIC(b_cont, update(b_cont, . ~ . + I(bmi^2))) ``` The quadratic term adds nothing (AIC 9885.9 against 9884.2 for the linear fit), so the linear fit stands. Cut-points chosen *after looking at the outcome* do more damage still: p-values and CIs for "optimal" cut-points are invalid without special corrections, and such cut-points fail to replicate across studies (Harrell 2015, §2.4.1). Categories remain fine for *descriptive* tables and for presentation-layer grouping — the mistake is putting them in the model when the continuous measurement exists. ## Presentation The presentation layer is independent of the coding decisions above. `labels` renames variables and levels; `reference_style = "annotation"` lifts the reference into the factor header for compact tables; `keep` / `drop` filter rows by regex after all statistics are computed (see the main vignette's *Display options* and *Filtering displayed coefficients*): ```{r presentation} table_regression( fit, labels = c(age = "Age (years)", employment_status = "Employment status"), reference_style = "annotation" ) ``` ## See also - `vignette("table-regression", package = "spicy")` for the shared table mechanics (variance estimators, nested layouts, filtering, output formats). - `vignette("table-regression-ordinal", package = "spicy")` for ordinal *outcomes* — the proportional-odds model, its tests and relaxations. - `vignette("table-regression-multinomial", package = "spicy")` for the scores-vs-dummies decision inside a multinomial model, where each freed coding costs one parameter per equation. ## References * Agresti, A. (2007). *An Introduction to Categorical Data Analysis* (2nd ed.). Wiley. * Gelman, A., Hill, J., & Vehtari, A. (2020). *Regression and Other Stories*. Cambridge University Press. * Grambsch, P. M., & O'Brien, P. C. (1991). The effects of transformations and preliminary tests for non-linearity in regression. *Statistics in Medicine*, 10(5), 697–709. * Harrell, F. E. (2015). *Regression Modeling Strategies* (2nd ed.). Springer. * Venables, W. N., & Ripley, B. D. (2002). *Modern Applied Statistics with S* (4th ed.). Springer.