6  Chương 5 - Mô hình Hồi quy cho Biến Phản hồi Đa phạm trù

Lời mở đầu

Thực tế kinh tế - xã hội hiếm khi chỉ đơn giản là những lựa chọn “có/không”. Con người và các tổ chức thường phải đối mặt với các quyết định có nhiều phương án: người tiêu dùng lựa chọn giữa nhiều thương hiệu khác nhau; nhà đầu tư phân bổ tài sản vào các kênh như cổ phiếu, trái phiếu hay bất động sản; các cơ quan xếp hạng tín dụng đánh giá doanh nghiệp theo các mức độ từ AAA đến D; hay khách hàng bày tỏ mức độ hài lòng của mình trên một thang đo nhiều bậc. Việc phân tích và dự báo các lựa chọn đa phạm trù này đòi hỏi những công cụ mô hình hóa chuyên biệt, vượt ra ngoài khuôn khổ của các mô hình hồi quy nhị phân.

Chương này sẽ mở rộng tầm nhìn của chúng ta từ thế giới nhị phân sang lĩnh vực đa dạng của các biến phản hồi có nhiều hơn hai kết quả. Một trong những sự phân biệt quan trọng nhất mà chúng ta sẽ học là giữa các biến phản hồi danh nghĩa (các lựa chọn không có thứ tự, ví dụ: chọn hãng xe nào) và các biến phản hồi thứ bậc (các lựa chọn có một trật tự tự nhiên, ví dụ: mức độ hài lòng). Sự phân biệt này là cốt lõi, vì nó dẫn đến hai họ mô hình hoàn toàn khác nhau.

Chúng ta sẽ khám phá chi tiết Mô hình Multinomial Logit, công cụ tiêu chuẩn để phân tích các lựa chọn danh nghĩa, và tìm hiểu về giả định quan trọng “IIA” của nó. Tiếp theo, chúng ta sẽ đi sâu vào Mô hình Ordinal Logit (hay Proportional Odds Model), một phương pháp thanh lịch và hiệu quả để mô hình hóa các kết quả có thứ tự, cùng với giả định “Tỷ lệ chênh tương ứng” (Proportional Odds) của nó.

Xuyên suốt chương, chúng ta sẽ không chỉ học lý thuyết mà còn tập trung vào kỹ năng thực hành: làm thế nào để xây dựng các mô hình này trong R, diễn giải các hệ số của chúng một cách chính xác, kiểm tra các giả định quan trọng, và đánh giá hiệu năng của mô hình. Thành thạo các kỹ thuật trong chương này sẽ cho phép bạn phân tích một lớp các vấn đề phức tạp và phổ biến trong nghiên cứu thị trường, quản trị, tài chính và chính sách công.

Mục tiêu chương

Sau khi hoàn thành chương này, người học sẽ có khả năng:

6.1 Phân biệt biến phản hồi đa phạm trù danh nghĩa và thứ bậc

Khi một biến định tính có nhiều hơn hai phạm trù, chúng ta gọi nó là biến đa phạm trù (polytomous hoặc multinomial). Bước đầu tiên và quan trọng nhất là xác định xem các phạm trù đó có một trật tự tự nhiên hay không.

  • Biến phản hồi danh nghĩa (Nominal Response Variable): Là một biến định tính có từ ba phạm trù trở lên mà giữa chúng không có một thứ tự hay xếp hạng vốn có. Các phạm trù chỉ đơn thuần là các nhãn khác nhau. Việc thay đổi thứ tự của các phạm trù không làm thay đổi bản chất của biến.
    • Ví dụ:
      • Lựa chọn nghề nghiệp: Kế toán, Marketing, Nhân sự, Tài chính.
      • Lựa chọn phương tiện đi làm: Xe buýt, Xe máy, Ô tô cá nhân, Tàu điện.
      • Nhà cung cấp dịch vụ Internet: Viettel, FPT, VNPT.
      • Ngành công nghiệp của một công ty: Sản xuất, Dịch vụ, Nông nghiệp.
  • Biến phản hồi thứ bậc (Ordinal Response Variable): Là một biến định tính có từ ba phạm trù trở lên mà giữa chúng có một trật tự hoặc xếp hạng tự nhiên. Chúng ta có thể nói phạm trù này “cao hơn”, “tốt hơn”, hoặc “nhiều hơn” phạm trù kia, nhưng khoảng cách giữa các bậc không nhất thiết phải bằng nhau.
    • Ví dụ:
      • Mức độ hài lòng của khách hàng: Rất không hài lòng < Không hài lòng < Bình thường < Hài lòng < Rất hài lòng.
      • Xếp hạng tín dụng: D < C < B < BB < BBB < A < AA < AAA.
      • Phản ứng với một chính sách: Phản đối mạnh mẽ < Phản đối < Trung lập < Ủng hộ < Ủng hộ mạnh mẽ.
      • Tình trạng sức khỏe tự đánh giá: Rất kém < Kém < Trung bình < Tốt < Rất tốt.

Bảng dưới đây tóm tắt sự khác biệt:

Đặc điểm Biến Phản hồi Danh nghĩa Biến Phản hồi Thứ bậc
Thứ tự phạm trù Không có ý nghĩa Có ý nghĩa và quan trọng
Thông tin chứa đựng Chỉ phân loại Phân loại và Xếp hạng
Ví dụ Lựa chọn thương hiệu, Ngành nghề Mức độ hài lòng, Xếp hạng tín dụng
Mô hình phù hợp Multinomial Logit Ordinal Logit (Proportional Odds)

