Sử dụng package Gamlss

Bài 2: Xác định quy luật phân phối cho mô hình GAM

1 Gamlss là một thư viện khổng lồ về quy luật phân phối

Trong bài trước, Nhi đã giới thiệu với các bạn package Gamlss, một bảo bối cho mô hình tuyến tính tổng quát. Một trong những ưu thế của gamlss là cấu trúc tinh tế của mô hình. Trước hết, GAM là những mô hình hồi quy bán tham số chuyên biệt về phân phối (distribution based semiparametric regression). Có nghĩa là trước tiên người dùng phải xác định 1 họ phân phối chuyên biệt cho biến kết quả ; sau đó tùy vào bản chất của phân phối này,từ 2 đến 4 mô hình khác nhau sẽ được xây dựng cho toàn bộ các tham số của phân phối.

Bốn tham số này được gọi tên là Mu, Sigma, Nu và Tau, tương ứng với 4 đặc tính : Vị trí trung tâm (Location = Mu), Scale (sai biệt, Sigma), Skewness (Nu) và Kurtosis (Tau). Như vậy Mu và Sigma xác định giá trị Trung bìhn và độ biến thiên của phân phối, trong khi Tau và Nu xác định kiểu hình của phân phối (hình dạng của đồ thị hàm pdf).

Thành tố thứ nhất của mô hình gamlss chính là họ phân phối và các tham số này. Thành tố thứ hai là các phần tử hồi quy (predictors) bao gồm hằng số (intercept), các biến độc lập, bậc đa thức và tương tác giữa chúng. Thành tố thứ 3 là phần tử bù trừ bao gồm tất cả các phương pháp hiệu chỉnh nhằm tăng khả năng (capacity) của mô hình, bao gồm smoothing splines, penalized splines, neural network, decision tree, piece wise…

Trong bài này chúng ta chỉ quan tâm đến Thành tố thứ nhất là Phân phối và tham số của nó.

Thống kê mô tả là một quy trình quen thuộc ta hay dùng để thăm dò dữ liệu, và mục tiêu của nó chính là nắm các đặc tính phân phối của dữ liệu trong mẫu. Nhi và các bạn vẫn hay tính giá trị trung bình, SD, skewness và kurtosis. Ta còn vẽ histogram để phân tích trực quan đặc tính này. Việc tính toán này mang lại thông tin khá đầy đủ về đặc tính phân phối của biến số. Giáo trình xác suất trong trường đại học đã mô tả cho chúng ta một số quy luật phân phối cổ điển như Gaussian, Poisson, Binomial, t và Chi2. Nhưng, để giúp các bạn liên hệ giữa biến số và mô hình GAMLSS, Nhi muốn chúng ta hình dung về biến số như 1 hàm toán học chứ không phải như 1 đại lượng. Thật vậy, bất cứ biến số nào cũng có thể xem như một hàm. Kết quả của hàm này là một con số (giá trị) và mỗi giá trị như vậy gắn với một xác suất mà giá trị đó có thể xuất hiện trong quần thể. Thí dụ, theo lý thuyết số lượng bạch cầu trong một mL máu tĩnh mạch có thể nhận bất cứ giá trị số nguyên dương nào trong khoảng từ 0 đến 10.000 nhưng xác suất quan sát được mỗi giá trị có thể khác nhau.

Hàm mật độ xác suất (probability density function – pdf ) cho phép ước tính xác suất cho mỗi giá trị bất kì trong 1 khoảng (thang đo) xác định. Hàm pdf có liên hệ rất gần gũi với biểu đồ tần suất (histogram). Trong thực hành, việc vẽ histogram kết hợp với đồ thị của hàm pdf (density curve) cho phép mô tả trực quan đặc tính phân phối của một biến số mà ta cần thăm dò. Từ histogram có thể ước tính hàm pdf nhưng bản chất của 2 thứ này hoàn toàn khác nhau : histogram là 1 phép thống kê (dựa vào chọn mẫu) trong khi PDF thuần túy là xác suất.

