Hướng dẫn vẽ đường hồi quy tuyến tính trong R

Dataset mtcars

Sử dụng dataset mtcars để đánh giá ảnh hưởng của khối lượng các dòng xe khác nhau, cột wt (đơn vị 1000 pound Mỹ, 1 pound Mỹ tương đương 0.45 kg) đến mức độ tiêu thụ xăng, cột mpg (đơn vị số dặm di chuyển trên mỗi gallon xăng, 1 gallon khoảng 3.78 lít).

# A data frame with 32 observations on 11 (numeric) variables.
# 
# [, 1] mpg     Miles/(US) gallon
# [, 2] cyl     Number of cylinders
# [, 3] disp    Displacement (cu.in.)
# [, 4] hp      Gross horsepower
# [, 5] drat    Rear axle ratio
# [, 6] wt      Weight (1000 lbs)
# [, 7] qsec    1/4 mile time
# [, 8] vs      Engine (0 = V-shaped, 1 = straight)
# [, 9] am      Transmission (0 = automatic, 1 = manual)
# [,10] gear    Number of forward gears
# [,11] carb    Number of carburetors

## In dataset
mtcars

##                      mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
## Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
## Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
## Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
## Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
## Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
## Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
## Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
## Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
## Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
## Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
## AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
## Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
## Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
## Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
## Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
## Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
## Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
## Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
## Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2

Công thức hồi quy

Linear regression | Hồi quy tuyến tính

Bước 1: Tính các thông số cho đường hồi quy tuyến tính

fit_lm <- lm(mpg ~ wt, data = mtcars) # Hàm lm() == Fitting Linear Models
coef_lm <- coef(lm(mpg ~ wt, data = mtcars)) # Tách hệ số b0 và b1 trong phương trình hồi quy
round(coef_lm, 4) -> ok # Làm tròn đến 4 chữ số sau dấu thập phân
round(unclass(summary(fit_lm))$r.squared, 4) -> r_ok # Tách hệ số multiple R-squared

## Kết quả sau khi fit theo linear regression
summary(fit_lm) 

## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Bước 2: Vẽ đồ thị

## Vẽ đồ thị scatter plot
plot(mpg ~ wt, data = mtcars, 
     pch = 19, col = 4, 
     xlim = c(1, 6), ylim = c(5, 35), 
     xlab = paste("wt =", "Khối lượng xe", "(\u00D7", "1000 pound Mỹ)"),
     ylab = paste("mpg =", "Mức độ tiêu thụ xăng", "(dặm/gallon)"))
title(main = "Tương quan giữa khối lượng xe và mức độ tiêu thụ xăng")     
mysubtitle <- "Linear regression | Trend line | Hồi quy tuyến tính"
mtext(line = 0.25, mysubtitle)

## Chèn đường trend line
abline(fit_lm, col = 3, lwd = 2)

## Chèn công thức toán
text(x = 4, y = max(mtcars$mpg)-3,  
     bquote(atop(Y == ~ .(ok[1]) + (.(ok[2])) ~ "\u00D7" ~ X, R^2 == .(r_ok))),
     cex = 1.2)

Exponential regression | Hồi quy phi tuyến theo log nepe

Bước 1: Tính các thông số cho đường hồi quy phi tuyến theo log nepe

fit_exponential <- lm(log(mpg) ~ wt, data = mtcars)
round(exp(fit_exponential$coefficients[1]), 3) -> he_so_a
round(fit_exponential$coefficients[2], 3) -> he_so_b
round(unclass(summary(fit_exponential))$r.squared, 4) -> r_ok_2

## Kết quả sau khi fit theo linear regression
summary(fit_exponential)

## 
## Call:
## lm(formula = log(mpg) ~ wt, data = mtcars)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.210346 -0.085932 -0.006136  0.061335  0.308623 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.83191    0.08396   45.64  < 2e-16 ***
## wt          -0.27178    0.02500  -10.87 6.31e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1362 on 30 degrees of freedom
## Multiple R-squared:  0.7976, Adjusted R-squared:  0.7908 
## F-statistic: 118.2 on 1 and 30 DF,  p-value: 6.31e-12

