Tensorflow Deep learning

Ứng dụng Hồi quy

1 Dẫn nhập: Mô hình giá trị tham chiếu

Trong thực hành lâm sàng, kết quả các xét nghiệm, thí dụ chức năng hô hấp (Pulmonary function test, PFT) chỉ là những con số vô nghĩa nếu không có giá trị tham chiếu (Reference values). Giá trị tham chiếu thường được ước lượng từ một mô hình hồi quy dựa trên dữ liệu từ quần thể người khỏe mạnh. Mô hình này thường có nguyên tắc chung là : dựa vào một hay nhiều biến số đầu vào (predictors) dễ đo đạc, độc lập với nhau nhưng có liên hệ chặt chẽ với kết quả cần dự báo (outcome), thí dụ Chủng tộc, Giới tính, các chỉ số nhân trắc (trọng lượng, chiều cao, BMI…).Nhiệm vụ của nhà thống kê là xác định hàm f cho phép ước lượng Outcome Y theo những predictors này.

\[Outcome Y \sim f(predictors...)\]

Khi áp dụng vào chẩn đoán, kết quả xét nghiệm (Observed, true value) của một cá thể bất kì sẽ được so sánh với giá trị dự báo (Predicted value) từ mô hình (hàm f), ý nghĩa của giá trị dự báo này là trung bình của Outcome Y trong một quần thể giả định có cùng đặc tính nhân trắc học, chủng tộc và giới tính với cá thể này. Do tuổi được dùng như predictor, mô hình f này đồng thời cũng đã mô phỏng sự thay đổi của Y dưới tác động của quy luật sinh lý bình thường (tăng trưởng, già lão, thoái hóa tự nhiên), nên người bác sĩ có thể loại trừ những yếu tố này khỏi biện luận chẩn đoán để có thể tập trung vào yếu tố nghi ngờ là có sự tăng/giảm bất thường của Y do nguyên nhân bệnh lý.

2 Nhược điểm của phương pháp thủ công

Trong hàng chục năm qua, Y giới tạo ra những mô hình giá trị tham chiếu theo cách hoàn toàn thủ công. Đầu tiên, xét nghiệm được thực hiện trên một mẫu với kích thước từ vài trăm đến hàng ngàn cá thể, được giả định là đại diện cho quần thể người khỏe mạnh. Dữ liệu sau đó phải được nhập thủ công vào máy tính và tích tụ theo thời gian cho đến khi đạt cỡ mẫu mong muốn. Thậm chí khi những nghiên cứu loại này được tổ chức ở cấp quốc gia, quá trình này thường rất chậm, và có thể kéo dài từ 1-3 năm tùy theo cỡ mẫu.

Sau khi có dữ liệu, một chuyên viên thống kê sẽ bắt đầu dựng mô hình. Đồng hành với lịch sử phát triển của ngành Thống kê, độ phức tạp của hàm f dần dần tăng lên – từ một hằng số trung bình đến mô hình tuyến tính đơn biến, rồi đa biến, rồi đến mô hình đa thức, phi tuyến tính. Vào năm 2018 này, đỉnh cao nhất về sự tinh tế của hàm f là mô hình 3 tham số LMS kết hợp với hàm spline và hàm đa thức. Tuy nhiên, có một điều không thay đổi, đó là việc dựng mô hình hoàn toàn thủ công, thường bằng phương pháp Step-wise. Quy trình này có thể kéo dài vài tháng hoặc lâu hơn.

Sau khi có mô hình, các tác giả sẽ công bố nó dưới dạng 1 bài báo khoa học trên tạp chí. Quy trình bình duyệt cho bài báo có thể dài từ 1 đến 6 tháng nữa cho đến khi nó thực sự được công bố. Sau khi xuất bản, bài báo sẽ được chuyển qua ngành công nghiệp thiết bị xét nghiệm, nơi các kỹ sư phần mềm sẽ đưa mô hình f vào code của phần mềm (PC software hay firmware). Công đoạn này có thể mất thêm vài tháng nữa.

Như vậy, từ khi ý tưởng về mô hình được để xuất cho đến khi nó thực sự được áp dụng trên lâm sàng thường phải mất từ 2-5 năm.

3 Ý tưởng về thiết bị xét nghiệm thông minh

Bây giờ, bạn thử tưởng tượng về một quy trình mới cho phép một cái máy xét nghiệm có khả năng tự thu thập dữ liệu một cách liên tục và ở thời gian thực, học từ dữ liệu này ngày này qua ngày khác và tự tìm ra quy luật (hàm f) tối ưu nhất cho phép ước lượng giá trị tham chiếu rồi áp dụng cho chính nó mà không cần một sự can thiệp nào từ con người ? Hơn nữa, cái máy còn có khả năng tự khắc phục những sai sót của mô hình khi tiếp xúc với dữ liệu mới và cập nhật mô hình để việc dự báo luôn chính xác nhất có thể ?

Đây có lẽ là hình thức sơ khai nhất của « trí tuệ nhân tạo », nhưng không phải là một chuyện viễn tưởng. Trong bài thực hành sau đây, Nhi sẽ thực hiện một thí nghiệm để minh họa cho ý tưởng này. Vào cuối bài, các bạn sẽ thấy rằng ý tưởng xây dựng giá trị tham chiếu tự động không những khả thi, mà còn tốt hơn so với cách làm mô hình thủ công, thậm chí với công cụ tốt nhất hiện nay là mô hình LMS.

4 Thí nghiệm minh họa

Trong thí dụ này, Nhi đặt vấn đề giả định như sau : Một khoa xét nghiệm chức năng hô hấp tại quốc gia C sở hữu một database kết quả thể tích phế nang (AV) và hệ số khuếch tán khí CO (DLCO) trên 627 người bình thường. Chưa có giá trị tham chiếu cho AV và DLCO ở người dân của quốc gia C. Trưởng khoa giao nhiệm vụ cho phòng Thống kê bệnh viện để dựng mô hình dự báo giá trị tham chiếu AV và DLCO cho riêng Nam và Nữ.Chuyên viên thống kê khá tự tin cho biết anh ta biết làm mô hình Lamnda Mu Sigma (LMS), một phương pháp tối tân nhất cho mô hình tham chiếu hiện được Hội hô hấp Châu Âu, Hoa kỳ và WHO áp dụng.

Cùng lúc đó, một cô gái làm ở phòng quản lý dữ liệu của bệnh viện muốn thử nghiệm một ý tưởng khác : Tôi sẽ áp dụng một Deep neural network (mạng thần kinh nhân tạo sâu) với Tensorflow để tạo ra một framework nhằm thực hiện cùng mục tiêu : Ước lượng giá trị AV và DLCO dựa vào dữ liệu đầu vào là 3 biến : Tuổi, Giới tính, Chiều cao. Hãy cho tôi cơ hội làm thử, rồi chúng ta sẽ so sánh hiệu quả giữa cách của tôi và cách của anh.

