Các dạng hồi quy
Có rất nhiều loại hồi quy như: Simple linear, Polynomial, Multiple linear, Multilevel, Multivariate, Logistic, Poisson, Cox proportional hazards, Time-series, Nonlinear, Nonparametric
Trong đó, 4 loại phổ biến hay được sử dụng là:
Thực hành hồi quy đơn biến
Như đã giới thiệu ở phần lý thuyết trước, hồi quy đơn biến được mô tả như sau: Là mô hình hồi quy với 1 biến phụ thuộc, 1 biến độc lập. Hàm hồi quy có dạng Y = a + bX
Y là biến phụ thuộc, X là biến độc lập, a là hệ số chặn (intercept), b là hệ số góc (coefficient)
Ví dụ: Với bộ dữ liệu women, có 2 biến là chiều cao, cân nặng của phụ nữ. Muốn dự báo cân nặng của phụ nữ dựa vào chiều cao của họ ta xây dựng mô hình hồi quy đơn biến. Bộ dữ liệu này có sẵn trong thư viện của RStudio, chúng ta có thể lấy ngay được bộ dữ liệu này:
data("women")
women %>% head() %>% datatable()women %>% ggplot(aes(x = height, y = weight)) +
geom_point( color = "red") +
xlab("Height (in inches)") + # 1 inch = 2.54 cm
ylab("Weight (in pounds)")+ # 1 pound = 0.45 kg
ggtitle("Graph to show the correlation")
Áp dụng hồi quy tuyến tính, sử dụng hàm lm(), phương pháp tính toán là MSE - Mean of Sum Square.
fit <- lm(weight ~ height, data = women )
summary(fit)##
## Call:
## lm(formula = weight ~ height, data = women)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7333 -1.1333 -0.3833 0.7417 3.1167
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -87.51667 5.93694 -14.74 1.71e-09 ***
## height 3.45000 0.09114 37.85 1.09e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-squared: 0.991, Adjusted R-squared: 0.9903
## F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14Như vậy phương trình hồi quy có dạng: \[\hat{weight}=-87.51667 + 3.45*height\] Các giá trị dự báo của mô hình này, tính bằng hàm fitted(), phần dư residual của mô hình, sử dụng hàm residuals(fit).
fitted(fit) %>% head(n=4)## 1 2 3 4
## 112.5833 116.0333 119.4833 122.9333women$predicton_value <- fitted(fit)
women$residual_value <- residuals(fit)Vẽ đồ thị các giá trị thực tế và các giá trị dự đoán từ mô hình
women %>% ggplot() +
geom_point(aes(x = height, y = weight), col = "blue", size =2, alpha = 0.8) +
geom_point(aes(x = height, y = predicton_value), col = "red", size = 2, alpha = 0.8) +
xlab("Height (in inches)") + # 1 inch = 2.54 cm
ylab("Weight (in pounds)")+ # 1 pound = 0.45 kg
ggtitle("Graph to show the correlation")
Nếu chúng ta muốn vẽ đường hồi quy và thêm khoảng tin cậy
ggplot(data = women, aes(x = height, y = weight)) +
geom_point(col = "red", shape =3, alpha = 0.7)+
geom_smooth(method = "lm", se = TRUE)+
xlab("Height (in inches)") + # 1 inch = 2.54 cm
ylab("Weight (in pounds)")+ # 1 pound = 0.45 kg
ggtitle("Graph to show the correlation")## `geom_smooth()` using formula 'y ~ x'
Thực hành hồi quy đa biến:
Mô hình hồi quy dạng: \[\hat{Y}=\beta_{0} +\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{3}X_{3}+...\] Trong phần này, chúng ta sử dụng bộ dữ liệu States trình bày các thông tin về 50 bang của nước Mỹ. Cụ thể chúng ta lấy dữ liệu, có tên là states.x77, bộ dữ liệu gồm 50 dòng, tương ứng với 50 bang của nước Mỹ và 8 cột biến với các nội dung: Population (tháng 1 năm 1975), Income (thu nhập bình quân đầu người - 1974), tỷ lệ mù chữ của dân số - Illiteracy, tỷ lệ phạm tội giết người và ngộ sát trên 100 nghìn dân Murder
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
states %>% head() %>% datatable()Xem phân tích tương quan:
chart.Correlation(states, histogram=TRUE) 
Chúng ta cùng xem mối tương quan giữa tỷ lệ tội phạm giết người và ngộ sát có liên quan gì đến các biến còn lại sử dụng hồi quy đa biến:
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
summary(fit)##
## Call:
## lm(formula = Murder ~ Population + Illiteracy + Income + Frost,
## data = states)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7960 -1.6495 -0.0811 1.4815 7.6210
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.235e+00 3.866e+00 0.319 0.7510
## Population 2.237e-04 9.052e-05 2.471 0.0173 *
## Illiteracy 4.143e+00 8.744e-01 4.738 2.19e-05 ***
## Income 6.442e-05 6.837e-04 0.094 0.9253
## Frost 5.813e-04 1.005e-02 0.058 0.9541
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.535 on 45 degrees of freedom
## Multiple R-squared: 0.567, Adjusted R-squared: 0.5285
## F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08Giá trị p-value Income và Frost > 0.05. Vì vậy, với mức ý nghĩa 0.05 có thể kết luận hệ số của 2 biến này bằng 0. Biến Frost ở đây là số ngày trung bình có nhiệt độ dưới mức đóng băng, là yếu tố về thời tiết do đó trong thực tế yếu tố này không có sự liên quan gì tới tỷ lệ tội phạm giết người. Biến Income cũng tương tự, thu nhập bình quân đầu người không có nhiều ảnh hưởng đến tỷ lệ phạm tội. Ở đây có 2 biến số quan trọng là dân số và tỷ lệ mù chữ ảnh hưởng rất nhiều đến tỷ lệ phạm tội giết người.
Tính phần dư và giá trị dự đoán:
states$pre <- fitted(fit)
states$res <- residuals(fit)Khoảng tin cậy của hệ số hồi quy
confint(fit)## 2.5 % 97.5 %
## (Intercept) -6.552191e+00 9.0213182149
## Population 4.136397e-05 0.0004059867
## Illiteracy 2.381799e+00 5.9038743192
## Income -1.312611e-03 0.0014414600
## Frost -1.966781e-02 0.0208304170Thay vì việc nhìn vào pvalue ở mô hình hồi quy, khoảng tin cậy của hệ số hồi quy cho phép đánh giá xem hệ số hồi quy có khác 0 hay không.
Khoảng tin cậy của Hệ số của Intercept, Income, Frost có cận dưới < 0, trong khi cận trên > 0 –> nhìn vào đây có thể biết được 3 hệ số này = 0 với mức ý nghĩa 0.05%
Kiểm định mô hình hồi quy
Dựa vào biểu đồ phần dư
par(mfrow=c(2,2))
plot(fit)
Biểu đồ 1: Vẽ tương quan giữa phần dư và kết quả dự báo, giá trị phần dư càng ở quanh mức 0 kết quả dự báo càng tốt. Như vậy, tỷ lệ giết người đang được dự báo quá cao thấp với thực tế tại 2 bang Rhode Island và Masschusetts, dự báo quá cao tại bang Nevada
Biểu đồ 2: Kiểm tra xem phần dư có phân phối chuẩn N(0,1) hay không, kết quả cho thấy phần dư tại 3 bang Nevada và Rhode Island và Alaska đang ko theo quy luật phân phối chuẩn
Biểu đồ 3: Đánh giá phương sai của phần dư có đồng nhất hay không
Biều đồ 4: Cho phép phát hiện ra các outliers trong phần dư
Dựa vào các kiểm định
Giả thiết 1: Sai số ngẫu nhiên có phân phối chuẩn
Giả thiết 2: Kỳ vọng của sai số ngẫu nhiên tại mỗi giá trị bằng 0
res <- residuals(fit)
t.test(res, mu = 0)##
## One Sample t-test
##
## data: res
## t = 6.5292e-16, df = 49, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.6903919 0.6903919
## sample estimates:
## mean of x
## 2.243128e-16p-value = 1: mô hình thỏa mãn giả thiết 2
Giả thiết 3: Phương sai của sai số ngẫu nhiên không đổi
H0: phương sai sai số không đổi Ha: Phương sai sai số thay đổi
ncvTest(fit)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 1.746514, Df = 1, p = 0.18632p-value = 0.18, có thể kết luận phương sai sai số không đổi với mức ý nghĩa < 0.18
Giả thiết 4: Giữa các biến độc lập không có mối quan hệ đa cộng tuyến hoàn hảo
Chú ý rằng ngoài hệ số VIF người ta còn lấy TOL = 1/VIF (nghịch đảo của VIF) là tiêu chí nhận định về hiện tượng đa cộng tuyến. Hiện tại, ngưỡng VIF (hay TOL, r23) bằng bao nhiêu để chỉ ra hiện tượng đa cộng tuyến vẫn là một vấn đề chưa có sự thống nhất giữa các nhà kinh tế lượng. Gujarati & Porter (2009) chỉ ra một số dấu hiệu của hiện tượng đa cộng tuyến trong mô hình khi:
- hệ số VIF ≥ 10
- hệ số tương quan r của bất kì cặp biến nào trong mô hình lớn hơn 0.8.
Trong khi đó Allison (1999) đưa ra tiêu chí chặt hơn khi chọn VIF > 2.5 (hay r > 0.775).
vif(fit)## Population Illiteracy Income Frost
## 1.245282 2.165848 1.345822 2.082547vif < 10 cho thấy các biến độc lập không có đa cộng tuyến
Lựa chọn mô hình
Sau khi kiểm định các giả thiết của mô hình MSE, thử nghiệm trường hợp bỏ 2 biến Income, Frost ra khỏi mô hình
fit <- lm(Murder ~ Population + Illiteracy, data = states)
summary(fit)##
## Call:
## lm(formula = Murder ~ Population + Illiteracy, data = states)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7652 -1.6561 -0.0898 1.4570 7.6758
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.652e+00 8.101e-01 2.039 0.04713 *
## Population 2.242e-04 7.984e-05 2.808 0.00724 **
## Illiteracy 4.081e+00 5.848e-01 6.978 8.83e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.481 on 47 degrees of freedom
## Multiple R-squared: 0.5668, Adjusted R-squared: 0.5484
## F-statistic: 30.75 on 2 and 47 DF, p-value: 2.893e-09par(mfrow=c(2,2))
plot(fit)
par(mfrow=c(1,1))