Hướng dẫn GAMLSS bài 7

Mô hình Lambda-Mu-Sigma (LMS)

1 Giới thiệu về mô hình Lambda-Mu-Sigma

Thân chào các bạn, trong bài thực hành số 7 cho package GAMLSS này, Nhi sẽ giới thiệu với các bạn một ứng dụng thiết thực và độc đáo của Gamlss vào Y học lâm sàng, đó là phương pháp mô hình 3 tham số (Lambda-Mu-Sigma, LMS)

Phương pháp LMS do 2 tác giả Rigby RA. và Stasinopoulos DM thiết kế vào năm 2001 như giải pháp chuyên dụng để ước lượng một biến số ngẫu nhiên, liên tục và không âm (thường là đại lượng sinh học) đặc trưng bởi quy luật biến thiên phi tuyến tính theo thời gian (cụ thể là Tuổi).

LMS trở thành một công cụ ưu thế trong Y học chính vì 3 lý do :

1. Tính linh hoạt : LMS không yêu cầu Y phải có phân phối bình thường (Gaussian).

2. Hữu ích cho chẩn đoán lâm sàng : Từ giá trị của tham số L,M,S chúng ta có thể ước tính chính xác 2 ngưỡng Upper limit và Lower limit của Y.

Cũng từ L,M,S, ta có thể ước tính Z-score (thang chuẩn hóa) cho Y. Ngưỡng ULN,LLN và Z-score được ứng dụng trong việc diễn giải kết quả xét nghiệm của nhiều chỉ số sinh học.

3. LMS Cho phép mô hình hóa một cách liên tục sự biến thiên theo quỹ đạo phức tạp, phi tuyến tính của Y.

Một số ứng dụng của LMS trên thực tế như :

• Mô hình tăng trưởng ở trẻ em (Cân nặng, chiều cao, BMI…).

• Mô hình ước tính giá trị « bình thường » của các xét nghiệm y khoa, thí dụ : dung tích phổi, hàm lượng Hb, huyết áp, các marker sinh lý, sinh hóa, huyết học, miễn dịch …

Tất cả những bài toán này đều có liên quan đến yếu tố Thời gian (Tuổi).

Hình sau đây trình bày đồ thị của 4x3 mô hình ước tính giá trị lý thuyết của 3 đại lượng sinh lý hô hấp là FEV1, FVC, FEV1/FVC riêng cho Nam/Nữ và cho 4 chủng tộc khác nhau. Tất cả chúng đều được dựng nên nhờ phương pháp LMS và sử dụng package gamlss, dựa vào hàng trăm ngàn kết quả hô hấp ký trên người bình thường từ 3-95 tuổi.

Các mô hình này do nhóm Global Lung function initiative (GLI) công bố vào năm 2012.

Trong một thí dụ khác, một nhóm tác giả Hà Lan khảo sát chỉ số Gripstrength ở nam và nữ, từ trẻ em đến người già:

LMS có bản chất là một thể đặc biệt thuộc họ mô hình GAM :

1. Nền tảng của LMS model là họ phân phối Box-Cox-Cole-Green (BCCG), một quy luật áp dụng để mô tả biến số ngẫu nhiên Y nằm trong khoảng từ 0 đến dương vô cùng. Họ phân phối này do Cole và Green thiết lập năm 1992 từ phương pháp hoán chuyển Box-Cox. Phân phối 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).

Hình sau đây minh họa đồ thị hàm PDF của một biến X có phân phối BCCG(mu=5,sigma=0.5 và nu=1):

Trước hết, giả sử ta có 1 biến số sau hoán chuyển Z, có phân phối chuẩn với Mu=0 và Sigma (SD) = 1,

Thì biến số ngẫu nhiên Y tuân theo quy luật BCCG nếu :

Khi Nu khác 0 :

\[Z = \frac{1}{\sigma \nu }\left [ \left ( \frac{Y}{\mu }\right )^{\nu } -1\right ]\] Khi Nu=0 :

\[Z = \frac{1}{\sigma }log\left ( \frac{Y}{\mu } \right )\] 2) 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, cân nặng trong mô hình dự báo dung tích phổi…

