Bài toán hồi quy

Keras v.s Gamlss

1 Giới thiệu

Trong bài thực hành ngắn này, Nhi muốn làm một thí nghiệm nhỏ để so sánh khả năng của 2 giải thuật khác nhau khi áp dụng cho cùng một bài toán hồi quy: một bên là Mô hình Deep neural network bằng Tensorflow (giao thức Keras), bên kia là mô hình Lambda-Mu-Sigma (LMS) thuộc phương pháp GAMLSS. Mỗi loại có thể xem như ứng cử viên mạnh nhất từ trường phái Whiteboxes và Blackboxes.

Tình huống đặt ra, đó là một nghiên cứu dịch tễ nhằm ước tính giá trị dung phổi bình thường (TLC) ở nam giới thuộc chủng tộc Caucasian.

2 Dữ liệu

library(tidyverse)

df=read.csv("https://raw.githubusercontent.com/kinokoberuji/R-Tutorials/master/LMSmodelGAMLSS.csv", sep=";")%>%as_tibble()

dfM=filter(df,Sex=="M")%>%dplyr::select(-Sex)

head(dfM)%>%knitr::kable()

TLC	Age	Height
7.017265	18	177.5
7.648100	18	181.0
6.737513	18	176.0
5.949827	22	171.0
7.389043	22	180.0
7.048908	23	185.0

3 Mô hình LMS (Gamlss)

Trước hết, ta sẽ dùng mô hình LMS từ package gamlss. Lý thuyết về mô hình LMS đã được Nhi giải thích cặn kẽ trong bài thực hành số 7 series Gamlss : http://rpubs.com/lengockhanhi/342110. Một cách ngắn gọn, LMS là một tập hợp 3 mô hình riêng biệt M,S,L cho 3 tham số : Mu (vị trị trung tâm), Sigma (Scale) và Lambda (Skewness), tương ứng với 3 tham số Mu (link function = log), Sigma (link=log) và Nu (link = identity) của quy luật phân phối Box-Cox-Cole-Green. Mỗi mô hình M,L,S là một mô hình GAM (generalized additive model), cho phép tích hợp thêm Smoothing splines (các hàm điều hòa/hiệu chỉnh) và hàm đa thức (Polynomial). Chính những hàm này cho phép đồ thị LMS model uốn lượn mềm mại theo những hình dạng phức tạp nhất có thể và tăng độ chính xác của mô hình.

Đầu tiên, ta sẽ dùng caret để chia dữ liệu thành 2 phần, 20% dành cho kiểm định và 80% để dựng mô hình:

library(caret)

set.seed(123)

idM=createDataPartition(y=dfM$TLC, p=0.8,list=FALSE)
trainM=dfM[idM,]
testM=dfM[-idM,]

Ta dựng mô hình LMS bằng quy trình Stepwise, trong đó ta sẽ xác định tổ hợp tối ưu cho 2 câu hỏi: Có cần mô hình cho Sigma/Nu hay không ? và: Nếu cần, thì nội dung mô hình là gì ?

Trong trường hợp này Nhi muốn chọn lọc mô hình L,M,S tối ưu dựa vào 3 giả định: hàm đa thức bậc 3 cho Height , bậc 2 cho Age, kèm hoặc không kèm Spline cho Age, công thức này sẽ lần lượt test cho Mu và Sigma, Nu được mặc định là hằng số.

library(gamlss)
nC<-detectCores()

m1=gamlss(data=trainM,
          TLC~1,
          sigma.formula = TLC~1,
          nu.formula = TLC~1,
          family=BCCG(mu.link="log"),
          trace=FALSE,
          parallel="multicore",
          ncpus = nC)

m2=stepGAICAll.A(m1, scope=list(lower=~1, 
                                upper=~poly(Age,2)+poly(Height,3)+pb(Age)),
                 sigma.scope = list(lower=~1, 
                                upper=~poly(Age,2)+poly(Height,3)+pb(Age)),
                 k=log(length(trainM)),
                 trace=FALSE,
                 parallel="multicore",
                 ncpus = nC
              )