Một cuộc thi đấu bắt đầu …

Trước hết, dữ liệu được chia ngẫu nhiên ra 2 phần: một phần dùng để dựng mô hình (n=500) chung cho cả 2 phương pháp LMS và Deep learning, một phần nhỏ hơn (n=127) sẽ được dùng làm dữ liệu kiểm định phẩm chất và so sánh 2 mô hình.

Lưu ý rằng dữ liệu kiểm định sẽ được dấu kín, cả 2 phía đều không biết gì về nội dung của dữ liệu này cho đến khi nó được dùng để so sánh mô hình của họ.

5 Thăm dò dữ liệu

library(tidyverse)

ndf<-read.csv("DLCOkeras.csv",sep=";")

ndf<-ndf%>%dplyr::select(-Group)%>%filter(DLCO>11.5,DLCO<55)

# original train/test data
idx=caret::createDataPartition(y=ndf$Sex, p=499/627,list=FALSE)
trainset=ndf[idx,]
testset=ndf[-idx,]

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

	Sex	Age	Height	AV	DLCO
1	F	20	173.0	6.00	29.2
2	F	22	165.0	4.76	26.3
3	F	22	168.6	5.40	27.9
4	F	23	164.0	5.36	25.4
5	F	24	170.0	5.77	29.3
8	F	25	155.0	4.61	26.6

Dữ liệu gồm 5 biến, 2 biến kết quả là AV và DLCO, 3 predictors là Sex, Age và Height. Trước hết, ta biết chắc chắn rằng AV hoặc DLCO phải liên hệ với Sex,Age và Height theo một kiểu nào đó, thí dụ: nữ giới thường có giá trị thấp hơn Nam giới, người già suy giảm giá trị do hiện tượng lão hóa sinh lý, dung tích phổi thì tỉ lệ thuận với chiều cao…

Tuy nhiên, vấn đề khó khăn đó là: những mối liên hệ này là phi tuyến tính như biểu đồ bên dưới:

# Data visualisation

trainset%>%gather(AV,DLCO,key="Outcome",value="Value")%>%
  ggplot(aes(x=Age,y=Value))+
  geom_point(aes(col=Sex),alpha=0.2)+
  geom_smooth(aes(fill=Sex,col=Sex),alpha=0.5)+
  theme_bw()+
  facet_grid(Outcome~Sex,scales = "free")+
  scale_fill_manual(values=c("red","blue"))+
  scale_color_manual(values=c("red3","blue3"))

trainset%>%gather(AV,DLCO,key="Outcome",value="Value")%>%
  ggplot(aes(x=Height,y=Value))+
  geom_point(aes(col=Sex),alpha=0.2)+
  geom_smooth(aes(fill=Sex,col=Sex),alpha=0.5)+
  theme_bw()+
  facet_grid(Outcome~Sex,scales = "free")+
  scale_fill_manual(values=c("red","blue"))+
  scale_color_manual(values=c("red3","blue3"))

Một mô hình hồi quy tốt phải có khả năng mô phỏng được quy luật phi tuyến tính này theo tuổi và/hoặc chiều cao. Ta sẽ theo dõi xem anh chuyên viên thống kê và chị chuyên viên data science xử lý vấn đề như thế nào nhé :

6 Mô hình LMS thủ công bằng Step-wise

Theo truyền thống, chúng ta cần phải dựng một mô hình riêng cho mỗi giới tính (M, F), do đó quy trình bắt đầu bằng việc chia dữ liệu gốc thành 2 phần cho riêng Nam và Nữ:

# trainset for LMS
trainLMS_M<-trainset%>%filter(Sex=="M")
trainLMS_F<-trainset%>%filter(Sex!="M")

# testset for LMS
testLMS_M<-testset%>%filter(Sex=="M")
testLMS_F<-testset%>%filter(Sex!="M")

Phương pháp mô hình 3 tham số (Lambda-Mu-Sigma, LMS) do 2 tác giả Rigby RA. và Stasinopoulos DM tạo ra vào năm 2001 để ước lượng một biến số ngẫu nhiên, liên tục và không âm, có quy luật biến thiên phi tuyến tính theo thời gian (Tuổi).LMS linh hoạt vì không dựa vào phân phối Gaussian như đa số mô hình thống kê khác. LMS cho phép tái hiện quỹ tích phức tạp của Y.

LMS có bản chất là một dạng đặc biệt thuộc họ mô hình GAM (Generalised additive model).Nền tảng của LMS model là họ phân phối Box-Cox-Cole-Green (BCCG), được xác định bởi 3 tham số, bao gồm vị trí trung tâm hoặc trung bình (Mu), Scale (Sigma) và Skewness (Nu).

LMS là một distributional model, tức 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ột mô hình hồi quy GAMLSS với biến kết quả Outcome là Y tuân theo quy luật BCCG :

\[Y \sim BCCG(\mu ,\sigma ,\nu )\]

Mỗi tham số theta (thí dụ Mu, Sigma, Nu, Zero_inflated,…) có thể được xem như một outcome thứ cấp, và có thể được mô hình hóa riêng biệt bởi một link function f với nội dung mô hình lần lượt như sau :

\[\mu \sim ln(\beta X + BSpline(x))\]

\[\sigma \sim ln(\beta X + BSpline(x))\]

\[\nu \sim identity(\beta X + BSpline(x))\]

Với : Mu, sigma và Nu lần lượt là 3 tham số của phân phối BCCG, Beta đại diện cho tham số hồi quy tương ứng với một model matrix X, là tập hợp của nhiều biến predictor trong mô hình, thí dụ Tuổi, Chiều cao trong mô hình dự báo dung AV hoặc DLCO…

Từ mô hình LMS, ta có thể ước tính : giá trị trung bình lý thuyết của Outcome, 2 ngưỡng cao nhất và thấp nhất.

\[LLN = Exp \left [ ln(\mu)+\frac{ln(1-1.645*\nu \sigma)}{\nu}\right ] \]

\[ULN = Exp \left [ ln(\mu)+\frac{ln(1+1.645*\nu \sigma)}{\nu}\right ]\] Mỗi mô hình M,L,S là một mô hình GAM (generalized additive model), nên matrix X 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 polynomial và smoothing 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 tính chính xác của mô hình. Tuy nhiên đây cũng có thể trở thành nhược điểm vì nếu quá tham lam, ta sẽ tạo ra mô hình bị overfit, nhất là khi dữ liệu trainset hỗn tạp và có số lượng lớn. Mô hình overfit chỉ hoạt động tốt trên chính dữ liệu dùng để tạo ra nó, nhưng sẽ phạm sai lầm khi áp dụng trên thực tế với dữ liệu mới.

