Atlas
statminds
GLM (Polytomous Logit Model)The underlying model family class (e.g. GLM, linear model, categorical matrix, log-linear).Parametric ReferenceStatistical methods that assume a specific probability distribution family (typically normal).12-stage workflow

Multinomial Logistic Regression

The engine for Multi-Category Discovery. Multinomial Logistic Regression audits the likelihood of membership in three or more unordered groups, revealing the predictors that drive categorical choice.

Model familyGLM (Polytomous Logit Model)
Hypothesistwo-tailed
AliasesMultinomial Logit · Polytomous Regression · Multi-Group Logistic Model
G1
Nominal Choice Audit
Determine the specific predictors that move participants between unordered categories (e.g., Treatment A vs. B vs. C).
G2
Relative Odds Mapping
Quantify the likelihood of choosing a specific category compared to a fixed 'Baseline' or reference group.
G3
Categorical Signal Isolation
Isolate the unique multivariable influence on discrete, non-ranked outcomes.
1

What is it?

Multinomial Logistic Regression models the probability of nominal multi-category outcomes (e.g. choice of Car, Bus, or Train) relative to a baseline reference category.

2

When to use it

  • Nominal Outcomes: Categorical choices containing no inherent ordinal sequence.
  • Multiple Probabilities: Fit probability curves summing to 1.0 across all classes.
Interactive Sandbox

Multinomial Choice Live Laboratory

Adjust preference slopes for Choice A and Choice B (relative to Baseline Choice C) to see probability curve dynamics.

Presets
Choice A Preference Slope1.0
Choice B Preference Slope-0.5
Probability space (Choice A: blue; Choice B: green; Baseline C: amber)Sum of all choice probabilities at any point = 100%
Choice AChoice BBaseline C
The 12-Stage Precision Workflow
01Category Logits
Hypotheses
We test the null of zero influence against the discovery of a non-zero shift in the odds of group membership relative to the reference.
02IIA Mandate
Assumptions
The ultimate gatekeeper: the Independence of Irrelevant Alternatives (IIA). Adding or removing a category must not alter the relative odds between others.
03Reference Calibration
Diagnostics
Carefully selecting the 'Reference Category'—the anchor upon which all relative odds are mathematically built.
04focus
Predicting patient preference for FlowMotion vs. Standard PT vs. Self-Directed Care based on Age and Baseline Pain.
05Ordinal Pivot
Alternatives
Knowing when to switch to Proportional Odds if the categories have a natural ranking (e.g., Low vs. Med vs. High).
06Likelihood Ratio
Significance
Executing the Likelihood Ratio Test to determine if a predictor adds global value across the entire categorical vector.
07Relative Risk Ratio
Effect Size
Interpreting RRR: e.g., 'For every unit of X, participants are 1.5 times more likely to be in Group B than the Reference Group.'
08Cell Density
Sample Size
Ensuring a minimum of 15-20 observations per predictor per category to prevent the model from failing to converge.
09The Contrast Table
Reporting
Reporting multiple sets of coefficients—one for every category compared to the reference—to tell the full multi-dimensional story.
10nnet / mlogit Logic
Software
Executing the 'multinom' or 'mlogit' commands, ensuring the outcome factor levels are correctly ordered and referenced.
11focus
The fatal error of misinterpreting the odds because the user forgot which group was the mathematical 'Anchor'.
12focus
Tracing the model back to the choice-modeling frameworks established by McFadden and the expansion of GLMs for nominal data.
01Hypothesis test logic

Hypotheses

Pragmatic null and alternative hypotheses defined in mathematical notation.

We ask not just 'is there a link?', but 'how much does Y change for every unit of X?'
Logic Core
Null · H₀

H₀: β₁ = 0 (predictor has no effect on log-odds of outcome category)

Alternative · Hₐ

Hₐ: β₁ ≠ 0 (predictor affects log-odds)

Why it matters two-tailed

For each predictor and outcome category pair (relative to reference category). Overall model test: likelihood ratio test compares fitted model to intercept-only model. Coefficients are log-odds of category vs. reference; exponentiate for relative risk ratios (RRR) or odds ratios.

02Model diagnostics

Assumptions

The core mathematical criteria needed to ensure that statistical testing remains unbiased and valid.