3. 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). Thí dụ nếu ta áp dụng hàm poly(Age,2) thì X sẽ bao gồm Age và Age^2).Hàm Spline phổ biến nhất là B_spline. Các bạn đọc lại những bài trước để rõ hơn về design matrix và các thành phần đa thức bậc cao này.

Chính những hàm polynomial và smoothing đó 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.

2 Minh họa: Ngưỡng giá trị bình thường của Dung tích phổi

Trong bài, Nhi sẽ tái hiện lại một quy trình dựng mô hình dự báo giá trị lý thuyết cho một đại lượng sinh lý (Dung tích phổi: Total Lung capacity) theo Tuổi và Chiều cao ở Nam giới, chủng tộc Caucassian.

Trước hết chúng ta tải dữ liệu vào R:

library(tidyverse)

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

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

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

Đây là một nghiên cứu dịch tễ học điển hình. Mục tiêu của chúng ta là tạo ra một mô hình tiên lượng cho phép ước tính : giá trị trung bình lý thuyết của TLC, 2 ngưỡng cao nhất và thấp nhất của khoảng TLC bình thường. Nếu ta biết giá trị TLC thực tế của đối tượng thì ta có thể tính Z-score.

\[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 ]\]

\[Z_{score}=\frac{\mu^{\nu -1}}{\nu \sigma }\]

\[\%predicted = \frac{100*Y}{\mu }\]

Giải thích thêm một chút:

Xét nghiệm chức năng hô hấp là một trong những tiêu chuẩn quan trọng để chẩn đoán xác định các bệnh lý về hô hấp. Trong thực hành lâm sàng, kết quả xét nghiệm đo dung tích phổi được diễn giải theo 3 cách:

1. Có giá trị TLC đo được (Y), bác sĩ so sánh giá trị này với giá trị trung bình dự báo (mean predicted)- chính là Mu trong mô hình LMS. Ngưỡng cắt chẩn đoán thường là 80%, nếu %pred<80, ta có thể kết luận TLC giảm bất thường.

2. Bác sĩ so sánh TLC đo được (Y) với 2 giá trị LLN và ULN mà mô hình LMS tính ra. Nếu Y < LLN cho thấy TLC giảm bất thường, nếu Y > ULN cho thấy TLC tăng cao bất thường.

3. Từ Y,L,M và S bác sĩ tính được Z-score, nếu Z-score < -1.645 cho thấy TLC giảm bất thường, nếu Z-score < -1.9 hoặc < -2 thì gần như chắc chắn là có sự suy giảm dung tích phổi bệnh lý.

Hình dưới đây minh họa cho một báo cáo xét nghiệm hô hấp ký với 3 chỉ số chính là FVC, FEV1 và FEV1/FVC, mô hình GLI 2012 cho phép tính được LLN, Z-score và %Pred cho bệnh nhân. Kết quả cho thấy FEV1/FVC < LLN và Z-score < -3, gợi ý hội chứng tắc nghẽn và bệnh lý nghi ngờ là Hen phế quản vì bệnh nhân có hồi phục sau khi dùng thuốc dãn phế quản.

Trở lại với thí dụ của chúng ta, ta sẽ bắt đầu bằng việc thăm dò dữ liệu

3 Bước 1: Thăm dò dữ liệu

Bước thăm dò này rất quan trọng, nội dung của nó gồm:

• Xác định phạm vi áp dụng của mô hình: Giới hạn Min/Max của mỗi predictor(Chiều cao, Tuổi) cũng chính là giới hạn của mô hình. Thí dụ nếu người cao tuổi nhất là 75 tuổi thì bạn không thể dùng mô hình để dự báo LLN, Mu, Z_score cho ông già 90 tuổi. Tương tự, nếu chiều cao thấp nhất trong data là 1m60 thì bạn không thể áp dụng mô hình cho người thấp 1m50. Khi bạn cố tình làm trái với quy ước này, đó là sự ngoại suy và kết quả không thể tin cậy. Phần cuối bài Nhi sẽ minh họa về ngoại suy.

• Phát hiện outliers

• Phác họa quỹ đạo của mô hình trên biểu đồ, và thăm dò hiệu ứng của các hàm đa thức

