Sử dụng Gamlss

Bài 1: Tại sao cần dùng Gamlss ?

1 Giới thiệu về Gamlss

Năm 2012, Nhóm nghiên cứu Philip H. Quanjer, Sanja Stanojevic, Tim J. Cole và Janet Stocks công bố mô hình GLI (Global Lung function initiative) cho phép ước tính giá trị bình thường cho xét nghiệm hô hấp ký của trẻ em và người lớn. Đó là mô hình phức tạp nhất tính đến thời điểm đó; Nhi rất ngạc nhiên khi nhìn thấy đồ thị của mô hình, vì nó thể hiện một cách hoàn hảo đường cong với nhiều đoạn uốn lượn mềm mại tương ứng với 3 giai đoạn : tăng trưởng nhanh của dung tích phổi từ 3-17 tuổi, ổn định 18-30 và sau đó bắt đầu thoái triển cho đến năm 80-90 tuổi.Để làm việc này, tác giả đã xác định phương trình cho cả 3 tham số Mu (vị trí trung bình), Sigma (coefficient of variation) và Lambda (skewness) của 1 phân phối Box-Cox Cole-Green cho biến kết quả, mỗi phương trình này là một mô hình đa thức bậc cao theo tuổi và chiều cao, kèm theo smoothing Spline. Đó là lần đầu tiên tôi nghe tới GAMLSS (khi đó tôi còn chưa biết R là gì).

Sau khi học R khoảng 1 năm và quen thuộc với hồi quy tuyến tính và mixed Model thì Nhi quyết định khám phá package Gamlss và nhận ra đây không phải là 1 package dễ học. Nhiều lần Nhi phải ngừng bước vì cấu trúc phức tạp bên trong gamlss và mãi đến năm 2016 mới luyện xong package này (mất 3 năm).

Có thể nói Gamlss là võ lâm chí tôn cho phái hồi quy. Nó là 1 framework phức tạp và tổng quát nhất cho mô hình hồi quy, với khả năng vô cùng mạnh mẽ và tinh tế như sau :

Gamlss hỗ trợ hơn 100 họ phân phối khác nhau, và hơn thế nữa, nó cho phép bạn tạo ra kiểu phân phối cho riêng mình do đó nó có khả năng biểu diễn tất cả quy luật phân phối nào dù kì dị đến đâu.

Gamlss cho phép mô hình hóa tất cả (từ 2-4) tham số của 1 họ phân phối bao gồm mô hình cho Mean (Mu, vị trí), Sigma (Scale), Tau và Nu (kiểu hình, shape). Thí dụ với phân phối Gaussian, bạn có thể dựng 1 mô hình cho Mean và 1 mô hình cho Sigma, với link function tùy chọn. Ngoài ra gamlss còn cho phép làm quantile regression.

Gamlss hỗ trợ tất cả các dạng Additive terms, bao gồm Smoothing với Penalized spline, P-splines, Cubic splines, Ridge và Lasso, GMRF, tensor product splines, thin plate splines; Loess spline, neural network, decision tree và MARS. Như vậy từ mô hình tuyến tính, bạn có thể biến hóa thành đủ loại mô hình phi tuyến tính, kết hợp hồi quy tuyến tính (cho linear terms) với Kernel, neural network và decision tree (cho smoothing term)

Gamlss cho phép dựng mô hình với Random Effects

Gamlss có cả tính năng Machine learning với K fold Cross-validation, bootstrap, Training và validation

2 Một thí dụ minh họa

Chúng ta sẽ sử dụng dataset về nồng độ glycosaminoglycans (GAG) ở trẻ em từ sơ sinh đến 18 tuổi. Dataset này có 314 quan sát. Mục tiêu của chúng ta là dựng một mô hình để khảo sát giá trị GAG theo tuổi.

df=read.csv("http://vincentarelbundock.github.io/Rdatasets/csv/MASS/GAGurine.csv")

age.cat <- function(x, lower = 0, upper, by = 10,
                    sep = "-", above.char = "+") {
  
  labs <- c(paste(seq(lower, upper - by, by = by),
                  seq(lower + by - 1, upper - 1, by = by),
                  sep = sep),
            paste(upper, above.char, sep = ""))
  
  cut(floor(x), breaks = c(seq(lower, upper, by = by), Inf),
      right = FALSE, labels = labs)
}


df$range=age.cat(x=df$Age,upper=17,by=2)

