NV1

PHẦN 1

TÓM TẮT CUỐN SÁCH “Generalized Linear Models With Examples in R” CỦA PETER K.DUNN VÀ GORDON K .SMYTH

Generalized Linear Models with Examples in R trình bày lý thuyết và ứng dụng mô hình GLM trong R, từ hồi quy tuyến tính, ANOVA đến các phần nâng cao về đại số ma trận và thuật toán. Cách trình bày linh hoạt, phù hợp cho nhiều đối tượng từ sinh viên đến chuyên gia thống kê.

CHƯƠNG 1: STATISTICAL MODELS - MÔ HÌNH THỐNG KÊ

Chương này bắt đầu với việc giới thiệu khái niệm mô hình thống kê, trong đó:

• y là biến kết quả(response variable)

• n là số quan sát

• \(x_1\), \(x_2\),…, \(x_p\) là các biến giải thích (explanatory variables)

Các biến giải thích có thể là định lượng (covariates) hoặc định tính (factors). Khi thực hiện phân tích data định lượng có thể dễ dàng hơn vì chúng được trình bài dưới dạng số học tuy nhiên dữ liệu định tính thì ngược lại. Để đưa các biến định tính vào mô hình, ta phải chuyển chúng thành dạng số thông qua các biến giả (dummy variables). Phương pháp mã hóa phổ biến là mã hóa treatment, với \(k-1\) biến giả cho một biến factor có k mức. Một cách đơn giản hơn, khi dữ liệu định tính gồm 3 mức độ đo lường sự hài lòng của khách hàng tuwg không thoải mái, bình thường và hài lòng ta sẽ đưa 2 biến giả vào mô hình.

Tại chương này cũng khuyến khích việc ứng dụng các biểu đồ để trình bày rõ hơn tại các bước đầu, tuy nhiên để đào sâu hơn vào cốt lõi của nó cần xây dựng các mô hình thống kê. Có 2 thành phần chính của mô hình thống kê:

• Thành phần hệ thống (systematic component): mô tả mối quan hệ giữa giá trị trung bình của biến kết quả với các biến giải thích hay mô tả kì vọng \(\mu_i = \mathbb{E}[y_i]\) của biến kết quả dựa trên các biến giải thích. Ví dụ, một mô hình đơn giản có dạng:

\[ \mu_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \beta_4 x_{4i} \]Trong đó:

\(x_{1i}\): tuổi của đối tượng thứ \(i\)

\(x_{2i}\): chiều cao (cm) của đối tượng thứ \(i\)

\(x_{3i}\): giới tính (biến giả; 1 = nam, 0 = nữ)

\(x_{4i}\): tình trạng hút thuốc (biến giả; 1 = có hút, 0 = không hút)

• Thành phần ngẫu nhiên (random component): mô tả sự biến thiên quanh giá trị trung bình \(\mu_i\). Có thể giả định phương sai không đổi \(\text{var}[y_i] = \sigma^2\), hoặc gỉa định phân phối chuẩn: \[ y_i \sim \mathcal{N}(\mu_i, \sigma^2) \]

Cần lưu ý rằng các giả định đơn giản như tuyến tính hay phương sai không đổi thường không phù hợp, vì dữ liệu có thể thể hiện mối quan hệ phi tuyến hoặc phương sai thay đổi theo giá trị trung bình. Do đó, cần cân nhắc các lựa chọn mô hình khác hoặc các phân phối khác ngoài phân phối chuẩn.

Mô hình thống kê là sự kết hợp giữa hai thành phần này để phản ánh đầy đủ đặc điểm của dữ liệu thực tế với mục đích cuối cùng là để dự báo và hiểu rõ mối quan hệ giữa các biến

Có 2 tiêu chí luôn cần được ưu tiên khi xây dựng mô hình hoàn chỉnh là đảm bảo được sự đơn giản và chính xác (parsimony and accuracy). Tức là mô hình đơn giản nhất nhưng vẫn mô tả đúng đặc điểm và sự biến thiên của dữ liệu, tránh việc giả thuyết chồng chéo và dày đặc các phương pháp khác nhau sẽ làm cho kết quả ước lượng khôn khớp với thực tế. Đây cũng chính là nguyên tắc “OCCAM’S RAZOR” trong kinh tế lượng.

Tuy nhiên, phương pháp thu thập dữ liệu ảnh hưởng rất lớn đến những kết luận có thể rút ra từ phân tích. . Trong nghiên cứu quan sát, dữ liệu chỉ được thu thập mà không can thiệp, nên mô hình chỉ thể hiện mối quan hệ tương quan, không đủ bằng chứng khẳng định được quan hệ nhân quả. Trong nghiên cứu thực nghiệm, khi người nghiên cứu kiểm soát biến giải thích, mô hình có thể hỗ trợ kết luận về quan hệ nhân quả.

Mặc dù các mô hình thống kê xử lý dữ liệu từ cả hai loại nghiên cứu này theo cách tương tự và kết luận thống kê có thể trông giống nhau, nhưng các kết luận khoa học rút ra từ thí nghiệm thường mạnh mẽ và đáng tin cậy hơn. Trong nghiên cứu quan sát, việc tốt nhất có thể làm là đo đạc và kiểm soát càng nhiều biến ngoại lai (extraneous variables) có thể ảnh hưởng đến biến phản hồi, nhằm điều chỉnh các ảnh hưởng không được kiểm soát trong thiết kế nghiên cứu.

Chương 2: Mô hình hồi quy tuyến tính

2.1. Mô hình cơ bản

Phản hồi \(y_i\) có phương sai \[ \mathrm{var}(y_i) = \sigma^2 w_i, \] với \(w_i\) là trọng số ưu tiên đã biết.

Kỳ vọng \[ E[y_i] = \mu_i \] tuyến tính theo biến giải thích: \[ \mu_i = \beta_0 + \sum_{j=1}^p \beta_j x_{ji}. \]

Các tham số \(\beta_j\) và \(\sigma^2\) cần ước lượng từ dữ liệu.

2.2 Các loại mô hình hồi quy tuyến tính

• Hồi quy tuyến tính đơn: chỉ có 1 biến giải thích.
• Hồi quy tuyến tính thường: tất cả \(w_i = 1\).
• Hồi quy tuyến tính có trọng số: các \(w_i\) khác nhau.
• Hồi quy tuyến tính bội: nhiều biến giải thích.
• Hồi quy tuyến tính chuẩn: giả định \[ y_i \sim N(\mu_i, \sigma^2 / w_i). \]

2.3 Ước lượng tham số

Dùng phương pháp bình phương tối thiểu có trọng số (WLS), tức là tối thiểu \[ S(\beta) = \sum_{i=1}^n w_i \left( y_i - \beta_0 - \sum_j \beta_j x_{ji} \right)^2. \]

Giải hệ phương trình đạo hàm để tìm \(\hat{\beta}\).

Ước lượng phương sai sai số: \[ s^2 = \frac{RSS}{n - p}, \quad \text{với } RSS = \sum w_i (y_i - \hat{\mu}_i)^2. \]

2.4 Độ không chắc chắn của ước lượng

Sai số chuẩn của \(\hat{\beta}_j\): \[ \mathrm{se}(\hat{\beta}_j) = s \sqrt{I_j^*}, \] trong đó \(I_j^*\) phản ánh mức độ độc lập của biến \(x_j\) với các biến khác.

Sai số chuẩn của giá trị dự báo \(\hat{\mu}_g\) tại điểm \(x_g\): \[ \mathrm{se}(\hat{\mu}_g) = \sqrt{x_g (X^T W X)^{-1} x_g^T \cdot s^2}. \]

2.5 Biểu diễn ma trận

Mô hình: \[ y = X \beta + \epsilon, \quad \mathrm{Var}(\epsilon) = \sigma^2 W^{-1}. \]

Ước lượng hệ số: \[ \hat{\beta} = (X^T W X)^{-1} X^T W y. \]

Ước lượng phương sai: \[ s^2 = \frac{(y - X \hat{\beta})^T W (y - X \hat{\beta})}{n - p}. \]

Ma trận phương sai–hiệp phương sai của \(\hat{\beta}\): \[ \mathrm{Var}(\hat{\beta}) = s^2 (X^T W X)^{-1}. \]

2.6 Suy luận thống kê: Kiểm định t

Giả định phân phối chuẩn cho sai số để sử dụng các kiểm định \(t\) và khoảng tin cậy.

Ước lượng \(\hat{\beta}_j\) có phân phối chuẩn với phương sai đã ước lượng.

Sai số chuẩn giúp đánh giá ý nghĩa thống kê và xây dựng khoảng tin cậy cho các hệ số.

CHƯƠNG 3: Linear Regression Models: Diagnostics and Model-Building

]Phát hiện và xử lý vi phạm giả định trong hồi quy tuyến tính bằng cách sử dụng phân tích chẩn đoán (diagnostics), đặc biệt thông qua phần dư (residuals).]{style=“color:blue”}

3.1 Types of Assumptions

1. Không có outlier (điểm ngoại lai) ảnh hưởng mạnh đến kết quả.
2. Quan hệ tuyến tính giữa biến phụ thuộc \(y\) và các biến giải thích.
3. Phương sai không đổi (homoscedasticity).
4. Độc lập giữa các quan sát.
5. Phân phối chuẩn của phần dư, đặc biệt trong hồi quy tuyến tính thường.

3.2: Các loại phần dư

• Phần dư thô (Raw residuals):

  \[ r_i = y_i - \hat{y}_i \]

• Phần dư chuẩn hóa (Standardized residuals):

  \[ r_i^* = \frac{r_i}{\hat{\sigma} \sqrt{1 - h_i}} \]

• Phần dư Student hóa (Studentized residuals):

  \[ t_i = \frac{r_i}{\hat{\sigma}_{(i)} \sqrt{1 - h_i}} \]

  Trong đó:

  • \(h_i\): leverage (đòn bẩy)
  • \(\hat{\sigma}_{(i)}\): ước lượng sai số chuẩn khi bỏ quan sát thứ \(i\)

3.3: Leverage (Đòn bẩy)

\[ h_i = \mathbf{x}_i^\top (X^\top X)^{-1} \mathbf{x}_i \]

• Đo lường mức độ ảnh hưởng của quan sát \(i\) đến dự đoán \(\hat{y}_i\).
• Nếu \(h_i\) lớn → quan sát có thể là điểm “đặc biệt” hoặc bất thường.

3.4: Kiểm tra giả định bằng biểu đồ phần dư

Giả định cần kiểm tra	Biểu đồ đề xuất
Tuyến tính	\(r_i^*\) vs. từng biến giải thích
Phương sai không đổi	\(r_i^*\) vs. \(\hat{y}_i\)
Phân phối chuẩn	Q-Q plot của phần dư
Độc lập giữa các quan sát	Phần dư theo thời gian hoặc vị trí

3.5: Phát hiện điểm ảnh hưởng lớn

• Cook’s Distance:

  \[ D_i = \frac{r_i^2}{p } \cdot \frac{h_i}{(1 - h_i)} \]

• DFFITS:

  \[ \text{DFFITS}_i = t_i \cdot \sqrt{\frac{h_i}{1 - h_i}} \]

• DFBETAS (với từng hệ số \(\beta_j\)):

  \[\text{DFBETAS}_{ij} = \frac{\hat{\beta}_j - \hat{\beta}_{j(i)}}{\text{SE}(\hat{\beta}_{j(i)})}\]

• CR:

  \[ CR = \frac{1}{1 - h} \cdot \left( \frac{n - p}{n - p + r^2} \right)^p \]

3.6: Khắc phục mô hình nếu có vi phạm giả định

• Vi phạm tính độc lập: sử dụng mô hình khác như:
  • GEE (Generalized Estimating Equations)
  • Mô hình hỗn hợp (Mixed Models)
• Phương sai thay đổi (Heteroscedasticity): biến đổi biến phụ thuộc \(y\):
  • Log: \(y' = \log(y)\)
  • Căn bậc hai: \(y' = \sqrt{y}\)
• Không tuyến tính:

  • Biến đổi biến giải thích

  • Thêm đa thức:

    \[ y = \beta_0 + \beta_1 x + \beta_2 x^2 + \dots + \beta_k x^k + \varepsilon \]

  • Hoặc sử dụng splines để mô hình hóa phi tuyến mềm dẻo hơn.

3.7: Xử lý outliers và điểm ảnh hưởng lớn

• Không nên tự động loại bỏ.
• Cần đánh giá lý do dữ liệu bất thường (do lỗi nhập liệu hay thực tế?).
• Có thể chạy mô hình 2 lần: có và không có điểm ảnh hưởng mạnh để so sánh kết quả.

3.8: Đa cộng tuyến (Collinearity)

• Xảy ra khi các biến giải thích có tương quan cao.
• Làm tăng phương sai ước lượng của hệ số \(\beta_j\), gây khó khăn trong diễn giải.

Phát hiện bằng:

• Variance Inflation Factor (VIF):

  \[ \text{VIF}_j = \frac{1}{1 - R_j^2} \]

  Trong đó \(R_j^2\) là hệ số xác định khi hồi quy \(x_j\) theo các biến giải thích còn lại.

CHƯƠNG 4: Beyond Linear Regression: The Method of Maximum Likelihood

]Giới thiệu phương pháp hợp lý cực đại (Maximum Likelihood - ML) để ước lượng mô hình, vượt ra ngoài hồi quy tuyến tính — chuẩn bị cho việc hiểu mô hình tuyến tính tổng quát (GLM).]{style=“color:blue”}

4.1 Khi Mô Hình Tuyến Tính Không Còn Phù Hợp

Những trường hợp mô hình tuyến tính thất bại:

• Tỉ lệ (0–1): dùng phân phối Binomial hoặc Bernoulli
• Số đếm: dùng Poisson hoặc Negative Binomial
• Dữ liệu dương, liên tục: dùng Gamma hoặc Inverse Gaussian

Phân phối chuẩn không phù hợp vì không đảm bảo đặc điểm miền giá trị và phương sai của dữ liệu.

4.2 Mô Hình Tuyến Tính Tổng Quát (GLM)

GLM là hệ thống mô hình mở rộng từ hồi quy tuyến tính với 3 thành phần chính:

1. Phân phối thành phần ngẫu nhiên: từ họ phân phối EDM
   
2. Thành phần hệ thống: dự báo tuyến tính
   \[ \eta = X \beta \]
3. Hàm liên kết:
   \[ g(\mu) = \eta \implies \mu = g^{-1}(\eta) \]

GLM thống nhất mô hình hóa cho nhiều loại dữ liệu.

4.3 Ước Lượng Hợp Lý Tối Đa (MLE)

MLE là tối đa hóa likelihood dựa trên dữ liệu quan sát.

Log-likelihood thường dùng hơn vì dễ tính, dễ đạo hàm.

Với phân phối chuẩn, MLE tương đương OLS.

Ví dụ với phân phối mũ:

\[ \hat{\theta} = \frac{1}{\bar{y}} \]

4.4 MLE Cho Một Tham Số

• Score function:
  \[ U(\zeta) = \frac{d\ell}{d\zeta} \]

• MLE: nghiệm của
  \[ U(\hat{\zeta}) = 0 \]

• Fisher Information:
  \[ I(\zeta) = \mathbb{E}[-\ell''(\zeta)] \]

• Sai số chuẩn:
  \[ SE(\hat{\zeta}) = \frac{1}{\sqrt{I(\hat{\zeta})}} \]