Sử dụng sai mô hình—ví dụ, dùng mô hình danh nghĩa cho dữ liệu thứ bậc—sẽ làm lãng phí thông tin về thứ tự và có thể làm giảm hiệu quả của mô hình. Ngược lại, áp dụng mô hình thứ bậc cho dữ liệu danh nghĩa là không hợp lệ vì nó áp đặt một cấu trúc thứ tự không tồn tại.

6.2 Mô hình Multinomial Logit (cho biến phản hồi danh nghĩa)

Mô hình Multinomial Logit, còn được gọi là mô hình hồi quy đa thức, là sự mở rộng trực tiếp của mô hình logistic nhị phân cho các biến phản hồi có nhiều hơn hai phạm trù danh nghĩa.

6.2.1 Mô hình Baseline-Category Logit: Xây dựng và lựa chọn phạm trù cơ sở

Ý tưởng cốt lõi của mô hình này là chọn một trong các phạm trù (\(J\)) làm phạm trù cơ sở (baseline category) và sau đó xây dựng \(J-1\) mô hình logistic nhị phân riêng biệt, mỗi mô hình so sánh một phạm trù khác với phạm trù cơ sở.

Giả sử biến phản hồi \(Y\) có \(J\) phạm trù. Chúng ta chọn phạm trù thứ \(J\) làm cơ sở. Mô hình sẽ ước lượng \(J-1\) phương trình logit sau: \[ \log\left(\frac{P(Y=1)}{P(Y=J)}\right) = \beta_{10} + \beta_{11}X_1 + ... + \beta_{1k}X_k \\ \log\left(\frac{P(Y=2)}{P(Y=J)}\right) = \beta_{20} + \beta_{21}X_1 + ... + \beta_{2k}X_k \\ ... \\ \log\left(\frac{P(Y=J-1)}{P(Y=J)}\right) = \beta_{(J-1)0} + \beta_{(J-1)1}X_1 + ... + \beta_{(J-1)k}X_k \] Mỗi phương trình có một bộ hệ số \(\beta\) riêng.

Từ các phương trình này, xác suất để một quan sát rơi vào mỗi phạm trù \(j\) có thể được tính toán: \[ P(Y=j) = \frac{\exp(\eta_j)}{1 + \sum_{h=1}^{J-1} \exp(\eta_h)} \quad \text{cho } j=1, ..., J-1 \\ P(Y=J) = \frac{1}{1 + \sum_{h=1}^{J-1} \exp(\eta_h)} \] trong đó \(\eta_j\) là bộ tiên đoán tuyến tính cho phạm trù \(j\).

Lựa chọn phạm trù cơ sở: Việc lựa chọn phạm trù cơ sở không ảnh hưởng đến xác suất dự báo hay độ phù hợp chung của mô hình, nhưng nó ảnh hưởng đến việc diễn giải các hệ số. Thường người ta chọn: * Phạm trù có tần suất lớn nhất. * Một phạm trù “tham chiếu” hoặc “status quo” có ý nghĩa (ví dụ: “không lựa chọn gì”). * R mặc định thường chọn phạm trù đầu tiên theo thứ tự alphabet hoặc thứ tự của factor.

6.2.2 Diễn giải ý nghĩa tham số (log-odds và tỷ số chênh tương đối)

Các hệ số \(\beta_{jk}\) trong mô hình multinomial logit biểu thị sự thay đổi trong log-odds của việc lựa chọn phạm trù \(j\) so với phạm trù cơ sở \(J\), khi biến độc lập \(X_k\) tăng một đơn vị.

Cách diễn giải hữu ích hơn là sử dụng Tỷ số Chênh Tương đối (Relative Risk Ratio - RRR), đôi khi cũng được gọi là Odds Ratio. \(RRR = \exp(\beta_{jk}) = \frac{Odds(Y=j|X_k+1)}{Odds(Y=j|X_k)} / \frac{Odds(Y=J|X_k+1)}{Odds(Y=J|X_k)}\)

Diễn giải: Khi \(X_k\) tăng 1 đơn vị, odds của việc chọn phạm trù \(j\) so với odds của việc chọn phạm trù cơ sở \(J\) sẽ thay đổi một lượng bằng \(\exp(\beta_{jk})\) lần, khi các biến khác không đổi.

6.2.3 Ước lượng (MLE) và suy diễn thống kê

  • Mô hình được ước lượng bằng MLE, dựa trên phân phối Đa thức (Multinomial).
  • Suy diễn thống kê cho các hệ số riêng lẻ (ví dụ, \(\beta_{jk}\)) vẫn dựa trên kiểm định Wald (z-test hoặc t-test trong một số output).
  • Để kiểm tra ý nghĩa tổng thể của một biến độc lập \(X_k\) trên tất cả các lựa chọn, chúng ta cần thực hiện một kiểm định Tỷ số hợp lý (LRT) so sánh mô hình có \(X_k\) và mô hình không có \(X_k\). Kiểm định này sẽ có \(J-1\) bậc tự do.

6.2.4 Kiểm định giả định IIA (Independence of Irrelevant Alternatives)

