Đánh giá mô hình phân loại Multiclass

1 Giới thiệu

Trong 2 bài trước, Nhi đã hướng dẫn các bạn cách tính thủ công tất cả các tiêu chí nhằm kiểm định mô hình phân loại nhị phân : https://rpubs.com/lengockhanhi/347941 và mô hình hồi quy : https://rpubs.com/lengockhanhi/445130 .Bài hôm nay Nhi sẽ tiếp tục giới thiệu một số tiêu chí đặc biệt dành cho mô hình phân loại multiclass.

Với bài toán Multiclass classification, mỗi trường hợp được dán 1 nhãn duy nhất trong tập hợp n nhãn giá trị được quy định từ trước. Như vậy các nhãn giá trị sẽ loại trừ lẫn nhau (khác với bài toán multilabel, khi một cá thể có thể nhận đồng thời nhiều nhãn).

2 Thí dụ minh họa

Trong thí dụ minh họa sau, ta có 400 bệnh nhân được xét nghiệm để đo nồng độ của 3 loại marker khác nhau là A,B,C. Mục tiêu của chúng ta là xây dựng một quy luật chẩn đoán để phân chia độ năng của một căn bệnh Y, với 4 mức độ gồm: Bình thường (Normal), bệnh nhẹ (Mild), bệnh trung bình (Moderate) và bệnh nặng (Severity), sử dụng chỉ số BMI và 3 biomarker A,B,C nói trên.

library(tidyverse)

df = read.csv("Multiclass.csv",sep = ',',dec = ".")
knitr::kable(head(df))

BMI	Marker_A	Marker_B	Marker_C	Severity
32.32	34.991	14.480	42.724	Moderate
31.31	86.754	89.415	102.693	Severe
27.40	54.385	74.213	63.430	Moderate
26.35	41.457	26.850	54.493	Moderate
23.24	16.252	12.590	33.832	Mild
27.70	45.826	27.428	55.275	Moderate

severity_col = c("#ffad2b","#ff2b55","#b1ef45","#a045ef")

df%>%gather(Marker_A,Marker_B,Marker_C,key="Marker",value="Score")%>%
  ggplot()+
  geom_density(aes(x=Score,fill=Severity),alpha=0.5)+
  theme_bw(8)+
  scale_fill_manual(values=severity_col)+
  scale_x_continuous(breaks=c(0,10,20,seq(from=30,to=300,by=20)))+
  facet_wrap(~Marker,
             ncol=1,scales = "free_y")

Trước kia, bài toán phân chia độ nặng của bệnh lý thường chỉ dựa vào 1 chỉ số lâm sàng hay cận lâm sàng duy nhất, thí dụ một marker X. Cách làm đơn giản nhất là phân chia thang đo của X thành nhiều phần dựa vào 1 quy luật phân bố nào đó, thí dụ tứ phân vị, trung vị, hoặc những ngưỡng cắt tùy tiện khác.

Nếu làm nghiêm túc hơn, người ta có thể dùng một mô hình multinomial logistic, hoặc dựa vào kết cục lâm sàng/tiên lượng sống còn … để xác định các ngưỡng cắt phù hợp cho từng độ nặng.Tuy nhiên, nếu quan sát phân bố của 3 marker trên biểu đồ, ta sẽ thấy rằng bất cứ biến số nào khi đứng một mình cũng không cho phép phân chia hiệu quả, ngay cả giữa 2 lớp cũng có sự chồng lắp đáng kể.

Ở thời đại mới, khi machine learning trở nên phổ biến, bài toán phân chia nhiều nhóm trở nên hết sức đơn giản, hơn nữa ta còn có thể đưa nhiều thông tin hơn vào mô hình. Một số algorithm xử lý bài toán multiclass một cách tự nhiên và chính xác, bao gồm mô hình cây, random forest hay neural network…

Để minh họa, Nhi sử dụng một mô hình RF, luyện trên 70% dữ liệu và kiểm định trên 30% còn lại. Nhi sẽ lần lượt giới thiệu các tiêu chí khác nhau, được dùng riêng cho loại mô hình multiclass này.

set.seed(1234)
idtest=caret::createDataPartition(y=df$Severity, p=0.7,list=FALSE)

trainset=df[idtest,]
testset=df[-idtest,]

table(trainset$Severity)

## 
##     Mild Moderate   Normal   Severe 
##       77       81       47       76

table(testset$Severity)

