Phân tích hồi qui tuyến tính đơn biến

Phân tích hồi qui tuyến tính đơn biến thực chất là một khai triển từ mô hình phân tích tương quan (correlation analysis) mà tôi đã giải thích trong phần trước. Phân tích tương quan cung cấp cho chúng ta hệ số tương quan (coefficient of correlation), phản ảnh mức độ liên hệ hay tương quan giữa hai biến. Phân tích hồi qui tuyến tính cung cấp cho chúng ta một mô hình để tiên lượng một biến số lâm sàng từ một yếu tố khác. Vì là “mô hình” (model) cho nên phải có tham số (parameter). Do đó, trong phân tích hồi qui tuyến tính, chúng ta còn phải ước tính các tham số của mô hình tiên lượng.

Mô hình hồi qui tuyến tính có lẽ là một trong những phương pháp phân tích thống kê phổ biến nhất, được áp dụng nhiều nhất (và cũng bị lạm dụng nhiều nhất) trong nghiên cứu y học. Mô hình này có một lịch sử khá lâu đời. Năm 1885, nhà khoa học gốc Anh, Francis Galton (một trong những nhà khoa học tiên phong trong di truyền học) giới thiệu khái niệm "regression" (hồi qui) trong một nghiên cứu mà trong đó ông chứng minh rằng chiều cao của những người con không có xu hướng tương quan với chiều cao của cha hay mẹ, mà có xu hướng tương quan với chiều cao trung bình của cha và mẹ. Ông gọi xu hướng này là hồi qui. Nhưng thật ra, Galton không phải là người đầu tiên phát triển, nhưng là người đầu tiên ứng dụng, mô hình hồi qui tuyến tính. Nhà toán học người Pháp thuộc loại hạng “huyền thoại” tên là Adrien Marie Legendre mới là người đầu tiên phát triển và công bố công trình nghiên cứu về hồi qui tuyến tính vào năm 1805 (nhưng lúc đó ông không dùng danh từ "regression"). Nhưng người có ý tưởng nguyên thủy về hồi qui tuyến tính là Carl Friedrich Gauss (một nhà toán học thuộc vào hạng huyền thoại khác), người từng đề cập đến khái niệm hồi qui vào đầu thế kỉ 19.

1. Tóm lược lí thuyết

Mô hình hồi qui tuyến tính (từ nay sẽ viết tắt là HQTT) phát biểu rằng: Gọi là đo lường của đối tượng i (i = 1, 2, 3, …, n) của một biến phụ thuộc, và là đo lường của một biến độc lập cũng của đối tượng i, mối liên hệ tuyến tính giữa hai biến có thể mô tả bằng phương trình với hai thông số a và b như sau:

Trong đó, a và b là hai tham số của mô hình hồi qui tuyến tính cần ước tính từ số liệu quan sát được, và  là phần dư, tức phần không thể tiên lượng bằng đo lường của biến số độc lập. Mô hình trên chỉ hợp lí khi các giả định sau đây đúng:

(i) Giá trị của x không chịu ảnh hưởng sai số đo lường (random error);

(ii)  tuân theo luật phân phối chuẩn với trung bình 0 và phương sai;

(iii)  không có tương quan gì với x; và

(iv) các giá trị nối tiếp của  (chẳng hạn như ) độc lập với nhau.

Với các giả định trên, và bởi vì hai tham số  là bất biến, cho nên, cho một đo lường của x chúng ta có thể ước tính trị số kì vọng (hay nói rõ hơn là số trung bình) của y như sau:

Vấn đề đặt ra là cho một loạt số liệu (x₁,y₁) , (x_2,y₂), . . . , (x_n,y_n), hai tham số  nên ước tính như thế nào. Phương pháp Phương pháp bình phương nhỏ nhất (còn gọi là least squares method) là phương pháp tốt nhất để ước tính hai tham số đó. Theo phương pháp này, chúng ta cần tìm hai ước số a và b (tương đương với  ) sao cho tổng số bình phương giữa giá trị quan sát (y_i) và giá trị tiên đoán () là thấp nhất, nói cách khác, chúng ta tối thiểu hóa:

Hóa ra, muốn tối thiểu hóa Q chúng ta chỉ cần giải hệ phương trình đơn giản sau đây:

Chú ý rằng, trong các phương trình trên,  là số trung bình của biến số x và y. Xin nhắc lại rằng chúng ta không biết được giá trị của a và b, mà chỉ có thể ước tính chúng, và ước số của hai tham số này chính là  . Thuật ngữ thống kê gọi a là intercept, và b là gradient hay slope. Như chúng ta thấy qua phương trình trên, intercept chính là giá trị của y khi x = 0.

Công thức [4] cho thấy ước số b chỉ đơn giản bằng hiệp biến của x và y chia cho phương sai của y. Tuy các công thức này mới nhìn qua có vẻ rắc rối, nhưng trong thực tế thì rất đơn giản, bạn đọc chỉ cần một máy tính cầm tay (calculator) hay tốt hơn nữa phần mềm Excel cũng có thể tính rất dễ dàng.

Ví dụ 1 (tiếp tục) – cân nặng và vòng eo: Trong phần trước (phân tích tương quan), chúng ta có số liệu về cân nặng và vòng eo của 15 đối tượng như sau (in lại để dễ theo dõi):

Bảng 1. Cân nặng và vòng eo của 15 đối tượng người Việt

Trọng lượng (weight; kg)	Vòng eo (waist; cm)
51.0	71.0
66.0	89.0
47.0	64.0
54.0	74.0
64.0	87.0
75.0	93.0
54.0	66.0
52.0	74.0
53.0	75.0
52.0	72.0
48.0	70.0
46.0	66.0
63.0	81.0
40.0	57.0
90.0	94.0

Gọi cân nặng là x và vòng eo là y. Với cách gọi này, chúng ta có ý muốn sử dụng cân nặng của để tiên đoán vòng eo của một đối tượng. Xin nhắc lại, trong bài trước, chúng ta đã có những kết quả tính toán sau đây:

Với n = 15 đối tượng, và dựa vào công thức [4] và [5] chúng ta có thể ước tính tham số b và a như sau:

Đến đây thì chúng ta đã có một công thức để tiên đoán vòng eo dựa vào trọng lượng của một đối tượng qua phương trình sau đây:

Chú ý rằng trong phương trình trên, biến số y có dấu mũ trên để nhắc nhở rằng đây là giá trị tiên đoán, để phân biệt với giá trị đo lường (thực tế) là y_i.

Ý nghĩa của phương trình này là gì? Ở đây, giá trị a = 30 không có ý nghĩa thực tế gì đáng kể, nhưng b = 0.80 có nghĩa là mỗi kg cân nặng tương quan với 0.8 cm vòng eo. Chẳng hạn như nếu chúng ta hỏi biết một bệnh nhân cân nặng là 60 kg, thì qua phương trình trên, có thể tiên đoán rằng vòng eo của bệnh nhân là:

nhưng nếu một bệnh nhân khác với cân nặng 62 kg, thì chúng ta có thể tiên đoán vòng eo của bệnh nhân là:

Bây giờ, chúng ta thử xem phương trình trên tiên đoán vòng eo chính xác ra sao, bằng cách sử dụng phương trình trên để ước tính vòng eo cho từng đối tượng trong Bảng 1 như sau:

Bảng 2. Tiên đoán vòng eo dựa vào cân nặng cho 15 đối tượng

Số ID (i)	Trọng lượng (xi)	Vòng eo thực tế (yi)	Vòng eo tiên đoán qua mô hình ()	Khác biệt giữa giá trị thực tế và tiên đoán
1	51.0	71.0	70.8	0.2
2	66.0	89.0	82.8	6.2
3	47.0	64.0	67.6	-3.6
4	54.0	74.0	73.2	0.8
5	64.0	87.0	81.2	5.8
6	75.0	93.0	90.0	3.0
7	54.0	66.0	73.2	-7.2
8	52.0	74.0	71.6	2.4
9	53.0	75.0	72.4	2.6
10	52.0	72.0	71.6	0.4
11	48.0	70.0	68.4	1.6
12	46.0	66.0	66.8	-0.8
13	63.0	81.0	80.4	0.6
14	40.0	57.0	62.0	-5.0
15	90.0	94.0	102.0	-8.0