Ngoài ra, hàm mật độ xác suất tích lũy (cumulative density function, cdf) là một cách khác để biểu diễn đặc tính phân phối của một biến. CDF biểu thị tổng mật độ xác suất của tất cả trường hợp Y tính đến vị trí Yi mà ta cần khảo sát. Thực ra cả pdf và cdf đều cung cấp thông tin như nhau, và có thể hoán chuyển .

Như Nhi đã nói ở trên, thông tin quan trọng nhất ta thường khai thác đó là vị trí trung tâm của dữ liệu nằm ở đâu, và sai biệt bao lớn (cách vị trí trung tâm bao xa) ? Vị trí trung tâm được xác định bằng Mean, Median, Mode…, sai biệt được ước tính bằng độ lệch chuẩn (sigma, sd) hay phương sai (variance, sigma^2). Ta thường hiểu Mean như là arithmetric mean (cho số thực hay số nguyên) hoặc geometric mean (cho tỉ lệ), và ta cũng quen báo cáo Mean +/- SD với giả định là biến số có phân phối chuẩn. Tuy nhiên vị trí trung tâm thực ra còn có thể xem (một cách tổng quát) như một mô hình tối giản có phương trình là Y ~ 1 (chỉ chứa hằng số Intercept) với link function = identity và họ phân phối tùy chọn. Lúc này thực chất ta ước tính Mu (trung bình, giá trị kì vọng) cho một phân phối dựa trên dữ liệu quan sát được. Ngoài Mean và Sd, ta còn quan tâm đến Skewness và Kurtosis, Skewness cho biết số liệu của chúng ta có hình dạng đối xứng hay bị lệch về 1 phía, Kurtosis cho biết hình dạng (đỉnh) của phân phối. Một cách tổng quát, Gamlss cung cấp 1 framework để mô hình hóa (ước tính giá trị của tham số đích tùy theo 1 hay nhiều hiệp biến số/yếu tố khác) cho cả Mu, sigma và tham số kiểu hình (Nu, Tau) của một phân phối bất kì.

Trước khi thực hành mô hình GAMLSS các bạn cần biết 1 sự thật gây shock, đó là : Trong thế giới thực, phân phối chuẩn không tồn tại (suy rộng ra là dữ liệu thực tế không bao giờ phù hợp tuyệt đối với bất cứ một dạng phân phối lý thuyết nào mà trường lớp dạy cho các bạn, bao gồm phân phối chuẩn). Do đó, thay vì bỏ công làm những kiểm định để kiểm chứng dữ liệu của mình có phân phối « bình thường » hay không, ta cần tìm một quy luật phân phối gần (phù hợp) nhất với dữ liệu đang có để dựng mô hình.

Bây giờ bắt đầu phần hay nhất của câu chuyện hôm nay, đó là Gamlss hỗ trợ tất cả (xin nhấn mạnh, tất cả) họ phân phối thông qua module gamlss.dist ; điều này có nghĩa là gamlss là một thư viện khổng lồ về lý thuyết xác suất. Bạn có thể bước vào, chọn bất kì phân phối nào trong danh mục hàng trăm họ phân phối khác nhau cho số thực, số nguyên dương, tỉ lệ, … và chơi đùa thỏa thích với hàm pdf, hàm cdf, mô phỏng với hàm random, ước tính quantiles, vặn vẹo density curve và histogram với 2-4 tham số,bạn có thể chặn 1 hay 2 đầu, thay đổi function link từ identity sang logarithm, logit, probit… vân vân, có thể trộn phân phối số nguyên và số thực với nhau, và tạo ra một quy luật phân phối của riêng mình.

Sau khi load package gamlss, theo mặc định packae gamlss.dist cũng được gọi lên đồng thời.

library(tidyverse)

library(gamlss)

nC<-detectCores()

Trong thí dụ minh họa sau, Nhi dùng hàm dNBI và pNBI trong gamlss.dist để mô phỏng hàm PDF và CDF của 1 phân phối Negative binomial Type I (NBI) có Mu=10, sigma=0.5