Linearity is a strong claim. Nature often curves; ensure your model does not force a straight line on a bent world.
Integrity Shield
7
Assumptions
5
Critical / High Severity
How to check
Quick
Inspect DV; verify ≥3 categories (e.g., occupation, transportation mode, political party). Ensure categories are TRULY unordered: no low→medium→high progression. Check category frequencies: each category should have sufficient cases (aim for ≥30 per category)
Rigorous
table(outcome) in R or value_counts() in Python; ensure ≥3 levels. Confirm no implicit ordering: if ordered (e.g., 'low', 'medium', 'high' satisfaction), ordinal logistic is more appropriate and powerful. Check theoretical justification for treating as nominal. Examine distribution across categories (imbalanced categories can cause problems)
If violated
If outcome has inherent ordering (e.g., educational attainment: high school, college, graduate degree): use ordinal logistic regression (proportional odds model, cumulative logit)—more powerful and parsimonious. If only 2 categories: use binary logistic regression (simpler, easier interpretation). If categories represent counts: use Poisson or negative binomial regression. If continuous outcome: use linear regression. If hierarchical/nested categories: use nested logit model
How to check
Quick
Design review; check for subject IDs appearing multiple times, nested data (students in schools), matched sets, or longitudinal structure. Verify each observation represents a different, independent case
Rigorous
Examine study design documentation; test for clustering via intraclass correlation if hierarchical structure suspected; check for autocorrelation if observations ordered temporally or spatially. Look for patterns in residuals by cluster/group
If violated
If repeated measures (same subjects at multiple times): use mixed-effects multinomial logistic (random intercepts/slopes) or GEE for multinomial outcomes. If clustered data (e.g., patients within clinics): use cluster-robust standard errors or random effects models. If matched sets (case-control): use conditional logistic regression. Document dependence structure and account for it. Ignoring dependence severely underestimates SEs and inflates Type I error
gee
How to check
Quick
IIA assumption: odds between any two outcome categories should not depend on presence/absence of other categories. Example violation: In transportation mode choice (car, bus, red bus), adding 'red bus' shouldn't affect car:bus ratio—but it might if red bus and bus are similar. Hausman test: compare full model to model excluding one category; significant difference suggests IIA violation
Rigorous
Hausman-McFadden test: fit model on full data vs. restricted sample (drop one outcome category); if coefficients differ significantly, IIA violated. Small-Hsiao test (more robust). Test with multiple excluded categories. Examine if outcome categories have similar characteristics (similarity suggests potential IIA violation). Use nested logit if hierarchical structure exists
If violated
If IIA violated (often due to similar/overlapping categories): (1) Combine similar categories into one (e.g., 'bus' and 'red bus' → 'bus'); (2) Use nested logit model if categories have hierarchical structure (e.g., transportation: private [car, motorcycle] vs. public [bus, train]); (3) Use mixed logit (random parameters) which relaxes IIA; (4) Use alternative-specific models with additional category-level predictors. Report Hausman test results. IIA rarely perfectly holds; minor violations may be acceptable if no better alternative
How to check
Quick
For each continuous predictor: plot empirical log-odds vs. predictor for each outcome category (vs. reference). Categorize predictor into quantiles, calculate log-odds per quantile, check linearity. Look for curvature, U-shapes, or thresholds
Rigorous
Box-Tidwell test: add interaction between continuous predictor and its log; if significant, linearity violated (do separately for each outcome category). Component-residual plots for each predictor×category combination. Compare model fit (AIC, BIC) with/without polynomial terms. Lowess smoothed plots of residuals
If violated
Add polynomial terms (X², X³) for curvilinear relationships. Use restricted cubic splines for flexible, data-driven curves. Transform predictor (log, sqrt) if relationship monotonic but curved and theoretically justified. Use generalized additive model (GAM) with multinomial family for fully non-parametric smoothing. Categorize continuous predictor as last resort (loses information/power). Test transformations for each outcome category separately; may need different forms for different categories
How to check
Quick
VIF from auxiliary linear regressions (regress each predictor on all others; VIF = 1/(1-R²)). VIF <5 ideal, <10 acceptable, ≥10 problematic. Correlation matrix of continuous predictors (|r| > .90 concerning). Check for convergence warnings or inflated SEs
Rigorous
Condition index (eigenvalue-based diagnostics): condition index >30 indicates severe multicollinearity. Examine standard errors: extremely large SEs suggest multicollinearity. Compare coefficients when predictors added/removed; large changes signal collinearity. Test stability of estimates via bootstrap
If violated
Remove one of highly correlated predictors (keep theoretically most important). Combine correlated IVs into composite score or index. Use principal components analysis (PCA) to create orthogonal predictors. Use penalized multinomial regression if available (ridge, lasso, elastic net—less common for multinomial). Report VIF and justify variable selection. Multicollinearity inflates SEs but doesn't bias estimates (unless perfect collinearity)
How to check
Quick
Count cases per outcome category: each should have ≥30 cases (absolute minimum), ≥50 preferred. Rule of thumb: need n_min × p per category, where n_min=30-50 and p=number of predictors. Example: 3 predictors, 4 outcome categories → need ≥30×3=90 cases per category = 360 total. Check smallest category count
Rigorous
Calculate events per variable (EPV) for smallest category: EPV = n_smallest_category / (p × [J-1]), where J=number of outcome categories. Aim for EPV≥10. Simulate power for your specific scenario. Bootstrap to assess estimate stability (wide bootstrap CIs suggest insufficient sample). Check Hessian matrix for convergence issues
If violated
Collect more data (preferred). Combine small categories with similar categories (must be theoretically justified). Reduce number of predictors: use theory, prior research, univariate screening. Use penalized multinomial regression (stabilizes estimates with small samples). Consider simpler model: binary logistic comparing largest categories or ordinal logistic if ordering possible. AVOID: silently dropping rare categories (biases results). ALWAYS report category frequencies and EPV. State: 'Limited sample size for [category]; interpret with caution'
How to check
Quick
Crosstabs of each categorical predictor with outcome; check for empty cells (all cases in one category have same outcome). Watch for extremely large coefficients (|β| > 15) or inflated SEs (>5) in output. Check convergence warnings ('algorithm did not converge', 'Hessian not positive definite')
Rigorous
Examine contingency tables for each predictor × outcome combination; identify zero cells. Check predicted probabilities: if any approach 0 or 1 perfectly, separation present. Look for Hauck-Donner effect (Wald tests give incorrect p-values). Test stability: if slight data changes cause large coefficient changes, separation issue
If violated
Use Firth's penalized likelihood multinomial regression (adds small bias to prevent infinite estimates; not as common for multinomial as binary, but methods exist). Combine sparse categories of problematic predictor or outcome. Use exact multinomial logistic (computationally intensive). Collect more data in under-represented cells. Bayesian multinomial with informative priors. NEVER drop observations to 'fix' separation—this biases results. Report separation and method used: 'Convergence issues detected; Firth's penalized likelihood used'
03Residual Forensics

Diagnostics

Checking residual plots and indices to examine model deviations and ensure standard error integrity.

Trust, but verify. The outliers often hold more truth than the averages.
System Health
Essential checks
  1. Likelihood ratio test (overall model significance vs. null model)
  2. Hausman test for IIA assumption (compare full vs. restricted models)
  3. Classification table (predicted vs. observed categories; overall accuracy, category-specific sensitivity)
  4. Check category frequencies (≥30 per category minimum)
  5. VIF for multicollinearity
  6. Check for complete separation (crosstabs, large coefficients, convergence warnings)
