mPlot: Vẽ Biểu Đồ Việc Chọn Lựa Mô Hình Cho Hồi Qui Tuyến Tính

0.1 Hồi qui tuyến tính

Trong trường hợp chúng ta muốn tìm ra một mô hình tuyến tính giữa biến số y và các biến số mô tả x(i) khác. Làm thế nào để chọn lựa được một mô hình tối ưu, mô tả được chính xác nhất sự biến thiên của y thông qua các biến số x(i).

Trước hết là chọn lựa ra được những biến số x(i) thật sự có liên quan với biến số y. Đó là trả lời câu hỏi: số lượng biến số và những biến số x(i) nào?

Tính ổn định của mô hình (model stability) là điều mà những nhà nghiên cứu quan tâm. Cụ thể, một mô hình tốt với dataset cụ thể, nhưng khi data thay đổi, dù nhỏ, đã làm cho mô hình thay đổi nhiều, tức là residual sum of squared tăng lên nhiều thì mô hình đó không ổn định. Người ta dùng thủ thuật boostrap để chọn lựa ra những mô hình ổn định nhất. Khi thử nghiệm nhiều mô hình, biến số nào có xác suất xuất hiện trong mô hình cao hơn sẽ được chọn lựa.

Có rất nhiều phương pháp được vận dụng để chọn lựa mô hình hồi qui, ở đây tôi muốn giới thiệu với các bạn chọn lựa mô hình với package mplot: Biểu đồ hóa việc chọn lựa mô hình để xem xét và quyết định chọn mô hình tối ưu.

0.2 Chuẩn bị dữ liệu

library(mplot)
data("diabetes")
head(diabetes)

##   age sex  bmi map  tc   ldl hdl tch    ltg glu   y
## 1  59   2 32.1 101 157  93.2  38   4 4.8598  87 151
## 2  48   1 21.6  87 183 103.2  70   3 3.8918  69  75
## 3  72   2 30.5  93 156  93.6  41   4 4.6728  85 141
## 4  24   1 25.3  84 198 131.4  40   5 4.8903  89 206
## 5  50   1 23.0 101 192 125.4  52   4 4.2905  80 135
## 6  23   1 22.6  89 139  64.8  61   2 4.1897  68  97

Bộ dữ liệu này gồm 10 biến số mô tả như age, sex, bmi, map và các số đo huyết tương như tc, ldl, hdl, tch, ltg, glu cùng biến số kết cục, y, đo lường sự tiến triển của bệnh tiểu đường. Chúng ta sẽ tìm mô hình hồi qui giữa y và các biến số còn lại.

0.3 Đi tìm mô hình: biến số nào có tương quan với y ?

Chúng ta khảo sát mối tương quan giữa y và các biến số còn lại.

library(corrplot)
M<-cor(diabetes)
# function for creating correlation matrix
cor.mtest <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat<- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
    }
  }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}
# matrix of the p-value of the correlation
p.mat <- cor.mtest(mtcars)
# Creating correlation plot
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(M, method="color", col=col(200),  
         type="upper", order="hclust", 
         addCoef.col = "black", # Add coefficient of correlation
         tl.col="black", tl.srt=45, #Text label color and rotation
         # Combine with significance
         p.mat = p.mat, sig.level = 0.01, insig = "blank", 
         # hide correlation coefficient on the principal diagonal
         diag=FALSE 
)

Trước hết, nhìn vào biến số y ở trung tâm của biểu đồ này. Chúng ta thấy y có tương quan cao hơn với các biến số bmi (0.59), ltg (0.57), map (0.44), tch (0.43), hdl (-0.39), glu(0.38)… Tuy rằng, những mối tương quan này không mạnh lắm.

Nếu một mô hình được chọn lựa thì các biến số trên, nhiều khả năng sẽ được đưa vào, sau khi đã loại trừ các mối tương quan lẫn nhau giữa chúng.

0.4 Package mplot

Các hàm của mplot:

vis(): tạo ra một object gồm các biến số được xem xét (variable inclusion)

af(): tạo ra một hàng rào giới hạn ngưỡng ý nghĩa của các biến số (adaptive fence)

0.4.1 Dùng Boostraped probability để so sánh các biến số

Chúng ta vẽ biểu đồ với Penalty (sum of squared coefficients) để kiểm soát các hệ số của các biến số được đưa vào trong mô hình và Boostraped probability là xác suất các biến số được chọn lựa trong những mô hình đang được máy tính đánh giá.

lm.d <- lm(y ~ ., data = diabetes)
vis.d <- vis(lm.d, B = 200, seed = 2018)
plot(vis.d, interactive = FALSE, which = "vip")

Đây là kết quả sau B = 200 lần boostraping data. Biểu đồ này lấy đường RV (redundent variable) làm ngưỡng ý nghĩa. Hai biến số glu và age không được xem là có ý nghĩa khi được đưa vào mô hình vì chúng nắm dưới đường RV.

bmi là biến số xuất hiện nhiều nhất trong các mô hình với xác suất 100%. Kế tiếp là các biến số ltg, map, sex…