plot(function(y) dNBI(y, mu = 10, sigma =0.5 ), from=0, to=50,n=50+1)%>%as_tibble()%>%ggplot(aes(x=x,y=y))+geom_area(alpha=0.5,fill="red",color="red")+theme_bw()+ggtitle("pdf of Negative Binomial type I;Mu=10;sigma=0.5")

cdf=stepfun(0:49, c(0, pNBI(0:49, mu=10, sigma=0.5 )),f=0)
p=plot(cdf,do.points=FALSE)

p$t=p$t[-1]

p=p%>%as_tibble()
p2 <- bind_rows(old = p,new = p %>% mutate(y = lag(y)),.id = "source")%>%arrange(t, source)

ggplot(p) + 
  geom_ribbon(aes(x =t, ymin = 0, ymax = y),data=p2,alpha=0.5,fill="red",color="red")+
  theme_bw()+scale_x_continuous("X")+
  scale_y_continuous("f(X)")+
  ggtitle("cdf of Negative Binomial type I;Mu=10;sigma=0.5")

Nhưng quan trọng nhất, đó là bạn có thể mang áp dụng bất kì họ phân phối nào trong danh mục cho mô hình của mình.

Cũng như trong phương pháp GLM cổ điển, tất cả 3 thành tố trong mô hình Gamlss đều có thể được lựa chọn, tinh chỉnh và tối ưu hóa nhờ vào quá trình chọn lọc, thử nghiệm, so sánh. Thí dụ, thành tố thứ nhất : phân phối (và tham số) có thể được chọn lọc ở nhiều cấp độ :

1. Trước công đoạn mô hình hóa : họ phân phối phù hợp được chỉ định dựa vào suy luận (bản chất biến số) và trực giác sau khi thăm dò số liệu (thống kê mô tả, histogram, density curve).

2. Quá trình thăm dò dữ liệu bằng thống kê mô tả có thể xem như tương đương với việc dựng 1 Mô hình tối giản : Thử nghiệm hàng loạt họ phân phối cho mô hình chỉ chứa intercept

3. Kết hợp với quy trình xây dựng mô hình : Model training (fitting), sử dụng các phương thức như Step-wise AIC, Likelihood ratio, K-fold cross validation, bootstrapping … nhằm chọn lọc mô hình chứa những biến số và họ phân phối phù hợp nhất.

4. Trong quá trình kiểm định mô hình : Khảo sát residual (sai số dự báo) của mô hình để xác nhận tính chính xác của họ phân phối đã chọn, so sánh hiệu năng mô hình (benchmark study)… Gamlss cung cấp gần như đầy đủ khả năng để làm tất cả 4 công đoạn nêu trên, kể cả K-folds cross validation (gamlss là 1 framework machine learning thật sự dành riêng cho algorithm GAM)

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

Nhi sẽ lấy 1 dataset đơn giản làm thí dụ. Đây là bộ số liệu về số lượng bạch cầu CD4 đo được trên 600 trẻ em sinh ra từ người mẹ bị nhiễm HIV. Mục tiêu của thí dụ là xây dựng mô hình dự báo giá trị lượng bạch cầu CD4 theo tuổi trong giới hạn từ sơ sinh đến 8 tuổi.

data(CD4)

psych::describe(CD4)

##     vars   n   mean     sd median trimmed    mad  min     max   range skew
## cd4    1 609 557.53 462.50 435.00  489.91 367.68 0.00 2327.00 2327.00 1.26
## age    2 609   2.62   1.31   2.35    2.49   1.20 0.25    8.18    7.93 1.08
##     kurtosis    se
## cd4     1.13 18.74
## age     1.55  0.05

Thăm dò số liệu cho ta thấy Biến cd4 có trung bình = 557.53, giá trị nhỏ nhất là 0, lớn nhất là 2327 (đơn vị tế bào/ml máu), skewness = 1.26, kurtosis = 1.13 và sd=462.5 Những giá trị này gợi ý về một phân phối điển hình của biến số đếm (count data), thực chất số lượng tế bào là 1 biến kiểu số và rời rạc chỉ nhận giá trị số nguyên dương.

3 Nhận định trực quan về đặc tính phân phối