Thăm dò số liệu cho ta thấy quan hệ giữa GAG và Age không thể nào được mô tả bằng một mô hình tuyến tính, ngoài ra phương sai của GAG cũng không đồng nhất giữa các phân nhóm tuổi

library(tidyverse)

mycolors=c("#ffe314","#ffb814","#ff7d14","#ff2f14","#ff2c28","#ff3251","#ff1499","#e714ff","#b814ff")


df%>%ggplot(aes(x=range,y=GAG))+geom_boxplot(aes(fill=range))+
  scale_fill_manual(values=mycolors)+
  theme_bw(10)+
  scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40,45,50,55))

df%>%ggplot(aes(x=Age,y=GAG))+
  stat_density2d(geom="polygon",aes(fill=df$range,alpha = ..level..),show.legend = F)+
  scale_fill_manual(values=mycolors)+
  geom_point(shape=21,color="black",aes(fill=df$range),alpha=0.8,show.legend = F)+
  theme_bw(10)+
  scale_x_continuous(breaks=c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18))+
  scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40,45,50))

3 So sánh Cú pháp Gamlss với hàm lm() và glm()

Đầu tiên, chúng ta sẽ làm quen với gamlss bằng 1 mô hình hồi quy tuyến tính đơn biến hết sức đơn giản, và như ta thấy, cú pháp gamlss hoàn toàn tương đồng với những hàm lm() và glm() trong base R. Nội dung mô hình (hệ số hồi quy) có thể xem như tương đương, tuy nhiên như đã đề cập, Gamlss cung cấp 2 mô hình cho cả Mu và Sigma thay vì chỉ dự báo Mu ; do đó độ tự do chỉ còn 311, vì 1 df đã được dùng để tính sigma như 1 hằng số)

l1=lm(data=df,GAG~Age)
summary(l1)

## 
## Call:
## lm(formula = GAG ~ Age, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.950 -4.217 -1.596  2.477 36.470 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 19.89381    0.52553   37.85   <2e-16 ***
## Age         -1.27253    0.07242  -17.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.386 on 312 degrees of freedom
## Multiple R-squared:  0.4974, Adjusted R-squared:  0.4958 
## F-statistic: 308.7 on 1 and 312 DF,  p-value: < 2.2e-16

library(gamlss)

g1=gamlss(data=df,formula=GAG~Age,family=NO,trace=FALSE)

summary(g1)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = GAG ~ Age, family = NO, data = df,  
##     trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 19.89381    0.52385   37.98   <2e-16 ***
## Age         -1.27253    0.07219  -17.63   <2e-16 ***
## ---
## 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.8509     0.0399   46.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  314 
## Degrees of Freedom for the fit:  3
##       Residual Deg. of Freedom:  311 
##                       at cycle:  2 
##  
## Global Deviance:     2053.451 
##             AIC:     2059.451 
##             SBC:     2070.699 
## ******************************************************************

Ta có thể tính Rsq, CI95% cho hệ số hồi quy và vẽ 1 số đồ thị để kiểm tra các giả định về phẩm chất mô hình. Như đã dự tính, mô hình tuyến tính cho ra kết quả rất tồi.

confint(g1)

##                 2.5 %    97.5 %
## (Intercept) 18.867081 20.920533
## Age         -1.414018 -1.131032

Rsq(g1,type="both")

## $CoxSnell
## [1] 0.4973689
## 
## $CraggUhler
## [1] 0.4977304

plot(g1)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  9.037173e-13 
##                        variance   =  1.003195 
##                coef. of skewness  =  2.093835 
##                coef. of kurtosis  =  10.03964 
## Filliben correlation coefficient  =  0.9140833 
## ******************************************************************

Để các bạn yên tâm, ta sẽ làm thử 1 mô hình hồi quy với phân phối Gamma sử dụng hàm gamlss và glm

g2=gamlss(data=df,formula=GAG~Age,family=GA,trace=FALSE)

summary(g2)

## ******************************************************************
## Family:  c("GA", "Gamma") 
## 
## Call:  gamlss(formula = GAG ~ Age, family = GA, data = df,  
##     trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.016866   0.028753  104.92   <2e-16 ***
## Age         -0.110826   0.003827  -28.96   <2e-16 ***
## ---
## 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.01510    0.03906  -25.99   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  314 
## Degrees of Freedom for the fit:  3
##       Residual Deg. of Freedom:  311 
##                       at cycle:  2 
##  
## Global Deviance:     1752.212 
##             AIC:     1758.212 
##             SBC:     1769.46 
## ******************************************************************