## 
##     Mild Moderate   Normal   Severe 
##       33       34       20       32

library(randomForest)

rfmod=randomForest(Severity~.,
                   data=trainset,
                   ntree=500,
                   mtry=2)

rfmod

## 
## Call:
##  randomForest(formula = Severity ~ ., data = trainset, ntree = 500,      mtry = 2) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 30.96%
## Confusion matrix:
##          Mild Moderate Normal Severe class.error
## Mild       50       19      8      0   0.3506494
## Moderate   24       42      0     15   0.4814815
## Normal     10        0     37      0   0.2127660
## Severe      0       11      0     65   0.1447368

3 Confusion matrix multiclass

Đầu tiên, ta áp dụng mô hình cho testset, kết quả prediction gồm probability cho mỗi nhãn giá trị trong 4 nhãn, và kết quả phân loại (response).

predict(rfmod,newdata=testset,type="prob")%>%
  as.data.frame()%>%
  mutate(response=predict(rfmod,newdata=testset))->pred_df

knitr::kable(head(pred_df))

Mild	Moderate	Normal	Severe	response
0.502	0.480	0.018	0.000	Mild
0.048	0.514	0.000	0.438	Moderate
0.152	0.814	0.000	0.034	Moderate
0.600	0.398	0.000	0.002	Mild
0.076	0.360	0.000	0.564	Severe
0.026	0.786	0.002	0.186	Moderate

Cũng như bài toán nhị phân, cách kiểm định đơn giản nhất là dựng một confusion matrix. Cấu trúc của matrix này cũng tương tự, chỉ khác về kích thước với nhiều hàng và cột hơn. Độ chính xác của mô hình được thể hiện qua đường chéo của matrix:

truth = testset$Severity
response = pred_df$response
pred = pred_df[-5]

conf.mat = table(truth,response)

conf.mat%>%knitr::kable()

	Mild	Moderate	Normal	Severe
Mild	25	6	2	0
Moderate	5	23	2	4
Normal	3	0	17	0
Severe	0	8	0	24

Quan sát cho thấy mô hình tương đối chính xác, nó hoạt động tốt cho 2 nhãn Normal và Severe tuy nhiên mô hình nhầm lẫn giữa 2 nhóm Mild và Moderate (hiện tượng này cũng tương tự như khi ta dùng ROC curve để xác định ngưỡng cắt cho 1 biến số duy nhất, 2 đầu thang đo và những giá trị cực điểm lúc nào cũng dễ dàng được nhận ra hơn các vùng trung gian)

4 Độ chính xác sau cân bằng

Balanced accuracy : được định nghĩa là trung bình giữa tỉ lệ xác định đúng (TPR) và tỉ lệ loại trừ đúng (TNR). Như vậy, có thể tính BAC như sau:

# bac

diag(conf.mat) 

##     Mild Moderate   Normal   Severe 
##       25       23       17       24

table(truth, truth)

##           truth
## truth      Mild Moderate Normal Severe
##   Mild       33        0      0      0
##   Moderate    0       34      0      0
##   Normal      0        0     20      0
##   Severe      0        0      0     32

bac = mean(diag(table(truth, response) / table(truth, truth)))

bac

## [1] 0.7585116

Kết quả cho ra BAC = 0.758

Diễn giải: Giá trị tốt nhất của BAC là 1, xấu nhất =0.

5 Tỉ lệ nhầm lẫn sau cân bằng

Balanced error rate (BER) có ý nghĩa ngược lại với BAC, nó là tỉ lệ nhầm lẫn trung bình cho tất cả các nhãn (bài toán nhị phân có 2 nhãn, ở đây ta có 4 nhãn), BER có thể được tính thủ công từ 2 vector truth và response, nhưng cũng có thể tính nhanh bằng công thức BER = 1 - BAC.

Do BER nghịch với BAC, nên giá trị tối ưu của nó =0,giá trị xấu nhất =1

ber = mean(diag(1 - (table(truth, response) / table(truth, truth))))

rbind(bac,ber,1-bac)%>%knitr::kable()

bac	0.7585116
ber	0.2414884
	0.2414884

6 Logloss

Logloss được định nghĩa bằng công thức: -mean(log(p_i)),với p_i là xác suất của nhãn giá trị thực sự ở trường hợp thứ i,

