Tìm Phân Phối Phù Hợp Cho Một Biến Số Trong Dataset Bằng fitdistrplus Package

Trong một số phép phân tích mô hình hồi qui biến số phụ thuộc (response variable) được yêu cầu xác định loại phân phối trước khi đưa vào mô hình hồi qui.

Bài viết này đề cập đến các functions đơn giản nhưng khá hiệu quả trong việc xác định phân phối khả dĩ phù hợp cho các biến số.

Mô tả bằng Biểu đồ phân phối của một biến số

Trước hết chúng ta phải biết biến số đó thuộc loại gì: liên tục, nhị phân, số đếm hay rời rạc. Việc có thể làm đầu tiên là biểu diễn biến số đó qua biểu đồ để xem hình dạng của nó giống với mô hình phân phối nào từ cái nhìn đầu tiên.

Tôi chọn dataset WDBC.data để minh họa các function này như sau:

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

Trước hết với V10, chúng ta tóm tắt để có cái nhìn tổng quan về biến số này

library(MASS)
library(survival)
library(lattice)
summary(df1$V10)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.02031 0.03350 0.04892 0.07400 0.20120
truehist(df1$V10,nbins = 30)

densityplot(~ df1$V10, plot.points = FALSE)

Đây là biến số liên tục, có các giá trị từ 0 đến 0.2012, lệch trái, không cân đối. Khi có thêm thông tin về skewness và kurtosis của V10 chúng ta có thể hình dung biến số này giống với phân phối nào. Chúng ta dùng hàm descdist() của package fitdistrplus

library(fitdistrplus)
descdist(df1$V10)

## summary statistics
## ------
## min:  0   max:  0.2012 
## median:  0.0335 
## mean:  0.04891915 
## estimated sd:  0.03880284 
## estimated skewness:  1.17118 
## estimated kurtosis:  4.066556

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ấ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.

Để chắc chắn hơn về kurtosis và skewness của V10, chúng ta dùng kĩ thuật boostraping với 1000 lần rút mẫu ngẫu nhiên sẽ có kết quả như sau:

descdist(df1$V10, boot = 1000)

## summary statistics
## ------
## min:  0   max:  0.2012 
## median:  0.0335 
## mean:  0.04891915 
## estimated sd:  0.03880284 
## estimated skewness:  1.17118 
## estimated kurtosis:  4.066556

Kết quả vẫn tương tự ở trên. Vùng màu vàng là biến thiên của V10 trong 1000 lần chọn mẫu boostrap.

Đánh giá trong mô hình phân phối

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. Vì Gamma nhận các giá trị (0,vô cực), cho nên khi có các giá trị bằng không mô hình sẽ báo lỗi. Trong phần summary ta thấy min V10 = 0. Để chạy được Gamma ta thay các giá trị 0 bằng 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.

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 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, 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.

Ngoài ra chúng ta còn có thể kiểm tra với phân phối khác, normal chẳng hạn, rồi so sánh chỉ số AIC để xem phân phối nào thực sự phù hợp hơn.

library(fitdistrplus)
f1gam <- fitdist(df1$V10, "gamma")
f1norm <- fitdist(df1$V10, "norm")
f1gam$aic
## [1] -2314.722
f1norm$aic
## [1] -2079.993

Chỉ số aic của gamma cho thấy thấp hơn (-2314.722) so với normal (-2079.993), vì thế phân phối gamma phù hợp hơn.

Cám ơn các bạn đã đọc bài viết. Các bạn nào có kinh nghiệm trong chủ đề này xin vui lòng góp ý để chúng ta hiểu rõ hơn về data distribution.