4.5 MLE Với Nhiều Tham Số

Mô hình tổng quát:

\[ \eta_i = \beta_0 + \beta_1 x_{i1} + \dots \]

Log-likelihood là tổng các log xác suất:

\[ \ell = \sum \log P(y_i; \mu_i) \]

Score equations:

\[ U(\beta_j) = \sum (y_i - \mu_i) x_{ij} \]

Giải hệ phương trình bằng các thuật toán lặp (Newton–Raphson, IRLS).

4.6 MLE Dùng Đại Số Ma Trận

• Hàm điểm số:
  \[ U(\beta) = X^T (y - \mu) \]

• Ma trận thông tin:
  \[ I(\beta) = X^T W X, \quad W = \text{diag}(\mu_i (1 - \mu_i)) \]

• Sai số chuẩn: là phần tử đường chéo của \(I^{-1}\).

4.7 Tính Chất Của MLE

Tính chất	Mô tả
Bất biến	Hàm của MLE là MLE của hàm đó
Không chệch tiệm cận	\(\mathbb{E}[\hat{\zeta}] \to \zeta\) khi \(n \to \infty\)
Hiệu quả tiệm cận	Không có ước lượng không chệch nào tốt hơn
Nhất quán	\(\hat{\zeta} \to \zeta\)
Phân phối chuẩn tiệm cận	\(\hat{\zeta} \sim \mathcal{N}(\zeta, 1/I(\zeta))\)

4.8 Kiểm Định Giả Thuyết Dựa Trên MLE

• Wald Test:

\[ W = \frac{(\hat{\zeta} - \zeta_0)^2}{\text{Var}(\hat{\zeta})} \sim \chi_1^2 \]

• Score Test:

\[ S = \frac{U(\zeta_0)^2}{I(\zeta_0)} \sim \chi_1^2 \]

• Likelihood Ratio Test (LRT):

\[ L = 2[\ell(\hat{\zeta}) - \ell(\zeta_0)] \sim \chi_1^2 \]

Cả ba kiểm định đều tiệm cận \(\chi^2\), dùng tốt khi \(n\) lớn.

4.9 Khoảng Tin Cậy (Confidence Intervals)

• Wald CI:

\[ \hat{\zeta} \pm z^* \cdot SE(\hat{\zeta}) \]

• Score CI và LRT CI: giải bất phương trình liên quan đến log-likelihood hoặc score.

Với nhiều tham số: miền tin cậy đa biến là hình elip trong không gian tham số.

4.10 So sánh mô hình không lồng nhau bằng AIC và BIC

Mục tiêu

So sánh các mô hình không lồng nhau bằng cách cân bằng giữa:

• Độ phù hợp (log-likelihood)
• Độ phức tạp (số tham số chưa biết)

→ Tránh overfitting khi chọn mô hình tốt nhất.

AIC – Akaike Information Criterion

\[ \mathrm{AIC} = -2 \ell(\hat{\zeta}_1, \ldots, \hat{\zeta}_p; y) + 2k \]

• \(\ell(\cdot)\): log-likelihood tại MLE
• \(k\): số tham số chưa biết

Phạt nhẹ cho mô hình phức tạp (hệ số phạt = 2).

AIC càng nhỏ → mô hình càng tốt

BIC – Bayesian Information Criterion

\[ \mathrm{BIC} = -2 \ell(\hat{\zeta}_1, \ldots, \hat{\zeta}_p; y) + k \log n \]

• \(n\): kích thước mẫu

Phạt nặng hơn AIC khi \(n\) lớn (hệ số phạt = \(\log n\)).

BIC càng nhỏ → mô hình càng tốt

CHƯƠNG 5: Generalized Linear Models: Structure

5.1 Giới thiệu và Tổng quan

Các chương trước (Chương 2 và 3) trình bày hồi quy tuyến tính và giả định phương sai không đổi — điều không phù hợp với nhiều loại dữ liệu thực tế (như đã thấy ở Chương 4).

GLMs mở rộng hồi quy tuyến tính bằng cách:

• Cho phép biến phụ thuộc theo nhiều phân phối xác suất khác nhau (không chỉ phân phối chuẩn).
• Thiết lập mối liên hệ phi tuyến giữa trung bình của biến phản hồi và tổ hợp tuyến tính của các biến giải thích thông qua hàm liên kết (link function).

Cấu trúc chương:

• 5.2: Hai thành phần chính của GLM
  
• 5.3–5.4: Thành phần ngẫu nhiên và mô hình phân tán hàm mũ (EDMs)
  
• 5.5: Thành phần hệ thống
  
• 5.6–5.7: Định nghĩa chính thức GLM và độ lệch (deviance)
  
• 5.8: So sánh GLM với hồi quy tuyến tính biến đổi

5.2 Hai Thành phần của GLM

(1) Thành phần ngẫu nhiên:

• Là phân phối xác suất của biến phản hồi.
• Lựa chọn dựa vào loại dữ liệu:
  • Nhị phân → Nhị thức (Binomial)
  • Đếm → Poisson
  • Liên tục dương → Gamma
• Mô tả mối quan hệ giữa trung bình và phương sai.

(2) Thành phần hệ thống:

• Biểu thức tuyến tính:
  \[ \eta = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p \]
• Trung bình của phản hồi \(\mu\) liên kết với \(\eta\) qua hàm liên kết \(g(\mu) = \eta\).

5.3 Thành phần Ngẫu nhiên: EDMs

• Liên tục: Normal, Gamma
  
• Rời rạc: Poisson, Binomial, Negative Binomial

Phân phối thuộc EDM có dạng:

\[ P(y; \theta, \phi) = a(y, \phi) \exp\left( \frac{y\theta - \kappa(\theta)}{\phi} \right) \]

Trong đó:

• \(\theta\): Tham số chuẩn tắc (canonical parameter)
  
• \(\kappa(\theta)\): Hàm tích lũy (cumulant function)
  
• \(\phi > 0\): Tham số phân tán
  
• \(a(y, \phi)\): Hàm chuẩn hóa

Ví dụ:

Phân phối	\(\theta\)	\(\phi\)	\(\kappa(\theta)\)
Normal	\(\mu\)	\(\sigma^2\)	\(\theta^2/2\)
Poisson	\(\log(\mu)\)	1	\(e^\theta = \mu\)
Binomial	\(\log\left( \frac{\mu}{1-\mu} \right)\)	\(1/m\)	—

Không phải mọi phân phối đều là EDM, ví dụ: Weibull (nếu \(\alpha \ne 1\)).

Hàm Sinh Moment (MGF) và Hàm Tích Lũy (CGF)

• MGF: \(M(t) = E[e^{ty}]\)
• CGF: \(K(t) = \log M(t)\)

Với EDM:

\[ K(t) = \frac{\kappa(\theta + t\phi) - \kappa(\theta)}{\phi} \]

Suy ra:

• Kỳ vọng: \(\mu = \kappa'(\theta) = E[y]\)
  
• Phương sai: \(\text{Var}(y) = \phi \kappa''(\theta)\)

Kết luận Tổng quát

\[ \mu = \frac{d\kappa(\theta)}{d\theta}, \quad \text{Var}(y) = \phi V(\mu) \]

Với:

\[ V(\mu) = \left( \frac{d\mu}{d\theta} \right)^{-1} \]

Ví dụ:

• Normal: \(V(\mu) = 1\) → phương sai không đổi
  
• Poisson: \(V(\mu) = \mu\) → phương sai bằng trung bình

5.4 Phân phối Phụ thuộc Trung bình

GLM gán phương sai theo trung bình:
\[ \text{Var}(y) = \phi V(\mu) \]
⇒ Khác với hồi quy tuyến tính, nơi phương sai không đổi.

5.5 Thành phần Hệ thống

5.5.1 Hàm Liên kết (Link Function)

• Giả định: \(g(\mu) = \eta\) với \(\eta = \beta_0 + \sum_{j=1}^p \beta_j x_j\)
• Hàm liên kết phải đơn điệu và khả vi.

Liên kết chuẩn tắc (canonical link):
\(\eta = \theta = g(\mu)\)

Chọn hàm liên kết phù hợp giúp:

• Giải thích dễ hơn
• Mô hình hóa hiệu quả hơn

5.5.2 Offset (Hệ số dịch chuyển)

• Là thành phần cố định không cần ước lượng.
• Ví dụ: điều chỉnh theo mức độ “phơi nhiễm” bằng \(\log(\text{dân số})\).
• Mô hình hóa:
  \[ \eta_i = o_i + \beta_0 + \sum_{j=1}^p \beta_j x_{ji} \]

5.6 Định nghĩa Chính thức của GLM

GLM bao gồm:

• Thành phần ngẫu nhiên:
  \[ y_i \sim \text{EDM}(\mu_i, \phi / w_i) \]

• Thành phần hệ thống:
  \[ \eta_i = o_i + \beta_0 + \sum_{j=1}^p \beta_j x_{ji} \]

• Hàm liên kết:
  \[ g(\mu_i) = \eta_i \]

Ký hiệu ngắn: glm(EDM; Link function)

5.7 Độ lệch Tổng (Deviance)

• Deviance đo độ phù hợp của mô hình:

  \[ D(y, \mu) = \sum w_i \cdot d(y_i, \mu_i) \]

• Scaled deviance:

  \[ D^*(y, \mu) = \frac{D(y, \mu)}{\phi} \]

• Nếu mô hình đúng, \(D^*\) xấp xỉ phân phối Chi-bình phương với \(n - p\) bậc tự do.

• Deviance liên quan đến log-likelihood.

5.8 Hồi Quy Xấp Xỉ GLM

Hồi quy tuyến tính dùng biến đổi (log, căn, …) để xử lý phương sai thay đổi nhưng:

• Chỉ là xấp xỉ
• Thiếu linh hoạt

GLM tốt hơn vì:

• Phân phối và hàm liên kết tách biệt
• Có thể làm việc trên thang đo gốc
• Dễ giải thích ảnh hưởng của biến giải thích đến trung bình \(\mu\)

CHƯƠNG 6: Generalized Linear Models: Estimation

6.1 Tính Toán Hàm Hợp Lý Cho β (Likelihood Calculations for β)

6.1.1 Đạo hàm hàm xác suất

Với một quan sát đơn lẻ
\(y \sim \text{EDM}(\mu, \phi / w)\), hàm log-xác suất có đạo hàm theo tham số tự nhiên \(\theta\) là:

\[ \frac{\partial \log P(y; \mu, \phi / w)}{\partial \theta} = \frac{w (y - \mu)}{\phi} \]

Dựa trên mối liên hệ giữa \(\theta\) và \(\mu\), đạo hàm theo \(\mu\) là:

\[ \frac{\partial \log P}{\partial \mu} = \frac{w (y - \mu)}{\phi V(\mu)} \]

Nếu hàm liên kết \(g(\mu) = \eta = \beta_0 + \sum_{j=1}^p \beta_j x_j\), thì đạo hàm theo \(\beta_j\) là:

\[ \frac{\partial \log P}{\partial \beta_j} = \frac{(y - \mu) w x_j}{\phi V(\mu)} \cdot \frac{d \mu}{d \eta} \]

6.1.2 Phương trình điểm và thông tin Fisher

Hàm log-likelihood tổng quát cho \(n\) quan sát:

\[ \ell(\beta_0, \ldots, \beta_p, \phi; y) = \sum_{i=1}^n \log P(y_i; \mu_i, \phi / w_i) \]

Phương trình điểm (score equations) cho từng \(\beta_j\) là:

\[ U(\beta_j) = \frac{1}{\phi} \sum_{i=1}^n W_i \frac{d \eta_i}{d \mu_i} (y_i - \mu_i) x_{ji} \]

trong đó trọng số làm việc (working weights) là:

\[ W_i = \frac{w_i}{V(\mu_i)} \left( \frac{d \mu_i}{d \eta_i} \right)^2 \]

Ma trận thông tin Fisher:

\[ I_{jk}(\beta) = \frac{1}{\phi} \sum_{i=1}^n W_i x_{ji} x_{ki} \]

---

6.2 Tính Toán Ước Lượng \(\beta\) (Computing Estimates of \(\beta\))

Thuật toán Fisher Scoring dùng để giải hệ phương trình điểm \(U(\beta_j) = 0\).
Mỗi bước lặp tương đương hồi quy bình phương nhỏ nhất có trọng số (weighted least squares) với biến phản hồi làm việc:

\[ z_i = \eta_i + \frac{d \eta_i}{d \mu_i} (y_i - \mu_i) \]

trên các biến giải thích \(x_{ji}\), dùng trọng số \(W_i\).

Thuật toán lặp quá trình cập nhật:

\[ \beta^{(r)} \to \eta^{(r)} \to \mu^{(r)} \to z^{(r)} \to \beta^{(r+1)} \]

Gọi thuật toán này là IRLS (Iteratively Reweighted Least Squares).

---

Điểm khởi đầu:

\[ \hat{\mu}_i^{(0)} = y_i \]

(có thể điều chỉnh nhỏ nếu \(y_i = 0\)).

6.3 Độ lệch phần dư (Residual Deviance)

Độ lệch đo lường phần biến thiên còn lại sau khi đã khớp mô hình:

\[ D(y, \hat{\mu}) = \sum_{i=1}^n w_i d(y_i, \hat{\mu}_i) \]

trong đó \(d(y, \mu)\) là độ lệch đơn vị (unit deviance).

Với phân phối Poisson:

\[ D(y, \hat{\mu}) = 2 \sum_{i=1}^n \left[ y_i \log\left(\frac{y_i}{\hat{\mu}_i}\right) - (y_i - \hat{\mu}_i) \right] \]

Nếu \(\phi \neq 1\), ta dùng scaled deviance:

\[ D^*(y, \hat{\mu}) = \frac{D(y, \hat{\mu})}{\phi} \]

Trong R, ta dùng deviance(fit) để lấy giá trị này từ mô hình glm.

6.4 Standard Errors for \(\hat{\beta}\)

Mục tiêu: Tính sai số chuẩn (standard errors - se) cho các ước lượng \(\hat{\beta}_j\) của hệ số hồi quy.

Phương pháp: Sai số chuẩn là căn bậc hai của các phần tử trên đường chéo ma trận nghịch đảo của ma trận thông tin (information matrix).

\[ se(\hat{\beta}_j) = \phi \, v_j \]

với \(v_j\) là căn bậc hai của phần tử trên đường chéo tương ứng trong ma trận nghịch đảo \((X^T W X)^{-1}\).

Nếu \(\phi\) (tham số phân tán) chưa biết, dùng giá trị ước lượng thay thế.

Ví dụ minh họa: Mô hình Poisson với dữ liệu noisy miner, kết quả trong R show estimate, std. error, z value.

6.5 Estimation of \(\beta\): Matrix Formulation

Công thức tổng quát cho vector score (đạo hàm log likelihood theo \(\beta\)):

\[ U = \frac{1}{\phi} X^T W M (y - \mu) \]

Ma trận thông tin Fisher:

\[ I = \frac{1}{\phi} X^T W X \]

Phương pháp Fisher scoring iteration để tìm \(\hat{\beta}\):

\[ \hat{\beta}^{(r+1)} = \hat{\beta}^{(r)} + I(\hat{\beta}^{(r)})^{-1} U(\hat{\beta}^{(r)}) \]

Biểu diễn lại bằng phương pháp Iteratively Weighted Least Squares (IWLS):

\[ \hat{\beta}^{(r+1)} = (X^T W X)^{-1} X^T W z \]

với

\[ z = \hat{\eta} + M (y - \hat{\mu}) \]

gọi là vector phản hồi công việc (working response).

Sau khi hội tụ, ma trận hiệp phương sai của \(\hat{\beta}\) được ước lượng bằng:

\[ \mathrm{var}[\hat{\beta}] = \phi (X^T W X)^{-1} \]

6.6 Estimation of GLMs Is Locally Like Linear Regression

Ý tưởng: Việc ước lượng GLM có thể coi là bài toán hồi quy tuyến tính trọng số, với vector phản hồi và trọng số cố định tại điểm hội tụ.

Các đại lượng như giá trị dự đoán \(\hat{\mu}\), phương sai \(\hat{\beta}_j\), độ nghiêng (leverages), residuals, Cook’s distance,… có thể tính tương tự như trong hồi quy tuyến tính cổ điển.

Điều này giúp tận dụng các kỹ thuật chuẩn của hồi quy tuyến tính để phân tích GLM.

6.7 Estimating \(\phi\)

6.7.1 Giới thiệu

Tham số \(\phi\) (tham số phân tán) không cần để ước lượng \(\beta\) nhưng cần thiết cho kiểm định giả thuyết và khoảng tin cậy.

Với các mô hình chuẩn (binomial, Poisson), \(\phi\) thường được gán cố định = 1.

6.7.2 Maximum Likelihood Estimator of \(\phi\)

Ước lượng ML cho \(\phi\) thường bị lệch, ví dụ trong hồi quy tuyến tính,

\[ \hat{\sigma}^2 = \frac{1}{n} \sum w_i (y_i - \hat{\mu}_i)^2 \]

bị lệch, do đó thường dùng ước lượng không lệch:

\[ s^2 = \frac{1}{n - p} \sum w_i (y_i - \hat{\mu}_i)^2 \]

6.7.3 Modified Profile Log-Likelihood Estimator of \(\phi\)

Phương pháp tinh vi hơn, dựa trên profile log-likelihood và có tính chất gần như không lệch và nhất quán.

Tuy nhiên, tính toán phức tạp, thường phải dùng phương pháp lặp.

6.7.4 Mean Deviance Estimator of \(\phi\)

Dựa trên trung bình deviance:

\[ \tilde{\phi} = \frac{D(y, \hat{\mu})}{n - p} \]

6.7.5 Pearson Estimator of \(\phi\)

Dựa trên thống kê Pearson:

\[ X^2 = \sum W_i (z_i - \hat{\eta}_i)^2 = \sum w_i \frac{(y_i - \hat{\mu}_i)^2}{V(\hat{\mu}_i)} \]

Ước lượng Pearson:

\[ \bar{\phi} = \frac{X^2}{n - p} \]

Ước lượng này rất phổ biến, có tính kháng (robust), được R sử dụng mặc định.

6.7.6 So sánh các ước lượng \(\phi\)

• ML estimator thường lệch.

• Modified profile log-likelihood tốt về mặt lý thuyết nhưng khó tính.

• Mean deviance và Pearson estimator dễ tính, Pearson estimator được dùng phổ biến vì ít giả định hơn.

Ví dụ với mô hình gamma và dữ liệu cây cherry được trình bày để minh họa.

6.8 Using R to Fit GLMs

Hàm chính trong R:

Cách sử dụng:

glm(formula, family = distribution(link = "link_function"), data = ...)

CHƯƠNG 7: Generalized Linear Models: Inference

7.1 Suy diễn cho hệ số khi φ đã biết

7.1.1 Kiểm định Wald cho từng hệ số hồi quy

• Giả thiết:
  \[ H_0: \beta_j = \beta_{j0} \quad (\text{thường là } 0) \]

• Thống kê kiểm định:
  \[ Z = \frac{\hat{\beta}_j - \beta_{j0}}{se(\hat{\beta}_j)} \]

• Khi \(H_0\) đúng, \(Z\) xấp xỉ phân phối chuẩn chuẩn hóa.

7.1.2 Khoảng tin cậy cho từng hệ số

• Khoảng tin cậy \(100(1-\alpha)\%\) cho \(\beta_j\):
  \[ \hat{\beta}_j \pm z_{\alpha/2} \times se(\hat{\beta}_j) \]

7.1.3 Khoảng tin cậy cho giá trị kỳ vọng \(\mu\)

• Dự đoán giá trị trung bình: \(\hat{\mu}\).

• Khoảng tin cậy được tính trên thang linear predictor \(\eta\), sau đó chuyển ngược về thang \(\mu\) bằng hàm link nghịch.

7.1.4 Kiểm định likelihood ratio để so sánh mô hình lồng nhau

• So sánh mô hình nhỏ (A) và mô hình lớn (B) với số tham số \(p_B > p_A\).

• Giả thiết:
  \[ H_0: \beta_{p_A+1} = \beta_{p_A+2} = \cdots = \beta_{p_B} = 0 \]

• Thống kê kiểm định:
  \[ L = \frac{D(y, \hat{\mu}^A) - D(y, \hat{\mu}^B)}{\phi} \]

  với \(D\) là deviance.

• Khi \(H_0\) đúng,
  \[ L \sim \chi^2_{p_B - p_A} \]

7.1.6 Kiểm định score

• Dùng để kiểm tra xem biến mới có nên thêm vào mô hình hay không.

• Thống kê score:
  \[ Z = \frac{\sum_i e(x_{p+1})_i \cdot e(y)_i}{\sqrt{\sum_i e(x_{p+1})_i^2}} \]

  với \(e(y)_i\) là residual làm việc,
  \(e(x_{p+1})_i\) là residual từ hồi quy trọng số của biến mới trên biến đã có.

7.1.7 Kiểm định score với ma trận

• Mở rộng kiểm định score để kiểm tra cùng lúc nhóm biến (k biến).

• Thống kê score dạng ma trận tương ứng.

7.2 Kết quả xấp xỉ lớn (Asymptotic Results)

• Nền tảng toán học cho các phân phối của kiểm định Wald, Score, và Likelihood Ratio.

• Giải thích các giả định mẫu lớn để các thống kê kiểm định xấp xỉ phân phối chuẩn hoặc \(\chi^2\).

7.3 Kiểm định Goodness-of-Fit

• Đánh giá liệu predictor tuyến tính có mô tả đủ xu hướng hệ thống trong dữ liệu không.

• Kiểm định dựa trên thống kê deviance hoặc Pearson chi-square.

• Phân phối các kiểm định goodness-of-fit dựa trên giả định phân tán nhỏ (small dispersion asymptotics).

7.4 Small Dispersion Asymptotics

• Mục đích: Giúp các kiểm định goodness-of-fit (deviance, Pearson) có hiệu quả ở mức quan sát đơn lẻ, không chỉ với mẫu lớn.

• Khác với large-sample asymptotics, small dispersion asymptotics xem xét trường hợp khi độ chính xác (precision) của từng quan sát tăng, tức là \(\phi \to 0\).

• Hai công cụ chính:

  • Saddlepoint approximation: tốt cho deviance statistics, hội tụ nhanh với tỉ lệ \(O(\phi)\).
    
  • Central Limit Theorem (CLT): áp dụng cho Pearson statistics, hội tụ chậm hơn với tỉ lệ \(O(\phi)\).

• Tiêu chí chính:
  \[ \tau = \phi \cdot V(y) \cdot (y - \text{boundary})^2 \leq \begin{cases} \frac{1}{3} & \text{cho Saddlepoint} \\ \frac{1}{5} & \text{cho CLT (nghiêm ngặt hơn)} \end{cases} \]

• Với các GLM phổ biến:

  • Binomial: Saddlepoint tốt nếu \(\min(m_i y_i) \geq 3\) và \(\min(m_i (1 - y_i)) \geq 3\); CLT tốt nếu \(\geq 5\).
    
  • Poisson: Saddlepoint tốt nếu \(\min(y_i) \geq 3\); CLT tốt nếu \(\geq 5\).
    
  • Gamma: Saddlepoint nếu \(\phi \leq \frac{1}{3}\); CLT nếu \(\phi \leq \frac{1}{5}\).

• Nếu không thỏa điều kiện, kiểm định deviance thường tốt hơn Pearson.

7.5 Suy diễn cho hệ số khi \(\phi\) chưa biết

7.5.1 Kiểm định Wald

• ước lượng \(\hat{\phi}\) để tính sai số chuẩn và thống kê Wald:
  \[ T = \frac{\hat{\beta}_j - \beta_{j0}}{se(\hat{\beta}_j)} \]

• Với mẫu lớn, \(T \sim N(0,1)\); với mẫu nhỏ hoặc vừa,
  \[ T \sim t_{n-p} \]

• R mặc định dùng ước lượng Pearson \(\hat{\phi}\).

7.5.2 Khoảng tin cậy cho \(\beta_j\)

• Khoảng tin cậy:
  \[ \hat{\beta}_j \pm t_{\alpha/2, n-p} \times se(\hat{\beta}_j) \]

• Hàm confint() trong R tính tự động.

7.5.3 Khoảng tin cậy cho giá trị dự đoán \(\mu\)

• Trên thang linear predictor:
  \[ \hat{\eta} \pm t_{\alpha/2, n-p} \times se(\hat{\eta}) \]

• Áp dụng hàm link nghịch để ra khoảng tin cậy trên \(\mu\).

7.5.4 Kiểm định likelihood ratio khi \(\phi\) chưa biết

• Dùng F-test:
  \[ F = \frac{\left[D(y, \hat{\mu}^A) - D(y, \hat{\mu}^B)\right] / (p_B - p_A)}{\hat{s}^2} \]

• \(\hat{s}^2\) là ước lượng phân tán \(\phi\) dựa trên mô hình lớn (B), có thể là: profile likelihood (\(\hat{\phi}^0\)), Pearson (\(\bar{\phi}\)), hoặc mean deviance (\(\tilde{\phi}\)).

• Với hồi quy tuyến tính (normal), phân phối F chính xác; với GLM khác, đây là xấp xỉ nhưng thường khá tốt.

• Nếu so sánh chỉ khác một hệ số, có thể dùng t-test với giá trị \(F^{1/2}\).

7.6 So sánh Wald, Score và Likelihood Ratio Tests

• Ba kiểm định phổ biến trong GLM: Wald test, Likelihood Ratio test (LRT), Score test.

• Wald test: thường dùng để kiểm định từng hệ số riêng lẻ, dễ hiểu, xuất hiện trong summary() của glm trong R.

• Likelihood Ratio test (LRT): dựa trên sự khác biệt của deviance, dùng trong hàm anova().

• Score test: ít dùng riêng, thường xuất hiện dưới dạng kiểm định goodness-of-fit Pearson. Nên dùng khi muốn kiểm tra thêm biến mới.

• Với hồi quy tuyến tính chuẩn, cả ba test có phân phối chính xác; với GLM, phân phối chỉ là xấp xỉ.

• Phân phối xấp xỉ của LRT và Score test thường tốt hơn Wald, nhất là với GLM Binomial, Poisson khi giá trị dự đoán gần ranh giới.

• Wald test không phù hợp khi hệ số ước lượng vô hạn; LRT vẫn chính xác hơn.

• Wald và Score test dùng cho kiểm định một hoặc hai phía; LRT truyền thống chỉ dùng kiểm định hai phía nhưng có thể áp dụng kiểm định một phía qua likelihood ratio dấu.

7.7 Chọn mô hình giữa các GLM không lồng nhau: AIC và BIC

• Các kiểm định trên chỉ áp dụng cho mô hình lồng nhau.

• Với mô hình không lồng nhau, dùng AIC và BIC để đánh giá.

• Không phải kiểm định thống kê chính thức, mà là tiêu chí chọn mô hình.

• Công thức khi \(\phi\) biết:
  \[ \mathrm{AIC} = -2 \log L + 2p \] \[ \mathrm{BIC} = -2 \log L + p \log n \]

• Khi \(\phi\) chưa biết:
  \[ \mathrm{AIC} = -2 \log L + 2(p+1) \] \[ \mathrm{BIC} = -2 \log L + (p+1) \log n \]

• Trong R: AIC(), BIC() tính theo công thức trên; extractAIC() tính AIC không hằng số cho hồi quy tuyến tính.

• Mô hình có AIC/BIC nhỏ hơn là tốt hơn.

7.8 Các phương pháp tự động chọn mô hình

• Các hàm R: drop1(), add1(), step().

• Dựa trên tiêu chí AIC để lựa chọn mô hình tối ưu.

• Cẩn trọng với GLM có ước lượng \(\phi\) vì \(\phi\) có thể khác nhau giữa các mô hình, khiến AIC chỉ là xấp xỉ.

• Phương pháp: loại dần (backward), thêm dần (forward), bước lặp (stepwise).

• Dù tự động, vẫn cần đánh giá lý thuyết và kiểm tra lại kết quả.

7.9 Sử dụng R để thực hiện các kiểm định trong GLM

• summary(fit): Bảng hệ số, sai số chuẩn, thống kê Wald, p-value, ước lượng \(\phi\), deviance, AIC, số lần lặp.

• deviance(fit): Trả về deviance của mô hình.

• df.residual(fit): Bậc tự do dư thừa.

• glm.scoretest(fit, x2): Kiểm định score cho biến mới (package statmod).

• anova(fit) hoặc anova(fit1, fit2, ...): So sánh mô hình lồng nhau, chọn test "F" (khi ước lượng \(\phi\)) hoặc "Chisq" (khi \(\phi\) biết).

• confint(fit, level=...): Khoảng tin cậy Wald cho hệ số.

• AIC(fit), BIC(fit): Trả về giá trị AIC, BIC.

• extractAIC(fit): Trả về AIC với penalty có thể chỉnh.

• drop1(fit), add1(fit): Loại hoặc thêm biến theo AIC, có thể chọn test "F" hoặc "Chisq".

• step(): Chọn mô hình tự động theo AIC, hỗ trợ các hướng forward, backward, both.

CHƯƠNG 8: Generalized Linear Models: Diagnostics

8.1 Giả định của Mô hình Tuyến tính Tổng quát (GLMs)

Các giả định chính khi xây dựng GLM gồm:

• Không có ngoại lệ (outliers): Tất cả quan sát phát sinh từ cùng một quá trình, phù hợp với cùng mô hình.
• Hàm liên kết (link function): Hàm liên kết được chọn đúng.
• Tính tuyến tính (linearity): Các biến giải thích quan trọng được đưa vào dưới dạng tuyến tính.
• Hàm phương sai (variance function): Hàm phương sai đúng được sử dụng.
• Tham số phân tán (dispersion parameter): Tham số \(\phi\) là hằng số.
• Tính độc lập (independence): Các phản hồi \(y_i\) độc lập.
• Phân phối (distribution): Các phản hồi tuân theo phân phối trong họ Exponential Dispersion Models (EDM).

Các giả định này được sắp xếp theo ảnh hưởng đến các moment của biến phản hồi: từ kỳ vọng (moment 1), phương sai (moment 2) đến phân phối đầy đủ. Mặc dù không luôn chính xác hoàn toàn, nhưng cần kiểm tra để đánh giá tính nhạy cảm của mô hình.

8.2 Residuals trong GLMs

8.2.1 Response Residuals không phù hợp cho GLM

Residual phản hồi:
\[ r_i = y_i - \hat{\mu}_i \]

Phương sai phụ thuộc vào kỳ vọng, nên residual này không đủ để đánh giá mô hình GLM.

8.2.2 Pearson Residuals

Để xử lý phương sai không đồng nhất, residual Pearson được định nghĩa:

\[ r_{P,i} = \frac{y_i - \hat{\mu}_i}{\sqrt{V(\hat{\mu}_i)/w_i}} \]

với \(V(\cdot)\) là hàm phương sai và \(w_i\) là trọng số. Residual này xấp xỉ chuẩn khi các điều kiện nhất định được thỏa mãn.

8.2.3 Deviance Residuals

Residual deviance là căn bậc hai có dấu của deviance đơn vị:

\[ r_{D,i} = \text{sign}(y_i - \hat{\mu}_i) \sqrt{w_i d(y_i, \hat{\mu}_i)} \]

với \(d(y_i, \hat{\mu}_i)\) là deviance đơn vị. Residual deviance cũng xấp xỉ chuẩn khi điều kiện phù hợp.

8.2.4 Quantile Residuals

Khi residual Pearson và deviance không chuẩn, đặc biệt với phân phối rời rạc, dùng residual quantile dựa trên hàm phân phối tích lũy (CDF):

\[ r_{Q,i} = \Phi^{-1}\big(F(y_i; \hat{\mu}_i, \phi)\big) \]

với \(\Phi^{-1}\) là hàm phân phối chuẩn ngược và \(F\) là CDF của phân phối phản hồi đã fit.

• Với phản hồi rời rạc, cần ngẫu nhiên hóa (randomized quantile residuals) để tạo tính liên tục, giúp residual này gần chuẩn hơn.

8.4 Leverages trong GLMs

8.4.1 Working Leverages

• GLM có thể xem như hồi quy tuyến tính tại từng điểm làm việc với working responses \(z_i\) và working weights \(W_i\) (phụ thuộc vào giá trị dự đoán \(\hat{\mu}_i\)).
• Leverage \(h_i\) thể hiện trọng số quan sát \(z_i\) trong tính toán giá trị tuyến tính dự đoán \(\hat{\eta}_i\).
• Nếu \(h_i\) nhỏ: nhiều quan sát cùng đóng góp cho \(\hat{\eta}_i\).
• Nếu \(h_i = 1\): giá trị dự đoán hoàn toàn do quan sát \(i\), tức \(\hat{\eta}_i = z_i\) và \(\hat{\mu}_i = y_i\).
• Phương sai residual làm việc \(e_i = z_i - \hat{\eta}_i\) xấp xỉ: \[ \mathrm{var}(e_i) \approx \phi V(\hat{\mu}_i)(1 - h_i) \]

8.4.2 Ma trận mũ (Hat Matrix)

\[ H = W^{1/2} X (X^T W X)^{-1} X^T W^{1/2} \]

• \(W\) là ma trận chéo trọng số tại vòng lặp cuối cùng.
• Leverages \(h_i\) là phần tử chéo của ma trận \(H\).
• Trong R, leverage tính bằng hàm hatvalues().

8.5 Residual Chuẩn hóa theo Leverage cho GLMs

• Residual thô gồm Pearson residual \(r_P\), deviance residual \(r_D\), và quantile residual \(r_Q\).
• Residual chuẩn hóa (standardized residuals) có phương sai gần bằng 1, được tính theo công thức: \[ r_P^{std} = \frac{r_P}{\sqrt{\phi(1 - h)}}, \quad r_D^{std} = \frac{r_D}{\sqrt{\phi(1 - h)}}, \quad r_Q^{std} = \frac{r_Q}{\sqrt{1 - h}} \] với \(h\) là leverage, \(\phi\) là tham số phân tán (dispersion parameter).

8.6 Khi nào nên dùng loại residual nào?

• Với phân phối chuẩn, các loại residual: Pearson, deviance, quantile đều gần chuẩn, có thể dùng thoải mái.
• Với phân phối rời rạc (discrete EDMs), nên dùng quantile residuals vì Pearson và deviance residuals thường không chuẩn, gây mẫu hình khó giải thích.
• Chuẩn hóa residuals giúp giảm biến đổi phương sai, nâng cao tính chính xác khi đánh giá mô hình.

8.7 Kiểm tra giả định mô hình

• Vẽ đồ thị residuals chuẩn hóa (standardized deviance hoặc quantile residuals) theo giá trị dự đoán \(\hat{\mu}\) hoặc theo từng biến giải thích để phát hiện xu hướng bất thường.
• Nếu xuất hiện xu hướng, cân nhắc chỉnh sửa mô hình: thay đổi hàm liên kết, thêm biến giải thích, hoặc biến đổi biến.
• Đồ thị working responses \(z_i\) theo predictor \(\hat{\eta}_i\) giúp kiểm tra tính đúng đắn của hàm liên kết.
• Sử dụng partial residual plots để kiểm tra xem biến giải thích đã được đưa vào mô hình đúng dạng tuyến tính chưa.

8.8 Outliers and Influential Observations (Quan sát ngoại lai và có ảnh hưởng)

Mục đích:
Phát hiện các điểm dữ liệu có ảnh hưởng lớn tới ước lượng mô hình.

Các chỉ số quan trọng:

• Deviance residuals \(r_i\): đánh giá sự khác biệt giữa quan sát và mô hình.
• Leverage \(h_i\): đại diện cho độ ảnh hưởng của điểm \(i\) trong không gian biến giải thích.
• Cook’s distance \(D_i\): đánh giá ảnh hưởng tổng thể của điểm \(i\) đến các hệ số hồi quy.

Công thức Cook’s distance:

\[ D_i = \frac{(r_i)^2 h_i}{p (1 - h_i)^2} \]

Trong đó, \(p\) là số tham số trong mô hình.

Cách dùng:
Điểm có Cook’s distance lớn (ví dụ \(D_i > \frac{4}{n}\)) có thể là điểm ảnh hưởng cần kiểm tra.

8.9 Residuals and Diagnostics (Phần dư và chẩn đoán)

Các loại phần dư:

• Raw residuals:
  \[ y_i - \hat{\mu}_i \]

• Pearson residuals:
  \[ r_i^{(P)} = \frac{y_i - \hat{\mu}_i}{\sqrt{\hat{V}(\hat{\mu}_i)}} \]

  với \(\hat{V}(\hat{\mu}_i)\) là phương sai ước lượng của \(y_i\).

• Deviance residuals:
  \[ r_i^{(D)} = \operatorname{sign}(y_i - \hat{\mu}_i) \sqrt{2 \left[ l(y_i; y_i) - l(y_i; \hat{\mu}_i) \right]} \]

  với \(l(\cdot)\) là hàm log-likelihood.

Mục đích:
Phát hiện sự bất thường trong dữ liệu, kiểm tra giả định mô hình.

8.10 Overdispersion and Quasi-likelihood (Phân tán vượt mức và quasi-likelihood)

Overdispersion xảy ra khi phương sai quan sát lớn hơn phương sai mô hình GLM dự đoán.

Ví dụ: Trong mô hình Poisson, giả định \(\operatorname{Var}(Y) = \mu\). Nếu dữ liệu có \(\operatorname{Var}(Y) > \mu\) gọi là overdispersion.

Giải pháp:
Sử dụng mô hình quasi-likelihood để điều chỉnh phương sai:

\[ \operatorname{Var}(Y_i) = \phi V(\mu_i) \]

với \(\phi > 1\) là hệ số phân tán (dispersion parameter).

Ước lượng \(\phi\):

\[ \hat{\phi} = \frac{1}{n-p} \sum_{i=1}^n \frac{(y_i - \hat{\mu}_i)^2}{V(\hat{\mu}_i)} \]

Sử dụng quasi-likelihood để điều chỉnh sai số chuẩn và kiểm định giả thuyết.

8.11 Collinearity (Đa cộng tuyến)

Khái niệm:
Các biến giải thích trong mô hình bị tương quan cao, gây khó khăn cho ước lượng.

Kiểm tra đa cộng tuyến:
Variance Inflation Factor (VIF):

\[ \text{VIF}_j = \frac{1}{1 - R_j^2} \]

Trong đó, \(R_j^2\) là hệ số xác định của hồi quy biến \(j\) trên các biến còn lại.

Ảnh hưởng:
VIF lớn (thường > 10) cho thấy biến đó bị đa cộng tuyến cao.

Ảnh hưởng tiêu cực đến độ chính xác ước lượng và kiểm định tham số.

Cách xử lý:
- Loại bỏ hoặc hợp nhất các biến có đa cộng tuyến cao.
- Sử dụng phương pháp hồi quy thu nhỏ như Ridge hoặc Lasso.

CHƯƠNG 9: GModels for Proportions: Binomial GLMs

9.1 Modelling Proportions

Mô hình hóa dữ liệu nhị phân hoặc dữ liệu tỷ lệ (ví dụ: số người hồi phục / tổng số).

• Dữ liệu nhị phân (binary): giá trị 0/1 (ví dụ, sống/chết).
• Dữ liệu tỷ lệ (proportion): số lần thành công / số lần thử.

Áp dụng mô hình tuyến tính tổng quát (GLM) với phân phối nhị thức (binomial), sử dụng các hàm liên kết như logit, probit, cloglog.

Cú pháp trong R:

glm(y/n ~ predictors, family = binomial)
# hoặc
glm(cbind(successes, failures) ~ x, family = binomial)

9.2 Link Functions

Các hàm liên kết phổ biến trong mô hình nhị thức:

• Logit: \[ \log\left(\frac{\pi}{1 - \pi}\right) \]

• Probit: dựa trên phân phối chuẩn tích lũy: \[ \Phi^{-1}(\pi) \]

• Complementary log-log (cloglog): \[ \log(-\log(1 - \pi)) \]

Lựa chọn hàm liên kết ảnh hưởng đến mô hình hóa và cách diễn giải kết quả.

9.3 Tolerance Distributions and the Probit Link

Hàm liên kết probit thích hợp khi mỗi cá thể có một ngưỡng phản ứng (tolerance) theo phân phối chuẩn.

• Áp dụng trong các lĩnh vực như sinh học, dược lý học, để mô hình hóa xác suất phản ứng theo liều lượng.

9.4 Odds, Odds Ratios and the Logit Link

• Odds (tỷ lệ cược): \[ \text{odds} = \frac{\pi}{1 - \pi} \]

• Odds ratio: tỉ số giữa odds của hai nhóm khác nhau, giúp so sánh mức độ ảnh hưởng.

• Logit link chuyển odds thành biểu thức tuyến tính: \[ \log\left(\frac{\pi}{1 - \pi}\right) \]

Hàm logit phổ biến do dễ diễn giải, đặc biệt trong thống kê y học và xã hội học.

9.5 Median Effective Dose, ED50

ED50 là liều lượng mà tại đó: \[ \pi = 0.5 \]

• Tức là 50% đối tượng có phản ứng với liều đó.
• Xác định từ mô hình bằng cách giải \(\pi = 0.5\).

Ứng dụng: đánh giá hiệu quả thuốc, nghiên cứu độc học.

9.6 The Complementary Log-Log Link in Assay Analysis

Hàm liên kết cloglog phù hợp khi xác suất thành công tăng nhanh ở mức liều cao (phân phối bất đối xứng):

\[ \log(-\log(1 - \pi)) \]

• Thường dùng trong các nghiên cứu về tử vong, phản ứng mạnh, hoặc dữ liệu thời gian sống (survival data).

9.7 Overdispersion

Quá phân tán (Overdispersion): khi phương sai lớn hơn kỳ vọng từ mô hình nhị thức.

Nguyên nhân:

• Thiếu biến giải thích.
• Sự phụ thuộc giữa các quan sát.
• Mẫu không đồng nhất.

Giải pháp:

• Sử dụng mô hình quasibinomial:

glm(y/n ~ x, family = quasibinomial)

9.8 When Wald Tests Fail

Wald test có thể không đáng tin cậy trong các trường hợp:

• Ước lượng nằm gần ranh giới (ví dụ: \(\pi \approx 0\) hoặc \(\pi \approx 1\)).
• Cỡ mẫu nhỏ.

Trong những tình huống này, Wald test có thể đánh giá sai ý nghĩa thống kê của hệ số.

Giải pháp: sử dụng kiểm định Likelihood Ratio Test (LRT), đáng tin cậy hơn vì:

• So sánh giữa mô hình đầy đủ và mô hình rút gọn.
• Dựa vào sự khác biệt về log-likelihood giữa hai mô hình.

9.9 No Goodness-of-Fit for Binary Responses

Với dữ liệu phản hồi nhị phân (0/1), không thể áp dụng các kiểm định goodness-of-fit truyền thống như Pearson hoặc deviance test vì:

• Không có nhiều quan sát lặp lại trong cùng một nhóm dự đoán (predictor combination).

Thay vào đó, có thể:

• Vẽ biểu đồ residuals (ví dụ: Pearson hoặc deviance residuals) để đánh giá sự phù hợp của mô hình.
• Gom nhóm dữ liệu theo predictor để tạo các nhóm có nhiều quan sát rồi mới kiểm định goodness-of-fit.
• Kiểm tra ảnh hưởng (influence):
  • Dùng influence.measures(), dfbeta, cooks.distance() trong R.
  • Xác định những điểm ngoại lai hoặc có ảnh hưởng mạnh đến mô hình.

Lưu ý: Việc gom nhóm cần thực hiện cẩn thận để tránh làm sai lệch dữ liệu gốc.

CHƯƠNG 10: Models for Counts: Poisson and Negative Binomial GLMs

10.1 Tổng quan về Poisson GLMs

Phân phối Poisson được định nghĩa bởi:

\[ Y_i \sim \text{Poisson}(\mu_i) \]

với \(\mu_i\) là kỳ vọng (mean) và cũng là phương sai.

Trong Poisson GLM:

Hàm liên kết chuẩn là logarithmic:

\[ \log(\mu_i) = \eta_i = \mathbf{x}_i^\top \boldsymbol{\beta} \]

Tức là:

\[ \mu_i = \exp(\mathbf{x}_i^\top \boldsymbol{\beta}) \]

Đặc điểm quan trọng:

• Mô hình đảm bảo \(\mu_i > 0\) (vì hàm mũ luôn dương).
• Phù hợp với các phản hồi dạng đếm không âm.
• Ước lượng các tham số bằng phương pháp tối đa hóa hợp lý (MLE) giống như trong các GLM khác.

Kiểm tra phân tán (dispersion):

• Nếu dữ liệu có phương sai lớn hơn kỳ vọng (overdispersion), mô hình Poisson có thể không còn phù hợp.
• Trong trường hợp đó, chuyển sang mô hình Negative Binomial là hợp lý hơn.

10.2 Mô hình hóa tỷ lệ (Modelling Rates)

Dữ liệu đếm thường được quan sát trên các khoảng thời gian hoặc không gian khác nhau, vì vậy cần điều chỉnh để so sánh công bằng.

Tỷ lệ (rate) được định nghĩa là:

\[ \text{Rate}_i = \frac{Y_i}{t_i} \]

trong đó \(t_i\) là đơn vị phơi nhiễm (exposure), chẳng hạn như thời gian, diện tích, hoặc số người.

Trong GLM, phơi nhiễm được đưa vào mô hình dưới dạng offset:

\[ \log(\mu_i) = \log(t_i) + \mathbf{x}_i^\top \boldsymbol{\beta} \]

hoặc tương đương:

\[ \log\left(\frac{\mu_i}{t_i}\right) = \mathbf{x}_i^\top \boldsymbol{\beta} \]

Điều này mô hình hóa \(\log(\text{rate})\) như một hàm tuyến tính của các biến giải thích.

Lưu ý:

• Offset là một thành phần cố định, không được ước lượng.
• Kỹ thuật này đảm bảo mô hình hóa chính xác các tỷ lệ xảy ra sự kiện trên các mức độ phơi nhiễm khác nhau.

10.3 Contigency Tables: Log-Linear Models

10.3.1 Giới thiệu chung về Bảng đối chiếu (Contingency Tables)

Bảng đối chiếu thường dùng để thể hiện dữ liệu đếm (count data) được phân loại theo các yếu tố phân loại (factors).

Bắt đầu với 2 yếu tố phân loại (bảng 2 chiều) rồi mở rộng đến 3 chiều và cao hơn.

Mục tiêu là phân tích mối quan hệ giữa các yếu tố phân loại qua các mô hình thống kê phù hợp.

10.3.2 Bảng 2 chiều: Thành phần hệ thống (Systematic Component)

Giả sử 2 yếu tố phân loại: \(A\) với \(I\) mức và \(B\) với \(J\) mức, ta có bảng \(I \times J\).

Ký hiệu:

• \(y_{ij}\): số đếm quan sát tại ô \((i,j)\)
• \(\mu_{ij}\): số đếm kỳ vọng (expected count) tại ô \((i,j)\)
• \(\pi_{ij}\): xác suất kỳ vọng quan sát thuộc ô \((i,j)\), với \(\mu_{ij} = m \pi_{ij}\), trong đó \(m\) là tổng số quan sát.
• \(m_{i \cdot}\): tổng hàng \(i\), \(m_{\cdot j}\): tổng cột \(j\).

Nếu \(A\) và \(B\) độc lập, ta có:

\[ \pi_{ij} = \pi_{i \cdot} \times \pi_{\cdot j} \]

Lấy log của kỳ vọng ta có công thức mô hình hệ thống:

\[ \log \mu_{ij} = \log m + \log \pi_{i \cdot} + \log \pi_{\cdot j} \]

Ví dụ cụ thể về bảng \(2 \times 2\) (về thái độ người Úc đối với thực phẩm biến đổi gen theo thu nhập) minh họa cách tạo và xử lý dữ liệu trong R.

10.3.3 Bảng 2 chiều: Thành phần ngẫu nhiên (Random Components)

10.3.3.1 Giới thiệu

Cách thức thu thập mẫu ảnh hưởng đến mô hình thống kê phù hợp.

Có ba kịch bản lấy mẫu phổ biến:

• Không cố định tổng hàng/cột (Không cố định các tổng biên): quan sát đến tự nhiên, không giới hạn.
• Tổng quan sát cố định: tổng số quan sát \(m\) cố định, phân bố theo multinomial.
• Tổng hàng hoặc tổng cột cố định: ví dụ như tổng số người thu nhập cao cố định, phân bố multinomial theo từng hàng hoặc cột.

10.3.3.2 Kịch bản không cố định tổng biên (No Marginal Totals Fixed)

Tổng quan sát \(m\) cũng là ngẫu nhiên (Poisson).

Mỗi ô \(y_{ij}\) có phân phối Poisson độc lập.

Hàm log-likelihood:

\[ l(\mu; y) = \sum_{i=1}^2 \sum_{j=1}^2 \left(-\mu_{ij} + y_{ij} \log \mu_{ij} \right) \]

Dữ liệu ví dụ có thể được phân tích bằng mô hình Poisson GLM với công thức:

\[ \log \hat{\mu}_{ij} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 \]

với biến giả \(x_1, x_2\) tương ứng với các mức của 2 yếu tố phân loại.

Mô hình tương tác cũng được xét để kiểm tra mối liên hệ giữa các yếu tố.

10.3.3.3 Kịch bản tổng quan sát cố định (Grand Total Fixed)

Tổng số quan sát \(m\) cố định.

Phân phối phù hợp là multinomial.

Hàm log-likelihood:

\[ l(\mu; y) = \sum_{i=1}^2 \sum_{j=1}^2 y_{ij} \log \mu_{ij} \]

Ràng buộc:

\[ \sum_{i=1}^2 \sum_{j=1}^2 \mu_{ij} = m \]

Mặc dù mô hình multinomial khác mô hình Poisson, nhưng nếu ta thêm hệ số hằng số (constant term) trong mô hình Poisson, mô hình Poisson có thể dùng để suy diễn như mô hình multinomial.

Điều này có ý nghĩa là mô hình Poisson với hằng số được coi như điều kiện trên tổng số quan sát \(m\).

10.3.3.4 Kịch bản tổng cột (hoặc tổng hàng) cố định (Column or Row Totals Fixed)

Giả sử tổng số người thuộc mỗi cột (hoặc hàng) cố định.

Phân phối multinomial áp dụng riêng cho từng cột (hoặc hàng).

Xác suất:

\[ P(\text{counts}) = \prod_{j=1}^2 \frac{m_{\cdot j}!}{\prod_i y_{ij}!} \prod_i \mu_{ij}^{y_{ij}} \]

Ràng buộc:

\[ \sum_i \mu_{ij} = m_{\cdot j} \quad \forall j \]

Mô hình Poisson có thể dùng lại nếu có các tham số đại diện cho các tổng hàng (hoặc cột) trong phần tử hệ thống (linear predictor).

10.3.4 Nghịch lý Simpson (Simpson’s Paradox)

Khi hợp nhất bảng đếm theo một yếu tố không hợp lý, có thể dẫn tới kết luận sai.

Ví dụ dữ liệu sỏi thận: nếu không phân tầng theo kích thước sỏi, kết luận cho rằng phương pháp B thành công hơn phương pháp A (thành công tổng thể cao hơn).

Tuy nhiên khi phân tầng riêng từng kích thước, lại thấy phương pháp A thành công hơn cả với sỏi nhỏ và sỏi lớn.

Đây chính là nghịch lý Simpson, do việc hợp nhất dữ liệu làm mất thông tin quan trọng về sự tương tác giữa các yếu tố.

10.3.5 Tương đương giữa GLM nhị thức và Poisson (Equivalence of Binomial and Poisson GLMs)

Trong một số trường hợp (số thử nghiệm lớn, tỷ lệ thành công nhỏ), phân phối nhị thức có thể được xấp xỉ bằng phân phối Poisson.

Ví dụ phân tích dữ liệu thái độ về GMO với mô hình nhị thức có logistic link và mô hình Poisson log-linear cho kết quả tương tự.

Ví dụ trên dữ liệu sỏi thận cũng được phân tích bằng mô hình nhị thức GLM (logistic) và cho kết luận tương đồng với mô hình bảng đếm.

10.3.6 Bảng đa chiều (Higher-Order Tables)

Với bảng có hơn 3 yếu tố, phân tích phức tạp hơn nhưng vẫn khả thi với R.

Ví dụ về nghiên cứu thanh thiếu niên với 4 yếu tố: tuổi, nhóm, giới tính, mức độ trầm cảm.

Mô hình phù hợp bao gồm các tương tác 2 chiều, 3 chiều (ví dụ tương tác giữa tuổi, trầm cảm, giới tính).

Phân tích thể hiện sự khác biệt trong tỷ lệ trầm cảm cao giữa các nhóm tuổi theo giới tính (ví dụ nữ giảm trầm cảm theo tuổi, nam giảm rồi tăng lại).

Việc đơn giản hóa bảng dữ liệu bằng cách gộp các nhóm có thể gây hiểu sai.

10.3.7 Giá trị bằng 0 cấu trúc trong bảng đếm (Structural Zeros in Contingency Tables)

Trong bảng đếm có thể xuất hiện ô có giá trị 0.

• Zero ngẫu nhiên (sampling zeros): do mẫu nhỏ, có thể xuất hiện giá trị 0 tạm thời, có thể dùng mô hình để dự đoán giá trị kỳ vọng.

• Zero cấu trúc (structural zeros): do bản chất không thể xảy ra, ví dụ nữ không thể bị ung thư tuyến tiền liệt.

Với structural zeros, phải xử lý đặc biệt, loại bỏ các ô này trước khi phân tích.

Ví dụ bảng ung thư theo giới tính ở Tây Úc: ung thư tuyến tiền liệt (nam), ung thư cổ tử cung (nữ) là zero cấu trúc; ung thư vú ở nam là zero ngẫu nhiên (rất hiếm).

10.4 Overdispersion trong Poisson GLMs

Overdispersion xảy ra khi phương sai của dữ liệu đếm lớn hơn giá trị trung bình (\(\mu\)), trái với giả định Poisson là \(\mathrm{Var}(y) = \mu\).

Nguyên nhân:

• Biến động ngẫu nhiên trong \(\mu\) dù các biến giải thích cố định.
• Các sự kiện đếm có tương quan dương, thường xảy ra khi các sự kiện xuất hiện thành cụm.

Hậu quả:

• Ước lượng tham số (\(\beta\)) có thể không bị ảnh hưởng nhiều, nhưng sai số chuẩn (standard errors) thường bị đánh giá thấp.
• Kết quả kiểm định có thể cho thấy ý nghĩa giả thuyết sai lệch (quá mức).

Phát hiện:

• Dùng kiểm định độ phù hợp mô hình (goodness-of-fit) như deviance hoặc Pearson chi-square.
• Khi giá trị kiểm định lớn hơn nhiều so với bậc tự do dư (residual df), hoặc sau khi loại bỏ biến giải thích và ngoại lệ, overdispersion là nguyên nhân.
• Khi giá trị đếm nhỏ, việc phát hiện overdispersion có thể khó do xấp xỉ chuẩn (asymptotic) không chính xác.

10.5.1 Negative Binomial GLMs

Một cách xử lý overdispersion là dùng mô hình phân phối âm bản (Negative Binomial), coi \(\mu\) cũng là biến ngẫu nhiên theo phân phối gamma.

• Phương sai:
  \[ \mathrm{Var}(y) = \mu + \psi \mu^2 \] trong đó \(\psi > 0\) là mức overdispersion.

• Tham số \(k = \frac{1}{\psi}\) được ước lượng cùng với các hệ số \(\beta\) bằng phương pháp Maximum Likelihood (ví dụ hàm glm.nb() trong R).

• So với Poisson, Negative Binomial thường cho sai số chuẩn lớn hơn và phù hợp hơn với dữ liệu có overdispersion.

10.5.2 Quasi-Poisson Models

• Quasi-Poisson giữ nguyên hàm phương sai dạng Poisson \(V(\mu) = \mu\) nhưng thêm hệ số tán xạ (dispersion parameter) \(\phi > 1\).
• Phương sai thực tế là:
  \[ \mathrm{Var}(y) = \phi \mu \]
• Phương pháp này dễ dùng, phổ biến, và vẫn cho ước lượng nhất quán và sai số chuẩn chính xác.
• Không dựa trên phân phối xác suất chuẩn, nên không thể dùng trực tiếp để tính likelihood, nhưng dùng để ước lượng và kiểm định phi tham số.

CHƯƠNG 11: Positive Continuous Data: Gamma and Inverse Gaussian GLMss

11.1 Mô hình hóa dữ liệu dương liên tục

Nhiều ứng dụng có biến phản hồi dương và liên tục, thường phân bố lệch phải vì ranh giới ở 0 giới hạn đuôi trái. Nếu biến này biến thiên nhiều bậc, độ lệch là tất yếu. Biến ràng buộc ở 0 cũng khiến phương sai thường tiến về 0 khi giá trị trung bình tiến về 0 (Mục 4.2). Dữ liệu dương liên tục thường có mối quan hệ phương sai–trung bình tăng dần.

Ngoài \[ V(\mu) = \mu \] đã quen với dữ liệu đếm, hai hàm phương sai tăng đơn giản nhất là \[ V(\mu) = \mu^2 \] và \[ V(\mu) = \mu^3 \] , tương ứng phân phối gamma và inverse Gaussian. Do đó, GLM dựa trên gamma và inverse Gaussian hữu ích với dữ liệu dương liên tục. Phân phối gamma ứng với dữ liệu tỷ lệ có hệ số biến thiên không đổi. Mô hình gamma trong R dùng family = Gamma(), inverse Gaussian dùng family = inverse.gaussian().

11.2 Phân phối Gamma

Hàm mật độ gamma:

\[ P(y; \alpha, \beta) = \frac{y^{\alpha - 1} e^{-y/\beta}}{\Gamma(\alpha) \beta^\alpha} \]

với \(y > 0\), \(\alpha > 0\) (tham số hình dạng), \(\beta > 0\) (tham số tỉ lệ), và

\[ E[y] = \alpha \beta, \quad Var[y] = \alpha \beta^2. \]

Viết lại theo \(\mu, \phi\):

\[ P(y; \mu, \phi) = \frac{y^{\frac{1}{\phi} - 1} e^{-\frac{y}{\phi \mu}}}{(\phi \mu)^{\frac{1}{\phi}} \Gamma\left(\frac{1}{\phi}\right)}. \]

Phương sai hàm gamma là

\[ V(\mu) = \mu^2. \]

Hệ số biến thiên (coefficient of variation) bằng tỉ lệ phương sai trên bình phương trung bình và giữ không đổi. Vì thế, gamma GLM phù hợp khi hệ số biến thiên xấp xỉ không đổi.

Ví dụ 11.1:
Dữ liệu sinh trưởng biomass cây Tilia cordata được mô tả trong bảng 11.1 (dữ liệu lime). Mối quan hệ được mô hình hóa giữa biomass lá và đường kính thân cây (DBH) cũng như tuổi cây. Các biểu đồ trong Fig. 11.1 cho thấy biomass luôn dương và phương sai tăng theo giá trị trung bình.

Ví dụ 11.2:
Phân nhóm theo tuổi và nguồn gốc, tính phương sai và trung bình nhóm, đồ thị \(\log(\text{phương sai})\) vs \(\log(\text{trung bình})\) cho thấy độ dốc xấp xỉ 2, tương ứng với

\[ V(\mu) \propto \mu^2, \]

phù hợp với gamma (Fig. 11.3).

Đơn deviance của gamma là:

\[ d(y, \mu) = 2 \left(-\log \frac{y}{\mu} + \frac{y - \mu}{\mu}\right). \]

Hàm liên kết chuẩn (canonical link) là nghịch đảo

\[ \eta = \frac{1}{\mu}, \]

nhưng trong thực tế thường dùng log-link để tránh các ràng buộc và cho phép giải thích theo nhân tử.

11.3 The Inverse Gaussian Distribution

Định nghĩa và hàm mật độ:

Phân phối Gaussian nghịch đảo (Inverse Gaussian) dùng để mô hình hóa dữ liệu liên tục dương, đặc biệt khi dữ liệu lệch phải mạnh hơn Gamma.

Hàm mật độ xác suất:

\[ P(y; \mu, \phi) = \frac{1}{\sqrt{2\pi y^3 \phi}} \exp\left( -\frac{(y - \mu)^2}{2\phi y \mu^2} \right) \]

Với \(y > 0\), \(\mu > 0\), và tham số phân tán \(\phi > 0\).

Phương sai:

\[ V(\mu) = \mu^3 \]

Link hàm chuẩn tắc (canonical link):

Canonical link:

\[ \eta = \mu^{-2} \]

Tuy nhiên, trong thực tế thường dùng log-link hoặc các hàm liên kết khác để đảm bảo \(\mu > 0\) và thuận tiện cho việc diễn giải.

Deviance:

Unit deviance:

\[ d(y, \mu) = \frac{(y - \mu)^2}{y \mu^2} \]

Residual deviance:

\[ D(y, \hat{\mu}) = \sum_{i=1}^{n} w_i d(y_i, \hat{\mu}_i) \]

Đặc biệt: deviance phân phối theo \(\chi^2_{n - p}\) một cách chính xác nhờ phép xấp xỉ saddlepoint – rất hiệu quả với phân phối inverse Gaussian.

Ứng dụng thực tế:

Chuyển động Brown (Brownian motion) với xu hướng dịch chuyển dương:

• Thời gian để đạt đến một khoảng cách cố định \(\delta\): phân phối inverse Gaussian.
• Khoảng cách sau một thời gian cố định \(T\): phân phối chuẩn.

Ví dụ: Nếu một hạt bắt đầu tại vị trí 0 và chuyển động Brown với độ trôi 0.5, thì: - Thời gian để vượt qua \(\delta = 5\) có phân phối inverse Gaussian. - Vị trí sau \(T = 20\) có phân phối chuẩn (normal).

11.4 Link Functions (Hàm liên kết)

Các hàm liên kết thường dùng cho Gamma & Inverse Gaussian:

• log (logarithmic): phổ biến nhất vì đảm bảo \(\mu > 0\)
• identity (đồng nhất)
• inverse (nghịch đảo)
• \(1/\mu^2\): là canonical link cho inverse Gaussian

Ví dụ (Example 11.4): Mô hình hóa dữ liệu lime với các link khác nhau

Dữ liệu: Foliage ~ Origin * log(DBH) với phân phối Gamma.

✔ Dùng log-link:

lime.log <- glm(Foliage ~ Origin * log(DBH),
                family = Gamma(link = "log"),
                data = lime)

11.5 Ước lượng Tham số Phân tán \(\phi\)

11.5.1 Ước lượng \(\phi\) cho Phân phối Gamma

• Trong mô hình GLM với phân phối Gamma, tham số phân tán \(\phi\) không có nghiệm chính xác (closed-form).
• Hàm digamma được định nghĩa là:

\[ \psi(x) = \frac{d}{dx} \log \Gamma(x) = \frac{\Gamma'(x)}{\Gamma(x)} \]

• Ước lượng hợp lý tối đa (MLE) là nghiệm của phương trình sau (giải bằng phương pháp số):

\[ D(y, \hat{\mu}) = -2 \sum_{i=1}^n w_i \log \phi - w_i \log w_i + w_i \psi\left(\frac{w_i}{\phi}\right) \]

11.5.2 Hai phương pháp ước lượng thay thế

1. Ước lượng Pearson:

\[ \bar{\phi} = \frac{1}{n - p} \sum_{i=1}^n w_i \left( \frac{y_i - \hat{\mu}_i}{\hat{\mu}_i} \right)^2 \]

2. Ước lượng từ độ lệch trung bình (mean deviance):

\[ \tilde{\phi} = \frac{D(y, \hat{\mu})}{n - p} \]

11.5.3 Ước lượng \(\phi\) cho Phân phối Inverse Gaussian

• Đối với phân phối Inverse Gaussian, ước lượng hợp lý tối đa (MLE) của \(\phi\) có công thức chính xác:

\[ \hat{\phi} = \frac{D(y, \hat{\mu})}{n} \]

Tuy nhiên, ước lượng này có độ chệch (biased) nên trong thực hành, thường sử dụng ước lượng từ độ lệch trung bình:

Mean Deviance Estimate (gần không chệch):

\[ \tilde{\phi} = \frac{D(y, \hat{\mu})}{n - p} \]

• Ước lượng này gần đúng với modified profile likelihood, thường ít bị chệch hơn.

Ước lượng Pearson (Pearson Estimate):

\[ \bar{\phi} = \frac{1}{n - p} \sum_{i=1}^{n} w_i \left( \frac{y_i - \hat{\mu}_i}{\hat{\mu}_i^{3/2}} \right)^2 \]

Ghi chú: \(\bar{\phi}\) có thể phù hợp hơn nếu dữ liệu đầu ra \(y_i\) bị làm tròn quá mức.

CHƯƠNG 12: Tweedie GLMs

12.1 Các EDM của Tweedie

12.1.1 Giới thiệu phân phối Tweedie

Các phân phối thuộc họ EDM (Exponential Dispersion Models) có dạng hàm phương sai:

• Normal: \(V(\mu) = \mu^0 = 1\)
• Poisson: \(V(\mu) = \mu^1\)
• Gamma: \(V(\mu) = \mu^2\)
• Inverse Gaussian: \(V(\mu) = \mu^3\)

Phân phối Tweedie tổng quát hoá tất cả các dạng trên với: \[ V(\mu) = \mu^\xi, \quad \text{với } \xi \in \mathbb{R} \setminus (0, 1) \] Trong đó, \(\xi\) được gọi là chỉ số Tweedie (index parameter).

\(\xi\)	Loại dữ liệu	Ví dụ
\(\xi < 0\)	Liên tục trên toàn trục số	Không dùng trong thực tế
\(\xi = 0\)	Phân phối chuẩn	Đã học Ch. 2–3
\(\xi = 1\)	Phân phối rời rạc (Poisson)	Chương 10
\(1 < \xi < 2\)	Liên tục có giá trị 0	Mô hình lượng mưa
\(\xi = 2\)	Gamma	Chương 11
\(\xi = 3\)	Inverse Gaussian	Chương 11
\(\xi > 2\)	Liên tục, dương, lệch phải mạnh	Phép biến đổi \(1/y\)

Lưu ý: Phân phối Tweedie thường không có dạng xác suất đóng, trừ các trường hợp đặc biệt. Do đó, cần sử dụng phương pháp gần đúng số học để đánh giá hàm xác suất.

12.1.2 Cấu trúc của Tweedie EDMs

Hàm liên kết và hàm tích phân log-chuẩn hóa như sau:

• Nếu \(\xi \neq 1\), thì: \[ \theta = \frac{\mu^{1-\xi} - 1}{1 - \xi} \]
• Nếu \(\xi = 1\), thì: \[ \theta = \log(\mu) \]

Tương tự, hàm tích phân log-chuẩn hoá \(\kappa(\theta)\):

• Nếu \(\xi \neq 2\): \[ \kappa(\theta) = \frac{\mu^{2 - \xi} - 1}{2 - \xi} \]
• Nếu \(\xi = 2\): \[ \kappa(\theta) = \log(\mu) \]

Độ lệch đơn vị (unit deviance):

\[ d(y, \mu) = \begin{cases} \text{Phức tạp – cần tính theo } y, \mu, \xi & \text{nếu } \xi \neq 1, 2 \\ 2 \left[ y \log\left(\frac{y}{\mu}\right) - (y - \mu) \right] & \text{nếu } \xi = 1 \\ 2 \left[ -\log\left(\frac{y}{\mu}\right) + \frac{y - \mu}{\mu} \right] & \text{nếu } \xi = 2 \\ \end{cases} \]

Tweedie Rescaling Identity:

\[ P_\xi(y; \mu, \phi) = c \cdot P_\xi(cy; c\mu, c^{2 - \xi} \phi) \]

12.1.3 EDMs Tweedie cho dữ liệu liên tục dương

Khi phân phối Gamma hoặc Inverse Gaussian không phù hợp (dữ liệu lệch mạnh, có đuôi dài), có thể sử dụng các phân phối Tweedie với \(\xi > 2\).

Ví dụ 12.1: Dữ liệu thời gian sống sót của động vật bị đầu độc:

• Gồm 3 loại độc tố × 4 điều trị → 48 quan sát.
• Biểu đồ sơ bộ (Hình 12.3): khác biệt rõ giữa nhóm điều trị.
• Đường hồi quy log(mean) vs. log(variance) có slope ≈ 3.95 → gợi ý \(\xi \approx 4\).

Bạn muốn mình tiếp tục viết phần mã R mô phỏng hoặc phân tích theo chương này không?

12.2 Tweedie GLMs

12.2.1 Giới thiệu

Mô hình tuyến tính tổng quát Tweedie có dạng: \[ \texttt{glm(Tweedie, } \xi \texttt{; link function)} \]

• Phạm vi chỉ số Tweedie \(\xi\) thường gặp:

  • \(1 < \xi < 2\): phân phối liên tục có giá trị bằng 0 (zero-inflated)
  • \(\xi > 2\): phân phối liên tục dương, không có giá trị 0

• Vì \(\mu > 0\), hàm liên kết phổ biến là: \[ \eta = \log(\mu) \]

• Tham số phân tán \(\phi\) thường được ước lượng bằng Pearson estimator.

• Giá trị \(\xi\) phải được xác định trước khi ước lượng mô hình GLM. Tuy nhiên:

  • Ước lượng \(\xi\) ít ảnh hưởng đến \(\hat{\beta}\)
  • Vì \(\hat{\xi}\) và \(\hat{\beta}\) thường có tương quan thấp

• Hồi quy Box–Cox có thể xấp xỉ Tweedie GLM, nhưng không hiệu quả khi có giá trị 0 hoặc gần 0 → Tweedie GLM phù hợp hơn.

12.2.2 Ước lượng chỉ số \(\xi\)

Phương pháp tuyến tính sơ bộ: Dựa trên mối quan hệ: \[ \log(\text{Var}[y]) = \log(\phi) + \xi \log(\mu) \]

• Chia dữ liệu thành các nhóm → tính \(\log(\text{variance})\) và \(\log(\text{mean})\)
• Hồi quy tuyến tính để ước lượng \(\xi\)

Ví dụ: Dữ liệu mưa tại Quilpie: - Phân nhóm theo giai đoạn thập kỷ và chỉ số SOI phase - Ước lượng cho thấy: \(\hat{\xi} \approx 1.55\) và \(\hat{\xi} \approx 1.95\)

Phương pháp nghiêm ngặt hơn: Ước lượng hợp lý tối đa (MLE)

Dùng profile likelihood:

1. Tạo lưới các giá trị \(\xi\)
2. Với mỗi \(\xi\):
   • Ước lượng \(\beta\) và \(\phi\) qua GLM
   • Tính log-likelihood
3. Vẽ profile likelihood theo \(\xi\)
4. \(\xi\) ứng với log-likelihood lớn nhất là MLE

Dùng trong R với package tweedie:

library(tweedie)

profile <- tweedie.profile(y ~ x1 + x2,
                           data = mydata,
                           p.vec = seq(1.1, 2, by = 0.05),
                           do.plot = TRUE)

profile$xi.max   # MLE của ξ
profile$phi.max  # MLE của φ
profile$ci       # Khoảng tin cậy 95% cho ξ

PHẦN 2

THỐNG KÊ MÔ TẢ DATA Supermarket Transactions.

2.1 Đọc Data

data <- read.csv("C:/Users/DELL/Downloads/Supermarket Transactions.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)

2.2 Khái quát data

2.2.1 Nội dung data

• Dữ liệu Supermarket Transactions ghi lại các giao dịch mua hàng tại chuỗi siêu thị, bao gồm thông tin về khách hàng, vị trí địa lý, cùng chi tiết các mặt hàng đã được mua.

• Bộ dữ liệu này giúp phân tích thói quen tiêu dùng của khách hàng, phân loại thị trường, cũng như đánh giá hiệu quả kinh doanh theo từng sản phẩm, khu vực và nhóm đối tượng khách hàng.

2.2.2 Danh sách các biến

• Tên các biến trong data bao gồm 16 biến định lượng và định tính

colnames(data)

##  [1] "X"                 "PurchaseDate"      "CustomerID"       
##  [4] "Gender"            "MaritalStatus"     "Homeowner"        
##  [7] "Children"          "AnnualIncome"      "City"             
## [10] "StateorProvince"   "Country"           "ProductFamily"    
## [13] "ProductDepartment" "ProductCategory"   "UnitsSold"        
## [16] "Revenue"

2.3 Thống kê mô tả cho các biến

2.3.1 Cấu trúc tổng quát

str(data)

## 'data.frame':    14059 obs. of  16 variables:
##  $ X                : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ PurchaseDate     : chr  "2007-12-18" "2007-12-20" "2007-12-21" "2007-12-21" ...
##  $ CustomerID       : int  7223 7841 8374 9619 1900 6696 9673 354 1293 7938 ...
##  $ Gender           : chr  "F" "M" "F" "M" ...
##  $ MaritalStatus    : chr  "S" "M" "M" "M" ...
##  $ Homeowner        : chr  "Y" "Y" "N" "Y" ...
##  $ Children         : int  2 5 2 3 3 3 2 2 3 1 ...
##  $ AnnualIncome     : chr  "$30K - $50K" "$70K - $90K" "$50K - $70K" "$30K - $50K" ...
##  $ City             : chr  "Los Angeles" "Los Angeles" "Bremerton" "Portland" ...
##  $ StateorProvince  : chr  "CA" "CA" "WA" "OR" ...
##  $ Country          : chr  "USA" "USA" "USA" "USA" ...
##  $ ProductFamily    : chr  "Food" "Food" "Food" "Food" ...
##  $ ProductDepartment: chr  "Snack Foods" "Produce" "Snack Foods" "Snacks" ...
##  $ ProductCategory  : chr  "Snack Foods" "Vegetables" "Snack Foods" "Candy" ...
##  $ UnitsSold        : int  5 5 3 4 4 3 4 6 1 2 ...
##  $ Revenue          : num  27.38 14.9 5.52 4.44 14 ...

Các biến trong data bao gồm:

Ký hiệu biến	Mô tả biến
X	Số thứ tự dòng, có thể bỏ qua
PurchaseDate	Ngày khách hàng mua hàng
CustomerID	Mã định danh khách hàng
Gender	Giới tính: “F” = nữ, “M” = nam
MaritalStatus	Tình trạng hôn nhân: “S” = độc thân, “M” = đã kết hôn
Homeowner	Có sở hữu nhà: “Y” = có, “N” = không
Children	Số con trong gia đình
AnnualIncome	Nhóm thu nhập hàng năm (ví dụ: “$30K - $50K”)
City	Thành phố nơi sinh sống hoặc mua hàng
StateorProvince	Bang hoặc tỉnh nơi sinh sống hoặc mua hàng
Country	Quốc gia nơi sinh sống hoặc mua hàng
ProductFamily	Nhóm sản phẩm chính (ví dụ: Food, Drink, …)
ProductDepartment	Phòng ban sản phẩm (ví dụ: Snack Foods, Produce, …)
ProductCategory	Danh mục sản phẩm cụ thể (ví dụ: Candy, Vegetables, …)
UnitsSold	Số lượng sản phẩm bán ra
Revenue	Doanh thu từ giao dịch (đơn vị USD)

2.3.2 Thống kê mô tả biến định lượng

library(psych)
describe(data[, c("UnitsSold", "Revenue", "Children")])

##           vars     n  mean   sd median trimmed  mad  min  max range  skew
## UnitsSold    1 14059  4.08 1.17   4.00    4.08 1.48 1.00  8.0  7.00  0.01
## Revenue      2 14059 13.00 8.22  11.25   12.05 7.40 0.53 56.7 56.17  1.13
## Children     3 14059  2.53 1.49   3.00    2.53 1.48 0.00  5.0  5.00 -0.02
##           kurtosis   se
## UnitsSold    -0.44 0.01
## Revenue       1.39 0.07
## Children     -1.03 0.01

2.3.2.1 Biến UnitsSold

Thống kê	Giá trị
N	14,059
Trung bình (Mean)	4.08
Độ lệch chuẩn (SD)	1.17
Trung vị (Median)	4.00
Trung bình cắt tỉa (Trimmed Mean)	4.08
Độ lệch tuyệt đối trung bình (MAD)	1.48
Giá trị nhỏ nhất (Min)	1.00
Giá trị lớn nhất (Max)	8.00
Khoảng giá trị (Range)	7.00
Độ lệch (Skewness)	0.012
Độ nhọn (Kurtosis)	-0.44
Sai số chuẩn (SE)	0.0099

Nhận xét:

• Trung bình là khoảng 4 đơn vị, cho thấy trung bình mỗi giao dịch khách hàng mua 4 sản phẩm.

• Độ lệch chuẩn 1.17 khá nhỏ so với trung bình, chứng tỏ số lượng bán ra ít biến động, dữ liệu tập trung tương đối gần nhau.

• Trung vị bằng 4, gần sát trung bình, nghĩa là phân phối dữ liệu khá cân đối, không lệch nhiều.

• Skewness = 0.012, gần 0, xác nhận phân phối dữ liệu gần như đối xứng.

• Kurtosis âm (-0.44) cho thấy phân phối hơi bằng phẳng (dữ liệu không có đỉnh nhọn hay ngoại lệ lớn).

• Khoảng giá trị từ 1 đến 8 đơn vị thể hiện số lượng bán ra trong giao dịch nằm trong khoảng hợp lý, không có giá trị bất thường.

Kết luận: biến UnitsSold có phân phối cân đối, ít biến động, phù hợp cho các phân tích thống kê chuẩn và không bị ảnh hưởng nhiều bởi ngoại lệ.

2.3.2.2 Biến Revenue

Thống kê	Giá trị
N	14,059
Trung bình (Mean)	13.00
Độ lệch chuẩn (SD)	8.22
Trung vị (Median)	11.25
Trung bình cắt tỉa (Trimmed Mean)	12.05
Độ lệch tuyệt đối trung bình (MAD)	7.40
Giá trị nhỏ nhất (Min)	0.53
Giá trị lớn nhất (Max)	56.70
Khoảng giá trị (Range)	56.17
Độ lệch (Skewness)	1.13
Độ nhọn (Kurtosis)	1.39
Sai số chuẩn (SE)	0.0693

Nhận xét:

• Trung bình khoảng 13 USD, với độ lệch chuẩn lớn (8.22), cho thấy doanh thu có sự phân hóa mạnh.

• Trung vị 11.25 USD thấp hơn trung bình, chứng tỏ phân phối lệch phải.

• Skewness = 1.13, khá cao, thể hiện doanh thu bị lệch phải mạnh, nghĩa là có một số giao dịch với doanh thu rất cao kéo trung bình lên.

• Kurtosis dương (1.39) cho thấy phân phối có đỉnh nhọn và nhiều ngoại lệ (outliers) ở phía đuôi phải.

• Giá trị doanh thu nhỏ nhất là 0.53 USD, lớn nhất lên đến 56.7 USD, khoảng cách rất rộng, thể hiện sự đa dạng về quy mô giao dịch.

Kết luận: biến Revenue phân phối lệch phải rõ ràng với nhiều giao dịch doanh thu nhỏ và một số ít giao dịch doanh thu rất lớn. Khi phân tích biến này cần chú ý đến ảnh hưởng của các giá trị ngoại lệ, có thể cân nhắc biến đổi log hoặc sử dụng phương pháp phân tích không phụ thuộc phân phối chuẩn.

2.3.2.3 Biến Children

Thống kê	Giá trị
N	14,059
Trung bình (Mean)	2.53
Độ lệch chuẩn (SD)	1.49
Trung vị (Median)	3.00
Trung bình cắt tỉa (Trimmed Mean)	2.53
Độ lệch tuyệt đối trung bình (MAD)	1.48
Giá trị nhỏ nhất (Min)	0.00
Giá trị lớn nhất (Max)	5.00
Khoảng giá trị (Range)	5.00
Độ lệch (Skewness)	-0.02
Độ nhọn (Kurtosis)	-1.03
Sai số chuẩn (SE)	0.0126

Nhận xét:

• rung bình là 2.53 con, khá phổ biến trong gia đình Việt Nam hiện đại.

• Độ lệch chuẩn 1.49 cho thấy sự phân tán tương đối lớn, nhưng không quá rộng.

• Trung vị là 3, tức một nửa gia đình có từ 3 con trở lên.

• Skewness âm nhẹ (-0.02), phân phối gần như đối xứng, hơi nghiêng trái nhẹ.

• Kurtosis âm (-1.03) cho thấy phân phối hơi phẳng, không có đỉnh nhọn, tức số lượng con phân bố đều, không tập trung cao vào một số ít giá trị.

• Số con tối thiểu là 0 (có thể có gia đình không có con), tối đa 5 con.

Kết luận: biến Children phân phối khá đều, không có sự lệch lớn hay ngoại lệ đáng kể. Phù hợp sử dụng trực tiếp trong các phân tích mà không cần xử lý thêm.

2.3.4 Thống kê mô tả biến định tính

Biến định tính thường chứa nhiều giá trị lặp lại (như giới tính “M”/“F”, tình trạng hôn nhân “S”/“M”…).

Bảng tần số cho biết mỗi mức độ xuất hiện bao nhiêu lần, giúp hiểu ngay cấu trúc phân phối của biến mà không cần nhìn qua hàng ngàn dòng dữ liệu.

2.3.4.1 Tạo bảng tần số

2.3.4.1.1 Biến Gender

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

## 
## Attaching package: 'ggplot2'

## The following objects are masked from 'package:psych':
## 
##     %+%, alpha

# Tóm tắt dữ liệu theo giới tính
summary_gender <- data %>%
  group_by(Gender) %>%
  summarise(
    Frequency = n(),
    Percentage = n() / nrow(data) * 100
  )

# In bảng tóm tắt
print(summary_gender)

## # A tibble: 2 × 3
##   Gender Frequency Percentage
##   <chr>      <int>      <dbl>
## 1 F           7170       51.0
## 2 M           6889       49.0

# Vẽ biểu đồ cột thể hiện tỷ lệ phần trăm theo giới tính
ggplot(summary_gender, aes(x = Gender, y = Percentage, fill = Gender)) +
  geom_col() +
   geom_text(aes(label = sprintf("%.1f%%", Percentage)), 
            hjust = -0.1, size = 3) +
  labs(title = "Tỷ lệ phần trăm theo giới tính",
       x = "Giới tính",
       y = "Tỷ lệ (%)") +
  theme_minimal()

Nhận xét:

• Mẫu dữ liệu có sự phân bố giới tính tương đối đồng đều, với nữ chiếm 50.99% và nam chiếm 49.00%.

• Sự chênh lệch chỉ khoảng 2%, cho thấy không có sự thiên lệch giới tính đáng kể trong dữ liệu.

2.3.4.1.2 Biến MaritalStatus

library(dplyr)


summary_marital <- data %>%
  filter(!is.na(MaritalStatus)) %>%  # Bỏ NA nếu có
  group_by(MaritalStatus) %>%
  summarise(
    Frequency = n(),
    Percentage = n() / nrow(data %>% filter(!is.na(MaritalStatus))) * 100
  )

print(summary_marital)

## # A tibble: 2 × 3
##   MaritalStatus Frequency Percentage
##   <chr>             <int>      <dbl>
## 1 M                  6866       48.8
## 2 S                  7193       51.2

library(ggplot2)

ggplot(summary_marital, aes(x = MaritalStatus, y = Percentage, fill = MaritalStatus)) +
  geom_col() +
   geom_text(aes(label = sprintf("%.1f%%", Percentage)), 
            hjust = -0.1, size = 3) +
  labs(title = "Tỷ lệ phần trăm theo tình trạng hôn nhân",
       x = "Tình trạng hôn nhân",
       y = "Tỷ lệ (%)") +
  theme_minimal()

Nhận xét:

• Tỷ lệ khách hàng độc thân chiếm 51.15%, cao hơn nhẹ so với tỷ lệ khách hàng đã kết hôn (48.85%).

• Mức chênh lệch là 2.3%, cho thấy hai nhóm này khá cân bằng về số lượng.

2.3.4.1.3 Biến Homeowner

library(ggplot2)

summary_homeowner <- data %>%
  filter(!is.na(Homeowner)) %>%
  group_by(Homeowner) %>%
  summarise(
    Frequency = n(),
    Percentage = n() / nrow(data %>% filter(!is.na(Homeowner))) * 100
  )

ggplot(summary_homeowner, aes(x = Homeowner, y = Percentage, fill = Homeowner)) +
  geom_col() +
   geom_text(aes(label = sprintf("%.1f%%", Percentage)), 
            hjust = -0.1, size = 3) +
  labs(title = "Tỷ lệ phần trăm theo tình trạng sở hữu nhà",
       x = "Tình trạng sở hữu nhà",
       y = "Tỷ lệ (%)") +
  theme_minimal()

Nhận xét:

• 60.06% khách hàng trong mẫu dữ liệu có sở hữu nhà, trong khi 39.94% không sở hữu nhà.

• Điều này cho thấy đa số khách hàng trong tập dữ liệu đang có điều kiện kinh tế tương đối ổn định (vì sở hữu nhà thường gắn với thu nhập và tài chính ổn định hơn).

2.3.4.1.4 Biến AnnualIncome

library(dplyr)
library(ggplot2)

# Tóm tắt dữ liệu theo AnnualIncome (loại bỏ NA nếu có)
summary_income <- data %>%
  filter(!is.na(AnnualIncome)) %>%
  group_by(AnnualIncome) %>%
  summarise(
    Frequency = n(),
    Percentage = n() / 14059 * 100
  )

# In bảng tóm tắt
print(summary_income)

## # A tibble: 8 × 3
##   AnnualIncome  Frequency Percentage
##   <chr>             <int>      <dbl>
## 1 $10K - $30K        3090      22.0 
## 2 $110K - $130K       643       4.57
## 3 $130K - $150K       760       5.41
## 4 $150K +             273       1.94
## 5 $30K - $50K        4601      32.7 
## 6 $50K - $70K        2370      16.9 
## 7 $70K - $90K        1709      12.2 
## 8 $90K - $110K        613       4.36

# Vẽ biểu đồ tròn
ggplot(summary_income, aes(x = "", y = Percentage, fill = AnnualIncome)) +
  geom_col(width = 1) +
   geom_text(aes(label = sprintf("%.1f%%", Percentage)), 
            hjust = -0.1, size = 3) +
  coord_polar(theta = "y") +
  labs(title = "Tỷ lệ phần trăm theo mức thu nhập hàng năm") +
  theme_void() +  # Xóa trục và lưới để giống biểu đồ tròn
  theme(legend.title = element_blank())

Nhận xét:

• Phân bố thu nhập trong data không đồng đều, phần lớn khách hàng thuộc các nhóm thu nhập trung bình.

• Nhóm chiếm tỷ lệ lớn nhất là $30K – $50K, với 32.73%, theo sau là nhóm $10K – $30K (21.98%) và $50K – $70K (16.86%).

• Các nhóm thu nhập cao hơn như $90K trở lên chiếm tỷ lệ rất thấp (chỉ khoảng 16% tổng cộng), cho thấy phần lớn khách hàng trong mẫu có thu nhập dưới $90K/năm.

• Phân bố thu nhập này phản ánh rõ tầng lớp thu nhập phổ thông, có thể hữu ích để xác định phân khúc thị trường mục tiêu cho các sản phẩm bình dân hoặc tầm trung.

2.3.4.1.5 Biến City

library(dplyr)
library(ggplot2)

# Loại bỏ các dòng NA trong City
filtered_data <- data %>% filter(!is.na(City))

# Tính tổng số dòng sau khi loại NA
total_rows <- nrow(filtered_data)

# Tóm tắt số lượng và phần trăm
summary_city <- filtered_data %>%
  group_by(City) %>%
  summarise(
    Frequency = n(),
    Percentage = n() / 14059 * 100
  )

# Hiển thị bảng
print(summary_city)

## # A tibble: 23 × 3
##    City          Frequency Percentage
##    <chr>             <int>      <dbl>
##  1 Acapulco            383      2.72 
##  2 Bellingham          143      1.02 
##  3 Beverly Hills       811      5.77 
##  4 Bremerton           834      5.93 
##  5 Camacho             452      3.22 
##  6 Guadalajara          75      0.533
##  7 Hidalgo             845      6.01 
##  8 Los Angeles         926      6.59 
##  9 Merida              654      4.65 
## 10 Mexico City         194      1.38 
## # ℹ 13 more rows

# Vẽ biểu đồ cột
ggplot(summary_city, aes(x = reorder(City, -Frequency), y = Percentage, fill = City)) +
  geom_col() +
   geom_text(aes(label = sprintf("%.1f%%", Percentage)), 
            hjust = -0.1, size = 3) +
  labs(title = "Tỷ lệ phần trăm theo từng thành phố",
       x = "Thành phố",
       y = "Tỷ lệ %") +
  theme_minimal() +
  theme(legend.position = "none") +
  coord_flip()

Nhận xét:

• Thành phố có số khách hàng nhiều nhất là Los Angeles, chiếm khoảng 6.59%, theo sau là Hidalgo (6.01%) và Bremerton (5.93%).

• Một số thành phố có tỷ lệ khách hàng rất nhỏ như Guadalajara (0.53%) hay Bellingham (1.02%), cho thấy phạm vi phân phối khách hàng tập trung tại một vài thành phố lớn.

• Việc các thành phố như Los Angeles, Beverly Hills và Mexico City xuất hiện có thể phản ánh các thị trường tiêu dùng năng động hoặc được nhắm mục tiêu cụ thể trong chiến dịch bán hàng.

2.3.4.1.6 Biến StateorProvince

library(dplyr)
library(ggplot2)

# Loại bỏ các dòng bị thiếu giá trị (NA) ở StateorProvince
filtered_data <- data %>% filter(!is.na(StateorProvince))

# Tính tổng số dòng sau khi lọc NA
total_rows <- nrow(filtered_data)

# Tóm tắt theo StateorProvince
summary_state <- filtered_data %>%
  group_by(StateorProvince) %>%
  summarise(
    Frequency = n(),
    Percentage = n() / 14059 * 100
  )

# Hiển thị bảng kết quả
print(summary_state)

## # A tibble: 10 × 3
##    StateorProvince Frequency Percentage
##    <chr>               <int>      <dbl>
##  1 BC                    809      5.75 
##  2 CA                   2733     19.4  
##  3 DF                    815      5.80 
##  4 Guerrero              383      2.72 
##  5 Jalisco                75      0.533
##  6 OR                   2262     16.1  
##  7 Veracruz              464      3.30 
##  8 WA                   4567     32.5  
##  9 Yucatan               654      4.65 
## 10 Zacatecas            1297      9.23

# Vẽ biểu đồ cột
ggplot(summary_state, aes(x = reorder(StateorProvince, -Percentage), y = Percentage, fill = StateorProvince)) +
  geom_col() +
   geom_text(aes(label = sprintf("%.1f%%", Percentage)), 
            hjust = -0.1, size = 3) +
  labs(title = "Số lượng theo Bang/Tỉnh",
       x = "Bang/Tỉnh",
       y = "Số lượng") +
  theme_minimal() +
  theme(legend.position = "none") +
  coord_flip()

Nhận xét:

• Bang WA (Washington) có tỷ lệ khách hàng cao nhất (32.48%), cho thấy đây là thị trường lớn nhất trong bộ dữ liệu, theo sau là CA (California) với 19.44% và OR (Oregon) với 16.09%.

• Một số bang/tỉnh có số lượng khách hàng thấp như Jalisco (0.53%) và Guerrero (2.72%), có thể là các khu vực phụ hoặc ít được tập trung.

• Sự tập trung khách hàng ở các bang như WA, CA, OR phản ánh rõ chiến lược thị trường có thể đang hướng mạnh về khu vực Tây Bắc và Bờ Tây Hoa Kỳ.

• Các bang/tỉnh Mexico như DF, Veracruz, Zacatecas cũng có mặt, tuy nhiên tỷ lệ thấp hơn, cho thấy có thể doanh nghiệp có hoạt động mở rộng nhưng chưa thực sự mạnh tại thị trường Mexico.

2.3.4.1.7 Biến Country

library(dplyr)
library(ggplot2)

# Bước 1: Lọc bỏ NA trong Country
filtered_data <- data %>% filter(!is.na(Country))

# Bước 2: Tính tổng số dòng hợp lệ
total_rows <- nrow(filtered_data)

# Bước 3: Tóm tắt số lượng và phần trăm
summary_country <- filtered_data %>%
  group_by(Country) %>%
  summarise(
    Frequency = n(),
    Percentage = n() / 14059 * 100
  )

# In bảng tóm tắt
print(summary_country)

## # A tibble: 3 × 3
##   Country Frequency Percentage
##   <chr>       <int>      <dbl>
## 1 Canada        809       5.75
## 2 Mexico       3688      26.2 
## 3 USA          9562      68.0

# Bước 4: Vẽ biểu đồ cột
ggplot(summary_country, aes(x = reorder(Country, -Percentage), y = Percentage, fill = Country)) +
  geom_col() +
   geom_text(aes(label = sprintf("%.1f%%", Percentage)), 
            hjust = -0.1, size = 3) +
  labs(title = "Tỷ lệ phần trăm theo quốc gia",
       x = "Quốc gia",
       y = "Tỷ lệ") +
  theme_minimal() +
  theme(legend.position = "none") +
  coord_flip()

Nhận xét:

• Dữ liệu cho thấy trong ba quốc gia, Mỹ chiếm tỷ lệ lớn nhất với 68% (9,562 đơn vị), Mexico đứng thứ hai với 26% (3,688 đơn vị), còn Canada chiếm tỷ lệ nhỏ nhất khoảng 5.75% (809 đơn vị). Điều này cho thấy phân bố dữ liệu không đồng đều, Mỹ là quốc gia chiếm ưu thế rõ rệt trong mẫu.

2.3.4.1.8 Biến ProductFamily

library(dplyr)
library(ggplot2)

# Bước 1: Loại bỏ giá trị NA trong ProductFamily
filtered_data <- data %>% filter(!is.na(ProductFamily))

# Bước 2: Tính tổng số dòng hợp lệ
total_rows <- nrow(filtered_data)

# Bước 3: Tóm tắt số lượng và phần trăm theo ProductFamily
summary_product <- filtered_data %>%
  group_by(ProductFamily) %>%
  summarise(
    Frequency = n(),
    Percentage = n() / total_rows * 100
  )

# Bước 4: Vẽ biểu đồ cột theo Percentage
ggplot(summary_product, aes(x = reorder(ProductFamily, -Percentage), y = Percentage, fill = ProductFamily)) +
  geom_col() +
   geom_text(aes(label = sprintf("%.1f%%", Percentage)), 
            hjust = -0.1, size = 3) +
  labs(title = "Tỷ lệ phần trăm theo dòng sản phẩm",
       x = "Dòng sản phẩm",
       y = "Tỷ lệ (%)") +
  theme_minimal() +
  theme(legend.position = "none") +
  coord_flip()

Nhận xét:

• Nhóm Food (thực phẩm) chiếm hơn 72% tổng số giao dịch, là mảng sản phẩm chủ lực và quan trọng nhất. → Cần ưu tiên các chính sách về danh mục sản phẩm, chất lượng, giá cả và chương trình khuyến mãi cho nhóm này.

• Non-Consumable (phi tiêu dùng, như đồ gia dụng, thiết bị điện tử…) chiếm gần 19%, là mảng phụ nhưng có tiềm năng phát triển. → Có thể tận dụng cơ hội upsell hoặc cross-sell với khách hàng đã mua nhóm Food.

• Drink (thức uống) chỉ chiếm khoảng 9%, là nhóm nhỏ hơn, nên tập trung vào những sản phẩm có biên lợi nhuận cao hoặc kết hợp bán kèm cùng nhóm Food để tăng giá trị.

2.3.4.1.9 Biến ProductDepartment

library(dplyr)
library(ggplot2)

# Bước 1: Lọc bỏ dòng bị thiếu ProductDepartment
filtered_data <- data %>% filter(!is.na(ProductDepartment))

# Bước 2: Tính tổng số dòng sau khi lọc NA
total_rows <- nrow(filtered_data)

# Bước 3: Tóm tắt số lượng và phần trăm
summary_department <- filtered_data %>%
  group_by(ProductDepartment) %>%
  summarise(
    Frequency = n(),
    Percentage = n() / total_rows * 100
  )

# Bước 4: Vẽ biểu đồ cột theo phần trăm
ggplot(summary_department, aes(x = reorder(ProductDepartment, -Percentage), y = Percentage, fill = ProductDepartment)) +
  geom_col() +
   geom_text(aes(label = sprintf("%.1f%%", Percentage)), 
            hjust = -0.1, size = 3) +
  labs(title = "Tỷ lệ phần trăm theo bộ phận sản phẩm",
       x = "Bộ phận sản phẩm",
       y = "Tỷ lệ (%)") +
  theme_minimal() +
  theme(legend.position = "none") +
  coord_flip()

Nhận xét:

• Nhóm Baked Goods (bánh nướng) chiếm khoảng 3%, còn Baking Goods (nguyên liệu làm bánh) chiếm gần 7.6%, cho thấy cả thành phẩm và nguyên liệu làm bánh đều được quan tâm. → Có thể phát triển các combo hoặc chương trình khuyến mãi kết hợp giữa hai nhóm này.

• Dairy (sữa và các sản phẩm từ sữa) chiếm khoảng 6.4%, và Canned Foods (thực phẩm đóng hộp) chiếm gần 7%, đều là những nhóm sản phẩm thiết yếu trong nhóm Food. → Nên duy trì chất lượng và đa dạng sản phẩm để giữ chân khách hàng.

• Các nhóm như Beverages (đồ uống) 4.8%, Alcoholic Beverages (đồ uống có cồn) 2.5%, Breakfast Foods (thực phẩm ăn sáng) 1.3% chiếm tỷ trọng nhỏ hơn, nhưng có thể là mảng phụ bổ trợ cho nhóm chính. → Có thể cân nhắc phát triển thêm các sản phẩm đặc thù hoặc các combo đa dạng.

• Các nhóm khác như Canned Products, Carousel, Checkout chiếm tỷ lệ rất nhỏ, có thể chưa phải là trọng tâm hiện tại.

2.3.4.1.10 Biến ProductCategory

library(dplyr)
library(ggplot2)

# Loại bỏ NA trong ProductCategory
filtered_data <- data %>% filter(!is.na(ProductCategory))

# Tổng số dòng sau khi lọc NA
total_rows <- nrow(filtered_data)

# Tóm tắt Frequency và Percentage
summary_category <- filtered_data %>%
  group_by(ProductCategory) %>%
  summarise(
    Frequency = n(),
    Percentage = n() / total_rows * 100
  )

print(summary_category)

## # A tibble: 45 × 3
##    ProductCategory   Frequency Percentage
##    <chr>                 <int>      <dbl>
##  1 Baking Goods            484      3.44 
##  2 Bathroom Products       365      2.60 
##  3 Beer and Wine           356      2.53 
##  4 Bread                   425      3.02 
##  5 Breakfast Foods         417      2.97 
##  6 Candles                  45      0.320
##  7 Candy                   352      2.50 
##  8 Canned Anchovies         44      0.313
##  9 Canned Clams             53      0.377
## 10 Canned Oysters           35      0.249
## # ℹ 35 more rows

# Vẽ biểu đồ cột theo phần trăm với nhãn phần trăm
ggplot(summary_category, aes(x = reorder(ProductCategory, -Percentage), y = Percentage, fill = ProductCategory)) +
  geom_col() +
  geom_text(aes(label = sprintf("%.1f%%", Percentage)), 
            hjust = -0.1, size = 3) +  # Hiển thị nhãn phần trăm bên ngoài cột
  labs(title = "Tỷ lệ phần trăm theo danh mục sản phẩm",
       x = "Danh mục sản phẩm",
       y = "Tỷ lệ (%)") +
  theme_minimal() +
  theme(legend.position = "none") +
  coord_flip() +
  ylim(0, max(summary_category$Percentage) * 1.15)  # Cách để có chỗ cho nhãn

Nhận xét:

• Nhóm Baking Goods, Bread và Breakfast Foods đều chiếm khoảng 3% mỗi nhóm, phản ánh nhu cầu ổn định cho thực phẩm cơ bản và nguyên liệu chế biến.

• Bathroom Products, Beer and Wine, Candy cũng chiếm tỷ lệ tương tự (~2.5%), cho thấy sự đa dạng trong sản phẩm tiêu dùng.

• Các sản phẩm đóng hộp như Canned Anchovies, Clams, Oysters chiếm tỷ lệ rất nhỏ dưới 0.4%, phù hợp với thị trường ngách hoặc sản phẩm đặc biệt.

2.3.4.2 Kết quả thống kê mô tả các biến định tính

# Load thư viện cần thiết
library(dplyr)
library(knitr)

# Danh sách các biến định tính
cat_vars <- c("Gender", "MaritalStatus", "Homeowner", "AnnualIncome",
              "City", "StateorProvince", "Country",
              "ProductFamily", "ProductDepartment", "ProductCategory")

# Hàm tóm tắt từng biến định tính
cat_summary <- function(df, var) {
  v <- trimws(df[[var]])         # Xử lý khoảng trắng dư
  v <- as.factor(v)              # Ép kiểu factor
  tbl  <- table(v, useNA = "ifany")
  prop <- prop.table(tbl)
  
  data.frame(
    Bien     = var,
    Muc      = names(tbl),
    So_luong = as.vector(tbl),
    Ty_le    = round(100 * as.vector(prop), 2),
    row.names = NULL,
    stringsAsFactors = FALSE
  )
}

# Gộp tất cả các biến lại thành 1 bảng
big_dt <- bind_rows(lapply(cat_vars, function(x) cat_summary(data, x)))

# Sắp xếp
big_dt <- big_dt %>%
  arrange(Bien, desc(So_luong))

# Tính thống kê mô tả
descriptive_stats_full <- big_dt %>%
  group_by(Bien) %>%
  summarise(
    Mean    = mean(So_luong, na.rm = TRUE),
    StdDev  = sd(So_luong, na.rm = TRUE),
    Min     = min(So_luong, na.rm = TRUE),
    Q1      = quantile(So_luong, 0.25, na.rm = TRUE),
    Median  = quantile(So_luong, 0.5, na.rm = TRUE),
    Q3      = quantile(So_luong, 0.75, na.rm = TRUE),
    Max     = max(So_luong, na.rm = TRUE)
  )

# Hiển thị bảng
kable(descriptive_stats_full,
      digits = 2,
      col.names = c("Biến", "Trung bình", "Độ lệch chuẩn", "Min", "Q1", "Trung vị", "Q3", "Max"),
      caption = "Thống kê mô tả đầy đủ số lượng các mức theo từng biến định tính")

Thống kê mô tả đầy đủ số lượng các mức theo từng biến định tính

Biến	Trung bình	Độ lệch chuẩn	Min	Q1	Trung vị	Q3	Max
AnnualIncome	1757.38	1511.35	273	635.50	1234.5	2550.00	4601
City	611.26	370.75	75	285.00	633.0	870.50	1386
Country	4686.33	4461.08	809	2248.50	3688.0	6625.00	9562
Gender	7029.50	198.70	6889	6959.25	7029.5	7099.75	7170
Homeowner	7029.50	2000.41	5615	6322.25	7029.5	7736.75	8444
MaritalStatus	7029.50	231.22	6866	6947.75	7029.5	7111.25	7193
ProductCategory	312.42	358.15	35	102.00	197.0	356.00	1728
ProductDepartment	639.05	569.43	59	190.50	390.5	958.50	1994
ProductFamily	4686.33	4786.18	1250	1953.00	2656.0	6404.50	10153
StateorProvince	1405.90	1393.40	75	511.50	812.0	2020.75	4567

• AnnualIncome: Trung bình thu nhập hàng năm của mẫu là 1,757.38 với độ lệch chuẩn cao 1,511.35, cho thấy sự phân bố thu nhập rất đa dạng trong dữ liệu. Mức thấp nhất là 273 và cao nhất đạt tới 4,601, thể hiện khoảng cách lớn giữa các cá nhân.

• City: Số lượng quan sát theo từng thành phố trung bình là 611.26, dao động từ 75 đến 1,386 với độ lệch chuẩn 370.75, cho thấy có sự chênh lệch đáng kể về kích thước mẫu ở các thành phố khác nhau, có thể do quy mô dân số hoặc mức độ tham gia khảo sát khác nhau.

• Country: Trung bình số lượng theo quốc gia khá lớn, 4,686.33 với độ lệch chuẩn rất cao 4,461.08, phạm vi từ 809 đến 9,562. Điều này cho thấy sự chênh lệch rõ ràng trong quy mô mẫu giữa các quốc gia, có thể phản ánh quy mô dân số hoặc mức độ thu thập dữ liệu.

• Gender, Homeowner, MaritalStatus: Ba biến này có trung bình số lượng xấp xỉ 7,000 và độ lệch chuẩn thấp hơn so với các biến khác, từ 198.7 đến 2,000. Điều này chỉ ra sự phân bố đều hơn về các nhóm giới tính, tình trạng sở hữu nhà và tình trạng hôn nhân trong mẫu dữ liệu.

• ProductCategory: Trung bình số lượng sản phẩm theo danh mục là 312.42, độ lệch chuẩn 358.15, dao động từ 35 đến 1,728. Phân bố này cho thấy một số danh mục sản phẩm rất phổ biến, trong khi những danh mục khác có mức độ quan tâm thấp hơn.

• ProductDepartment: Có trung bình 639.05 quan sát theo phòng ban sản phẩm, với phạm vi từ 59 đến 1,994 và độ lệch chuẩn 569.43, cho thấy sự đa dạng và chênh lệch lớn về số lượng trong từng phòng ban sản phẩm.

• ProductFamily: Số liệu về nhóm sản phẩm cho thấy sự phân bố rất rộng với trung bình 4,686.33 và độ lệch chuẩn lớn 4,786.18, từ 1,250 đến 10,153. Điều này cho thấy có những nhóm sản phẩm rất phổ biến, chiếm phần lớn số liệu, trong khi nhóm khác ít được đại diện hơn.

• StateorProvince: Số lượng theo khu vực hành chính có trung bình 1,405.90 và độ lệch chuẩn 1,393.40, dao động từ 75 đến 4,567, thể hiện sự phân bố không đồng đều về quy mô mẫu theo khu vực.