## --------------------------------------------------- 
## Start:  AIC= 546.27 
##  TLC ~ 1 
## 
## --------------------------------------------------- 
## Start:  AIC= 437.74 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 435.34 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 435.34 
##  ~pb(Age) 
## 
## --------------------------------------------------- 
## Start:  AIC= 435.34 
##  TLC ~ poly(Height, 3) + pb(Age) 
## 
## ---------------------------------------------------

Kết quả Stepwise như sau:

Mô hình cho Mu có công thức: Mu ~ poly(Height,3)+Age+ MuSpline.

Mô hình cho Sigma chỉ gồm Sigma ~ Age + SigmSplin.

Nu = hằng số.

summary(m2)

## ******************************************************************
## Family:  c("BCCG", "Box-Cox-Cole-Green") 
## 
## Call:  
## gamlss(formula = TLC ~ poly(Height, 3) + pb(Age), sigma.formula = ~pb(Age),  
##     nu.formula = ~1, family = BCCG(mu.link = "log"),  
##     data = trainM, trace = FALSE, parallel = "multicore",      ncpus = nC) 
## 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.9041367  0.0222386  85.623  < 2e-16 ***
## poly(Height, 3)1  1.1785602  0.1076800  10.945  < 2e-16 ***
## poly(Height, 3)2 -0.3753507  0.1063796  -3.528 0.000527 ***
## poly(Height, 3)3 -0.1604175  0.1038782  -1.544 0.124224    
## pb(Age)           0.0009752  0.0005334   1.828 0.069093 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.568038   0.155515  -16.51   <2e-16 ***
## pb(Age)      0.006680   0.003356    1.99    0.048 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.05722    0.59892   0.096    0.924
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  195 
## Degrees of Freedom for the fit:  9.53894
##       Residual Deg. of Freedom:  185.4611 
##                       at cycle:  5 
##  
## Global Deviance:     424.8601 
##             AIC:     443.938 
##             SBC:     475.1589 
## ******************************************************************

4 Mô hình Deep neural network với Keras

Keras là một nền tảng chuyên dụng cho Deep learning. Bản thể của Keras dùng ngôn ngữ Python, nhưng có giao thức để sử dụng trong R (package keras), bạn cần install trước một gói Python thí dụ Anaconda để có thể sử dụng keras trong R.

Giao thức Keras tiếp nhận một cấu trúc dữ liệu đặc biệt gọi là các tensor. Ta có thể hình dung về tensors như khái niệm tổng quát để gọi tên những đối tượng dữ liệu mà ta từng biết trong R theo số chiều thông tin (dimension). Trong thí dụ này, ta có thể gọi vector biến kết quả là là 1D tensor, còn matrix hay dataframe gồm các features được gọi là 2D tensor.

Trong R, tensor được lưu dưới dạng array, do đó trước hết ta cần hoán chuyển những dataframe hiện có thành array. Điểm khác biệt thứ hai đó là Keras tiếp nhận riêng tensors cho tập features và cho biến kết quả nên ta phải cắt 2 phần này riêng lẻ từ tập train và test

library(keras)

train_data<-trainM%>%dplyr::select(-1)%>%as.matrix()%>%as.array.default(dimnames=NULL)  # Convert dataframe to array
dimnames(train_data) <- NULL

train_targets<-trainM%>%.$TLC%>%as.array() 

test_data<-testM%>%dplyr::select(-1)%>%as.matrix()%>%as.array.default(dimnames=NULL)  # Convert dataframe to array
dimnames(train_data) <- NULL

test_targets<-testM%>%.$TLC%>%as.array() 

str(train_data)

##  num [1:195, 1:2] 18 18 18 22 23 23 24 25 25 26 ...

str(test_data)

##  num [1:48, 1:2] 22 27 32 33 34 34 35 47 51 61 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:2] "Age" "Height"

str(train_targets)

##  num [1:195(1d)] 7.02 7.65 6.74 5.95 7.05 ...

str(test_targets)

