Logistic Regression — Part 2: Two Predictors
Linguistics / Cognitive Science — Beginner Version (Multifactorial)
K.Osolin
2026-06-22
1 What’s new compared with Part 1
You already know the workflow from the one-predictor template (Part 1). This version changes only one thing: the model now contains two predictors at once.
That has one important consequence for how you read the results:
Each effect is reported while holding the other predictor constant. The effect of Modality below means its effect after accounting for the length of the recipient — and vice versa. Statisticians call this an adjusted effect. This is the whole reason to put predictors in the same model.
Run the blocks in order, top to bottom. STUDENT: marks the few things you may need to change. The real focus of this template is the section “How to read your regression table” — that is where the interpretation happens.
2 Setup and data
# First time only: remove the # and run, then put the # back.
# install.packages(c("tidyverse", "broom", "gtsummary", "ggeffects"))
library(tidyverse)
library(broom)
library(gtsummary)
library(ggeffects)STUDENT: Put your
.csvin the same folder as this file and set its name below.
data <- read_csv("Sem1_bresnan_et_al_2008_dative.csv")
glimpse(data) # <chr> = category variable, <dbl>/<int> = numeric variable## Rows: 3,263
## Columns: 16
## $ Speaker <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Modality <chr> "written", "written", "written", "written", "wr…
## $ Verb <chr> "feed", "give", "give", "give", "offer", "give"…
## $ SemanticClass <chr> "t", "a", "a", "a", "c", "a", "t", "a", "a", "a…
## $ LengthOfRecipient <dbl> 1, 2, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 1, 2, 1, 1,…
## $ AnimacyOfRec <chr> "animate", "animate", "animate", "animate", "an…
## $ DefinOfRec <chr> "definite", "definite", "definite", "definite",…
## $ PronomOfRec <chr> "pronominal", "nonpronominal", "nonpronominal",…
## $ LengthOfTheme <dbl> 14, 3, 13, 5, 3, 4, 4, 1, 11, 2, 3, 3, 5, 2, 8,…
## $ AnimacyOfTheme <chr> "inanimate", "inanimate", "inanimate", "inanima…
## $ DefinOfTheme <chr> "indefinite", "indefinite", "definite", "indefi…
## $ PronomOfTheme <chr> "nonpronominal", "nonpronominal", "nonpronomina…
## $ RealizationOfRecipient <chr> "NP", "NP", "NP", "NP", "NP", "NP", "NP", "NP",…
## $ AccessOfRec <chr> "given", "given", "given", "given", "given", "g…
## $ AccessOfTheme <chr> "new", "new", "new", "new", "new", "new", "new"…
## $ HeavinessNum <dbl> 13, 1, 12, 4, 1, 2, 2, 0, 10, 1, 1, 1, 4, 0, 7,…3 Choose your two predictors
The outcome is always RealizationOfRecipient. You pick
the two predictors.
STUDENT: Edit
selected_vars(exact spelling fromglimpse()) and give each predictor a friendly label invariable_labels. For the clearest teaching example, use one categorical and one numeric predictor (as below). Two categorical predictors also work — just list two categorical column names.
selected_vars <- c("AnimacyOfTheme", "LengthOfTheme")
variable_labels <- c(
AnimacyOfTheme = "Animacy of Theme (inanimate vs animate)",
LengthOfTheme = "Length of Theme (words)"
)
# Keep outcome + chosen predictors, drop rows with missing values,
# and make sure R treats text columns as categorical factors.
data <- data %>%
select(RealizationOfRecipient, all_of(selected_vars)) %>%
drop_na() %>%
mutate(across(where(is.character), as.factor))R always predicts the probability of the second outcome level; the first is the reference (baseline).
# The SECOND name in this output is the outcome being modelled (here: PP).
# Keep this in mind when interpreting every number in the results.
levels(data$RealizationOfRecipient)## [1] "NP" "PP"4 A quick look before modelling
# How many observations per outcome category?
data %>% count(RealizationOfRecipient)# Descriptive summary per predictor:
# - numeric predictor → mean and SD split by outcome
# - categorical predictor → cross-tabulation with outcome
for (var in selected_vars) {
cat("\n=====", var, "=====\n")
if (is.numeric(data[[var]])) {
print(
data %>%
group_by(RealizationOfRecipient) %>%
summarise(Mean = round(mean(.data[[var]]), 2),
SD = round(sd(.data[[var]]), 2),
.groups = "drop")
)
} else {
print(table(data[[var]], data$RealizationOfRecipient))
}
}##
## ===== AnimacyOfTheme =====
##
## NP PP
## animate 22 52
## inanimate 2392 797
##
## ===== LengthOfTheme =====
## # A tibble: 2 × 3
## RealizationOfRecipient Mean SD
## <fct> <dbl> <dbl>
## 1 NP 4.75 4.71
## 2 PP 2.9 2.725 Fit the model and produce the tables
~ . means “predict the outcome from every other column
in data” — i.e. your two predictors.
family = binomial tells R to fit a logistic (not linear)
regression.
model_full <- glm(RealizationOfRecipient ~ ., data = data, family = binomial)# Main results table: odds ratios (OR) with 95% confidence intervals.
# This is the table you will interpret row by row in the next section.
broom::tidy(model_full, exponentiate = TRUE, conf.int = TRUE) %>%
mutate(across(where(is.numeric), ~ round(.x, 3))) %>%
rename(Predictor = term,
OR = estimate,
SE = std.error,
z = statistic,
p = p.value,
CI_low = conf.low,
CI_high = conf.high)# The same results in a clean, report-ready layout.
tbl_regression(model_full, exponentiate = TRUE)| Characteristic | OR | 95% CI | p-value |
|---|---|---|---|
| AnimacyOfTheme | |||
| animate | — | — | |
| inanimate | 0.16 | 0.09, 0.26 | <0.001 |
| LengthOfTheme | 0.85 | 0.82, 0.88 | <0.001 |
| Abbreviations: CI = Confidence Interval, OR = Odds Ratio | |||
6 How to read your regression table
This is the part that matters. Below is an illustrative table (realistic but made-up numbers) of the kind your code just produced. We will read it one row at a time. Your own numbers will differ, but the way you read them is identical.
| Predictor | B | SE | z | p | OR | CI_low | CI_high |
|---|---|---|---|---|---|---|---|
| (Intercept) | −0.95 | 0.21 | −4.52 | <.001 | 0.39 | 0.26 | 0.58 |
| Modality written | 0.71 | 0.18 | 3.94 | <.001 | 2.03 | 1.43 | 2.89 |
| LengthOfRecipient | 0.28 | 0.05 | 5.60 | <.001 | 1.32 | 1.20 | 1.46 |
Before reading any row, fix two things in your mind:
- What is being predicted? From
levels()the second level is PP, so every number in the table is about the odds of PP (versus NP). - OR is the column you interpret. B expresses the same effect on the hard-to-read log-odds scale; the odds ratio (OR) is B made readable. Quick guide: OR = 1 → no effect; OR > 1 → raises odds of PP; OR < 1 → lowers them; (OR − 1) × 100% converts OR to a percentage change in odds.
6.1 Row 1 —
(Intercept): the baseline
The intercept is the model’s starting point: the predicted odds of PP when every categorical predictor is at its reference level (here: Modality = spoken) and every numeric predictor equals zero (Length = 0 words). Its OR of 0.39 means PP is less likely than NP at that particular baseline. You normally do not interpret the intercept substantively — it is just the anchor the other rows are measured from. Move on.
6.2 Row 2 — a
categorical predictor (Modality)
First, the row name: R writes a categorical row as
the variable name joined to the level it describes —
Modalitywritten means “written, compared with the
baseline level spoken.” R never prints a separate row for the
baseline level; it is “hidden” inside the intercept.
Now read it, holding length of recipient constant:
- Direction (OR vs 1). OR = 2.03 > 1, so being written increases the odds of PP relative to spoken.
- Size. The odds of PP are about 2.03 times higher for written than for spoken — i.e. roughly +103% ((2.03 − 1) × 100).
- Is it significant? Two checks that must both pass: the p-value is below .05 (here < .001), and the 95% CI [1.43, 2.89] does not contain 1. Both say: significant. (If a CI straddles 1, the effect is not significant regardless of what the OR looks like.)
In one sentence: holding length constant, written texts have about twice the odds of a PP realization as spoken texts, and this effect is significant.
6.3 Row 3 — a
numeric predictor (LengthOfRecipient)
A numeric predictor has no “levels”, so its OR is read per one-unit increase, holding modality constant:
- Direction. OR = 1.32 > 1 → each extra word increases the odds of PP.
- Size. Every additional word multiplies the odds of PP by 1.32 (about +32% per word). Note that effects are multiplicative: two extra words ≈ 1.32 × 1.32.
- Significant? p < .001 and CI [1.20, 1.46] excludes 1 → yes.
In one sentence: holding modality constant, each additional word of the recipient raises the odds of a PP realization by about 32%, significantly so.
6.4 Putting it together (the “adjusted” idea)
Both predictors are significant at the same time, each adjusted for the other. So the Modality effect is not an artefact of written texts happening to have longer recipients — the model has already removed that overlap, and vice versa. That is exactly what a multifactorial model buys you over two separate one-predictor models.
A 30-second reading recipe for any row (except the intercept): (1) row name → which predictor / which level; (2) OR vs 1 → direction; (3) how far OR is from 1 → size; (4) p < .05 and CI excludes 1 → significance. Always end with “…holding the other predictor constant.”
7 Does each predictor add value beyond the other?
We use nested model comparison (likelihood-ratio / chi-square difference test) to ask whether each predictor improves model fit above and beyond what the other predictor already explains.
The logic is a three-step ladder for each predictor:
null model → model with predictor A only → full model (A + B)- Step 1 (null → A only): does predictor A help at all on its own?
- Step 2 (A only → full): does predictor B add anything on top of A?
We then repeat with the order swapped to evaluate A’s unique contribution beyond B. The p-value in step 2 is the one that answers each research question.
pred1 <- selected_vars[1] # first predictor, e.g. "Modality"
pred2 <- selected_vars[2] # second predictor, e.g. "LengthOfRecipient"
# Build all four models on the same dataset.
model_null <- glm(RealizationOfRecipient ~ 1, data = data, family = binomial)
model_pred1 <- glm(RealizationOfRecipient ~ get(pred1), data = data, family = binomial)
model_pred2 <- glm(RealizationOfRecipient ~ get(pred2), data = data, family = binomial)
# model_full was already fitted in section 5 above.
# ── Test 1: does pred2 add value beyond pred1? ────────────────────────────
# Read the SECOND chi-square row: null→pred1 (does pred1 help on its own?),
# then pred1→full (does pred2 add anything on top of pred1?).
cat("=== Does", pred2, "add value beyond", pred1, "? ===\n")## === Does LengthOfTheme add value beyond AnimacyOfTheme ? ===anova(model_null, model_pred1, model_full, test = "Chisq")# ── Test 2: does pred1 add value beyond pred2? ────────────────────────────
# Same logic, order reversed: read the SECOND row for pred1's unique contribution.
cat("\n=== Does", pred1, "add value beyond", pred2, "? ===\n")##
## === Does AnimacyOfTheme add value beyond LengthOfTheme ? ===anova(model_null, model_pred2, model_full, test = "Chisq")How to read the output of each anova()
call:
The table has two chi-square rows (the null model itself has no test row):
| Row | Comparison | What it tests |
|---|---|---|
| 1 | null → single-predictor model | Does this predictor help at all on its own? |
| 2 | single-predictor → full model | Does the other predictor add unique value? |
Focus on row 2: if p < .05 there, the predictor being added in that step earns its place after the other predictor is already accounted for. This is the adjusted test — it mirrors exactly the logic of the coefficients in the regression table.
APA reporting template for each test (use row 2 values):
“Adding [pred2] to a model already containing [pred1] significantly improved model fit, χ²([Δdf]) = [Δdeviance], p = [p], indicating that [pred2] explains variance in recipient realization beyond what [pred1] alone accounts for.”
If p ≥ .05: “Adding [pred2] to a model already containing [pred1] did not significantly improve model fit, χ²([Δdf]) = [Δdeviance], p = [p].”
8 Effect plots
One plot per predictor: the predicted probability of PP across that predictor’s range, holding the other predictor constant at its mean (numeric) or most frequent level (categorical). Shaded bands / error bars are 95% CIs.
for (var in selected_vars) {
eff <- ggpredict(model_full, terms = var)
label <- if (var %in% names(variable_labels)) variable_labels[[var]] else var
print(plot(eff) + ggtitle(paste("Predicted probability of PP by", label)))
}
9 Ready-made APA sentences
Write one sentence per predictor. Copy B, SE, z, p and the OR with its CI directly from the odds-ratio table in section 5. Round to 2 decimal places; report p exactly or as p < .001.
Overall model fit (use values from
anova(model_null, model_full, test = "Chisq"))
A logistic regression predicted recipient realization (NP vs. PP) from [predictor 1] and [predictor 2]. The full model significantly improved on a null model, χ²([Δdf]) = [Δdeviance], p [= / < .001].
Significant categorical predictor
Holding [other predictor] constant, [Predictor] significantly predicted recipient realization, B = [B], SE = [SE], z = [z], p = [p]; the odds of a PP realization were [OR] times [higher / lower] for [level] than for [baseline], OR = [OR], 95% CI [[CI_low], [CI_high]].
Significant numeric predictor
Holding [other predictor] constant, [Predictor] significantly predicted recipient realization, B = [B], SE = [SE], z = [z], p = [p]; each additional [unit] [increased / decreased] the odds of a PP realization by a factor of [OR], OR = [OR], 95% CI [[CI_low], [CI_high]].
Non-significant predictor
Once [other predictor] was accounted for, [Predictor] did not significantly predict recipient realization, B = [B], SE = [SE], z = [z], p = [p].
Nested model comparison (use values from row 2 of
the relevant anova() output)
Adding [pred2] to a model already containing [pred1] significantly improved model fit, χ²([Δdf]) = [Δdeviance], p = [p], indicating that [pred2] explains unique variance in recipient realization beyond [pred1].