0.4.2 Dùng -2*Log-LikeliHood (residual sum of squared: -2LL) để so sánh các biến số

plot(vis.d, interactive = FALSE, which = "lvk")

Một cách khái quát khi số lượng biến số được đưa vào mô hình nhiều hơn thì -2LL giảm đi, nói cách khác các mô hình chính xác hơn. Chương trình tự động đưa biến số age vào để so sánh.

Tuy nhiên, biểu đồ cho thấy không có sự khác biệt đáng kể về -2LL trong những mô hình có biến số age và không age. Khi chọn từ 2-3 biến số đẩ đưa vào mô hình thì age càng không thể xuất hiện.

Tương tự, chúng ta có thể so sánh vai trò của 2 biến số bmi và hdl trong các mô hình qua biểu đồ sau.

plot1 <- plot(vis.d, interactive = FALSE, highlight = "bmi", which = "lvk")
plot2 <- plot(vis.d, interactive = FALSE, highlight = "hdl", which = "lvk")
library(gridExtra)
grid.arrange(plot1, plot2, ncol = 2)

Vai trò của hdl không khác mấy với biến số age.

Ngược lại, bmi có vai trò quan trọng trong các mô hình. Biểu đồ bên trái, khi có bmi trong mô hình (màu đỏ) thì -2LL giảm đáng kể so với những mô hình không có bmi (màu xanh).

0.4.3 Kết hợp -2LL và Boostraped probability để so sánh các mô hình trên biểu đồ

plot1.1 <- plot(vis.d, interactive = FALSE, which = "boot", max.circle = 10, highlight = "bmi") 
plot1.2 <- plot(vis.d, interactive = FALSE, which = "boot", max.circle = 10, highlight = "hdl") 
grid.arrange(plot1.1, plot1.2, ncol = 2)

Kết quả được thể hiện tương tự trong sự kết hợp hai thông số trên.

Những mô hình có sự hiện diện của bmi cho -2LL thấp hơn đáng kể trong sự so sanh. Vòng tròn to cho biết mô hình đó có xác suất chọn lựa cao (gần bằng 1). Những mô hình có từ 7 biến số trở lên không cải thiện hơn -2LL. Mô hình có 5 - 6 biến số đã đạt được -2LL thấp tương đối, với xác suất chọn lựa cao (vòng tròn) khá lớn.

0.4.4 Xem bảng xác suất chọn mô hình và LogLikelihood

print(vis.d)

##                                     name prob logLikelihood
##                                      y~1 1.00      -2547.17
##                                    y~bmi 0.66      -2454.02
##                                    y~ltg 0.34      -2461.86
##                                y~bmi+ltg 0.99      -2411.20
##                            y~bmi+map+ltg 0.66      -2402.61
##                         y~bmi+map+tc+ltg 0.40      -2397.48
##                        y~bmi+map+hdl+ltg 0.34      -2397.71
##                    y~sex+bmi+map+hdl+ltg 0.69      -2390.13
##                 y~sex+bmi+map+tc+ldl+ltg 0.40      -2387.30
##  y~sex+bmi+map+tc+ldl+hdl+tch+ltg+glu+RV 0.30      -2385.39

Xác suất chọn lựa càng cao càng tốt đồng thời loglikelihood càng thấp càng tốt. Vậy mô hình:

y ~ sex + bmi + map + hdl + ltg

Với P = 69% và LogLikelihood = - 2390.13 có thể là tối ưu để chọn lựa.

0.4.5 Sử dụng hàm af

Trong những mô hình được so sánh, mô hình nào có xác suất chọn lựa cao nhất với số lượng biến số ít nhất sẽ là tối ưu. Vì vậy hãy tìm đỉnh xác suất đầu tiên để phát hiện mô hình tốt nhất.

af.d <- af(lm.d, B = 200, n.c = 100, c.max = 100, seed = 2018)
plot(af.d, interactive = FALSE, best.only = TRUE, legend.position = "right")

Đỉnh xác suất đầu tiên đạt được với những mô hình gồm 5 biến số kể trên trong biểu đồ adaptive fence. Chúng ta có thêm bằng chứng để chọn lựa mô hình này.

y ~ sex + bmi + map + hdl + ltg

Đi đến kết luận: Tiến triển của bệnh tiểu đường có thể tiên lượng thông qua các biến số như giới tính của bệnh nhân, cộng với bmi, áp suất động mạch (map), hdl và ltg.

0.4.6 Interactive plots

Sẽ rất thú vị với những biểu đồ interactive, trong đó các bạn có thể chọn lựa các thông số để biểu đồ thay đổi tức thời các kết quả với code:

mplot(lm.d, vis.d, af.d)

0.5 Kết luận

Với sự hỗ trợ của package mplot, chúng ta có thêm một công cụ để phân tích và chọn lựa mô hình linear regression.

Lưu ý mplot cũng hoạt động tốt với các mô hình glm (hồi qui logistic - classification).

0.6 Đọc thêm về mplot package

https://www.jstatsoft.org/article/view/v083i09