Cột số 4 cho thấy giá trị vòng eo được tiên đoán qua phương trình , và cột số 5 cung cấp cho chúng ta độ chính xác của tiên đoán (lấy vòng eo đo lường trừ cho vòng eo tiên đoán). Như có thể thấy qua bảng này, phương trình tiên đoán khá chính xác vòng eo. Chẳng hạn như với đối tượng 1, vòng eo thực tế là 71 cm, nhưng mô hình tiên đoán là 70.8 cm. Điều này có nghĩa gì? Nó có nghĩa là giá trị tiên đoán 70.8 cm chính là vòng eo trung bình của tất cả những người với cân nặng 51 kg.

Tuy nhiên nhìn qua bảng trên, chúng ta cũng thấy ở một số đối tượng, phương trình tiên đoán không mấy tốt. Chẳng hạn như đối tượng 2, phương trình tiên đoán thấp hơn thực tế đến 6.2 cm, nhưng với đối tượng 15, phương trình tiên đoán cao hơn thực tế đến 8 cm! Chúng ta sẽ quay lại để thẩm định chất lượng và độ chính xác của mô hình này trong một phần tiếp theo.

Ngoài ra, cách tốt hơn là vẽ biểu đồ so sánh giữa giá trị đo lường vòng eo và giá trị tiên đoán như Biểu đồ 1 dưới đây thể hiện. Biểu đồ này cho thấy, mô hình trên tiên đoán khá chính xác vòng eo ở những người có trọng lượng dưới 60 kg, nhưng ở những người có trọng lượng cao hơn ngưỡng này, thì giá trị tiên đoán không mấy chính xác.

2. Kiểm định giả thuyết về a và b

Quay lại với mô hình  , chúng ta có thể rút ra vài nhận xét như sau: Nếu, thì phương trình đơn giản thành , tức là không có mới tương quan nào giữa x và y; nhưng nếu  (tức có thể âm hay dương) thì mối liên hệ giữa x và y hiện hữu. Do đó, kiểm định mô hình hồi qui tuyến tính tập trung vào kiểm định giá thuyết . Để kiểm định giả thuyết này, chúng ta cần tính toán phương sai của b (vì nên nhớ rằng b là ước số của ).

Chúng ta biết rằng  là số trung bình của vòng eo, và phương sai của vòng eo, trước khi biết cân nặng, có thể ước tính như sau:

Công thức trên còn được gọi là phương sai vô điều kiện, Nhưng trong mô hình [1], như vừa đề cập trên,  chính là vòng eo trung bình với điều kiện x_i. Chính vì thế mà phương sai của y (kí hiệu s²) với điều kiện biết cân nặng được ước tính bằng cách thay thế :

(Chú ý rằng mẫu số là n – 2, chứ không phải n – 1, vì đây là phương sai được ước tính với 2 tham số a và b, cho nên nphải trừ cho 2).

Gọi e_i là độ khác biệt giữa vòng eo thực tế và vòng eo tiên đoán:

Trong thuật ngữ thống kê học, e_i còn được gọi là “residual”. Phương sai trong phương trình [6] có thể viết lại như sau:

s² chính là ước số của  trong mô hình [1]

Sau vài thao tác đại số, có thể chứng minh rằng phương sai của b và a có thể viết như sau:

và:

Khi n (số đối tượng tương đối lớn), b tuân theo luật phân phối chuẩn với số trung bình là  và phương sai như trình bày trong [8]. Do đó, kiểm định giả thuyết = 0 có thể dựa vào tỉ số T_b sau đây:

(Chú ý, s_x chính là độ lệch chuẩn của x). Nếu  = 0 thì T_b tuân theo luật phân phối t với bậc tự do là n – 2.

Ngoài ra, kiểm định giả thuyết  = 0 cũng có thể tính toán qua tỉ số T_a như sau:

Nếu  = 0 thì t tuân theo luật phân phối t với bậc tự do là n – 2.

Ví dụ 1 (tiếp tục): Chúng ta tiếp tục sử dụng số liệu của ví dụ 1 để minh họa cho các tính toán trên. Để ước tính phương sai s² theo công thức [6], chúng ta cần tính toán bình phương của khác biệt giữa y_i và  và tổng số như sau:

Bảng 3. Tính toán phương sai của y