Quy trình thăm dò, mô tả và thử nghiệm có thể được thực thi thông qua 2 hàm trong gamlss : truehist, histDist

Hàm truehist chỉ đơn giản vẽ ra 1 histogram, t có cũng thể làm việc này dễ dàng với base R hay cầu kì hơn với ggplot2

truehist(CD4$cd4)

hist(CD4$cd4)

CD4%>%ggplot(aes(x=cd4))+
  geom_density(alpha=.8,fill="#c507ff")+
  geom_histogram(aes(y=..density..,fill=..density..),colour="black",alpha=0.7,show.legend = F)+
  theme_bw()+scale_fill_gradient(low="#ff9c07",high="#ff072f")+
  ggtitle("Density + Histogram")

cdfdf=cbind.data.frame(cd4=unique(CD4$cd4),ecdf = ecdf(CD4$cd4)(unique(CD4$cd4)))

ggplot(cdfdf) + 
  geom_area(aes(x=cd4, y=ecdf),color="red",fill="#ff0741",alpha=0.8)+
  theme_bw()+scale_x_continuous("CD4 count")+
  scale_y_continuous("ecdf")+
  ggtitle("ECDF of cd4")

CD4%>%ggplot(aes(x=age,y=cd4))+
  geom_point(shape=21,aes(color=age,fill=age),show.legend = F)+
  geom_smooth(se=F,color="darkviolet")+
  theme_bw()+
  scale_fill_gradient2(high="#8e00b2",mid="#ff9c07",low="#ef9b09")+
  scale_color_gradient2(high="#8e00b2",mid="#ff9c07",low="#ef9b09")

4 Thăm dò phân phối dựa vào mô hình tối giản

Hàm histDist cho phép thăm dò 1 họ phân phối bất kì mà gamlss.dist hỗ trợ, thí dụ Negative binomial, bản chất của nó chính là một mô hình tối giản chỉ chứa Intercept với 1 họ phân phối xác định. Mô hình này cũng được lượng giá bằng BIC, AIC, Likelihood ratio… để dựa vào đó ta có thể so sánh các họ phân phối với nhau.

histDist(data=CD4,cd4,family=NBI())

## 
## Family:  c("NBI", "Negative Binomial type I") 
## Fitting method: "nlminb" 
## 
## Call:  gamlssML(formula = cd4, family = "NBI", data = CD4) 
## 
## Mu Coefficients:
## [1]  6.324
## Sigma Coefficients:
## [1]  -0.1869
## 
##  Degrees of Freedom for the fit: 2 Residual Deg. of Freedom   607 
## Global Deviance:     8909.14 
##             AIC:     8913.14 
##             SBC:     8921.96

fitDist là 1 hàm “thần kì”, khi áp dụng nó cho 1 biến số Y ta chỉ cần cung cấp thông tin về bản chất của thang đo (thí dụ số thực, số thực mở rộng, hỗn hợp, số đếm, xác suất…), gamlss sẽ sử dụng brute-force để kiểm tra toàn bộ (từ 20-40 loại) những họ phân phối tiềm năng cho biến số này, và tự động xác định phân phối tối ưu nhất.

Trong thí dụ này, hàm fitDist cho biết phân phối thích hợp nhất cho CD4 là Zero inflated negative binomial (ZINBI). Lưu ý là tiêu chuẩn lựa chọn ở đây là BIC (AIC hiệu chỉnh với k=log(cỡ mẫu))

Object xuất ra là 1 mô hình tối giản với phân phối ZINBI, ta có thể khảo sát nó : Ta có thể xâm nhập vào object để xem xét BIC của 16 họ phân phối mà mô hình có thể converged, ZINBI đứng đầu danh sách còn Poisson xếp chót bảng

set.seed=123

explore=fitDist(cd4,data=CD4,type="counts",k=log(nrow(CD4)),parallel="multicore",ncpus = nC)

explore

