Một số mô hình hồi quy và áp dụng đối với tập dữ liệu dự báo suy thận mạn
Sinh viên thực hiện: Vũ Đức Trung
2026-01-11
1 ĐẶT VẤN ĐỀ
Trong thời đại mới, sự giao thoa giữa y học lâm sàng và khoa học dữ liệu đã mở ra những hướng đi đột phá trong việc chăm sóc sức khỏe. Việc kết hợp hai lĩnh vực này trong dự báo đóng vai trò then chốt giúp giảm thiểu rủi ro và tối ưu hóa quy trình điều trị. Khả năng sử dụng các mô hình toán học để nhận diện sớm nguy cơ dựa trên các chỉ số sinh học không chỉ mang ý nghĩa hàn lâm mà còn là công cụ hỗ trợ đắc lực cho các bác sĩ trong việc can thiệp kịp thời, đặc biệt với các bệnh lý diễn tiến âm thầm.
Suy thận mạn tính (Chronic Kidney Disease - CKD) là một trong những nguyên nhân hàng đầu gây gánh nặng bệnh tật trên toàn thế giới. Các yếu tố nguy cơ dẫn đến suy giảm chức năng thận thường đan xen phức tạp giữa đặc điểm nhân khẩu học (tuổi tác, thể trạng) và các chỉ số sinh hóa máu, nước tiểu cùng tiền sử bệnh lý (tăng huyết áp, tiểu đường). Vì vậy, việc định lượng mức độ suy giảm chức năng thận và dự báo xác suất mắc bệnh là một bài toán thực tiễn vô cùng cấp thiết.
Tiểu luận này được thực hiện nhằm áp dụng các kiến thức về xác suất thống kê nâng cao và các mô hình hồi quy để giải quyết bài toán trên qua tập dữ liệu “Kidney Disease Dataset”. Cụ thể, bài viết tập trung vào hai phương pháp cốt lõi:
Hồi quy tuyến tính (Linear Regression): Được sử dụng để phân tích sự phụ thuộc của mức lọc cầu thận ước tính (eGFR) – chỉ số vàng đánh giá chức năng thận – vào các yếu tố như độ tuổi, nồng độ Creatinine, Cystatin C và các chỉ số sinh hóa khác. Mô hình này giúp xác định mức độ ảnh hưởng của từng yếu tố lên sự suy giảm chức năng thận định lượng.
Hồi quy Logistic (Logistic Regression): Đây là công cụ quan trọng nhất để dự báo biến nhị phân (phân loại một cá nhân là “Khỏe mạnh” hay “Có nguy cơ mắc bệnh”). Mô hình này giúp thiết lập mối quan hệ giữa các biến độc lập và xác suất xảy ra biến cố bệnh lý, từ đó đưa ra những cảnh báo sớm dựa trên ngưỡng xác suất cụ thể.
Thông qua việc sử dụng phần mềm R và định dạng R Markdown, em sẽ thực hiện phân tích tập dữ liệu “Kidney Disease Dataset” để hiện thực hóa quy trình từ khai phá dữ liệu, xây dựng mô hình đến kiểm định giả thuyết toán học. Mục tiêu cuối cùng là tìm ra một mô hình dự báo đủ tốt để nâng cao khả năng cảnh báo sớm nguy cơ suy thận dựa trên các dữ liệu thống kê tin cậy.
2 CƠ SỞ TOÁN HỌC
Chương này trình bày các nguyên lý toán học nền tảng của các mô hình hồi quy. Các định nghĩa, định lý và chứng minh được tổng hợp từ tài liệu chuyên khảo của Sheldon M. Ross [1], giáo trình bài giảng của TS. Đào Huy Cường [3] và phương pháp ứng dụng của GS.Nguyễn Văn Tuấn [2].
2.1 Mô hình hồi quy tuyến tính
Nhiều vấn đề kỹ thuật và khoa học quan tâm đến việc xác định mối quan hệ giữa một tập hợp các biến số. Trong nhiều tình huống, có một biến đầu ra đơn lẻ \(Y\), còn được gọi là biến phụ thuộc, phụ thuộc vào giá trị của một tập hợp các biến đầu vào, còn được gọi là biến giải thích, \(x_1, \dots, x_r\). Loại quan hệ đơn giản nhất giữa biến phụ thuộc \(Y\) và các biến đầu vào \(x_1, \dots, x_r\) là mối quan hệ tuyến tính. Nghĩa là, với các hằng số \(\beta_0, \beta_1, \dots, \beta_r\), phương trình sau sẽ được thỏa mãn:
\[ Y = \beta_0 + \beta_1 x_1 + \cdots + \beta_r x_r \tag{1} \]
Nếu đây là mối quan hệ giữa \(Y\) và các \(x_i, i = 1, \dots, r\), thì (khi các \(\beta_i\) đã được xác định) chúng ta có thể dự đoán chính xác phản hồi cho bất kỳ tập hợp giá trị đầu vào nào. Tuy nhiên, trên thực tế, độ chính xác như vậy gần như không bao giờ đạt được, và điều tốt nhất người ta có thể mong đợi là phương trình (1) sẽ có thêm phần sai số ngẫu nhiên. Điều này có nghĩa là mối quan hệ rõ ràng được biểu diễn như sau:
\[ Y = \beta_0 + \beta_1 x_1 + \cdots + \beta_r x_r + \epsilon \tag{2} \]
trong đó \(e\), đại diện cho sai số ngẫu nhiên, được giả định là một biến ngẫu nhiên có trung bình bằng 0. Ngoài ra, một cách khác để biểu diễn phương trình (2) là:
\[ E[Y|\mathbf{x}] = \beta_0 + \beta_1 x_1 + \cdots + \beta_r x_r \]
trong đó \(\mathbf{x} = (x_1, \dots, x_r)\) là tập hợp các biến giải thích, và \(E[Y|\mathbf{x}]\) là kỳ vọng thu được với biến đầu vào \(\mathbf{x}\).
Phương trình (2) được gọi là phương trình hồi quy tuyến tính. Chúng ta nói rằng nó mô tả sự hồi quy của \(Y\) trên tập hợp các biến giải thích \(x_1, \dots, x_r\). Các đại lượng \(\beta_0, \beta_1, \dots, \beta_r\) được gọi là hệ số hồi quy (regression coefficients), và thường phải được ước lượng từ một tập dữ liệu.
- Một phương trình hồi quy chứa một biến giải thích duy nhất — tức là, trường hợp \(r = 1\) — được gọi là phương trình hồi quy đơn (simple regression equation).
- Trong khi đó, phương trình chứa nhiều biến giải thích được gọi là phương trình hồi quy đa biến (multiple regression equation).
Từ đây, ta có thể định nghĩa một mô hình hồi quy tuyến tính đơn giản như sau.
Định nghĩa ([1, Mục 9.1, Trang 358]): Giả sử ta có \(n\) cặp dữ liệu \((x_1, Y_1), (x_2, Y_2), \dots, (x_n, Y_n)\), trong đó \(x_i\) là biến giải thích (hay biến đầu vào) và \(Y_i\) là biến phụ thuộc (hay biến đầu ra). Mô hình hồi quy tuyến tính giả định mối quan hệ: \[ Y_i = \alpha + \beta x_i + \epsilon_i, \quad i = 1, \dots, n \tag{3} \] Trong đó:
\(\alpha, \beta\): Là các tham số của mô hình (lần lượt là hệ số chặn và hệ số góc).
\(\epsilon_i\): Là các sai số ngẫu nhiên độc lập, có phân phối chuẩn với kỳ vọng \(E[\epsilon_i] = 0\) và phương sai không đổi \(Var(\epsilon_i) = \sigma^2\).
2.1.1 Mô hình hồi quy tuyến tính đơn biến
Theo Sheldon M. Ross [1, Trang 357], mục tiêu của phân tích hồi quy là sử dụng các dữ liệu quan sát để ước lượng hàm số mô tả mối quan hệ giữa biến đầu vào \(x\) và biến phản hồi \(Y\).
Định nghĩa ([3, trang 28]). Giả sử ứng với mỗi giá trị của biến đầu vào \(x\) (còn gọi là biến giải thích) ta có biến đầu ra \(Y\) (còn gọi là biến phụ thuộc), và ta muốn mô hình hóa giá trị trung bình \(E[Y]\). Mô hình hồi quy tuyến tính đơn biến (simple linear regression model) giả định rằng \(E[Y]\) có quan hệ tuyến tính với \(x\), cụ thể là \[E[Y] = \alpha +\beta x\], trong đó \(\alpha\) và \(\beta\) được gọi là các hệ số hồi quy
Tiếp theo, ta quan tâm đến các ước lượng hệ số hồi quy.
2.1.1.1 Ước lượng các hệ số hồi quy
Để các suy diễn thống kê từ mô hình có giá trị, Sheldon Ross [1, Trang 356] thiết lập các giả định tiên quyết đối với sai số trong (3):
Quay lại với (3), ta giả sử \(\hat{\alpha}\) là ước lượng của hệ số \(\alpha\) và \(\hat{\beta}\) là ước lượng của hệ số \(\beta\), thì giá trị ước lượng của biến đầu ra tại \(x_i\) sẽ là \[\hat{y}_i = \hat{\alpha} + \hat{\beta}x_i. \tag{4} \]
Vì giá trị quan sát được của \(Y\) tại \(x_i\) là \(Y_i\) nên ta có các phần dư (residuals) giữa giá trị quan sát được và giá trị ước lượng của biến đầu ra là \[ r_i = Y_i - \hat{y}_i, \quad i = 1, 2, \dots, n. \tag{4} \]
Ta gọi \[ \text{SSR} = \displaystyle\sum_{i=1}^{n} r_i^2 \tag{5} \] là tổng bình phương các phần dư (sum of squares of the residuals)} giữa giá trị quan sát được và giá trị ước lượng của biến đầu ra.
Phương pháp bình phương tối thiểu (least squares estimator)} là phương pháp chọn các ước lượng \(\hat{\alpha}\) và \(\hat{\beta}\) sao cho tổng bình phương \(\text{SSR}\) trong (5) đạt giá trị nhỏ nhất.
Một số đặc trưng của \(\epsilon\) được Sheldon M. Ross [1, Trang 357] chỉ ra như sau
Kỳ vọng bằng không: \(E[\epsilon] = 0\). Điều này dẫn đến \(E[Y|x] = \alpha + \beta x\).
Phương sai không đổi: \(Var(\epsilon) = \sigma^2\) với mọi giá trị của \(x\).
Tính độc lập và phân phối chuẩn: Các sai số \(\epsilon_i\) là các biến ngẫu nhiên độc lập và tuân theo phân phối chuẩn \(\epsilon \sim N(0, \sigma^2)\).
Các ước lượng trên xoay quanh phương pháp ước lượng bình phương tối thiếu như sau.
\(\hat{\alpha}=\overline{Y}-\dfrac{S_{xY}}{S_{xx}}\overline{x}, \quad \hat{\beta}=\dfrac{S_{xY}}{S_{xx}}\) (2.5)
trong đó\(S_{xY}=\displaystyle\sum_{i=1}^{n}(x_{i}-\overline{x})(Y_{i}-\overline{Y}), \quad S_{xx}=\displaystyle\sum_{i=1}^{n}(x_{i}-\overline{x})^{2}.\) (2.6)
\(\dfrac{\partial SSR}{\partial \hat{\alpha}} = -2 \displaystyle\sum_{i=1}^{n}(Y_{i} - \hat{\alpha} - \hat{\beta}x_{i})\)
\(\dfrac{\partial SSR}{\partial \hat{\beta}} = -2 \displaystyle\sum_{i=1}^{n}x_{i}(Y_{i} - \hat{\alpha} - \hat{\beta}x_{i})\)
Giải hệ phương trình bằng cách cho các đạo hàm bằng 0:\(\displaystyle\sum_{i=1}^{n}Y_{i} = n\hat{\alpha} + \hat{\beta}\sum_{i=1}^{n}x_{i}\)
\(\displaystyle\sum_{i=1}^{n}x_{i}Y_{i} = \hat{\alpha}\displaystyle\sum_{i=1}^{n}x_{i} + \hat{\beta}\sum_{i=1}^{n}x_{i}^{2}\)
Từ phương trình đầu ta có \(\hat{\alpha} = \overline{Y} - \hat{\beta}\overline{x}\). Thay vào phương trình thứ hai và biến đổi đại số, ta thu được kết quả như trong (2.5) và (2.6). Ma trận Hess của \(SSR\) là một ma trận xác định dương nên điểm dừng này chính là điểm cực tiểu.
\(y=\hat{\alpha}+\hat{\beta}x\) (2.7)
là đường thẳng hồi quy ước lượng.
Vậy là ta đã xây dựng được các ước lượng hồi quy dựa vào phương pháp bình phương tối thiếu. Khi đó ta quan tâm thêm đến các phân phối của chúng có đặc trưng gì.
2.1.1.2 Phân phối các ước lượng hệ số hồi quy
Định lý 2.2. ([Định lý 3.2, trang 31]) Giả sử \(Y_{i}\sim\mathcal{N}(\alpha+\beta
x_{i},\sigma^{2}), i=1,2,...,n\) là các biến ngẫu nhiên độc lập.
Khi đó:
(a) \(\hat{\beta}\sim\mathcal{N}(\beta,\dfrac{\sigma^{2}}{S_{xx}})\)
(b) \(\hat{\alpha}\sim\mathcal{N}(\alpha,\dfrac{\sigma^{2}\sum
x_{i}^{2}}{nS_{xx}})\)
\(\hat{\beta} = \dfrac{\sum (x_i - \overline{x})Y_i}{S_{xx}}\)
Vì các \(Y_i \sim \mathcal{N}(\alpha +
\beta x_i, \sigma^2)\) độc lập, nên \(\hat{\beta}\) là tổng các biến ngẫu nhiên
phân phối chuẩn, do đó cũng có phân phối chuẩn. Kỳ vọng \(E[\hat{\beta}] = \beta\) và phương sai
\(Var(\hat{\beta}) =
\dfrac{\sigma^2}{S_{xx}}\).
(b) Tương tự, \(\hat{\alpha} = \overline{Y} -
\hat{\beta}\overline{x}\) cũng là tổ hợp tuyến tính của các \(Y_i\), do đó có phân phối chuẩn với kỳ vọng
\(\alpha\) và phương sai \(\dfrac{\sigma^2 \sum x_i^2}{n
S_{xx}}\).
Định lý 2.3. ([Định lý 3.3, trang 32]) Giả sử \(Y_{i}\sim\mathcal{N}(\alpha+\beta
x_{i},\sigma^{2}), i=1,2,...,n\) là các biến ngẫu nhiên độc lập.
Khi đó:
(a) \(SSR=S_{YY}-\dfrac{(S_{xY})^{2}}{S_{xx}}.\)
(b) \(\dfrac{SSR}{\sigma^{2}}\sim\chi_{n-2}^{2}.\)
Chứng minh. (Chứng minh chi tiết được trình bày trong [3, trang 32-34]), ta tóm tắt lại cách chứng minh như sau.
- Thay \(\hat{\alpha} = \overline{Y} -
\hat{\beta}\overline{x}\) vào \(SSR =
\displaystyle\sum_{i=1}^n (Y_i - \hat{\alpha} -
\hat{\beta}x_i)^2\), ta được:
\(SSR = \displaystyle\sum_{i=1}^n [(Y_i - \overline{Y}) - \hat{\beta}(x_i - \overline{x})]^2 = \displaystyle\sum (Y_i - \overline{Y})^2 - 2\hat{\beta}\displaystyle\sum (x_i - \overline{x})(Y_i - \overline{Y}) + \hat{\beta}^2 \displaystyle\sum (x_i - \overline{x})^2\)
\(SSR = S_{YY} - 2\hat{\beta} S_{xY} + \hat{\beta}^2 S_{xx}\)
Thay \(\hat{\beta} = S_{xY}/S_{xx}\) vào ta được \(SSR = S_{YY} - \dfrac{(S_{xY})^2}{S_{xx}}\).
(b) Bằng việc sử dụng biến đổi trực giao hoặc phân tích vector, có thể chứng minh được \(\dfrac{SSR}{\sigma^2} = \displaystyle\sum_{i=1}^n Z_i^2 - Z_1^2 - Z_2^2\) với \(Z_i \sim \mathcal{N}(0,1)\) độc lập, dẫn đến \(\dfrac{SSR}{\sigma^2} \sim \chi^2_{n-2}\).
Nhận xét 2.1. ([3, trang 34]) Vì \(\dfrac{SSR}{\sigma^{2}}\sim\chi_{n-2}^{2}\) nên ta suy ra \(E[SSR/(n-2)]=\sigma^{2}\) do đó \(SSR/(n-2)\) là một ước lượng “không chệch” cho phương sai chung \(\sigma^{2}\) của các sai số \(\epsilon_{i}.\)
\(\hat{\sigma}=\sqrt{\dfrac{SSR}{n-2}}\) (2.12)
là sai số chuẩn phần dư (residual standard error).
\(\dfrac{\hat{\alpha}-\alpha}{\hat{\sigma}\sqrt{\displaystyle\sum x_{i}^{2}}/\sqrt{nS_{xx}}}\sim t_{n-2}, \quad \dfrac{\hat{\beta}-\beta}{\hat{\sigma}/\sqrt{S_{xx}}}\sim t_{n-2}.\) (2.13)
\(r_{i}\sim\mathcal{N}(0,\sigma^{2}(1-1/n-(x_{i}-\overline{x})^{2}/S_{xx}))\) \(i=1,...,n\) (2.14)
\(r_i = \epsilon_i - \overline{\epsilon} - \dfrac{(x_i - \overline{x})\displaystyle\sum (x_j - \overline{x})\epsilon_j}{S_{xx}}\)
Vì đây là tổ hợp tuyến tính của các \(\epsilon_i \sim \mathcal{N}(0, \sigma^2)\), nên \(r_i\) có phân phối chuẩn. Kỳ vọng \(E[r_i] = 0\). Bằng cách tính \(Var(r_i) = Var(Y_i) - Var(\hat{y}_i)\) (vì \(Cov(r_i, \hat{y}_i) = 0\)), ta thu được phương sai như trong công thức (2.14).
Nhận xét 2.3. ([3, trang 36]) Các phần mềm thống kê thường cung cấp biểu đồ Residuals vs Fitted của các phần dư \(r_{i}\) với các giá trị ước lượng của biến đầu ra \(\hat{y_{i}}\) để giúp người làm thống kê đánh giá trực quan các giả định của mô hình hồi quy tuyến tính.
\(\dfrac{r_{i}}{\hat{\sigma}\sqrt{1-1/n-(x_{i}-\overline{x})^{2}/S_{xx}}}, \quad i=1,...,n\) (2.15)
là các phần dư chuẩn hóa (standardized residuals).
\(\dfrac{r_{i}}{\hat{\sigma}\sqrt{1-1/n-(x_{i}-\overline{x})^{2}/S_{xx}}}\sim t_{n-2},\) \(i=1,...,n.\) (2.16)
Việc xây dựng các mô hình hồi quy tuyến tính luôn cần quan tân đến các sai số. Để phản ánh một cách có cơ sở về tính đúng đắn của một mô hình hồi quy, ta quan tâm đến hệ số xác định.
2.1.1.3 Hệ số xác định
Theo TS. Đào Huy Cường[3], khi xét mô hình hồi quy tuyến tính đơn biến (2.1). Ta biết rằng đại lượng \(S_{YY} = \displaystyle\sum_{i=1}^{n}(Y_{i}-\overline{Y})^{2}\) đo lường độ biến thiên của các \(Y_{i}\). Chẳng hạn, nếu \(S_{YY}=0\) thì các \(Y_{i}\) đều bằng nhau. Hơn nữa, ta có:
\(\displaystyle\sum_{i=1}^{n}(Y_{i}-\overline{Y})^{2}=\sum_{i=1}^{n}(Y_{i}-\hat{y_{i}})^{2}+\sum_{i=1}^{n}(\hat{y_{i}}-\overline{Y})^{2}\) (2.17)
\(S_{YY}=SSR+SSE\) (2.18)
\(SSE=\displaystyle\sum_{i=1}^{n}(\hat{y_{i}}-\overline{Y})^{2}\) (2.19)
Như vậy, độ biến thiên của các \(Y_{i}\) được giải thích bởi 2 yếu tố:
- Yếu tố 1 (SSE): được giải thích bởi các giá trị ước lượng \(\hat{y}_{i}\) từ mô hình hồi quy.
- Yếu tố 2 (SSR): tổng bình phương các phần dư \(r_{i}\).
Ta gọi SSE (explained sum of squares) là độ biến thiên của \(Y\) được giải thích bởi mô hình hồi quy tuyến tính đơn biến (2.1). Như vậy, tỉ số \(\dfrac{SSE}{S_{YY}}\) cho biết có bao nhiêu phần trăm độ biến thiên của \(Y\) được giải thích bởi mô hình hồi quy tuyến tính đơn biến (2.1).
Từ đây, ta xây dựng đại lượng như sau:
\(R^{2}=\dfrac{SSE}{S_{YY}}=1-\dfrac{SSR}{S_{YY}}\) (2.20)
được gọi là hệ số xác định của mô hình hồi quy tuyến tính đơn biến.
Nhận xét 2.5. ([3, trang 38]) Giá trị của \(R^2\) là một chỉ báo về mức độ phù hợp toàn cục(goodness of fit) của một mô hình hồi quy tuyến tính đơn biến đối với tapapj số liệu cặp đôi \((x_i, Y_i)\). Giá trị của \(R^2\) luôn nằm giữa 0 và 1. Với giá trị gần 1 cho biết mức độ phù hợp tốt; với giá trị gầ 0 cho biết mức độ phù hợp kém.
Ta thấy rằng \(R^2\) luôn tăng khi thêm biến đầu vào, điều này dễ gây ảo tưởng tằng mô hình hồi quy đã phù hợp. Để khắc phcuj, người ta hiệu chỉnh \(R^2\) theo bậc tự do của SSR và \(S_{YY}\).
\(R_{adj}^{2}=1-\dfrac{SSR/(n-2)}{S_{YY}/(n-1)}\) (2.21)
được gọi là hệ số xác định hiệu chỉnh (Adjusted R-squared) của mô hình hồi quy tuyến tính đơn biến.
\(F=\dfrac{SSE/1}{SSR/(n-2)}\sim F_{1,n-2}.\) (2.22)
Như vậy kiểm định F cũng được sử dụng để kiểm định sự phù hợp toàn cục của mô hình hồi quy tuyến tính đơn biến.
Nhận xét 2.7. ([3, trang 39]) Hệ số xác định \(R^{2}\) bằng bình phương hệ số tương quan mẫu \(r_{x,Y}\): \(R^{2}=(r_{x,Y})^{2}.\)
2.1.2 Mô hình hồi quy tuyến tính đa thức
Định nghĩa ([2, mục 10.4, trang 174]). mô hình hối qui đa thức tìm mối liên hệ giữa biến phụ thuộc \(Y\) và biến độc lập \(X\) theo những hàm số sau đây: \[ Y = \alpha + \beta_1 X + \beta_2 X^2 + \beta_3 X^3 + .. + \beta_p X^p + \epsilon. \] Trong đó các thông số \(\beta_j (j = 1, 2, 3, ... p)\) là hệ số đo lường mối liên hệ giữa \(Y\) và \(X\); và \(\epsilon\) là phẩn dư của mô hình, với giả định \(\epsilon\) tuân theo luật phán phối chuẩn với trung bình 0 và phương sai \(\sigma^2\). Cho một dãy cập số \((Y_1,X_1), (Y_2,X_2), (Y_3,X_3)), (Y_n, X_n),\) chúng ta có thể áp dụng phương pháp bình phương nhô nhất để ước tính \(\beta\) và \(\sigma^2\)
Về phương pháp, mô hình này cũng sử dụng phương pháp bình phương tối thiểu.
Theo Sheldon Ross [2, mục 9.9], các ước lượng bình phương tối thiểu \(\beta_0, \beta_1, \ldots, \beta_r\) là những giá trị tối thiểu hóa:
\[ \sum_{i=1}^{n} \left(Y_i - \beta_0 - \beta_1 x_i - \beta_2 x_i^2 - \cdots - \beta_r x_i^r\right)^2 \]
Để xác định các ước lượng này, ta lấy đạo hàm riêng theo \(\beta_0, \ldots, \beta_r\) và cho bằng 0 để tìm các giá trị cực tiểu. Sắp xếp lại các phương trình thu được dẫn tới hệ phương trình chuẩn:
\[ \sum_{i=1}^{n} Y_i = \beta_0 n + \beta_1 \sum_{i=1}^{n} x_i + \beta_2 \sum_{i=1}^{n} x_i^2 + \cdots + \beta_r \sum_{i=1}^{n} x_i^r \]
\[ \sum_{i=1}^{n} x_i Y_i = \beta_0 \sum_{i=1}^{n} x_i + \beta_1 \sum_{i=1}^{n} x_i^2 + \beta_2 \sum_{i=1}^{n} x_i^3 + \cdots + \beta_r \sum_{i=1}^{n} x_i^{r+1} \]
\[ \sum_{i=1}^{n} x_i^2 Y_i = \beta_0 \sum_{i=1}^{n} x_i^2 + \beta_1 \sum_{i=1}^{n} x_i^3 + \cdots + \beta_r \sum_{i=1}^{n} x_i^{r+2} \]
\[ \vdots \]
\[ \sum_{i=1}^{n} x_i^r Y_i = \beta_0 \sum_{i=1}^{n} x_i^r + \beta_1 \sum_{i=1}^{n} x_i^{r+1} + \cdots + \beta_r \sum_{i=1}^{n} x_i^{2r} \]
Khi khớp một đa thức cho một tập dữ liệu các cặp \((x, y)\), ta thường có thể xác định bậc cần thiết của đa thức bằng cách quan sát biểu đồ phân tán. Nên dùng bậc thấp nhất có thể nhưng vẫn mô tả dữ liệu một cách đầy đủ. Ta cũng không nên dùng đa thức đã khớp để dự đoán ở các mức đầu vào quá xa so với vùng dữ liệu đã dùng để khớp, vì đa thức có thể chỉ hợp lệ trong vùng lân cận các mức đầu vào đã sử dụng.
Một số nguyên tắc trong xây dựng mô hình đó là nguyên tắc chọn số mũ càng nhỏ càng tốt, do đó trong thực tế ta sẽ bắt đầu với bậc 1 rồi tăng dần lên.
2.1.3 Mô hình hồi quy tuyến tính đa biến
Định nghĩa ([3, trang 39]). Giả sử ứng với mỗi bộ giá trị của các biến đầu vào \(x_1, x_2,..., x_k\) (còn gọi là biến giải thích) ta có biến đầu ra \(Y\) (còn gọi là biến phụ thuộc), và ta muoonc mô hình hóa giá trị trung bình \(E[Y]\). Mô hình hồ quy tuyến tính đa biến giả đinh rằng \(E[Y]\) có quan hệ tuyến tính với \(x_2, x_2,..., x_k\): \[ E[Y]= \beta_0 +\beta_1 x_1 +\beta_2 x_2 +...+ \beta_k x_k \tag{2.23} \]
2.1.3.1 Ước lượng các hệ số hồi quy
Giả sử \(Y_i, i=\overline{1,n}\) là \(n\) giá trị quan sát được với bộ đầu vào \(x_{ij}, j= \overline{1,k}\). Khi đó ta xây dựng được các đại lượng sau:
\[ SSR= \displaystyle\sum_{i=1}^n r_i^2, \tag{2.24} \]
trong đó \(r_i- Y_i-\hat{y_i}, i\overline{1,n}\) là các phần dư giữa giá rị quna sát được và giá trị ước lượng của biến \(Y\) \[ \hat{y_i}= \hat{\beta_0} +\hat{\beta_1} x_{i1} +\hat{\beta_2} x_{i2} +...+ \hat{\beta_k} x_{ik}, \tag{2.26} \]\(\hat{\beta}=(X^{T}X)^{-1}X^{T}Y\) (2.27)
nếu \(rank(X) = k + 1\), trong đó\(\hat{\beta}=\begin{bmatrix}\hat{\beta}_{0}\\ \hat{\beta}_{1}\\ \vdots\\ \hat{\beta}_{k}\end{bmatrix}, \quad X=\begin{bmatrix}1&x_{11}&x_{12}&\cdot\cdot\cdot&x_{1k}\\ 1&x_{21}&x_{22}&\cdot\cdot\cdot&x_{2k}\\ \vdots&\vdots&\vdots&\vdots\\ 1&x_{n1}&x_{n2}&\cdot\cdot\cdot&x_{nk}\end{bmatrix}, \quad Y=\begin{bmatrix}Y_{1}\\ Y_{2}\\ \vdots\\ Y_{n}\end{bmatrix}.\) (2.28)
\(SSR=(Y-X\hat{\beta})^{T}(Y-X\hat{\beta}) = Y^{T}Y-2\hat{\beta}^{T}X^{T}Y+\hat{\beta}^{T}X^{T}X\hat{\beta}.\)
Vector gradient của \(SSR\) là\(\nabla SSR=-2X^{T}Y+2X^{T}X\hat{\beta}.\)
Vì \(rank(X)=k+1\) nên \(X^{T}X\) khả nghịch, do đó ta tìm được duy nhất một điểm dừng \(\hat{\beta}=(X^{T}X)^{-1}X^{T}Y\). Hơn nữa, Hessian của \(SSR\) là\(\nabla^{2}SSR=2X^{T}X\)
đó là ma trận xác định dương. Vậy \(SSR\) đạt giá trị nhỏ nhất tại \(\hat{\beta}=(X^{T}X)^{-1}X^{T}Y\).
Như vậy, ta đã có xây dựng một mô hình hồi quy đa biến. Sau đây ta quan tâm đên scacs đặc trưng về phân phối của các ước lượng hệ số hồi quy.
2.1.3.2 Phân phối của các ước lượng hệ số hồi quy
\(Y_{i}\sim\mathcal{N}(\beta_{0}+\beta_{1}x_{i1}+\beta_{2}x_{i2}+\cdot\cdot\cdot+\beta_{k}x_{ik},\sigma^{2}), \quad i=1,2,...,n\)
là các biến ngẫu nhiên độc lập. Khi đó\(\hat{\beta}_{j}\sim\mathcal{N}(\beta_{j},\sigma^{2}[(X^{T}X)^{-1}]_{j+1,j+1}) \quad j=0,1,...,k.\) (2.29)
\(\epsilon_{i}=Y_{i}-(\beta_{0}+\beta_{1}x_{i1}+\beta_{2}x_{i2}+\cdot\cdot\cdot+\beta_{k}x_{ik}),\)
thì \(\epsilon_{i}\sim\mathcal{N}(0,\sigma^2)\) là các biến ngẫu nhiên độc lập. Từ (2.27) và (2.30) ta có\(\hat{\beta}=(X^{T}X)^{-1}X^{T}Y = (X^{T}X)^{-1}X^{T}(X\beta+\epsilon) = \beta+(X^{T}X)^{-1}X^{T}\epsilon,\)
trong đó\(\beta=\begin{bmatrix}\beta_{0}\\ \beta_{1}\\ \vdots\\ \beta_{k}\end{bmatrix}, \quad \epsilon=\begin{bmatrix}\epsilon_{1}\\ \epsilon_{2}\\ \vdots\\ \epsilon_{n}\end{bmatrix}\) (2.30)
Do đó\(Cov(\hat{\beta})=Cov((X^{T}X)^{-1}X^{T}\epsilon) = E[(X^{T}X)^{-1}X^{T}\epsilon((X^{T}X)^{-1}X^{T}\epsilon)^{T}]\)
\(= (X^{T}X)^{-1}X^{T}E[\epsilon\epsilon^{T}]((X^{T}X)^{-1}X^{T})^{T} = (X^{T}X)^{-1}X^{T}(\sigma^{2}I)X(X^{T}X)^{-1}\)
\(= \sigma^{2}(X^{T}X)^{-1}X^{T}X(X^{T}X)^{-1} = \sigma^{2}(X^{T}X)^{-1}.\)
Vậy \(\hat{\beta}_{j}\sim\mathcal{N}(\beta_{j},\sigma^{2}[(X^{T}X)^{-1}]_{j+1,j+1})\) với mọi \(j=0,1,...,k\).
\(Y_{i}\sim\mathcal{N}(\beta_{0}+\beta_{1}x_{i1}+\beta_{2}x_{i2}+\cdot\cdot\cdot+\beta_{k}x_{ik},\sigma^{2}), \quad i=1,2,...,n\)
là các biến ngẫu nhiên độc lập. Khi đó:
(a) \(SSR=Y^{T}(I-H)Y\) trong đó \(H=X(X^{T}X)^{-1}X^{T}.\)
(b) \(\dfrac{SSR}{\sigma^{2}}\sim\chi_{n-k-1}^{2}.\)
\(SSR=(Y-X\hat{\beta})^{T}(Y-X\hat{\beta}) = Y^{T}Y-2\hat{\beta}^{T}X^{T}Y+\hat{\beta}^{T}X^{T}X\hat{\beta}\)
\(= Y^{T}Y-2\hat{\beta}^{T}X^{T}Y+\hat{\beta}^{T}X^{T}X(X^{T}X)^{-1}X^{T}Y = Y^{T}Y-\hat{\beta}^{T}X^{T}Y\)
\(= Y^{T}Y-((X^{T}X)^{-1}X^{T}Y)^{T}X^{T}Y = Y^{T}Y-Y^{T}X(X^{T}X)^{-1}X^{T}Y\)
\(= Y^{T}(I-X(X^{T}X)^{-1}X^{T})Y = Y^{T}(I-H)Y.\)
- Dễ thấy \(H\) và \(I-H\) lần lượt là phép chiếu trực giao lên
không gian \(col(X)\) và \(col(X)^{\perp}\). Do đó ta có thể chọn một
ma trận trực giao \(U=[U_{1},U_{2}]\)
trong đó \(U_{1}\) gồm các cột là cơ sở
trực chuẩn của không gian \(col(X)\),
\(U_{2}\) gồm các cột là cơ sở trực
chuẩn của không gian \(col(X)^{\perp}\)
sao cho phép chiếu \(I-H\) khi viết
trong cơ sở \(U\) trở thành ma trận
\(\begin{bmatrix}I_{n-k-1}&0\\
0&0_{k+1}\end{bmatrix}\). Đặt \(W=U^{T}\epsilon\). Vì \(U^{T}U=I\) nên ta có \(Cov(W)=\sigma^{2}I\). Suy ra \(W_{i}\sim\mathcal{N}(0,\sigma^{2})\), \(i=1,...,n\) và là các biến độc lập. Ta có
\(SSR=Y^{T}(I-H)Y = (X\beta+\epsilon)^{T}(I-H)(X\beta+\epsilon) = \epsilon^{T}(I-H)\epsilon\)
\(= W^{T}\begin{bmatrix}I_{n-k-1}&0\\ 0&0_{k+1}\end{bmatrix}W = \displaystyle\sum_{i=1}^{n-k-1}W_{i}^{2}.\)
Vậy \(\dfrac{SSR}{\sigma^{2}}=\displaystyle\sum_{i=1}^{n-k-1}\left(\dfrac{W_{i}}{\sigma}\right)^{2}\sim\chi_{n-k-1}^{2}.\)
Để thuận tiện hơn khi làm việc với các sai số, ta xây dựng sai số chuẩn phần phần dư như sau:
\(\hat{\sigma}=\sqrt{\dfrac{SSR}{n-k-1}}\) (2.31)
là sai số chuẩn phần dư (residual standard error).
\(\dfrac{\hat{\beta}_{j}-\beta_{j}}{\hat{\sigma}\sqrt{[(X^{T}X)^{-1}]_{j+1,j+1}}}\sim t_{n-k-1}, \quad j=0,1,...,k\) (2.32)
Kết quả này thường được các phần mềm thống kê tóm tắt trong bảng hệ số hồi quy dưới đây, đồng thời giúp người làm thống kê sử dụng kiểm định \(t\) để đưa ra kết luận về các giả thuyết \(H_{0}:\beta_{j}=0\) với \(j=1,...,k\). Nếu \(P\) value của hệ số ước lượng \(\hat{\beta}_{j}\) nhỏ \((<0.05)\) thì ta đủ cơ sở để bác bỏ \(H_{0},\) tức là biến \(x_{j}\) có ảnh hưởng đáng kể đến \(Y\).
\(r_{i}\sim\mathcal{N}(0,\sigma^{2}(1-H_{i,i})), \quad i=1,...,n\) (2.33)
\(r = (I-H)(X\beta + \epsilon) = (I-H)\epsilon\).
Do đó \(r\) là biến đổi tuyến tính của vector ngẫu nhiên phân phối chuẩn \(\epsilon\), nên \(r\) có phân phối chuẩn với \(E[r] = (I-H)E[\epsilon] = 0\). Ma trận hiệp phương sai của \(r\) là\(Cov(r) = (I-H)Cov(\epsilon)(I-H)^T = \sigma^2(I-H)(I-H) = \sigma^2(I-H)\).
Suy ra \(Var(r_i) = \sigma^2(I-H)_{i,i} = \sigma^2(1-H_{i,i})\). Vậy \(r_i \sim \mathcal{N}(0, \sigma^2(1-H_{i,i}))\).
\(\dfrac{r_{i}}{\hat{\sigma}\sqrt{1-H_{i,i}}}, \quad i=1,...,n\) (2.34)
là các phần dư chuẩn hóa (standardized residuals).
\(\dfrac{r_{i}}{\hat{\sigma}\sqrt{1-H_{i,i}}}\sim t_{n-k-1}, \quad i=1,...,n.\) (2.35)
Kết quả này xấp xỉ phân phối chuẩn chính tắc \(\mathcal{N}(0,1)\) khi \(n \to \infty\).
Trong mô hình đa biến khi xuất hiện quá nhiều biến sẽ xảy ra vấn đề các biến thuộc tính chồng lập với nhau. Theo GS. Nguyễn Văn Tuấn [2, mục 10.5, trang 184] hiện tượng này gọi là đa cộng tuyến và gây ra các vấn đề sau:
Thông tin trùng lặp: Các biến có tương quan cao cung cấp một lượng thông tin như nhau để tiên đoán biến phụ thuộc, dẫn đến việc đưa cả hai vào mô hình là không cần thiết.
Làm mất ý nghĩa thống kê của biến: Khi một biến đã có mặt trong mô hình, ảnh hưởng của biến khác (có tương quan cao với nó) có thể không còn ý nghĩa thống kê (P>0.05), ngay cả khi biến đó thực sự có ảnh hưởng khi xét đơn lẻ. Điều này khiến việc diễn dịch vai trò của từng biến trở nên khó khăn.
Khi thực hiện trong R ta có thể dùng chỉ số VIF để đánh giá mo hình có đa cộng tuyến không. Với \(VIF<5\) thì mô hình tốt và đa cộng tuyến không ảnh hưởng nhiều. Với \(5 \leq VIF \leq 10\) thì ta cần xem xét cẩn thận và đánh giá lại. Với \(VIF > 10\) thì hiện tượng đa cộng tuyến ảnh hưởng mạnh tới mô hình và có biến nào đó bị phản ảnh hoàn toàn theo các biến còn lại trong mô hình.
Tương tự như mô hình hồi quy đơn biến, ta cũng xây dựng hệ số xác định của mô hình đa biến để đánh giá sự phù hợp của mô hình khi áp dụng.
2.1.3.3 Hệ số xác định
Trong mô hình hồi quy tuyến tính đa biến trên, ta chứng minh được
\[ S_{YY}=SSR+SSE \tag{2.36} \]
trong đó \(SSE= \displaystyle\sum_{i=1}^n (\hat{\beta_0}+\hat{\beta_1}x_{i1}+...+\hat{\beta_k}x_{ik}-\overline{Y})^2\)
\(R^{2} = 1 - \dfrac{SSR}{S_{YY}}\) (2.38)
được gọi là hệ số xác định (coefficient of determination) của mô hình hồi quy tuyến tính đa biến.
Tương tự như mô hình đơn biến, để xử lý vấn đề \(R^2\) tăng khi thêm biến đầu vào, ta xây dựng đại lượng hiệu chỉnh sau
\(R_{adj}^{2} = 1 - \dfrac{SSR/(n-k-1)}{S_{YY}/(n-1)}\) (2.39)
được gọi là hệ số xác định hiệu chỉnh (Adjusted R-squared).
\(\dfrac{S_{YY}}{\sigma^2}\sim \chi_{n-1}^2, \dfrac{SSE}{\sigma^2} \sim \chi_k^2\) (2.40)
\(F = \dfrac{(S_{YY}-SSR)/k}{SSR/(n-k-1)} \sim F_{k, n-k-1}\) (2.41)
Kiểm định này giúp đánh giá sự phù hợp toàn cục của mô hình hồi quy.
2.2 Mô hình hồi quy Logistic
Định nghĩa ([3, trang 46]). Giả sử ứng với mỗi bộ giá trị của biến đầu vào \(x_1, x_2,..., x_k\) ta có biến đầu ra \(Y\) chỉ nhận một trong hai giá trị 0 hoặc 1, và ta muốn mohinhf hóa xác xuất
\[ P(Y=1)=p \tag{2.42} \] Mô hình hồi quy logistic giả định rằng logit (log odds unit) của xác suất \(p\) có quan hệ tuyến tính với \(x_1, x_2, ..., x_k\)
\[ ln \dfrac{p}{1-p}= \beta_0+\beta_1x_1+...+\beta_k x_k \tag{2.43} \] Suy ra công thức của xác suất \(p\) là \[ p=\dfrac{1}{1+e^{-(\beta_0+\beta_1x_1+...+\beta_k x_k)}} \tag{2.44} \] Người ta còn gọi công thức (2.44) là hàm logistic.
2.2.1 Ước lượng các hệ số hồi quy
Giả sử ta quan sát được bảng giá trị:
| \(x_{1}\) | \(x_{2}\) | \(\dots\) | \(x_{k}\) | \(Y\) |
|---|---|---|---|---|
| \(x_{11}\) | \(x_{12}\) | \(\dots\) | \(x_{1k}\) | \(Y_{1}\) |
| \(x_{21}\) | \(x_{22}\) | \(\dots\) | \(x_{2k}\) | \(Y_{2}\) |
| \(\vdots\) | \(\vdots\) | \(\ddots\) | \(\vdots\) | \(\vdots\) |
| \(x_{n1}\) | \(x_{n2}\) | \(\dots\) | \(x_{nk}\) | \(Y_{n}\) |
và sử dụng bảng giá trị này để ước lượng các hệ số \(\beta_{0}, \beta_{1}, \dots, \beta_{k}\) trong mô hình hồi quy logistic:
\(ln\dfrac{p}{1-p} = \beta_{0} + \beta_{1}x_{1} + \dots + \beta_{k}x_{k}\)
Giả sử \(\hat{\beta}_{0}, \hat{\beta}_{1}, \dots, \hat{\beta}_{k}\) là các ước lượng của các hệ số \(\beta_{0}, \beta_{1}, \dots, \beta_{k}\), thì giá trị ước lượng của xác suất \(p\) tại bộ giá trị \((x_{i1}, \dots, x_{ik})\) của các biến đầu vào là:
\(\hat{p}_{i} = \dfrac{1}{1+e^{-(\hat{\beta}_{0} + \hat{\beta}_{1}x_{i1} + \dots + \hat{\beta}_{k}x_{ik})}}, \quad i=1,2,\dots,n\) (2.45)
Vì giá trị quan sát được của \(Y\) tại \((x_{i1}, \dots, x_{ik})\) là \(Y_{i}\) nên ta có hàm khả năng (likelihood) của các ước lượng \(\hat{\beta}_{0}, \hat{\beta}_{1}, \dots, \hat{\beta}_{k}\) là:
\(L(\hat{\beta}_{0}, \hat{\beta}_{1}, \dots, \hat{\beta}_{k}) = \displaystyle\prod_{i=1}^{n} \hat{p}_{i}^{Y_{i}}(1-\hat{p}_{i})^{1-Y_{i}}\) (2.46)
Phương pháp ước lượng khả năng cực đại (maximum likelihood estimator) là phương pháp chọn các ước lượng \(\hat{\beta}\) sao cho hàm khả năng đạt giá trị lớn nhất.
\(X^{T}(Y-\hat{p}) = 0\) (2.47)
trong đó:\(X = \begin{bmatrix} 1 & x_{11} & x_{12} & \dots & x_{1k} \\ 1 & x_{21} & x_{22} & \dots & x_{2k} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \dots & x_{nk} \end{bmatrix}, \hat{\beta} = \begin{bmatrix} \hat{\beta}_{0} \\ \hat{\beta}_{1} \\ \vdots \\ \hat{\beta}_{k} \end{bmatrix}, Y = \begin{bmatrix} Y_{1} \\ Y_{2} \\ \vdots \\ Y_{n} \end{bmatrix}, \hat{p} = \begin{bmatrix} \hat{p}_{1} \\ \hat{p}_{2} \\ \vdots \\ \hat{p}_{n} \end{bmatrix}\) (2.48)
Chứng minh. Lấy logarit tự nhiên của hàm khả năng (2.46) ta được hàm log-likelihood: \[\mathcal{L}(\hat{\beta}) = \displaystyle\sum_{i=1}^{n} (Y_{i} ln \hat{p}_{i} + (1 - Y_{i}) ln(1 - \hat{p}_{i}))\] \[= \displaystyle\sum_{i=1}^{n} (Y_{i} ln \dfrac{\hat{p}_{i}}{1 - \hat{p}_{i}} + ln(1 - \hat{p}_{i}))\] \[= \displaystyle\sum_{i=1}^{n} (Y_{i} X_{i, \bullet} \hat{\beta} - ln(1 + e^{X_{i, \bullet} \hat{\beta}}))\]
Vector gradient của \(\mathcal{L}(\hat{\beta})\) là:\(\nabla \mathcal{L}(\hat{\beta}) = \displaystyle\sum_{i=1}^{n} (Y_{i} - \hat{p}_{i}) X_{i, \bullet}^{T} = X^{T} (Y - \hat{p})\) (2.49)
Do đó điểm dừng \(\hat{\beta}\) là nghiệm của hệ phương trình \(X^{T}(Y - \hat{p}) = 0\). Hessian của \(\mathcal{L}(\hat{\beta})\) là:\(\nabla^{2} \mathcal{L}(\hat{\beta}) = -\displaystyle\sum_{i=1}^{n} \hat{p}_{i} (1 - \hat{p}_{i}) X_{i, \bullet}^{T} X_{i, \bullet} = -X^{T} \hat{W} X\) (2.50)
trong đó \(\hat{W} = diag(\hat{p}_{1}(1 - \hat{p}_{1}), \dots, \hat{p}_{n}(1 - \hat{p}_{n}))\). Vì \(X^{T} \hat{W} X\) là ma trận nửa xác định dương nên mọi điểm dừng (nếu có) đều là điểm cực đại toàn cục.
Nhận xét 2.12. ([3, trang 48]) Hệ phương trình điểm dừng (2.47) là phi tuyến nên người ta thường tìm nghiệm gần đúng bằng các phương pháp số như phương pháp Newton:\(\hat{\beta}^{(n+1)} = \hat{\beta}^{(n)} - [\nabla^{2} \mathcal{L}(\hat{\beta}^{(n)})]^{-1} \nabla \mathcal{L}(\hat{\beta}^{(n)})\) (2.51)
2.2.2 Phân phối của các ước lượng hệ số hồi quy
\(\hat{\beta_{j}}\sim\mathcal{N}(\beta_{j},[(X^{T}\hat{W}X)^{-1}]_{j+1,j+1}), \quad j=0,1,...,k\) (2.52)
\(\dfrac{\hat{\beta}_{j}-\beta_{j}}{\sqrt{[(X^{T}\hat{W}X)^{-1}]_{j+1,j+1}}}\sim\mathcal{N}(0,1), \quad j=0,1,...,k\) (2.53)
Kết quả này thường được các phần mềm thống kê tóm tắt trong bảng hệ số hồi quy, đồng thời giúp người làm thống kê sử dụng kiểm định Z để đưa ra kết luận về các giả thuyết \(H_{0}:\beta_{j}=0\) với \(j=1,...,k\). Nếu P value của hệ số ước lượng \(\hat{\beta}_{j}\) nhỏ \((<0.05)\) thì ta đủ cơ sở để bác bỏ \(H_{0}\) tức là biến \(x_{j}\) có ảnh hưởng đáng kể đến xác suất \(Y=1\).
\(d_{i}=sign(Y_{i}-\hat{p_{i}})\sqrt{-2(Y_{i}ln\hat{p_{i}}+(1-Y_{i})ln(1-\hat{p_{i}}))}, \quad i=1,...,n\) (2.54)
là các phần dư độ lệch (deviance residuals).
Nhận xét 2.14. ([3, trang 50]) Nếu mô hình dự đoán tốt thì \(\hat{p}_i \approx Y_i\), khi đó \(d_i \approx 0\). Nếu mô hình dự đoán sai nghiêm trọng thì \(\hat{p}_i \approx 0\) nhưng \(Y_i = 1\) hoặc ngược lại, khi đó \(|d_i|\) rất lớn. Dấu \(sign(Y_i - \hat{p}_i)\) giúp giữ thông tin “mô hình dự đoán thừa hay thiếu”, giống như các \(r_i\) của mô hình hồi quy tuyến tính. Như vậy ta có thể hiểu các \(d_i\) là các phần dư dựa trên hàm log-likelihood.
\(D = \displaystyle\sum_{i=1}^{n} d_i^2\) (2.55)
là độ lệch dư (residual deviance).
\(D = \displaystyle\sum_{i=1}^{n} d_i^2 = -2 \displaystyle\sum_{i=1}^{n} (Y_i \ln \hat{p}_i + (1 - Y_i) \ln(1 - \hat{p}_i)) = -2\mathcal{L}(\hat{\beta})\) (2.56)
\(D_0 = -2 \displaystyle\sum_{i=1}^{n} (Y_i \ln \hat{p} + (1 - Y_i) \ln(1 - \hat{p}))\) (2.57)
trong đó \(\hat{p} = \dfrac{1}{n} \displaystyle\sum_{i=1}^{n} Y_i\), là độ lệch dư của mô hình rỗng (null deviance).
Nhận xét 2.16. ([3, trang 51]) Độ độ lệch dư \(D_0\) dùng để đo mức độ sai khác của mô hình rỗng (mô hình chỉ có hệ số \(\beta_0\), không có biến giải thích) so với mô hình hoàn hảo. Ta có thể xem \(D_0\) như là mức “phù hợp tệ nhất” vì nó không dùng biến giải thích nào. Vì lý do này nên người ta thường dùng nó để làm mức xuất phát để so sánh với độ lệch dư của các mô hình có biến giải thích.
\(AIC = -2\mathcal{L}(\hat{\beta}) + 2(k + 1)\) (2.59)
là tiêu chuẩn thông tin Akaike (Akaike Information Criterion).
Nhận xét 2.17. (3, trang 51) Vì \(AIC = D + 2(k + 1)\) nên nó cũng là một đại lượng đo mức độ phù hợp dữ liệu của mô hình đang xét nhưng có tính đến hình phạt (penalty) nếu sử dụng quá nhiều biến giải thích. Như vậy AIC cân bằng giữa độ phù hợp và sự đơn giản. Nó thường được dùng khi so sánh nhiều mô hình trên cùng một bộ dữ liệu, mô hình nào có AIC nhỏ hơn thì được xem là tương đối tốt hơn.
Ta đã tìm hiểu các mô hình hồi quy thường dược áp dụng, tuy nhiên, khi ứng dụng vào thực tiễn ta cần xem thử mô hình ta xây dựng được có tốt, có hợp lý không. Phần cơ sở lý thuyết sau đây sẽ định hướng cho phân tích mẫu số liệu của chúng ta.
2.3 Kiểm định tính hợp lý
Theo TS. Đào Huy Cường [3, trang 39], các suy diễn thống kê và độ tin cậy của các hệ số ước lượng \(\hat{\beta}\) dựa trên giả định rằng sai số ngẫu nhiên \(\epsilon\) có kỳ vọng bằng 0 và phương sai \(\sigma^2\) không đổi. GS. Nguyễn Văn Tuấn [2, trang 268] nhấn mạnh rằng một mô hình có các hệ số mang ý nghĩa thống kê (\(p < 0.05\)) vẫn có thể là một mô hình không hợp lý nếu vi phạm các điều kiện về phân phối của phần dư. Do đó ta có thể thấy tính hợp lý của mô hình phụ thuộc khá nhiều vào yếu tố phần du. Mục này ta sẽ nghiên cứu về các đặc trưng của phần dư và cách kiểm định của chúng.
2.3.1 Kiểm định chuẩn của phần dư
Giả định về tính chuẩn của phần dư là một trong những điều kiện tiên quyết quan trọng nhất để các kết luận thống kê (như kiểm định \(t\) cho các hệ số hay kiểm định \(F\) cho toàn bộ mô hình) có giá trị.
1. Định nghĩa phần dư:
Dựa trên bài giảng của TS. Đào Huy Cường [3, Trang 29], với mỗi giá trị quan sát \(Y_i\) và giá trị dự báo tương ứng \(\hat{y}_i\), phần dư (residual) \(r_i\) được xác định bởi công thức: \[r_i = Y_i - \hat{y}_i\] Theo mô hình lý thuyết, ta giả định các sai số ngẫu nhiên \(\epsilon_i\) (và do đó là các phần dư \(r_i\)) phải tuân theo phân phối chuẩn với kỳ vọng bằng 0: \(r_i \sim N(0, \sigma^2)\).
2. Phương pháp kiểm tra trực quan bằng biểu đồ Q-Q (Quantile-Quantile Plot):
Theo hướng dẫn thực hành của GS. Nguyễn Văn Tuấn [2, mục 10.2.3], cách phổ biến và hiệu quả nhất để kiểm tra tính chuẩn là sử dụng biểu đồ Q-Q Plot.
Nguyên lý: Biểu đồ này so sánh các phân vị (quantiles) của dữ liệu phần dư thực tế với các phân vị của một phân phối chuẩn lý thuyết.
Cách sử dụng: Nếu các điểm dữ liệu trên biểu đồ nằm tập trung sát hoặc dọc theo đường thẳng 45 độ, ta có thể kết luận phần dư tuân theo phân phối chuẩn. Ngược lại, nếu các điểm dữ liệu uốn cong hoặc lệch xa đường thẳng, mô hình đã vi phạm giả định về tính chuẩn [3, trang 37].
3. Kiểm định thống kê Shapiro-Wilk:
Để bổ sung cho phương pháp hình ảnh, ta sử dụng kiểm định
Shapiro-Wilk (thường được thực hiện qua hàm shapiro.test()
trong R) nhằm đưa ra kết luận mang tính định lượng [2 trang 142].
Nguyên lý giả thuyết không (\(H_0\)): Phần dư có phân phối chuẩn.
Quy tắc sử dụng: Nếu trị số \(P\) (p-value) thu được lớn hơn mức ý nghĩa \(\alpha = 0.05\), ta không có đủ bằng chứng để bác bỏ \(H_0\), đồng nghĩa với việc giả định về tính chuẩn được thỏa mãn.
4. Hệ quả của việc vi phạm giả định:
Sheldon Ross [1, Mục 9.3, trang 361] nhấn mạnh rằng nếu phần dư không có phân phối chuẩn, các khoảng tin cậy và các trị số \(P\) của các hệ số hồi quy sẽ không còn chính xác, dẫn đến những kết luận sai lệch về mối quan hệ giữa các biến độc lập (như tuổi, BMI) và nguy cơ suy thận.
2.3.2 Kiểm định phương sai sai số không đổi
Giả định về phương sai không đổi (Homoscedasticity) đóng vai trò then chốt trong việc đảm bảo tính hiệu quả của các ước lượng bình phương tối thiểu (LSE).
1. Định nghĩa toán học: Theo Sheldon Ross [1, trang 361], trong mô hình hồi quy, ta giả định rằng phương sai của sai số ngẫu nhiên \(\epsilon_i\) là hằng số đối với mọi giá trị của biến độc lập \(x_i\): \[Var(\epsilon_i) = \sigma^2, \quad \forall i = 1, \dots, n\] Điều này có nghĩa là mức độ biến động (độ phân tán) của các quan sát xung quanh đường hồi quy phải đồng đều trên toàn bộ dải dữ liệu của các biến dự báo (như tuổi hay chỉ số BMI).
2. Phương pháp kiểm tra bằng biểu đồ Residuals vs Fitted: Dựa trên Nhận xét 2.3 đã trình bày, ta sử dụng biểu đồ giữa phần dư (Residuals) và giá trị dự báo (Fitted values) để kiểm chứng giả định này.
Trường hợp lý tưởng (Homoscedasticity): Các điểm phần dư phân bố một cách ngẫu nhiên và đồng đều trong một dải nằm ngang bao quanh đường \(y=0\). Không có bất kỳ hình mẫu hoặc xu hướng rõ rệt nào xuất hiện.
Trường hợp vi phạm (Heteroscedasticity): Nếu các điểm dữ liệu tạo thành hình phễu (phình ra hoặc thắt lại khi giá trị dự báo tăng lên), đây là dấu hiệu cho thấy phương sai thay đổi [1, Mục 9.9].
3. Tầm quan trọng và hệ quả của việc vi phạm: Sheldon Ross [1, Mục 9.6] chỉ ra rằng khi hiện tượng phương sai thay đổi (Heteroscedasticity) xảy ra, mặc dù các ước lượng hệ số \(\hat{\beta}\) vẫn không chệch, nhưng chúng không còn là các ước lượng có phương sai nhỏ nhất.
Hệ quả trực tiếp là:
Các sai số chuẩn (Standard Errors) của các hệ số bị tính toán sai.
Các kiểm định \(t\) và khoảng tin cậy không còn đáng tin cậy, dẫn đến việc kết luận sai lệch về ý nghĩa của các yếu tố rủi ro gây suy thận.
4. Cách khắc phục sơ bộ: Theo GS. Nguyễn Văn Tuấn [2, trang 452], nếu phát hiện phương sai thay đổi mạnh, nhà nghiên cứu có thể cân nhắc sử dụng các phép biến đổi dữ liệu (như lấy logarit biến phụ thuộc) hoặc sử dụng phương pháp bình phương tối thiểu có trọng số (Weighted Least Squares) để điều chỉnh mô hình.
2.4 Các chỉ số đánh giá hiệu năng mô hình
2.4.1 Đánh giá mức độ phù hợp của mô hình hồi quy (Goodness-of-Fit)
Để đánh giá mức độ giải thích của mô hình đối với sự biến thiên của dữ liệu thực tế, ta sử dụng các chỉ số về độ phù hợp (Goodness-of-fit).
Hệ số xác định và hệ số xác định hiệu chỉnh là những chỉ số cơ bản nhất và đã được trình bày trong 2.1.1.3 và 2.1.2.3
Tiêu chuẩn thông tin Akaike (AIC - Akaike Information Criterion) đánh giá mức độ phù hợp dữ liệu của mô hình đã được trình bày ở Định nghĩa 2.13 và Nhận xét 2.17. Ngoài ra [2, trang 184] còn chỉ ra cách thức sử dụng AIC hiệu quả đó là bắt đầu với một mô hình đầy đủ chứa tất cả các biến đã qua bước 1 và 2, sau đó dùng hàm
step()với thuật toán backward để loại dần các biến không đóng góp vào việc giảm chỉ số AIC.Sai số chuẩn của ước lượng (Residual Standard Error - RSE) được nói đến trong Định nghĩa 2.2 là một chỉ số quan trọng, chỉ số này đại diện cho sai số không thể giải thích được của mô hình. RSE càng nhỏ, mô hình dự báo càng chính xác.
2.4.2 Chỉ số cho mô hình phân loại (Classification Metrics)
Trong tiểu luận này, chúng ta nghiên cứu về dự báo y khoa, việc mô hình có “khớp” với dữ liệu (thể hiện qua AIC) là chưa đủ. Ta cần đánh giá khả năng phân loại chính xác giữa nhóm có bệnh và nhóm không bệnh. Theo GS. Nguyễn Văn Tuấn [2], các chỉ số sau đây là tiêu chuẩn để đánh giá hiệu năng mô hình.
1. Ma trận chẩn đoán: Ma trận chận đoán là một bảng thống kê chéo giữa giá trị thực tế quan sát được và giá trị do mô hình dự báo. Từ bảng này, ta xác định được 4 trạng thái cơ bản:
Dương tính thật (True Positive - TP): Người thực sự bị bệnh và mô hình dự báo đúng.
Âm tính thật (True Negative - TN): Người thực sự không bị bệnh và mô hình dự báo đúng.
Dương tính giả (False Positive - FP): Người không bị bệnh nhưng mô hình dự báo là có (Sai lầm loại I).
Âm tính giả (False Negative - FN): Người thực sự bị bệnh nhưng mô hình dự báo là không (Sai lầm loại II).
2. Độ chính xác tổng thể (Accuracy): Là tỉ lệ giữa số ca dự báo đúng trên tổng số quan sát. \[Accuracy = \frac{TP + TN}{TP + TN + FP + FN}\] Tuy nhiên, GS. Nguyễn Văn Tuấn [2] lưu ý rằng nếu tập dữ liệu bị mất cân bằng (ví dụ: số người không bị suy thận nhiều áp đảo số người bị), thì độ chính xác cao không đồng nghĩa với việc mô hình tốt.
3. Độ nhạy (Sensitivity) [2, mục 21.2.8]: Độ nhạy đo lường khả năng của mô hình trong việc nhận diện chính xác những người thực sự bị bệnh. \[Sensitivity = \frac{TP}{TP + FN}\]
- Ý nghĩa y khoa: Trong dự báo bị suy thận, đây là chỉ số quan trọng nhất. Một mô hình có độ nhạy thấp sẽ bỏ sót các bệnh nhân đang gặp nguy hiểm (tăng tỉ lệ FN), dẫn đến hậu quả nghiêm trọng về sức khỏe.
4. Độ đặc hiệu (Specificity) [2, mục 21.2.8] Độ đặc hiệu đo lường khả năng của mô hình trong việc nhận diện chính xác những người không bị bệnh. \[Specificity = \frac{TN}{TN + FP}\]
- Ý nghĩa: Chỉ số này giúp giảm thiểu việc chẩn đoán nhầm cho người khỏe mạnh, tránh gây lo lắng và chi phí điều trị không cần thiết (giảm tỉ lệ FP).
5. Giá trị tiên đoán (Predictive Values):
Giá trị tiên đoán dương (PPV - Precision): Xác suất một người thực sự bị bệnh khi mô hình báo “Có”: \(PPV = TP / (TP + FP)\).
Giá trị tiên đoán âm (NPV): Xác suất một người thực sự không bệnh khi mô hình báo “Không”: \(NPV = TN / (TN + FN)\).
2.4.3 Đường cong ROC và chỉ số AUC
Trong bài toán dự báo suy thận, việc lựa chọn một ngưỡng cắt (threshold) tối ưu để phân loại bệnh nhân là rất quan trọng. Đường cong ROC và chỉ số AUC là các công cụ mạnh mẽ nhất để đánh giá năng lực phân loại tổng quát của mô hình hồi quy Logistic mà không phụ thuộc vào một ngưỡng cắt cố định.
1. Đường cong ROC (Receiver Operating Characteristic)[2, mục 21.2.9] Theo GS. Nguyễn Văn Tuấn, đường cong ROC là một đồ thị biểu diễn mối tương quan giữa:
Trục tung (Y): Độ nhạy (Sensitivity), hay còn gọi là Tỉ lệ dương tính thật (True Positive Rate - TPR).
Trục hoành (X): 1 - Độ đặc hiệu, hay còn gọi là Tỉ lệ dương tính giả (False Positive Rate - FPR).
Đường cong này được tạo ra bằng cách thay đổi liên tục ngưỡng xác suất \(p\) (từ 0 đến 1) của mô hình Logistic để xác định khả năng phân biệt giữa nhóm bị bệnh và nhóm không bị bệnh.
Một mô hình lý tưởng sẽ có đường cong áp sát góc trên bên trái của đồ thị.
Đường chéo nối từ góc (0,0) đến (1,1) đại diện cho một mô hình dự báo ngẫu nhiên (không có giá trị sử dụng).
2. Chỉ số AUC (Area Under the Curve)[2, mục 21.2.9] Chỉ số AUC chính là diện tích nằm dưới đường cong ROC. Về mặt xác suất, AUC đại diện cho khả năng mô hình sẽ xếp hạng một cá thể được chọn ngẫu nhiên từ nhóm “Có bệnh” có xác suất rủi ro cao hơn một cá thể chọn ngẫu nhiên từ nhóm “Không bệnh”.
3. Thang đo đánh giá hiệu năng (theo GS. Nguyễn Văn Tuấn [2]): Để kết luận về độ tin cậy của mô hình dự báo bị bệnh, ta dựa vào giá trị AUC như sau:
AUC = 0.5: Mô hình không có giá trị phân loại (tương đương tung đồng xu).
0.7 \(\le\) AUC < 0.8: Mô hình có khả năng phân loại ở mức Khá.
0.8 \(\le\) AUC < 0.9: Mô hình có khả năng phân loại Tốt.
AUC \(\ge\) 0.9: Mô hình có khả năng phân loại Xuất sắc.
4. Ý nghĩa trong bài toán dự báo suy thận: Việc đạt được chỉ số AUC cao khẳng định rằng các biến độc lập như Tuổi, BMI, nồng độ đường huyết và tình trạng tăng huyết áp đã được mô hình Logistic kết hợp một cách hiệu quả để tạo ra một công cụ dự báo có giá trị lâm sàng cao. Chỉ số này cũng giúp nhà nghiên cứu so sánh các mô hình khác nhau để chọn ra phương án tối ưu nhất cho việc tầm soát bệnh lý.
2.4.4 Kiểm định Likelihood Ratio
Trong hồi quy Logistic, để đánh giá xem việc đưa các biến độc lập (như tuổi, BMI, huyết áp) vào mô hình có thực sự giúp cải thiện khả năng dự báo suy thận so với một mô hình rỗng hay không, ta sử dụng kiểm định Likelihood Ratio (LR).
1. Cơ sở toán học và hàm Log-khả năng đã được trình bày ở mục 2.2.1.
2. Độ hiệu quả của mô hình. Để so sánh hiệu quả của hai mô hình, ta sử dụng giá trị độ lệch dư được trình bày ở Định nghĩa 2.11 và Nhận xét 2.15:
3. Ứng dụng trong R [2, mục 12.4] Trong thực hành,
kiểm định này thường được thực hiện thông qua hàm
anova(model_null, model_full, test="Chisq").
2.5 Nguyên tắc chọn biến đầu vào
- Ba tiêu chuẩn của một mô hình tối ưu Theo GS. Nguyễn Văn Tuấn [2, trang 180, trang 251-253], một mô hình được gọi là tối ưu khi đáp ứng đầy đủ ba tiêu chuẩn sau:
Đơn giản (Simple): Mô hình nên có ít biến số độc lập nhất có thể để việc diễn dịch kết quả trở nên dễ dàng và tránh rườm rà. Đây còn được gọi là Nguyên lý tiết kiệm (Principle of Parsimony) – một mô hình với 3 biến giải thích dữ liệu tốt tương đương mô hình 5 biến thì mô hình 3 biến sẽ được chọn.
Đầy đủ (Adequate): Mô hình phải có khả năng mô tả dữ liệu một cách thỏa đáng, nghĩa là các giá trị do mô hình tiên đoán phải càng gần với giá trị thực tế quan sát được càng tốt.
Có ý nghĩa thực tế: Mô hình phải được yểm trợ bằng lý thuyết chuyên ngành (như ý nghĩa sinh học trong y khoa, ý nghĩa lâm sàng…). Một biến số có thể có ý nghĩa về mặt toán học nhưng nếu không có ý nghĩa thực tế thì mô hình đó chỉ là một trò chơi con số.
- Các công cụ định lượng để lựa chọn biến
Để hiện thực hóa các nguyên tắc trên, GS. Nguyễn Văn Tuấn đề xuất sử dụng các chỉ số thống kê hiện đại thay vì chỉ dựa vào trị số P truyền thống:
Tiêu chuẩn thông tin Akaike (AIC) [2, trang 180]: Đây là thước đo dùng để so sánh các mô hình. AIC giúp dung hòa giữa độ chính xác (tổng bình phương phần dư thấp) và số lượng biến (số thông số ít),. Một mô hình có trị số AIC càng thấp thì được xem là càng tối ưu,.
Bayesian Model Averaging (BMA) [2, trang 186, trang 257]: Đây là phương pháp hiện đại hơn để đánh giá sự bất định trong việc xây dựng mô hình. Thay vì chỉ chọn một mô hình sau cùng, BMA xem xét tất cả các tổ hợp biến khả dĩ và tính toán xác suất một biến thực sự có ảnh hưởng đến kết quả,. Những biến xuất hiện với xác suất cao (nhất quán) trong nhiều mô hình sẽ là những ứng viên tốt nhất để đưa vào mô hình cuối cùng.
- Lưu ý về mối tương quan giữa các biến độc lập
Khi chọn biến, cần xem xét mối tương quan giữa các biến độc lập với nhau. Nếu hai biến độc lập có tương quan rất mạnh (ví dụ như hai chỉ số sinh học cùng phản ánh một quá trình), chúng cung cấp lượng thông tin như nhau. Trong trường hợp này, không cần đưa cả hai vào mô hình vì sẽ làm mô hình phức tạp hơn mà không tăng khả năng tiên đoán đáng kể.
3 GIẢI QUYẾT VẤN ĐỀ
3.1 Nạp dữ liệu và quan sát tổng quan
Đầu tiên, ta cài dặt môi trường làm việc với một số gói thông dụng và thư viện tương ứng của chúng.
## package 'tidyverse' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\win\AppData\Local\Temp\RtmpCoaQ1i\downloaded_packages## package 'janitor' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\win\AppData\Local\Temp\RtmpCoaQ1i\downloaded_packages## package 'lubridate' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\win\AppData\Local\Temp\RtmpCoaQ1i\downloaded_packages## package 'ggplot2' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\win\AppData\Local\Temp\RtmpCoaQ1i\downloaded_packages## package 'caret' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\win\AppData\Local\Temp\RtmpCoaQ1i\downloaded_packages## package 'e1071' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\win\AppData\Local\Temp\RtmpCoaQ1i\downloaded_packageslibrary(tidyverse) # Thư viện tổng hợp phân tích dữ liệu
library(janitor) # Thư viện làm sạch dữ liệu
library(lubridate) # Thư viên ngày giờ
library(dplyr) # Thư viện thao tác dữ liệu
library(ggplot2)Tiếp theo, ta quan sát sơ bộ dữ liệu.
# Ta nạp dữ liệu vào (kn là viết tắt của kidney)
kn <- read.csv("D:/R/Tieu luan xsnc/kidney_disease_dataset.csv")
# Ta xem sơ qua dữ liệu
View(kn)
# Ta
str(kn)## 'data.frame': 400 obs. of 26 variables:
## $ id : int 0 1 2 3 4 5 6 7 8 9 ...
## $ age : num 48 7 62 48 51 60 68 24 52 53 ...
## $ bp : num 80 50 80 70 80 90 70 NA 100 90 ...
## $ sg : num 1.02 1.02 1.01 1 1.01 ...
## $ al : num 1 4 2 4 2 3 0 2 3 2 ...
## $ su : num 0 0 3 0 0 0 0 4 0 0 ...
## $ rbc : chr "" "" "normal" "normal" ...
## $ pc : chr "normal" "normal" "normal" "abnormal" ...
## $ pcc : chr "notpresent" "notpresent" "notpresent" "present" ...
## $ ba : chr "notpresent" "notpresent" "notpresent" "notpresent" ...
## $ bgr : num 121 NA 423 117 106 74 100 410 138 70 ...
## $ bu : num 36 18 53 56 26 25 54 31 60 107 ...
## $ sc : num 1.2 0.8 1.8 3.8 1.4 1.1 24 1.1 1.9 7.2 ...
## $ sod : num NA NA NA 111 NA 142 104 NA NA 114 ...
## $ pot : num NA NA NA 2.5 NA 3.2 4 NA NA 3.7 ...
## $ hemo : num 15.4 11.3 9.6 11.2 11.6 12.2 12.4 12.4 10.8 9.5 ...
## $ pcv : chr "44" "38" "31" "32" ...
## $ wc : chr "7800" "6000" "7500" "6700" ...
## $ rc : chr "5.2" "" "" "3.9" ...
## $ htn : chr "yes" "no" "no" "yes" ...
## $ dm : chr "yes" "no" "yes" "no" ...
## $ cad : chr "no" "no" "no" "no" ...
## $ appet : chr "good" "good" "poor" "poor" ...
## $ pe : chr "no" "no" "no" "yes" ...
## $ ane : chr "no" "no" "yes" "yes" ...
## $ classification: chr "ckd" "ckd" "ckd" "ckd" ...Giải thích tên biến thuộc tính (tên cột)
| Tên biến gốc | Tên rút gọn | Ý nghĩa | Đơn vị / Loại dữ liệu |
|---|---|---|---|
| id | id | Mã số bệnh nhân | Số nguyên (Bỏ qua khi PT) |
| age | age | Tuổi của bệnh nhân | Năm (Số thực) |
| bp | bp | Huyết áp tâm thu | mm/Hg (Số thực) |
| sg | sg | Tỷ trọng nước tiểu | 1.005, 1.010, … (Số thực) |
| al | al | Albumin niệu | 0, 1, 2, 3, 4, 5 (Định mức) |
| su | su | Đường niệu | 0, 1, 2, 3, 4, 5 (Định mức) |
| rbc | rbc | Hồng cầu trong nước tiểu | normal / abnormal |
| pc | pc | Tế bào mủ trong nước tiểu | normal / abnormal |
| pcc | pcc | Cụm tế bào mủ | present / notpresent |
| ba | ba | Vi khuẩn | present / notpresent |
| bgr | glucose | Đường huyết ngẫu nhiên | mg/dl (Số thực) |
| bu | urea | Ure trong máu | mg/dl (Số thực) |
| sc | creatinine | Creatinine huyết thanh | mg/dl (Số thực) |
| sod | sodium | Nồng độ Natri | mEq/L (Số thực) |
| pot | potassium | Nồng độ Kali | mEq/L (Số thực) |
| hemo | hemo | Nồng độ Huyết sắc tố | gms (Số thực) |
| pcv | pcv | Dung tích hồng cầu | % (Số thực) |
| wc | wc | Số lượng bạch cầu | cells/cumm (Số nguyên) |
| rc | rc | Số lượng hồng cầu | millions/cumm (Số thực) |
| htn | htn | Tiền sử tăng huyết áp | yes / no |
| dm | dm | Tiền sử đái tháo đường | yes / no |
| cad | cad | Bệnh mạch vành | yes / no |
| appet | appetite | Sự thèm ăn | good / poor |
| pe | edema | Phù chân | yes / no |
| ane | anemia | Tình trạng thiếu máu | yes / no |
| classification | target | Kết luận bệnh (CKD) | ckd / notckd (Biến mục tiêu) |
| (Tính thêm) | egfr | Mức lọc cầu thận | ml/min/1.73m2 (Số thực) |
Đánh giá sơ bộ:
Quy mô mẫu: Tập dữ liệu bao gồm 400 quan sát bệnh nhân thực tế với 26 biến thuộc tính lâm sàng.
Các biến số cho Hồi quy tuyến tính: Tập trung vào các chỉ số sinh hóa
quan trọng như: bp (Huyết áp), sc (Creatinine
huyết thanh), hemo (Nồng độ huyết sắc tố). Đặc biệt, biến
phụ thuộc chính là eGFR (Mức lọc cầu thận) sẽ được tính
toán từ age và sc theo công thức y khoa để đánh giá chính
xác mức độ suy giảm chức năng thận.
Biến mục tiêu cho Hồi quy Logistic: Biến classification đóng vai trò
là biến mục tiêu, phân loại bệnh nhân thành hai nhóm: ckd
(Mắc bệnh thận mãn tính) và notckd (Không mắc bệnh). Biến
này sẽ được mã hóa về dạng nhị phân (0 và 1) để phù hợp với mô hình.
Các biến định tính: Các yếu tố nguy cơ quan trọng như
htn (Tăng huyết áp), dm (Đái tháo đường),
cad (Bệnh mạch vành) và các dấu hiệu lâm sàng như
ane (Thiếu máu), pe (Phù chân) sẽ được chuyển
đổi sang dạng số để đưa vào mô hình dự báo.
- Đặc điểm dữ liệu khuyết thiếu: Điểm cần lưu ý đặc biệt ở mẫu số liệu này là sự hiện diện của các giá trị khuyết thiếu ở nhiều thuộc tính.
3.2 Làm sạch và xử lý ban đầu dữ liệu
Như đã nhận xét ở trên, ta thấy có 26 biến thuộc tính (cột), ta sẽ lọc bỏ các biến không quan trọng sau dựa trên cơ sở y học và dựa trên nguyên tắc giảm nhiễu trong phân tích dữ liệu.
3.2.1 Xử lý rút gọn
Nhóm 1: Loại bỏ do dư thừa và đa cộng tuyến
+) pcv (Dung tích hồng cầu) và rc (Số lượng
hồng cầu)
Lý do y khoa: Trong y học về thận, thiếu máu là hệ quả của việc thiếu hormone Erythropoietin (do thận sản xuất). Cả
pcv,rcvàhemođều cùng đo lường một tình trạng sinh lý này.Lý do trong phân tích dữ liệu:
hemođược coi là tiêu chuẩn vàng lâm sàng để đánh giá thiếu máu. Việc giữ lại cả 3 biến sẽ gây hiện tượng đa cộng tuyến (các biến giải thích trùng lặp nhau), làm sai lệch hệ số hồi quy.rccũng là biến có tỷ lệ khuyết thiếu rất cao (~32%) nên đáng để loại bỏ khỏi phân tích của chúng ta.
Nhóm 2: Loại bỏ do dữ liệu khuyết thiếu quá lớn (>30%)
+) rbc (Hồng cầu trong nước tiểu)
Lý do y khoa: Sự hiện diện của hồng cầu trong nước tiểu có độ nhạy thấp và không đặc hiệu cho suy thận mạn. Nó có thể do sỏi, nhiễm trùng nhất thời hoặc thậm chí là chu kỳ sinh lý ở nữ giới.
Lý do trong phân tích dữ liệu: Với tỷ lệ khuyết thiếu lên đến hơn 38%, việc điền khuyết cho biến này sẽ tạo ra sai số lớn, làm giảm độ tin cậy của mô hình.
+) wc (Số lượng bạch cầu)
Lý do y khoa: Bạch cầu tăng phản ánh tình trạng viêm hoặc nhiễm trùng toàn thân, không phản ánh trực tiếp mức độ suy giảm chức năng của đơn vị thận (Nephron).
Lý do trong phân tích dữ liệu: Biến này khuyết thiếu khoảng 26% và thường dao động mạnh do các tác nhân ngoài thận, gây nhiễu cho dự báo dài hạn.
Nhóm 3: Loại bỏ do tính chủ quan và sai số quan sát
+) appet (Sự thèm ăn) và pe (Phù chân)
Lý do y khoa: Đây là các triệu chứng lâm sàng mang tính cảm tính. “Chán ăn” hay “Phù” phụ thuộc vào cảm nhận của bệnh nhân và sự đánh giá bằng mắt của bác sĩ.
Lỹ do trong phân tích dữ liệu: Trong mô hình toán học nâng cao, ta ưu tiên các chỉ số định lượng (số liệu máy đo) như creatinine hay albumin hơn là các biến định tính chủ quan để đảm bảo tính khách quan của kết quả.
+) ane (Tình trạng thiếu máu - dạng Yes/No)
Lý do y khoa: Đây là biến định danh dựa trên quan sát lâm sàng (nhìn niêm mạc mắt nhợt…).
Lý do trong phân tích dữ liệu: Chúng ta đã giữ lại biến số thực
hemo(Huyết sắc tố). Biếnhemocung cấp thông tin chính xác và chi tiết hơn nhiều so với việc chỉ kết luận đơn thuần là “có” hay “không” thiếu máu.
Nhóm 4: Loại bỏ do biến động thấp (Low Variance/Acute phase)
+) pcc (Cụm tế bào mủ) và ba (Vi khuẩn)
Lý do y khoa: Hai biến này thường chỉ xuất hiện trong các đợt nhiễm trùng đường tiết niệu cấp tính (UTI). Đối với bệnh thận mạn tính (CKD), phần lớn bệnh nhân sẽ có kết quả là “notpresent”.
Lý do trong phân tích dữ liệu: Do phần lớn dữ liệu tập trung ở một giá trị (notpresent), các biến này không giúp ích nhiều trong việc phân loại mức độ bệnh thận, dễ làm mô hình bị lệch.
Nhóm 5: Biến định danh không mang giá trị sinh học
+) id (Mã bệnh nhân)
- Lý do: Đây chỉ là số thứ tự định danh hành chính, không có liên hệ nhân quả với tình trạng bệnh lý.
Sau đây ta phân loại và định nghĩa các biến thuộc tính theo các nhóm chức năng dưới đây. Các tên rút gọn sẽ được sử dụng thống nhất trong suốt quá trình xử lý mã nguồn R.
3.2.2 Xử lý tên biến
Tổng kết lại ta giữ lại các biến thuộc tính sau
Nhân khẩu học: age
Sinh hiệu: bp
Nước tiểu: sg, al, su, pc
Máu: bgr, bu, sc, sod, pot, hemo
Bệnh lý: htn, dm, cad
Mục tiêu: target (và egfr cho hồi quy tuyến tính)
Các biến này có tên đã được rút ngắn sắn nên ta sẽ không thay đổi gì thêm. Thông tin chi tiết về các biens thuộc tính này được trình bày ở mục “Nạp dữ liệu và quan sát tổng quan”
3.2.3 Xử lý khuyết thiếu
Trong xử lý dữ liệu y tế, việc giải quyết giá trị khuyết thiếu (NA) là cực kỳ quan trọng vì nó ảnh hưởng trực tiếp đến độ tin cậy của mô hình hồi quy.
Kiểm tra khuyết thiếu
Ta kiểm tra sự khuyết thiếu theo các biến thuộc tính
library(tidyverse)
# 1. Đọc dữ liệu và chọn lọc 16 biến chính
df <- read.csv("kidney_disease_dataset.csv", stringsAsFactors = FALSE) %>%
select(age, bp, sg, al, su, pc, bgr, bu, sc, sod, pot, hemo, htn, dm, cad, classification)
# 2. Quan sát cấu trúc dữ liệu (Biến nào là số, biến nào là chữ)
glimpse(df)## Rows: 400
## Columns: 16
## $ age <dbl> 48, 7, 62, 48, 51, 60, 68, 24, 52, 53, 50, 63, 68, 68, …
## $ bp <dbl> 80, 50, 80, 70, 80, 90, 70, NA, 100, 90, 60, 70, 70, 70…
## $ sg <dbl> 1.020, 1.020, 1.010, 1.005, 1.010, 1.015, 1.010, 1.015,…
## $ al <dbl> 1, 4, 2, 4, 2, 3, 0, 2, 3, 2, 2, 3, 3, NA, 3, 3, 2, NA,…
## $ su <dbl> 0, 0, 3, 0, 0, 0, 0, 4, 0, 0, 4, 0, 1, NA, 2, 0, 0, NA,…
## $ pc <chr> "normal", "normal", "normal", "abnormal", "normal", "",…
## $ bgr <dbl> 121, NA, 423, 117, 106, 74, 100, 410, 138, 70, 490, 380…
## $ bu <dbl> 36, 18, 53, 56, 26, 25, 54, 31, 60, 107, 55, 60, 72, 86…
## $ sc <dbl> 1.2, 0.8, 1.8, 3.8, 1.4, 1.1, 24.0, 1.1, 1.9, 7.2, 4.0,…
## $ sod <dbl> NA, NA, NA, 111.0, NA, 142.0, 104.0, NA, NA, 114.0, NA,…
## $ pot <dbl> NA, NA, NA, 2.5, NA, 3.2, 4.0, NA, NA, 3.7, NA, 4.2, 5.…
## $ hemo <dbl> 15.4, 11.3, 9.6, 11.2, 11.6, 12.2, 12.4, 12.4, 10.8, 9.…
## $ htn <chr> "yes", "no", "no", "yes", "no", "yes", "no", "no", "yes…
## $ dm <chr> "yes", "no", "yes", "no", "no", "yes", "no", "yes", "ye…
## $ cad <chr> "no", "no", "no", "no", "no", "no", "no", "no", "no", "…
## $ classification <chr> "ckd", "ckd", "ckd", "ckd", "ckd", "ckd", "ckd", "ckd",…# 3. Kiểm tra số lượng giá trị thiếu (NA) cho từng biến
missing_data <- colSums(is.na(df))
print("So luong gia tri khuyet thieu (NA) moi bien:")## [1] "So luong gia tri khuyet thieu (NA) moi bien:"## age bp sg al su
## 9 12 47 46 49
## pc bgr bu sc sod
## 0 44 19 17 87
## pot hemo htn dm cad
## 88 52 0 0 0
## classification
## 0## age bp sg al
## Min. : 2.00 Min. : 50.00 Min. :1.005 Min. :0.000
## 1st Qu.:42.00 1st Qu.: 70.00 1st Qu.:1.010 1st Qu.:0.000
## Median :55.00 Median : 80.00 Median :1.020 Median :0.000
## Mean :51.48 Mean : 76.47 Mean :1.017 Mean :1.017
## 3rd Qu.:64.50 3rd Qu.: 80.00 3rd Qu.:1.020 3rd Qu.:2.000
## Max. :90.00 Max. :180.00 Max. :1.025 Max. :5.000
## NA's :9 NA's :12 NA's :47 NA's :46
## su pc bgr bu
## Min. :0.0000 Length:400 Min. : 22 Min. : 1.50
## 1st Qu.:0.0000 Class :character 1st Qu.: 99 1st Qu.: 27.00
## Median :0.0000 Mode :character Median :121 Median : 42.00
## Mean :0.4501 Mean :148 Mean : 57.43
## 3rd Qu.:0.0000 3rd Qu.:163 3rd Qu.: 66.00
## Max. :5.0000 Max. :490 Max. :391.00
## NA's :49 NA's :44 NA's :19
## sc sod pot hemo
## Min. : 0.400 Min. : 4.5 Min. : 2.500 Min. : 3.10
## 1st Qu.: 0.900 1st Qu.:135.0 1st Qu.: 3.800 1st Qu.:10.30
## Median : 1.300 Median :138.0 Median : 4.400 Median :12.65
## Mean : 3.072 Mean :137.5 Mean : 4.627 Mean :12.53
## 3rd Qu.: 2.800 3rd Qu.:142.0 3rd Qu.: 4.900 3rd Qu.:15.00
## Max. :76.000 Max. :163.0 Max. :47.000 Max. :17.80
## NA's :17 NA's :87 NA's :88 NA's :52
## htn dm cad classification
## Length:400 Length:400 Length:400 Length:400
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## Có 3 vấn đề chính chúng ta cần xử lý trước khi tính toán:
Biến số bị thiếu (NA): Đặc biệt là Natri (
sod), Kali (pot) và Huyết sắc tố (hemo) thiếu khá nhiều (~50-88 mẫu).Biến chữ: Các biến như
htn,dm,classificationđang ở dạng chữ “yes”/“no”. Máy tính không cộng trừ nhân chia chữ được, nên phải đổi về 0 và 1.Lỗi nhập liệu ẩn
Xử lý
# 1. Điền giá trị thiếu (NA) bằng Trung vị (Median) cho các cột. Cách này giúp giữ nguyên 400 mẫu mà không làm lệch dữ liệu
df_clean <- df %>%
mutate(across(where(is.numeric), ~ifelse(is.na(.), median(., na.rm = TRUE), .)))
# 2. Làm sạch chữ và chuyển về dạng số (0/1)
df_clean <- df_clean %>%
mutate(
# Xóa khoảng trắng thừa
across(where(is.character), str_trim),
# Chuyển yes/no thành 1/0
htn = ifelse(htn == "yes", 1, 0),
dm = ifelse(dm == "yes", 1, 0),
cad = ifelse(cad == "yes", 1, 0),
# Chuyển pc (tế bào mủ): abnormal = 1, normal = 0, ở đây lưu ý là các ô trống trong biến `pc` không coi là N/A mà là chuỗi trống nên kết quả trong cơ sở dữ liệu mới sẽ là số 0
pc = ifelse(pc == "abnormal", 1, 0)
)
# Xem lại kết quả sau khi làm sạch
head(df_clean)## age bp sg al su pc bgr bu sc sod pot hemo htn dm cad classification
## 1 48 80 1.020 1 0 0 121 36 1.2 138 4.4 15.4 1 1 0 ckd
## 2 7 50 1.020 4 0 0 121 18 0.8 138 4.4 11.3 0 0 0 ckd
## 3 62 80 1.010 2 3 0 423 53 1.8 138 4.4 9.6 0 1 0 ckd
## 4 48 70 1.005 4 0 1 117 56 3.8 111 2.5 11.2 1 0 0 ckd
## 5 51 80 1.010 2 0 0 106 26 1.4 138 4.4 11.6 0 0 0 ckd
## 6 60 90 1.015 3 0 0 74 25 1.1 142 3.2 12.2 1 1 0 ckd## age bp sg al su
## Min. : 2.00 Min. : 50.00 Min. :1.005 Min. :0.0 Min. :0.000
## 1st Qu.:42.00 1st Qu.: 70.00 1st Qu.:1.015 1st Qu.:0.0 1st Qu.:0.000
## Median :55.00 Median : 80.00 Median :1.020 Median :0.0 Median :0.000
## Mean :51.56 Mean : 76.58 Mean :1.018 Mean :0.9 Mean :0.395
## 3rd Qu.:64.00 3rd Qu.: 80.00 3rd Qu.:1.020 3rd Qu.:2.0 3rd Qu.:0.000
## Max. :90.00 Max. :180.00 Max. :1.025 Max. :5.0 Max. :5.000
## pc bgr bu sc
## Min. :0.00 Min. : 22.0 Min. : 1.50 Min. : 0.400
## 1st Qu.:0.00 1st Qu.:101.0 1st Qu.: 27.00 1st Qu.: 0.900
## Median :0.00 Median :121.0 Median : 42.00 Median : 1.300
## Mean :0.19 Mean :145.1 Mean : 56.69 Mean : 2.997
## 3rd Qu.:0.00 3rd Qu.:150.0 3rd Qu.: 61.75 3rd Qu.: 2.725
## Max. :1.00 Max. :490.0 Max. :391.00 Max. :76.000
## sod pot hemo htn
## Min. : 4.5 Min. : 2.500 Min. : 3.10 Min. :0.0000
## 1st Qu.:135.0 1st Qu.: 4.000 1st Qu.:10.88 1st Qu.:0.0000
## Median :138.0 Median : 4.400 Median :12.65 Median :0.0000
## Mean :137.6 Mean : 4.577 Mean :12.54 Mean :0.3675
## 3rd Qu.:141.0 3rd Qu.: 4.800 3rd Qu.:14.62 3rd Qu.:1.0000
## Max. :163.0 Max. :47.000 Max. :17.80 Max. :1.0000
## dm cad classification
## Min. :0.0000 Min. :0.000 Length:400
## 1st Qu.:0.0000 1st Qu.:0.000 Class :character
## Median :0.0000 Median :0.000 Mode :character
## Mean :0.3425 Mean :0.085
## 3rd Qu.:1.0000 3rd Qu.:0.000
## Max. :1.0000 Max. :1.0003.2.4 Chuẩn hóa biến mục tiêu (biến đầu ra)
Mô hình hồi quy tuyến tính
Đối với hồi quy tuyến tính ta có biến mục tiêu chính là chỉ số eGPR. Chỉ số này chưa có trong cơ sở dữ liệu nên ta sẽ tính toán thêm.
# Tính eGFR (Biến phụ thuộc cho hồi quy tuyến tính)
# Công thức MDRD: 186 * (Creatinine^-1.154) * (Tuổi^-0.203)
df_clean <- df_clean %>%
mutate(egfr = 186 * (sc^-1.154) * (age^-0.203))Mô hình hồi quy tuyến tính Logistic
Vì máy tính (đặc biệt là các hàm hồi quy trong R như lm
hay glm) chỉ làm việc với các con số, nên ta thực hiện
chuyển đổi các cột thuộc tính dạng chữ (“ckd”, “notckd”) thành các giá
trị số.
3.3 Xây dựng mô hình hồi quy
Đê xây dựng mô hình hồi quy, trước hết, ta chia cơ sở dữ liệu của chúng ta làm 2 phần:
Phần huấn luyện mô hình (train_data) gồm 80% quan sát của cơ sở dữ liệu.
Phần kiểm tra mô hình (test_data) gồm 20% quan sát của cơ sở dữ liệu.
# 1. Thiết lập seed để kết quả chia ngẫu nhiên có thể lặp lại được
set.seed(1000)
# 2. Xác định số lượng dòng cho tập huấn luyện (80% dữ liệu)
n <- nrow(df_clean)
sample_size <- floor(0.8 * n)
# 3. Lấy ngẫu nhiên chỉ số các dòng
train_indices <- sample(seq_len(n), size = sample_size)
# 4. Chia dữ liệu thành 2 phần
train_data <- df_clean[train_indices, ]
test_data <- df_clean[-train_indices, ]
# 5. Kiểm tra lại kích thước của 2 tập dữ liệu
cat("Số lượng quan sát tập Huấn luyện (Train):", nrow(train_data), "\n")## Số lượng quan sát tập Huấn luyện (Train): 320## Số lượng quan sát tập Kiểm tra (Test): 80Tiếp theo, ta sẽ sử dụng các nguyên tắc được tìm hiểu trong mục 2.5. Giữa mô hình hồi quy tuyến tính và mô hình hồi quy logistic có sự khác biệt cơ bản về biến đầu ra. Do đó quy trình và cách thức chọn biến đầu vào của hai mô hình sẽ có sự khác biết. Trước hết, ta lựa chọn biến đầu vào cho mô hình hồi quy tuyến tính.
3.3.1 Biến đầu vào cho mô hình hồi quy tuyến tính
1. Phân tích tương quan
Bước thực hiện này là quan trọng để hiểu được cấu trúc của cơ sở dữ liệu. Từ đó lựa chọn ra những biến thuộc tính để xây dựng mô hình hồi quy sao cho không bị trùng lặp với nhau và lựa chọn được ít biến thuộc tính nhất cho mô hình, đảm bảo nguyên tắc đơn giản. Bước làm này còn rất quan trọng khi làm việc với cơ sở dữ liệu cho nhiều biến thuộc tính như này. (xem thêm Nhận xét 2.7)
library(ggcorrplot)
library(dplyr)
library(stringr)
# 2. Tính ma trận tương quan
corr_matrix <- cor(df_clean)
# 3. Vẽ Heatmap
ggcorrplot(corr_matrix,
hc.order = TRUE,
type = "lower",
lab = TRUE,
lab_size = 3,
colors = c("#D73027", "white", "#4575B4"),
title = "Ma trận tương quan - Dữ liệu thực tế UCI",
ggtheme = theme_minimal())
Ta thực hiện lọc các biến thuộc tính phù hợp để làm đầu vào dựa trên nguyên tắc được trình bày ở mục 2.5.
Dựa vào heatmap, ta thấy các biến có tương quan mạnh với
egfrlàbu,sc,htn,hemo,dm,age,sg(không tính biếntargetvì nếu bị bệnh thì chắc chắnegfrsẽ cao).Tiếp theo ta loại
dmvìdm(tiểu đường) vàhtn(tăng huyết áp) có mối tương quan cao, trong y khoa thì 2 bệnh này hầu như luôn đi cùng nhau nên ta chọn 1 trong 2. Ta chọnhtnvì nó có tương quan vớiegfrcao hơn.Ta loại
buvìbu(ure máu) vàsc(creatinine) có tương quan mạnh với nhau. Mặc dùbucó tuong quan mạnh vớiegfrhơn nhưng ta vẫn chọnbuvì trong y học Creatinine là là chỉ số đặc trung hơn gắn liền với bệnh thận, còn ure máu có thể phát sinh từ nhiều loại bệnh khác nhau nữa.Ta không loại cả hai biến
hemo(hemoglobin) vàhtn(tăng huyết áp) mặc dù chúng có tương quan nghịch với nhau vì chúng là 2 khía cạnh khác nhau trong y khoa (một bên là thiếu máu, một bên là tăng huyết áp), sự có mặt cả 2 chỉ số có thể là một gọi ý quan trọng để nhận diện bệnh thận với các bệnh khác.
Để chắc chắn hơn ta sử dụng tối ưu hóa bằng chỉ số AIC được trình bày
ở mục 2.4.1. Ta bắt đầu với một mô hình đầy đủ chứa tất cả 14 biến (trừ
egfr và target), sau đó dùng hàm
step() với thuật toán backward để loại dần các biến không
đóng góp vào việc giảm chỉ số AIC.
# Xây dựng mô hình đầy đủ
full_model <- lm(egfr ~ sg + hemo + pot + bu + sc + su + bgr + htn + dm + age + cad + bp + al + pc, data = train_data)
# Tự động lựa chọn biến tối ưu
optimal_lm <- step(full_model, direction = "backward", trace = 0)
summary(optimal_lm)##
## Call:
## lm(formula = egfr ~ sg + hemo + bu + sc + htn + age, data = train_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -102.044 -24.323 -6.249 16.715 139.174
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1352.5036 459.6863 -2.942 0.00350 **
## sg 1403.3011 455.8593 3.078 0.00227 **
## hemo 4.8952 1.1427 4.284 2.45e-05 ***
## bu -0.1861 0.0601 -3.096 0.00214 **
## sc -1.1304 0.4962 -2.278 0.02339 *
## htn -19.3745 5.9524 -3.255 0.00126 **
## age -0.8732 0.1403 -6.222 1.57e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37.88 on 313 degrees of freedom
## Multiple R-squared: 0.5459, Adjusted R-squared: 0.5372
## F-statistic: 62.72 on 6 and 313 DF, p-value: < 2.2e-16Ta thấy giữa chọn thủ công và chọn bằng step() có sự
tương đồng nhau. Đến đây ta chọn được mô hình hồi quy tuyến tính đa biến
cho egfr sẽ có biến đầu vào là sg,
hemo, bu, sc, htn,
age. Vậy ta có thể thấy việc sử dụng kiến thức y khoa vào
chọn biến đầu vào có thể xử lý nhanh bằng công cụ R.
Để chắc chắn hơn nữa, ta dùng chỉ số VIF để đánh giá. Đây là chỉ số đặc trung của mô hình hồi quy tuyến tính đa biến, tránh vấn đề đa cộng tuyến khi có nhiều biến.
## package 'car' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\win\AppData\Local\Temp\RtmpCoaQ1i\downloaded_packages## sg hemo bu sc htn age
## 1.308469 2.116659 2.045404 1.657281 1.845321 1.2422593.3.2 Biến đầu vào cho mô hình logistic
Với mô hình này ta trực tiếp sử dụng công cụ R để thực hiện chọn biến đầu vào nhanh.
# 1. Xây dựng mô hình Logistic đầy đủ (Full Model) (bỏ biến egfr)
full_logic <- glm(target ~ sg + hemo + pot + bu + sc + su + bgr + htn + dm + age + cad + bp + al + pc, data = train_data, family = "binomial")
# 2. Sử dụng Stepwise Backward để tối ưu AIC
optimal_logic <- step(full_logic, direction = "backward", trace = 0)
# 3. Xem các biến được giữ lại
summary(optimal_logic)##
## Call:
## glm(formula = target ~ sg + hemo + sc + bgr + al, family = "binomial",
## data = train_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.416e+04 1.810e+06 0.013 0.989
## sg -2.264e+04 1.701e+06 -0.013 0.989
## hemo -1.230e+02 8.016e+03 -0.015 0.988
## sc 1.477e+02 8.378e+03 0.018 0.986
## bgr 2.697e+00 1.626e+02 0.017 0.987
## al 7.761e+01 6.920e+03 0.011 0.991
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 4.2237e+02 on 319 degrees of freedom
## Residual deviance: 1.0724e-06 on 314 degrees of freedom
## AIC: 12
##
## Number of Fisher Scoring iterations: 25Vậy là theo chỉ số AIC, các biến đầu vào được chọn là
sg, hemo, sc, bgr,
al. Tuy nhiên vấn đề gặp phải ở đây là hiện tượng khớp hoàn
toàn. Nghĩa là có một chỉ số nào đó trong đây dự báo đúng 100% bị bệnh
hay không.
# Kiểm tra xem có giá trị nào của Albumin mà 100% thuộc 1 nhóm không
table(train_data$al, train_data$target)##
## 0 1
## 0 119 80
## 1 0 33
## 2 0 32
## 3 0 34
## 4 0 21
## 5 0 1##
## 0 1
## 1.005 0 4
## 1.01 0 68
## 1.015 0 60
## 1.02 68 59
## 1.025 51 10##
## 0 1
## 3.1 0 1
## 5.5 0 1
## 5.6 0 1
## 5.8 0 1
## 6.1 0 1
## 6.2 0 1
## 6.3 0 1
## 6.8 0 1
## 7.1 0 1
## 7.3 0 1
## 7.5 0 1
## 7.6 0 1
## 7.7 0 1
## 7.9 0 3
## 8 0 2
## 8.1 0 3
## 8.2 0 1
## 8.3 0 3
## 8.4 0 1
## 8.5 0 1
## 8.6 0 2
## 8.7 0 2
## 8.8 0 1
## 9 0 1
## 9.1 0 4
## 9.2 0 1
## 9.3 0 2
## 9.4 0 2
## 9.5 0 2
## 9.6 0 3
## 9.7 0 3
## 9.8 0 6
## 9.9 0 3
## 10 0 3
## 10.1 0 2
## 10.2 0 2
## 10.3 0 5
## 10.4 0 1
## 10.5 0 2
## 10.6 0 1
## 10.7 0 1
## 10.8 0 4
## 10.9 0 7
## 11 0 3
## 11.1 0 6
## 11.2 0 3
## 11.3 0 6
## 11.4 0 2
## 11.5 0 2
## 11.6 0 3
## 11.7 0 2
## 11.8 0 2
## 11.9 0 2
## 12 0 6
## 12.1 0 2
## 12.2 0 3
## 12.3 0 2
## 12.4 0 4
## 12.5 0 4
## 12.6 0 3
## 12.65 5 35
## 12.7 0 3
## 12.8 0 1
## 12.9 0 1
## 13 3 4
## 13.1 0 1
## 13.2 2 1
## 13.4 1 1
## 13.5 3 1
## 13.6 4 1
## 13.7 3 0
## 13.8 3 1
## 13.9 4 1
## 14 3 1
## 14.1 3 0
## 14.2 2 0
## 14.3 2 2
## 14.4 3 0
## 14.5 3 0
## 14.6 2 0
## 14.7 2 0
## 14.8 2 0
## 14.9 2 0
## 15 9 2
## 15.1 2 0
## 15.2 2 1
## 15.3 3 0
## 15.4 4 0
## 15.5 4 0
## 15.6 1 1
## 15.7 3 0
## 15.8 2 0
## 15.9 3 0
## 16 2 0
## 16.1 3 1
## 16.2 3 0
## 16.3 2 0
## 16.4 3 0
## 16.5 2 0
## 16.6 2 0
## 16.7 1 0
## 16.8 1 0
## 16.9 2 0
## 17 3 0
## 17.1 1 0
## 17.2 2 0
## 17.3 1 0
## 17.4 1 0
## 17.5 1 0
## 17.6 1 0
## 17.7 1 0
## 17.8 2 0##
## 0 1
## 0.4 1 0
## 0.5 19 1
## 0.6 13 2
## 0.7 15 4
## 0.8 7 6
## 0.9 13 6
## 1 11 6
## 1.1 12 6
## 1.2 23 5
## 1.3 5 14
## 1.4 0 7
## 1.5 0 7
## 1.6 0 6
## 1.7 0 7
## 1.8 0 7
## 1.9 0 5
## 2 0 5
## 2.1 0 5
## 2.2 0 8
## 2.3 0 1
## 2.4 0 3
## 2.5 0 5
## 2.6 0 1
## 2.7 0 3
## 2.8 0 6
## 2.9 0 2
## 3 0 2
## 3.2 0 4
## 3.25 0 1
## 3.3 0 4
## 3.4 0 2
## 3.6 0 1
## 3.8 0 1
## 3.9 0 4
## 4 0 1
## 4.1 0 2
## 4.3 0 1
## 4.4 0 2
## 4.6 0 2
## 5.2 0 2
## 5.3 0 2
## 5.6 0 2
## 5.9 0 1
## 6 0 2
## 6.1 0 1
## 6.3 0 2
## 6.5 0 2
## 6.7 0 1
## 6.8 0 1
## 7.1 0 1
## 7.2 0 2
## 7.3 0 4
## 7.5 0 1
## 7.7 0 1
## 9.2 0 1
## 9.3 0 1
## 9.6 0 1
## 9.7 0 1
## 10.2 0 1
## 11.8 0 1
## 12 0 1
## 12.8 0 1
## 13 0 1
## 13.3 0 1
## 13.4 0 1
## 13.5 0 1
## 13.8 0 1
## 14.2 0 1
## 15 0 1
## 16.4 0 1
## 16.9 0 1
## 18.1 0 1
## 24 0 1
## 32 0 1
## 76 0 1##
## 0 1
## 22 0 1
## 70 2 2
## 74 1 2
## 76 0 3
## 78 1 2
## 79 2 0
## 80 1 1
## 81 3 0
## 82 3 0
## 83 2 0
## 84 0 1
## 85 1 0
## 86 1 1
## 87 1 0
## 88 2 1
## 89 1 0
## 90 0 1
## 91 2 1
## 92 1 3
## 93 2 4
## 94 2 2
## 95 2 2
## 96 2 0
## 97 3 0
## 98 1 1
## 99 6 3
## 100 5 3
## 101 0 2
## 102 2 1
## 103 0 3
## 104 2 3
## 105 2 1
## 106 2 1
## 107 3 4
## 108 0 2
## 109 3 2
## 110 1 1
## 111 1 1
## 112 2 1
## 113 2 1
## 114 2 2
## 116 1 0
## 117 3 3
## 118 1 1
## 119 2 0
## 120 1 1
## 121 7 29
## 122 2 1
## 123 2 3
## 124 2 2
## 125 4 0
## 127 2 1
## 128 1 1
## 129 1 2
## 130 4 1
## 131 4 1
## 132 4 1
## 133 3 0
## 134 1 0
## 137 2 0
## 138 0 1
## 139 0 3
## 140 3 2
## 141 0 2
## 143 0 1
## 144 0 2
## 146 0 1
## 148 0 1
## 150 0 2
## 153 0 2
## 156 0 1
## 157 0 1
## 158 0 2
## 159 0 1
## 160 0 1
## 162 0 1
## 163 0 1
## 165 0 1
## 169 0 1
## 171 0 1
## 172 0 2
## 173 0 1
## 176 0 1
## 182 0 1
## 184 0 1
## 192 0 2
## 201 0 1
## 203 0 1
## 204 0 1
## 207 0 2
## 208 0 2
## 210 0 1
## 213 0 2
## 214 0 3
## 219 0 3
## 220 0 1
## 224 0 1
## 226 0 1
## 230 0 1
## 233 0 1
## 234 0 1
## 238 0 1
## 239 0 1
## 241 0 1
## 242 0 1
## 246 0 1
## 250 0 1
## 251 0 1
## 252 0 1
## 253 0 1
## 255 0 1
## 261 0 1
## 263 0 1
## 264 0 1
## 269 0 1
## 270 0 1
## 273 0 1
## 280 0 1
## 288 0 1
## 294 0 1
## 295 0 1
## 297 0 1
## 303 0 1
## 307 0 1
## 308 0 1
## 309 0 1
## 323 0 1
## 341 0 1
## 352 0 1
## 360 0 2
## 380 0 1
## 410 0 1
## 423 0 1
## 424 0 2
## 425 0 1
## 447 0 1Nhìn vào các bảng ta thấy trong cơ sở dữ liệu này.
Nếu
al(Albumin) trong máu lớn hơn 0 thì đều sẽ bị bệnh.Những ai có nồng độ Creatinine từ ngường
scbằng 1.4 đến 76 đều bị bệnh, chỉ có 10 ngưỡng đầu từ 0.4 đến 1.3 có sự phân hóa giữa bị bênh và không bệnh.Trong 5 mốc khảo sát
sg(tỷ trọng nước tiểu) thì 3 mốc đầu đều ứng với người bị bệnh.
Đây là những nguyên nhân gây ra sự khớp hoàn hảo. Nếu muốn dự đoán thật chính xác ta sẽ dùng mô hình chứa 3 biến này luôn.
Nếu muốn giảm còn ít biến nhất và tránh việc R không phân tích được
mô hình, ta có thể xây dựng mô hình với chỉ 2 biến hemo và
bgr.
# 1. Xây dựng mô hình Logistic với 2 biến hemo và bgr
model_logic <- glm(target ~ hemo + bgr , data = train_data, family = "binomial")
# xem summary
summary(model_logic)##
## Call:
## glm(formula = target ~ hemo + bgr, family = "binomial", data = train_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 24.455508 3.634715 6.728 1.72e-11 ***
## hemo -2.138509 0.296430 -7.214 5.42e-13 ***
## bgr 0.036730 0.009285 3.956 7.62e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 422.37 on 319 degrees of freedom
## Residual deviance: 108.14 on 317 degrees of freedom
## AIC: 114.14
##
## Number of Fisher Scoring iterations: 8Ta thấy với 2 biến này mô hình được thiết lập thoát khỏi sự khớp hoàn hảo và R ra được 2 biến đều có ý nghĩa cao (\(p<0.001\)).
3.3.3 Phân phối của dữ liệu
Dựa trên heatmap vừa vẽ ta quan tâm sâu hơn về phân phối trong cơ sở dữ liệu.
library(ggplot2)
library(patchwork) # Để ghép các biểu đồ lại với nhau
# 1. Biểu đồ cột cho biến mục tiêu (Target)
p1 <- ggplot(df_clean, aes(x = factor(target), fill = factor(target))) +
geom_bar() +
scale_x_discrete(labels = c("Khỏe mạnh", "Bệnh thận")) +
labs(title = "Tỷ lệ mắc bệnh trong mẫu", x = "Trạng thái", y = "Số lượng") +
theme_minimal() + theme(legend.position = "none")
# 2. Biểu đồ phân phối của eGFR (Biến phụ thuộc hồi quy tuyến tính)
p2 <- ggplot(df_clean, aes(x = egfr)) +
geom_histogram(bins = 30, fill = "skyblue", color = "white") +
geom_density(aes(y = ..count.. * 5), color = "red") + # Thêm đường cong mật độ
labs(title = "Phân bố mức lọc cầu thận (eGFR)", x = "eGFR (ml/min/1.73m2)", y = "Tần suất") +
theme_minimal()
# 3. Boxplot so sánh Hemoglobin, bgr (Biến tương quan mạnh nhất)
p3 <- ggplot(df_clean, aes(x = factor(target), y = hemo, fill = factor(target))) +
geom_boxplot() +
scale_x_discrete(labels = c("Khỏe mạnh", "Bệnh thận")) +
labs(title = "Hemoglobin theo trạng thái bệnh", x = "", y = "Hemo (gms)") +
theme_minimal() + theme(legend.position = "none")
p4 <- ggplot(df_clean, aes(x = factor(target), y = bgr, fill = factor(target))) +
geom_boxplot() +
scale_x_discrete(labels = c("Khỏe mạnh", "Bệnh thận")) +
labs(title = "Đường huyết theo trạng thái bệnh", x = "", y = "BGR (mg/dl)") +
theme_minimal() + theme(legend.position = "none")
# 4. So sánh sg, hemo, bu, sc, htn, age (Biến quan trọng tính eGFR)
p5 <- ggplot(train_data, aes(x = sg, y = egfr)) +
geom_point(alpha = 0.5, color = "steelblue") +
theme_minimal()
p6 <- ggplot(train_data, aes(x = hemo, y = egfr)) +
geom_point(alpha = 0.5, color = "steelblue") +
theme_minimal()
p7 <- ggplot(train_data, aes(x = bu, y = egfr)) +
geom_point(alpha = 0.5, color = "steelblue") +
theme_minimal()
p8 <- ggplot(train_data, aes(x = sc, y = egfr)) +
geom_point(alpha = 0.5, color = "steelblue") +
theme_minimal()
p9 <- ggplot(train_data, aes(x = htn, y = egfr)) +
geom_point(alpha = 0.5, color = "steelblue") +
theme_minimal()
p10 <- ggplot(train_data, aes(x = age, y = egfr)) +
geom_point(alpha = 0.5, color = "steelblue") +
theme_minimal()
# Ghép 4 biểu đồ lại thành một khung hình duy nhất
(p1 | p2) / (p3 | p4) / (p5/p6)/(p7/p8)/(p9/p10)
- Quan sát các biểu đồ ta thấy có một số vấn đề về giá trị ngoại lệ, từ đây đặt ra việc cần làm là nên làm sạch thêm dữ liệu.
3.3.3.1 Làm sạch y khoa
Ta thấy rằng có 2 tương quan giữa sc và
egfr có sự bất thường. Về y học, nồng độ Creatinine trên 20
mg/dl thể hiện chức năng thận hầu như không còn (mức lọc cầu thận bằng
0) nên ta sẽ loại chúng ra. Tương tự nồng độ Ure trong máu
(bu) hơn 200 mg/dl là biểu hiện của nhiễm độc nghiêm trọng
nên ta không cần tính toán thêm mà phải có phương án xử lý trực
tiếp.
Vậy với vấn đề này ta sẽ cần điều chỉnh lại tập
train_data và test_data của chúng ta. Về
nguyên tắc thống kê ta sẽ làm sạch trên tập data_clean rồi
chia thành 2 tập lại từ đâu.
library(tidyverse)
# 1. Lọc bỏ ngoại lệ cho bu và sc trên toàn bộ dữ liệu sạch
df_final <- df_clean %>%
filter(sc <= 15, bu <= 200) # Ngưỡng đề xuất dựa trên biểu đồ bạn đã xem
# 2. Thực hiện chia dữ liệu (80/20)
set.seed(1000) # Giữ seed để kết quả ổn định
train_indices <- sample(1:nrow(df_final), 0.8 * nrow(df_final))
train_data <- df_final[train_indices, ]
test_data <- df_final[-train_indices, ]
# 3. Kiểm tra lại kích thước của 2 tập dữ liệu
cat("Số lượng quan sát tập Huấn luyện (Train):", nrow(train_data), "\n")## Số lượng quan sát tập Huấn luyện (Train): 306## Số lượng quan sát tập Kiểm tra (Test): 773.3.3.2 Đánh giá biến
Kiểm tra tính chuẩn : Qua biểu đồ Histogram, ta thấy chỉ số eGFR có phân phối lệch phải. Điều này cho thấy sự hiện diện của một số giá trị ngoại lai (outliers) ở nhóm bệnh nhân có chức năng thận cao. Sự lệch chuẩn này có thể ảnh hưởng đến giả định về phân phối của phần dư trong mô hình hồi quy tuyến tính, do đó cần được lưu ý khi diễn giải các khoảng tin cậy và giá trị \(p\)-value
Xác nhận giả thuyết:
Nhìn vào biểu đồ Boxplot của Hemoglobin theo trạng thái bệnh, ta thấy nhóm “Bệnh thận” có dải màu thấp hơn hẳn nhóm “Khỏe mạnh”. Điều này chứng minh cho con số \(-0.73\) trong Heatmap là có thật và rất rõ ràng rằng nếu bị bênh thì Hemoglobin sẽ thấp. Tương tự, nhóm bệnh thận có mức đường huyết cao hơn và biến động lớn hơn. Do đó, việc chọn 2 biến này để làm mô dình hồi quy logistic là phù hợp.
Nhìn vào các biểu đồ điểm ta thấy khó có thể xây dựng một mô hình hồi quy tuyến tính đơn biến nào với biến đầu ra là
egfrvà biến đầu vào là 1 trong các biếnsg,hemo,bu,sc,htn,age. Tuy nhiên, ta thấy giữaegfrvàsccó vẻ có mối liên hệ đặc biệt và có thể xây dựng một mô hình hồi quy đa thức giữa chúng.
3.4 Thực hiện xây dựng mô hình hồi quy
Như đã phân tích, ta sẽ có 3 mô hình
Mô hình hồi quy tuyến tính đa thức giữa
egfrvàsc.Mô hình hồi quy tuyến tính đa biến giữa
egfrvàsg,hemo,bu,sc,htn,ageMô hình hồi quy Logistic giữa
targetvàhemo,bgr
3.4.1 Mô hình hồi quy tuyến tính đa thức
Ta thự hiện xây dựng mô hình hồi quy đa thức với biến đầu vào là
sc (Chỉ số Creatinine) và biến đầu ra là egpr
(Chỉ số mức độ lọc cầu thận). Khác với các mô hình khác, mô hình đa thức
cần thử nghiệm để chọn ra bậc phù hợp, do đó, ta cần chia tiếp tập
train_data thành 2 phần.
# 1. Thiết lập seed để kết quả chia ngẫu nhiên có thể lặp lại được
set.seed(1000)
# 2. Xác định số lượng dòng cho tập huấn luyện (80% dữ liệu)
n <- nrow(train_data)
sample_polysize <- floor(0.8 * n)
# 3. Lấy ngẫu nhiên chỉ số các dòng
train_polyindices <- sample(seq_len(n), size = sample_polysize)
# 4. Chia dữ liệu thành 2 phần
train_polydata <- train_data[train_polyindices, ]
test_polydata <- train_data[-train_polyindices, ]
# 5. Kiểm tra lại kích thước của 2 tập dữ liệu
cat("Số lượng quan sát tập Huấn luyện đa thức (Train):", nrow(train_polydata), "\n")## Số lượng quan sát tập Huấn luyện đa thức (Train): 244## Số lượng quan sát tập Kiểm tra đa thức (Test): 62Sau khi đã có tập huấn luyện và tập kiểm tra bậc tối ưu cho mô hình hồi quy đa thức, ta cây dựng mô hình.
library(tidyverse)
# 1. Khởi tạo dataframe để lưu kết quả các bậc
poly_comparison <- data.frame(
Bac = 1:5,
R2_Hieu_Chinh = NA,
RMSE_Tap_Kiem_Tra = NA
)
# 2. Vòng lặp chạy thử từ bậc 1 đến bậc 5
for (i in 1:5) {
# Xây dựng mô hình đa thức bậc i
# raw = TRUE để dễ dàng đọc hệ số beta tương ứng với x^i
model_poly <- lm(egfr ~ poly(sc, i, raw = TRUE), data = train_polydata)
# Lấy Adjusted R-squared
poly_comparison$R2_Hieu_Chinh[i] <- summary(model_poly)$adj.r.squared
# Dự báo trên tập test để tính RMSE (Root Mean Squared Error)
preds <- predict(model_poly, newdata = test_polydata)
poly_comparison$RMSE_Tap_Kiem_Tra[i] <- sqrt(mean((test_polydata$egfr - preds)^2))
}
# 3. Hiển thị bảng kết quả so sánh
print("BẢNG SO SÁNH ĐỘ CHÍNH XÁC QUA CÁC BẬC:")## [1] "BẢNG SO SÁNH ĐỘ CHÍNH XÁC QUA CÁC BẬC:"## Bac R2_Hieu_Chinh RMSE_Tap_Kiem_Tra
## 1 1 0.3866725 53.62300
## 2 2 0.6175932 44.34420
## 3 3 0.7898801 37.58419
## 4 4 0.8721495 35.63644
## 5 5 0.9225536 37.20720# 4. Tìm bậc có RMSE thấp nhất (đây thường là bậc tối ưu)
best_degree <- poly_comparison$Bac[which.min(poly_comparison$RMSE_Tap_Kiem_Tra)]
cat(" best degree:", best_degree, "\n")## best degree: 4Vậy bậc tối ưu cho mô hình đa thức là bậc 4. Ta xây dựng được mô hình hồi quy tuyến tính đa thức như sau:
# 1. Xây dựng mô hình đa thức bậc 4
# Lưu ý: Sử dụng tập train_data đã được lọc nhiễu (sc <= 15/20)
model_poly4 <- lm(egfr ~ poly(sc, 4, raw = TRUE), data = train_data)
# 2. Xem tóm tắt chi tiết kết quả
summary(model_poly4)##
## Call:
## lm(formula = egfr ~ poly(sc, 4, raw = TRUE), data = train_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -50.914 -16.714 -6.228 14.733 105.679
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.329e+02 4.992e+00 46.66 <2e-16 ***
## poly(sc, 4, raw = TRUE)1 -1.644e+02 6.146e+00 -26.75 <2e-16 ***
## poly(sc, 4, raw = TRUE)2 3.866e+01 2.003e+00 19.30 <2e-16 ***
## poly(sc, 4, raw = TRUE)3 -3.531e+00 2.328e-01 -15.17 <2e-16 ***
## poly(sc, 4, raw = TRUE)4 1.088e-01 8.638e-03 12.60 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.76 on 301 degrees of freedom
## Multiple R-squared: 0.8459, Adjusted R-squared: 0.8438
## F-statistic: 413.1 on 4 and 301 DF, p-value: < 2.2e-16Đánh giá
Residuals (Phần dư): Phần dư là sự chênh lệch giữa giá trị thực tế của bệnh nhân và giá trị mà mô hình dự báo được.
- Min (-50.91) & Max (105.68)
- Ý nghĩa: Dữ liệu vẫn còn khoảng lệnh tương đối nhiều. Trung vị (Median) là -6.228, cho thấy đa số các điểm dự báo chỉ lệch rất ít so với thực tế. Tuy nhiên, vẫn có những trường hợp cá biệt có eGFR thực tế cao hơn dự báo khoảng 105 đơn vị.
Coefficients (Hệ số hồi quy): Đây là các tham số xác định hình dạng đường cong của mô hình.
Estimate (Hệ số ước lượng):
- Intercept (232.9): Đây là giá trị eGFR dự báo khi Creatinine (sc) bằng 0. Trong y tế, điều này giúp thiết lập điểm bắt đầu của đường cong ở mức lọc cầu thận lý tưởng cao nhất.
- poly(sc, 4, raw = TRUE)1 (-164.4): Hệ số của bậc 1 (tuyến tính), cho
thấy xu hướng giảm mạnh của eGFR khi
scbắt đầu tăng. - Các bậc cao (bậc 2, 3, 4): Các hệ số và tạo nên độ cong uốn lượn, giúp mô hình khớp với dữ liệu sinh học thực tế thay vì chỉ là một đường thẳng cứng nhắc.
Std. Error (Sai số chuẩn): Các chỉ số này đều rất nhỏ so với Estimate (ví dụ: bậc 4 chỉ có 0.0086). Điều này chứng tỏ các hệ số hồi quy được ước lượng với độ tin cậy cao.
t value: Tất cả các giá trị tuyệt đối của t đều rất lớn (từ 12.60 đến 46.66), vượt xa ngưỡng chuẩn > 2, khẳng định sức mạnh của các biến dự báo.
Pr(>|t|) (p-value): Tất cả đều nhỏ hơn 2e-16 (gần như bằng 0). Kết luận: Mọi bậc của đa thức (
sc) đều có đóng góp cực kỳ quan trọng và có ý nghĩa thống kê vào việc giải thích eGFR.
Các chỉ số đánh giá sự hợp lý
Residual standard error (21.76): Trung bình dự báo của mô hình chỉ còn lệch khoảng 21.76 đơn vị eGFR.
Multiple R-squared (0.8459): Mô hình giải thích được 84.59% sự biến thiên của dữ liệu.
Adjusted R-squared (0.8438): Kết quả 84.38% là rất tốt đối với dữ liệu y sinh.
F-statistic (Kiểm định F): F-statistic (413.1) & p-value (< 2.2e-16). Chỉ số F rất cao và p-value siêu nhỏ xác nhận toàn bộ mô hình đa thức này có ý nghĩa toán học rất cao.
Phương trình dự báo:
\[ egfr = 232.9 - 164.4(sc) + 38.66(sc^2) - 3.531(sc^3) + 0.1088(sc^4) \]
Nhận xét mô hình: Mô hình hồi quy đa thức bậc 4 đã phản ánh 84,38% mối quan hệ phi tuyến giữa Creatinine và eGFR. Việc loại bỏ các ngoại lệ và áp dụng bậc 4 đã tạo ra một mô hình dự báo có độ tin cậy cao, đủ tiêu chuẩn để làm cơ sở cho các phân tích y khoa chuyên sâu hơn.
3.4.2 Mô hình hồi quy tuyến tính đa biến
Ta thực hiện xây dựng mô hình hồi quy tuyến tính đa biến.
# 1. Xây dựng mô hình hồi quy đa biến
# Biến đầu ra: egfr
# Biến đầu vào: sg, hemo, bu, sc, htn, age
model2 <- lm(egfr ~ sg + hemo + bu +htn + age, data = train_data)
# 2. Xuất kết quả chi tiết
summary_model2 <- summary(model2)
print(summary_model2)##
## Call:
## lm(formula = egfr ~ sg + hemo + bu + htn + age, data = train_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -120.858 -22.546 -6.146 12.461 135.629
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.336e+03 4.586e+02 -2.914 0.00384 **
## sg 1.416e+03 4.556e+02 3.109 0.00206 **
## hemo 3.639e+00 1.186e+00 3.067 0.00236 **
## bu -4.458e-01 7.266e-02 -6.136 2.68e-09 ***
## htn -1.499e+01 6.112e+00 -2.452 0.01476 *
## age -9.858e-01 1.359e-01 -7.251 3.54e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37.85 on 300 degrees of freedom
## Multiple R-squared: 0.5354, Adjusted R-squared: 0.5277
## F-statistic: 69.15 on 5 and 300 DF, p-value: < 2.2e-16## 2.5 % 97.5 %
## (Intercept) -2238.9174911 -433.9282496
## sg 519.8577544 2312.8460770
## hemo 1.3039371 5.9731978
## bu -0.5888316 -0.3028422
## htn -27.0189238 -2.9617553
## age -1.2532997 -0.7182317Đánh giá
Residuals (Phần dư): Phần dư cho biết mô hình đoán sai bao nhiêu.
Min (-120.86) và Max (135.63): Dải sai số ở mô hình đa biến vẫn còn lại đáng kể do các biến khác còn xuất hiện nhiều giastrij ngoại lệ.
Ý nghĩa: Việc thêm nhiều biến vào đã giúp mô hình kiểm soát được các giá trị ngoại lệ tốt hơn. Tuy nhiên, dải sai số vẫn lệch (đuôi dương 146 dài hơn đuôi âm -107), xác nhận dữ liệu \(eGFR\) vẫn còn tính chất lệch phải.
Coefficients (Hệ số hồi quy): Đây là phần giúp giải thích cơ chế tác động của các yếu tố lâm sàng đến chức năng lọc của thận (\(eGFR\)).
Intercept (-1336): Đây là hệ số chặn mang giá trị âm rất lớn. Về mặt toán học, đây là điểm xuất phát của mô hình, tuy nhiên về y học nó không có ý nghĩa thực tế vì các chỉ số như sg hay hemo không bao giờ bằng 0.
sg(1416): Đây là biến có tác động mạnh nhất trong mô hình với trị số \(p = 0.002\) (\(p < 0.01\)). Hệ số dương cực lớn (1416) cho thấy tỷ trọng nước tiểu là một chỉ báo quan trọng: khi tỷ trọng nước tiểu tăng (thể hiện khả năng cô đặc của thận tốt), chỉ số \(eGFR\) cũng tăng mạnh tương ứng.hemo(3.64): Có ý nghĩa thống kê mạnh (\(p = 0.002 < 0.01\)). Với mỗi 1 g/dL Hemoglobin tăng thêm, \(eGFR\) dự báo tăng thêm 3.64 đơn vị. Điều này khẳng định mối liên hệ chặt chẽ giữa tình trạng thiếu máu và sự suy giảm chức năng thận.bu(-0.45): Biến này có ý nghĩa thống kê cực kỳ mạnh (\(p < 0.001\)). Hệ số âm cho thấy cứ mỗi đơn vị Ure máu tăng lên, thận lọc kém đi 0.45 đơn vị. Ure là chất độc mà thận cần đào thải, nên nồng độ này trong máu tăng tỷ lệ nghịch với sức khỏe của thận.htn(-14.99): Biến định tính về tăng huyết áp có ý nghĩa thống kê (\(p = 0.014 < 0.05\)). Một người bị tăng huyết áp sẽ có chỉ số \(eGFR\) thấp hơn khoảng 15 đơn vị so với người không bị, khi các yếu tố khác được giữ nguyên. Đây là tác động tiêu cực rất đáng kể của huyết áp lên màng lọc cầu thận.age(-0.99): Có ý nghĩa thống kê rất mạnh (\(p < 0.001\)). Cứ mỗi năm già đi, mức lọc cầu thận giảm trung bình khoảng 0.99 đơn vị (gần 1 đơn vị mỗi năm). Điều này phản ánh sự lão hóa tự nhiên của chức năng thận theo thời gian.
Đánh giá mức độ hợp lý:
Residual standard error (37.85): Sai số trung bình của mô hình là 37.85 đơn vị. Con số này cho thấy độ lệch trung bình giữa giá trị dự báo và giá trị thực tế đã được cải thiện (giảm xuống) so với các phiên bản trước đó.
Adjusted R-squared (0.5277): Mô hình đa biến này giải thích được 52.77% sự biến thiên của \(eGFR\). Đây là một mức độ giải thích khá ổn định cho thấy sự kết hợp của các biến lâm sàng có khả năng bao quát hơn một nửa các yếu tố ảnh hưởng đến chức năng thận.
Kiểm định F (F-statistic): F-statistic (69.15) với p-value < 2.2e-16:
- Ý nghĩa: Chỉ số F-statistic tăng lên 69.15 và p-value vẫn duy trì ở mức cực thấp. Điều này khẳng định mạnh mẽ rằng mô hình đa biến này có giá trị thực sự về mặt thống kê.
Khoảng tin cậy (Confint - 95%):
bu(-0.588 đến -0.302) vàage(-1.253 đến -0.718): Cả hai biến này có khoảng tin cậy rất hẹp và đều nằm hoàn toàn về phía âm. Điều này cho thấy tác động tiêu cực của nồng độ Ure máu và tuổi tác lên \(eGFR\) là cực kỳ ổn định và đáng tin cậy trong mẫu nghiên cứu.htn(-27.018 đến -2.961): Toàn bộ khoảng tin cậy đều mang giá trị âm. Ta có thể chắc chắn 95% rằng sự hiện diện của tình trạng cao huyết áp sẽ kéo chỉ số \(eGFR\) xuống từ 2.96 đến tận 27.02 đơn vị.hemo(1.303 đến 5.973): Toàn bộ khoảng tin cậy đều nằm ở miền dương. Điều này xác nhận rằng với mức tin cậy 95%, việc tăng nồng độ Hemoglobin trong máu chắc chắn sẽ cải thiện chỉ số chức năng thận \(eGFR\).sg(519.86 đến 2312.85): Đây là biến có khoảng tin cậy rộng nhất nhưng hoàn toàn dương, cho thấy tỷ trọng nước tiểu là một nhân tố dự báo có sức tác động rất lớn đến \(eGFR\)..
Nhận xét mô hình: Mô hình đa biến mới với các chỉ số được cập nhật cho thấy tính logic cao hơn. Việc loại bỏ các biến không cần thiết đã giúp tăng chỉ số F-statistic, đồng thời các khoảng tin cậy đều không chứa giá trị 0, khẳng định vai trò quan trọng của từng biến còn lại trong việc dự báo chức năng thận.
3.4.3 Mô hình hồi quy Logistic
Khác với hồi quy tuyến tính (dự báo chỉ số), hồi quy Logistic giúp ta dự báo xác suất (0% đến 100%) một bệnh nhân có bị bệnh thận mạn (CKD) hay không.
# 1. Xây dựng mô hình Logistic
# Biến phụ thuộc: target (1: ckd, 0: notckd)
# Biến độc lập: hemo, bgr
model3 <- glm(target ~ hemo + bgr, data = train_data, family = binomial)
# 2. Xem kết quả chi tiết
summary(model3)##
## Call:
## glm(formula = target ~ hemo + bgr, family = binomial, data = train_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 21.444360 3.238971 6.621 3.57e-11 ***
## hemo -1.872531 0.259124 -7.226 4.96e-13 ***
## bgr 0.031669 0.008539 3.709 0.000208 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 409.86 on 305 degrees of freedom
## Residual deviance: 125.94 on 303 degrees of freedom
## AIC: 131.94
##
## Number of Fisher Scoring iterations: 8# 3. Tính Odds Ratio (Tỷ số số chênh) - Chỉ số "vàng" trong y tế
odds_ratios <- exp(coef(model3))
print("Odds Ratios:")## [1] "Odds Ratios:"## (Intercept) hemo bgr
## 2.056682e+09 1.537341e-01 1.032176e+00Đánh giá
Hệ số hồi quy (Coefficients): Giải thích cơ chế tác động đến xác suất mắc bệnh .
- Estimate (Hệ số ước lượng):
Intercept (21.444): Đây là giá trị hằng số trong phương trình logit.
hemo(-1.8725): Hệ số âm cho thấy mối tương quan nghịch. Khi chỉ số Hemoglobin tăng 1 đơn vị, xác suất mắc bệnh thận giảm xuống đáng kể. Điều này phù hợp với y khoa vì thiếu máu (hemo thấp) là biểu hiện đặc trưng của suy thận.bgr(0.0317): Hệ số dương cho thấy mối tương quan thuận. Khi đường huyết (bgr) tăng, xác suất mắc bệnh thận cũng tăng theo.
- Estimate (Hệ số ước lượng):
Tỷ số số chênh (Odds Ratio - ): :
Đối với
hemo: (0.15) Nghĩa là khi tăng 1 đơn vị Hemoglobin, cơ hội (odds) mắc bệnh chỉ còn bằng 0.15 lần so với trước (giảm 85%).Đối với
bgr: (1.032) Nghĩa là khi tăng 1 đơn vị đường huyết, cơ hội mắc bệnh tăng thêm khoảng 3.2%.
Pr(>|z|) (p-value): Cả hai biến đều có p-value cực nhỏ (
hemo: vàbgr: ), nhỏ hơn rất nhiều so với ngưỡng . Kết luận: Cả Hemoglobin và Đường huyết đều là những nhân tố dự báo quan trọng cho bệnh thận.Đánh giá mức độ phù hợp của mô hình:
Null deviance (409.86) & Residual deviance (125.94): Giá trị Deviance đã giảm mạnh từ 409.86 xuống còn 125.94 (giảm khoảng 284 đơn vị) sau khi thêm các biến dự báo. Điều này chứng tỏ mô hình có khả năng giải thích dữ liệu rất tốt.
AIC (131.94): Chỉ số cho thấy mô hình vừa đủ gọn nhẹ (ít biến) nhưng vẫn giữ được độ chính xác cao.
Fisher Scoring iterations (8): Thuật toán hội tụ chỉ sau 8 lần lặp. Đây là minh chứng cho thấy hiện tượng khớp hoàn toàn đã biến mất, mô hình hiện tại ổn định về mặt toán học.
Nhận xét mô hình: Mô hình Logistic rút gọn với hai
biến hemo và bgr là một lựa chọn tốt. Nó đảm
bảo được tính đơn giản nhưng lại mang lại hiệu quả dự báo mạnh và các
biến đều có ý nghĩa lâm sàng rõ rệt. Kết quả này cho thấy chức năng thận
chịu ảnh hưởng lớn từ tình trạng thiếu máu và kiểm soát đường huyết của
bệnh nhân.
4 KẾT QUẢ
Trong phần này, chúng ta đưa ra kết quả từ ba mô hình hồi quy đã xây dựng trên tập huấn luyện (80% dữ liệu) và đánh giá khả năng dự báo của chúng dựa trên tập kiểm tra (20% dữ liệu).
4.1 Mô hình 1: Hồi quy tuyến tính đa thức
Mô hình này xem xét tác động trực tiếp của nồng độ Creatinine
(sc) lên mức lọc cầu thận (egfr).
Biểu đồ
library(ggplot2)
# Xây dựng mô hình với bậc tối ưu
final_model <- lm(egfr ~ poly(sc, best_degree, raw = TRUE), data = train_polydata)
# Vẽ biểu đồ đường cong hồi quy
ggplot(train_data, aes(x = sc, y = egfr)) +
geom_point(alpha = 0.4, color = "steelblue") +
stat_smooth(method = "lm", formula = y ~ poly(x, best_degree),
color = "darkred", size = 1.2) +
labs(title = paste("Mô hình hồi quy đa thức bậc", best_degree, "giữa sc và egfr"),
subtitle = paste("R2 Hiệu chỉnh:", round(summary(final_model)$adj.r.squared, 4)),
x = "Serum Creatinine (sc)", y = "eGFR") +
theme_minimal() 
Về biểu đồ ta dễ thấy phần lớn dữ phiệu phân tán nằm xung quanh đường hồi quy.
Kiểm định toán học
library(ggplot2)
# Thiết lập hiển thị 4 biểu đồ trên cùng 1 khung hình
par(mfrow = c(2, 2))
# Vẽ biểu đồ chẩn đoán cho mô hình đa thức bậc 4
plot(model_poly4, main = "Kiểm định giả định mô hình bậc 4") #lệnh này tự động tạo ra 4 biểu đồ quan trọng nhất
1. Biểu đồ Residuals vs Fitted (Kiểm tra tính tuyến tính)
Quan sát: Đường màu đỏ uốn lượn rất nhẹ và gần như bám sát đường đứt đoạn \(0\).
Nhận xét: Điều này chứng tỏ mô hình đa thức bậc 4 đã nắm bắt được hầu hết cấu trúc của dữ liệu. Đường màu đỏ có dạng chữ U thể hiện sự phi tuyến giữa biến đầu vào và biến đầu ra. Các điểm số 173, 258, 338 được hệ thống đánh dấu là những điểm có phần dư lớn cần chú ý.
2. Biểu đồ Normal Q-Q (Kiểm tra phân phối chuẩn)
Quan sát: Phần lớn các điểm dữ liệu nằm sát đường chéo đứt đoạn. Tuy nhiên, ở phía đầu bên phải (giá trị cao), các điểm dữ liệu có xu hướng lệch lên trên khỏi đường thẳng.
Nhận xét: Phần dư của mô hình tuân theo phân phối chuẩn ở dải giá trị trung tâm nhưng có hiện tượng lệch ở các giá trị \(eGFR\) cao. Điều này phản ánh thực tế dữ liệu y sinh thường có các ca bệnh đặc biệt mà mô hình toán học chưa thể khớp hoàn toàn.
3. Biểu đồ Scale-Location (Kiểm tra tính đồng nhất của phương sai)
Quan sát: Đường màu đỏ có xu hướng đi lên nhẹ khi giá trị dự báo (\(Fitted \ values\)) tăng lên.
Nhận xét: Có dấu hiệu nhẹ của hiện tượng phương sai thay đổi (heteroscedasticity). Cụ thể, khi mức lọc cầu thận dự báo càng cao, độ biến động của sai số cũng tăng lên một chút. Tuy nhiên, mức độ này không quá nghiêm trọng để làm hỏng mô hình.
4. Biểu đồ Residuals vs Leverage (Kiểm tra điểm ảnh hưởng)
Quan sát: Không có điểm dữ liệu nào vượt qua đường giới hạn Cook’s distance (đường đỏ đứt đoạn ở ngưỡng 0.5 và 1.0). Điểm số 81 có đòn bẩy (\(Leverage\)) khá cao nhưng vẫn nằm trong vùng an toàn.
Nhận xét: Việc ta loại bỏ các giá trị \(sc > 20\) đã phát huy tác dụng. Mô hình không bị chi phối bởi các số liệu đơn lẻ, đảm bảo tính tổng quát hóa cho toàn bộ tập dữ liệu.
Kiểm tra với test_data
# 1. Dự báo trên tập Test bằng mô hình bậc 4
preds_poly4 <- predict(model_poly4, newdata = test_data)
# 2. Tính sai số RMSE (Root Mean Squared Error)
rmse_poly4 <- sqrt(mean((test_data$egfr - preds_poly4)^2))
cat("RMSE trên tập Test:", rmse_poly4, "\n")## RMSE trên tập Test: 17.56918Chỉ số RMSE (Root Mean Squared Error): 17.57Ý nghĩa: Chỉ số này đo lường độ lệch trung bình giữa giá trị \(eGFR\) dự báo và giá trị thực tế của bệnh nhân trong tập dữ liệu kiểm tra độc lập. Với kết quả 17.57, mô hình cho thấy khả năng dự báo mức lọc cầu thận khá chính xác trên những dữ liệu mà nó chưa từng thấy trước đó.
Một điểm rất tốt là chỉ số RMSE trên tập Test (17.57) thấp hơn so với Residual Standard Error trên tập huấn luyện (21.76).
4.2 Mô hình 2: Hồi quy tuyến tính đa biến
Kiểm định toán học
1.Biểu đồ Residuals vs Fitted (Tính tuyến tính)
Quan sát: Đường màu đỏ có xu hướng dốc xuống nhẹ và các điểm thặng dư phân tán rộng hơn so với mô hình đa thức (trục Tung từ -150 đến 100).
Nhận xét: Mô hình đa biến có sai số lớn hơn đáng kể so với mô hình đa thức. Các điểm 125, 338, 161 là những điểm có sai số dự báo lớn nhất. Điều này giải thích tại sao \(R^2\) hiệu chỉnh chỉ đạt 52.77%.
2.Biểu đồ Normal Q-Q (Phân phối chuẩn):
Quan sát: Các điểm thặng dư bám khá sát đường chéo ở đoạn giữa nhưng lệch rất mạnh ở hai đầu, đặc biệt là phía đuôi bên phải.
Nhận xét: Giả định về phân phối chuẩn của sai số bị vi phạm nhẹ do hiện tượng lệch phải. Điều này thường gặp trong dữ liệu y tế khi có những ca bệnh nặng với chỉ số lâm sàng bất thường vượt xa mức trung bình.
3.Biểu đồ Scale-Location (Tính đồng nhất phương sai):
Quan sát: Đường màu đỏ có xu hướng đi lên rõ rệt khi giá trị dự báo tăng.
Nhận xét: Có hiện tượng phương sai thay đổi (Heteroscedasticity). Khi mức lọc cầu thận dự báo càng cao, sai số của mô hình càng có xu hướng biến động mạnh hơn.
4.Biểu đồ Residuals vs Leverage (Điểm ảnh hưởng):
Quan sát: Các điểm thặng dư hầu hết nằm trong vùng an toàn. Không có điểm nào vượt qua đường giới hạn Cook’s distance (0.5 hoặc 1.0).
Nhận xét: Mặc dù xuất hiện một số điểm có đòn bẩy cao như 173, 125, 720, chúng không gây ra sự thay đổi quá lớn làm lệch lạc các hệ số hồi quy của toàn mô hình.
Kiểm tra với test_data
# Tính RMSE cho mô hình đa biến
preds_multi <- predict(model2, newdata = test_data)
rmse_multi <- sqrt(mean((test_data$egfr - preds_multi)^2))
cat("RMSE đa biến:", rmse_multi)## RMSE đa biến: 35.90962Chỉ số này đo lường độ lệch trung bình giữa giá trị \(eGFR\) dự báo và thực tế trên tập dữ liệu
kiểm tra độc lập cho mô hình kết hợp 5 biến lâm sàng (sg,
hemo, bu, htn, age).
Với kết quả 35.91, sai số dự báo của mô hình đa biến nằm ở mức trung
bình và phản ánh đúng sai số đã quan sát được trên tập huấn luyện (\(RSE = 37.85\)).
Sai số của mô hình đa biến cao gấp hơn 2 lần so với mô hình đa thức
bậc 4. Điều này khẳng định rằng mặc dù việc kết hợp nhiều biến lâm sàng
giúp mang lại cái nhìn tổng quát về bệnh lý, nhưng để dự báo chính xác
chỉ số \(eGFR\), mô hình hóa mối quan
hệ phi tuyến của riêng Creatinine (sc) vẫn đem lại hiệu quả
hơn hẳn.
4.3 Mô hình 3: Hồi quy Logistic
- Các biến dự báo: Qua quá trình lựa chọn biến, mô hình tối ưu bao gồm
hai biến: Hemoglobin (
hemo) và Đường huyết (bgr). - Phương trình Logit:
\[ \ln\left(\dfrac{P}{1-P}\right) = 21.444 - 1.872 \cdot hemo + 0.032 \cdot bgr \]
Kiểm định toán học và kiểm tra với
test_data
## Installing package into 'C:/Users/win/AppData/Local/R/win-library/4.4'
## (as 'lib' is unspecified)## package 'pROC' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\win\AppData\Local\Temp\RtmpCoaQ1i\downloaded_packages## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var# 1. Tính toán xác suất dự báo trên tập Test
prob_test <- predict(model3, newdata = test_data, type = "response")
# 2. Tạo đối tượng roc
roc_obj <- roc(test_data$target, prob_test)## Setting levels: control = 0, case = 1## Setting direction: controls < cases# 3. Vẽ biểu đồ ROC
plot(roc_obj, col = "blue", main = paste("Đường cong ROC (AUC =", round(auc(roc_obj), 4), ")"))
abline(a = 0, b = 1, lty = 2, col = "red") # Đường 50/50 ngẫu nhiên
## Chi so AUC cua mo hinh la: 0.9794326Đường cong ROC: Biểu đồ cho thấy đường cong bám sát góc trên bên trái, thể hiện khả năng phân loại rất tốt.
Chỉ số AUC (Area Under the Curve): 0.9794. Ý nghĩa: Với AUC xấp xỉ 0.98, mô hình có xác suất tới 98% phân loại đúng một bệnh nhân có bệnh và một người khỏe mạnh. Đây là một chỉ số phân loại thuộc nhóm xuất sắc trong dự báo y sinh.
Ma trận nhầm lẫn với test_data
library(caret)
# 1. Đảm bảo biến thực tế là factor với 2 levels rõ ràng
actual_factor <- factor(test_data$target, levels = c("0", "1"))
# 2. Dự báo xác suất
probs <- predict(model3, newdata = test_data, type = "response")
# 3. Chuyển đổi xác suất thành nhãn (0 hoặc 1)
preds_numeric <- ifelse(probs > 0.5, 1, 0)
# 4. Quan trọng: Ép nhãn dự báo về factor với ĐÚNG 2 levels như actual_factor
preds_factor <- factor(preds_numeric, levels = c("0", "1"))
# 5. Chạy lại ma trận nhầm lẫn
library(caret)
confusionMatrix(preds_factor, actual_factor)## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 29 2
## 1 1 45
##
## Accuracy : 0.961
## 95% CI : (0.8903, 0.9919)
## No Information Rate : 0.6104
## P-Value [Acc > NIR] : 6.288e-13
##
## Kappa : 0.9186
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9667
## Specificity : 0.9574
## Pos Pred Value : 0.9355
## Neg Pred Value : 0.9783
## Prevalence : 0.3896
## Detection Rate : 0.3766
## Detection Prevalence : 0.4026
## Balanced Accuracy : 0.9621
##
## 'Positive' Class : 0
## Ta có một số nhận xét từ mô hình như sau:
Dự báo đúng: 74 trường hợp (29 ca khỏe mạnh và 45 ca bệnh).
Dự báo sai: Chỉ có 3 trường hợp (1 ca bị chẩn đoán nhầm là có bệnh và 2 ca bị bỏ sót bệnh).
Accuracy (Độ chính xác tổng quát): 96.1%. Mô hình phân loại đúng tình trạng bệnh cho hơn 96% bệnh nhân trong tập kiểm tra. Đây là mức độ chính xác rất cao, khẳng định sự kết hợp giữa Hemoglobin và Đường huyết mang lại khả năng dự báo rất tốt.
Sensitivity (Độ nhạy): 96.67%. Khả năng phát hiện chính xác những người khỏe mạnh là 96.67%.
Specificity (Độ đặc hiệu): 95.74%. Khả năng xác định đúng những người thực sự mắc bệnh thận (CKD) là 95.74%. Trong y tế, chỉ số này cực kỳ quan trọng vì nó giúp đảm bảo rất ít bệnh nhân bị bỏ quên trong quá trình sàng lọc.
Kappa: 0.9186. Chỉ số Kappa > 0.9 cho thấy sự đồng thuận giữa dự báo của mô hình và thực tế lâm sàng đạt mức gần như hoàn hảo (Almost Perfect Agreement), loại trừ yếu tố ngẫu nhiên.
Nhận xét tổng quát: Mô hình hồi quy Logistic với
hai biến hemo và bgr có ý nghĩa thống kê mạnh
và mang lại giá trị thực tiễn cao trong sàng lọc bệnh. Kết quả cho thấy
tình trạng thiếu máu và kiểm soát đường huyết là hai chỉ báo quan trọng
nhất để phân loại bệnh nhân trong mẫu nghiên cứu này.
5 BÀN LUẬN
5.1 So sánh mô hình đa thức với đa biến
Mô hình hồi quy đa thức bậc 4 đã thể hiện sự vượt trội về khả năng dự
báo so với cả mô hình tuyến tính đơn biến ban đầu và mô hình tuyến tính
đa biến. Cụ thể, mô hình đa thức chỉ sử dụng duy nhất biến Creatinine
(sc) nhưng đạt chỉ số \(R^2\) hiệu chỉnh lên tới 84.38%, cao hơn
đáng kể so với mức 52.77% của mô hình đa biến kết hợp 5 yếu tố lâm sàng
khác.
Sự khác biệt này phản ánh hai vấn đề cốt lõi về mặt thống kê và y
khoa: Bản chất phi tuyến tính của dữ liệu: Mối quan hệ giữa nồng độ
Creatinine và mức lọc cầu thận (eGFR) vốn dĩ mang tính
nghịch biến phi tuyến tính mạnh. Khi sử dụng mô hình đa biến dưới dạng
đường thẳng (Linear), chúng ta đã vô tình bỏ qua độ cong của dữ liệu,
dẫn đến sai số RMSE trên tập kiểm tra vọt lên mức 35.91.
Ngược lại, việc nâng lên bậc 4 giúp mô hình uốn lượn linh hoạt để
khớp với quy luật chức năng thận, giảm sai số RMSE xuống chỉ còn 17.57.
Vai trò của việc xử lý ngoại lệ: Việc chủ động loại bỏ các điểm dữ liệu
nhiễu (sc > 20) đã giúp mô hình đa thức đạt được trạng
thái ổn định tối ưu.
5.2 Mối liên hệ đa yếu tố
Mặc dù mô hình đa thức bậc 4 chiếm ưu thế về khả năng dự báo con số \(eGFR\), nhưng mô hình hồi quy đa biến lại cung cấp một cái nhìn toàn diện về các yếu tố bệnh lý đang tác động lên bệnh nhân.
Kết quả cho thấy các chỉ số xét nghiệm máu và nước tiểu cơ bản có mối liên hệ mật thiết với chức năng thận:
Tầm quan trọng của Hemoglobin (
hemo): Cả mô hình đa biến và Logistic đều xác nhận Hemoglobin là một chỉ báo quan trọng với hệ số dương (3.64) và có ý nghĩa thống kê mạnh (\(p = 0.002\)). Về mặt lâm sàng, điều này phản ánh tình trạng thiếu máu do thận giảm sản xuất hormone Erythropoietin khi bị tổn thương. Kết quả này khẳng định rằng cải thiện tình trạng thiếu máu có liên quan trực tiếp đến việc nâng cao mức lọc cầu thận của bệnh nhân.Tác động tiêu cực của Ure máu (
bu) và tăng huyết áp (htn): Ure máu mang hệ số âm (-0.45) với mức ý nghĩa cực cao (\(p < 0.001\)), minh chứng cho việc tích tụ chất độc nitrogen trong máu khi thận mất khả năng đào thải. Đồng thời, tình trạng tăng huyết áp (htn) làm giảm chỉ số \(eGFR\) trung bình tới 14.99 đơn vị. Khoảng tin cậy của biếnhtn(-27.02 đến -2.96) hoàn toàn nằm ở miền âm, khẳng định huyết áp cao là tác nhân trực tiếp gây xơ hóa màng lọc cầu thận.Chỉ số tỷ trọng nước tiểu (
sg) và tuổi tác (age): Hệ số dương rất lớn củasg(1416) cho thấy khả năng cô đặc nước tiểu là dấu hiệu vàng cho một quả thận khỏe mạnh. Ngược lại, biếnagemang dấu âm (-0.99) phản ánh quy luật lão hóa tự nhiên, khi cứ mỗi năm già đi, chức năng thận sẽ suy giảm khoảng 1 đơn vị \(eGFR\).Việc loại bỏ biến Creatinine (sc) ra khỏi mô hình đa biến cuối cùng không có nghĩa là nó không quan trọng, mà thực tế là do các biến nhưbuvàhemođã phản ánh một phần thông tin tương đồng, và mô hình đa biến này được tối ưu hóa để làm nổi bật các yếu tố nguy cơ lâm sàng khác ngoài chỉ số Creatinine truyền thống.
5.3 Ý nghĩa của mô hình Logistic
Nếu mô hình đa thức và đa biến đóng vai trò quan trọng trong việc theo dõi tiến triển và giải thích bệnh lý, thì mô hình hồi quy Logistic lại thể hiện giá trị thực tiễn trong công tác sàng lọc cộng đồng.
Với chỉ số AUC đạt 0.9794 và độ chính xác thực tế trên tập kiểm tra lên tới 96.1%, đây là một công cụ phân loại có độ tin cậy rất cao. Sức mạnh của mô hình này nằm ở tính tối giản và hiệu quả:
Khả năng phân loại xuất sắc: Chỉ số AUC tiệm cận 1.0 cùng với độ nhạy (96.67%) và độ đặc hiệu (95.74%) rất cao cho thấy mô hình cực kỳ nhạy bén trong việc nhận diện bệnh nhân mắc bệnh thận mạn tính, đồng thời hạn chế tối đa tỷ lệ chẩn đoán nhầm cho người khỏe mạnh. Điều này đặc biệt có ý nghĩa trong y tế dự phòng, giúp phát hiện sớm các ca bệnh tiềm ẩn để can thiệp kịp thời.
Tính ứng dụng từ các chỉ số cơ bản: Mô hình Logistic chỉ sử dụng hai biến là Hemoglobin (
hemo) và Đường huyết (bgr). Đây là hai chỉ số xét nghiệm máu phổ thông, rẻ tiền và có kết quả nhanh. Việc xác định rằng cứ tăng 1 đơn vị Hemoglobin thì cơ hội mắc bệnh giảm tới 85% (\(OR \approx 0.15\)) cung cấp một quy tắc lâm sàng đơn giản nhưng mạnh mẽ cho các bác sĩ trong việc đánh giá nguy cơ bệnh nhân ngay từ tuyến cơ sở.Mối liên hệ giữa tiểu đường và bệnh thận: Hệ số dương của biến
bgr(\(OR \approx 1.03\)) một lần nữa khẳng định vai trò của đường huyết cao như một yếu tố nguy cơ trực tiếp. Kết quả này cho thấy mô hình không chỉ là một thuật toán dự báo mà còn phản ánh đúng thực trạng bệnh lý: kiểm soát đường huyết tốt là chìa khóa để ngăn ngừa sự hình thành của bệnh thận mạn tính.
Tóm lại, mô hình Logistic đóng vai trò là sàng lọc đầu tiên trong quy trình chẩn đoán. Sự kết hợp giữa khả năng phân loại chính xác của mô hình Logistic và khả năng dự báo chỉ số chi tiết của mô hình đa thức bậc 4 sẽ tạo nên một bộ công cụ hỗ trợ ra quyết định lâm sàng toàn diện.
6 TÀI LIỆU THAM KHẢO
Sheldon M. Ross (2014), Introduction to Probability and Statistics for Engineers and Scientists, 5th Edition, Elsevier Academic Press.
Nguyễn Văn Tuấn (2014), Phân tích dữ liệu với R, Nhà xuất bản Tổng hợp Thành phố Hồ Chí Minh.
Đào Huy Cường (2025), Bài giảng Xác suất thống kê nâng cao, Khoa Toán - Tin học, Trường Đại học Sư phạm TP.HCM.