##  num [1:48(1d)] 7.39 7.78 8.69 8.04 6.75 ...

Một điểm cần lưu ý tiếp theo đó là Neural network chỉ hoạt động tối ưu khi dữ liệu đầu vào được chuẩn hóa bằng hàm scale, sử dụng trung bình của tập train:

mu <- apply(train_data, 2, mean)
sigma <- apply(train_data, 2, sd)
train_data <- scale(train_data, center = mu, scale = sigma)

test_data <- scale(test_data, center = mu, scale = sigma)

head(train_data)%>%knitr::kable()

-1.616520	0.0889939
-1.616520	0.5362570
-1.616520	-0.1026903
-1.365097	-0.7416377
-1.302241	1.0474149
-1.302241	-1.7000587

Trong thí dụ này, do bài toán rất đơn giản nên ta chỉ cần dùng một mạng neuron có cấu trúc 2 lớp hidden, mỗi lớp gồm 64 đơn vị neurons. Đầu ra của mạng lưới là 1 neuron duy nhất (và không kèm activation function như trong bài toán tiên lượng xác suất) – đủ để xuất kết quả là một con số (scalar) và tuyến tính – chính là giá trị cần ước lượng.

Mạng lưới được đóng gói với loss function là mean squared error (sai số bình phương trung bình) vì đây là bài toán hồi quy. Trong quá trình huấn luyện, tiêu chí được sử dụng để cải thiện hiệu năng của mô hình là MSE (sai số tuyệt đối trung bình).

build_model <- function() {
  model <- keras_model_sequential() %>% 
    layer_dense(units = 64, activation = "relu", 
                input_shape = dim(train_data)[[2]]) %>% 
    layer_dense(units = 64, activation = "relu") %>% 
    layer_dense(units = 1) 
  
  model %>% compile(
    optimizer = "rmsprop", 
    loss = "mse", 
    metrics = c("mae")
  )
}

Tiếp theo, ta thực hiện một quy trình kiểm chứng chéo 5 blocks lặp lại 100 lần

set.seed(123)

k <- 5
indices <- sample(1:nrow(train_data))
folds <- cut(1:length(indices), breaks = k, labels = FALSE) 

num_epochs <- 100

all_mae_histories <- NULL
for (i in 1:k) {
  cat("processing fold #", i, "\n")
  
  # Prepare the validation data: data from partition # k
  val_indices <- which(folds == i, arr.ind = TRUE)
  val_data <- train_data[val_indices,]
  val_targets <- train_targets[val_indices]
  
  # Prepare the training data: data from all other partitions
  partial_train_data <- train_data[-val_indices,]
  partial_train_targets <- train_targets[-val_indices]
  
  # Build the Keras model (already compiled)
  model <- build_model()
  
  # Train the model (in silent mode, verbose=0)
  history <- model %>% fit(
    partial_train_data, partial_train_targets,
    validation_data = list(val_data, val_targets),
    epochs = num_epochs, batch_size = 1, verbose = 0
  )
  mae_history <- history$metrics$val_mean_absolute_error
  all_mae_histories <- rbind(all_mae_histories, mae_history)
}

## processing fold # 1 
## processing fold # 2 
## processing fold # 3 
## processing fold # 4 
## processing fold # 5

median_mae_history <- data.frame(
  epoch = seq(1:ncol(all_mae_histories)),
  validation_mae = apply(all_mae_histories, 2, median)
)

median_mae_history%>%ggplot(aes(x = epoch, y = validation_mae))+geom_line(col="red",size=1)+
  theme_bw()+
  geom_hline(yintercept = median(median_mae_history$validation_mae),linetype=2)

ggplot(median_mae_history, aes(x = epoch, y = validation_mae))+geom_smooth(col="red",fill="red")+
  theme_bw()+
  geom_hline(yintercept = median(median_mae_history$validation_mae),linetype=2)

Cuối cùng, ta dựng một mô hình chính thức với số epochs = 50 và batch_size=16

set.seed(123)

model <- build_model()