## 
## Family:  c("ZINBI", "Zero inflated negative binomial type I") 
## Fitting method: "nlminb" 
## 
## Call:  gamlssML(formula = y, family = DIST[i], parallel = "multicore",  
##     ncpus = ..2, data = sys.parent()) 
## 
## Mu Coefficients:
## [1]  6.335
## Sigma Coefficients:
## [1]  -0.3108
## Nu Coefficients:
## [1]  -4.478
## 
##  Degrees of Freedom for the fit: 3 Residual Deg. of Freedom   606 
## Global Deviance:     8875.95 
##             AIC:     8881.95 
##             SBC:     8895.19

plot(explore)

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  0.008977322 
##                        variance   =  0.9817547 
##                coef. of skewness  =  -0.1344454 
##                coef. of kurtosis  =  3.020631 
## Filliben correlation coefficient  =  0.9966616 
## ******************************************************************

explore$fits

##      ZINBI      ZANBI       NBII        NBI       GEOM     SICHEL 
##   8895.188   8895.188   8921.963   8921.963   8927.554   8928.375 
##        DEL     WARING      ZIPIG        PIG       ZALG       YULE 
##   8928.375   8933.966   9222.426   9454.245   9931.132  14307.884 
##       ZIP2        ZAP        ZIP         PO 
## 213819.609 213819.609 213819.609 221587.433

Có 1 package khác trong R là fitdistrplus cũng cho phép làm điều tương tự, chỉ khác là nó không phổ quát như gamlss.dist, vì số họ phân phối mà hàm descdist hỗ trợ chỉ đếm trên đầu ngón tay. Ta thử khảo sát biến cd4 bằng package fitdistrplus như sau:

Dựa vào 1 bootstrap 10 ngàn lượt và Biểu đồ Cullen&Frey , ta có thể khẳng định phân phối của cd4 KHÔNG thể là Poisson, mà cũng chẳng phải là Negative binomial và chắc chắn không thể là normal

library(fitdistrplus)

## Warning: package 'fitdistrplus' was built under R version 3.4.1

## Loading required package: survival

descdist(CD4$cd4,discrete = TRUE,boot=10000,boot.col="gold",obs.pch=16,obs.col="red4")

## summary statistics
## ------
## min:  0   max:  2327 
## median:  435 
## mean:  557.5337 
## estimated sd:  462.4983 
## estimated skewness:  1.263977 
## estimated kurtosis:  4.163153

Gamlss hỗ trợ 24 họ phân phối cho số đếm (count data) , đặc tính của chúng nằm trong bảng sau đây: Phân phối ZINBI có 3 tham số Mu (link function= log), Sigma (link function = log) và Nu (link function = logit).

Ta có thể tìm hiểu về bản chất của họ phân phối bằng hàm gamlss.family và show.link

gamlss.family(ZINBI)

## 
## GAMLSS Family: ZINBI Zero inflated negative binomial type I 
## Link function for mu   : log 
## Link function for sigma: log 
## Link function for nu   : logit

show.link(ZINBI)

## $mu
## c("inverse", "log", "identity")
## 
## $sigma
## c("inverse", "log", "identity")
## 
## $nu
## c("logit", "probit", "cloglog", "log", "own")

Module gamlss.tr cho phép người dùng tạo ra 1 phân phối của riêng mình bằng cách chặn 2 đầu 1 phân phối cho trước, thí dụ ta có thể tạo ra 1 phân phối mới có tên là ZINBI1000 với giá trị dao động trong khoảng từ 0-1000

Từ lúc này, họ phân phối mới có thể được đưa vào hàm fitDist và cả hàm gamlss để sử dụng như bất cứ họ phân phối nào được gamlss chính thức hỗ trợ.

library(gamlss.tr)
gen.trun(par=c(0,1000),family="ZINBI", name="1000", type="right",varying=T)

## A truncated family of distributions from ZINBI has been generated 
##  and saved under the names:  
##  dZINBI1000 pZINBI1000 qZINBI1000 rZINBI1000 ZINBI1000 
## The type of truncation is right 
##  and the truncation parameter is 0 1000

gamlss.family(ZINBI1000)

## 
## GAMLSS Family: ZINBI1000 right truncated Zero inflated negative binomial type I 
## Link function for mu   : log 
## Link function for sigma: log 
## Link function for nu   : logit