Đây là giả định quan trọng và cũng là điểm yếu tiềm tàng của mô hình multinomial logit. * Nội dung giả định: Tỷ số xác suất của việc lựa chọn giữa hai phương án bất kỳ chỉ phụ thuộc vào các đặc điểm của hai phương án đó, và không bị ảnh hưởng bởi sự có mặt hay vắng mặt của các phương án (“không liên quan”) khác. \(\frac{P(Y=j)}{P(Y=h)}\) không phụ thuộc vào các phương án còn lại. * Vấn đề “Red bus/Blue bus”: Giả sử bạn chọn giữa đi làm bằng ô tô và xe buýt màu đỏ. Tỷ lệ odds là 2:1. Bây giờ, một lựa chọn mới xuất hiện: xe buýt màu xanh, giống hệt xe buýt đỏ. Theo giả định IIA, tỷ lệ odds giữa ô tô và xe buýt đỏ vẫn phải là 2:1. Nhưng thực tế, nhiều người sẽ chia sự lựa chọn “xe buýt” của mình ra, và tỷ lệ odds giữa ô tô và xe buýt đỏ có thể trở thành 2:0.5 (tức 4:1). Giả định IIA bị vi phạm. * Khi nào IIA có thể bị vi phạm: Khi các phương án lựa chọn có sự tương đồng hoặc có thể thay thế cho nhau. * Kiểm định: * Hausman-McFadden test: Ý tưởng là ước lượng mô hình đầy đủ, sau đó loại bỏ một phạm trù và ước lượng lại mô hình trên tập con các phạm trù còn lại. Nếu IIA đúng, các hệ số không nên thay đổi một cách có hệ thống. * Small-Hsiao test: Tương tự nhưng ổn định hơn với mẫu nhỏ. * Trong R, các kiểm định này có thể được thực hiện bằng gói mlogit. * Nếu IIA bị vi phạm: Cần xem xét các mô hình thay thế như Nested Logit (nhóm các lựa chọn tương tự nhau lại) hoặc Multinomial Probit (phức tạp hơn nhiều).

6.3 Mô hình cho Biến Phản hồi Thứ bậc (Ordinal Response Models)

Khi biến phản hồi có các phạm trù được sắp xếp theo thứ tự, chúng ta nên sử dụng thông tin về thứ tự này. Mô hình phổ biến nhất là mô hình Ordinal Logit.

6.3.1 Mô hình Ordinal Logit (Proportional Odds Logistic Regression - POLR)

Mô hình này không mô hình hóa xác suất của một phạm trù cụ thể, mà là xác suất tích lũy (cumulative probability).

6.3.1.1 5.3.1.1. Xây dựng mô hình dựa trên logit tích lũy (Cumulative logits)

Cho một biến phản hồi thứ bậc \(Y\) có \(J\) phạm trù được sắp xếp \(1 < 2 < ... < J\). Mô hình xem xét xác suất để \(Y\) nằm ở hoặc dưới một mức \(j\) nào đó, \(P(Y \le j)\). Mô hình định nghĩa \(J-1\) logit tích lũy: \[ \text{logit}[P(Y \le j)] = \log\left(\frac{P(Y \le j)}{P(Y > j)}\right) = \alpha_j - (\beta_1 X_1 + ... + \beta_k X_k) \] cho \(j = 1, 2, ..., J-1\). (Lưu ý: Một số phần mềm, bao gồm polr trong R, dùng công thức \(\alpha_j + \beta_1 X_1 + ...\). Dấu của \(\beta\) sẽ ngược lại nhưng ý nghĩa không đổi. Chúng ta sẽ theo công thức này cho tiện với R.) \[ \text{logit}[P(Y \le j)] = \alpha_j + \beta_1 X_1 + ... + \beta_k X_k \]

Các thành phần của mô hình:

  • \(\alpha_j\) (Thresholds/Cut-points): Đây là các hệ số chặn (intercepts), có \(J-1\) giá trị khác nhau. Chúng xác định ngưỡng để chuyển từ phạm trù này sang phạm trù khác trên thang đo biến tiềm ẩn. \(\alpha_1 < \alpha_2 < ... < \alpha_{J-1}\).
  • \(\beta_k\) (Coefficients): Đây là các hệ số hồi quy cho các biến độc lập. Điểm mấu chốt là chỉ có một bộ hệ số \(\beta\) duy nhất cho tất cả \(J-1\) phương trình logit.

6.3.1.2 5.3.1.2. Diễn giải ý nghĩa tham số (cumulative odds ratios)

Hệ số \(\beta_k\) biểu thị sự thay đổi trong log-odds của việc \(Y \le j\) (so với \(Y > j\)) khi \(X_k\) tăng một đơn vị. Khi lấy hàm mũ, \(\exp(\beta_k)\), chúng ta có Tỷ số chênh tích lũy (Cumulative Odds Ratio). * Diễn giải: Khi \(X_k\) tăng 1 đơn vị, odds của việc kết quả nằm ở hoặc dưới một mức bất kỳ \(j\) (so với trên mức \(j\)) sẽ thay đổi một lượng bằng \(\exp(\beta_k)\) lần, khi các biến khác không đổi. * Nếu \(\beta_k > 0\), \(\exp(\beta_k)>1\): Tăng \(X_k\) làm tăng odds rơi vào các phạm trù cao hơn của \(Y\). * Nếu \(\beta_k < 0\), \(\exp(\beta_k)<1\): Tăng \(X_k\) làm tăng odds rơi vào các phạm trù thấp hơn của \(Y\).

6.3.1.3 5.3.1.3. Ước lượng (MLE) và suy diễn thống kê

  • Mô hình được ước lượng bằng MLE.
  • Suy diễn thống kê cho các hệ số \(\beta\) và các ngưỡng \(\alpha\) dựa trên kiểm định Wald (t-test trong polr summary).