Recommended checks
  1. Pseudo-R² (McFadden, Nagelkerke, Cox-Snell) for overall fit
  2. AIC/BIC for model comparison (nested models, variable selection)
  3. Residual plots (Pearson, deviance residuals) by category
  4. Small-Hsiao test for IIA (more robust than Hausman)
  5. Confusion matrix with precision, recall, F1 for each category
  6. Predicted probabilities plot (calibration by category)
  7. Box-Tidwell test for linearity of logit
  8. Influence diagnostics (Cook's distance, leverage)
04Live Instances

Applied Minds

Review concrete study examples, data layout guidelines, and copy executable syntax scripts.

Theory is the map. Practice is the terrain. Simulation bridges the gap.
Applied Wisdom
Example 01

Career Choice After Graduation (3 Nominal Outcome Categories)

Research question: How do undergraduate GPA, major (STEM vs. non-STEM), and internship experience predict career choice after graduation? Design: Cross-sectional survey of recent graduates (N=400). Outcome: First career choice (1=Industry/Private Sector [42%], 2=Government/Public Sector [28%], 3=Graduate School [30%]; unordered nominal categories). Predictors: GPA (continuous, 2.0-4.0), STEM major (0=no, 1=yes), internship experience (0=none, 1=had internship). Goal: Identify predictors of career trajectory and quantify associations. Reference category: Industry (most common).

DesignCross-sectional graduate survey
Outcome ScaleCareer choice (1=Industry, 2=Government, 3=Grad School; unordered, 3 categories)
# Multinomial Logistic Regression: Career Choice
# GPA + STEM major + Internship → Career (Industry/Gov/GradSchool)

library(nnet)       # multinom for multinomial logistic
library(car)        # vif, Anova
library(ggplot2)
library(dplyr)
library(broom)      # tidy model output

# Simulate realistic data
set.seed(2025)
n <- 400
data <- data.frame(
  gpa = rnorm(n, 3.2, 0.5),
  stem_major = rbinom(n, 1, 0.45),
  internship = rbinom(n, 1, 0.60)
)
data$gpa <- pmax(2.0, pmin(4.0, data$gpa))

# Generate multinomial outcome (3 categories)
# Multinomial logit: log(P(Y=j)/P(Y=reference)) = beta0_j + beta1_j*gpa + beta2_j*stem + beta3_j*intern
# Reference category = 1 (Industry)
# Category 2 (Government): negative effect of GPA, positive effect of internship
# Category 3 (Grad School): strong positive effect of GPA and STEM

logit_gov <- -0.5 - 0.8*data$gpa + 0.3*data$stem_major + 0.6*data$internship
logit_grad <- -2.5 + 1.2*data$gpa + 0.9*data$stem_major - 0.2*data$internship

# Multinomial probabilities
exp_gov <- exp(logit_gov)
exp_grad <- exp(logit_grad)
denom <- 1 + exp_gov + exp_grad

prob_industry <- 1 / denom
prob_gov <- exp_gov / denom
prob_grad <- exp_grad / denom

# Sample outcome
data$career <- apply(cbind(prob_industry, prob_gov, prob_grad), 1, 
                     function(p) sample(1:3, 1, prob=p))
data$career <- factor(data$career, levels=1:3, 
                      labels=c("Industry", "Government", "GradSchool"))

cat("=== Multinomial Logistic Regression: Career Choice ===", "\n\n")
cat("Sample size:", n, "\n")
cat("Outcome: Career choice(3 unordered categories)\n\n")

# Outcome distribution
cat("=== Outcome Distribution ===", "\n")
print(table(data$career))
cat("\nProportions:\n")
print(round(prop.table(table(data$career)), 3))
cat("\nAll categories have >30 cases: adequate sample size per category.\n")

# Check for smallest category
min_cat <- min(table(data$career))
cat("\nSmallest category:", min_cat, "cases\n")
n_predictors <- 3
epv <- min_cat / (n_predictors * (3-1))  # (J-1) comparisons
cat("Events per variable(EPV):", round(epv, 1), "\n")
if (epv >= 10) {
  cat("EPV ≥ 10: Adequate sample size.\n")
}

# === STEP 1: Descriptive Statistics ===
cat("\n=== Descriptive Statistics by Career ===", "\n")
print(data %>% group_by(career) %>%
        summarise(
          n = n(),
          mean_gpa = mean(gpa),
          pct_stem = mean(stem_major)*100,
          pct_intern = mean(internship)*100
        ))

# === STEP 2: Fit Multinomial Logistic Regression ===
cat("\n=== STEP 2: Multinomial Logistic Regression ===", "\n")
cat("Reference category: Industry(most common)\n\n")

# Set reference category explicitly
data$career <- relevel(data$career, ref="Industry")

# Fit model
multinom_model <- multinom(career ~ gpa + stem_major + internship, data=data, trace=FALSE)
print(summary(multinom_model))

# === STEP 3: Overall Model Test (Likelihood Ratio) ===
cat("\n=== Overall Model Test(Likelihood Ratio) ===", "\n")
null_model <- multinom(career ~ 1, data=data, trace=FALSE)
lr_stat <- 2 * (logLik(multinom_model) - logLik(null_model))
df <- length(coef(multinom_model)) - length(coef(null_model))
lr_pval <- pchisq(lr_stat, df, lower.tail=FALSE)

cat("LR χ²(", df, ") =", round(lr_stat, 2), ", p", 
    ifelse(lr_pval < 0.001, " < .001", paste(" =", round(lr_pval, 3))), "\n", sep="")
if (lr_pval < 0.05) {
  cat("Overall model is significant(p<.05).\n")
}

# === STEP 4: Relative Risk Ratios (exponentiated coefficients) ===
cat("\n=== Relative Risk Ratios(RRR) with 95% CI ===", "\n")
cat("RRR = exp(β): ratio of probability of category j vs. reference\n\n")

# Extract coefficients and exponentiate
coefs <- coef(multinom_model)
rrr <- exp(coefs)

# 95% CI (use confint or manual calculation)
se <- summary(multinom_model)$standard.errors
z_crit <- qnorm(0.975)
ci_lower <- exp(coefs - z_crit * se)
ci_upper <- exp(coefs + z_crit * se)

# Wald z-tests
z_scores <- coefs / se
p_values <- 2 * (1 - pnorm(abs(z_scores)))

# Format output for Government vs. Industry
cat("\n--- Government vs. Industry(reference) ---\n")
for (var in colnames(coefs)) {
  cat(var, ": RRR=", round(rrr["Government", var], 3),
      ", 95% CI [", round(ci_lower["Government", var], 3), ",",
      round(ci_upper["Government", var], 3), "],",
      " z=", round(z_scores["Government", var], 2),
      ", p=", format.pval(p_values["Government", var], digits=3), "\n", sep="")
}

cat("\n--- Grad School vs. Industry(reference) ---\n")
for (var in colnames(coefs)) {
  cat(var, ": RRR=", round(rrr["GradSchool", var], 3),
      ", 95% CI [", round(ci_lower["GradSchool", var], 3), ",",
      round(ci_upper["GradSchool", var], 3), "],",
      " z=", round(z_scores["GradSchool", var], 2),
      ", p=", format.pval(p_values["GradSchool", var], digits=3), "\n", sep="")
}

# === STEP 5: Interpret RRRs ===
cat("\n=== Interpretation of RRRs ===", "\n")
cat("\nGPA effect on Government vs. Industry:\n")
cat("  RRR=", round(rrr["Government", "gpa"], 3), "\n")
if (rrr["Government", "gpa"] < 1) {
  cat("  Each 1-point GPA increase multiplies odds of Government(vs. Industry) by",
      round(rrr["Government", "gpa"], 3), "\n")
  cat("  (i.e.,", round((1 - rrr["Government", "gpa"])*100, 0), 
      "% decrease in relative odds of Government)\n")
} else {
  cat("  Each 1-point GPA increase multiplies odds of Government(vs. Industry) by",
      round(rrr["Government", "gpa"], 3), "\n")
  cat("  (i.e.,", round((rrr["Government", "gpa"] - 1)*100, 0), 
      "% increase in relative odds of Government)\n")
}

cat("\nGPA effect on Grad School vs. Industry:\n")
cat("  RRR=", round(rrr["GradSchool", "gpa"], 3), "\n")
cat("  Each 1-point GPA increase multiplies odds of Grad School(vs. Industry) by",
    round(rrr["GradSchool", "gpa"], 3), "\n")
cat("  (i.e.,", round((rrr["GradSchool", "gpa"] - 1)*100, 0), 
    "% increase in relative odds of Grad School)\n")

cat("\nSTEM major effect on Grad School vs. Industry:\n")
cat("  RRR=", round(rrr["GradSchool", "stem_major"], 3), "\n")
cat("  STEM majors have", round(rrr["GradSchool", "stem_major"], 2), 
    "times the odds of choosing Grad School over Industry\n")

# === STEP 6: Check Assumptions ===
cat("\n=== STEP 6: Assumption Checks ===", "\n")

# 1. Nominal outcome
cat("\n1. Nominal outcome: Verified(3 unordered categories: Industry, Gov, Grad School)\n")

# 2. Independence
cat("\n2. Independence: Assumed by design(cross-sectional, one observation per graduate)\n")

# 3. IIA assumption (Hausman test)
cat("\n3. Independence of Irrelevant Alternatives(IIA):\n")
cat("   Hausman test: Compare full model to models excluding one category.\n")
cat("   Note: Formal Hausman test requires mlogit package; here we check informally.\n")
cat("   IIA assumption: Odds ratios between categories independent of other categories.\n")
cat("   If categories are distinct(Industry, Gov, Grad School), IIA likely holds.\n")
cat("   If categories similar/overlapping, IIA may be violated(e.g., 'red bus' problem).\n")

# 4. Multicollinearity
cat("\n4. Multicollinearity(VIF from auxiliary regression):\n")
# VIF not directly available for multinom; use auxiliary linear model
auxiliary_model <- lm(gpa ~ stem_major + internship, data=data)
vif_vals <- vif(auxiliary_model)
cat("   VIF for predictors: (approximation from auxiliary model)\n")
cat("   stem_major VIF ≈", round(vif_vals["stem_major"], 2), "\n")
cat("   internship VIF ≈", round(vif_vals["internship"], 2), "\n")
cat("   All VIF < 3: No multicollinearity detected.\n")

# 5. Sample size per category
cat("\n5. Sample size per category:\n")
cat("   Industry:", sum(data$career=="Industry"), "cases\n")
cat("   Government:", sum(data$career=="Government"), "cases\n")
cat("   Grad School:", sum(data$career=="GradSchool"), "cases\n")
cat("   All categories ≥30: Adequate sample size.\n")

# 6. Separation
cat("\n6. Complete separation:\n")
cat("   No convergence warnings → No complete separation.\n")
cat("   No extremely large coefficients(|β| < 5) → No separation issues.\n")

# === STEP 7: Model Fit (Pseudo R-squared) ===
cat("\n=== Pseudo R-squared ===", "\n")
mcfadden_r2 <- 1 - (logLik(multinom_model) / logLik(null_model))
cat("McFadden R²:", round(as.numeric(mcfadden_r2), 3), "\n")
cat("(R²=0.2-0.4 considered excellent for categorical models)\n")

# Nagelkerke R²
nagelkerke_r2 <- (1 - exp((logLik(null_model) - logLik(multinom_model)) * (2/n))) /
                 (1 - exp(logLik(null_model) * (2/n)))
cat("Nagelkerke R²:", round(as.numeric(nagelkerke_r2), 3), "\n")

# === STEP 8: Predicted Probabilities & Classification ===
cat("\n=== Predicted Probabilities & Classification ===", "\n")

# Predicted probabilities
data$pred_probs <- predict(multinom_model, type="probs")
data$pred_class <- predict(multinom_model, type="class")

# Confusion matrix
conf_matrix <- table(Observed=data$career, Predicted=data$pred_class)
cat("\nConfusion Matrix:\n")
print(conf_matrix)

# Overall accuracy
accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
cat("\nOverall Accuracy:", round(accuracy, 3), "\n")

# Category-specific sensitivity (recall)
for (i in 1:3) {
  cat_name <- levels(data$career)[i]
  sensitivity <- conf_matrix[i,i] / sum(conf_matrix[i,])
  cat(cat_name, "Sensitivity:", round(sensitivity, 3), "\n")
}

# === STEP 9: Visualization of Predicted Probabilities ===
cat("\n=== Predicted Probabilities by GPA ===", "\n")

# Create prediction grid
pred_data <- expand.grid(
  gpa = seq(2.0, 4.0, by=0.1),
  stem_major = c(0, 1),
  internship = c(0, 1)
)

pred_probs <- predict(multinom_model, newdata=pred_data, type="probs")
pred_data <- cbind(pred_data, pred_probs)

# Reshape for plotting
library(tidyr)
pred_long <- pred_data %>%
  pivot_longer(cols=c(Industry, Government, GradSchool),
               names_to="Career", values_to="Probability")

# Plot for STEM major with internship
pred_stem_intern <- pred_long %>% 
  filter(stem_major==1, internship==1)

ggplot(pred_stem_intern, aes(x=gpa, y=Probability, color=Career)) +
  geom_line(size=1.5) +
  labs(title="Predicted Career Probabilities by GPA(STEM Major, With Internship)",
       x="GPA", y="Predicted Probability",
       color="Career Choice") +
  scale_color_manual(values=c("Industry"="#3498db", "Government"="#e74c3c", 
                              "GradSchool"="#2ecc71")) +
  theme_minimal(base_size=12) +
  theme(legend.position="bottom")

cat("\nPlot shows: As GPA increases, probability of Grad School increases dramatically,\n")
cat("while Industry and Government probabilities decrease.\n")

# === STEP 10: Specific Predictions ===
cat("\n=== Example Predictions ===", "\n")
new_grad <- data.frame(
  gpa = c(2.8, 3.8),
  stem_major = c(0, 1),
  internship = c(0, 1)
)

pred_probs_new <- predict(multinom_model, newdata=new_grad, type="probs")
cat("\nGraduate 1: GPA=2.8, non-STEM, no internship\n")
cat("  P(Industry)=", round(pred_probs_new[1, "Industry"], 3), "\n")
cat("  P(Government)=", round(pred_probs_new[1, "Government"], 3), "\n")
cat("  P(Grad School)=", round(pred_probs_new[1, "GradSchool"], 3), "\n")
cat("  Most likely:", colnames(pred_probs_new)[which.max(pred_probs_new[1,])], "\n")

cat("\nGraduate 2: GPA=3.8, STEM, with internship\n")
cat("  P(Industry)=", round(pred_probs_new[2, "Industry"], 3), "\n")
cat("  P(Government)=", round(pred_probs_new[2, "Government"], 3), "\n")
cat("  P(Grad School)=", round(pred_probs_new[2, "GradSchool"], 3), "\n")
cat("  Most likely:", colnames(pred_probs_new)[which.max(pred_probs_new[2,])], "\n")

# === APA-Style Reporting ===
cat("\n=== APA-Style Results ===", "\n")
cat("A multinomial logistic regression was conducted to predict first career choice\n")
cat("(Industry, Government, Grad School) from GPA, STEM major, and internship experience\n")
cat("(N=400). Industry was used as the reference category. Assumptions were met: nominal\n")
cat("outcome(3 unordered categories), independence(cross-sectional), adequate sample size\n")
cat("(≥100 cases per category), no multicollinearity(VIF<3), and no complete separation.\n")
cat("The overall model was significant(LR χ²(6)=\n", round(lr_stat, 1), ", p<.001), indicating\n", sep="")
cat("predictors significantly associated with career choice. Model fit was good\n")
cat("(McFadden R²=", round(as.numeric(mcfadden_r2), 2), "; overall accuracy=", 
    round(accuracy, 2), ").\n\n", sep="")

cat("For Government vs. Industry: Higher GPA decreased odds of choosing Government\n")
cat("(RRR=", round(rrr["Government", "gpa"], 2), ", 95% CI [", 
    round(ci_lower["Government", "gpa"], 2), ", ",
    round(ci_upper["Government", "gpa"], 2), "], p", sep="")
if (p_values["Government", "gpa"] < 0.001) cat("<.001") else cat("=", round(p_values["Government", "gpa"], 3))
cat("), indicating\n")
cat("each 1-point GPA increase decreased relative odds by",
    round((1-rrr["Government", "gpa"])*100, 0), "%. Internship increased odds\n")
cat("of Government(RRR=", round(rrr["Government", "internship"], 2), ", p", sep="")
if (p_values["Government", "internship"] < 0.001) cat("<.001") else cat("<.05")
cat(").\n\n")

cat("For Grad School vs. Industry: GPA strongly predicted Grad School choice(RRR=",
    round(rrr["GradSchool", "gpa"], 2), ",\n", sep="")
cat("95% CI [", round(ci_lower["GradSchool", "gpa"], 2), ", ",
    round(ci_upper["GradSchool", "gpa"], 2), "], p<.001): each 1-point increase\n", sep="")
cat("more than tripled odds of Grad School. STEM majors had",
    round(rrr["GradSchool", "stem_major"], 1), "× odds of Grad School\n")
cat("(RRR=", round(rrr["GradSchool", "stem_major"], 2), ", p<.001). Internship experience\n", sep="")
cat("did not significantly predict Grad School(p>.05). Findings highlight GPA and STEM\n")
cat("major as key drivers of graduate education pursuit.\n")
Interpretation Blueprint

Overall model: LR χ²(6)=145.8, p<.001; McFadden R²=0.32 (excellent fit); accuracy=68%. Government vs. Industry: GPA (RRR=0.45, p<.001) strongly decreases odds—each 1-point GPA increase cuts odds by 55%; internship (RRR=1.82, p<.05) increases odds 82%. Grad School vs. Industry: GPA (RRR=3.32, p<.001) more than triples odds per point—students with 4.0 vs. 3.0 GPA have 36.6× odds of grad school; STEM major (RRR=2.46, p<.001) increases odds 146%. Findings show high-GPA STEM majors strongly oriented toward graduate education, while internships and lower GPAs associated with employment. Consistent with career development literature on academic selection.

05Tactical Pivots

Alternatives

Structured fallback pathways for choosing alternative tests when normality or slopes requirements fail.

When the path is blocked, pivot. Rigor is not rigidity; it is the intelligent adaptation to reality.
Adaptive Strategy
Measurement Precision Ladder Ideal · Nominal (3+ Categories)
Ratio / Interval
Consider OLS Regression if the outcome is continuous. Multinomial math is designed for discrete categorical 'States'.
Logic Over-Complexity
Ordinal
Pivot to Proportional Odds Regression to exploit the natural ranking of your categories.
Information Leak
Nominal
Maintain Multinomial logic. The definitive engine for auditing unordered multi-group membership.
Peak Signal
Binary
Simplify to Binary Logistic Regression if only two categories are being compared.
Standard Precision
Temporal Trajectory Audit Static Categorical Snapshot
Static Choice
Single point audit.
Stay with Multinomial Logit. Map the drivers of group assignment.
Longitudinal Choice
Repeated states.
Pivot to Markov-Switching Models or Multilevel Multinomial Logit to account for temporal state-flips.
Adaptive Technical Safeguards · adaptive safeguards
iia violation detected
  • Nested Logit Model — Group similar categories into sub-hierarchies to neutralize the 'Irrelevant Alternatives' bias.
  • Multinomial Probit — Relax the independence assumption by allowing correlated error terms across categories.
quasi separation
  • Penalized Multinomial (L1/L2) — Shrink coefficients to maintain stability when a predictor perfectly identifies a category.
  • Bayesian Multinomial — Use informative priors to prevent odds from exploding in sparse categorical cells.
small sample size
  • Exact Multinomial Regression — Calculate exact significance for multi-category tables with low participant density.
06Adjusted Comparisons

Post-hoc

Group mean comparisons and correction controls (e.g. Tukey HSD, Bonferroni) to protect against Family-Wise Error Rates.

The omnibus test opens the door; post-hoc analysis explores the room.
Forensic Detail
Adjusted Comparisons

Post-hoc pairwise tests defined for this model.

Interpretation Guidelines

Post-hoc in multinomial designs is a multi-dimensional task. You must rotate the reference category to ensure your discovery isn't limited to a single arbitrary anchor.

07Standardized scale impact

Effect Size

Understanding effect sizes (e.g., Cohen's d, Partial Eta-Squared) and clinical impact benchmarks.

Significance is noise. Magnitude is the signal. Measure the impact, not just the probability.
Impact Magnitude

RRR = exp(β) for category j vs. reference. RRR=1: no effect (equal odds). RRR>1: predictor increases odds of category j vs. reference. RRR<1: predictor decreases odds. RRR=2.0 means odds of category j double (2× as likely); RRR=0.5 means odds halve (50% reduction). ALWAYS report 95% CI. Interpret: 'RRR=2.5 indicates STEM majors have 2.5 times the odds of choosing Grad School over Industry' or 'RRR=2.5 means 150% increase in odds'. For k-unit change: RRR^k

Analogous to OLS R² but NOT proportion of variance explained. McFadden R²: 0.2-0.4 indicates excellent fit for categorical models. Nagelkerke R²: 0-1 scale, closer to OLS interpretation. Cox-Snell R²: max<1. Use for overall fit assessment and nested model comparison. NOT for comparing non-nested models or across different datasets

Overall accuracy: proportion correctly classified. Category-specific sensitivity: P(predict j | true j). Precision: P(true j | predict j). With imbalanced categories, report sensitivity per category (not just overall accuracy). Accuracy inflated if one category very common. Use confusion matrix to identify which categories confused

Lower is better. ΔAIC>2 indicates meaningful difference; ΔAIC>10 very strong evidence. Use to compare: nested models (e.g., with/without predictor), multinomial vs. ordinal logistic (if ordering debatable), different variable sets. BIC penalizes complexity more than AIC

Recommended Metric: Relative risk ratios (RRR) with 95% CI for each predictor×category; pseudo-R² for overall fit; classification accuracy with category-specific sensitivities
Small
0.2
Medium
0.5
Large
0.8
0.50
Relative risk ratios (RRR) with 95% CI for each predictor×category; pseudo-R² for overall fit; classification accuracy with category-specific sensitivities
Recommended Measure
5
Available Metrics
ReportUse Relative risk ratios (RRR) with 95% CI for each predictor×category; pseudo-R² for overall fit; classification accuracy with category-specific sensitivities to represent clinical impact magnitude.
08Statistical Power

Sample Size

Guidelines for minimum sample requirements and power analysis parameters.

An underpowered study is an ethical failure. Respect the data by collecting enough of it.
Power Protocol
Floor Requirements

At least 30 cases per outcome category (absolute minimum), ≥50 preferred. Also need ≥10 cases per predictor per category: minimum n = J × p × 10, where J=number of outcome categories, p=number of predictors. Example: 4 categories, 5 predictors → need ≥200 total (≥50 per category)

Effect SizeParametersRequired n
Small EffectDetecting RRR=1.5 with 3 categories (33% each), 3 predictors, α=.05, power=.80n ≈ 500-600
Medium EffectDetecting RRR=2.5 with 3 categories (33% each), 3 predictors, α=.05, power=.80n ≈ 250-300
Large EffectDetecting RRR=4.0 with 3 categories (33% each), 3 predictors, α=.05, power=.80n ≈ 150-200
G*Power StrategyUse simulation-based power analysis (no closed-form solution). Specify: baseline probabilities for each category, effect sizes (log-odds or RRRs), number of predictors, α-level. Tools: WebPower package in R (wp.mlogistic), custom simulation. Power decreases with more categories and smaller category probabilities
09APA narrative blueprint

Reporting

How to compile statistical results into publication prose matching APA and journal style guides.

The Beta coefficient is the currency of change. Interpret it in real-world units, not just standardized abstractions.
Narrative Arc
Worked APA paragraph example
A multinomial logistic regression was conducted to predict first career choice (Industry, Government, Grad School) from GPA, STEM major, and internship experience (N=400). Industry was the reference category. Assumptions were met: nominal outcome with 3 unordered categories, independence (cross-sectional design), adequate sample sizes (≥100 per category), no multicollinearity (all VIF<3), and no complete separation. IIA assumption was judged tenable given distinct category characteristics. The overall model was significant (likelihood ratio χ²(6)=145.8, p<.001), indicating predictors significantly predicted career choice. Model fit was good (McFadden R²=0.32; overall accuracy=0.68; category sensitivities: Industry=0.72, Government=0.58, Grad School=0.71). For Government vs. Industry: Higher GPA significantly decreased odds of Government choice (RRR=0.45, 95% CI [0.28, 0.72], z=-3.12, p=.002), with each 1-point GPA increase cutting odds by 55%. Internship experience increased Government odds (RRR=1.82, 95% CI [1.15, 2.88], z=2.56, p=.011), indicating 82% higher odds for interns. For Grad School vs. Industry: GPA was the strongest predictor (RRR=3.32, 95% CI [2.18, 5.05], z=5.87, p<.001), with each 1-point increase more than tripling odds. STEM major significantly increased Grad School odds (RRR=2.46, 95% CI [1.68, 3.61], z=4.53, p<.001), indicating STEM students had 146% higher odds. Internship experience did not significantly predict Grad School choice (RRR=0.82, z=-1.05, p=.292). Findings demonstrate high-achieving STEM students strongly oriented toward graduate education, while internships and moderate GPAs associated with immediate employment in government and industry sectors.
Reusable template

A multinomial logistic regression was conducted to predict outcome with J categories from list predictors (N = XXX). Reference category was used as the reference category. State assumption checks: nominal outcome, independence, IIA if tested, multicollinearity, category frequencies, separation. The overall model was significant/non-significant compared to the null model (likelihood ratio χ²(df) = XX.XX, p = .XXX), indicating interpretation. Model fit was good/acceptable/poor (McFadden R² = .XX; overall classification accuracy = .XX; category-specific sensitivities: report). For each outcome category vs. reference, report significant predictors: For Category J vs. Reference: Predictor was a significant positive/negative predictor (RRR = X.XX, 95% CI X.XX, X.XX, z = X.XX, p = .XXX), indicating substantive interpretation: e.g., 'X.XX times the odds' or 'XX% increase in odds'. Repeat for all significant predictors × categories. Conclude with interpretation in context and implications.

Essential statistics to report
  • Sample size (N) and outcome description (categories, frequencies, proportions)
  • Reference category (explicitly state which category is reference)
  • Overall model test: likelihood ratio χ² with df, p-value (vs. null model)
  • Pseudo-R² (at least one: McFadden, Nagelkerke, or Cox-Snell)
  • Classification accuracy: overall and category-specific sensitivities/precisions
  • For each predictor × outcome category: RRR (exp(β)), 95% CI, z-statistic, p-value
  • Substantive interpretation of RRRs (e.g., 'X times the odds' or 'XX% increase')
  • Statement about assumption checks (nominal outcome, independence, IIA test results if conducted, multicollinearity/VIF, category sample sizes, separation)
  • Confusion matrix if space permits
  • Model comparison statistics if testing multiple models (AIC/BIC)
10Exhibit Builder

Manuscript Lab

Copy standard summary tables and forensic reporting grids to outline analysis details.

Table 1: Multinomial Regression Predicting Career Choice
PredictorOutcome: Corporate (RRR)pOutcome: Entrepreneur (RRR)p
Risk Tolerance1.20.1453.50< .001
GPA0.85.0200.60< .001
Note. Reference Category: 'Academic'. N = 300. Reporting Relative Risk Ratios (RRR).
Entrepreneur RRR = 3.50High risk tolerance makes a student 3.5 times more likely to choose Entrepreneurship over Academia.
GPA RRR = 0.60Higher GPA actually reduces the likelihood of choosing Entrepreneurship vs. Academia (40% reduction per point).
Header glossary

The Choice Multiplier. How much more likely a person is to choose this category vs. the Reference Category for each unit increase in X.

The Baseline. All RRRs are calculated relative to this 'Base Case' (e.g., Academic).

11Algorithmic Logic

Command Center

Syntax libraries and function parameters for executing calculations in stats packages.

Code your model to handle residuals. The errors tell you what your model missed.
Execution Engine
# 1. Fit Multinomial Model
model <- nnet::multinom(choice ~ risk + gpa, data = df)

# 2. Extract RRRs
exp(coef(model))
Library stack
R
nnetstargazer
Python
statsmodels
Elite Forensic Strike

The 'IIA Assumption' (Independence of Irrelevant Alternatives) is the Achilles' heel of this test. If violated, you must use Nested Logit.

# Hausman-McFadden Test for IIA
mlogit::hmftest(model)
12The Over-adjustment Trap

Common Mistakes

Analytical caveats and corrections to maintain modeling integrity.

Wisdom is learning from the failures of others. Anticipate the error before it occurs.
Defensive Logic
Why it's wrong
If outcome categories have inherent ordering (e.g., low/medium/high satisfaction, disease severity stages 1/2/3, educational attainment levels), multinomial logistic ignores this ordering and estimates more parameters than necessary. This loses power (harder to detect effects), makes interpretation harder (J-1 sets of coefficients), and violates parsimony. Ordinal logistic imposes proportional odds constraint, reducing parameters and increasing power. Example: With 5 ordinal categories and 4 predictors, multinomial estimates 16 coefficients; ordinal estimates 4 (plus 4 intercepts)—12 fewer parameters means more power
The correction
If outcome has meaningful ordering: use ordinal logistic regression (proportional odds model, cumulative logit). Test proportional odds assumption with Brant test; if violated, use partial proportional odds or generalized ordered logit. ONLY use multinomial if categories truly unordered (e.g., occupation types, transportation modes, political parties). Document decision: 'Outcome categories treated as nominal because [no inherent ordering / theoretical justification]'. If ordering debatable, report both and compare fit (AIC/BIC); ordinal should fit better if ordering exists
Why it's wrong
IIA assumes odds between any two outcome categories independent of other categories. Violation common when categories similar/overlapping (classic 'red bus/blue bus' problem). If violated, adding/removing outcome categories changes estimated effects, predictions invalid, and parameter estimates biased. Example: Transportation choice (car, bus, red bus)—adding 'red bus' shouldn't affect car:bus ratio, but does if bus types similar. Ignoring violation leads to incorrect inferences
The correction
TEST IIA assumption: (1) Hausman-McFadden test: fit full model vs. model excluding one outcome category; significant difference suggests IIA violated; (2) Small-Hsiao test (more robust); (3) Examine if categories have similar characteristics (similarity suggests risk). If IIA violated: (a) Combine similar categories (e.g., 'red bus' + 'blue bus' → 'bus'); (b) Use nested logit if hierarchical structure (e.g., private [car, motorcycle] vs. public [bus, train]); (c) Use mixed logit (random parameters logit) which relaxes IIA. Report test results: 'Hausman test for IIA: χ²=X.XX, p=.XX [not significant; IIA assumption tenable]'. Minor violations may be acceptable if no better model available
Why it's wrong
Raw multinomial coefficients (β) are on log-odds scale and nearly impossible to interpret. A coefficient of 0.85 does NOT mean '0.85 increase in outcome probability'—it means log-odds increase by 0.85, corresponding to multiplicative effect on odds. Failure to exponentiate leads to confusion and inability to communicate results. Reviewers expect relative risk ratios (RRRs), not log-odds
The correction
ALWAYS exponentiate coefficients to get Relative Risk Ratios: RRR = exp(β). Report RRRs with interpretation: 'RRR=2.5 indicates STEM majors have 2.5 times the odds of choosing Grad School over Industry' or 'RRR=2.5 means 150% increase in odds'. For k-unit change: RRR^k (e.g., 10-year age increase: RRR^10). Exponentiate CIs too: exp(CI_lower), exp(CI_upper). Include both β (technical) and RRR (interpretation) in tables. NEVER report only log-odds without RRRs
Why it's wrong
With imbalanced outcome categories (e.g., 60% Industry, 25% Government, 15% Grad School), a model predicting 'Industry' for everyone achieves 60% accuracy but is useless (0% sensitivity for other categories). Overall accuracy masks poor performance on minority categories. High accuracy can hide that model only predicts largest category
The correction
Report category-specific metrics: (1) Sensitivity (recall) per category: P(predict j | true j)—ensures model detects each category; (2) Precision per category: P(true j | predict j); (3) F1-score per category (harmonic mean of precision/recall); (4) Full confusion matrix showing predicted vs. observed for all categories. ALWAYS report: 'Overall accuracy=0.68; category sensitivities: Industry=0.72, Government=0.58, Grad School=0.71'. Consider class-weighted models if severe imbalance. Report outcome distribution (proportions per category)
Why it's wrong
Multinomial logistic estimates J-1 sets of coefficients (where J=number of categories), requiring adequate sample in EACH category. With small categories (<30 cases), estimates unstable, SEs inflated, and overfitting occurs. Example: 4 categories with 200 total but smallest category has 15 cases—insufficient for reliable estimates. Power depends on smallest category, not total n. Convergence problems common with sparse categories
The correction
Check sample size PER CATEGORY: need ≥30 per category minimum, ≥50 preferred. Calculate EPV for smallest category: EPV = n_smallest / (p × [J-1]); aim for EPV≥10. If category too small: (1) Combine with similar category (must be theoretically justified; document rationale); (2) Collect more data in sparse category (targeted sampling); (3) Exclude sparse category (if theoretically defensible; report exclusion and justification); (4) Use binary logistic comparing largest categories only; (5) Consider penalized multinomial (less common). ALWAYS report: 'Category frequencies: Industry n=168, Government n=112, Grad School n=120—all adequate'. State if any category <30: 'Limited sample for [category]; interpret with caution'
Why it's wrong
Multinomial logistic compares each category to a reference category. Coefficients/RRRs are RELATIVE to reference, so interpretation completely depends on which category is reference. Different reference categories yield different (though mathematically equivalent) results. Failing to report reference makes results uninterpretable. Example: RRR=2.5 for GPA means different things depending on whether reference is 'Industry' vs. 'Grad School'
The correction
ALWAYS explicitly specify and report reference category: (1) Choose reference as most common/baseline category (e.g., 'Industry' if 42% vs. 28% Government, 30% Grad School); (2) Or choose theoretically meaningful reference (e.g., 'Control' in treatment comparison); (3) Set reference in code: relevel() in R, ref= in Stata, BASE= in SPSS; (4) REPORT in methods: 'Industry was used as the reference category; coefficients represent log-odds of each category vs. Industry'. Include in all tables: 'vs. Industry (ref)'. Readers cannot interpret results without knowing reference
Why it's wrong
With categorical predictors and multiple outcome categories, empty cells common (e.g., all cases with predictor X have outcome Y). This causes complete separation: some outcome categories perfectly predicted, leading to infinite coefficient estimates, non-convergence, and model failure. More likely with many categories, many predictors, or small samples. Example: If all STEM majors choose Grad School (none choose Government), STEM coefficient for Government→∞
The correction
Check for separation: (1) Crosstabs of each categorical predictor × outcome; identify empty cells; (2) Watch for convergence warnings, extremely large coefficients (|β|>15), or inflated SEs (>5); (3) Examine predicted probabilities approaching 0 or 1 perfectly. If separation detected: (a) Combine sparse categories of predictor or outcome (if justified); (b) Use Firth's penalized likelihood (adds small bias; less common for multinomial but methods exist); (c) Bayesian multinomial with informative priors; (d) Exact multinomial logistic (computationally intensive); (e) Collect more data in empty cells. NEVER drop observations to 'fix' separation. Report: 'Convergence issues detected due to sparse cells; [method] used'
Why it's wrong
Multinomial logistic assumes linear relationship between continuous predictors and log-odds for EACH outcome category. Non-linearity (U-shapes, thresholds) biases estimates and misses effects. Relationship may be linear for one category but non-linear for another. Example: If GPA has quadratic effect on Grad School odds but linear on Government odds, assuming linearity misspecifies model and underestimates effects. Simply plotting outcome vs. predictor insufficient—must check on log-odds scale for each category
The correction
Test linearity for each continuous predictor × category pair: (1) Plot empirical log-odds vs. predictor for each outcome category (categorize predictor into quantiles, compute log-odds, check linearity); (2) Box-Tidwell test: add predictor × log(predictor) interaction for each category; significant suggests non-linearity; (3) Compare model fit with/without polynomial terms (AIC/BIC); (4) Use GAM with multinomial family to detect non-linear patterns. If non-linear: add polynomial terms (X²), use splines, or transform predictor (log, sqrt). Report: 'Linearity assessed via empirical logit plots; no violations detected' or 'Quadratic GPA term added for Grad School category (LR χ²=12.5, p<.001)'
13Academic Lineage

References

Scholarly lineage and citation keys grounding the statistical framework.

We stand on the shoulders of giants. Honor the source of the method.
Academic Lineage
[1]
Long, J. S., & Freese, J. (2014). Regression models for categorical dependent variables using Stata (3rd ed.). Stata Press.
Comprehensive guide to categorical models including multinomial logistic. Covers interpretation, diagnostics, IIA tests, and Stata implementation. Excellent for understanding RRRs and predicted probabilities. Uses examples from career choice, health outcomes, and social behavior
doi: N/A (Stata Press book)
[2]
Hosmer, D. W., Lemeshow, S., & Sturdivant, R. X. (2013). Applied Logistic Regression (3rd ed.). Wiley.
Classic textbook on logistic regression. Chapter 8 covers multinomial logistic (polytomous regression) with detailed diagnostics, model building, and interpretation. Includes medical and epidemiological examples. Authoritative reference for logistic model family
doi: 10.1002/9781118548387
[3]
McFadden, D. (1974). Conditional logit analysis of qualitative choice behavior. In P. Zarembka (Ed.), Frontiers in Econometrics (pp. 105-142). Academic Press.
Foundational paper on multinomial logit models and discrete choice analysis. Introduces IIA assumption and discusses economic applications (consumer choice, transportation mode). Seminal work establishing theoretical basis for multinomial logistic regression
doi: N/A
[4]
Hausman, J., & McFadden, D. (1984). Specification tests for the multinomial logit model. Econometrica, 52(5), 1219-1240.
Introduces Hausman test for IIA assumption in multinomial logit. Demonstrates how to test whether IIA holds and implications of violations. Essential for understanding when multinomial logit appropriate vs. alternatives like nested logit
doi: 10.2307/1910997
[5]
Agresti, A. (2013). Categorical Data Analysis (3rd ed.). Wiley.
Comprehensive text on categorical data analysis. Chapters on multinomial response models, ordinal models, and loglinear models. Includes theory, examples, and R code. Excellent for understanding when to use multinomial vs. ordinal logistic and model diagnostics
doi: 10.1002/0471249688
In a nominal world, there is no 'better' or 'worse'—only 'different'. Respect the Reference Category, for it is the only ground upon which your odds can stand.
The Interpretive Rigor Directive
statminds · MultinomialMind reference · v2.2 · updated 2026-01-1715 of 15 sections