Số ID (i)	Trọng lượng (xi )	Vòng eo thực tế ( yi)	Vòng eo tiên đoán qua mô hình ()	Khác biệt giữa giá trị thực tế và tiên đoán
1	51.0	71.0	70.8	0.2	0.04
2	66.0	89.0	82.8	6.2	38.44
3	47.0	64.0	67.6	-3.6	12.96
4	54.0	74.0	73.2	0.8	0.64
5	64.0	87.0	81.2	5.8	33.64
6	75.0	93.0	90.0	3.0	9.00
7	54.0	66.0	73.2	-7.2	51.84
8	52.0	74.0	71.6	2.4	5.76
9	53.0	75.0	72.4	2.6	6.76
10	52.0	72.0	71.6	0.4	0.16
11	48.0	70.0	68.4	1.6	2.56
12	46.0	66.0	66.8	-0.8	0.64
13	63.0	81.0	80.4	0.6	0.36
14	40.0	57.0	62.0	-5.0	25.00
15	90.0	94.0	102.0	-8.0	64.00
	= 57.0, sx= 12.8	=75.5 sy=11.1			Tổng cộng: 251.7

Bảng trên cho thấy s²= 251.7 / 13 = 19.36. Do đó, phương sai của b là:

và độ lệch chuẩn của b là: . Do đó, kiểm định giả thuyết  = 0 là:

Nếu  = 0 thì T_b dao động từ -2 đến +2, nhưng ở đây, chúng ta thấy T_b = 8.69, tức cao hơn 4 lần so với giá trị kì vọng, cho nên chúng ta có thể kết luận rằng mối liên hệ giữa vòng eo và cân nặng có ý nghĩa thống kê.

3. Phân tích phương sai

Một trong những mục đích của phân tích hồi qui tuyến tính là tìm hiểu xem biến độc lập có thể giải thích bao nhiêu phần trăm độ biến thiên của biến phụ thuộc. Trong ví dụ này, cụ thể là chúng ta muốn biết bao nhiêu phần trăm của độ biến thiên (hay khác biệt) giữa các cá nhân về vòng eo có thể giải thích bằng cân nặng. Cụm từ căn bản ở đây là “biến thiên” (thuật ngữ thống kê học là variation). Chú ý rằng mỗi cá nhân có 3 giá trị: vòng eo thực tế y_i, vòng eo tiên đoán bằng cân nặng  , và vòng eo trung bình của quần thể . Mối liên hệ giữa ba giá trị này có thể mô tả như sau:

Nói cách khác, độ khác biệt giữa vòng eo của một cá nhân và số trung bình là tổng số khác biệt của: (a) giữa giá trị tiên đoán và số trung bình (), và (b) giữa giá trị thực tế và giá trị tiên đoán .

Do đó, một chỉ số có thể đo lường độ biến thiên của một biến là tổng bình phương của biến đó. Nói cách khác, nếu y_i là vòng eo của từng cá nhân i và  là vòng eo trung bình, thì tổng bình phương là (kí hiệu SST):

Trong mô hình [1], chúng ta có là giá trị tiên đoán của y_i sau khi “điều chỉnh” cho cân nặng (x_i), cho nên, cùng một logic như trên, có thể nói rằng tổng biến thiên của vòng eo mà mô hình [1] có thể giải thích được là (kí hiệu SSR):

Và phần còn lại, tức phần biến thiên của vòng eo không thể giải thích bằng mô hình [1] là (kí hiệu SSE):

Do đó, qua mối liên hệ [9], chúng ta có thể chứng minh rằng:

SST = SSR + SSE

hay:

Chúng ta thấy để ước tính SST, chúng ta phải “tiêu” ra một tham số (số trung bình), cho nên bậc tự do (degrees of freedom) của SST là n – 1. Do đó, số trung bình bình phương (mean squares) là MST = SST / (n – 1). Để tính SSE, chúng ta phải cần đến hai tham số (a và b), cho nên bậc tự do của SSE là n – 2; do đó, số trung bình bình phương là MSE = SSE/ (n – 2). Các chỉ số này có thể tóm lược trong một bảng phân tích phương sai (analysis of variance) như sau:

Bảng 4. Phân tích phương sai cho mô hình hồi qui tuyến tính

Nguồn biến thiên	Bậc tự do (degrees of freedom)	Tổng bình phương (sum of squares)	Trung bình bình phương (mean squares)
Hồi qui (regression)	1		MSR = SSR / 1
Phần dư (residual)	n – 2		MSE = SSE / (n – 2)
Tổng số biến thiên	n – 1