Mô hình LMS được dựng bằng package gamlss trong R. Ta sẽ sử dụng một chiến thuật đơn giản đó là Backward Step-wise (triệt thoái tự động, dựa vào BIC). Thí dụ ta có đến 3 mô hình con là L,M,S; và mỗi mô hình dùng 1 đa thức bậc cao bao gồm smoothing spline, và tất cả tổ hợp có thể giữa chúng phải được kiểm tra bằng câu hỏi: Có cần mô hình bậc 2,3,4 cho Age hoặc Height không cho Mu/Sigma/Nu ?

Các biến trong formula sẽ bị triệt tiêu dần đến khi mô hình có AIC tối ưu hiện ra.

Sau đây lần lượt là 4 mô hình cho Nam/nữ, AV và DLCO:

# LMS

library(gamlss)
nC<-detectCores()

# DLCO for Male

mm1=gamlss(data=trainLMS_M[,-1],
          DLCO~1,
          sigma.formula = DLCO~1,
          nu.formula = DLCO~1,
          family=BCCG(mu.link="log"),
          trace=FALSE,
          parallel="multicore",
          ncpus = nC)

mm2=stepGAICAll.A(mm1,scope=list(lower=~1, 
                                upper=~poly(Age,3)+poly(Height,3)+pb(Age)),
                 sigma.scope = list(lower=~1, 
                                    upper=~poly(Age,3)+poly(Height,3)+pb(Age)),
                 nu.scope = list(lower=~1,upper=~poly(Age,3)+poly(Height,3)),
                 k=log(length(trainLMS_M[,-1])),
                 trace=FALSE,
                 parallel="multicore",
                 method=RS(1000),
                 ncpus = nC
)

## --------------------------------------------------- 
## Start:  AIC= 1665.44 
##  DLCO ~ 1 
## 
## --------------------------------------------------- 
## Start:  AIC= 1548.02 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 1546.01 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 1536.22 
##  ~poly(Height, 3) 
## 
## --------------------------------------------------- 
## Start:  AIC= 1536.22 
##  DLCO ~ poly(Age, 3) + poly(Height, 3) 
## 
## ---------------------------------------------------

modDL_M=mm2

summary(modDL_M)

## ******************************************************************
## Family:  c("BCCG", "Box-Cox-Cole-Green") 
## 
## Call:  
## gamlss(formula = DLCO ~ poly(Age, 3) + poly(Height, 3), sigma.formula = ~poly(Height,  
##     3), nu.formula = ~poly(Height, 3) + poly(Age, 3), family = BCCG(mu.link = "log"),  
##     data = trainLMS_M[, -1], trace = FALSE, parallel = "multicore",  
##     ncpus = nC) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.533025   0.008146 433.696  < 2e-16 ***
## poly(Age, 3)1    -1.533083   0.121116 -12.658  < 2e-16 ***
## poly(Age, 3)2    -0.542061   0.112982  -4.798 2.78e-06 ***
## poly(Age, 3)3    -0.071340   0.101436  -0.703    0.483    
## poly(Height, 3)1  0.838180   0.124046   6.757 1.01e-10 ***
## poly(Height, 3)2 -0.098194   0.082379  -1.192    0.234    
## poly(Height, 3)3 -0.054692   0.078729  -0.695    0.488    
## ---
## 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)      -1.89192    0.04446 -42.558  < 2e-16 ***
## poly(Height, 3)1 -0.97254    0.70711  -1.375   0.1703    
## poly(Height, 3)2  1.83223    0.70711   2.591   0.0101 *  
## poly(Height, 3)3  2.87401    0.70711   4.064 6.46e-05 ***
## ---
## 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.7268     0.3211   2.264 0.024451 *  
## poly(Height, 3)1  17.1680     4.7256   3.633 0.000341 ***
## poly(Height, 3)2  44.6022     4.4419  10.041  < 2e-16 ***
## poly(Height, 3)3  13.1565     3.5889   3.666 0.000302 ***
## poly(Age, 3)1     21.4322     4.7096   4.551 8.41e-06 ***
## poly(Age, 3)2    -16.6489     4.7888  -3.477 0.000600 ***
## poly(Age, 3)3     19.7044     4.6986   4.194 3.83e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  253 
## Degrees of Freedom for the fit:  18
##       Residual Deg. of Freedom:  235 
##                       at cycle:  20 
##  
## Global Deviance:     1511.266 
##             AIC:     1547.266 
##             SBC:     1610.867 
## ******************************************************************

# DLCO for female

mf1=gamlss(data=trainLMS_F[,-1],
           DLCO~1,
           sigma.formula = DLCO~1,
           nu.formula = DLCO~1,
           family=BCCG(mu.link="log"),
           trace=FALSE,
           parallel="multicore",
           ncpus = nC)

mf2=stepGAICAll.A(mf1,scope=list(lower=~1, 
                                 upper=~poly(Age,3)+poly(Height,3)+pb(Age)),
                  sigma.scope = list(lower=~1, 
                                     upper=~poly(Age,3)+pb(Age)),
                  nu.scope = list(lower=~1,upper=~Age),
                  k=log(length(trainLMS_F[,-1])),
                  trace=FALSE,
                  parallel="multicore",
                  method=RS(10000),
                  ncpus = nC
)

## --------------------------------------------------- 
## Start:  AIC= 1549.39 
##  DLCO ~ 1 
## 
## --------------------------------------------------- 
## Start:  AIC= 1375.57 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 1375.57 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 1375.57 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 1375.57 
##  DLCO ~ poly(Age, 3) + poly(Height, 3) 
## 
## ---------------------------------------------------

modDL_F=mf2

summary(modDL_F)

## ******************************************************************
## Family:  c("BCCG", "Box-Cox-Cole-Green") 
## 
## Call:  
## gamlss(formula = DLCO ~ poly(Age, 3) + poly(Height, 3), sigma.formula = ~1,  
##     nu.formula = ~1, family = BCCG(mu.link = "log"), data = trainLMS_F[,  
##         -1], trace = FALSE, parallel = "multicore", ncpus = nC) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.191868   0.009995 319.337  < 2e-16 ***
## poly(Age, 3)1    -2.103494   0.167392 -12.566  < 2e-16 ***
## poly(Age, 3)2    -0.760809   0.159128  -4.781 3.04e-06 ***
## poly(Age, 3)3    -0.014776   0.159795  -0.092   0.9264    
## poly(Height, 3)1  0.835098   0.170029   4.911 1.67e-06 ***
## poly(Height, 3)2 -0.138109   0.157733  -0.876   0.3821    
## poly(Height, 3)3  0.261908   0.158388   1.654   0.0995 .  
## ---
## 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) -1.84908    0.04499   -41.1   <2e-16 ***
## ---
## 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.2648     0.3056   0.866    0.387
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  247 
## Degrees of Freedom for the fit:  9
##       Residual Deg. of Freedom:  238 
##                       at cycle:  4 
##  
## Global Deviance:     1363.093 
##             AIC:     1381.093 
##             SBC:     1412.677 
## ******************************************************************