Nhi sử dụng package tableone để làm thống kê mô tả: Ở nam giới, phạm vi của chiều cao là 158-200 cm, còn Tuổi là 18-85.

library(tableone)

vars=c("Age","Height","TLC")

sumtab=CreateContTable(vars, 
                      strata = "Sex",
                      funcNames=c("n","mean","median","sd","min","max"),
                      test=T,
                      data = df)

summary(sumtab,
        digits = 3)

## Sex: F
##          n   mean median    sd    min    max
## Age    239  44.75   43.0 17.53  18.00  87.00
## Height 239 163.97  163.0  7.17 147.00 182.00
## TLC    239   5.31    5.3  0.86   3.18   7.58
## -------------------------------------------------------- 
## Sex: M
##          n   mean median     sd    min    max
## Age    243  43.51  43.00 16.083  18.00  85.00
## Height 243 176.72 176.00  7.515 158.00 200.00
## TLC    243   7.07   7.05  0.969   4.61   9.82
## 
## p-values
##             pNormal   pNonNormal
## Age    4.176566e-01 6.521788e-01
## Height 1.023777e-60 1.036753e-49
## TLC    1.859500e-70 2.072124e-55
## 
## Standardize mean differences
##            1 vs 2
## Age    0.07387375
## Height 1.73666068
## TLC    1.92412053

Density curves cho thấy biến kết quả TLC có phân phối gần như bình thường ở cả Nam và Nữ:

df%>%gather(TLC:Height,key="Variable",value="Value")%>%
  ggplot()+
  geom_density(aes(x=Value,fill=Sex),alpha=0.6)+
  facet_wrap(~Variable,ncol=2,scales="free")+
  theme_bw()

Ta tách riêng ra một dataset dành cho Nam giới:

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

TLC có liên hệ phi tuyến tính với Chiều cao và tuổi, như đồ thị hàm loess smoothing cho thấy:

dfM%>%gather(Height,Age,key="Predictor",value="Value")%>%
  ggplot(aes(x=Value,y=TLC))+
  geom_point(shape=21,col="red",fill="red",alpha=0.3)+
  geom_smooth(alpha=0.5,fill="red",col="red4")+
  facet_wrap(~Predictor,ncol=2,scales="free")+
  theme_bw()

Nếu ta ước tính TLC theo Age và Height bằng một mô hình tuyến tính,đồ thị và mặt phẳng hồi quy sẽ như thế này:

dfM%>%gather(Height,Age,key="Predictor",value="Value")%>%
  ggplot(aes(x=Value,y=TLC))+
  geom_point(shape=21,col="black",fill="grey80",alpha=0.4)+
  geom_smooth(alpha=0.5,fill="blue",col="blue",method="lm")+
  facet_wrap(~Predictor,ncol=2,scales="free")+
  theme_bw()

library(mgcv)

fit <- gam(TLC~Age+Height,
           data=dfM)

vis.gam(fit,view=c("Age","Height"),theta=50,phi=30,n.grid=30,color="topo",ticktype = "detailed") 

Nếu áp dụng một hàm đa thức bậc 3 cho Height và bậc 2 cho Age,một linkfunction=log cho TLC, đồ thị và mặt phẳng hồi quy sẽ như sau:

dfM%>%ggplot(aes(x=Height,y=log(TLC)))+
  geom_point(alpha=0.3)+
  geom_smooth(method=glm,
              formula = y ~ poly(x,3),alpha=0.3,fill="violet",col="purple")+
  theme_bw()+ggtitle("ln(TLC)~Height^3+Height^2+Height")

dfM%>%ggplot(aes(x=Age,y=log(TLC)))+
  geom_point(alpha=0.3)+
  geom_smooth(method=glm,
              formula = y ~ poly(x,2),alpha=0.3,fill="violet",col="purple")+
  theme_bw()+ggtitle("ln(TLC)~Age^2+Age")

fit <- gam(log(TLC)~Age+poly(Age,2)+poly(Height,3),
           data=dfM)

vis.gam(fit,view=c("Age","Height"),theta=60,phi=25,n.grid=30,color="topo",ticktype = "detailed") 