6.3.1.4 5.3.1.4. Kiểm định giả định Proportional Odds (Parallel Lines Assumption)

Đây là giả định cốt lõi của mô hình POLR. * Nội dung giả định: Tác động của mỗi biến độc lập (hệ số \(\beta\)) là như nhau (không đổi) qua tất cả các điểm cắt (cut-points) \(j = 1, ..., J-1\). Điều này có nghĩa là mối quan hệ giữa \(X\) và \(Y\) không phụ thuộc vào việc bạn đang so sánh \(Y \le 1\) với \(Y > 1\), hay \(Y \le 2\) với \(Y > 2\), v.v. * Kiểm định: * Brant test: Đây là kiểm định phổ biến nhất. Nó ước lượng các mô hình logistic nhị phân riêng biệt tại mỗi điểm cắt và so sánh các hệ số \(\beta\) tương ứng. * \(H_0\): Giả định Proportional Odds được thỏa mãn. * \(H_a\): Giả định Proportional Odds bị vi phạm. * Một p-value lớn (ví dụ > 0.05) là điều mong muốn, cho thấy không có bằng chứng để bác bỏ giả định này. * Trong R, gói brant cung cấp hàm brant().

6.3.2 Các mô hình khác cho biến phản hồi thứ bậc (giới thiệu)

Nếu giả định Proportional Odds bị vi phạm, bạn có thể cân nhắc các mô hình thay thế: * Partial Proportional Odds Model: Cho phép một vài biến độc lập vi phạm giả định (có hệ số \(\beta\) khác nhau ở các điểm cắt), trong khi các biến khác vẫn tuân thủ. * Adjacent-Categories Logit Model: Mô hình hóa logit của việc rơi vào phạm trù \(j\) so với phạm trù ngay liền kề nó (\(j+1\)). * Continuation-Ratio Logit Model: Mô hình hóa logit của việc rơi vào phạm trù \(j\) với điều kiện đã biết đối tượng nằm ở hoặc trên phạm trù \(j\).

6.4 Đánh giá và lựa chọn mô hình cho biến phản hồi đa phạm trù

  1. Các tiêu chí thông tin (AIC, BIC): Vẫn là công cụ hữu ích để so sánh các mô hình khác nhau (ví dụ, mô hình có các bộ biến độc lập khác nhau) trên cùng một bộ dữ liệu. Mô hình có AIC/BIC thấp hơn được ưa thích.
  2. Độ chính xác dự báo (Confusion matrix):
    • Sau khi xây dựng mô hình, ta có thể dự báo phạm trù có xác suất cao nhất cho mỗi quan sát.
    • So sánh phạm trù dự báo với phạm trù thực tế để tạo ra một ma trận nhầm lẫn (confusion matrix) \(J \times J\).
    • Từ đó, tính được độ chính xác tổng thể (overall accuracy) và độ chính xác cho từng phạm trù.

6.5 Ứng dụng R trong xây dựng và đánh giá mô hình hồi quy đa phạm trù

6.5.1 Mô hình Multinomial Logit (gói nnet và mlogit)

Hàm chính là multinom() trong gói nnet (đi kèm với R).

# Ví dụ: Dự đoán loài hoa Iris dựa trên các đặc điểm
data(iris)
str(iris) # Species là biến phản hồi danh nghĩa 3 mức

# Cài đặt và tải gói nnet
library(nnet)

# Xây dựng mô hình, chọn 'setosa' làm mức cơ sở
# relevel() để thay đổi mức cơ sở
iris$Species <- relevel(iris$Species, ref = "setosa")
iris_multinom_model <- multinom(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
                                data = iris)
summary(iris_multinom_model)

# Diễn giải hệ số:
# Output có 2 bộ hệ số: một cho 'versicolor' và một cho 'virginica' (so với 'setosa')
# Tính RRR (exp(coef))
exp(coef(iris_multinom_model))

# Đánh giá độ chính xác
pred_species <- predict(iris_multinom_model, newdata = iris, type = "class")
confusion_matrix_multi <- table(Predicted = pred_species, Actual = iris$Species)
print(confusion_matrix_multi)

accuracy_multi <- sum(diag(confusion_matrix_multi)) / sum(confusion_matrix_multi)
cat("Overall Accuracy:", round(accuracy_multi, 4), "\n")

6.5.2 Mô hình Ordinal Logit (gói MASS và ordinal)

Hàm polr() (Proportional Odds Logistic Regression) trong gói MASS là hàm cổ điển.

# Ví dụ: Đánh giá sự hài lòng về nhà ở dựa trên các yếu tố
# Dữ liệu `housing` trong gói MASS
library(MASS)
data(housing)
str(housing) # Sat là biến phản hồi thứ bậc: Low < Medium < High

# Xây dựng mô hình
housing_polr_model <- polr(Sat ~ Infl + Type + Cont, data = housing, Hess = TRUE)
# Hess=TRUE để tính ma trận Hessian, cần cho summary()
summary(housing_polr_model)

# Diễn giải hệ số:
# Coefs là log-odds, Intercepts là các ngưỡng
# Tính OR
exp(coef(housing_polr_model))

# Kiểm tra giả định Proportional Odds
# install.packages("brant")
library(brant)
brant_test_result <- brant(housing_polr_model)
print(brant_test_result)
# Xem p-value của dòng "Omnibus". Nếu lớn, giả định không bị vi phạm.

