Keywords: Mô hình hồi quy logistic, LASSO, Dữ liệu về bệnh nhân mắc bệnh tim mạch, thử nghiệm dự báo, đánh giá mô hình.

1    Introduction

Năm 2020, Phó giáo sư lâm sàng nội khoa giáo dục tin học Robert Hoyt, MD (Đại học California tại Irvine - UCI) đã công bố một bộ dữ liệu Dự đoán Bệnh tim (Heart Disease Prediction - HDP), bao gồm 270 nghiên cứu trường hợp của các cá nhân được phân loại là mắc hoặc không mắc bệnh tim dựa trên kết quả của các phương pháp chẩn đoán y học về tim mạch. Tập dữ liệu này cung cấp cho ta một cái nhìn quan trọng trong việc đánh giá mối quan hệ giữa các yếu tố nguy cơ và sức khỏe tim mạch.

Trong lĩnh vực Khoa học thống kê, có nhiều tập dữ liệu mà ta cần quan tâm đến mối quan hệ giữa các loại dữ liệu với nhau. Như vậy, người ta cũng cố gắng khảo sát xem liệu một bệnh nhân bị mắc bệnh tim có mối quan hệ gì tới một số yếu tố như tuổi tác, giới tính, các yếu tố lâm sàng, chẩn đoán tim mạch… hay không. Người ta sẽ thu thập, quan sát và phân tích tập dữ liệu về các yếu tố đó của những người bị mắc và không bị mắc bệnh. Một phương pháp thống kê được sử dụng phổ biến là phân tích hồi quy (regression analysis). Đây là một phương pháp phân tích thống kê dùng để xác định mối quan hệ giữa một biến phụ thuộc (biến kết quả) với một hoặc nhiều biến độc lập (biến dự đoán).

Từ đó, mục tiêu của ta là sử dụng các phương pháp phân tích hồi quy để khảo sát tập dữ liệu của Robert Hoyt để thu lại được các kết quả ý nghĩa về mặt thống kê cũng như về mặt thực tiễn. Trong bài tiểu luận này, tác giả sử dụng một phần của tập dữ liệu, trong đó chọn ra 100 dữ liệu ngẫu nhiên và mỗi bệnh nhân (mỗi dữ liệu) được ghi nhận số liệu về 8 yếu tố lâm sàng, cũng là 9 biến giải thích, đó là các biến về độ tuổi (Age), giới tính (Sex), chỉ số huyết áp (BP), mức cholesterol (C), kết quả điện tâm đồ (EKG), nhịp tim cực đại (HR), Đoạn ST chênh xuống (ST) và số lượng mạch máu được thấy trên hình ảnh huỳnh quang (number of vessels fluroscopy - NVF). Dưới đây là toàn bộ bộ dữ liệu ta sẽ xét trong bài tiểu luận này.

Tập dữ liệu về rủi ro mắc bệnh tim và một số chẩn đoán lâm sàng
No.  Age Sex BP Cholesterol EKG Max.HR ST.depression Num.of.ves.flu Heart.Disease
1 67 0 115 564 2 160 1.6 0 0
2 67 1 100 299 2 125 0.9 2 1
3 68 0 120 211 2 115 1.5 0 0
4 34 1 118 182 2 174 0.0 0 0
5 62 0 138 294 0 106 1.9 3 1
6 51 1 140 298 0 122 4.2 3 1
7 46 1 150 231 0 147 3.6 0 1
8 67 1 125 254 0 163 0.2 2 1
9 50 1 129 196 0 163 0.0 0 0
10 42 1 120 240 0 194 0.8 0 0
11 56 0 134 409 2 150 1.9 2 1
12 41 1 110 172 2 158 0.0 0 1
13 42 0 102 265 2 122 0.6 0 0
14 53 1 130 246 2 173 0.0 3 0
15 43 1 130 315 0 162 1.9 1 0
16 56 1 132 184 2 105 2.1 1 1
17 52 1 108 233 0 147 0.1 3 0
18 62 0 140 394 2 157 1.2 0 0
19 70 1 160 269 0 112 2.9 1 1
20 54 1 140 239 0 160 1.2 0 0
21 70 1 145 174 0 125 2.6 0 1
22 54 1 108 309 0 156 0.0 0 0
23 35 1 126 282 2 156 0.0 0 1
24 48 1 124 255 0 175 0.0 2 0
25 55 0 135 250 2 161 1.4 0 0
26 58 0 100 248 2 122 1.0 0 0
27 54 0 110 214 0 158 1.6 0 0
28 69 0 140 239 0 151 1.8 2 0
29 77 1 125 304 2 162 0.0 3 1
30 68 1 118 277 0 151 1.0 1 0
31 69 1 140 254 2 146 2.0 3 1
32 69 1 160 234 2 131 0.1 1 0
33 65 1 138 282 2 174 1.4 1 1
34 45 0 138 236 2 152 0.2 0 0
35 34 0 118 210 0 192 0.7 0 0
36 57 1 132 207 0 168 0.0 0 0
37 64 1 145 212 2 132 2.0 2 1
38 59 1 138 271 2 182 0.0 0 0
39 50 1 140 233 0 163 0.6 1 1
40 51 1 125 213 2 125 1.4 1 0
41 54 1 192 283 2 195 0.0 1 1
42 53 1 123 282 0 95 2.0 2 1
43 52 1 112 230 0 160 0.0 1 1
44 40 1 110 167 2 114 2.0 0 1
45 58 1 132 224 2 173 3.2 2 1
46 41 0 112 268 2 172 0.0 0 0
47 41 1 112 250 0 179 0.0 0 0
48 50 0 120 219 0 158 1.6 0 0
49 58 1 125 300 2 171 0.0 2 1
50 54 0 108 267 2 167 0.0 0 0
51 51 0 130 256 2 149 0.5 0 0
52 46 0 105 204 0 172 0.0 0 0
53 55 1 140 217 0 111 5.6 0 1
54 45 1 128 308 2 170 0.0 0 0
55 56 1 120 193 2 162 1.9 0 0
56 66 0 178 228 0 165 1.0 2 1
57 38 1 120 231 0 182 3.8 0 1
58 62 0 150 244 0 154 1.4 0 1
59 55 1 130 262 0 155 0.0 0 0
60 58 1 128 259 2 130 3.0 2 1
61 43 1 110 211 0 161 0.0 0 0
62 64 0 180 325 0 154 0.0 0 0
63 50 0 110 254 2 159 0.0 0 0
64 53 1 130 197 2 152 1.2 0 0
65 64 0 130 303 0 122 2.0 2 0
66 60 1 125 258 2 141 2.8 1 1
67 51 1 140 299 0 173 1.6 0 1
68 55 1 160 289 2 145 0.8 1 1
69 52 1 120 325 0 172 0.2 0 0
70 68 1 180 274 2 150 1.6 0 1
71 39 1 140 321 2 182 0.0 0 0
72 53 0 130 264 2 143 0.4 0 0
73 62 0 140 268 2 160 3.6 2 1
74 51 0 140 308 2 142 1.5 1 0
75 60 1 130 253 0 144 1.4 1 1
76 65 1 110 248 2 158 0.6 2 1
77 65 0 155 269 0 148 0.8 0 0
78 60 1 140 185 2 155 3.0 0 1
79 60 1 145 282 2 142 2.8 2 1
80 54 1 120 188 0 113 1.4 1 1
81 44 1 130 219 2 188 0.0 0 0
82 44 1 112 290 2 153 0.0 1 1
83 47 1 138 257 2 156 0.0 0 0
84 51 1 110 175 0 123 0.6 0 0
85 71 0 160 302 0 162 0.4 2 0
86 61 1 150 243 0 137 1.0 0 0
87 55 1 132 353 0 132 1.2 1 1
88 64 1 140 335 0 158 0.0 0 1
89 43 1 150 247 0 171 1.5 0 0
90 58 0 120 340 0 172 0.0 0 0
91 60 1 130 206 2 132 2.4 2 1
92 58 1 120 284 2 160 1.8 0 1
93 49 1 130 266 0 171 0.6 0 0
94 48 1 110 229 0 168 1.0 0 1
95 52 1 172 199 0 162 0.5 0 0
96 44 1 120 263 0 173 0.0 0 0
97 56 0 140 294 2 153 1.3 0 0
98 57 1 140 192 0 148 0.4 0 0
99 59 1 150 212 0 157 1.6 0 0
100 56 1 125 249 2 144 1.2 1 1