Kiểm định F

Để kiểm định ý nghĩa thống kê của giả thuyết  = 0, chúng ta đã làm quen với tỉ số T_b. Nhưng còn một kiểm định tương đương khác là kiểm định F, với công thức sau đây:

Thật ra, F như định nghĩa trong công thức [15] chính là bình phương của T_b.

Hệ số xác định (coefficient of determination)

Bởi vì tổng số biến thiên là SST, và trong số này, biến thiên có thể giải thích qua mô hình hồi qui tuyến tính [1] là SSR, cho nên chúng ta có thể ước tính số phần trăm mà mô hình có thể giải thích tổng biến thiên của y. Hệ số này được gọi là hệ số xác định và kí hiệu là R²:

Nhìn qua mối liên hệ [13], chúng ta dễ dàng thấy R² có giá trị từ 0 đến 1. Nếu R² = 0, mô hình hồi qui tuyến tính coi như vô dụng, vì không giải thích phần trăm nào biến thiên của y. Nếu R²= 1 hay gần 1, mô hình hồi qui tuyến tính có thể tiên đoán chính xác giá trị của y.

Tuy nhiên cần phải nhấn mạnh rằng một mô hình với R² cao không có nghĩa là rằng mô hình tốt. Thật vậy, R² có thể cao nếu b cao hay biến độc lập có range (dãy số) lớn. Ngoài ra, R² có thể cao khi mô hình tuyến tính được áp dụng cho một mối liên hệ phi tuyến tính.

Ví dụ 1 (tiếp tục): Chúng ta tiếp tục sử dụng số liệu của ví dụ 1 để minh họa cho các tính toán trên.

Bảng 5. Chi tiết phân tích phương sai của mô hình hồi qui tuyến tính

Số ID (i)	yi
1	71.0	70.8	20.25	0.04
2	89.0	82.8	182.25	38.44
3	64.0	67.6	132.25	12.96
4	74.0	73.2	2.25	0.64
5	87.0	81.2	132.25	33.64
6	93.0	90.0	306.25	9.00
7	66.0	73.2	90.25	51.84
8	74.0	71.6	2.25	5.76
9	75.0	72.4	0.25	6.76
10	72.0	71.6	12.25	0.16
11	70.0	68.4	30.25	2.56
12	66.0	66.8	90.25	0.64
13	81.0	80.4	30.25	0.36
14	57.0	62.0	342.25	25.00
15	94.0	102.0	342.25	64.00
Tổng cộng:			1715.75	251.7

Qua tính toán trình bày trong bảng trên, chúng ta có SST = 1715.75, và SSE = 251.7. Do đó, SSR = SST – SSE =1715.75 – 251.7 = 1463.95. Bảng 6 sau đây tóm lược các chỉ số trên:

Bảng 6. Phân tích phương sai (ví dụ 1)

Nguồn biến thiên	Bậc tự do (degrees of freedom)	Tổng bình phương (sum of squares)	Trung bình bình phương (mean squares)
Hồi qui (regression)	1	SSR = 1463.95	1463.95
Phần dư (residual)	13	SSE = 251.7	19.37
Tổng số biến thiên	14	SST = 1715.75

Hệ số xác định bội, theo công thức [13], có thể ước tính là:  . Nói cách khác, mô hình tuyến tính với cân nặng là biến độc lập có thể giải thích khoảng 85% tổng biến thiên (hay những khác biệt về) vòng eo giữa các cá nhân.

4. Phân tích phần dư và kiểm tra giả định

Như đề cập trên, phần dư (residual) là độ khác biệt giữa giá trị thực tế và giá trị tiên đoán của biếnphụ thuộc: . Phần dư rất quan trọng cho việc thẩm định tính hợp lí và độ chính xác của mô hình tiên đoán. Xin nhắc lại rằng mô hình hồi qui tuyến tính mà chúng ta áp dụng trong Ví dụ 1 chỉ có giá trị khoa học nếu mô hình này đáp ứng những giả định sau đây:

1. Mối liên hệ giữa cân nặng và vòng eo (biến độc lập và biến phụ thuộc, hay x và y) phải là mối liên hệ tuyến tính, tức tuân thủ theo một đường thẳng;

