Hiểu về đường cong ROC và AUC dưới giả định Binormal

Phần II

1 Giới thiệu

Đường cong ROC (receiver operating characteristic) là một công cụ cơ bản để đánh giá hiệu suất của một chỉ dấu chẩn đoán.

Trong bài giảng này, chúng ta xây dựng trực giác từng bước, bắt đầu từ phân loại dựa trên ngưỡng và kết thúc với biểu diễn ROC dưới giả định binormal, cũng như cách hiểu xác suất của diện tích dưới đường cong (AUC).

Xét một chỉ dấu liên tục X được đo trên hai quần thể:

• nhóm không bệnh
  
• nhóm có bệnh

Giá trị X càng lớn thì càng có bằng chứng cho bệnh.

2 Phân loại theo ngưỡng

Với một ngưỡng r, ta định nghĩa quy tắc:

\text{dương tính nếu } X \ge r.

Khi đó, ta có hai đại lượng:

\mathrm{FPF}(r) = P(X \ge r \mid \text{không bệnh}),

\mathrm{TPF}(r) = P(X \ge r \mid \text{có bệnh}).

Mỗi ngưỡng r cho một điểm:

(\mathrm{FPF}(r), \mathrm{TPF}(r)).

3 Định nghĩa đường cong ROC

Đường cong ROC là tập hợp các điểm khi thay đổi r:

\mathrm{ROC} = \left\{ (\mathrm{FPF}(r), \mathrm{TPF}(r)) : r \in \mathbb R \right\}.

Đây là cách hiểu trực quan nhất: thay đổi ngưỡng và tính toán các xác suất.

4 ROC như một hàm theo false-positive rate

Ta có thể viết lại ROC theo biến

t = \mathrm{FPF}(r).

Định nghĩa

\mathrm{ROC}(t) = \mathrm{TPF}(r_t),

với r_t thỏa

\mathrm{FPF}(r_t) = t.

Khi đó,

\mathrm{ROC} = \left\{ (t, \mathrm{ROC}(t)) : t \in (0,1) \right\}.

Dạng này hữu ích khi cần tích phân.

5 Diện tích dưới đường cong ROC

Diện tích dưới đường cong là

\mathrm{AUC} = \int_0^1 \mathrm{ROC}(t)\,dt.

Nó biểu diễn giá trị trung bình của true-positive rate khi xét toàn bộ false-positive rate.

6 AUC dưới dạng xác suất

Ta sẽ chứng minh:

\mathrm{AUC} = P(X_1 > X_0),

trong đó:

• X_1: giá trị marker của một người bệnh được chọn ngẫu nhiên
  
• X_0: giá trị marker của một người không bệnh được chọn ngẫu nhiên

6.1 Suy luận từng bước

Bắt đầu từ

\mathrm{AUC} = \int_0^1 \mathrm{ROC}(t)\,dt.

Ta có

t = \mathrm{FPF}(r) = P(X_0 \ge r) = 1 - F_0(r),

và

\mathrm{ROC}(t) = \mathrm{TPF}(r) = P(X_1 \ge r).

Suy ra

\mathrm{AUC} = \int_{-\infty}^{+\infty} \mathrm{TPF}(r)\, dF_0(r).

Thay vào:

\mathrm{AUC} = \int_{-\infty}^{+\infty} P(X_1 \ge r)\, dF_0(r).

Ý nghĩa của dF_0(r) là ta đang lấy trung bình theo phân phối của X_0.

Do đó:

\mathrm{AUC} = E_{X_0}\bigl[ P(X_1 \ge X_0 \mid X_0) \bigr].

Vì X_0 và X_1 độc lập:

E_{X_0}\bigl[ P(X_1 \ge X_0 \mid X_0) \bigr] = P(X_1 \ge X_0).

Suy ra:

\mathrm{AUC} = P(X_1 > X_0).

6.2 Diễn giải hình học

Xét mặt phẳng 2 chiều:

• trục hoành: X_0
  
• trục tung: X_1

Đường X_1 = X_0 chia mặt phẳng thành hai vùng.

Vùng phía trên:

X_1 > X_0

Do đó:

\mathrm{AUC} = \text{xác suất nằm phía trên đường chéo}.

Đây là một cách biểu diễn khác với ROC, nhưng cùng một bản chất.

7 Mô hình Binormal

Giả sử:

X_0 \sim N(\mu_0, \sigma_0^2), \qquad X_1 \sim N(\mu_1, \sigma_1^2).

7.1 False-positive fraction

\mathrm{FPF}(r) = 1 - \Phi\left(\frac{r - \mu_0}{\sigma_0}\right).

Đặt

t = \mathrm{FPF}(r).

Gọi Z_t sao cho

1 - \Phi(Z_t) = t.

Suy ra:

r = \mu_0 + \sigma_0 Z_t.

7.2 True-positive fraction

\mathrm{TPF}(r) = 1 - \Phi\left(\frac{r - \mu_1}{\sigma_1}\right).

Thay vào:

\mathrm{TPF}(t) = 1 - \Phi\left( \frac{\mu_0 + \sigma_0 Z_t - \mu_1}{\sigma_1} \right).

Đặt

a = \frac{\mu_1 - \mu_0}{\sigma_1}, \qquad b = \frac{\sigma_0}{\sigma_1}.

Ta được:

\mathrm{ROC}(t) = 1 - \Phi(bZ_t - a).

7.3 Dạng cuối cùng

\left[ t,\, 1 - \Phi(bZ_t - a) \right], \quad t \in (0,1).

8 Diễn giải tham số

• a: độ tách biệt giữa hai nhóm
  
• b: so sánh độ biến thiên giữa hai nhóm

Nếu \sigma_0 = \sigma_1 thì b = 1.

9 Kết luận chính

• ROC được xây dựng bằng cách thay đổi ngưỡng
  
• AUC là diện tích dưới đường cong
  
• AUC cũng có thể hiểu là:

\mathrm{AUC} = P(X_1 > X_0)

• Dưới giả định binormal, ROC có dạng:

\mathrm{ROC}(t) = 1 - \Phi(bZ_t - a)

Các cách biểu diễn này là tương đương, nhưng cung cấp những góc nhìn khác nhau về khả năng phân biệt của marker.

10 Dữ liệu minh họa: glucose máu và diabetes

Ví dụ này sử dụng dữ liệu glucose máu để minh họa cách ước lượng Bayesian cho diện tích dưới đường cong ROC, tức AUC.

Biến quan sát chính là y_i, biểu diễn giá trị glucose máu của người thứ i. Mỗi người thuộc một trong hai nhóm:

• d_i = 1: người có diabetes
• d_i = 0: người không có diabetes

Trong dữ liệu này, vector y chứa các giá trị glucose máu, còn vector d chứa chỉ báo nhóm bệnh tương ứng.

Về mặt mô hình, ta giả định giá trị glucose trong mỗi nhóm có phân phối chuẩn. Cụ thể,

y_i \sim N(\mu_i, \sigma^2_{d_i}).

Mean của từng người được viết dưới dạng mô hình hồi quy tuyến tính đơn giản:

\mu_i = \beta_1 + \beta_2 d_i.

Do đó, nếu người đó không có diabetes, tức d_i = 0, thì

\mu_i = \beta_1.

Vì vậy, \beta_1 là mean glucose của nhóm không diabetes.

Nếu người đó có diabetes, tức d_i = 1, thì

\mu_i = \beta_1 + \beta_2.

Vì vậy, \beta_2 là độ chênh lệch mean glucose giữa nhóm diabetes và nhóm không diabetes.

Mô hình cũng cho phép hai nhóm có phương sai khác nhau:

\sigma_0^2

là phương sai glucose trong nhóm không diabetes, và

\sigma_1^2