# Train it on the entirety of the data.
model %>% fit(train_data, train_targets,
              epochs = 50, batch_size = 16, verbose = 0)

summary(model)

## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #     
## ===========================================================================
## dense_16 (Dense)                 (None, 64)                    192         
## ___________________________________________________________________________
## dense_17 (Dense)                 (None, 64)                    4160        
## ___________________________________________________________________________
## dense_18 (Dense)                 (None, 1)                     65          
## ===========================================================================
## Total params: 4,417
## Trainable params: 4,417
## Non-trainable params: 0
## ___________________________________________________________________________

result <- model %>% evaluate(test_data, test_targets)

result

## $loss
## [1] 0.6035802
## 
## $mean_absolute_error
## [1] 0.6450254

5 So sánh hiệu năng 2 mô hình

Truth=testM$TLC
DNN=predict(model,test_data)
LMS=predict(m2,newdata=testM,type="response")

## new prediction 
## New way of prediction in pb()  (starting from GAMLSS version 5.0-3)

pdf=cbind(Truth,DNN,LMS)%>%as_data_frame()

colnames(pdf)=c("Truth","DeepNN","LMS")

pdf$Age=testM$Age

library(mlr)

regr.task= mlr::makeRegrTask(id = "dfM", data=testM, target = "TLC")
regr.lrn = makeLearner("regr.glm")

dummy=mlr::train(regr.lrn,regr.task)

dumpredDNN=predict(dummy,regr.task)
dumpredLMS=predict(dummy,regr.task)

dumpredDNN$data$response<-DNN
dumpredLMS$data$response<-LMS

Đây là hiệu năng mô hình Deep neural net với Tensorflo (Keras)

mets=list(mse,mae,medae,rmse)
performance(dumpredDNN,measures =mets)

##       mse       mae     medae      rmse 
## 0.6035802 0.6450254 0.5570394 0.7769043

Còn đây là hiệu năng mô hình LMS của gamlss

performance(dumpredLMS,measures =mets)

##       mse       mae     medae      rmse 
## 0.5183984 0.5819141 0.5102973 0.7199989

Màu xanh là mô hình LMS, màu đỏ là mô hình Deep neural network

pdf%>%ggplot()+geom_point(aes(x=Age,y=Truth),col="black",alpha=0.5)+
  geom_smooth(aes(x=Age,y=DNN),col="red",fill="red",alpha=0.2)+
  geom_smooth(aes(x=Age,y=LMS),col="blue",fill="blue",alpha=0.2)+theme_bw()+scale_y_continuous("TLC")

## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'

Trung bình, mô hình Deep neural network chỉ tiên lượng dung tích phổi sai khoảng 50 mL

mean(pdf$DeepNN-pdf$Truth)

## [1] -0.02100793

6 Nhận xét

Deep neural network (còn gọi là Deep learning) là một giải pháp vô cùng mạnh mẽ cho bài toán hồi quy, thậm chí nó tương đương với mô hình phi tuyến tính với hàm Spline bù trừ vốn là công cụ tối ưu nhất hiện nay cho các mô hình tăng trưởng trong y học lâm sàng.

Tuy nhiên, mô hình LMS có 2 ưu điểm mà Deep NN không thể thay thế, thứ nhất là tính tường minh: nội dung mô hình có thể được tính thủ công từ lookup table, thứ hai là khả năng ước tính được các ngưỡng giới hạn trên và dưới của giá trị bình thường trong quần thể từ 3 tham số L,M và S, cho phép diễn giải kết quả về mặt lâm sàng.Trong khi đó mô hình Deep NN không thể diễn giải được và chỉ có thể được áp dụng dưới dạng hard code trong firmware của thiết bị xét nghiệm.

Tuy nhiên, có thể dự báo trong tương lai giải pháp Deep learning sẽ cho phép những thiết bị xét nghiệm có khả năng tự tạo ra mô hình ước lượng giá trị tham chiếu cho chính nó mà không cần sự can thiệp của con người, và đây là một ứng dụng không thể xem thường.