fitDist(cd4,data=CD4,type="counts",k=log(nrow(CD4)),extra="ZINBI1000",parallel="multicore",ncpus = nC)

## 
## Family:  c("ZINBI", "Zero inflated negative binomial type I") 
## Fitting method: "nlminb" 
## 
## Call:  gamlssML(formula = y, family = DIST[i], parallel = "multicore",  
##     ncpus = ..2, data = sys.parent()) 
## 
## Mu Coefficients:
## [1]  6.335
## Sigma Coefficients:
## [1]  -0.3108
## Nu Coefficients:
## [1]  -4.478
## 
##  Degrees of Freedom for the fit: 3 Residual Deg. of Freedom   606 
## Global Deviance:     8875.95 
##             AIC:     8881.95 
##             SBC:     8895.19

5 Chọn lọc họ phân phối trong quá trình xây dựng mô hình: Vuong test và stepwise GAIC

Bây giờ Nhi sẽ thử cách làm khác, đó là kiểm tra họ phân phối TRONG KHI dựng mô hình, tức là mô hình sẽ chứa đầy đủ predictor chứ không chỉ có Intercept. Cách làm thủ công đơn giản nhất đó là so sánh bắt cặp 2 mô hình với 2 họ phân phối khác nhau (thí dụ Poisson và ZINBI), sau đó sử dụng 2 kiểm định Vương-Clarke để so sánh độ phù hợp dữ liệu của 2 mô hình (ta không thể dùng Likelihood ratio test thông thường do 2 mô hình này non-nested)

mPO=gamlss(formula=cd4~pb(age),data=CD4,family=PO(),
          parallel="multicore",
          ncpus = nC)

mZINBI=gamlss(formula=cd4~pb(age),data=CD4,family=ZINBI(),
              parallel="multicore",
              ncpus = nC)

Kết quả của Vuong test cho ra chỉ số thống kê Z rất thấp, còn Clarke test thì cho ra giá trị p rất thấp, cho thấy mô hình ZINBI tốt hơn so với Poisson

VC.test(mPO, mZINBI, sig.lev = 0.05)

##  Vuong's test: -16.973 model mZINBI is preferred over mPO 
## Clarke's test: 83 p-value= 0 mZINBI is preferred over mPO

Ta cũng có thể dùng tiêu chuẩn AIC và BIC để lựa chọn mô hình : Kết quả của cả BIC và AIC đều ủng hộ cho mô hình ZINBI

GAIC(mPO,mZINBI)

##              df        AIC
## mZINBI 6.657584   8672.814
## mPO    6.951758 147256.463

GAIC(mPO,mZINBI,k=log(nrow(CD4)))

##              df        AIC
## mZINBI 6.657584   8702.186
## mPO    6.951758 147287.132

Ta có thể thực hiện quy trình này cho hàng loạt biến số nhưng sử dụng trực tiếp tiêu chí AIC hoặc BIC để quyết định lựa chọn phân phối nào. Trong thí dụ sau Nhi so sánh 4 mô hình với phân phối Poisson, Negative binomial, Zero inflated poisson, và zero inflated negative binomial.

#GAIC

mPO=gamlss(cd4~pb(age),
            data=CD4,
            family=PO(), 
            rand=rand,
            parallel="multicore",
            ncpus = nC)

mZINBI=gamlss(cd4~pb(age),
            data=CD4,
            family=ZINBI(), 
            rand=rand,
            parallel="multicore",
            ncpus = nC)

mNBI=gamlss(cd4~pb(age),
          data=CD4,
          family=NBII(), 
          rand=rand,
          parallel="multicore",
          ncpus = nC)

mZIP=gamlss(cd4~pb(age),
            data=CD4,
            family=ZIP(), 
            rand=rand,
            parallel="multicore",
            ncpus = nC,method=RS(500))

GAIC(mPO,mZIP,mZINBI,mNBI,k=log(nrow(CD4)))