Bước 2: Vẽ đồ thị

## Vẽ đồ thị scatter plot
plot(mpg ~ wt, data = mtcars,
     pch = 19, col = 4,
     xlim = c(1, 6), ylim = c(5, 35),
     xlab = paste("wt =", "Khối lượng xe", "(\u00D7", "1000 pound Mỹ)"),
     ylab = paste("mpg =", "Mức độ tiêu thụ xăng", "(dặm/gallon)"))
title(main = "Tương quan giữa khối lượng xe và mức độ tiêu thụ xăng")
mysubtitle <- "Exponential regression | Hồi quy phi tuyến theo log nepe"
mtext(line = 0.25, mysubtitle)

## Chèn đường exponential regression curve
curve(he_so_a*exp(he_so_b * x), 0, 6, add = TRUE, col = 3, lwd = 2)

## Chèn công thức toán
text(x = 4,
     y = max(mtcars$mpg)-3,
     bquote(atop(Y == ~ .(he_so_a) ~ "\u00D7" ~ e ^(.(he_so_b) ~ "\u00D7" ~ X), R^2 == .(r_ok_2))),
     cex = 1.2)

Polynomial regression | Hồi quy đa thức

Bước 1: Tính các thông số cho đường hồi quy đa thức (bậc 3)

fit_polynomial <- lm(mpg ~ poly(wt, 3, raw = TRUE), data = mtcars)
fit_polynomial$coefficients[1] -> he_so_poly_b0
fit_polynomial$coefficients[2] -> he_so_poly_b1
fit_polynomial$coefficients[3] -> he_so_poly_b2
fit_polynomial$coefficients[4] -> he_so_poly_b3

## Làm tròn số
round(he_so_poly_b0, 4) -> he_so_poly_b0
round(he_so_poly_b1, 4) -> he_so_poly_b1
round(he_so_poly_b2, 4) -> he_so_poly_b2
round(he_so_poly_b3, 4) -> he_so_poly_b3
round(unclass(summary(fit_polynomial))$r.squared, 4) -> r_ok_3

## Kết quả sau khi fit theo linear regression
summary(fit_polynomial)

## 
## Call:
## lm(formula = mpg ~ poly(wt, 3, raw = TRUE), data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.506 -1.999 -0.768  1.490  6.188 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)               48.40370   15.58379   3.106  0.00431 **
## poly(wt, 3, raw = TRUE)1 -11.82598   15.46346  -0.765  0.45081   
## poly(wt, 3, raw = TRUE)2   0.68938    4.74034   0.145  0.88541   
## poly(wt, 3, raw = TRUE)3   0.04594    0.45070   0.102  0.91954   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.697 on 28 degrees of freedom
## Multiple R-squared:  0.8191, Adjusted R-squared:  0.7997 
## F-statistic: 42.27 on 3 and 28 DF,  p-value: 1.585e-10

Bước 2: Vẽ đồ thị

## Vẽ đồ thị scatter plot
plot(mpg ~ wt, data = mtcars,
     pch = 19, col = 4,
     xlim = c(1, 6), ylim = c(5, 35),
     xlab = paste("wt =", "Khối lượng xe", "(\u00D7", "1000 pound Mỹ)"),
     ylab = paste("mpg =", "Mức độ tiêu thụ xăng", "(dặm/gallon)"))
title(main = "Tương quan giữa khối lượng xe và mức độ tiêu thụ xăng")
mysubtitle <- "Polynomial regression | Hồi quy đa thức (bậc 3)"
mtext(line = 0.25, mysubtitle)

## Chèn đường polynomial regression curve
curve(he_so_poly_b0 + he_so_poly_b1*x^1 + he_so_poly_b2*x^2 + he_so_poly_b3*x^3, 0, 6, 
      add = TRUE, col = 3, lwd = 2)