Với bộ dữ liệu này, chúng ta có cơ hội đánh giá cách các đặc điểm này tương tác với nhau để xác định mức độ nguy cơ của một cá nhân trong việc phát triển các vấn đề tim mạch dẫn đến suy tim hoặc đột quỵ. Với kiến thức này, chúng ta có thể tạo ra các chiến lược phòng ngừa vượt ra ngoài những gì điều trị y tế truyền thống có thể làm bằng cách xác định sớm những người có nguy cơ và hỗ trợ các chuyên gia y tế trong việc điều trị họ tốt hơn. Bằng cách phân tích kết hợp các biến lâm sàng được giải thích ở đây, chúng ta có một công cụ mạnh mẽ trong tay để cố gắng chống lại bệnh tim mạch trước khi nó có cơ hội hình thành.

Bài tiểu luận này sẽ bao gồm năm phần chính:

  • Phần đầu tiên giới thiệu bộ dữ liệu được sử dụng để phân tích hồi quy.
  • Phần thứ hai trình bày về các cơ sở lý thuyết về các mô hình hồi quy được sử dụng trong bài tiểu luận như hồi quy logistic và hồi quy LASSO.
  • Phần thứ ba tập trung mô tả tập dữ liệu và áp dụng các mô hình hồi quy để phân tích mối quan hệ của biến đầu ra với các biến giải thích. Có đoạn code R tương ứng.
  • Phần thứ tư tổng kết một số kết quả chính thu được qua việc vận dụng mô hình hồi quy, thực hiện đánh giá các mô hình hồi quy và đánh giá sai số của mỗi mô hình.
  • Phần cuối cùng bàn về một số nhận xét của riêng tác giả và tài liệu tham khảo.

2   Preliminaries

2.1   Multiple Linear Regression

Đầu tiên, ta phát biểu một số ký hiệu được sử dụng xuyên suốt trong bài tiểu luận này.

Trong một bộ dữ liệu, ta ký hiệu mỗi bộ biến đầu vào \(x = (x_1, x_2, \ldots, x_n)\) (trong đó \(n\) là số lượng biến đầu vào). Giả sử ứng với mỗi biến đầu vào \(x\), ta thu được một biến đầu ra \(Y\) (còn gọi là biến phụ thuộc). Mối quan hệ giữa \(Y\) và \(x_1, x_2, \ldots, x_p\) có thể được xấp xỉ bởi một mô hình hồi quy dạng \[Y = f(x_1, x_2, \ldots, x_n ) + \varepsilon,\] trong đó \(\varepsilon\) là biến ngẫu nhiên đại diện cho sai số của mô hình.

Khi hàm \(f\) biểu thị mối quan hệ tuyến tính giữa \(Y\) và \(x_1, x_2, \ldots, x_p\), ta có thể viết mô hình hồi quy tuyến tính dưới dạng \[Y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \cdots + \beta_nx_n + \varepsilon,\] trong đó \(\beta_0, \beta_1, \ldots, \beta_n\) được gọi là các hệ số hồi quy.