# VA for Male

mm1=gamlss(data=trainLMS_M[,-1],
          AV~1,
          sigma.formula = AV~1,
          nu.formula = AV~1,
          family=BCCG(mu.link="log"),
          trace=FALSE,
          parallel="multicore",
          ncpus = nC)

mm2=stepGAICAll.A(mm1, scope=list(lower=~1, 
                                upper=~poly(Age,3)+poly(Height,3)+pb(Age)),
                 sigma.scope = list(lower=~1, 
                                    upper=~poly(Age,3)+poly(Height,3)+pb(Age)),
                 nu.scope = list(lower=~1, 
                                    upper=~poly(Age,3)+poly(Height,3)+pb(Age)),
                 k=log(length(trainLMS_M[,-1])),
                 trace=FALSE,
                 parallel="multicore",
                 method=RS(10000),
                 ncpus = nC
)

## --------------------------------------------------- 
## Start:  AIC= 694.71 
##  AV ~ 1 
## 
## --------------------------------------------------- 
## Start:  AIC= 559.96 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 559.38 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 549.6 
##  ~pb(Age) 
## 
## --------------------------------------------------- 
## Start:  AIC= 549.42 
##  AV ~ poly(Height, 3) + pb(Age) 
## 
## ---------------------------------------------------

modAV_M=mm2

summary(modAV_M)

## ******************************************************************
## Family:  c("BCCG", "Box-Cox-Cole-Green") 
## 
## Call:  
## gamlss(formula = AV ~ poly(Height, 3) + pb(Age), sigma.formula = ~1,  
##     nu.formula = ~poly(Height, 3) + pb(Age), family = BCCG(mu.link = "log"),  
##     data = trainLMS_M[, -1], trace = FALSE, parallel = "multicore",  
##     ncpus = nC) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.9021994  0.0144106 132.000  < 2e-16 ***
## poly(Height, 3)1  1.4046132  0.0783293  17.932  < 2e-16 ***
## poly(Height, 3)2 -0.0772345  0.0362332  -2.132  0.03403 *  
## poly(Height, 3)3 -0.3997339  0.0397826 -10.048  < 2e-16 ***
## pb(Age)           0.0009241  0.0003098   2.983  0.00314 ** 
## ---
## 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.26702    0.04446  -50.99   <2e-16 ***
## ---
## 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)        2.08298    1.31278   1.587   0.1139    
## poly(Height, 3)1  55.36866    7.46891   7.413 1.94e-12 ***
## poly(Height, 3)2 -15.75267    7.29528  -2.159   0.0318 *  
## poly(Height, 3)3  68.13854    7.30389   9.329  < 2e-16 ***
## pb(Age)           -0.04297    0.02934  -1.465   0.1443    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## 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 may not be reliable. 
## ------------------------------------------------------------------
## No. of observations in the fit:  253 
## Degrees of Freedom for the fit:  13.11931
##       Residual Deg. of Freedom:  239.8807 
##                       at cycle:  12 
##  
## Global Deviance:     531.2316 
##             AIC:     557.4702 
##             SBC:     603.8259 
## ******************************************************************

# AV for Female

mf1=gamlss(data=trainLMS_F[,-1],
           AV~1,
           sigma.formula = AV~1,
           nu.formula = AV~1,
           family=BCCG(mu.link="log"),
           trace=FALSE,
           parallel="multicore",
           ncpus = nC)

mf2=stepGAICAll.A(mm1, scope=list(lower=~1, 
                                  upper=~poly(Age,3)+poly(Height,3)+pb(Age)),
                  sigma.scope = list(lower=~1, 
                                     upper=~poly(Age,3)+poly(Height,3)+pb(Age)),
                  nu.scope = list(lower=~1, 
                                  upper=~poly(Age,3)+poly(Height,3)+pb(Age)),
                  k=log(length(trainLMS_F[,-1])),
                  trace=FALSE,
                  parallel="multicore",
                  method=RS(10000),
                  ncpus = nC
)

## --------------------------------------------------- 
## Start:  AIC= 694.71 
##  AV ~ 1 
## 
## --------------------------------------------------- 
## Start:  AIC= 559.96 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 559.38 
##  ~1 
## 
## --------------------------------------------------- 
## Start:  AIC= 549.6 
##  ~pb(Age) 
## 
## --------------------------------------------------- 
## Start:  AIC= 549.42 
##  AV ~ poly(Height, 3) + pb(Age) 
## 
## ---------------------------------------------------

modAV_F=mf2

summary(modAV_F)

## ******************************************************************
## Family:  c("BCCG", "Box-Cox-Cole-Green") 
## 
## Call:  
## gamlss(formula = AV ~ poly(Height, 3) + pb(Age), sigma.formula = ~1,  
##     nu.formula = ~poly(Height, 3) + pb(Age), family = BCCG(mu.link = "log"),  
##     data = trainLMS_M[, -1], trace = FALSE, parallel = "multicore",  
##     ncpus = nC) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.9021994  0.0144106 132.000  < 2e-16 ***
## poly(Height, 3)1  1.4046132  0.0783293  17.932  < 2e-16 ***
## poly(Height, 3)2 -0.0772345  0.0362332  -2.132  0.03403 *  
## poly(Height, 3)3 -0.3997339  0.0397826 -10.048  < 2e-16 ***
## pb(Age)           0.0009241  0.0003098   2.983  0.00314 ** 
## ---
## 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.26702    0.04446  -50.99   <2e-16 ***
## ---
## 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)        2.08298    1.31278   1.587   0.1139    
## poly(Height, 3)1  55.36866    7.46891   7.413 1.94e-12 ***
## poly(Height, 3)2 -15.75267    7.29528  -2.159   0.0318 *  
## poly(Height, 3)3  68.13854    7.30389   9.329  < 2e-16 ***
## pb(Age)           -0.04297    0.02934  -1.465   0.1443    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## 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 may not be reliable. 
## ------------------------------------------------------------------
## No. of observations in the fit:  253 
## Degrees of Freedom for the fit:  13.11931
##       Residual Deg. of Freedom:  239.8807 
##                       at cycle:  12 
##  
## Global Deviance:     531.2316 
##             AIC:     557.4702 
##             SBC:     603.8259 
## ******************************************************************

