Ví dụ về Hồi quy tuyến tính bội
H.V.Ha
11/12/2020
Mô tả dữ liệu
Tập tin marketing.csv dùng để phân tích ảnh hưởng của các hình thức quảng cáo lên doanh thu. Dữ liệu bao gồm 4 biến: \(youtube\), \(facebook\) và \(newspaper\) là số tiền chi cho quảng cáo (Đơn vị: \(1000\) usd), và \(sales\) là biến doanh thu. Thí nghiệm được lặp lại \(200\) lần.
Yêu cầu: xây dựng mô hình hồi quy bội (multiple regression) để phân tích tác động của các hình thức quảng cáo lên doanh thu và diễn giải kết quả.
Đọc dữ liệu, gán vào data.frame marketing:
marketing <- read.csv("marketing.csv", header = T, sep = ";")
marketing[1:20,] # Hien thi 20 hang dau tien cua bo du lieu ## youtube facebook newspaper sales
## 1 276.12 45.36 83.04 26.52
## 2 53.40 47.16 54.12 12.48
## 3 20.64 55.08 83.16 11.16
## 4 181.80 49.56 70.20 22.20
## 5 216.96 12.96 70.08 15.48
## 6 10.44 58.68 90.00 8.64
## 7 69.00 39.36 28.20 14.16
## 8 144.24 23.52 13.92 15.84
## 9 10.32 2.52 1.20 5.76
## 10 239.76 3.12 25.44 12.72
## 11 79.32 6.96 29.04 10.32
## 12 257.64 28.80 4.80 20.88
## 13 28.56 42.12 79.08 11.04
## 14 117.00 9.12 8.64 11.64
## 15 244.92 39.48 55.20 22.80
## 16 234.48 57.24 63.48 26.88
## 17 81.36 43.92 136.80 15.00
## 18 337.68 47.52 66.96 29.28
## 19 83.04 24.60 21.96 13.56
## 20 176.76 28.68 22.92 17.521) Xây dựng mô hình
Mô hình hồi quy bội bao gồm:
- Biến phụ thuộc: doanh thu \(sales\)
- Biến dư báo (biến độc lập): \(youtube\), \(facebook\), \(newspaper\) lần lượt là số tiền chi cho quảng cáo trên Youtube, Facebook và trên báo chí
Mô hình được biểu diễn như sau: \[ sales = \beta_0 + \beta_1\times youtube + \beta_2\times facebook + \beta_3\times newspaper + \epsilon. \]
Ta có thể vẽ đồ thị phân tán của từng cặp biến để nhận xét mối quan hệ tuyến tính giữa chúng:
par(mfrow = c(1,3))
plot(marketing$youtube, marketing$sales, pch = 16, col = 'blue', main = "Youtube vs Sales")
plot(marketing$facebook, marketing$sales, pch = 16, col = 'blue', main = "Facebook vs Sales")
plot(marketing$newspaper, marketing$sales, pch = 16, col = 'blue', main = "Newspaper vs Sales")
Từ các đồ thị phân tán, ta có thể nhận xét rằng các biến \(youtube\) và \(facebook\) có mối quan hệ tuyến tính với biến \(sales\) trong khi biến \(newspaper\) không có. Mặt khác, mối quan hệ giữa \(youtube\) và \(sales\), nói một cách chính xác hơn là phi tuyến tính hơn là tuyến tính.
Ta ước lượng các hệ số \(\beta_i\), \(i=0,\ldots,3\) sử dụng lệnh lm():
model.marketing <- lm(sales ~ youtube + facebook + newspaper, data = marketing)
summary(model.marketing)##
## Call:
## lm(formula = sales ~ youtube + facebook + newspaper, data = marketing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.5932 -1.0690 0.2902 1.4272 3.3951
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.526667 0.374290 9.422 <2e-16 ***
## youtube 0.045765 0.001395 32.809 <2e-16 ***
## facebook 0.188530 0.008611 21.893 <2e-16 ***
## newspaper -0.001037 0.005871 -0.177 0.86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.023 on 196 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8956
## F-statistic: 570.3 on 3 and 196 DF, p-value: < 2.2e-16Từ kết quả phân tích, ta thu được \(\hat \beta_0 = 3.526667\), \(\hat\beta_1 = 0.045765\), \(\hat\beta_2 = 0.188530\) và \(\hat\beta_3 = -0.001037\). Như vậy, đường thẳng hồi quy ước lượng cho bởi phương trình sau:
\[ \widehat{sales} = 3.526667 + 0.045765\times youtube + 0.188530\times facebook -0.001037\times newspaper. \]
2) Phân tích kết quả
Trước hết, ta thấy rằng \(p\)-value tương ứng với thống kê \(F\) bé hơn \(2.2\times 10^{-16}\), có ý nghĩa rất cao. Điều này chỉ ra rằng, ít nhất một biến dự báo trong mô hình có ý nghĩa giải thích rất cao cho biến doanh thu \(sales\).
Để xét ảnh hưởng cụ thể của từng biến độc lập, ta xét trọng số (hệ số \(\beta_i\)) và \(p\)-value tương ứng. Ta thấy rằng \(p\)-value tương ứng với hai biến \(youtube\) và \(facebook\) đều bé hơn \(2\times 10^{-16}\), điều này nói lên rằng ảnh hưởng của hai biến có ý nghĩa rất cao lên biến doanh thu \(sales\). Mặt khác, hệ số hồi quy \(\beta_i\) của một biến dự báo cũng có thể được xem như ảnh hưởng trung bình lên biến phụ thuộc doanh thu khi tăng một đơn vị của biến dự báo đó, giả sử rằng các biến dự báo khác không đổi. Cụ thể, \(\hat\beta_1 = 0.045765\) thì với mỗi \(1000\) usd gia tăng trên chi phí quảng cáo trên Youtube, ta có thể kỳ vọng doanh thu sẽ tăng lên \(0.045765\times 1000 \approx 46\) đơn vị về mặt trung bình (giả sử rằng chi phí quảng cáo trên Facebook và báo chí không đổi). Tương tự, với mỗi \(1000\) usd gia tăng trên chi phí quảng cáo trên Facebook, ta có thể kỳ vọng rằng doanh thu sẽ tăng lên \(0.188530\times 1000 \approx 189\) đơn vị về mặt trung bình.
Ta cũng nhận thấy rằng \(newspaper\) không có nhiều ảnh hưởng lên doanh thu (\(\hat\beta_3 = -0.001037\)). \(p\)-value của bài toán kiểm định giả thuyết \(H_0:\beta_3 = 0\) bằng \(0.86\) cho thấy rằng biến này không có ý nghĩa đối với mô hình hồi quy bội.
Vì biến \(newspaper\) không có ý nghĩa trong mô hình hồi quy bội được xây dựng, ta có thể loại bỏ biến này và xây dựng lại mô hình hồi quy bội chỉ với hai biến \(youtube\) và \(facebook\):
model2.marketing <- lm(sales ~ youtube + facebook, data = marketing)
summary(model2.marketing) ##
## Call:
## lm(formula = sales ~ youtube + facebook, data = marketing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.5572 -1.0502 0.2906 1.4049 3.3994
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.50532 0.35339 9.919 <2e-16 ***
## youtube 0.04575 0.00139 32.909 <2e-16 ***
## facebook 0.18799 0.00804 23.382 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.018 on 197 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8962
## F-statistic: 859.6 on 2 and 197 DF, p-value: < 2.2e-16Như vậy, mô hình hồi quy bội về ảnh hưởng của chi phí quảng cáo theo các hình thứ khác nhau (Youtube, Facebook) lên doanh thu được cho bởi:
\[ \widehat{sales} = 3.50532 + 0.04575\times youtube + 0.18799\times facebook. \] Hệ số \(R^2\) hiệu chỉnh bằng \(0.8962\) nghĩa là \(89.62\%\) sự biến thiên trong doanh thu \(sales\) được giải thích bởi các biến \(youtube\) và \(facebook\) tức là chi phí chi cho việc quảng cáo trên Facebook và Youtube.
Đề tìm khoảng tin cậy cho các hệ số hồi quy, ta sử dụng hàm confint():
confint(model2.marketing)## 2.5 % 97.5 %
## (Intercept) 2.80841159 4.20222820
## youtube 0.04301292 0.04849671
## facebook 0.17213877 0.20384969Khoảng tin cậy \(95\%\) cho các hệ số hồi quy cho bởi: \[ 2.808 \le \beta_0 \le 4.202, \\ 0.043 \le \beta_1 \le 0.048, \\ 0.172 \le \beta_2 \le 0.204. \]
3) Kiểm tra các giả định của mô hình
Nhắc lại các giả định của mô hình hồi quy: \(Y_i = \beta_0 + \beta_1X_1 + \cdots + \beta_pX_p + \epsilon_i,\quad i = 1,\ldots, n\)
- Tính tuyến tính của dữ liệu: mối quan hệ giữa biến dự báo \(X\) và biến phụ thuộc \(Y\) được giả sử là tuyến tính
- Sai số có phân phối chuẩn
- Phương sai của các sai số là hằng số: \(\epsilon_i \sim \mathcal{N}(0, \sigma^2)\)
- Các sai số \(\epsilon_1, \ldots, \epsilon_n\) thì độc lập với nhau
Ta thực hiện phân tích thặng dư để kiểm tra các giả định của mô hình:
par(mfrow = c(2, 2))
plot(model2.marketing)
Đồ thị thứ 1 (Residuals vs Fitted) vẽ các giá trị dự báo (\(\hat y_i\)) với các giá trị thặng dư (sai số) tương ứng, dùng để kiểm tra tính tuyến tính của dữ liệu (giả định 1) và tính đồng nhất của các phương sai sai số (giả định 3). Nếu như giả định về tính tuyến tính của dữ liệu KHÔNG thỏa, ta sẽ quan sát thấy rằng các điểm thặng dư (residuals) trên đồ thị sẽ phân bố theo một hình mẫu (pattern) đặc trưng nào đó (ví dụ parabol). Nếu đường màu đỏ trên đồ thị phân tán là đường thẳng nằm ngang mà không phải là đường cong, thì giả định tính tuyến tính của dữ liệu được thỏa mãn. Để kiểm tra giả định thứ 3 (phương sai đồng nhất) thì các điểm thặng dự phải phân tán đều nhau xung quanh đường thẳng \(y = 0\).
Đồ thị thứ 2 (Normal Q-Q) cho phép kiểm tra giả định về phân phối chuẩn của các sai số. Nếu các điểm thặng dư nằm trên cùng 1 đường thẳng thì điều kiện về phân phối chuẩn được thỏa.
Đồ thị thứ 3 (Scale - Location) vẽ căn bậc hai của các giá trị thặng dư được chuẩn hóa với các giá trị dự báo, được dùng để kiểm tra giả định thứ 3 (phương sai của các sai số là hằng số). Nếu như đường màu đỏ trên đồ thị là đường thẳng nằm ngang và các điểm thặng dư phân tán đều xung quanh đường thẳng này thì giả định thứ 3 được thỏa. Nếu như đường màu đỏ có độ dốc (hoặc cong) hoặc các điểm thặng dư phân tán không đều xung quanh đường thẳng này, thì giả định thứ 3 bị vi phạm.
Đồ thị thứ 4 (Residuals vs Leverage) cho phép xác định những điểm có ảnh hưởng cao (influential observations), nếu chúng có hiện diện trong bộ dữ liệu. Những điểm có ảnh hưởng cao này có thể là các điểm outliers, là những điểm có thể gây nhiều ảnh hưởng nhất khi phân tích dữ liệu. Nếu như ta quan sát thấy một đường thẳng màu đỏ đứt nét (Cook’s distance), và có một số điểm vượt qua đường thẳng khoảng cách này, nghĩa là các điểm đó là các điểm có ảnh hưởng cao. Nếu như ta chỉ quan sát thấy đường thẳng khoảng cách Cook ở góc của đồ thị và không có điểm nào vượt qua nó, nghĩa không có điểm nào thực sự có ảnh hưởng cao.
Nhận xét:
- Đồ thị Normal Q-Q cho thấy giả định sai số có phân phối chuẩn được thỏa mãn.
- Đồ thị thứ 1 (Residuals vs Fitted) cho thấy giả định về tính tuyến tính của dữ liệu hơi bị vi phạm, ta có thể thấy rằng sự vi phạm này bởi vì mối quan hệ giữa \(sales\) và \(youtube\) là phi tuyến tính.
- Đồ thị thứ 1 và thứ 3 (Scale - Location) cho ta thấy rằng giả định về tính đồng nhất của phương sai cũng hơi bị vi phạm. Tuy nhiên, ta cũng thấy này sự vi phạm này tương đối nhỏ và có thể chấp nhận được.
- Đồ thị thứ tư chỉ ra có các quan trắc thứ 6, 36 và 131 có thể là các điểm có ảnh hưởng cao trong bộ dữ liệu.
plot(model2.marketing, 5)
Tuy nhiên ta cũng quan sát thấy rằng các điểm này chưa vượt qua đường thẳng khoảng cách Cook (đường thẳng đứt nét màu đỏ Cook’s distance). Do vậy, các điểm này chưa thực sự là các điểm có ảnh hưởng cao trong bộ dữ liệu. Do đó ta không cần phải loại bỏ chúng khi phân tích.
Để vẽ từng đồ thị, ta dùng các lệnh sau:
plot(model2.marketing, 1)
plot(model2.marketing, 2)
plot(model2.marketing, 3)
plot(model2.marketing, 5)