Rsq(g2)

## [1] 0.7003374

plot(g2)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -0.003817336 
##                        variance   =  1.002205 
##                coef. of skewness  =  1.132927 
##                coef. of kurtosis  =  5.525881 
## Filliben correlation coefficient  =  0.9595208 
## ******************************************************************

l2=glm(data=df,GAG~Age,family=Gamma(link="log"))

summary(l2)

## 
## Call:
## glm(formula = GAG ~ Age, family = Gamma(link = "log"), data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2949  -0.2697  -0.1167   0.0969   1.5039  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.016891   0.034434   87.61   <2e-16 ***
## Age         -0.110831   0.004745  -23.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.1750716)
## 
##     Null deviance: 134.215  on 313  degrees of freedom
## Residual deviance:  42.132  on 312  degrees of freedom
## AIC: 1758.3
## 
## Number of Fisher Scoring iterations: 6

Ta thấy rằng Rsq đã có cải thiện đáng kể ngay sau khi thay Gaussian bằng Gamma, tuy nhiên đây chưa phải là kết quả tối ưu, như đồ thị QQ plot cho thấy

4 Hàm gamlss với các họ phân phối Log normal, Box-Cox T và Box-Cox Cole Green

Các bạn đã quen với cú pháp hàm gamlss rồi, bây giờ ta chơi 1 trò chơi nhỏ, đó là lần lượt thử 3 họ phân phối khác nhau là Logarithmic normal, Box-Cox T và Box Cox Cole Green.

Các bạn có thể nhận ra là với những họ phân phối phức tạp hơn, Gamlss bắt đầu dựng mô hình cho nhiều tham số hơn, thí dụ Sigma, Tau và Nu ; bạn có thể đặt formula cho từng tham số tùy thích. Đây là thứ mà glm() không thể làm được.

g3=gamlss(data=df,formula=GAG~Age,sigma.fo=~Age,family=BCT,trace=FALSE)

summary(g3)

## ******************************************************************
## Family:  c("BCT", "Box-Cox t") 
## 
## Call:  gamlss(formula = GAG ~ Age, sigma.formula = ~Age, family = BCT,  
##     data = df, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 15.23596    0.44151   34.51   <2e-16 ***
## Age         -0.76879    0.03159  -24.34   <2e-16 ***
## ---
## 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.025102   0.073839 -13.883   <2e-16 ***
## Age         -0.001339   0.010173  -0.132    0.895    
## ---
## 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.4833     0.1506  -3.209  0.00147 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.3199     0.4016   5.777 1.86e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  314 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  308 
##                       at cycle:  7 
##  
## Global Deviance:     1802.902 
##             AIC:     1814.902 
##             SBC:     1837.398 
## ******************************************************************

Rsq(g3,type ="both")

## $CoxSnell
## [1] 0.6364217
## 
## $CraggUhler
## [1] 0.6371652

plot(g3)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  0.004530074 
##                        variance   =  1.002391 
##                coef. of skewness  =  -0.04952054 
##                coef. of kurtosis  =  3.200192 
## Filliben correlation coefficient  =  0.9941185 
## ******************************************************************

g4=gamlss(data=df,formula=GAG~Age,family=LOGNO2(),trace=FALSE,method=RS(1000))

summary(g4)

## ******************************************************************
## Family:  c("LOGNO2", "Log Normal 2") 
## 
## Call:  gamlss(formula = GAG ~ Age, family = LOGNO2(), data = df,  
##     method = RS(1000), trace = FALSE) 
## 
## Fitting method: RS(1000) 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.888570   0.029988   96.33   <2e-16 ***
## Age         -0.104180   0.004117  -25.30   <2e-16 ***
## ---
## 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.03240    0.04089  -25.25   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  314 
## Degrees of Freedom for the fit:  3
##       Residual Deg. of Freedom:  311 
##                       at cycle:  135 
##  
## Global Deviance:     1727.618 
##             AIC:     1733.618 
##             SBC:     1744.866 
## ******************************************************************

Rsq(g4,type ="both")

## $CoxSnell
## [1] 0.7150416
## 
## $CraggUhler
## [1] 0.7158736

plot(g4)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  0.07329526 
##                        variance   =  0.9978055 
##                coef. of skewness  =  0.4928808 
##                coef. of kurtosis  =  4.648521 
## Filliben correlation coefficient  =  0.9784179 
## ******************************************************************

