Đánh giá mô hình hồi quy

1 Giới thiệu

Cách đây khoảng 1 năm, Nhi đã viết một bài tổng quan về các tiêu chí để đánh giá hiệu năng của một mô hình phân loại: https://rpubs.com/lengockhanhi/347941 ; Xét rằng phân tích hồi quy cũng là một vấn đề phổ biến trong nghiên cứu y học, hôm nay Nhi cũng sẽ thực hiện một bài viết tương tự cho mô hình Hồi quy.

Như các bạn đã biết, trong phân tích hồi quy chúng ta đi tìm một quy luật (mô hình) cho phép tiên lượng giá trị của một biến kết quả từ dữ liệu đầu vào. Khi làm việc này, chính là ta đã xác định một hàm để rút gọn không gian dữ liệu rộng lớn và tất cả những quan hệ phức tạp giữa các biến số trong không gian này thành một con đường đơn giản, hay nói cách khác, ta đang giản lược hóa thế giới thực thành những quy luật đơn giản. Như vậy mô hình hồi quy có thể được dùng để giải quyết mục tiêu tiên lượng lẫn diễn dịch (suy diễn thống kê về mối tương quan, hiệu ứng, so sánh,… ).Trong cả 2 mục tiêu, điều kiện quan trọng nhất đảm bảo giá trị cho kết quả tiên lượng và suy diễn thống kê, đó là mô hình phải chính xác. Nhưng làm thế nào ta có thể đánh giá được độ chính xác của mô hình ? Hay so sánh phẩm chất giữa nhiều mô hình với nhau ?

Trong bài này, Nhi sẽ giới thiệu với các bạn tất cả những chỉ số cho phép đánh giá phẩm chất của một mô hình hồi quy.

Trước hết, Nhi sẽ dùng một thí dụ minh họa đơn giản với dữ liệu DLCO mà Nhi từng dùng trước kia. Mục tiêu là tiên lượng giá trị của DLCO, một đại lượng sinh lý hô hấp, dựa vào 3 biến là Giới tính, Tuổi và chiều cao.

library(tidyverse)
library(ggpubr)
library(rsample)

df = read.csv("https://raw.githubusercontent.com/kinokoberuji/R-Tutorials/master/DLCOkeras1.csv",sep=";")%>%
  dplyr::select(Sex,Age,Height,DLCO)

head(df)

##   Sex Age Height DLCO
## 1   F  20  173.0 29.2
## 2   F  22  165.0 26.3
## 3   F  22  168.6 27.9
## 4   F  23  164.0 25.4
## 5   F  24  170.0 29.3
## 6   F  25  170.5 24.8

Nhi chia dữ liệu thành 2 phần: Trainset (n=511) và Testset (n=126)

set.seed(123)
idx=caret::createDataPartition(y=df$Age, p=0.8,list=FALSE)
trainset=df[idx,]
testset=df[-idx,]

Sau đó, Nhi dựng một mô hình Polynomial đơn giản bằng hàm glm, có nội dung: DLCO ~ Sex + poly(Age, 2) + Height.

model = glm(data=trainset,
            formula=DLCO ~ Sex + poly(Age,2)+Height)

summary(model)

## 
## Call:
## glm(formula = DLCO ~ Sex + poly(Age, 2) + Height, data = trainset)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -30.8986   -2.8045   -0.4108    3.0483   19.9276  
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -11.57376    5.40141  -2.143   0.0326 *  
## SexM            6.35204    0.59877  10.609  < 2e-16 ***
## poly(Age, 2)1 -67.11159    5.37940 -12.476  < 2e-16 ***
## poly(Age, 2)2 -21.83199    5.18266  -4.213 2.99e-05 ***
## Height          0.22215    0.03274   6.785 3.25e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 25.83874)
## 
##     Null deviance: 33355  on 510  degrees of freedom
## Residual deviance: 13074  on 506  degrees of freedom
## AIC: 3118.8
## 
## Number of Fisher Scoring iterations: 2