Giả sử ta quan sát được bảng giá trị có \(k \in \mathbb{N}\) bộ dữ liệu đầu vào

\(x_1\) \(x_2\) \(\ldots\) \(x_k\) \(Y\)
\(x_{11}\) \(x_{12}\) \(\ldots\) \(x_{1k}\) \(Y_1\)
\(x_{21}\) \(x_{22}\) \(\ldots\) \(x_{2k}\) \(Y_2\)
\(\vdots\) \(\vdots\) \(\ddots\) \(\vdots\) \(\vdots\)
\(x_{n1}\) \(x_{n2}\) \(\ldots\) \(x_{nk}\) \(Y_n\)

Như vậy, theo bảng giá trị, ta có một hệ các phương trình tuyến tính \[Y_i = \sum_{j=1}^n x_{ij} \beta_j + \varepsilon_i,\] trong đó \(\beta = (\beta_1, \beta_2, \ldots, \beta_n)^T\) là các tham số cần tìm; \(x_{ij}\) là quan sát thứ \(i\) của biến thứ \(j\) với \(j = \overline{1, n}\); \(\varepsilon_i\) là sai số ngẫu nhiên của mô hình với các giả thuyết OLS: \[\begin{align*} E(\varepsilon_i) &= 0; \ \forall i. \\ Var(\varepsilon_i) & = \sigma^2; \ \forall i. \\ Cov(\varepsilon_i, \varepsilon_t) &= 0; \ \forall i \ne t. \end{align*}\] Mục tiêu của bài toán là hồi quy là tìm mô hình ước lượng \(\hat{y_i} = \displaystyle \sum_{j=1}^n x_{ij} \hat{\beta}_j\), trong đó các hệ số \(\hat{\beta}_j\); \(j = \overline{1, n}\) được gọi là các hệ số ước lượng của các tham số \(\hat{\beta}_j\); \(j = \overline{1, n}\). Với một mẫu dữ liệu gồm \(k\) quan sát, các hệ số \(\hat{\beta}_j\) được ước tính dựa trên bài toán tìm cực trị (giá trị nhỏ nhất) của hàm phần dư SSR bằng phương pháp OLS (Ordinary Least Squares - Bình phương tối thiểu) sau \[\begin{align*} SSR = \sum_{i=1}^k \left(Y_i - \hat{\beta}_0 - \sum_{j=1}^n \hat{\beta}_jx_{ij}\right). \end{align*}\]

Bằng tính toán ma trận sử dụng vectơ gradient của \(SSR\), ta có thể tìm được \(\hat{\beta}\) để hàm phần dư đạt giá trị nhỏ nhất bằng phương pháp OLS.