Lưu ý là logloss càng cao càng xấu, giá trị tối ưu của nó =0, nhưng không có giá trị tối đa, do đó khi tính logloss ta phải đặt ra 1 ngưỡng cho xác suất tối thiểu, thí dụ 1e-20, để nó không thể = 0 (không thể tính log của 0).

# logloss

pred_mat = pred

pred_mat[pred_mat > 1 - (1e-20)] = 1 - (1e-20)
pred_mat[pred_mat < (1e-20)] = 1e-20
  
c_truth = match(as.character(truth), colnames(pred_mat))
inds = cbind(BBmisc::seq_row(pred_mat), c_truth)
  
logloss = -1 * mean(log(pred[inds]))

logloss

## [1] 0.6772009

7 Kappa

Hệ số Kappa của Cohen là một trị số thống kê dùng để đánh giá mức độ tương hợp giữa 2 vector, được xác định theo công thức: kappa = 1 - (1 - p0) / (1 - pe).Với p0 = tần suất tương hợp trên thực tế và pe = tần suất tương hợp giả định (nếu 2 vector là độc lập với nhau)

Ta có thể tính thủ công kappa từ 2 vector truth và response như sau:

# conf mat
  
conf.mat = table(truth, response)
conf.mat = conf.mat / sum(conf.mat)
p0 = sum(diag(conf.mat))
rowsum = rowSums(conf.mat)
colsum = colSums(conf.mat)
pe = sum(rowsum * colsum) / sum(conf.mat)^2
# kappa
kappa = 1 - (1 - p0) / (1 - pe)

round(cbind(p0,pe,kappa),5) %>%knitr::kable()

p0	pe	kappa
0.7479	0.25867	0.65994

8 Kappa hiệu chỉnh trọng số

Đặc biệt cho những trường hợp mà vector kết quả có tính thứ bậc (thí dụ độ nặng của bệnh lý), ta có hệ số kappa hiểu chỉnh bởi một matrix trọng số

Dưới đây là cách tính thủ công kappa với trọng số:

#Wkappa
conf.mat = table(truth, response)
conf.mat = conf.mat / sum(conf.mat)
# get expected probs under independence
rowsum = rowSums(conf.mat)
colsum = colSums(conf.mat)
expected.mat = rowsum %*% t(colsum)

# get weights
class.values = seq_along(levels(truth)) - 1L
k_weights = outer(class.values, class.values, FUN = function(x, y) (x - y)^2)

# calculate weighted kappa
wkappa = 1 - sum(k_weights * conf.mat) / sum(k_weights * expected.mat)


k_weights 

##      [,1] [,2] [,3] [,4]
## [1,]    0    1    4    9
## [2,]    1    0    1    4
## [3,]    4    1    0    1
## [4,]    9    4    1    0

round(rbind(kappa,wkappa),5) %>%knitr::kable()

kappa	0.65994
wkappa	0.73811

Cả 2 hệ số kappa này dao động từ -1 đến 1, giá trị tối ưu =1, xấu nhất =-1

9 Chỉ số Brier multiclass

Brier score là một tiêu chí đo lường nguy cơ sai lầm của mô hình, nó dao động từ 0 đến 2, giá trị tốt nhất =0 và xấu nhất =2. brier score được tính bằng công thức:

Brier score = (1/n) sum_i sum_j (y_ij - p_ij)^2, với y_ij = 1 nếu trường hợp i nhận nhãn j (trái lại, sẽ =0), còn p_ij là xác suất tiên lượng của mô hình cho nhãn j tại trường hợp i

  # multiclass brier
  
truth_f = factor(truth, levels = colnames(pred))
mat01 = mlr::createDummyFeatures(truth_f)
mc_brier = mean(rowSums((pred - mat01)^2))
  
mc_brier

## [1] 0.3735652

Brier score của mô hình = 0.37, một giá trị không tệ.

10 AU1P

Diện tích dưới đường cong ROC là một tiêu chí đo lường hiệu năng của một quy luật nhị phân, thí dụ sự phân biệt giữa 2 nhãn A và B. Một bài toán multiclass có C nhãn, do đó ta có thể dựng lên tất cả c(c-1) đường cong ROC để phân biệt từng cặp nhãn với nhau.

Tiêu chí AU1P là giá trị AUC trung bình của c(c - 1) đường cong ROC,nhưng có xét thêm tỉ lệ phân bố của mỗi class trong dữ liệu.