4 Bước 2: Dựng mô hình LMS

Bây giờ chúng ta qua bước 2, dựng mô hình LMS trong gamlss. Trước hết, 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,]

4.1 Chiến lược 1: Triệt thoái tự động dựa vào BIC và AIC

gamlss cung cấp 2 giải pháp rất đơn giản để thăm dò và chọn lọc mô hình tối ưu. Cả 2 đều dựa vào tiêu chí AIC hoặc BIC, nhưng khác nhau về cách tái chọn mẫu / kiểm định.

Cách thứ nhất không dùng tái chọn mẫu, các mô hình được dựng và kiểm tra tuần tự trên CÙNG một mẫu duy nhất (thí dụ trainM). Các biến số trong formula sẽ bị triệt tiêu dần đến khi mô hình có AIC tối ưu hiện ra.

Quy trình Stepwise này thực ra phức tạp hơn ta tưởng, vì ở đây ta có đến 3 mô hình con là L,M,S; và một công thức với nhiều biến predictor, bao gồm smoothing spline, như vậy tất cả tổ hợp có thể giữa chúng phải được kiểm tra, bao gồm 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 
## ******************************************************************

Rsq(m2,type="both")

## $CoxSnell
## [1] 0.45431
## 
## $CraggUhler
## [1] 0.484216

plot(m2)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -0.000254358 
##                        variance   =  1.005155 
##                coef. of skewness  =  -0.01853313 
##                coef. of kurtosis  =  3.01777 
## Filliben correlation coefficient  =  0.9963444 
## ******************************************************************

Dù R2 khá thấp, nhưng mô hình tối ưu mà Stepwise chọn được không cho thấy vấn đề gì nghiêm trọng.

4.2 Chiến lược 2: Kiểm chứng chéo K block

Cách chọn lọc mô hình thứ hai phức tạp hơn, thường dùng để kiểm định kết quả của Stepwise. Nó có bản chất là một quy trình K-folds cross validation.

Nguyên tắc lựa chọn mô hình vẫn dựa trên tiêu chí AIC và BIC, tuy nhiên ta phải xác định trước những mô hình nào cần được kiểm tra, thí dụ ta có 3 mô hình M1, M2, M3. Quy trình crossvalidation sẽ được áp dụng cho mỗi mô hình, và tính AIC trung bình. Sau đó ta so sánh AIC này giữa 3 mô hình, mô hình nào có AIC thấp nhất sẽ được chọn.

Mỗi mô hình sẽ được dựng trên (k-1) block dữ liệu và kiểm định độc lập trên block còn lại, quy trình này lặp lại K lần, như vậy mỗi một block trong số K block đều được dùng làm tập kiểm định 1 lần và đều tham gia vào việc huấn luyện mô hình.

Thí dụ, Nhi đặt mô hình g1 với công thức tối ưu như trên:

TLC_{poly(Height,3)+Age+pb(Age), sigma.formula=TLC}Age+pb(Age),

Sau đó thay đổi công thức để tạo ra thêm 3 mô hình g2,g3,g4 đơn giản hơn, thí dụ giảm bậc đa thức của Height xuống còn 2, hay tăng bậc đa thức cho Age từ 1 lên 2, hay bỏ hẳn Age và chỉ giữ lại Height.

rand1 <- sample (5,nrow(trainM),replace=TRUE)

g1=gamlssCV(rand=rand1,
            data=trainM,
            TLC~poly(Height,3)+Age+pb(Age),
            sigma.formula=TLC~Age+pb(Age),
            family=BCCG(mu.link="log"),
            trace=FALSE,
            parallel="multicore",
            ncpus = nC
            )

## fold 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) 
## fold 2
## 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) 
## fold 3
## 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) 
## fold 4
## 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) 
## fold 5
## 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)

g2=gamlssCV(rand=rand1,
            data=trainM,
            TLC~poly(Height,3)+poly(Age,2)+pb(Age),
            sigma.formula=TLC~poly(Age,2)+pb(Age),
            family=BCCG(mu.link="log"),
            trace=FALSE,
            parallel="multicore",
            ncpus = nC
)