# Đánh giá độ chính xác
pred_satisfaction <- predict(housing_polr_model, newdata = housing, type = "class")
confusion_matrix_ord <- table(Predicted = pred_satisfaction, Actual = housing$Sat)
print(confusion_matrix_ord)

accuracy_ord <- sum(diag(confusion_matrix_ord)) / sum(confusion_matrix_ord)
cat("Overall Accuracy:", round(accuracy_ord, 4), "\n")

6.6 Tóm tắt chương

Chương 5 đã trang bị cho chúng ta các phương pháp mô hình hóa chuyên biệt cho các biến phản hồi có nhiều hơn hai phạm trù, một tình huống rất phổ biến trong thực tế.

  1. Phân loại biến phản hồi: Nền tảng của chương là sự phân biệt cốt lõi giữa biến phản hồi danh nghĩa (không có thứ tự) và thứ bậc (có thứ tự tự nhiên), dẫn đến việc lựa chọn các họ mô hình khác nhau.
  2. Mô hình Multinomial Logit: Chúng ta đã tìm hiểu sâu về mô hình Baseline-Category Logit để phân tích các lựa chọn danh nghĩa. Điểm mấu chốt là mô hình này tạo ra \(J-1\) bộ hệ số so sánh với một phạm trù cơ sở, và việc diễn giải chúng dựa trên Tỷ số chênh tương đối (RRR). Giả định quan trọng IIA và các vấn đề liên quan cũng đã được thảo luận.
  3. Mô hình Ordinal Logit (POLR): Đối với dữ liệu thứ bậc, chúng ta đã học về mô hình Proportional Odds, một phương pháp hiệu quả tận dụng thông tin về thứ tự. Mô hình này ước lượng một bộ hệ số \(\beta\) duy nhất cho các biến độc lập và các ngưỡng \(\alpha_j\) riêng biệt. Diễn giải dựa trên Tỷ số chênh tích lũy. Giả định nền tảng Proportional Odds và cách kiểm định nó bằng Brant test đã được trình bày chi tiết.
  4. Đánh giá và Thực hành: Chúng ta đã học cách đánh giá các mô hình này thông qua các tiêu chí thông tin (AIC/BIC) và độ chính xác dự báo (ma trận nhầm lẫn). Quan trọng hơn, chúng ta đã thực hành xây dựng, diễn giải và kiểm định các mô hình này trong R bằng các gói nnet (cho multinom), MASS (cho polr), và brant (để kiểm tra giả định).

Với kiến thức từ chương này, người học có thể tự tin tiếp cận các bài toán phân tích lựa chọn phức tạp, từ dự báo thương hiệu mà khách hàng sẽ chọn, đến mô hình hóa các thang điểm đánh giá trong khảo sát.

6.7 Case studies

Case Study 5.1: Dự báo Lựa chọn Phương tiện đi lại

  • Bối cảnh: Một cơ quan quy hoạch đô thị muốn hiểu các yếu tố ảnh hưởng đến việc người dân lựa chọn phương tiện đi lại chính (mode: Car, Bus, Train). Các yếu tố bao gồm income (thu nhập), time (thời gian di chuyển), và cost (chi phí di chuyển).

  • Nhiệm vụ:

    1. Xác định loại biến phản hồi và mô hình phù hợp.
    2. Xây dựng mô hình Multinomial Logit trong R, với “Car” là phạm trù cơ sở.
    3. Diễn giải RRR của biến time đối với lựa chọn “Bus” so với “Car”.
    4. Dự báo lựa chọn có xác suất cao nhất cho một người có income, time, cost cụ thể.

  • Phân tích bằng R:

    # Giả lập dữ liệu
set.seed(501)
n_mode <- 500
df_mode <- data.frame(
  income = runif(n_mode, 10, 100),
  time_car = rnorm(n_mode, 30, 10),
  cost_car = rnorm(n_mode, 5, 2),
  time_bus = rnorm(n_mode, 50, 15),
  cost_bus = rnorm(n_mode, 2, 0.5),
  time_train = rnorm(n_mode, 40, 12),
  cost_train = rnorm(n_mode, 4, 1)
)
# Mô phỏng lựa chọn dựa trên 'utility'
utility_car = -0.1*df_mode$time_car - 0.5*df_mode$cost_car + 0.02*df_mode$income
utility_bus = -0.1*df_mode$time_bus - 0.5*df_mode$cost_bus
utility_train = -0.1*df_mode$time_train - 0.5*df_mode$cost_train

max_utility <- pmax(utility_car, utility_bus, utility_train)
df_mode$mode <- ifelse(max_utility == utility_car, "Car",
                       ifelse(max_utility == utility_bus, "Bus", "Train"))
df_mode$mode <- factor(df_mode$mode, levels = c("Car", "Bus", "Train"))

# Cần định dạng lại dữ liệu cho multinom (dạng 'wide' sang 'long' là tốt nhất, 
# nhưng để đơn giản, ta dùng mô hình với các đặc điểm cá nhân)
# Giả sử ta chỉ có đặc điểm cá nhân
df_choice <- data.frame(
    mode = df_mode$mode,
    income = df_mode$income,
    household_size = sample(1:5, n_mode, replace=TRUE)
)

library(nnet)
mode_model <- multinom(mode ~ income + household_size, data = df_choice)
summary(mode_model)

# Diễn giải RRR
rrr_mode <- exp(coef(mode_model))
print(rrr_mode)
# Diễn giải: RRR của income cho lựa chọn Bus (so với Car) là...