## Chèn công thức toán
text(x = 4, 
     y = max(mtcars$mpg)-5,  
     bquote(atop(Y == ~ .(he_so_poly_b0) + (.(he_so_poly_b1)) ~ "\u00D7" ~ X + 
                     .(he_so_poly_b2) ~ "\u00D7" ~ X^2 + .(he_so_poly_b3) ~ "\u00D7" ~ X^3, 
                 R^2 == .(r_ok_3))),
     cex = 1.1
     )

LOWESS locally-weighted polynomial regression | Scatter plot smoothing

Bước 1: Tính thông số cho đường LOWESS smoother

lowess(mtcars$mpg ~ mtcars$wt)

## $x
##  [1] 1.513 1.615 1.835 1.935 2.140 2.200 2.320 2.465 2.620 2.770 2.780 2.875
## [13] 3.150 3.170 3.190 3.215 3.435 3.440 3.440 3.440 3.460 3.520 3.570 3.570
## [25] 3.730 3.780 3.840 3.845 4.070 5.250 5.345 5.424
## 
## $y
##  [1] 31.50709 30.62375 28.72820 27.88054 26.18179 25.69430 24.72257 23.47379
##  [9] 22.32223 21.29130 21.24459 20.71497 19.58833 19.47259 19.34050 19.17743
## [17] 17.88756 17.85655 17.85655 17.85655 17.73252 17.38098 17.10126 17.10126
## [25] 16.24386 15.98592 15.68318 15.65831 14.82813 11.98084 11.70630 11.47829

Bước 2: Vẽ đồ thị

## Vẽ đồ thị scatter plot
plot(mpg ~ wt, data = mtcars,
     pch = 19, col = 4,
     xlim = c(1, 6), ylim = c(5, 35),
     xlab = paste("wt =", "Khối lượng xe", "(\u00D7", "1000 pound Mỹ)"),
     ylab = paste("mpg =", "Mức độ tiêu thụ xăng", "(dặm/gallon)"))
title(main = "Tương quan giữa khối lượng xe và mức độ tiêu thụ xăng")
mysubtitle <- "LOWESS smoother | Hồi quy đa thức có trọng số cục bộ"
mtext(line = 0.25, mysubtitle)

## Chèn đường scatter plot smoothing
lines(lowess(mtcars$mpg ~ mtcars$wt), col = 3, lwd = 2)

LOESS locally polynomial regression | Hồi quy đa thức cục bộ

Cách 1

Bước 1: Tính thông số cho đường LOESS smoother

Thay đổi mức độ smooth từ 50% đến 90% ở tham số span trong function loess() sẽ có các đường fitting gần sát với bộ dữ liệu.

## Sắp xếp dataset theo cột `wt`
mtcars[order(mtcars$wt), ] -> mtcars_order_by_wt

## Đường fitting 1
loess50 <- loess(mpg ~ wt, data = mtcars_order_by_wt, span= 0.5)
smooth50 <- predict(loess50)  

## Đường fitting 2
loess75 <- loess(mpg ~ wt, data = mtcars_order_by_wt, span= 0.75)
smooth75 <- predict(loess75) 

## Đường fitting 3
loess90 <- loess(mpg ~ wt, data = mtcars_order_by_wt, span= 0.9)
smooth90 <- predict(loess90) 

Bước 2: Vẽ đồ thị

## Vẽ đồ thị scatter plot
plot(mpg ~ wt, data = mtcars_order_by_wt,
     pch = 19, col = 4,
     xlim = c(1, 6), ylim = c(5, 35),
     xlab = paste("wt =", "Khối lượng xe", "(\u00D7", "1000 pound Mỹ)"),
     ylab = paste("mpg =", "Mức độ tiêu thụ xăng", "(dặm/gallon)"))
title(main = "Tương quan giữa khối lượng xe và mức độ tiêu thụ xăng")
mysubtitle <- "LOESS smoother | Hồi quy đa thức cục bộ"
mtext(line = 0.25, mysubtitle)