Định lý 1.  Ước lượng của các hệ số hồi quy trong mô hình hồi quy tuyến tính đa biến bằng phương pháp OLS là \[\begin{align*} \hat{\beta} = (\textbf{X}^T\textbf{X})^{-1} \textbf{X}^T\textbf{Y}, \end{align*}\] nếu \(rank(X) = n+1\), trong đó \[\textbf{X} = \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1n} \\ 1 & x_{21} & x_{22} & \cdots & x_{2n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{k1} & x_{k2} & \cdots & x_{kn} \end{pmatrix}, \quad \beta = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_n \end{pmatrix}, \quad \textbf{Y} = \begin{pmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_k \end{pmatrix}.\]

Chứng minh Ta có \[\begin{align*} SSR &= (\textbf{Y} - \textbf{X}\hat{\beta})^T(\textbf{Y} - \textbf{X}\hat{\beta}) \\ &= \textbf{Y}^T\textbf{Y} - 2\hat{\beta}^T \textbf{X}^T \textbf{Y} + \hat{\beta}^T \textbf{X}^T \textbf{X}\hat{\beta}. \end{align*}\] Sử dụng vectơ gradient và để ý rằng \(\textbf{X}^T \textbf{X}\) khả nghịch nên ta suy ra \(SSR\) đạt giá trị nhỏ nhất tại \[\hat{\beta} = (\textbf{X}^T \textbf{X})^{-1}\textbf{X}^T \textbf{Y}.\]

Định nghĩa 2.  Ta gọi đại lượng \[\hat{\sigma} = \sqrt{\dfrac{SSR}{k-n-1}}\] là sai số chuẩn phần dư (residual standard error).

Định nghĩa 3.  Đại lượng \[R^2 = \dfrac{SSE}{S_{YY}} = 1 - \dfrac{SSR}{S_{YY}}\] được gọi là hệ số xác định của mô hình hồi quy tuyến tính đa biến. Trong đó \[S_{YY} = SSR + SSE,\] với \[SSE = \sum_{i=1}^k (\hat{\beta}_0 + \hat{\beta}+1 x_{i1} + \cdots + \hat{\beta}_n x_{in} - \overline{Y})^2.\]

2.2   Logistic Regression

Giả sử ứng với mỗi bộ giá trị của các biến giải thích \(x_1, x_2, \ldots, x_k\) ta có biến đầu ra \(Y\) chỉ nhận một trong hai giá trị là \(1\) hoặc \(0\), và ta muốn mô hình hóa xác suất \[P(Y=1)= p.\] Đầu tiên, ta phát biểu khái niệm về hàm sigmoid, tỉ lệ chênh lệch (odds ratio) và logit.

Hàm sigmoid \[S(x) := \dfrac{e^x}{1+e^{x}},\] trong đó \(x\) là một tổ hợp tuyến tính của các biến độc lập, có đồ thị được vẽ dưới dạng hình chữ S (S-shaped), giới hạn giá trị đầu ra \(S(x)\) nằm trong khoảng từ \(0\) đến \(1\). Vì vậy, mô hình hồi quy logistic có thể được thiết lập trong tình huống giả định rằng giá trị của xác suất sao cho \(Y\) nhận giá trị \(1\) tại các biến \((X_1, X_2, \ldots, X_k) = (x_1, x_2, \ldots, x_k)\) là \[\begin{align}\label{logisticfunc} p &= P(Y = 1 \mid X_1 = x_1, \ldots, X_p = x_p) \nonumber \\ &= \dfrac{e^{\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p}}{ 1 + e^{\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p}}. \end{align}\] Hàm trên được gọi là hàm số hồi quy logistic.

Nếu \(p\) là xác suất để một sự kiện diễn ra thì ta gọi tỉ lệ \(\dfrac{p}{1-p}\) là của biến cố đó. Vì \[\begin{align*} 1 - p = P(Y=0 \mid X_1 = x_1, X_2 = x_2, \ldots, X_p = x_p) = \dfrac{1}{1+ e^{\beta_0 + \beta_1x_1 + \cdots + \beta_p x_p}}, \end{align*}\] nên \[\begin{align}\label{oddrat} \dfrac{p}{1-p} = e^{\beta_0 + \beta_1x_1 + \cdots + \beta_p x_p}. \end{align}\] Lấy logarithm tự nhiên hai vế, ta thu được \[\begin{align}\label{logit} \ln \left(\dfrac{p}{1-p}\right) = \beta_0 + \beta_1x_1 + \cdots + \beta_p x_p. \end{align}\] Logarithm tự nhiên của tỉ lệ chênh lệch được gọi là một logit. Vậy theo công thức trên, ta thấy rằng mô hình hồi quy logistic là một mô hình dạng \[\begin{align*} g(E[Y]) = \beta_0 + \beta_1x_1 + \cdots + \beta_p x_p, \end{align*}\] trong đó \(g\) là logit của kỳ vọng \(E[Y]\).

Ta sẽ đánh giá mô hình bằng hàm mất mát (Loss Function). Đây là một hàm số được sử dụng để đo lường mức độ lỗi mà mô hình tạo ra khi dự đoán các kết quả từ dữ liệu đầu vào. Trong bài toán Hồi quy logistic, chúng ta sử dụng hàm mất mát Cross-Entropy (còn gọi là Log Loss) để đánh giá hiệu năng của mô hình. \[\begin{align*} L(\beta) = \sum_{i=1}^n y_i \ln p_i + (1-y_i)\ln (1-p_i), \end{align*}\] trong đó

  1. \(\beta = (\beta_0, \beta_1, \ldots, \beta_k)^T\) là các ước lượng hệ số của mô hình.
  2. \(y_i\) là giá trị thực tế của đầu ra thứ \(i\).
  3. \(p_i\): xác suất dự đoán của biến đầu ra thứ \(i\).

Hàm Cross-Entropy đo lường khoảng cách giữa hai phân phối xác suất \(y_i\) và \(p_i\). Khi mô hình dự đoán chính xác, tức là nếu \(y_i = 1\) thì \(p_i\) càng gần \(1\) và nếu \(y_i = 0\) thì \(p_i\) càng gần \(0\), từ đó hàm mất mát sẽ tiến gần về \(0\).

Bằng phương pháp ước lượng khả năng cực đại, ước lượng của các hệ số hồi quy trong mô hình logistic là nghiệm (nếu có) của hệ phương trình phi tuyến: \[\begin{align*} \textbf{X}^T(\textbf{Y} - \textbf{p}) = \textbf{0}, \end{align*}\] trong đó \[\textbf{X} = \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1k} \\ 1 & x_{21} & x_{22} & \cdots & x_{2k} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{nk} \end{pmatrix}, \quad \beta = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_k \end{pmatrix}, \quad \textbf{Y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}, \quad \textbf{p} = \begin{pmatrix} p_1 \\ p_2 \\ \vdots \\ p_n \end{pmatrix}.\]

Định nghĩa 4.  Ta gọi các đại lượng \[d_i = \mathrm{sign} (y_i - p_i)\sqrt{-2\left(y_1 \ln p_i + (1-y_i)\ln (1-p_i)\right)}, \ \quad i= 1, 2, \ldots, n\] là các phần dư độ lệch (deviance residuals).

Định nghĩa 5.  Ta gọi đại lượng tổng bình phương các phần dư độ lệch \[D = \sum_{i=1}^n d_i^2\] là độ lệch dư (residual deviance).

Định nghĩa 6.  Ta gọi đại lượng \[AIC = -2\mathcal{L}(\hat{\beta}) + 2(n+1)\] là tiêu chuẩn thông tin Akaike (Akaike Information Criterion).

Các mô hình tuyến tính tổng quát nói chung được sử dụng phổ biến hơn vì sự đơn giản trong việc sử dụng, cài đặt và diễn giải kết quả. Tuy nhiên, để có thể sử dụng được mô hình, ta cần phải có những giả định như sai số có phân bố chuẩn, dữ liệu quan hệ tuyến tính để có được những kết quả dự báo hợp lý. Mặt khác, mô hình hồi quy tuyến tính gặp nhiều trở ngại như lỗi dự báo cao hoặc gặp khó khăn nếu như ta gặp phải những dữ liệu phức tạp như:

  1. Có số liệu trống (missing value).
  2. Số liệu không phải dạng số.
  3. Số lượng biến gấp nhiều lần so với số lượng mẫu (hay \(k >>n)\), điều này xảy ra trong trường hợp chúng ta có một số lượng lớn các biến giải thích được cho là quan trọng nhưng lại có ít quan sát.