# Dự báo
new_person <- data.frame(income = 60, household_size = 3)
predict(mode_model, newdata = new_person, type = "probs") # xác suất
predict(mode_model, newdata = new_person, type = "class") # lựa chọn

Case Study 5.2: Mô hình hóa Xếp hạng Trái phiếu Doanh nghiệp

  • Bối cảnh: Một quỹ đầu tư muốn xây dựng mô hình dự báo xếp hạng tín dụng của trái phiếu doanh nghiệp (rating: AAA, AA, A, BBB, BB, B). Các biến dự báo bao gồm profit_margin (tỷ suất lợi nhuận), leverage_ratio (tỷ lệ đòn bẩy), và company_size (quy mô công ty: Small, Medium, Large).

  • Nhiệm vụ:

    1. Xác định loại biến phản hồi và mô hình phù hợp.
    2. Xây dựng mô hình Ordinal Logit.
    3. Kiểm tra giả định Proportional Odds.
    4. Diễn giải ý nghĩa của hệ số leverage_ratio.

  • Phân tích bằng R:

    set.seed(502)
n_bond <- 300
rating_levels <- c("B", "BB", "BBB", "A", "AA", "AAA")
profit_margin <- runif(n_bond, 0.05, 0.25)
leverage_ratio <- runif(n_bond, 0.2, 0.8)
company_size <- factor(sample(c("Small", "Medium", "Large"), n_bond, replace = TRUE),
                       levels = c("Small", "Medium", "Large"), ordered = TRUE)

# Tạo biến tiềm ẩn
latent_rating <- 2*profit_margin - 5*leverage_ratio + 1*as.numeric(company_size) + rnorm(n_bond, 0, 1)
# Chia thành các hạng
cutpoints <- quantile(latent_rating, probs = c(0.15, 0.3, 0.55, 0.75, 0.9))
rating <- cut(latent_rating, breaks = c(-Inf, cutpoints, Inf), labels = rating_levels)
bond_data <- data.frame(rating, profit_margin, leverage_ratio, company_size)

library(MASS)
rating_model <- polr(rating ~ profit_margin + leverage_ratio + company_size, 
                     data = bond_data, Hess = TRUE)
summary(rating_model)

library(brant)
brant(rating_model) # p-value lớn cho thấy giả định được thỏa mãn

# Diễn giải leverage_ratio
# Hệ số âm -> khi leverage_ratio tăng, odds rơi vào hạng thấp hơn (B, BB) sẽ tăng
# exp(coef(rating_model)["leverage_ratio"]) -> OR

Case Study 5.3: Phân tích các yếu tố ảnh hưởng đến Mức độ Hài lòng công việc

  • Bối cảnh: Phòng nhân sự muốn biết các yếu tố nào (salary, work_life_balance_score 1-10, promotion_opportunity Yes/No) ảnh hưởng đến mức độ hài lòng công việc của nhân viên (job_satisfaction: Low, Medium, High).
  • Nhiệm vụ: Xây dựng mô hình Ordinal Logit và diễn giải kết quả. Đánh giá độ chính xác của mô hình.
  • Phân tích bằng R: Tương tự Case Study 5.2, sử dụng polr() và sau đó tạo ma trận nhầm lẫn để tính độ chính xác.

Case Study 5.4: Kiểm định Giả định IIA trong lựa chọn sản phẩm

  • Bối cảnh: Một công ty nghiên cứu thị trường có dữ liệu về lựa chọn smartphone của người dùng (brand: Apple, Samsung, Google, Other). Họ xây dựng mô hình multinomial logit và muốn kiểm tra giả định IIA.

  • Nhiệm vụ:

    1. Xây dựng mô hình multinomial logit.
    2. Sử dụng gói mlogit để thực hiện kiểm định Hausman-McFadden.

  • Phân tích bằng R:

    # Dữ liệu cần có dạng "long" cho mlogit
# Giả sử đã có dữ liệu dạng wide
# set.seed(504)
# ... tạo dữ liệu brand_choice_data ...

# install.packages("mlogit")
# library(mlogit)
# Dữ liệu cần được định dạng lại bằng mlogit.data()
# mlogit_data <- mlogit.data(brand_choice_data, choice = "brand", shape = "wide")

# Xây dựng mô hình
# mlogit_model <- mlogit(brand ~ 1 | income + age, data = mlogit_data)

# Kiểm định IIA (ví dụ, loại bỏ "Other")
# mlogit_model_subset <- mlogit(brand ~ 1 | income + age, data = mlogit_data, 
#                              alt.subset = c("Apple", "Samsung", "Google"))
# hmftest(mlogit_model, mlogit_model_subset)
# (Ghi chú: Phần này có tính minh họa, vì mlogit yêu cầu cấu trúc dữ liệu phức tạp hơn)

Case Study 5.5: Khi Giả định Proportional Odds bị Vi phạm

  • Bối cảnh: Trong Case Study 5.2, giả sử kiểm định Brant cho thấy giả định Proportional Odds bị vi phạm, đặc biệt là đối với biến company_size.

  • Nhiệm vụ:

    1. Diễn giải kết quả của kiểm định Brant.
    2. Đề xuất các hướng xử lý.
    3. Thử xây dựng các mô hình logistic nhị phân riêng biệt để xem hệ số của company_size thay đổi như thế nào ở các điểm cắt khác nhau.

  • Phân tích bằng R:

    # Giả sử brant(rating_model) cho p-value nhỏ