7 Mô hình mạng thần kinh nhân tạo sâu (Deep learning)

Cùng lúc này, cô chuyên viên Data science đang giải quyết vấn đề theo một hướng đi hoàn toàn khác với cách làm truyền thống. Sau đây là những phân tích của cá nhân cô ấy:

1. Chúng ta cắt bài toán hồi quy này ra thành nhiều mảnh nhỏ (mô hình riêng cho Nam, Nữ, AV, DLCO…) để làm gì cho khó khăn ? Tại sao ta không tạo ra 1 mô hình DUY NHẤT cho phép ước lượng đồng thời cả 2 outcomes (AV, DLCO),và cho tất cả mọi cá thể bất kể giới tính ?

2. Ta đã xác định về dữ liệu đầu vào, sẽ gồm 3 biến: Sex, Age, Height, ta cũng đã xác định sẽ có 2 outcomes là AV và DLCO. Như vậy mục tiêu của ta là tạo ra một nền tảng cho phép xuất ra 2 giá trị AV và DLCO gần nhất so với giá trị thực tế. Ta còn muốn rằng nền tảng này sẽ cho phép cái máy tự nó thích nghi với dữ liệu đầu vào mới, mà không cần phải tháo dỡ hệ thống ra rồi lắp ráp lại từ đầu như phương pháp Step-wise thủ công. Điều này giống như việc thiết kế một bệnh viện để nó hoạt động trong 10-20 năm, đầu vào tiếp nhận bệnh nhân, đầu ra bệnh nhân lành bệnh. Hệ thống bên trong bệnh viện (phòng ốc, hành lang, y cụ, nhân viên y tế) là cố định, nhưng có khả năng thích ứng với nhiều bệnh lý khác nhau, chữa nhiều bệnh đồng thời. Như vậy khi gặp bệnh khó chữa, vấn đề là hệ thống phải thích nghi chứ ta không thể đập bỏ bệnh viện và xây lại một cái mới.

Hệ thống đó chính là một mạng neuron nhân tạo có cấu trúc gồm 3 phần:

Lớp đầu tiên gồm 3 perceptron, có chức năng nhận dữ liệu đầu vào, tương ứng 3 biến là Sex (giá trị nhị phân), Age và Height.

Tín hiệu đi qua mỗi neuron sẽ được kích hoạt bằng một hàm để điều chỉnh trọng số của tín hiệu xuất ra qua những neuron tiếp theo nằm trong lớp sâu (hidden layers). Mạng neuron này sẽ có 2 hay 3 lớp sâu như vậy, mỗi lớp chứa vài chục đến vài trăm nodes. Trước khi dẫn đến lớp cuối cùng, tín hiệu sẽ đi qua một lớp có vai trò hiệu chính sử dụng 1 hàm l2 regularisation.

Lớp cuối cùng là đầu ra của mạng neuron, sẽ gồm 2 nodes tương ứng với 2 outcomes là DLCO và AV.

Cấu trúc này sẽ được “huấn luyện” bằng từng đoạn dữ liệu trong trainset và kiểm định qua hàng trăm lượt, cho đến khi đạt trạng thái tối ưu. Quá trình huấn luyện này có mục tiêu là giảm thiểu một hàm mất mát (loss function) và tiêu chí về sai số chẩn đoán, thí dụ MAE hay RMSE…

Sau khi huấn luyện xong, ta có trong tay 1 mô hình Deep neural network có khả năng ước lượng đồng thời AV và DLCO từ 3 biến đầu vào.

Công cụ được sử dụng là package keras, giao thức cho phép R tương tác với Tensofflow của Python (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). Thí dụ vector là 1D tensor, còn matrix gồm nhiều features, nhiều instances đượ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. Keras tiếp nhận riêng tensors cho tập features và cho biến kết quả nên ta phải tách riêng 2 phần này từ tập train và test

Cần lưu ý là Neural network chỉ hoạt động tối ưu khi dữ liệu đầu vào được chuẩn hó, và nó chỉ tiếp nhận dữ liệu số - nên các predictor như Sex phải được chuyển thành dummy variable hay binary variable (giá trị 0/1)

Việc chuẩn hóa này được thực hiện trên cả 2 tập trainset và testset, tuy nhiên chỉ dựa vào mean và sd của trainset để đảm bảo tính bí mật của testset. Package recipes cho phép ta tạo 1 hàm chuẩn hóa như vậy

library(recipes)

# Recipe for standardisation

dlcorec<-recipe(trainset,AV+DLCO~Sex+Age+Height)%>%
  add_role(Sex,new_role = "stratify")%>%
  step_center(all_predictors())%>%
  step_scale(all_predictors())

# Standardised trainset+testset

std_func<-prep(dlcorec,training=trainset,retain=F)
std_test<-bake(std_func, newdata=testset)
std_train<-bake(std_func, newdata=trainset)

# Predictors to tensors

train_data<-std_train%>%
  dplyr::select(-DLCO,-AV)%>%
  mutate(Sex=as.numeric(trainset$Sex)-1)%>%
  as.matrix()%>%
  as.array.default(dimnames=NULL)  

test_data<-std_test%>%
  dplyr::select(-DLCO,-AV)%>%
  mutate(Sex=as.numeric(testset$Sex)-1)%>%
  as.matrix()%>%
  as.array.default(dimnames=NULL)  

# Outcomes to tensors
train_targets<-trainset%>%
  dplyr::select(AV,DLCO)%>%
  as.matrix()%>%
  as.array.default(dimnames=NULL)  

test_targets<-testset%>%
  dplyr::select(AV,DLCO)%>%
  as.matrix()%>%
  as.array.default(dimnames=NULL)

Ta thiết kế mạng neuron trong keras: mạng này có 1 lớp đầu vào , 4 lớp sâu dùng hàm rectified linear unit (relu), 1 lớp hiệu chỉnh với L2 reg, cuối cùng là lớp kết quả gồm 2 neuron. Hàm mất mát là MAE, tiêu chí huấn luyện là MSE

# Design a deep neural net

# ANN

library(keras)

build_model <- function() {
  model <- keras_model_sequential() %>% 
    layer_dense(units = 128, activation = "relu", 
                input_shape = dim(train_data)[[2]]) %>% 
    layer_dense(units = 128, activation = "relu") %>%
    layer_dense(units = 128, activation = "relu") %>%
    layer_dense(units =128, kernel_regularizer = regularizer_l2(0.001))%>%   
    layer_dense(units = 2)
    
  model %>% compile(
    optimizer = "rmsprop",
    loss = "mae",
    metric= "msle"
  )
}