Ngay cả trong trường hợp \(n>>k\), để giảm số lượng biến giải thích trong mô hình hồi quy, người ta thường phải trải qua các bước ước lượng bình phương tối thiểu, sau đó dựa vào P-value để thực hiện kiểm định giả thuyết sau \[\begin{align*} \begin{cases} H_0: \beta_j = 0, \\ H_1: \beta_j \ne 0. \end{cases} \ \text{ Hoặc } \ \begin{cases} H_0: \beta_m = \beta_{m+1} = \cdots = \beta_{m+p} = 0, \\ H_1: \overline{H}_0. \end{cases} \end{align*}\] Từ đó có thể xác định và loại trừ những biến giải thích không ảnh hưởng đến biến phụ thuộc \(Y\). Sau đó, thử lại mô hình mới với các biến giải thích còn lại, dựa vào các hệ số như \(R^2, AIC,\ldots\) để so sánh và lựa chọn tập hợp con của các biến giải thích mà có ảnh hưởng đến biến phụ thuộc \(Y\).

Nhưng theo các bài báo [1] và [2], việc lựa chọn tập hợp các biến giải thích có ảnh hưởng đến biến phụ thuộc bằng kiểm định dựa trên P-value không còn chính xác, do đó cần phải có một phương pháp hồi quy tốt hơn thay thế cho hồi quy tuyến tính đa biến OLS. Bài tiểu luận này giới thiệu một phương pháp hồi quy hữu dụng mới là phương pháp hồi quy LASSO.

2.3   LASSO Regression

Đầu tiên, ta giới thiệu về một hiện tượng thường gặp trong một số dữ liệu đang muốn phân tích hồi quy, đó là hiện tượng đa cộng tuyến (multicolinearity)

Định nghĩa 7.  Trong một số tập dữ liệu, nếu xảy ra hiện tượng các biến độc lập trong mô hình hồi quy phụ thuộc tuyến tính lẫn nhau thì ta gọi hiện tượng đó là đa cộng tuyến.

Đa cộng tuyến thường xảy ra khi bộ dữ liệu chứa các biến phụ thuộc biểu diễn cùng một khái niệm, chẳng hạn như sức mạnh tay và lực nắm. Hoặc khi ta tạo ra một biến phụ thuộc mới (biến giả tạo) từ biến phụ thuộc đã có như tạo ra biến mới BMI từ biến cân nặng và chiều cao. Trong đó, một nguyên nhân quan trọng dễ dẫn đến hiện tượng này là khi các dòng dữ liệu quá ít so với số lượng các biến giải thích (\(k>>n\)) điều này dễ làm ma trận \(X^TX\) trở nên gần như suy biến.

Một hiện tượng khác trong quá trình xử lý dữ liệu xảy ra là hiện tượng quá khớp (overfitting).

Định nghĩa 8.  Hiện tượng quá khớp là một hiện tượng không mong muốn, xảy ra khi mô hình quá cồng kềnh, phức tạp, hoặc quá nhiều tham số, dẫn đến mô hình học quá kỹ dữ liệu huấn luyện đến mức học luôn cả nhiễu, sai số ngẫu nhiên và đặc điểm không quan trọng. Kết quả là dự đoán trên tập huấn luyện rất chính xác nhưng dự đoán trên dữ liệu mới lại rất kém.

Người ta nhận ra một mô hình bị quá khớp khi training error (\(R^2\) training) rất thấp nhưng test error lại rất cao (\(R^2\) test).

Từ hai vấn đề đã nêu, ta sẽ sử dụng phương pháp hồi quy chính quy hóa LASSO để giải quyết vấn đề \(X^TX\) không khả nghịch và vấn đề quá khớp do mô hình quá phức tạp.

Định nghĩa 9.  LASSO (Least Absolute Shrinkage Selection Sperator) là một phương pháp để ước lượng các tham số của mô hình hồi quy tuyến tính được đề xuất bởi Tibshirani (1996). Trong đó các hệ số \(\hat{\beta}_j, j = \overline{1, n}\) được ước tính dựa trên bài toán tìm cực trị của hàm \[S(\hat{\beta}) = \sum_{i=1}^p \left(y_i - \sum_{j=1}^n x_{ij} \beta_j\right)^2 \ \text{ với điều kiện ràng buộc: } \ \|\hat{\beta}\|_1 \le t,\] trong đó \(\|\hat{\beta}\|_1 = \displaystyle \sum_{j=1}^n |\hat{\beta}_j|\) là chuẩn \(\mathcal{l}_1\) của vectơ \(\hat{\beta}\) và \(t\) là hằng số dương.

Bài toán cực trị có điều kiện tương đương bài toán Lagrange sau: \[\begin{align*} L(\hat{\beta}, \lambda) = \sum_{i=1}^p \left(y_i - \sum_{j=1}^n x_{ij}\hat{\beta}_j)\right)^2 + \lambda \sum_{j=1}^n |\hat{\beta}_j|, \end{align*}\] trong đó \(\lambda\) là nhân tử Lagrange dùng để điều chỉnh mô hình.

Để xác định tham số \(\lambda\), ta chọn \(\lambda\) sao cho trung bình bình phương của training error và validation error nhỏ nhất, nghĩa là \[MSE = \dfrac{\displaystyle \sum_{i=1}^p e_i^2}{n},\] trong đó \(e_i\) là chênh lệch giữa giá trị dự báo và giá trị thực.

Mối quan hệ giữa \(\lambda\) và \(t\) được cho bởi công thức: \[\begin{align*} & \hat{\beta} = (X^TX)^{-1}(X^TY) \\ \Longrightarrow \quad & Y'X(X'X + \lambda I)^{-1} (X'X + \lambda I)^{-1}(X'Y) = t. \end{align*}\]