## fold 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) 
## fold 2
## 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) 
## fold 3
## 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) 
## fold 4
## 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) 
## fold 5
## 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)

g3=gamlssCV(rand=rand1,
            data=trainM,
            TLC~poly(Age,2)+pb(Age),
            sigma.formula=TLC~Age+pb(Age),
            family=BCCG(mu.link="log"),
            trace=FALSE,
            parallel="multicore",
            ncpus = nC)

## fold 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) 
## fold 2
## 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) 
## fold 3
## 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) 
## fold 4
## 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) 
## fold 5
## 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)

g4=gamlssCV(rand=rand1,
            data=trainM,
            TLC~poly(Height,3)+pb(Age),
            sigma.formula=TLC~Age+pb(Age),
            family=BCCG(mu.link="log"),
            trace=FALSE,
            parallel="multicore",
            ncpus = nC
)

## fold 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) 
## fold 2
## 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) 
## fold 3
## 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) 
## fold 4
## 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) 
## fold 5
## 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)

CV(g1,g2,g3,g4)

##    val[o.val]
## g1   462.3335
## g4   462.3335
## g2   470.9797
## g3   601.6367

Kết quả của CV khẳng định rằng mô hình G1 vẫn là tối ưu

5 Bước 3: Biện luận mô hình

Sau khi có mô hình tối ưu, đây là lúc ta khai thác nó, trước tiên Nhi muốn kiểm tra lại mô hình này trên testset, một cách độc lập, và tính những chỉ số như R2,rmse; mse…

library(mlr)

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

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

dumpred=predict(dummy,regr.task)
dumpred$data$response<-predict(m2,newdata=testM,type="response")

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

mets=list(rsq,mse,rmse)

performance(dumpred,measures =mets)

##       rsq       mse      rmse 
## 0.3998107 0.5183984 0.7199989

Tiếp theo, Nhi dùng hàm predictAll để ước tính cả 3 tham số L,M,S trong mô hình. Chỗ lắt léo duy nhất ở đây đó là Nu trong BCCG chính là Lambda trong method LMS:

lms=predictAll(m2,newdata=dfM,type="response")

## 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)

pdf=dfM%>%mutate(Mu=lms$mu,
                 LLN=exp(log(lms$mu)+(log(1-1.645*lms$nu*lms$sigma))/lms$nu),
                 ULN=exp(log(lms$mu)+(log(1+1.645*lms$nu*lms$sigma))/lms$nu),
                 Z_score=((TLC/lms$mu)^(lms$nu)-1)/(lms$nu*lms$sigma),
                 pcPred=100*TLC/lms$mu,
                 cutof80=0.8*lms$mu)

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

TLC	Age	Height	Mu	LLN	ULN	Z_score	pcPred	cutof80
7.017265	18	177.5	6.831976	5.922563	7.871907	0.3096632	102.71210	5.465580
7.648100	18	181.0	7.129516	6.180497	8.214737	0.8135273	107.27377	5.703613
6.737513	18	176.0	6.692085	5.801294	7.710723	0.0782444	100.67883	5.353668
5.949827	22	171.0	6.276478	5.419907	7.259549	-0.6007951	94.79563	5.021183
7.389043	22	180.0	7.128793	6.155904	8.245360	0.4040904	103.65068	5.703035
7.048908	23	185.0	7.491865	6.463038	8.673720	-0.6803729	94.08750	5.993492

Sau đó, Nhi tạo ra 2 biến mới là ToleranceZscore và Tolerance%pred, để kiểm tra liệu trong 243 người bình thường này, liệu mô hình có phân loại nhầm ai đó thành bệnh lý hay không ? Như vậy Z-score không được thấp hơn -1.645 còn %Pred không được thấp hơn 80%

pdf$ToleranceZscore=if_else(pdf$Z_score<(-1.645),
                           "Unacceptable", 
                           "Acceptable")

pdf$TolerancePcPred=if_else(pdf$pcPred<80,
                           "Unacceptable", 
                           "Acceptable")

Ta sẽ nhìn kết quả này một cách trực quan như sau:

pdf%>%ggplot(aes(x=Age))+
  geom_point(aes(y=TLC),color="grey80")+
  geom_smooth(aes(y=ULN),se=F,color="blue",linetype=2)+
  geom_smooth(aes(y=LLN),se=F,color="red",linetype=2)+
  geom_smooth(aes(y=Mu),se=F,color="black")+
  geom_smooth(aes(y=cutof80),se=F,color="orange")+
  geom_point(aes(y=Mu,color=ToleranceZscore),size=2,alpha=0.6)+
  scale_color_brewer(palette = "Set1",direction = -1)+
  theme_bw()+ggtitle("Mu,%pred,LLN and ULN for Males")

Trên hình vẽ: màu đen là Mu, Xanh là ULN, đỏ là LLN, màu cam là ngưỡng 80% của Mu

pdf%>%ggplot(aes(x=Height))+
  geom_point(aes(y=TLC),color="grey80")+
  geom_smooth(aes(y=ULN),se=F,color="blue",linetype=2)+
  geom_smooth(aes(y=LLN),se=F,color="red",linetype=2)+
  geom_smooth(aes(y=Mu),se=F,color="black")+
  geom_smooth(aes(y=cutof80),se=F,color="orange")+
  geom_point(aes(y=Mu,color=ToleranceZscore),size=2,alpha=0.6)+
  scale_color_brewer(palette = "Set1",direction = -1)+
  theme_bw()+ggtitle("Mean predicted,LLN and ULN for Males")

Ta có thể kết luận rằng mô hình đã thể hiện khá trung thành khuynh hướng biến thiên của TLC theo Age và Height. Chỉ có 1 số ít cá thể bị phân loại nhầm trên thực tế.

Khi chồng lắp 2 phân phối của Mu (màu cam) và của TLC thực tế (màu đỏ),ta thấy mô hình cho phép ước tính chính xác vị trí trung tâm của TLC trong mẫu, nhưng có một phần nào đó TLC đã bị đánh giá thấp hơn thực tế.

pdf%>%ggplot()+
  geom_density(aes(x=Mu),color="orange",linetype=1,fill="orange",alpha=0.5)+
  geom_density(aes(x=TLC),color="red",linetype=2,fill="red",alpha=0.2)+
  theme_bw()

Khi chuyển cả 2 thành thang đo Z-score, thì Z-score ước tính bởi mô hình hoàn toàn tương hợp với Z-score của TLC trên thực tế.

pdf%>%ggplot()+
  geom_density(aes(x=Z_score),color="red",linetype=1,fill="red",alpha=0.3)+
  geom_density(aes(x=scale(TLC)),color="blue",linetype=1,fill="blue",alpha=0.3)+
  theme_bw()+ggtitle("Predicted (Red) vs Truth (Blue)")

6 Bước 4: Ngoại suy

Cuối cùng, Nhi tò mò muốn biết liệu ta có thể ngoại suy mô hình này cho những người già > 85 tuổi và những thiếu niên < 18 tuổi hay không ? Để làm việc này Nhi tạo ra dữ liệu mô phỏng:

oldmendf=filter(dfM,Age==max(dfM$Age))
boydf=filter(dfM,Age==min(dfM$Age))

extdf=data_frame(Height=rnorm(oldmendf$Height-5,n=200,sd=10),
                 Age=rnorm(92,5,n=200))

expred=predictAll(m2,
                  newdata=extdf,
                  type="response")

## 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)

extdf=extdf%>%mutate(Mu=expred$mu,
                 LLN=exp(log(expred$mu)+
                           (log(1-1.645*expred$nu*expred$sigma))/expred$nu),
                 ULN=exp(log(expred$mu)+
                           (log(1+1.645*expred$nu*expred$sigma))/expred$nu),
                 cutof80=0.8*expred$mu)

youngdf=data_frame(Height=rnorm(mean(boydf$Height)-10,n=200,sd=10),
                 Age=rnorm(15,3,n=200))

ypred=predictAll(m2,
                  newdata=youngdf,
                  type="response")

## 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)