GAIC(g1,g2,g3,g4,k=0)

##    df      AIC
## g4  3 1727.618
## g2  3 1752.212
## g3  6 1802.902
## g1  3 2053.451

g5=gamlss(GAG~pb(Age),sigma.fo=~Age,nu.fo=~Age,family=BCCG,data=df,method=RS(1000))

## GAMLSS-RS iteration 1: Global Deviance = 1602.824 
## GAMLSS-RS iteration 2: Global Deviance = 1589.356 
## GAMLSS-RS iteration 3: Global Deviance = 1589.539 
## GAMLSS-RS iteration 4: Global Deviance = 1589.772 
## GAMLSS-RS iteration 5: Global Deviance = 1589.859 
## GAMLSS-RS iteration 6: Global Deviance = 1589.898 
## GAMLSS-RS iteration 7: Global Deviance = 1589.908 
## GAMLSS-RS iteration 8: Global Deviance = 1589.906 
## GAMLSS-RS iteration 9: Global Deviance = 1589.908 
## GAMLSS-RS iteration 10: Global Deviance = 1589.906 
## GAMLSS-RS iteration 11: Global Deviance = 1589.907

summary(g5)

## ******************************************************************
## Family:  c("BCCG", "Box-Cox-Cole-Green") 
## 
## Call:  gamlss(formula = GAG ~ pb(Age), sigma.formula = ~Age,  
##     nu.formula = ~Age, family = BCCG, data = df, method = RS(1000)) 
## 
## Fitting method: RS(1000) 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 16.46134    0.29701   55.42   <2e-16 ***
## pb(Age)     -0.94372    0.02686  -35.14   <2e-16 ***
## ---
## 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.36105    0.05672 -23.996   <2e-16 ***
## Age          0.02002    0.00789   2.538   0.0117 *  
## ---
## 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.05033    0.19041   0.264    0.792
## Age         -0.03250    0.02324  -1.399    0.163
## 
## ------------------------------------------------------------------
## 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:  314 
## Degrees of Freedom for the fit:  13.69461
##       Residual Deg. of Freedom:  300.3054 
##                       at cycle:  11 
##  
## Global Deviance:     1589.907 
##             AIC:     1617.296 
##             SBC:     1668.643 
## ******************************************************************

Rsq(g5,type ="both")

## $CoxSnell
## [1] 0.8154966
## 
## $CraggUhler
## [1] 0.8164492

plot(g5)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  0.00756865 
##                        variance   =  1.003139 
##                coef. of skewness  =  -0.03483846 
##                coef. of kurtosis  =  5.315206 
## Filliben correlation coefficient  =  0.984145 
## ******************************************************************

5 Sức mạnh của gamlss

Tiếp theo, Nhi sẽ biểu diễn sức mạnh của gamlss bằng 1 mô hình phức tạp hơn, gồm đa thức bậc 3 cho Mu, đa thức bậc 2 cho Sigma, kèm theo Spline bù trừ theo Age cho cả 2, với phân phối Gamma

Hàm GAIC( ) cho phép ta so sánh AIC của nhiều mô hình với nhau , kết quả cho thấy mô hình số 6 là tối ưu với AIC = 1611 và R2 = 0.82

g6=gamlss(GAG~poly(Age,3)+pb(Age),sigma.fo=~poly(Age,2)+pb(Age),family=GA,data=df,method=RS(1000))

## GAMLSS-RS iteration 1: Global Deviance = 1597.256 
## GAMLSS-RS iteration 2: Global Deviance = 1597.094 
## GAMLSS-RS iteration 3: Global Deviance = 1597.094

summary(g6)

## ******************************************************************
## Family:  c("GA", "Gamma") 
## 
## Call:  
## gamlss(formula = GAG ~ poly(Age, 3) + pb(Age), sigma.formula = ~poly(Age,  
##     2) + pb(Age), family = GA, data = df, method = RS(1000)) 
## 
## Fitting method: RS(1000) 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.41000    0.01702 141.630  < 2e-16 ***
## poly(Age, 3)1 -9.70609    0.40250 -24.115  < 2e-16 ***
## poly(Age, 3)2  2.82411    0.37473   7.536 5.47e-13 ***
## poly(Age, 3)3 -2.16491    0.33848  -6.396 5.94e-10 ***
## ---
## 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.25573    0.03932 -31.936  < 2e-16 ***
## poly(Age, 2)1  2.56608    0.69098   3.714 0.000243 ***
## poly(Age, 2)2  2.64892    0.68115   3.889 0.000123 ***
## ---
## 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 maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  314 
## Degrees of Freedom for the fit:  7.005903
##       Residual Deg. of Freedom:  306.9941 
##                       at cycle:  3 
##  
## Global Deviance:     1597.094 
##             AIC:     1611.106 
##             SBC:     1637.374 
## ******************************************************************