Khi \(\lambda\) đủ lớn, sẽ có một tham số hồi quy tiến dần về \(0\), do đó chúng không đóng vai trò gì trong mô hình hồi quy, ta có thể loại ra khỏi mô hình đang xét.

Hình dưới là minh hoạ cho mô hình LASSO trong trường hợp hai chiều. Như vậy, tương tự phương pháp OLS, LASSO tìm các ước lượng tham số của một mô hình hồi quy bằng cách cực tiểu hóa hàm mục tiêu là tổng bình phương các phần dư. Tuy nhiên, phương pháp ước lượng LASSO áp đặt ràng buộc tổng giá trị tuyệt đối của các tham số không vượt quá một tham số điều chỉnh.

Dựa vào cơ sở lý thuyết ở trên, ta sẽ vận dụng các mô hình hồi quy đó để khảo sát bộ dữ liệu về rủi ro mắc bệnh tim của các bệnh nhân đựa chọn thông qua đánh giá mối quan hệ giữa việc có mắc bệnh tim hay không mắc bệnh tim và các chẩn đoán lâm sàng.

3   Methods

3.1   Dataset

Bây giờ, ta sẽ tiến hành phân tích bộ dữ liệu về rủi ro mắc bệnh tim và các yếu tố lâm sàng.

Dưới đây là tổng quan về các thông tin của các biến trong bộ dữ liệu.

df <- read.csv("HDP100.csv")

summary(df)
##       Age             Sex             BP         Cholesterol         EKG   
##  Min.   :34.00   Min.   :0.00   Min.   :100.0   Min.   :167.0   Min.   :0  
##  1st Qu.:49.75   1st Qu.:0.00   1st Qu.:120.0   1st Qu.:219.0   1st Qu.:0  
##  Median :55.00   Median :1.00   Median :130.0   Median :254.0   Median :1  
##  Mean   :54.73   Mean   :0.71   Mean   :131.3   Mean   :257.7   Mean   :1  
##  3rd Qu.:61.25   3rd Qu.:1.00   3rd Qu.:140.0   3rd Qu.:283.2   3rd Qu.:2  
##  Max.   :77.00   Max.   :1.00   Max.   :192.0   Max.   :564.0   Max.   :2  
##      Max.HR      ST.depression   Num.of.ves.flu Heart.Disease 
##  Min.   : 95.0   Min.   :0.000   Min.   :0.0    Min.   :0.00  
##  1st Qu.:142.8   1st Qu.:0.000   1st Qu.:0.0    1st Qu.:0.00  
##  Median :156.5   Median :1.000   Median :0.0    Median :0.00  
##  Mean   :152.7   Mean   :1.116   Mean   :0.7    Mean   :0.44  
##  3rd Qu.:167.2   3rd Qu.:1.650   3rd Qu.:1.0    3rd Qu.:1.00  
##  Max.   :195.0   Max.   :5.600   Max.   :3.0    Max.   :1.00
str(df)
## 'data.frame':    100 obs. of  9 variables:
##  $ Age           : int  67 67 68 34 62 51 46 67 50 42 ...
##  $ Sex           : int  0 1 0 1 0 1 1 1 1 1 ...
##  $ BP            : int  115 100 120 118 138 140 150 125 129 120 ...
##  $ Cholesterol   : int  564 299 211 182 294 298 231 254 196 240 ...
##  $ EKG           : int  2 2 2 2 0 0 0 0 0 0 ...
##  $ Max.HR        : int  160 125 115 174 106 122 147 163 163 194 ...
##  $ ST.depression : num  1.6 0.9 1.5 0 1.9 4.2 3.6 0.2 0 0.8 ...
##  $ Num.of.ves.flu: int  0 2 0 0 3 3 0 2 0 0 ...
##  $ Heart.Disease : int  0 1 0 0 1 1 1 1 0 0 ...

Biến phụ thuộc được xét trong bộ dữ liệu này là biến về các ghi nhận các cá nhân có mắc bệnh tim hay không (Heart Disease). Đây là một biến nhị phân, tức chỉ nhận giá trị là 1 (nếu bệnh nhân có mắc bệnh về tim mạch) hoặc 0 (nếu bệnh nhân không mắc bệnh). Trong bộ dữ liệu này có 44 người mắc bệnh tim và 56 người không mắc bệnh tim, do đó có thể nói xác suất để một người mắc bệnh tim là khoảng 0.44, khả năng mắc bệnh hoặc không mắc bệnh có thể coi như là như nhau.

Bộ dữ liệu có 8 biến đầu vào (tức biến giải thích) được xét đến. Đầu tiên là biến về độ tuổi (Age). Tổng quan, những người được xét là những người trưởng thành (tuổi từ 34 đến 77).

Biến thứ hai là biến về giới tính (Sex), đây là một biến nhị phân, trong đó một người được ghi nhận là 1 nếu họ là nam và được ghi nhận là 0 nếu họ là nữ.

Biến thức ba là biến huyết áp (BP), biến này ghi nhận huyết áp tâm thu (đơn vị: mmHg) của bệnh nhân. Một người sẽ được coi là huyết áp bình thường trong khoảng từ 90/60 - 120/80 (mmHg).

Biến thứ tư (Cholesterol) ghi nhận lượng cholesterol của bệnh nhân trong máu (đơn vị: mg/dL), trong đó một người sẽ được coi là bình thường nếu có lượng cholesterol nhỏ hơn 200 mg/dL và nguy cơ cao nếu lớn hơn 240 mg/dL.