# Train ANN on whole data.

model <- build_model()

model %>% fit(train_data, train_targets,
              epochs = 200, batch_size = 30,verbose=0,
              validation_split = 0.1) -> dnnmod

summary(model)

## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #     
## ===========================================================================
## dense_1 (Dense)                  (None, 128)                   512         
## ___________________________________________________________________________
## dense_2 (Dense)                  (None, 128)                   16512       
## ___________________________________________________________________________
## dense_3 (Dense)                  (None, 128)                   16512       
## ___________________________________________________________________________
## dense_4 (Dense)                  (None, 128)                   16512       
## ___________________________________________________________________________
## dense_5 (Dense)                  (None, 2)                     258         
## ===========================================================================
## Total params: 50,306
## Trainable params: 50,306
## Non-trainable params: 0
## ___________________________________________________________________________

plot(dnnmod)+theme_bw()

Sau khi huấn luyện 200 lượt với tỉ lệ kiểm định 30%, sử dụng ngẫu nhiên 30 trường hợp mỗi lượt, ta có kết quả sau cùng trong object dnnmod.

Phần lớn thời gian được sử dụng là để thiết kế cấu trúc mạng neuron nhưng công việc này tương đối nhẹ nhàng so với quy trình step wise thủ công mà phòng thống kê đang làm.

8 So sánh hiệu quả 2 phương pháp

Đây là thời điểm quyết định, hai mô hình LMS (phái thống kê) và Deep neural network (phái Machine learning) đã hoàn tất. Ta sẽ kiểm tra độ chính xác của chúng trên tập testset.

Bên cạnh đó, một mô hình tuyến tính đơn giản cho 2 chỉ số DLCO và AV được ERS đề xuất năm 2017 sẽ được dùng làm nhóm chứng,chúng ta hy vọng rằng có sự khác biệt giữa 2 mô hình mới và mô hình tuyến tính, cũng như có sự tương phản rõ nét về hiệu năng giữa 2 phương pháp Deep learning và LMS theo GAMLSS.

# Validating DLCO prediction on testset

y_pred=predict(model,test_data)

library(mlr)

# Metrics to be estimated
mets=list(mse,mae,medae,rmse)

dlco.task= mlr::makeRegrTask(id = "DLCO", data=testset, target = "DLCO")
dlco.lrn = makeLearner("regr.glm")
dummydlco=mlr::train(dlco.lrn,dlco.task)
dumpredDNNdlco=predict(dummydlco,dlco.task)

# Implement y_pred to dummy
dumpredDNNdlco$data$response<-y_pred[,2]


# Validating VA prediction on testset

va.task= mlr::makeRegrTask(id = "VA", data=testset, target = "AV")
va.lrn = makeLearner("regr.glm")
dummyva=mlr::train(va.lrn,va.task)
dumpredDNNva=predict(dummyva,va.task)

# Implement y_pred to dummy
dumpredDNNva$data$response<-y_pred[,1]

# Validating DLCO by LMS

pred_DL_M=predict(modDL_M,newdata=testLMS_M[,-1],type="response")
pred_DL_F=predict(modDL_F,newdata=testLMS_F[,-1],type="response")
pred_AV_M=predict(modAV_M,newdata=testLMS_M[,-1],type="response")

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

pred_AV_F=predict(modAV_F,newdata=testLMS_F[,-1],type="response")

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

test_DL_LMS=rbind(testLMS_M,testLMS_F)%>%
  mutate(predDL=c(pred_DL_M,pred_DL_F),
         predAV=c(pred_AV_M,pred_AV_F))

# Dummy LMS
dumpredLMS_DL=predict(dummydlco,dlco.task)
dumpredLMS_AV=predict(dummyva,va.task)

# Implement

dumpredLMS_DL$data$response<-test_DL_LMS$predDL
dumpredLMS_AV$data$response<-test_DL_LMS$predAV

            
# ERS 2017 linear model:

# DLCO
# Male: (0.3*HEIGHT)-(0.002*(AGE^2))+(5.47*1)-17.75    RSD=4.4377
# Female: (0.3*HEIGHT)-(0.002*(AGE^2))+(5.47*0)-17.75  RSD=4.7377

#AV
# male: 0.082*HEIGHT+0.72*1-8.2   RSD=0.766
# Female: 0.082*HEIGHT+0.72*0-8.2 RSD=0.766

predERSmod=function(newdata){
  Sex=as.numeric(newdata$Sex)-1
  Age=newdata$Age
  Height=newdata$Height
  rsdDL=4.4377
  rsdAV=0.766
  predERS=list(Mu_DL=(0.3*Height)-(0.002*(Age^2))+(5.47*Sex)-17.75,
               LLN_DL=(0.3*Height)-(0.002*(Age^2))+(5.47*Sex)-17.75-1.645*rsdDL,
               ULN_DL=(0.3*Height)-(0.002*(Age^2))+(5.47*Sex)-17.75+1.645*rsdDL,
               Mu_AV=0.082*Height+0.72*Sex-8.2,
                LLN_AV=0.082*Height+0.72*Sex-8.2-1.645*rsdAV,
                ULN_AV=0.082*Height+0.72*Sex-8.2+1.645*rsdAV)
return(predERS)
}

predERS=predERSmod(testset)

# Dummy ERS
dumpredERS_DL=predict(dummydlco,dlco.task)
dumpredERS_AV=predict(dummyva,va.task)

# Implement

dumpredERS_DL$data$response<-predERS$Mu_DL
dumpredERS_AV$data$response<-predERS$Mu_AV

# Compare DLCO
pe1=performance(dumpredLMS_DL,measures =mets)
pe2=performance(dumpredDNNdlco,measures =mets)
pe3=performance(dumpredERS_DL,measures =mets)

# Compare AV
pe4=performance(dumpredLMS_AV,measures =mets)
pe5=performance(dumpredDNNva,measures =mets)
pe6=performance(dumpredERS_AV,measures =mets)

verdict=data_frame(Outcome=c("DLCO","DLCO","DLCO","AV","AV","AV"),
                   Method=rep(c("Step_wise_LMS","Deep_NN","ERS_linear"),2),
                   MSE=c(pe1[[1]],pe2[[1]],pe3[[1]],pe4[[1]],pe5[[1]],pe6[[1]]),
                   MAE=c(pe1[[2]],pe2[[2]],pe3[[2]],pe4[[2]],pe5[[2]],pe6[[2]]),
                   MEDAE=c(pe1[[3]],pe2[[3]],pe3[[3]],pe4[[3]],pe5[[3]],pe6[[3]]),
                   RMSE=c(pe1[[4]],pe2[[4]],pe3[[4]],pe4[[4]],pe5[[4]],pe6[[4]])
                   )