2. Phương sai của y không thay đổi tùy theo giá trị của x; hay nói cách khác, phần dư e_i không biến chuyển một cách có hệ thống với x_i;

3. Biến phụ thuộc y (hay e_i) tuân theo luật phân phối chuẩn; hay một cách tương đương;

4. Các giá trị của e_i không liên quan nhau.

Để kiểm định những giả định trên, một số biểu đồ sau đây có thể áp dụng:

o                                Biểu đồ e_i và  để tìm xem có các giá trị “outlier” (tức những giá trị mà mô hình không tiên đoán chính xác), và xem mô hình có vi phạm giả định 1 và 2 hay không.

o                                Biểu đồ so sánh giá trị quan sát của e_i và giá trị kì vọng (dựa vào luật phân phối chuẩn) của e_i để kiểm tra xem giả định 3 có đáp ứng hay không.

o                                Biểu đồ e_i và x để kiểm tra xem có cần hoán chuyển x hay không.

Tuy nhiên, phương sai của e_i không cố định. Do đó, các nhà nghiên cứu khuyến cáo chúng ta nên sử dụng phần dư chuẩn hóa (thuật ngữ thống kê là standardised residuals). Phần dư chuẩn hóa (kí hiệu r_i) được định nghĩa như sau: lấy phần dư e_i chia cho độ lệch chuẩn của mô hình:

Với cách hoán chuyển này, r_i sẽ có giá trị trung bình là 0 và phương sai bằng 1. Chúng ta có thể sử dụng r_i để kiểm tra các giả định của mô hình hồi qui tuyến tính.

(Còn tiếp)

Chú thích kĩ thuật:

Các mã R sau đây đã được sử dụng cho phân tích vừa trình bày.

# Mô phỏng cho biểu đồ 1d

zn1 <- rnorm(100)

zn2 <- rnorm(100)

r = 0.80

# mô phỏng x với trung bình 5 và độ lệch chuẩn 1.5

x <- 5 + 1.5*zn1

# mô phỏng y với trung bình 10 và độ lệch chuẩn 2 và r=0.80

y <- 10 + r*2.0*zn1 + 2.0*sqrt(1-r^2)*zn2

# vẽ biểu đồ

plot(y ~ x, pch=16)

# Mô phỏng cho biểu đồ 1e

zn1 <- rnorm(100)

zn2 <- rnorm(100)

r = -0.80

# mô phỏng x với trung bình 5 và độ lệch chuẩn 1.5

x <- 5 + 1.5*zn1

# mô phỏng y với trung bình 10 và độ lệch chuẩn 2 và r=0.80

y <- 10 + r*2.0*zn1 + 2.0*sqrt(1-r^2)*zn2

# vẽ biểu đồ

plot(y ~ x, pch=16)

# Mô phỏng cho biểu đồ 1f

zn1 <- rnorm(100)

zn2 <- rnorm(100)

r = 0.001

# mô phỏng x với trung bình 5 và độ lệch chuẩn 1.5

x <- 5 + 1.5*zn1

# mô phỏng y với trung bình 10 và độ lệch chuẩn 2 và r=0.001

y <- 10 + r*2.0*zn1 + 2.0*sqrt(1-r^2)*zn2

# vẽ biểu đồ

plot(y ~ x, pch=16)

# nhập số liệu trọng lượng và vòng eo

weight <- c(51,66,47,54,64,75,54,52,53,52,48,46,63,40,90)

waist <- c(71,89,64,74,87,93,66,74,75,72,70,66,81,57,94)

# trung bình

mean(weight)

mean(waist)

# độ lệch chuẩn

sd(weight)

sd(waist)

# hiệp biến và r

cov(weight, waist)

cov(weight, waist) / (sd(weight)*sd(waist))

# hay sử dụng hàm cor trong R

cor(weight, waist)

# vẽ biểu đồ hồi qui tuyến tính

model <- lm(waist ~ weight)

plot(waist ~ weight, pch=16, xlab="Weight", ylab="Waist circumference")

abline(model)

par(mfrow=c(2,2))

plot(model)

predicted <- 29.9582+0.7996*weight

residuals <- waist-predicted

plot(residuals~weight, pch=16)

height <- c(169,153,153,166,165,176,154,151,149,152,152,147,167,150,165)

hip <- c(87,98,85,84,100,98,90,94,85,88,92,88,93,67,106)