Biến thứ năm (EKG) là kết quả điện tâm đồ, nhận ba giá trị 0, 1, 2, các giá trị này đề cập đến độ lớn (biên độ) của sóng và đoạn, nơi 0 mm (hoặc không có sóng/đoạn, đường đẳng điện) là bình thường, 1 mm có thể là dấu hiệu dẹt/âm nhẹ, còn 2 mm (hoặc lớn hơn) thường báo hiệu một vấn đề cần chú ý như ST chênh hoặc T cao nhọn, tùy thuộc vào bối cảnh cụ thể (vị trí chuyển đạo, triệu chứng bệnh nhân).

Biến thứ sáu là biến về nhịp tim cực đại (Max HR) của bệnh nhân (đơn vị: nhịp/phút), nhịp tim này được đo khi bệnh nhân không vận động hay làm gì gắng sức.

Biến thứ bảy đo độ chênh xuống của đoạn ST trong kết quả điện tâm đồ.Thường đoạn ST chênh xuống trên 0.5mm mới có ý nghĩa trong chẩn đoán. Đoạn ST chênh xuống thường gặp trong các bệnh lý: Thiếu máu cục bộ cơ tim, nhồi máu cơ tim thành sau, hạ Kali máu, phì đại thất phải, đau thắt ngực điển hình và bệnh mạch vành…

Biến đầu ra thứ tám (biến cuối cùng) là biến đo số lượng mạch vành chính được nhìn thấy qua kỹ thuật huỳnh quang (fluoroscopy). Trong đó mang các giá trị 0, 1, 2, 3, 4. Nếu bệnh nhân nhận được kết quả là 0, tức là bệnh nhân không có mạch vành lớn nào được nhìn thấy rõ (không được tô màu), từ đó có thể suy ra mạch vành bị tắc nghẽn nặng, thuốc cản quang không vào được, lòng mạch nhỏ hoặc dòng chảy yếu, hoặc dòng máu bất thường. Nếu bệnh nhân nhận được kết quả 1, nghĩa là có 1 mạch vành lớn được nhìn thấy rõ (tô màu). Vậy có một mạch vành thông suốt đủ để thuốc cản quang đi qua và hiện hình. Hai hoặc ba mạch vành còn lại có thể hẹp, tắc hoặc không hiển thị rõ, nguy cơ bệnh tim thấp hơn giá trị 0, nhưng vẫn cao nếu mạch còn lại hẹp nặng.

Tiếp theo, ta sẽ khảo sát mối tương quan tuyến tính giữa các cặp biến dạng số (không phải biến nhị phân) với nhau. Bên dưới là bảng liệt kê các hệ số góc của đường thẳng hồi quy của từng cặp biến thông qua khảo sát mô hình hồi quy tuyến tính đơn biến, quan sát biểu đồ phân tán (scatter plot) và kèm theo đoạn code R tương ứng.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
## 
##     group_rows
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
num_df <- df %>% 
  select(Age, BP, Cholesterol, Max.HR = Max.HR, ST.depression)
library(corrplot)
## corrplot 0.95 loaded
corrplot(cor(num_df, use="complete.obs"),
         method = "color",
         addCoef.col = "black",
         tl.col = "black")

Từ đây, dựa vào hệ số góc của đường thẳng hồi quy, ta thử khảo sát mô hình hồi quy tuyến tính đơn \[\text{Cholesterol} = \beta_0 + \beta_1 \text{Age} + \varepsilon.\]

library(ggplot2)

ggplot(df, aes(x = Age, y = as.numeric(Cholesterol))) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(title = "Liên hệ giữa độ tuổi và huyết áp")
## `geom_smooth()` using formula = 'y ~ x'

  summary(lm(as.numeric(Cholesterol) ~ Age, data = df))
## 
## Call:
## lm(formula = as.numeric(Cholesterol) ~ Age, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -104.357  -38.026   -3.553   25.110  289.707 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 183.5257    33.0638   5.551 2.43e-07 ***
## Age           1.3547     0.5957   2.274   0.0251 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 54.92 on 98 degrees of freedom
## Multiple R-squared:  0.05012,    Adjusted R-squared:  0.04043 
## F-statistic: 5.171 on 1 and 98 DF,  p-value: 0.02514

Dựa trên biểu đồ phân tán này, ta có thể thấy rằng, đường hồi quy có độ dốc dương, nghĩa là tuổi càng cao thì huyết áp có khả năng cao. Vùng xám (confidence interval) hẹp ở giữa, cho thấy rằng mô hình dự đoán tốt hơn. Ở hai biên (tuổi quá thấp hoặc quá cao), vùng này rộng hơn nên độ bất định tăng. Vậy tuổi tác cao cũng có mối liên hệ với việc một người có huyết áp cao và dữ liệu ủng hộ xu hướng này.

Tuy nhiên, ta thấy rằng phần dư của mô hình lệch khá lớn, do đó mô hình này dự đoán không tốt tuổi có ảnh hưởng đến Cholesterol, đại lượng \(R^2\) nhỏ và sai số cao. Tuy nhiên, do p-value nhỏ (vì số lượng mẫu nhiều) nên mô hình này có ý nghĩa thống kê.

Dưới đây ta áp dụng mô hình hồi quy tuyến tính đơn cho hai biến huyết áp và cholesterol và cũng nhận được kết quả tương tự.

## `geom_smooth()` using formula = 'y ~ x'

## 
## Call:
## lm(formula = as.numeric(Cholesterol) ~ BP, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -86.740 -38.492  -2.431  25.702 310.007 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 228.1288    41.2790   5.527  2.7e-07 ***
## BP            0.2249     0.3113   0.722    0.472    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 56.2 on 98 degrees of freedom
## Multiple R-squared:  0.005297,   Adjusted R-squared:  -0.004854 
## F-statistic: 0.5218 on 1 and 98 DF,  p-value: 0.4718