verdict%>%knitr::kable()

Outcome	Method	MSE	MAE	MEDAE	RMSE
DLCO	Step_wise_LMS	123.9756040	9.8277355	9.607005	11.1344333
DLCO	Deep_NN	21.2717903	3.5261171	2.852706	4.6121351
DLCO	ERS_linear	22.9205832	4.0352126	3.838000	4.7875446
AV	Step_wise_LMS	2.8022101	1.4087842	1.304354	1.6739803
AV	Deep_NN	0.5301520	0.5686747	0.454293	0.7281154
AV	ERS_linear	0.5583332	0.6060008	0.536000	0.7472170

Kết quả kiểm định cho thấy một điều bất ngờ, đó là một mô hình phức tạp như LMS lại có mức độ chính xác kém hơn so với một mô hình tuyến tính đơn giản (ERS). Mô hình Deep learning có hiệu quả tương đương, thậm chí có thể tốt hơn so với mô hình của ERS. Kết quả được trình bày bằng biểu đồ:

verdict%>%gather(MSE:RMSE,key="Metric",value="Score")%>%
  ggplot(aes(x=Method,y=Score,fill=Method))+
  geom_bar(stat="identity",position="dodge",alpha=0.6,col="black")+
  theme_bw()+coord_flip()+
  facet_grid(Metric~Outcome,scales="free")

Độ chính xác của mô hình Deep learning và ERS có thể được quan sát trên sự chồng lắp biểu đồ tuyến kí giữa dữ liệu ước tính bởi mô hình và dữ liệu thực tế tron tập kiểm định như sau:

# LLN and ULN

# rsd for DLCO

(trainset$DLCO-as.vector(y_pred[,2]))%>%
  density()%>%
  plot(col="red")

rsd_DL=sd(trainset$DLCO-as.vector(y_pred[,2]))

# rsd for AV

(trainset$AV-as.vector(y_pred[,1]))%>%
  density()%>%
  plot(col="blue")

rsd_AV=sd(trainset$AV-as.vector(y_pred[,1]))


# LLN and ULN by LMS

pred_DL_M_all=predictAll(modDL_M,newdata=testLMS_M[,-1])
pred_DL_F_all=predictAll(modDL_F,newdata=testLMS_F[,-1])
pred_AV_M_all=predictAll(modAV_M,newdata=testLMS_M[,-1])

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

pred_AV_F_all=predictAll(modAV_F,newdata=testLMS_F[,-1])

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

# DLCO
MuDL=c(pred_DL_M_all$mu , pred_DL_F_all$mu)
sigmaDL=c(pred_DL_M_all$sigma , pred_DL_F_all$sigma)
NuDL=c(pred_DL_M_all$nu , pred_DL_F_all$nu)
LLN_DL_LMS=exp(log(MuDL)+(log(1-1.645*NuDL*sigmaDL))/NuDL)
ULN_DL_LMS=exp(log(MuDL)+(log(1+1.645*NuDL*sigmaDL))/NuDL)

# AV
MuAV=c(pred_AV_M_all$mu , pred_AV_F_all$mu)
sigmaAV=c(pred_AV_M_all$sigma , pred_AV_F_all$sigma)
NuAV=c(pred_AV_M_all$nu , pred_AV_F_all$nu)
LLN_AV_LMS=exp(log(MuAV)+(log(1-1.645*NuAV*sigmaAV))/NuAV)
ULN_AV_LMS=exp(log(MuAV)+(log(1+1.645*NuAV*sigmaAV))/NuAV)

test_all<-testset%>%mutate(Mu_AV_ANN=as.vector(y_pred[,1]),
                          LLN_AV_ANN=qnorm(p=0.2,
                                           mean=as.vector(y_pred[,1]),
                                           sd=rsd_AV),
                          ULN_AV_ANN=qnorm(p=0.8,
                                           mean=as.vector(y_pred[,1]),
                                           sd=rsd_AV),
                          Mu_DL_ANN=as.vector(y_pred[,2]),
                          LLN_DL_ANN=qnorm(p=0.2,
                                           mean=as.vector(y_pred[,2]),
                                           sd=rsd_DL),
                          ULN_DL_ANN=qnorm(p=0.8,
                                           mean=as.vector(y_pred[,2]),
                                           sd=rsd_DL),
                          Mu_DL_LMS=MuDL,
                          LLN_DL_LMS=LLN_DL_LMS,
                          ULN_DL_LMS=ULN_DL_LMS,
                          Mu_AV_LMS=MuAV,
                          LLN_AV_LMS=LLN_AV_LMS,
                          ULN_AV_LMS=ULN_AV_LMS,
                          Mu_DL_ERS=predERS$Mu_DL,
                          LLN_DL_ERS=predERS$LLN_DL,
                          ULN_DL_ERS=predERS$ULN_DL,
                          Mu_AV_ERS=predERS$Mu_AV,
                          LLN_AV_ERS=predERS$LLN_AV,
                          ULN_AV_ERS=predERS$ULN_AV
                          )
                          
p1=test_all%>%ggplot(aes(x=Age))+
  geom_point(aes(y=AV),alpha=0.3,col="black")+
  geom_smooth(aes(y=Mu_AV_ANN),fill="red",col="red3",alpha=0.2)+
  geom_smooth(aes(y=AV),fill="green",col="green4",alpha=0.2)+
  geom_smooth(aes(y=LLN_AV_ANN),fill="red",col="red3",se=F,linetype=2)+
  geom_smooth(aes(y=ULN_AV_ANN),fill="red",col="red3",se=F,linetype=2)+
  theme_bw()+ggtitle("Deep neural network")+scale_y_continuous(limits = c(2,9))

p2=test_all%>%ggplot(aes(x=Age))+
  geom_point(aes(y=AV),alpha=0.3,col="black")+
  geom_smooth(aes(y=Mu_AV_LMS),fill="blue",col="blue3",alpha=0.2)+
  geom_smooth(aes(y=AV),fill="green",col="green3",alpha=0.2)+
  geom_smooth(aes(y=LLN_AV_LMS),fill="blue",col="blue3",se=F,linetype=2)+
  geom_smooth(aes(y=ULN_AV_LMS),fill="blue",col="blue3",se=F,linetype=2)+
  theme_bw()+ggtitle("LMS")+scale_y_continuous(limits = c(2,9))