Có rất nhiều algorithm có thể được dùng để tạo ra quy luật (mô hình hồi quy) mà ta cần, từ đơn giản như một mô hình tuyến tính, đến những công cụ kì lạ khác như gamlss, random forest, gaussian process hay thậm chí neural network; nhưng các bạn chỉ cần hình dung về chúng như một hàm F, cho phép ước tính kết quả y_hat từ dữ liệu đầu vào (X):

fi ~ F(X)

Dĩ nhiên, chúng ta mong muốn rằng kết quả (fi) phải gần với quan sát thực tế (yi) nhất có thể, càng gần càng tốt. Do đó, một cách tự nhiên, ta nhìn vào khác biệt giữa giá trị tiên lượng và giá trị thực tế và gọi đó là sai số :

err = fi - yi

Đây chính là nền tảng của mọi tiêu chí kiểm định mô hình mà ta sẽ thấy trong phần tiếp theo:

Nếu nghiên cứu có mục tiêu diễn dịch, ta thường kiểm định mô hình trên chính dữ liệu gốc (trainset) hoặc có sử dụng một hình thức tái chọn mẫu như bootstrap, cross-validation… , nhưng cho mục tiêu tiên lượng, ta bắt buộc phải kiểm định mô hình trên một quần thể độc lập khác (testset).

2 Kiểm định trực quan bằng biểu đồ

Trước hết, ta hoàn toàn có thể kiểm tra phẩm chất của mô hình một cách trực quan bằng hình ảnh mà không cần đến các chỉ số thống kê. Cách làm này sẽ giúp bạn gây ấn tượng với khán giả/độc giả khi thuyết trình, in poster hay công bố báo chí, và thực ra, chúng chuyển nhiều thông tin hơn bất cứ một con số đơn độc nào. Do đó, Nhi sẽ nói về phần này trước:

Nhi tạo ra 1 dataframe mới từ testset, với 3 cột truth = giá trị thực, predicted = giá trị từ mô hình, và error = sai biệt giữa 2 giá trị này.

pdf = data_frame(sex=testset$Sex,
                 age = testset$Age,
                 height = testset$Height,
                 truth = testset$DLCO,
                 predicted = predict(model,testset),
                 error= predicted - truth)

head(pdf)

## # A tibble: 6 x 6
##   sex     age height truth predicted  error
##   <fct> <int>  <dbl> <dbl>     <dbl>  <dbl>
## 1 F        22   165   26.3      27.8  1.53 
## 2 F        22   169.  27.9      28.6  0.732
## 3 F        26   149   19.2      24.2  4.97 
## 4 F        26   159   24.8      26.4  1.60 
## 5 F        28   171.  27.1      29.0  1.92 
## 6 F        30   173   30.2      29.3 -0.901

2.1 So sánh mật độ phân bố giữa thực tế và tiên lượng:

Tất cả các mô hình hồi quy đều hoạt động dựa trên một giả định về phân bố của biến kết quả ngẫu nhiên, nếu mô hình chính xác,hình ảnh phân bố kết quả của nó sẽ đồng dạng và trùng lắp với phân bố thực tế của đại lượng mà ta muốn ước tính: Tuy nhiên, hầu hết mô hình chỉ có khả năng giảm thiểu được sai biệt giữa fi và vị trí trung tâm của y, do đó kết quả tiên lượng thường tập trung quanh trung vị hay trung bình của y, nhưng không bao giờ chồng lắp một cách hoàn hảo và tái hiện được phân bố thực.

Trên biểu đồ, phẩm chất của mô hình được đánh giá bằng mức độ đồng dạng giữa 2 phân bố, vị trí trung tâm, trong khi sai sót của mô hình thể hiện qua phần diện tích không chồng lắp.