ydf=youngdf%>%mutate(Mu=ypred$mu,
                     LLN=exp(log(ypred$mu)+
                               (log(1-1.645*ypred$nu*ypred$sigma))/ypred$nu),
                     ULN=exp(log(ypred$mu)+
                               (log(1+1.645*ypred$nu*ypred$sigma))/ypred$nu),
                     cutof80=0.8*ypred$mu)

extdf=extdf%>%mutate(Mu=expred$mu,
                     LLN=exp(log(expred$mu)+
                               (log(1-1.645*expred$nu*expred$sigma))/expred$nu),
                     ULN=exp(log(expred$mu)+
                               (log(1+1.645*expred$nu*expred$sigma))/expred$nu),
                     cutof80=0.8*expred$mu)

extradf=rbind(extdf,ydf,
              dplyr::select(pdf,colnames(extdf)))

extradf=extradf%>%mutate(Extra=if_else(Age<18 | Age>85,"Extrapolation","True"))

extradf%>%ggplot(aes(x=Age))+
  geom_smooth(aes(y=ULN),se=F,color="blue",linetype=2)+
  geom_smooth(aes(y=LLN),se=F,color="red",linetype=2)+
  geom_smooth(aes(y=Mu),se=F,color="black")+
  geom_smooth(aes(y=cutof80),se=F,color="orange")+
  geom_point(aes(y=Mu,col=Extra,fill=Extra),shape=21,size=2,alpha=0.5,show.legend = F)+
  geom_vline(xintercept = c(18,85),col=c("blue4","red4"),linetype=1)+
  scale_y_continuous("TLC")+
  scale_color_brewer(palette = "Set1")+
  scale_fill_brewer(palette = "Set1")+
  theme_bw()+ggtitle("Extrapolation for <18 and >85")

Kết quả mô phỏng cho phép dự báo kết quả ngoại suy mô hình cho những đối tượng ngoài giới hạn Tuổi trong mẫu khảo sát. Đây chỉ là kết quả mô phỏng, với giả định chiều cao không biến đổi quá lớn và có phân phối chuẩn. Trên thực tế những thiếu niên từ 14-17 tuổi có thể tăng trưởng chiều cao nhiều hơn, tương tự người già > 85 tuổi có khuynh hướng giảm chiều cao.

7 Tổng kết

Bài thực hành đến đây là kết thúc. Nhi đã chuyển đến các bạn toàn bộ nguyên lý đằng sau mô hình LMS và những bí quyết để xây dựng nó trong gamlss. Những điều này chưa bao giờ được giải thích đầy đủ bởi các tổ chức như GLI và các nhóm nghiên cứu khác. Bây giờ thì họ không còn độc quyền nữa. Nếu có đủ dữ liệu tốt trong tay, các bạn hoàn toàn có thể làm được những mô hình tăng trưởng và mô hình ước tính ngưỡng giá trị chẩn đoán cho xét nghiệm sinh lý, sinh hóa cho người Việt.

Cảm ơn các bạn và hẹn gặp lại trong một tutorial khác.

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. Stanojevic S, Wade A, Stocks J, et al. Reference ranges for spirometry across all ages. A new approach. Am J Respir Crit Care Med 2008; 177: 253–260.

2. Quanjer PH, Stanojevic S, Cole TJ, et al. and the ERS Global Lung Function Initiative. Multi-ethnic reference values for spirometry for the 3-95 years age range: the Global Lung Function 2012 equations. Eur Respir J 2012; 40: 1324-1343.

3. Sandercock G., Voss C., Cohen D., Taylor M., and Stasinopoulos D. M. (2012) Centile curves and normative values for the twenty metre shuttle-run test in English schoolchildren. Journal of Sports Sciences, DOI:10.1080/02640414.2012.660185

4. WHO Multicentre Growth Reference Study Group (2007) WHO Child Growth Standards: Head circumference-for-age, arm circumference-for-age, triceps circumference-for-age and subscapular skinford-for-age: Methods and development. Geneva: World Health Organization. Available from: http://www.who.int/childgrowth/standards/second_set/technical_report_2/en/.

5. Rigby RA, Stasinopoulos DM. Smooth centile curves for skew and kurtotic data modelled using the Box-Cox power exponential distribution. Stat Med. 2004;23:3053-76.