##               df        AIC
## mZINBI  6.657584   8702.186
## mNBI   10.950083   8834.964
## mZIP   23.446697 133872.144
## mPO     6.951758 147287.132

6 Chọn lọc quy luật phân phối dựa vào K folds cross validation

Kiểm chứng chéo lặp lại là 1 cách làm khác. Ngoài quy trình stepwise, gamlss còn cho phép dùng K fold cross validation để so sánh độ phù hợp của hàng loạt họ phân phối, thí dụ Poisson, Negative binomial, Zero inflated poisson, zero inflated negative binomial. Chúng ta sẽ thực hiện 1 quy trình như vậy với 10 folds cross validation, chia ngẫu nhiên dataset CD4 ra thành 10 blocks, dùng 9 blocks để dựng mô hình rồi kiểm định trên 1 block còn lại, cứ thế lặp lại 10 lần.

set.seed(123)
rand=sample (10, nrow(CD4), replace=TRUE)

mcvPO=gamlssCV(cd4~pb(age),
               data=CD4,
               family=PO(), 
               rand=rand,
               parallel="multicore",
               ncpus = nC)

mcvNBI=gamlssCV(cd4~pb(age),
            data=CD4,
            family=NBII(), 
            rand=rand,
            parallel="multicore",
            ncpus = nC)

mcvZINBI=gamlssCV(cd4~pb(age),
            data=CD4,
            family=ZINBI(), 
            rand=rand,
            parallel="multicore",
            ncpus = nC)

mcvZIP=gamlssCV(cd4~pb(age),
            data=CD4,
            family=ZIP(), 
            rand=rand,
            parallel="multicore",
            ncpus = nC)

CV(mcvZIP,mcvZINBI,mcvNBI,mcvPO)

##          val[o.val]
## mcvZINBI   8671.850
## mcvNBI     8813.045
## mcvZIP   142360.871
## mcvPO    150765.188

7 Chọn lọc biến số cho mô hình của từng tham số, sử dụng StepwiseGAIC

Một khi đã chắc chắn về họ phân phối, ta có thể đi bước tiếp theo, đó là xác định xem liệu có cần mô hình hóa cho tất cả 3 tham số không ? và mô hình phức tạp (sâu) đến mức nào ? Hay nói cách khác: Ta sẽ phân phối predictor vào mỗi tham số Mu, Sigma, Nu như thế nào cho hợp lý ? Vì theo mặc định, nếu không khai báo formula cho Sigma và Nu thì gamlss ưu tiên dựng mô hình cho Mu (trung bình dự báo), còn sigma, Nu và Tau chỉ được mô hình bằng hằng số Intercept. Câu hỏi này có thể được giải quyết bằng hàm stepGAICAll.A; hàm này sẽ thực hiện 1 quy trình lựa chọn biến số cho từng tham số Mu, Sigma, Nu (và Tau, nếu cần) theo quy tắc stepwise, tức là loại bỏ dần dần biến số, dựa vào tiêu chuẩn AIC hoặc BIC. Trong thí dụ dưới đây, chúng ta đặt ra 2 ngưỡng: thấp nhất là 1 mô hình chỉ chứa Intercept, cao nhất là mô hình đa thức bậc 5 cho biến số age kèm theo penalize B spline cho age, tiêu chuẩn lựa chọn là BIC.

Đặc biệt, gamlss cho phép thực hiện quy trình song song và tận dụng tối đa là 4 cores trong máy tính để tăng tốc độ xử lý (tăng đến 50% cho 1 máy tính Core i5)

m1<-gamlss(cd4~1, data=CD4, family=ZINBI(), trace=FALSE)
m2<-stepGAICAll.A(m1, 
                  scope=list(lower=~1, 
                             upper=~poly(age,5)+pb(age)
                             ),
                  k=log(nrow(CD4)),
                  parallel="multicore",
                  ncpus = nC
                  )

Kết quả không bất ngờ lắm: quy trình stepwiseAIC cho thấy chỉ cần pbspline cho Mu, và Intercept cho 2 tham số Sigma và Nu.

summary(m2)