pdf%>%gather(truth,predicted,key="Y",value="DLCO")%>%
  ggplot(aes(x=DLCO,fill=Y))+
  geom_density(alpha=0.3)+
  scale_fill_manual(values=c("blue","red"))+
  facet_wrap(~sex,ncol=1)+
  theme_bw()

Một hình thức đơn giản hơn là biểu đồ boxplot, dạng biểu đồ này không thể hiện được tất cả hình ảnh của phân bố, nhưng nó trình bày rất tốt những mốc so sánh như: trung vị, min, max, tứ phân vị, và những giá trị outliers. Khi đặt song song 2 boxplot của thực tế và tiên lượng, ta có thể hình dung về giới hạn ứng dụng của mô hình, và khuynh hướng của nó (ước lượng quá cao, quá thấp, hay phù hợp):

pdf%>%gather(truth,predicted,key="Y",value="DLCO")%>%
  ggplot(aes(Y,DLCO,fill=Y,col=Y))+
  geom_jitter(alpha=0.5)+
  geom_boxplot(alpha=0.5)+
  coord_flip()+
  stat_compare_means(method="t.test",paired = TRUE,label.y = 45)+
  scale_color_manual(values=c("blue","red"))+
  scale_fill_manual(values=c("blue","red"))+
  facet_wrap(~sex,ncol=1)+
  theme_bw()

2.2 Tương quan tuyến tính giữa giá trị thực tế và tiên lượng:

Áp dụng một biểu đồ tán xạ và một đồ thị tuyến tính giữa 2 vectors: predicted và truth cho phép khảo sát tương quan tuyến tính giữa kết quả của mô hình và quan sát thực tế.Nếu bạn thích, có thể ghi chú thêm giá trị của hệ số tương quan lên biểu đồ.

pdf%>%ggplot(aes(truth,predicted,color=sex))+
  geom_jitter(alpha=0.5)+
  geom_smooth(aes(fill=sex),alpha=0.2,method="lm")+
  stat_cor(method="spearman",label.x = 15, label.y = 45)+
  scale_color_manual(values=c("blue","red"))+
  scale_fill_manual(values=c("blue","red"))+
  facet_wrap(~sex,ncol=2)+
  theme_bw()

2.3 Khảo sát trực tiếp sai biệt giữa thực tế và tiên lượng

Một hướng đi khác đó là ta chỉ quan tâm đến sai số của mô hình, thí dụ biểu đồ mật độ phân bố có thể được dùng để kiểm tra các giả định về sai số : phân bố chuẩn, trung bình = 0, đối xứng…

pdf%>%
  ggplot(aes(x=error,fill=sex))+
  geom_density(alpha=0.5)+
  scale_fill_manual(values=c("pink","skyblue"))+
  theme_bw()

Một biểu đồ tán xạ trình bày giá trị tuyệt đối của sai số, và liên hệ với thang đo của kết quả cần tiên lượng: Hình vẽ cho biết khuynh hướng sai biệt của mô hình (over hay underestimation), và sai biệt này là ngẫu nhiên hay có hệ thống: (Dạng biểu đồ này tương tự như Bland-Altman plot)

pdf%>%mutate(est=if_else(.$error>0,"Over","Under"))%>%
  ggplot(aes(x=truth,y=error,col=est))+
  geom_jitter(alpha=0.5)+
  geom_hline(yintercept = 0,linetype=2,col="red")+
  scale_color_manual(values=c("red","blue"))+
  facet_wrap(~sex,ncol=2)+
  theme_bw()

2.4 Kiểm tra tính hợp lý của nội dung mô hình

Bạn có thể sử dụng dạng biểu đồ sau đây để kiểm tra đồng thời : liệu mô hình có chính xác không ? Quy luật bên trong của mô hình có phù hợp với quan sát thực tế hay không ? Thí dụ, mô hình hiện thời sử dụng một hàm đa thức bậc 2 để diễn tả sự biến thiên của DLCO tùy theo tuổi - ta muốn kiểm tra xem quy luật này có áp dụng được cho một quần thể độc lập hay không ?

