Giới thiệu
Các phần tử bất thường (Outliers) có thể là nguyên nhân gây ra sự vi phạm các giả thiết của mô hình hồi quy tuyến tính và làm méo mó kết quả dự báo của mô hình.
Có hai tiêu chuẩn được sử dụng để đánh giá Outliers. Tiêu chuẩn thứ nhất (như phần mềm SPSS sử dụng) là bất cứ quan sát nào nằm ngoài 2 lần độ lệnh chuẩn so với mean của phần dư sẽ là các quan sát bất thường.
Tiêu chuẩn thứ hai là căn cứ vào khoản cách Cook - Cook’s Distance.
Kết quả chỉ ra rằng:
Nếu loại Outliers theo tiêu chuẩn thứ nhất thì mức độ giải thích của mô hình tăng gần 5 lần.
Nếu loại Outliers theo tiêu chuẩn thứ hai thì mức độ giải thích của mô hình tăng xấp xỉ 7 lần.
Số liệu sử dụng ở đây là bộ dữ liệu ceosal1 được sử dụng ở chương 2 cuốn Introductory Econometrics: A Modern Approach của Jeffrey Wooldridge.
# Load dữ liệu:
rm(list = ls())
require(readstata13)
ceosal1 <- read.dta13("http://fmwww.bc.edu/ec-p/data/wooldridge/ceosal1.dta")
# Xem qua:
library(tidyverse)
ceosal1 %>% head()## salary pcsalary sales roe pcroe ros indus finance consprod utility
## 1 1095 20 27595.0 14.1 106.4 191 1 0 0 0
## 2 1001 32 9958.0 10.9 -30.6 13 1 0 0 0
## 3 1122 9 6125.9 23.5 -16.3 14 1 0 0 0
## 4 578 -9 16246.0 5.9 -25.7 -21 1 0 0 0
## 5 1368 7 21783.2 13.8 -3.0 56 1 0 0 0
## 6 1145 5 6021.4 20.0 1.0 55 1 0 0 0
## lsalary lsales
## 1 6.998509 10.225389
## 2 6.908755 9.206132
## 3 7.022868 8.720281
## 4 6.359574 9.695602
## 5 7.221105 9.988894
## 6 7.043160 8.703075#----------------------------------------------
# Version 1: OLS Model với đầy đủ quán sát
#----------------------------------------------
ols1 <- ceosal1 %>% lm(salary ~ roe, data = .)
ols1 %>% summary()##
## Call:
## lm(formula = salary ~ roe, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1160.2 -526.0 -254.0 138.8 13499.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 963.19 213.24 4.517 1.05e-05 ***
## roe 18.50 11.12 1.663 0.0978 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1367 on 207 degrees of freedom
## Multiple R-squared: 0.01319, Adjusted R-squared: 0.008421
## F-statistic: 2.767 on 1 and 207 DF, p-value: 0.09777theme_set(theme_minimal())
ceosal1 %>%
ggplot(aes(roe, salary)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm")
# Bằng mắt thường chúng ta có thể thấy có một số quan sát
# có mức lương rất cao. Chúng là những outliers. Nhưng ta
# đã biết outliers có ảnh hưởng rất mạnh đến kết quả của
# OLS. Dưới đây chúng ta sẽ tìm cách loại bỏ các outliers
# theo một số tiêu chuẩn thống kê và thực hiện lại OLS.
# Trước hết chúng ta chẩn đoán các lỗi có thể có của mô
# hình cũng như một số quan sát bất thường:
library("ggfortify")
autoplot(ols1, which = 1:6, ncol = 3, label.size = 1)
#------------------------------------------------------
# Version 2: Loại bỏ outliers theo tiêu chuẩn
# nằm ngoài 2 lần độ lệch chuẩn so với mean phần dư
#-------------------------------------------------------
ct <- mean(ols1$residuals) + 2*sd(ols1$residuals)
cd <- mean(ols1$residuals) - 2*sd(ols1$residuals)
library(magrittr)
ceosal1 %<>%
mutate(phan_du = ols1$residuals,
out = case_when(phan_du > ct ~ "Outlier",
phan_du < cd ~ "Outlier",
phan_du <= ct & phan_du >= cd ~ "Normal"))
# Các quan sát là outliers theo tiêu chuẩn này:
ceosal1 %>% filter(out == "Outlier")## salary pcsalary sales roe pcroe ros indus finance consprod utility
## 1 4143 -10 2678.4 14.4 -12.8 42 0 1 0 0
## 2 6640 4 8946.0 10.2 87.6 34 0 1 0 0
## 3 11233 17 6047.9 22.9 9.9 74 0 0 1 0
## 4 14822 1 2159.2 19.4 -37.3 17 0 0 1 0
## lsalary lsales phan_du out
## 1 8.329175 7.892975 2913.392 Outlier
## 2 8.800867 9.098962 5488.097 Outlier
## 3 9.326612 8.707466 9846.132 Outlier
## 4 9.603868 7.677493 13499.886 Outlier# Nếu chúng ta loại bỏ chúng ra và thực hiện lại OLS thì:
ols2 <- ceosal1 %>%
filter(out != "Outlier") %>%
lm(salary ~ roe, data = .)
ols2 %>% summary()##
## Call:
## lm(formula = salary ~ roe, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -990.8 -372.2 -104.5 266.9 2798.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 853.563 87.866 9.714 < 2e-16 ***
## roe 15.870 4.576 3.468 0.000639 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 560.4 on 203 degrees of freedom
## Multiple R-squared: 0.05595, Adjusted R-squared: 0.0513
## F-statistic: 12.03 on 1 and 203 DF, p-value: 0.000639# Đường hồi quy của mô hình này cũng như minh họa đủ mọi thứ trong
# đó biểu diễn các outliers, đường màu blue là đường hồi quy đủ,
# đường green là đường hồi quy đã bỏ đi 4 (hay 1.9% quan sát) là outliers:
library(ggrepel)
ceosal1 %>%
ggplot(aes(roe, salary)) +
geom_point(data = ceosal1 %>% filter(out == "Outlier"),
color = "red", size = 2) +
geom_smooth(method = "lm", se = FALSE) +
geom_point(data = ceosal1 %>% filter(out != "Outlier"),
color = "purple", alpha = 0.3) +
geom_smooth(data = ceosal1 %>% filter(out != "Outlier"),
method = "lm", se = FALSE, color = "green") +
geom_text_repel(aes(label = out),
data = ceosal1 %>% filter(out == "Outlier"), force = 19)
# Có thể thấy là sau khi loại bỏ 4 quan sát này thì kết quả OLS
# thay đổi đáng kể. Thứ nhất là hệ số góc:
ols1$coefficients # Mô hình ban đầu. ## (Intercept) roe
## 963.19134 18.50119ols2$coefficients # Mô hình loại bỏ 4 quan sát bất thường. ## (Intercept) roe
## 853.56290 15.87006# Đặc biệt là mức độ giải thích của mô hình thay đổi rất lớn.
# Cụ thể mô hình sau khi bỏ outliers có khả năng giải thích tăng 5 lần:
ols1 %>% summary() -> u
u$r.squared # Mô hình ban đầu. ## [1] 0.01318862ols2 %>% summary() -> v
v$r.squared # Mô hình bỏ bốn outliers. ## [1] 0.05594679# Nghĩa là mô hình 2 nếu được sử dụng để dự báo mức lương
# phụ thuộc ra sao vào roe - hiệu quả điều hành công ti của
# các CEO thì sẽ dự báo chính xác hơn so với mô hình 1.
#---------------------------------------------------------------
# Version 3: Loại bỏ outlier theo
# khoảng cách Cook (Cook’s Distance)
# Đọc thêm: https://en.wikipedia.org/wiki/Cook%27s_distance
#---------------------------------------------------------------
# Tính khoản cách Cook:
cut_off <- 1 / nrow(ceosal1)
# Xếp loại quan sát theo tiêu chuẩn khoảng cách Cook:
ceosal1 %<>% mutate(d = cooks.distance(ols1),
out_d = case_when(d > cut_off ~ "Outlier",
d <= cut_off ~ "Normal"))
# Các quan sát là bất thường theo tiêu chuẩn Cook:
ceosal1 %>%
filter(out_d == "Outlier") %>%
select(-indus, -finance, -phan_du, -d)## salary pcsalary sales roe pcroe ros consprod utility lsalary
## 1 3844 61 20604.0 12.1 -5.9 18 0 0 8.254269
## 2 4143 -10 2678.4 14.4 -12.8 42 0 0 8.329175
## 3 6640 4 8946.0 10.2 87.6 34 0 0 8.800867
## 4 612 39 3618.9 33.3 -10.5 68 0 0 6.416732
## 5 3142 15 10236.3 35.7 27.2 155 1 0 8.052615
## 6 1130 1 1686.3 44.5 3.4 99 1 0 7.029973
## 7 11233 17 6047.9 22.9 9.9 74 1 0 9.326612
## 8 3646 -11 3921.5 7.8 -64.2 29 1 0 8.201385
## 9 14822 1 2159.2 19.4 -37.3 17 1 0 9.603868
## lsales out out_d
## 1 9.933241 Normal Outlier
## 2 7.892975 Outlier Outlier
## 3 9.098962 Outlier Outlier
## 4 8.193925 Normal Outlier
## 5 9.233695 Normal Outlier
## 6 7.430292 Normal Outlier
## 7 8.707466 Outlier Outlier
## 8 8.274229 Normal Outlier
## 9 7.677493 Outlier Outlier# Theo tiêu chuẩn này thì số các quán sát bất thường là chặt hơn: có
# 9 quan sát tức tương đương 4.3 % kích thước mẫu ban đầu. Thực hiện
# hồi quy sau khi loại outlier theo tiêu chuẩn Cook:
ols3 <- ceosal1 %>%
filter(out_d != "Outlier") %>%
lm(salary ~ roe, data = .)
# Kết quả của mô hình OLS này cho thấy mức độ giải thích của mô
# hình (nếu căn cứ vào R2) cao gấp gần 7 lần (0.09001 so với 0.01318862):
ols3 %>% summary()##
## Call:
## lm(formula = salary ~ roe, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -973.68 -340.29 -71.92 289.09 1893.04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 785.79 77.08 10.195 < 2e-16 ***
## roe 18.10 4.09 4.425 1.59e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 475.7 on 198 degrees of freedom
## Multiple R-squared: 0.09001, Adjusted R-squared: 0.08541
## F-statistic: 19.58 on 1 and 198 DF, p-value: 1.588e-05# Có thể hình ảnh hóa các điểm với khoảng cách Cook tương ứng:
autoplot(ols1, which = 4, ncol = 3, label.size = 1) # Cách 1. 
# Tại cột biến N chỉ thị số thứ tự của quan sát:
u <- ceosal1 %>%
mutate(N = 1:nrow(.))
u %>%
ggplot(aes(N, d)) +
geom_point(alpha = 0.2) +
geom_point(data = u %>% filter(d > cut_off), color = "red") +
geom_text_repel(data = u %>% filter(d > cut_off),
aes(label = N), force = 19) +
labs(x = NULL, y = "Cook's Distance",
title = "Outliers based on Cook's Distance") # Cách 2. 
#-------------------------------------------
# Mở rộng
#-------------------------------------------
# Tất cả các bước dài dòng trên có thể được làm ngắn gọn như sau
# để so sánh các mô hình. Trước hết tạo ra ba bộ số liệu như sau:
df1 <- ceosal1 %>%
select(salary, roe) %>%
mutate(Model = "All") # Bộ số liệu nguyên bản.
df2 <- ceosal1 %>%
filter(out != "Outlier") %>%
select(salary, roe) %>%
mutate(Model = "Remove4") # Bộ bỏ 4.
df3 <- ceosal1 %>%
filter(out_d != "Outlier") %>%
select(salary, roe) %>%
mutate(Model = "Remove9") # Bộ bỏ 9.
# Hợp nhất cả ba bộ dữ liệu:
all_df <- bind_rows(df1, df2, df3)
library(purrr)
# So sánh kết quả hồi quy của ba phiên bản:
all_df %>%
split(.$Model) %>%
map(function(df) lm(salary ~ roe, data = df)) %>%
map(summary)## $All
##
## Call:
## lm(formula = salary ~ roe, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1160.2 -526.0 -254.0 138.8 13499.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 963.19 213.24 4.517 1.05e-05 ***
## roe 18.50 11.12 1.663 0.0978 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1367 on 207 degrees of freedom
## Multiple R-squared: 0.01319, Adjusted R-squared: 0.008421
## F-statistic: 2.767 on 1 and 207 DF, p-value: 0.09777
##
##
## $Remove4
##
## Call:
## lm(formula = salary ~ roe, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -990.8 -372.2 -104.5 266.9 2798.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 853.563 87.866 9.714 < 2e-16 ***
## roe 15.870 4.576 3.468 0.000639 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 560.4 on 203 degrees of freedom
## Multiple R-squared: 0.05595, Adjusted R-squared: 0.0513
## F-statistic: 12.03 on 1 and 203 DF, p-value: 0.000639
##
##
## $Remove9
##
## Call:
## lm(formula = salary ~ roe, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -973.68 -340.29 -71.92 289.09 1893.04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 785.79 77.08 10.195 < 2e-16 ***
## roe 18.10 4.09 4.425 1.59e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 475.7 on 198 degrees of freedom
## Multiple R-squared: 0.09001, Adjusted R-squared: 0.08541
## F-statistic: 19.58 on 1 and 198 DF, p-value: 1.588e-05# So sánh hệ số R2:
all_df %>%
split(.$Model) %>%
map(function(df) lm(salary ~ roe, data = df)) %>%
map(summary) %>%
map_dbl("r.squared")## All Remove4 Remove9
## 0.01318862 0.05594679 0.09000992# Hay so sánh các hệ số hồi quy:
all_df %>%
split(.$Model) %>%
map(function(df) lm(salary ~ roe, data = df)) %>%
map(summary) %>%
map("coefficients")## $All
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 963.19134 213.24026 4.516930 1.053352e-05
## roe 18.50119 11.12325 1.663289 9.776775e-02
##
## $Remove4
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 853.56290 87.865950 9.714376 1.458024e-18
## roe 15.87006 4.575527 3.468466 6.390450e-04
##
## $Remove9
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 785.79549 77.078456 10.19475 6.907394e-20
## roe 18.10064 4.090106 4.42547 1.587915e-05