VARIABLE DISTRIBUTIONS

FITTING A DISTRIBUTION: Những Hạn Chế

1 . Giới thiệu

Trong khi đi tìm phân phối khả dĩ cho một biến số chúng ta có thể cậy nhờ đến một vài công cụ, đó là những hàm (R functions) có thể giúp chỉ ra những phân phối khả dĩ cho biến số mà ta đang khảo sát.

Các hàm như fitDist() trong Gamlss, descdistr() trong fitdistrplus đều có những mặt hữu ích nhưng cũng có những hạn chế, đôi khi làm chúng ta nhầm lẫn những kết quả và khó có lời giải thích hợp lí.

Chúng ta lần lượt đi qua các hàm này để xem có những vấn đề gì cần làm sáng tỏ.

2 . fitDist() của Gamlss

Để minh họa, tôi lấy dataset WDBC.data để phân tích.

library(readr)
df<-read.table(url("https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data"),header = FALSE, sep = ",")
df1<-df[,c(1,2,4,7,10,14)]

Tôi chọn V10 làm response variable (biến số Y) cho một mô hình regression từ các biến số V4, V7 và V14. Như vậy trước hết tôi sẽ tìm một phân phối khả dĩ cho V10.

Trong phần dưới đây, sẽ có thử phân tích V10 với phân phối Gamma (chỉ nhận những giá trị dương: lớn hơn 0). Vì thế, để chạy được Gamma tôi thay các giá trị V10 = 0 bằng V10 = 0.0001, tức là rất nhỏ để không làm thay đổi V10, nhưng đủ để V10 chạy trong phân phối Gamma.

library(gamlss)
library(MASS)
df1$V10[df1$V10==0]<-0.0001
m0 <- fitDist(V10, data=df1, type="realline")
m0$fits

##     SHASH    SHASHo   SHASHo2       SN2    exGAUS      EGB2      SEP3 
## -2333.057 -2330.291 -2330.291 -2326.165 -2317.881 -2316.819 -2311.261 
##       ST2       ST1       SN1      JSUo       JSU      SEP4       SST 
## -2308.943 -2308.878 -2307.068 -2306.747 -2306.747 -2304.786 -2304.230 
##       ST3       ST5        RG       ST4      SEP1        LO        TF 
## -2304.230 -2291.493 -2250.227 -2162.998 -2111.418 -2095.586 -2095.347 
##       TF2      SEP2        PE       PE2        NO        GT        GU 
## -2095.347 -2092.622 -2080.633 -2080.633 -2079.993 -1943.216 -1814.237

Kết quả cho ra 28 phân phối có thể có cho V10 được sắp xếp theo thứ tự từ tốt nhất đến xấu nhất dựa trên chỉ số AIC ở bên dưới mỗi loại phân phối. Ở trên đây, phân phối SHASH (Sinh-Arshinh) là tốt nhất với AIC = - 2332.432 (thấp nhất).

Tuy nhiên theo các Tác giả của cuốn sách “Flexible Regression and Smoothing Using GAMLSS in R” thì các phân phối trên được gọi là marginal distributions, phù hợp khi phân tích V10 đứng riêng một mình. Khi có hiện diện của các biến số khác trong một mô hình hồi qui (như V4, V7 và V14 ở trên) thì đôi khi cho thấy một phân phối khác phù hợp hơn, được gọi là conditional distributions.

Để chứng minh cho điều này tôi lập một mô hình hồi qui và phân tích xem mô hình đó tốt nhất trong dạng phân phối nào của V10. Câu trả lời có còn là phân phối SHASH nữa hay không.

library(gamlss)
df1$V10[df1$V10==0]<-0.0001
m_norm <- gamlss(V10 ~ V4+V7+V14, family=NO, data=df1, trace=FALSE)
m_gamma <- gamlss(V10 ~ V4+V7+V14, family=GA, data=df1, trace=FALSE)
m_SHASH <- gamlss(V10 ~ V4+V7+V14, family=SHASH, data=df1, trace=FALSE)
m_exGAU <- gamlss(V10 ~ V4+V7+V14, family=exGAUS, data=df1, trace=FALSE)
GAIC(m_norm,m_gamma,m_SHASH,m_exGAU)