## ******************************************************************
## Family:  c("ZINBI", "Zero inflated negative binomial type I") 
## 
## Call:  gamlss(formula = cd4 ~ pb(age), family = ZINBI(), data = CD4,  
##     trace = FALSE, sigma.formula = ~1, nu.formula = ~1) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.67254    0.06711  114.32   <2e-16 ***
## pb(age)     -0.43241    0.02284  -18.93   <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) -0.62324    0.05441  -11.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.4559     0.3808   -11.7   <2e-16 ***
## ---
## 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:  609 
## Degrees of Freedom for the fit:  6.657584
##       Residual Deg. of Freedom:  602.3424 
##                       at cycle:  5 
##  
## Global Deviance:     8659.499 
##             AIC:     8672.814 
##             SBC:     8702.186 
## ******************************************************************

model<-gamlss(cd4~pb(age),data=CD4, family=ZINBI(), trace=FALSE)

Chúng ta kết thúc bài hôm nay bằng 2 hình vẽ đẹp mắt biểu diễn đồ thị của 1 mô hình zeroinflated negatitive binomial kèm theo spline bù trừ nhằm ước tính số lượng tế bào bạch cầu CD4 theo tuổi ở trẻ nhỏ.

summary(model)

## ******************************************************************
## Family:  c("ZINBI", "Zero inflated negative binomial type I") 
## 
## Call:  gamlss(formula = cd4 ~ pb(age), family = ZINBI(), data = CD4,  
##     trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.67254    0.06711  114.32   <2e-16 ***
## pb(age)     -0.43241    0.02284  -18.93   <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) -0.62324    0.05441  -11.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.4559     0.3808   -11.7   <2e-16 ***
## ---
## 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:  609 
## Degrees of Freedom for the fit:  6.657584
##       Residual Deg. of Freedom:  602.3424 
##                       at cycle:  5 
##  
## Global Deviance:     8659.499 
##             AIC:     8672.814 
##             SBC:     8702.186 
## ******************************************************************

plot(model)

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  0.02955598 
##                        variance   =  0.9040999 
##                coef. of skewness  =  -0.3523449 
##                coef. of kurtosis  =  2.879732 
## Filliben correlation coefficient  =  0.9942961 
## ******************************************************************

CD4$predict=predict(model,type="response")
CD4%>%ggplot()+
  geom_point(shape=21,aes(x=age,y=cd4,color=age,fill=age),show.legend = F)+
  geom_smooth(aes(x=age,y=predict),se=F,color="darkviolet",linetype=2)+
  theme_bw()+
  scale_fill_gradient2(high="#8e00b2",mid="#ff9c07",low="#ef9b09")+
  scale_color_gradient2(high="#8e00b2",mid="#ff9c07",low="#ef9b09")

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

library(gamlss.util)

## Warning: package 'gamlss.util' was built under R version 3.4.1

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(viridisLite)

plotSimpleGamlss(y=cd4,x=age,formula=cd4~pb(age),family=ZINBI(),data=CD4,cols=viridis(option="A",n=100),x.val=c(1,2,3,4,5,6,7,8),val=200,N=1000,ylim=c(0,2300))

## new prediction

8 Kết luận

Trong bài thực hành thứ 2 của series về Gamlss này, Nhi muốn chuyển đến các bạn 3 thông điệp như sau:

1. Khái niệm về thành phần phân phối trong cấu trúc mô hình GAM

2. Cách thăm dò và chọn lọc quy luật phân phối phù hợp cho mô hình

3. Quy trình Stepwise và Kfold cross validation

Gamlss là một thư viện khổng lồ về lý thuyết xác suất. Tuy nhiên chúng ta chỉ mới hé mở cánh cửa này. Hy vọng với sự tò mò của mình, các bạn có thể tự mình khám phá hơn 100 họ phân phối còn lại cho biến số thực, biến liên tục dương, biến liên tục giới hạn, tỉ lệ, xác suất… Trong những bài tiếp theo chúng ta sẽ tiếp tục khám phá cấu trúc bên trong mô hình GAM bao gồm thành tố hồi quy và các yếu tố bù trừ.

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