Hướng dẫn cách chạy hồi quy tuyến tính kiểu R trong Python một cách dễ dàng

Python rất phổ biến trong phân tích thống kê nhờ vào số lượng thư viện khổng lồ. Một trong những phép tính thống kê phổ biến nhất là hồi quy tuyến tính (linear regression). Thư viện statsmodels cung cấp một số công cụ mạnh mẽ cho hồi quy và phân tích phương sai. Dưới đây là cách bắt đầu với các mô hình tuyến tính.

Xem thêm: Hướng dẫn cách chỉnh sửa biến môi trường (Environment Variables) trên Windows 10 hoặc 11

1. Thư viện statsmodels là gì?

statsmodels là một thư viện Python dùng để chạy các phép kiểm thống kê phổ biến. Nó đặc biệt hướng đến phân tích hồi quy, nhất là loại phân tích bạn thường thấy trong kinh tế lượng, nhưng bạn không cần phải là nhà kinh tế học để sử dụng nó.

Thư viện này có một đường cong học tập nhất định, nhưng một khi đã quen, bạn sẽ thấy nó linh hoạt hơn nhiều so với các hàm hồi quy mà bạn tìm thấy trong các chương trình bảng tính như Excel. Tuy nhiên, nó sẽ không tự vẽ biểu đồ cho bạn. Nếu bạn muốn tạo biểu đồ phân tán cổ điển với đường hồi quy được vẽ lên trên, bạn nên sử dụng một thư viện như Seaborn.

Một lợi thế của việc sử dụng statsmodels là nó được kiểm tra chéo về độ chính xác với các gói phần mềm thống kê khác như R, Stata và SAS. Vì vậy, đây có thể là gói phù hợp với bạn nếu bạn đang làm nghiên cứu chuyên nghiệp hoặc học thuật.

2. Hồi quy tuyến tính đơn giản (Simple Linear Regression)

Nếu bạn chỉ muốn xác định mối quan hệ của một biến phụ thuộc (y), hay còn gọi là biến nội sinh (endogenous variable) trong ngôn ngữ kinh tế lượng và statsmodels, với biến ngoại sinh (x), hay biến độc lập (independent variable), bạn có thể dễ dàng làm điều này với statsmodels.

Nếu bạn có kinh nghiệm với R hoặc muốn một cách nhanh chóng để tạo một hồi quy với statsmodels bằng cách sử dụng một Pandas DataFrame, bạn có thể sử dụng các công thức kiểu R.

Đầu tiên, bạn cần import statsmodels và API công thức của nó:

import statsmodels.formula.api as smf
import statsmodels.api as sm

Tôi sẽ sử dụng mô hình quen thuộc của mình, thực chất là từ Seaborn, về tiền tip và tổng hóa đơn trong một nhà hàng ở New York City.

import seaborn as sns
sns.set_theme()
tips = sns.load_dataset('tips')

Tôi có thể xem qua tập dữ liệu với phương thức head:

tips.head()

Tôi muốn tìm hiểu mối quan hệ giữa tiền tip so với tổng hóa đơn. Tôi có thể vẽ mối quan hệ này bằng biểu đồ phân tán trong Seaborn:

sns.relplot(x='total_bill',y='tip',data=tips)

Trông có vẻ có một mối quan hệ tuyến tính dương. Tôi có thể thử một biểu đồ hồi quy:

sns.regplot(x='total_bill',y='tip',data=tips)

Quả thực có một mối quan hệ tuyến tính dương, như được chỉ ra bởi đường thẳng. Với hóa đơn lớn hơn, tiền tip dường như cũng tăng lên. Biểu đồ này sẽ không cung cấp cho tôi các hệ số mà tôi cần cho đường thẳng này. statsmodels sẽ làm điều đó.

Đầu tiên, tôi sẽ fit (khớp) đường thẳng.

results = smf.ols('tip ~ total_bill',data = tips).fit()

Điều này sử dụng API công thức của statsmodels. Nó sử dụng module patsy, mượn các quy ước được thiết lập bởi ngôn ngữ R. “ols” là viết tắt của “Ordinary Least Squares” (Bình phương tối thiểu thông thường), phương pháp được sử dụng để tạo hồi quy. Phương pháp này giảm thiểu bình phương khoảng cách từ đường thẳng đến các điểm dữ liệu, còn gọi là phần dư (residual), nhiều nhất có thể.

Lệnh trên thiết lập cột tip là biến nội sinh, và cột total_bill là biến ngoại sinh. Dấu ngã ($~$) trong ngữ cảnh này giống như một dấu bằng cách điệu. Nó trông lạ lùng, nhưng là một cách ngắn gọn để chỉ định mối quan hệ. Tham số data = tips cho statsmodels biết sử dụng DataFrame tips. Phương thức .fit() khớp các điểm dữ liệu vào một đối tượng kết quả, trong trường hợp này được gọi là “results.”

Bạn có thể xem tóm tắt theo một trong hai cách. Bạn có thể sử dụng lệnh print của Python nếu bạn đang viết một script Python:

print(results.summary())

Nếu bạn đang sử dụng một phiên tương tác, chẳng hạn như IPython hoặc Jupyter notebook, bạn có thể chỉ cần gõ tên đối tượng:

results.summary()

Điều này sẽ hiển thị kết quả của hồi quy.

Bạn cũng có thể sử dụng mảng NumPy để thực hiện hồi quy. Bạn sẽ tạo một ma trận “thiết kế” của biến độc lập cộng với một cột độ chặn (intercept) gồm các giá trị 1.

Hãy tạo các giá trị $x$ và $y$ bằng cách tạo các số ngẫu nhiên với NumPy:

x = np.linspace(-10,10,100)
y = 2*x - 3

statsmodels bao gồm một hàm để tạo độ chặn này với phương thức add_intercept. Điều này sẽ tạo ma trận thiết kế.

X = sm.add_intercept(x)

Việc tạo mô hình tương tự như sử dụng công thức kiểu R.

model = sm.OLS(y,X)
results = model.fit()
print(results.summary())

Lưu ý rằng khi sử dụng API statsmodels thay vì API công thức, hàm OLS được viết hoa toàn bộ.

3. Hồi quy tuyến tính đa biến (Multiple Linear Regression)

Thật dễ dàng để mở rộng hồi quy đơn giản sang nhiều hơn một biến ngoại sinh (multiple variables). Điều này hữu ích để xem liệu có thể có một biến gây nhiễu (confounding variable) dẫn đến một tương quan sai lầm hay không. Thay vì khớp một đường thẳng trên một tập hợp các điểm dữ liệu trên mặt phẳng, bạn đang khớp một mặt phẳng trên các điểm dữ liệu trong không gian ba chiều, và một siêu mặt phẳng (hyperplane) với nhiều hơn hai biến ngoại sinh và do đó nhiều hơn ba chiều.

Điều này dễ dàng đạt được với các công thức kiểu R. Giả sử chúng ta muốn thêm kích thước bữa tiệc (party size) vào hồi quy tiền tip của chúng ta. Chúng ta có thể làm điều đó:

results = smf.ols('tip ~ total_bill + size',data = tips).fit()
results.summary()

Để thêm các biến hồi quy (regressors), bạn chỉ cần thêm chúng vào công thức bằng dấu cộng (“+”).

Bạn cũng có thể sử dụng điều này để mô hình hóa các mối quan hệ phi tuyến tính, chẳng hạn như phương trình bậc hai. Hãy tạo một mô hình khác bằng cách sử dụng NumPy và tạo một pandas DataFrame:

x = np.linspace(-10,10,100)
y = 3*x**2 + 2*x + 5
df = pd.DataFrame({'x':x,'y':y})

Chúng ta có thể hình dung điều này bằng một biểu đồ phân tán khác:

sns.relplot(x='x',y='y',data=df)

Biểu đồ dường như cho thấy một parabol bậc hai cổ điển.

Chúng ta có thể tạo một công thức khác để mô hình hóa điều này với statsmodels bằng cách thêm bình phương các giá trị $x$ làm biến hồi quy:

results = smf.ols('y ~ x + I(x**2)',data=df).fit()
results.summary()

Dấu “I()” xung quanh $x^2$ là để báo cho statsmodels biết rằng đây không phải là một biến hồi quy riêng biệt, mà là một phép toán trên cột $x$ trong pandas DataFrame.

4. Giải thích kết quả hồi quy (Interpreting Regression Results)

Khi bạn xem thông tin results.summary(), bạn sẽ thấy một bảng với một số giá trị. Hãy cùng xem qua nó.

R-squared (Hệ số xác định) ở góc trên bên phải là bình phương của hệ số tương quan ($r$). Nó cho bạn biết mức độ tương quan giữa các biến độc lập hoặc ngoại sinh với các biến phụ thuộc hoặc nội sinh, và do đó, mức độ phù hợp của mô hình với dữ liệu. $R^2$ càng gần 1, thì càng có nhiều tương quan giữa các biến. Dữ liệu thực tế hiếm khi tương quan hoàn hảo.

Adjusted R-squared (Hệ số xác định điều chỉnh) hữu ích cho hồi quy đa biến và điều chỉnh cho hồi quy sai lầm, cho bạn một hệ số tương quan chính xác hơn.

Bảng bên dưới hiển thị độ chặn (intercept), các hệ số (coefficients) của mỗi biến hồi quy, và mức độ phù hợp của chúng. Với hồi quy tuyến tính đơn giản, phương trình khớp với dạng $y = mx + b$. Giá trị của hệ số biến độc lập sẽ là $m$, và độ chặn là $b$. Điều này tự nhiên mở rộng cho nhiều biến hồi quy hơn.

Các cột bên phải của bảng hiển thị khoảng tin cậy nơi các giá trị thực có thể nằm.

Cột “std err” (sai số chuẩn) đo lường sai số chuẩn, là khoảng cách từ các điểm dữ liệu đến đường thẳng được khớp. Số này càng thấp, mô hình càng phù hợp.

