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: //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["//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.8Nhi 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: 2Có 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].
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 ## ## 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.901Tấ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]:
Á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[]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[]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[]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].
Trong nhóm này, ta có:
- MAE : Trung bình của sai biệt tuyệt đối
\[MAE = \frac{\sum abs[f_{i}-y_{i}]}{n}\]
- MEDAE: Trung vị của sai biệt tuyệt đối
\[MEDAE = Median[abs[f_{i}-y_{i}]]\]
- SAE: Tổng sai biệt tuyệt đối
\[SAE = \sum_{i}^{n}abs[f_{i}-y_{i}]\]
- 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à