colAUC = function(samples, truth, maximum = TRUE) {
  y = as.factor(truth)
  X = as.matrix(samples)
  if (nrow(X) == 1)
    X = t(X)
  nr = nrow(X)
  nc = ncol(X)
  ny = table(y)
  ul = as.factor(rownames(ny))
  nl = length(ny)
  if (nl <= 1)
    stop("colAUC: List of labels 'y' have to contain at least 2 class labels.")
  if (!is.numeric(X))
    stop("colAUC: 'X' must be numeric")
  if (nr != length(y))
    stop("colAUC: length(y) and nrow(X) must be the same")
  per = t(utils::combn(1:nl, 2))
  np = nrow(per)
  auc = matrix(0.5, np, nc)
  rownames(auc) = paste(ul[per[, 1]], " vs. ", ul[per[, 2]], sep = "")
  colnames(auc) = colnames(X)
  # Wilcoxon AUC
  idxl = vector(mode = "list", length = nl)
  for (i in 1:nl) idxl[[i]] = which(y == ul[i])
  for (j in 1:nc) {
    for (i in 1:np) {
      c1 = per[i, 1]
      c2 = per[i, 2]
      n1 = as.numeric(ny[c1])
      n2 = as.numeric(ny[c2])
      if (n1 > 0 & n2 > 0) {
        r = rank(c(X[idxl[[c1]], j], X[idxl[[c2]], j]))
        auc[i, j] = (sum(r[1:n1]) - n1 * (n1 + 1) / 2) / (n1 * n2)
      }
    }
  }
  if (maximum == TRUE) {
    auc = pmax(auc, 1 - auc)
  }
  return(auc)
}

# au1p
m = colAUC(as.matrix(pred), truth)

weights = table(truth) / length(truth)

m = m * matrix(rep(weights, each = nrow(m)), ncol = length(weights))

m%>%knitr::kable()

	Mild	Moderate	Normal	Severe
Mild vs. Moderate	0.2449333	0.2291826	0.1398313	0.2224120
Mild vs. Normal	0.2296218	0.2666667	0.1591546	0.1861981
Mild vs. Severe	0.2723214	0.1904762	0.1567673	0.2647059
Moderate vs. Normal	0.1598616	0.2789916	0.1643599	0.2479486
Moderate vs. Severe	0.2346183	0.2351628	0.1190219	0.2429560
Normal vs. Severe	0.2214154	0.2174107	0.1680672	0.2686975

weights%>%knitr::kable()

truth	Freq
Mild	0.2773109
Moderate	0.2857143
Normal	0.1680672
Severe	0.2689076

Tính AU1P:

c = c(combn(1:nlevels(truth), 2))

au1p = sum(m[cbind(rep(seq_len(nrow(m)), each = 2), c)]) / (nlevels(truth) - 1)

au1p

## [1] 0.9193849

11 AU1U

Tương tự, chỉ số AU1U cũng tính tring bình AUC của c(c - 1) đường cong ROC nhằm phân biệt bắt cặp tuần tự giữa các nhãn, nhưng khác với AU1P, AU1U giả định các nhãn có phân bố đồng dạng.

#au1u
  
m = colAUC(as.matrix(pred), truth)

m%>%knitr::kable()

	Mild	Moderate	Normal	Severe
Mild vs. Moderate	0.8832442	0.8021390	0.8319964	0.8270945
Mild vs. Normal	0.8280303	0.9333333	0.9469697	0.6924242
Mild vs. Severe	0.9820076	0.6666667	0.9327652	0.9843750
Moderate vs. Normal	0.5764706	0.9764706	0.9779412	0.9220588
Moderate vs. Severe	0.8460478	0.8230699	0.7081801	0.9034926
Normal vs. Severe	0.7984375	0.7609375	1.0000000	0.9992188

c = c(combn(1:nlevels(truth), 2))
au1u = mean(m[cbind(rep(seq_len(nrow(m)), each = 2), c)])

au1u

## [1] 0.9255799

12 AUNP

Tiêu chí AUNP cũng tính AUC trung bình, nhưng của 4 đường cong ROC tương ứng với 4 quy luật nhị phân, nhằm phân biệt 1 nhãn với 3 nhãn còn lại, có xét đến tỉ lệ phân bố thực tế của các nhãn.

aunp=sum(BBmisc::vnapply(1:nlevels(truth),function(i) mean(truth == levels(truth)[i]) * (colAUC(as.matrix(pred)[, i], truth == levels(truth)[i]))))

aunp