Các p-value (giá trị $p$) xác định ý nghĩa thống kê của các biến. statsmodels thực hiện kiểm định $t$ trên mỗi biến của giả thuyết không (null hypothesis) rằng độ dốc là 0. Nếu giá trị $p$ thấp hơn ngưỡng xác định trước (thường là tối thiểu $0.05$), kết quả là có ý nghĩa.

Phần dưới cùng của biểu đồ hiển thị phân phối của phần dư (residuals), hoặc sự khác biệt giữa các điểm dữ liệu và đường thẳng. Nó phải càng gần với phân phối chuẩn càng tốt để có sự phù hợp tốt. Skew (độ xiên) và Kurtosis (độ nhọn) đo lường mức độ phần dư khác biệt so với phân phối chuẩn.

5. Phân tích phương sai một chiều (One-way Analysis of Variance – ANOVA)

Bạn có thể tự hỏi làm thế nào để so sánh một giá trị số trên một biến phân loại. Để làm điều đó, bạn cần chuyển sang phân tích phương sai, hay ANOVA. Thật dễ dàng đạt được điều này với statsmodels. Vì chúng ta chỉ sử dụng một loại, nên điều này được gọi là ANOVA một chiều.

Hãy tải một tập dữ liệu khác, lần này là về chim cánh cụt ở Nam Cực:

penguins = sns.load_dataset('penguins')
penguins.head()

Hãy xem liệu loài (species) có phải là một yếu tố dự đoán đáng kể về chiều dài mỏ (bill length) hay không. Chúng ta có thể sử dụng một mô hình tuyến tính tương tự như hồi quy, nhưng với một biến phân loại:

penguin_lm = smf.ols('bill_length_mm ~ species',data=penguins).fit()

Chúng ta có thể đưa điều này vào hàm anova_lm của statsmodels:

results = sm.stats.anova_lm(penguin_lm)

Chúng ta có thể xem kết quả:

print(results)

Chúng ta sẽ thấy một bảng kết quả, nhưng con số cần chú ý là giá trị $p$. Vì nó rất thấp (ở dạng ký hiệu khoa học âm), nên loài là một yếu tố dự đoán đáng kể về chiều dài mỏ.

6. Phân tích phương sai đa chiều (Multi-way Analysis of Variance)

Chúng ta có thể mở rộng điều này với một biến phân loại khác bằng ANOVA hai chiều. Thật dễ dàng để sửa đổi công thức ban đầu của chúng ta và thêm đảo (island) vào hỗn hợp để xem liệu điều đó cũng có thể dự đoán chiều dài mỏ hay không:

penguin_multi_lm = smf.ols('bill_length_mm ~ species * island',data=penguins).fit()
results = sm.stats.anova_lm(penguin_multi_lm)
print(results)

Điều này tương tự như hồi quy đa biến, nhưng chúng ta sử dụng phép nhân ($*$) cho công thức thay vì phép cộng.

Hồi quy đa biến và ANOVA như những gì bạn vừa thấy gần như không thể thực hiện thủ công, với lượng dữ liệu lớn và các công thức phức tạp. statsmodels làm cho việc mô hình hóa tinh vi trở nên sẵn có trong tầm tay bạn.

7. Kết luận

Vậy là bạn đã chính thức làm chủ hồi quy tuyến tính kiểu R ngay trong Python – từ y ~ x1 + x2 quen thuộc đến kết quả summary() chi tiết, tất cả chỉ với vài dòng code ngắn gọn, mạnh mẽ và dễ đọc nhờ statsmodels và cú pháp patsy. Không còn phải chuyển qua lại giữa hai ngôn ngữ, bạn giờ đây có thể phân tích dữ liệu chuyên nghiệp, trực quan hóa mô hình và đưa ra dự báo chính xác chỉ trong một môi trường duy nhất: Python.

Đây chính là bước tiến để bạn biến những dòng dữ liệu thô thành quyết định thông minh – dù là nghiên cứu khoa học, phân tích kinh doanh hay học máy cơ bản.

Xem thêm: Hướng dẫn cách sao lưu và khôi phục tệp, ứng dụng, cài đặt bằng OneDrive trên Windows 11

Nếu bạn đang cần máy trạm cấu hình cao cho Data Science, laptop AI, phần mềm bản quyền (Python, R, Jupyter, VS Code…) với giá tốt và hỗ trợ cài đặt sẵn môi trường lập trình, hãy đến ngay COHOTECH – cửa hàng công nghệ chuyên sâu cho dân IT, Data Analyst & Developer. Chúng tôi không chỉ bán máy, mà còn tặng kèm hướng dẫn tối ưu hiệu suất để bạn làm việc mượt mà nhất!

Bạn đã chạy thử mô hình hồi quy kiểu R trong Python chưa? Kết quả có giống với R không? Hãy bình luận bên dưới chia sẻ output của bạn, hoặc hỏi thêm nếu cần hỗ trợ code! Và đừng quên chia sẻ bài viết này đến cộng đồng Data, giúp mọi người cùng tiến bộ nhé!