# 1. Diễn giải: Tác động của ít nhất một biến (ví dụ: company_size) lên
# odds của việc có hạng cao hơn là không nhất quán trên toàn bộ thang đo xếp hạng.

# 2. Hướng xử lý:
# - Sử dụng Partial Proportional Odds Model (gói `VGAM` với hàm `vglm`).
# - Sử dụng mô hình multinomial logit (bỏ qua thông tin thứ bậc).
# - Gộp các phạm trù của biến phản hồi lại.

# 3. Minh họa bằng các mô hình nhị phân
bond_data$rating_le_A <- ifelse(bond_data$rating <= "A", 1, 0)
bond_data$rating_le_BBB <- ifelse(bond_data$rating <= "BBB", 1, 0)

model_cut1 <- glm(rating_le_BBB ~ profit_margin + leverage_ratio + company_size,
                  data = bond_data, family = binomial)
model_cut2 <- glm(rating_le_A ~ profit_margin + leverage_ratio + company_size,
                  data = bond_data, family = binomial)

cat("Hệ số tại ngưỡng BBB:\n")
print(coef(model_cut1))
cat("\nHệ số tại ngưỡng A:\n")
print(coef(model_cut2))
# So sánh hệ số của 'company_size' trong hai mô hình để thấy sự khác biệt.

6.8 Bài tập

6.8.1 Bài tập lý thuyết

  1. Hãy cho 3 ví dụ về biến phản hồi danh nghĩa và 3 ví dụ về biến phản hồi thứ bậc trong lĩnh vực tài chính.
  2. Giải thích ý tưởng cơ bản của mô hình Baseline-Category Logit. Tại sao chúng ta cần \(J-1\) phương trình cho một biến phản hồi có \(J\) phạm trù?
  3. Trong mô hình multinomial logit, \(\exp(\beta_{jk})\) được diễn giải là gì? (với \(j\) là một phạm trù, \(k\) là một biến độc lập, và phạm trù cơ sở là \(J\)).
  4. Giả định IIA (Independence of Irrelevant Alternatives) là gì? Tại sao nó lại quan trọng và có thể là một vấn đề?
  5. Giải thích ý tưởng cơ bản của mô hình Ordinal Logit (Proportional Odds). Nó mô hình hóa đại lượng nào?
  6. Giả định Proportional Odds có nghĩa là gì? Tại sao nó còn được gọi là giả định “Parallel Lines”?
  7. Trong mô hình Ordinal Logit, ý nghĩa của các tham số \(\alpha_j\) (ngưỡng) và \(\beta_k\) (hệ số) là gì? Tại sao chỉ có một bộ \(\beta\)?
  8. Diễn giải ý nghĩa của \(\exp(\beta_k)\) trong mô hình Ordinal Logit.
  9. Nếu kiểm định Brant cho p-value = 0.02, bạn sẽ kết luận gì về mô hình Ordinal Logit của mình? Nêu ít nhất 2 hướng xử lý.
  10. So sánh hai cách tiếp cận: (1) sử dụng mô hình multinomial logit cho một biến phản hồi thứ bậc, và (2) sử dụng mô hình ordinal logit. Cách nào thường tốt hơn và tại sao?