## [1] 0.9156638

13 AUNU

Tương tự, AUNU cũng xét 4 quy luật nhị phân như AUNP, nhưng với giả định các phân nhóm có phân bố đồng dạng

aunu=mean(BBmisc::vnapply(1:nlevels(truth),function(i) colAUC(as.matrix(pred)[, i],truth == levels(truth)[i])))

aunu

## [1] 0.9222792

14 Tổng kết

Bây giờ, ta tạo ra 1 hàm để tính đồng thời 10 tiêu chí kể trên từ matrix xác suất của C nhãn, vector nhãn thực tế và vector nhãn predict.

multiclass_perfomance = function(pred,truth,response){
  if (length(unique(truth)) != nlevels(truth)){
    warning("Error: Missing label(s) detected")
    return(NA_real_)
  }
  #aunu
  aunu=mean(BBmisc::vnapply(1:nlevels(truth),
                       function(i) colAUC(as.matrix(pred)[, i],
                                                   truth == levels(truth)[i])))
  # aunp
  aunp=sum(BBmisc::vnapply(1:nlevels(truth), 
                           function(i) mean(truth == levels(truth)[i]) * (colAUC(as.matrix(pred)[, i], 
                                                                                          truth == levels(truth)[i]))))
  # au1p
  m = colAUC(as.matrix(pred), truth)
  weights = table(truth) / length(truth)
  m = m * matrix(rep(weights, each = nrow(m)), ncol = length(weights))
  c = c(combn(1:nlevels(truth), 2))
  au1p = sum(m[cbind(rep(seq_len(nrow(m)), each = 2), c)]) / (nlevels(truth) - 1)
  
  #au1u
  
  m = colAUC(as.matrix(pred), truth)
  c = c(combn(1:nlevels(truth), 2))
  mean(m[cbind(rep(seq_len(nrow(m)), each = 2), c)])
  au1u = mean(m[cbind(rep(seq_len(nrow(m)), each = 2), c)])
  
  # conf mat
  
  conf.mat = table(truth, response)
  conf.mat = conf.mat / sum(conf.mat)
  p0 = sum(diag(conf.mat))
  rowsum = rowSums(conf.mat)
  colsum = colSums(conf.mat)
  pe = sum(rowsum * colsum) / sum(conf.mat)^2
  # kappa
  kappa = 1 - (1 - p0) / (1 - pe)
  
  #Wkappa
  
  conf.mat = table(truth, response)
  conf.mat = conf.mat / sum(conf.mat)
  # get expected probs under independence
  rowsum = rowSums(conf.mat)
  colsum = colSums(conf.mat)
  expected.mat = rowsum %*% t(colsum)
  
  # get weights
  class.values = seq_along(levels(truth)) - 1L
  k_weights = outer(class.values, class.values, FUN = function(x, y) (x - y)^2)
  
  # calculate weighted kappa
  wkappa = 1 - sum(k_weights * conf.mat) / sum(k_weights * expected.mat)
  
  # multiclass brier
  
  truth_f = factor(truth, levels = colnames(pred))
  mat01 = mlr::createDummyFeatures(truth_f)
  mc_brier = mean(rowSums((pred - mat01)^2))
  
  # bac
  bac = mean(diag(table(truth, response) / table(truth, truth)))
  
  # ber
  ber = mean(diag(1 - (table(truth, response) / table(truth, truth))))
  
  # logloss
  
  pred[pred > 1 - (1e-20)] = 1 - (1e-20)
  pred[pred < (1e-20)] = 1e-20
  
  c_truth = match(as.character(truth), colnames(pred))
  inds = cbind(BBmisc::seq_row(pred), c_truth)
  
  logloss = -1 * mean(log(pred[inds]))
  
  result = rbind(bac,ber,logloss,mc_brier,
                 aunu,aunp,au1u,au1p,
                 kappa,wkappa)
  
  return(result)
  }

Áp dụng hàm này để lấy kết quả của 10 tiêu chí.

multiclass_perfomance(truth=truth,
                      pred=pred,
                      response=response)%>%knitr::kable()

bac	0.7585116
ber	0.2414884
logloss	0.6772009
mc_brier	0.3735652
aunu	0.9222792
aunp	0.9156638
au1u	0.9255799
au1p	0.9193849
kappa	0.6599352
wkappa	0.7381062

Bài thực hành đến đây là kết thúc; hẹn gặp lại các bạn lần tới.