pdf%>%gather(truth,predicted,key="Y",value="DLCO")%>%ggplot()+
  geom_jitter(aes(age,DLCO),alpha=0.5)+
  geom_smooth(aes(age,DLCO,col=Y,fill=Y),alpha=0.2)+
  scale_color_manual(values=c("blue","red"))+
  scale_fill_manual(values=c("blue","red"))+
  facet_wrap(~sex,ncol=2)+
  ggtitle("Test_set")+
  theme_bw()

Cho bài toán diễn dịch, bạn có thể áp dụng biểu đồ cho dữ liệu gốc (trainset):

trainset%>%mutate(truth=trainset$DLCO,predicted=predict(model,trainset))%>%
  gather(truth,predicted,key="Y",value="DLCO")%>%ggplot()+
  geom_jitter(aes(Age,DLCO),alpha=0.2)+
  geom_smooth(aes(Age,DLCO,col=Y,fill=Y),alpha=0.5)+
  scale_color_manual(values=c("blue","red"))+
  scale_fill_manual(values=c("blue","red"))+
  facet_wrap(~Sex,ncol=2)+
  ggtitle("Train_set")+
  theme_bw()

## Warning: attributes are not identical across measure variables;
## they will be dropped

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Bạn cũng có thể thay biểu đồ tuyến kí bằng 2D density plot:

pdf%>%gather(truth,predicted,key="Y",value="DLCO")%>%
  ggplot(aes(age,DLCO))+
  stat_density_2d(geom = "polygon", aes(col = Y,fill=Y, alpha = ..level..))+
  geom_jitter(alpha=0.5)+
  scale_color_manual(values=c("blue4","red4"))+
  scale_fill_manual(values=c("blue","red"))+
  facet_wrap(~sex,ncol=2)+
  theme_bw()

3 Các tiêu chí kiểm định mô hình hồi quy

Trong phần tiếp theo, Nhi sẽ liệt kê hơn 14 tiêu chí có thể sử dụng để kiểm tra phẩm chất của mô hình hồi quy. Như đã nói ở trên, tất cả các tiêu chí này đều nhắm đến mục tiêu đo lường kích thước sai biệt giữa thực tế (truth, yi) và tiên lượng (predicted, fi).

3.1 Nhóm 1: Các tiêu chí dựa vào giá trị sai biệt tuyệt đối

Trong nhóm này, ta có:

1. MAE : Trung bình của sai biệt tuyệt đối

\[MAE = \frac{\sum abs(f_{i}-y_{i})}{n}\]

2. MEDAE: Trung vị của sai biệt tuyệt đối

\[MEDAE = Median(abs(f_{i}-y_{i}))\]

3. SAE: Tổng sai biệt tuyệt đối

\[SAE = \sum_{i}^{n}abs(f_{i}-y_{i})\]

4. MAPE Mean absolute percentage error

\[MAPE = \frac{\sum {(\frac{abs(y_{i}-f_{i})}{y_{i}})}}{n}\]

4 tiêu chí này rất dễ hiểu, cả 3 đều đo lường sai biệt giữa thực tế (yi) và tiên lượng của mô hình (fi),giá trị tuyệt đối được dùng để tránh sai lầm trong trường hợp mô hình đồng thờicó nguy cơ đánh giá quá cao và quá thấp, dẫn đến việc sai số > 0 và <0 triệt tiêu lẫn nhau.

Tiêu chí thứ 4: MAPE đo lường sai biệt theo tỉ lệ % , dùng cho những trường hợp mà biến kết quả có đơn vị quá thấp hoặc quá cao

3.2 Nhóm 2:Các tiêu chí dựa vào bình phương sai số

Các tiêu chí trong nhóm này không sử dụng hàm abs, nhưng dựa vào bình phương của sai số; ta lần lượt có:

5. MSE : trung bình bình phương sai số