6.8.2 Bài tập thực hành với R

  1. Phân tích dữ liệu iris:

    1. Sử dụng dữ liệu iris. Xây dựng mô hình multinom để dự báo Species dựa trên Sepal.Length và Sepal.Width.
    2. Chọn virginica làm phạm trù cơ sở và chạy lại mô hình. So sánh các hệ số với mô hình ở câu a.
    3. Tính toán và diễn giải RRR cho biến Sepal.Length đối với setosa so với virginica.

  2. Đánh giá mô hình multinom: Sử dụng mô hình từ bài 11a.

    1. Dự báo phạm trù cho toàn bộ tập dữ liệu iris.
    2. Tạo ma trận nhầm lẫn (confusion matrix).
    3. Tính độ chính xác tổng thể của mô hình.

  3. Phân tích dữ liệu housing của MASS:

    1. Sử dụng dữ liệu housing từ gói MASS. Xây dựng mô hình polr để dự báo Sat (mức độ hài lòng) dựa trên Infl (ảnh hưởng của ban quản lý) và Type (loại nhà cho thuê).
    2. In ra summary() của mô hình. Các biến độc lập có ý nghĩa thống kê không?
    3. Diễn giải ý nghĩa của các ngưỡng (intercepts).

  4. Diễn giải và Kiểm định mô hình polr: Sử dụng mô hình từ bài 13.

    1. Tính và diễn giải Odds Ratio cho InflMedium và InflHigh.
    2. Thực hiện kiểm định Brant cho mô hình này. Giả định Proportional Odds có được thỏa mãn không?

  5. Dự báo với polr: Sử dụng mô hình từ bài 13.

    1. Dự báo mức độ hài lòng có xác suất cao nhất cho những người thuê nhà loại “Apartment” với mức ảnh hưởng “High”.
    2. Dự báo xác suất cho mỗi mức độ hài lòng (Low, Medium, High) cho trường hợp trên. (Gợi ý: predict(..., type = "probs")).

  6. Mô hình hóa dữ liệu khảo sát: Sử dụng bộ dữ liệu survey từ gói MASS. Biến Smoke là thói quen hút thuốc (Heavy, Regul, Occas, Never), có thể coi là thứ bậc.

    1. Chuyển Smoke thành một factor có thứ tự.
    2. Xây dựng mô hình ordinal logit dự báo Smoke dựa trên Sex và Age.
    3. Kiểm tra giả định Proportional Odds.
    4. Dựa trên kết quả, bạn kết luận gì về mối liên hệ giữa tuổi, giới tính và thói quen hút thuốc?

  7. So sánh mô hình Ordinal và Multinomial: Vẫn sử dụng dữ liệu survey và biến phản hồi Smoke.

    1. Xây dựng một mô hình multinomial logit (multinom) để dự báo Smoke dựa trên Sex và Age (coi Smoke như biến danh nghĩa).
    2. So sánh AIC của mô hình này với mô hình ordinal logit từ bài 16. Mô hình nào tốt hơn? Tại sao?

  8. Mô hình hóa lựa chọn kinh tế: Sử dụng dữ liệu Heating từ gói mlogit. Dữ liệu này ở dạng “long”.

    1. Khám phá dữ liệu, tìm hiểu các biến. Biến depvar cho biết lựa chọn hệ thống sưởi.
    2. Xây dựng mô hình multinomial logit đơn giản bằng multinom để dự báo lựa chọn hệ thống sưởi dựa trên income.

  9. Tạo và phân tích dữ liệu thứ bậc:

    1. Tạo một data frame giả lập với biến phản hồi là “Mức độ rủi ro đầu tư” (Low, Medium, High) và các biến dự báo là tuoi_nha_dau_tu và kien_thuc_tc (thang điểm 1-10).
    2. Xây dựng mô hình polr.
    3. Vẽ biểu đồ các hiệu ứng (effects plot) bằng gói effects để trực quan hóa tác động của các biến độc lập lên xác suất của mỗi phạm trù. (Gợi ý: plot(allEffects(model))).

  10. Ma trận nhầm lẫn chi tiết: Sử dụng mô hình và các dự báo từ bài 14.

    1. Sử dụng hàm confusionMatrix() từ gói caret (install.packages("caret")) để có một ma trận nhầm lẫn chi tiết hơn, bao gồm độ nhạy, độ đặc hiệu cho từng lớp.
    2. Dựa vào kết quả, mô hình dự báo tốt nhất cho mức độ hài lòng nào?

  11. Mô hình hóa biến đếm (ôn tập GLM): Sử dụng dữ liệu housing. Xây dựng mô hình Poisson (glm) để dự báo Freq dựa trên Sat, Infl, Type, Cont. So sánh AIC của mô hình này với các mô hình polr đã xây dựng. (Lưu ý: đây là hai cách tiếp cận khác nhau, một bên coi Freq là biến phản hồi, một bên coi Sat là biến phản hồi).

  12. So sánh các mức trong multinom: Sử dụng mô hình iris_multinom_model từ phần lý thuyết. Làm thế nào để bạn tìm ra RRR của versicolor so với virginica? (Gợi ý: bạn cần tính lại mô hình với virginica làm cơ sở, hoặc thực hiện một phép biến đổi tuyến tính trên các hệ số).

  13. Xây dựng mô hình Adjacent-Categories (Nâng cao):

    1. Sử dụng gói VGAM. Tìm hiểu hàm vglm() với family = acat().
    2. Áp dụng nó cho dữ liệu housing để dự báo Sat.
    3. So sánh kết quả với mô hình polr.

  14. Trường hợp IIA bị vi phạm (Tư duy): Giả sử bạn đang mô hình hóa lựa chọn của người tiêu dùng giữa các loại nước giải khát: Cola, Pepsi, Nước cam. Bạn nghi ngờ rằng Cola và Pepsi là các sản phẩm thay thế gần gũi và có thể vi phạm IIA.

    1. Mô tả cách bạn sẽ kiểm tra giả định này.
    2. Nếu giả định bị vi phạm, bạn sẽ đề xuất mô hình nào thay thế? (Gợi ý: Nested Logit). Mô tả cấu trúc của mô hình đó (ví dụ, các “tổ” lựa chọn).

  15. Diễn giải hệ số của biến factor trong polr: Sử dụng mô hình housing_polr_model.

    1. Diễn giải hệ số của TypeTownhouse. Nó so sánh với mức cơ sở nào?
    2. Tính Odds Ratio cho nó. Diễn giải ý nghĩa của OR này.

6.9 Tài liệu tham khảo

  • Agresti, A. (2010). Analysis of Ordinal Categorical Data (2nd ed.). Wiley.
  • Agresti, A. (2018). An Introduction to Categorical Data Analysis (3rd ed.). Wiley.
  • Fox, J., & Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). Sage.
  • Long, J. S., & Freese, J. (2014). Regression Models for Categorical Dependent Variables Using Stata (3rd ed.). Stata Press. (Nhiều khái niệm có thể áp dụng cho R).
  • Powers, D. A., & Xie, Y. (2008). Statistical Methods for Categorical Data Analysis (2nd ed.). Emerald Group Publishing.
  • Venables, W. N., & Ripley, B. D. (2002). Modern Applied Statistics with S (4th ed.). Springer. (Tài liệu gốc cho gói nnet và MASS).
  • Christensen, R. H. B. (2019). Analysis of Ordinal Data with R. https://cran.r-project.org/web/packages/ordinal/vignettes/ordinal_tutorial.pdf (Vignette cho gói ordinal).
  • Các tài liệu hướng dẫn của gói R: nnet, MASS, ordinal, brant, mlogit, VGAM, effects.