Rsq(g6,type ="both")

## $CoxSnell
## [1] 0.8171525
## 
## $CraggUhler
## [1] 0.8180771

plot(g6)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -0.001029199 
##                        variance   =  1.003389 
##                coef. of skewness  =  0.4691101 
##                coef. of kurtosis  =  4.942875 
## Filliben correlation coefficient  =  0.9775264 
## ******************************************************************

wp(object=g6,xvar=df$Age)

## number of missing points from plot= 0  out of  78

## number of missing points from plot= 0  out of  79

## number of missing points from plot= 1  out of  79

## number of missing points from plot= 0  out of  78

GAIC(g1,g2,g3,g4,g5,g6)

##           df      AIC
## g6  7.005903 1611.106
## g5 13.694613 1617.296
## g4  3.000000 1733.618
## g2  3.000000 1758.212
## g3  6.000000 1814.902
## g1  3.000000 2059.451

So Sánh trực quan giữa 3 mô hình số 2, số 5 và số 6 :

plot(GAG~ Age, data = df, cex.lab=1.5,pch=16,col=grey(0.7))
lines(df$Age[order(df$Age)], fitted(g2)[order(df$Age)],
      lty=1, lwd=2.5, col = "orange")

lines(df$Age[order(df$Age)], fitted(g5)[order(df$Age)],
      lty=1, lwd=2.5, col = "red3")

lines(df$Age[order(df$Age)], fitted(g6)[order(df$Age)],
      lty=1, lwd=2.5, col = "darkviolet")

legend("topright",legend=c("Gamma","BCCG+Spline","GA-Cubic-PS"),
       col=c("orange","red3","violet"), lwd=2, title="3 Models")

Đây là đồ thị của mô hình số 6

df=df%>%mutate(.,predict=predict(g6,type="response"))

df%>%ggplot()+
  stat_density2d(geom="polygon",aes(x=Age,y=GAG,fill=df$range,alpha = ..level..),show.legend = F)+
  scale_fill_manual(values=mycolors)+
  geom_point(shape=21,color="black",aes(x=Age,y=GAG,fill=df$range),alpha=0.2,show.legend = F)+
  geom_smooth(aes(x=Age,y=predict),se=F,color="black",linetype=2,show.legend = F)+
  theme_bw(10)+
  scale_x_continuous(breaks=c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18))+
  scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40,45,50))

Gamlss còn cho phép ta vẽ 1 đồ thị cho bách phân vị của giá trị dự báo :

centiles.fan(g6, xvar=df$Age, cent=c(2.5,5,10,25,50,75,90,95,97.5),colors = "heat",points=T,col=gray(0.3))

Hoặc 1 đồ thị tuyệt đẹp khác như dưới đây

library(gamlss.util)
library(viridisLite)

plotSimpleGamlss(GAG,Age,model=g6, data=df,cols=viridis(option="A",n=100),x.val=seq(1,18,2),val=6,N=500,ylim=c(0,50))

## new prediction 
## new prediction

6 Kết luận

Bài đầu tiên này Nhi chỉ muốn « khoe » với các bạn một số tính năng mà gamlss có thể làm được, nhằm giải thích tại sao chúng ta cần phải học về Gamlss và sử dụng nó. Đây là bài đầu tiên trong một series mà nhóm MLM sẽ tiến hành trong vài tháng tới, với mục tiêu hướng dẫn sử dụng thành thạo package này. Chúng tôi sẽ giải thích cú pháp của các hàm trong gamlss, những quy trình, các họ phân phối, các loại smothing splines… thông qua case study.

Sau khi luyện đến cảnh giới của package này, các bạn có thể hoàn toàn kiểm soát các quy luật phi tuyến tính thí dụ mô hình tăng trưởng cho chỉ số nhân trắc học, sinh lý, sinh hóa theo tuổi một cách tối ưu nhất.

Xin cảm ơn và hẹn gặp lại.