\[MSE = \frac{\sum (f_{i}-y_{i})^{2}}{N}\]

6. MEDSE: Trung vị bình phương sai số

\[MEDSE = Median ((f_{i}-y_{i})^{2})\]

7. SSE: tổng bình phương sai số

\[SSE = \sum (f_{i}-y_{i})^{2}\]

8. RMSE : Căn bậc 2 của trung bình bình phương sai số

\[RMSE = \sqrt{\frac{\sum_{i}^{n} (f_{i}-y_{i})^{2}}{n}}\]

9. MSLE Mean squared logarithmic error

\[MSLE = \frac{\sum (log(f_{i}+1)- log(y_{i}+1))^{2}}{n}\]

10. RMSLE: Root mean squared logarithmic errors

\[RMSLE = \sqrt{\frac{\sum (log(f_{i}+1)- log(y_{i}+1))^{2}}{n}}\] ## Nhóm 3: Khảo sát sai biệt tương đối của mô hình

2 tiêu chí trong nhóm 3 này không đo lường sai biệt của mô hình một cách độc lập và tuyệt đối, nhưng đánh giá tương đối tỉ lệ giữa sai biệt của mô hình và sai biệt nội tại của biến kết quả bên trong quần thể

11. RRSE Root relative squared error

\[RRSE = \sqrt{\frac{\sum(f_{i}-y_{i})^{2}}{\sum(y_{i} - \bar{y})^{2}}}\]

12. RAE : Relative absolute error

\[RAE = \frac{\sum (abs(f_{i}-y_{i}))}{\sum (abs(y_{i}-\bar{y}))}\]

3.3 Nhóm 4 :

13. Rsq: hệ số xác định : coefficient of determination (R2)

Bản thân hệ số này không cung cấp thông tin về độ chính xác của mô hình nhưng cho biết mô hình giải thích được bao nhiêu phần phương sai của biến kết quả trong mẫu. Thông tin này khá quan trọng cho mô hình diễn dịch.

Hệ số R2 có liên hệ với những thành phần phương sai mà ta đã biết trong ANOVA cổ điển, bao gồm tổng phương sai (SSE), phương sai của mô hình (mà mô hình có khả năng giải thích, MSE) và phương sai tồn lưu (residual sum of squares)

\[R^{2} = \frac{SSM}{SSE} = 1 - \frac{SSR}{SSE} = \frac{\sum (f_{i}-\bar{y})^{2}}{\sum (y_{i} - \bar{y})^{2}} = 1- \frac{\sum (y_{i}-f_{i})^{2}}{\sum (y_{i} - \bar{y})^{2}}\]

Hình: SSE = vùng màu đỏ, SSR = vùng màu xanh, f = mô hình

3.4 Nhóm 5 : Các tiêu chí đánh giá mối tương quan tuyến tính giữa giá trị thực và tiên lượng:

Bản chất của 3 tiêu chí này chính là 3 phương pháp tính hệ số tương quan trong thống kê cổ điển mà ta từng biết, gồm

14. Kendall’s Tau (thứ hạng)

\[\tau = \frac{số cặp trùng hợp - số cặp không trùng hợp}{n(n-1)/2}\]

15. Rho của Spearman (phi tham số)

\[\rho = \frac{cov(rankY,rankY_{hat})}{\sigma_{rankY} \sigma_{rankY_{hat}}}\]

16. Pearson’r

\[r = \frac{cov(rankY,rankY_{hat})}{\sigma_{Y} \sigma_{Y_{hat}}}\]

4 Viết một hàm tính thủ công tất cả tiêu chí kiểm định mô hình hồi quy