## Chèn đường scatter plot smoothing
lines(smooth50, x = mtcars_order_by_wt$wt, col='red', lwd = 2, lty = 2)
lines(smooth75, x = mtcars_order_by_wt$wt, col='darkgreen', lwd = 2, lty = 2)
lines(smooth90, x = mtcars_order_by_wt$wt , col='purple', lwd = 2, lty = 2)

## Chú thích
legend('bottomleft', 
       legend= c('50%', '75%', '90%'),
       col=c('red', 'darkgreen', 'purple'), lwd = 2, lty = 2, 
       title = 'Mức độ smooth')

Cách 2

Bước 1: Tính thông số cho đường LOESS smoother

## Sắp xếp dataset theo cột `wt`
mtcars[order(mtcars$wt), ] -> mtcars_order_by_wt

## Đường fitting 
plx <- predict(loess(mtcars_order_by_wt$mpg ~ mtcars_order_by_wt$wt), se = TRUE)

Bước 2: Vẽ đồ thị

## Vẽ đồ thị scatter plot
plot(mpg ~ wt, data = mtcars_order_by_wt,
     pch = 19, col = 4,
     xlim = c(1, 6), ylim = c(5, 35),
     xlab = paste("wt =", "Khối lượng xe", "(\u00D7", "1000 pound Mỹ)"),
     ylab = paste("mpg =", "Mức độ tiêu thụ xăng", "(dặm/gallon)"))
title(main = "Tương quan giữa khối lượng xe và mức độ tiêu thụ xăng")
mysubtitle <- "LOESS smoother với khoảng tin cậy CI = 95%"
mtext(line = 0.25, mysubtitle)

## Chèn đường scatter plot smoothing
lines(mtcars_order_by_wt$wt, plx$fit, col = 3, lwd = 2, lty = 1)
# CI là 95%, do đó 2 tail sẽ là 2.5%
lines(mtcars_order_by_wt$wt, plx$fit - qt(0.975, plx$df)*plx$se, col = 3, lwd = 2, lty = 2)
lines(mtcars_order_by_wt$wt, plx$fit + qt(0.975, plx$df)*plx$se, col = 3, lwd = 2, lty = 2)

Tài liệu tham khảo

• https://www.r-bloggers.com/2011/05/sab-r-metrics-basics-of-loess-regression/

• https://educationalresearchtechniques.com/2018/01/05/local-regression-in-r/

• http://www.sthda.com/english/wiki/scatter-plots-r-base-graphs/

• https://www.biolinfo.com/curved-lines-in-r/

• https://stackoverflow.com/questions/22717930/how-to-get-the-confidence-intervals-for-lowess-fit-using-r

Sơ kết

Trên đây là hướng dẫn cách vẽ đồ thị thể hiện đường hồi quy trong R sử dụng base graphic. Để học R bài bản từ A đến Z, kính mời Bạn tham gia khóa học “HDSD R để xử lý dữ liệu” để có nền tảng vững chắc về R nhằm tự tay làm các câu chuyện dữ liệu của riêng mình!

Nội dung khóa học: www.tuhocr.com

Hành trình ngàn dặm bắt đầu từ bước chân đầu tiên.

ĐĂNG KÝ NGAY: https://www.tuhocr.com/register

Session info

sessionInfo()

## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.31   R6_2.5.1        lifecycle_1.0.3 jsonlite_1.8.4 
##  [5] magrittr_2.0.3  evaluate_0.19   highr_0.10      stringi_1.7.8  
##  [9] cachem_1.0.6    rlang_1.0.6     cli_3.5.0       rstudioapi_0.14
## [13] jquerylib_0.1.4 bslib_0.4.2     vctrs_0.5.1     rmarkdown_2.19 
## [17] tools_4.2.2     stringr_1.5.0   glue_1.6.2      xfun_0.36      
## [21] yaml_2.3.6      fastmap_1.1.0   compiler_4.2.2  htmltools_0.5.4
## [25] knitr_1.41      sass_0.4.4