Từ đây, ta sẽ tiến hành chạy các mô hình hồi quy để dự báo khả năng mắc bệnh tim dựa trên các biến đầu vào đã được phân tích ở trên.

3.2   Regression analysis models for the response varible

Ta sẽ chạy mô hình logistic cho bộ dữ liệu trên.

model_logit <- glm(Heart.Disease ~ ., data = df, family = binomial)
summary(model_logit)
## 
## Call:
## glm(formula = Heart.Disease ~ ., family = binomial, data = df)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -6.965201   3.849953  -1.809 0.070425 .  
## Age             0.008795   0.037672   0.233 0.815399    
## Sex             2.895653   0.866967   3.340 0.000838 ***
## BP              0.016416   0.016905   0.971 0.331505    
## Cholesterol     0.004777   0.005300   0.901 0.367439    
## EKG             0.346052   0.285532   1.212 0.225530    
## Max.HR         -0.010632   0.017131  -0.621 0.534824    
## ST.depression   1.158690   0.359020   3.227 0.001249 ** 
## Num.of.ves.flu  0.927022   0.336786   2.753 0.005913 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 137.186  on 99  degrees of freedom
## Residual deviance:  80.201  on 91  degrees of freedom
## AIC: 98.201
## 
## Number of Fisher Scoring iterations: 5

Đồ thị

library(ggplot2)
library(tidyr)
library(dplyr)

df$Heart.Disease <- as.numeric(df$Heart.Disease)

vars <- names(df)[!names(df) %in% c("Heart.Disease")]

for (v in vars) {
  p <- ggplot(df, aes_string(x = v, y = "Heart.Disease")) +
    geom_point(alpha = 0.6) +
    geom_smooth(method = "glm",
                method.args = list(family = "binomial"),
                se = TRUE) +
    labs(title = paste("Hồi quy logistic: Heart.Disease ~", v),
         x = v, y = "Xác suất mắc bệnh tim") +
    theme_minimal(base_size = 14)
  print(p)
}

residuals_std <- rstandard(model_logit)

qqnorm(residuals_std)
qqline(residuals_std, col = "red", lwd = 2)

Tiếp theo, ta sẽ chạy mô hình LASSO cho dữ liệu trên.

library(glmnet)
library(tidyverse)
library(conflicted)
glimpse(df)
## Rows: 100
## Columns: 9
## $ Age            <int> 67, 67, 68, 34, 62, 51, 46, 67, 50, 42, 56, 41, 42, 53,…
## $ Sex            <int> 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1…
## $ BP             <int> 115, 100, 120, 118, 138, 140, 150, 125, 129, 120, 134, …
## $ Cholesterol    <int> 564, 299, 211, 182, 294, 298, 231, 254, 196, 240, 409, …
## $ EKG            <int> 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 0, 2, 0, 2, 0…
## $ Max.HR         <int> 160, 125, 115, 174, 106, 122, 147, 163, 163, 194, 150, …
## $ ST.depression  <dbl> 1.6, 0.9, 1.5, 0.0, 1.9, 4.2, 3.6, 0.2, 0.0, 0.8, 1.9, …
## $ Num.of.ves.flu <int> 0, 2, 0, 0, 3, 3, 0, 2, 0, 0, 2, 0, 0, 3, 1, 1, 3, 0, 1…
## $ Heart.Disease  <dbl> 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1…
df$Heart.Disease <- as.numeric(df$Heart.Disease)
y <- df$Heart.Disease
X <- df%>% select(-Heart.Disease) %>% as.matrix()


lasso_model <- glmnet(X, y, alpha = 1)
# Xem đường cong LASSO
plot(lasso_model, xvar = "lambda", label = TRUE)

# Tìm lambda tối ưu bằng cross-validation
set.seed(123)
cv_lasso <- cv.glmnet(X, y, alpha = 1)

# Vẽ CV-curve
plot(cv_lasso)

# Lambda tối ưu
best_lambda <- cv_lasso$lambda.min
best_lambda
## [1] 0.01565683
coef(cv_lasso, s = "lambda.min")
## 9 x 1 sparse Matrix of class "dgCMatrix"
##                   lambda.min
## (Intercept)    -2.494892e-01
## Age             2.424950e-03
## Sex             3.098400e-01
## BP              7.672145e-04
## Cholesterol     2.122087e-05
## EKG             4.560985e-02
## Max.HR         -5.022959e-04
## ST.depression   1.505977e-01
## Num.of.ves.flu  1.336462e-01
y_pred <- predict(cv_lasso, newx = X, s = "lambda.min")
head(y_pred)
##      lambda.min
## [1,]  0.2649893
## [2,]  0.7371519
## [3,]  0.2713029
## [4,]  0.2410128
## [5,]  0.6468030
## [6,]  1.2699259
rmse <- sqrt(mean((y - y_pred)^2))
rmse
## [1] 0.3732274

4   Results

test <- read.csv("HDP80.csv")
prob_logit <- predict(model_logit, test, type = "response")
pred_logit <- ifelse(prob_logit > 0.5, 1, 0)

table(test$Heart.Disease, pred_logit)
##    pred_logit
##      0  1
##   0 43  4
##   1  8 25
library(Metrics)
logit_acc <- mean(pred_logit == test$Heart.Disease)
logit_acc
## [1] 0.85
x_test <- model.matrix(Heart.Disease ~ ., test)[, -1]
y_test <- test$Heart.Disease
prob_lasso <- predict(lasso_model, newx = x_test, type = "response")
pred_lasso <- ifelse(prob_lasso > 0.5, 1, 0)



lasso_acc <- mean(pred_lasso == y_test)
lasso_acc
## [1] 0.8057377

5   Discussion

6   References