scoring = function(methods,truth,predicted){
  
  SSE = function(truth, predicted) {
    sum((predicted - truth)^2)
  }
  
  # MSE
  
  MSE = function(truth, predicted) {
    mean((predicted - truth)^2)
  }
  
  # RMSE
  
  RMSE = function(truth, predicted) {
    sqrt(MSE(truth, predicted))
  }
  
  # MEDSE
  
  MEDSE = function(truth, predicted) {
    median((predicted - truth)^2)
  }
  
  # SAE
  
  SAE = function(truth, predicted) {
    sum(abs(predicted - truth))
  }
  
  # MAE
  
  MAE = function(truth, predicted) {
    mean(abs(predicted - truth))
  }
  
  # MEDAE
  
  MEDAE = function(truth, predicted) {
    median(abs(predicted - truth))
  }
  
  # RSQ
  
  RSQ = function(truth,  predicted) {
    rss = SSE(truth,  predicted)
    ess = sum((truth - mean(truth))^2L)
    if (ess == 0){
      warning("Error: all truth values are equal")
      return(NA_real_)
    }
    1 - rss / ess
  }

  # RRSE Root relative squared error
  
  RRSE = function(truth, predicted){
    tss = sum((truth - mean(truth))^2L)
    if (tss == 0){
      warning("Error: all truth values are equal.")
      return(NA_real_)
    }
    sqrt(SSE(truth, predicted) / tss)
  }
  
  # RAE : Relative absolute error
  
  RAE = function(truth, predicted){
    meanad = sum(abs(truth - mean(truth)))
    if (meanad == 0){
      warning("Error:  all truth values are equal.")
      return(NA_real_)
    }
    return(SAE(truth, predicted) / meanad)
  }
  
  # MAPE Mean absolute percentage error
  
  MAPE = function(truth, predicted){
    if (any(truth == 0)){
      warning("Error: truth value is equal to 0.")
      return(NA_real_)
    }
    return(mean(abs((truth - predicted) / truth)))
  }
  
  # Mean squared logarithmic error
  
  MSLE = function(truth, predicted) {
    if (any(truth < -1))
      stop("All truth values must be greater or equal -1")
    if (any(predicted < -1))
      stop("All predicted values must be greater or equal -1")
    
    mean((log(predicted + 1) - log(truth + 1))^2)
  }
  
  # Root mean squared logarithmic error
  
  RMSLE = function(truth, predicted) {
    sqrt(MSLE(truth, predicted))
  }
  
  # Kendall Tau
  
  KendallTau = function(truth, predicted) {
    cor(truth, predicted, use = "na.or.complete", method = "kendall")
  }
  
  # rho
  
  SpearmanRho = function(truth, predicted) {
    cor(truth, predicted, use = "na.or.complete", method = "spearman")
  }
  
  # Pearson r
  
  PearsonR = function(truth, predicted) {
    cor(truth, predicted, use = "na.or.complete", method = "pearson")
  }
  
  scores = data_frame(Score=methods,Value=rep(NA,length(methods)))
  
  for (i in 1:length(methods)){
    scoring_func = get(methods[i])
    scores$Value[i]= scoring_func(truth,predicted)
  }
  
  return(scores)
}

Áp dụng hàm này cho kết quả trên testset: Chú ý danh pháp viết tắt của các tiêu chí:

scores = scoring(methods=c('MAE','MAPE','MEDAE','MEDSE','MSE','MSLE',
                  'RAE','RMSE','RMSLE','RRSE','RSQ','SAE','SSE',
                  'PearsonR','KendallTau','SpearmanRho'),
        truth=pdf$truth,
        predicted=pdf$predicted)

knitr::kable(scores)

Score	Value
MAE	4.0078789
MAPE	0.1316461
MEDAE	3.4934013
MEDSE	12.2039131
MSE	25.9979765
MSLE	0.0234722
RAE	0.6218403
RMSE	5.0988211
RMSLE	0.1532062
RRSE	0.6241336
RSQ	0.6104573
SAE	504.9927456
SSE	3275.7450392
PearsonR	0.7832422
KendallTau	0.6043231
SpearmanRho	0.7972955

Bài thực hành đến đây là hết. Tạm biệt các bạn và hẹn gặp lại :)