##         df       AIC
## m_gamma  5 -2559.335
## m_exGAU  6 -2476.131
## m_SHASH  7 -2469.193
## m_norm   5 -2386.096

Kết quả cho thấy mô hình tốt nhất là m_gamma với giá trị AIC thấp nhất (-2559.335). Phù hợp với V10 là phân phối gamma chứ không phải là SHASH như từ kết quả của FitDist() trong package Gamlss. Gamma (GA) thậm chí còn không nằm trong các phân phối tốt nhất mà fitDist() gợi ý cho chúng ta.

Nói cách khác, FitDist() cho ta những marginal distributions tốt nhất cho V10 khi nó được phân tích riêng một mình. Nhưng khi đưa vào mô hình hồi qui với sự hiện diện của các biến số khác, V4, V7 và V14, thì đã có những thay đổi. Lúc này gamma là phân phối phù hợp nhất (conditrional distribution) cho V10 trong mô hình hồi qui trên.

3 . descdist() của fitdistrplus

library(fitdistrplus)
descdist(df1$V10,discrete = FALSE,boot = 1000)

## summary statistics
## ------
## min:  1e-04   max:  0.2012 
## median:  0.0335 
## mean:  0.04892143 
## estimated sd:  0.03879996 
## estimated skewness:  1.171546 
## estimated kurtosis:  4.067015

Nhìn vào biểu đồ trên với các chú thích như sau:

*: là vị trí tương đương kurtosis ~3, skewness = 0, đó là phân phối normal

Tam giác: là phân phối gần normal

+: là phân phối logistic

Dãi bôi xám là phân phối beta

Dấu tròn màu xanh chỉ v10 nằm trong vùng của phân phối beta và gần với gamma (đường vạch rời thô). Tuy nhiên vì phân phối beta nhận các giá trị trong khoảng (0-1), cho nên nhiều khả năng V10 của chúng ta gần với gamma.

Sau khi xác định V10 có nhiều khả năng có phân phối Gamma, chúng ta kiểm tra lại các thông số của V10 trong mô hình Gamma.

df1$V10[df1$V10==0]<-0.0001
library(fitdistrplus)
library(logspline)
f1gam <- fitdist(df1$V10, "gamma")
plot(f1gam)

summary(f1gam)

## Fitting of the distribution ' gamma ' by maximum likelihood 
## Parameters : 
##        estimate Std. Error
## shape  1.302095 0.06950869
## rate  26.615079 1.72472738
## Loglikelihood:  1159.361   AIC:  -2314.722   BIC:  -2306.035 
## Correlation matrix:
##           shape      rate
## shape 1.0000000 0.8237657
## rate  0.8237657 1.0000000

Biểu đồ so sánh các giá trị thật của V10 so với phân phối Gamma trên lí thuyết. Ngoại trừ biểu đồ Q-Q Plot có sai biệt ở phần cuối, thì ba biểu đồ về density, probability và CDF đều khá tương đồng với lí thuyết.

Vì thế chúng ta chấp nhận gamma là phân phối phù hợp của biến số V10 để đưa vào các mô hình phân tích.

Tuy nhiên, hạn chế của hàm descdidt() là chỉ cho ta biết V10 có thể là một vài phân phối cơ bản ở trên, không có nhiều phân phối phù hợp với rất nhiều kiểu data khác thường trong thực tế.

4 . Tóm lại

Như một người bạn nhiều kinh nghiệm của tôi đã nói: “Trên thực tế chẳng bao giờ tìm được một biến số nào phân phối chuẩn cả”.

Vì thế việc xác định phân phối của một biến số trong data thực tế chúng ta cần có cái nhìn hoài nghi về kết quả của các công cụ máy tính. Phải dựa vào các đặc tính cơ bàn của biến số như liên tục, số đếm, rời rạc, factor, hay nhị phân, cũng như các biểu đồ histogram, density để có sự cảm nhận, phán đoán ban đầu về biến số.

Cám ơn các bạn đã đọc bài viết này. Có điều gì chưa đúng xin vui lòng chỉ giúp.