là phương sai glucose trong nhóm diabetes.

Dưới giả định binormal, AUC có thể được tính từ các tham số của hai phân phối chuẩn. Nếu nhóm diabetes có mean lớn hơn nhóm không diabetes, thì

\mathrm{AUC} = \Phi\left( \frac{\beta_2}{\sqrt{\sigma_0^2+\sigma_1^2}} \right),

trong đó \Phi là hàm phân phối tích lũy của chuẩn tắc.

Công thức này có thể hiểu như sau: AUC là xác suất rằng một người có diabetes được chọn ngẫu nhiên có giá trị glucose cao hơn một người không có diabetes được chọn ngẫu nhiên.

\mathrm{AUC} = P(X_1 > X_0),

với X_1 là glucose của người có diabetes và X_0 là glucose của người không có diabetes.

11 Stan code

pacman::p_load(cmdstanr, posterior, HDInterval, gt, tidyverse)

string_stan<- "
data {
  int<lower=1> N;
  vector[N] y;
  array[N] int<lower=0, upper=1> d;
}

parameters {
  vector[2] beta;
  vector<lower=0>[2] sigma;
}

transformed parameters {

 real auc;
 real a;
 real b;
 
 auc = Phi(beta[2] / sqrt(square(sigma[1]) + square(sigma[2])));
 a = beta[2]/sigma[2];
 b = sigma[1]/sigma[2];
}

model {
  // Priors
  beta ~ normal(0, 100);

  // Weakly informative priors for group-specific standard deviations
  sigma ~ normal(0, 10);

  // Likelihood
  for (i in 1:N) {
    y[i] ~ normal(beta[1] + beta[2] * d[i], sigma[d[i] + 1]);
  }
}
"

y <- c(
  123,129,115,131,119,111,129,127,118,111,
  131,118,126,130,122,112,122,128,123,119,
  132,118,126,136,118,122,119,117,129,120,
  125,115,131,123,130,113,128,138,119,118,
  124,127,139,120,122,120,114,114,122,127,
  123,118,131,130,139,125,135,121,124,
  109,106,100,88,106,108,110,111,112,94,
  122,110,113,106,114,101,99,128,106
)

d <- c(
  rep(1, 59),
  rep(0, 19)
)



stan_data <- list(
  N = length(y),
  y = y,
  d = d
)

mod <- cmdstan_model(write_stan_file(string_stan))

fit <- mod$sample(
  data = stan_data,
  chains = 4,
  parallel_chains = 4,
  iter_warmup = 4000,
  iter_sampling = 1000,
  seed = 123,
  show_messages = F
)

# ----
fit$summary(c("beta", "sigma", "auc"), mean, hdi)|>
 mutate(variable = c("$\\beta_0$", "$\\beta_1$", "$\\sigma_0$", "$\\sigma_1$", "AUC"))|>
 mutate(across(mean:upper, ~round(.,3)))|>
 gt::gt()|>
 gt::fmt_markdown(columns = everything())|>
 gt::tab_spanner(label = "Highest density interval", column = c("lower", "upper"))

variable	mean	Highest density interval
		lower	upper
\beta_0	107.457	103.389	111.874
\beta_1	16.181	11.626	20.774
\sigma_0	9.419	6.625	12.611
\sigma_1	7.112	5.921	8.474
AUC	0.91	0.832	0.975

tt = seq(0.001, 0.99, length = 100)
d<-fit$draws(c("a", "b"), format = "matrix")
z = apply(d,1, \(x){
 a = x[1]; b = x[2]
 1 - pnorm(b*qnorm(1-tt) - a)
})|> 
 apply(1, \(k){
  c(mean = mean(k), hdi(k))
 })
 

plot(tt, z["lower",], type = "l", xlim = c(0.001, 0.99), xlab = "t", ylab = "ROC(t)")
lines(tt, z["upper",])
lines(tt, z[ "mean",], col = "red")

Đường cong ROC