p3=test_all%>%ggplot(aes(x=Age))+
  geom_point(aes(y=AV),alpha=0.3,col="black")+
  geom_smooth(aes(y=Mu_AV_ERS),fill="purple",col="blue3",alpha=0.2)+
  geom_smooth(aes(y=AV),fill="green",col="green3",alpha=0.2)+
  geom_smooth(aes(y=LLN_AV_ERS),fill="purple",col="purple",se=F,linetype=2)+
  geom_smooth(aes(y=ULN_AV_ERS),fill="purple",col="purple",se=F,linetype=2)+
  theme_bw()+ggtitle("ERS")+scale_y_continuous(limits = c(2,9))

p4=test_all%>%ggplot(aes(x=Age))+
  geom_point(aes(y=DLCO),alpha=0.3,col="black")+
  geom_smooth(aes(y=DLCO),fill="green",col="green4",alpha=0.2)+
  geom_smooth(aes(y=Mu_DL_ANN),fill="red",col="red3",alpha=0.2)+
  geom_smooth(aes(y=LLN_DL_ANN),fill="red",col="red3",se=F,linetype=2)+
  geom_smooth(aes(y=ULN_DL_ANN),fill="red",col="red3",se=F,linetype=2)+
  theme_bw()+scale_y_continuous(limits = c(5,50))+ggtitle("Deep neural network")

p5=test_all%>%ggplot(aes(x=Age))+
  geom_point(aes(y=DLCO),alpha=0.3,col="black")+
  geom_smooth(aes(y=DLCO),fill="green",col="green4",alpha=0.2)+
  geom_smooth(aes(y=Mu_DL_LMS),fill="blue",col="blue3",alpha=0.2)+
  geom_smooth(aes(y=LLN_DL_LMS),fill="blue",col="blue3",se=F,linetype=2)+
  geom_smooth(aes(y=ULN_DL_LMS),fill="blue",col="blue3",se=F,linetype=2)+
  theme_bw()+scale_y_continuous(limits = c(5,50))+ggtitle("LMS")

p6=test_all%>%ggplot(aes(x=Age))+
  geom_point(aes(y=DLCO),alpha=0.3,col="black")+
  geom_smooth(aes(y=DLCO),fill="green",col="green4",alpha=0.2)+
  geom_smooth(aes(y=Mu_DL_ERS),fill="violet",col="purple",alpha=0.2)+
  geom_smooth(aes(y=LLN_DL_ERS),fill="violet",col="purple",se=F,linetype=2)+
  geom_smooth(aes(y=ULN_DL_ERS),fill="violet",col="purple",se=F,linetype=2)+
  theme_bw()+scale_y_continuous(limits = c(5,50))+ggtitle("ERS")

gridExtra::grid.arrange(p1,p2,p3,ncol=2)

gridExtra::grid.arrange(p4,p5,p6,ncol=2)

9 Bàn luận

Sự phổ biến và phát triển của các phương pháp Machine learning mà giải thuật hiện đại nhất là Deep learning hay mạng thần kinh nhân tạo đã cung cấp một giải pháp hoàn toàn mới cho bài toán xây dựng giá trị tham chiếu cho xét nghiệm y khoa. Trong thí dụ ngắn này, Nhi đã chứng tỏ rằng Deep learning với cấu trúc mạng neuro phù hợp sẽ mạnh hơn nhiều so với những mô hình tuyến tính cổ điển, thậm chí những mô hình phức tạp nhất như LMS với Splines.

Ta cùng nhìn lại 2 giải pháp để so sánh giữa chúng:

1. Những ưu thế của Deep learning so với mô hình LMS:

• Deep neural network cho phép giải quyết cùng lúc nhiều outcomes trong một mô hình duy nhất, LMS không thể làm được điều này.

• Cấu trúc của mạng neuron chỉ cần thiết kế 1 lần, nhưng có thể tái sử dụng nhiều lần , và việc thiết kế có tính thứ bậc. Trái lại, phương pháp step-wise luôn phá bỏ cấu trúc mô hình cũ và thay thế bằng cấu trúc mới khi thay đổi dữ liệu đầu vào. Việc thăm dò, thử, sai và thử lại những tổ hợp của bậc đa thức trong mô hình LMS bằng Stepwise cũng khó khăn hơn nhiều so với việc dựng mạng neuron.

2. Những nhược điểm của Deep learning so với LMS:

• Mô hình Deep neural network là một hộp đen (blackbox), ta không thể nhìn thấy cơ chế bên trong của nó. Tuy nhiên, điều này không thực sự quan trọng. Trên thực tế các mô hình LMS cũng là 1 dạng blackbox, vì không có bác sĩ nào có khả năng tính toán thủ công từ mô hình, nhất là hàm splines sẽ tạo ra những matrices phức tạp và hoàn toàn ẩn. Một khi đã đưa vào software, người ta không quan tâm đến nội dung của model nữa.

• Mô hình Deep neural net chưa cho phép ước tính phân vị một cách chính xác như LMS. Trong thí dụ này, Nhi ước tính LLN và ULN dựa vào giả định là residual error có phân bố chuẩn, nên LLN được tính từ hàm quantile với Mean= predicted và sd = sd của residual error.

Một câu hỏi thú vị nữa: Tại sao ta không dùng Deep learning để đưa ra chẩn đoán luôn (phân biệt các bệnh lý với người bình thường) mà lại phải đi vòng qua mô hình hồi quy ? Câu trả lời đó là : (1) bài toán multiclass, multilabel classification luôn phức tạp hơn nhiều so với 1 bài toán hồi quy, và (2) Thu thập dữ liệu của bệnh nhân với chẩn đoán xác định thì khó khăn hơn rất nhiều so với chỉ thu thập dữ liệu trên người khỏe mạnh, (3) Đọc kết quả xét nghiệm dựa vào giá trị tham chiếu vẫn còn là phương pháp thực hành lâm sàng quy ước, các bác sĩ không thấy dễ chịu với ý tưởng cái máy thay thế vai trò của họ hoàn toàn.

Có thể trong một tương lai gần, máy tính sẽ thay thế cho vai trò của chuyên viên thống kê, vì trí tuệ nhân tạo mang lại ưu thế quá lớn so với các làm việc thủ công. Thay vì phải chờ 3-5 năm để có một mô hình, một cái máy có thể làm việc này hàng tuần, thậm chí hằng ngày, với độ chính xác cao hơn nhiều.

Tài liệu tham khảo

1. Rigby RA, Stasinopoulos DM. Generalized additive models for location, scale and shape (with discussion). Appl Statist 2005; 54: 507-554.

2. Francois Chollet, J. J. Allaire. Deep Learning with R. Publisher: Manning Publications. 2018. ISBN: 9781617295546