Sử dụng package Gamlss

Bài 3: Suy diễn thống kê với mô hình Gamlss

1 Chúng ta dựng mô hình hồi quy để làm gì ?

Thân chào các bạn, hôm nay chúng ta lại gặp nhau trong một bài thực hành mới của series về package. Chủ đề của bài hôm nay là Suy diễn thống kê với mô hình Gamlss. Khi viết bài này, Nhi giả định rằng các bạn đều đã nắm rõ những kiến thức cơ bản về suy diễn thống kê (hay diễn giải) kết quả của một mô hình hồi quy, như chúng ta vẫn hay làm với bất cứ package hồi quy nào khác trong R, do đó bài này không phức tạp chút nào, thậm chí có thể sẽ nhàm chán với các bạn đã có thừa kinh nghiệm về GLM. Tuy nhiên Nhi hy vọng sẽ mang đến cho các bạn ít nhất là một vài phút giải trí thú vị.

Không biết có bao giờ các bạn thắc mắc tự hỏi mình : Trong phân tích dữ liệu nói chung và nghiên cứu y học nói riêng, chúng ta dựng mô hình hồi quy để làm gì ? Mô hình hồi quy có 3 ứng dụng chính như sau :

1. Để thăm dò dữ liệu khi chưa có giả thuyết hay ý tưởng nào cụ thể rõ ràng. Thí dụ đơn giản nhất là khảo sát sự tương quan giữa các biến số, nhận diện, khám phá một khuynh hướng hay quy luật mới nào đó. Các bạn thực hành data science chắc từng trải qua cảm giác tò mò khi phân tích một bộ dữ liệu theo cách hoàn toàn khác so với mục đích nguyên thủy của tác giả (thí dụ dataset iris nổi tiếng có thể dùng để minh họa cho rất nhiều phương pháp thống kê). Mô hình còn cho phép ta nhận diện những biến số quan trọng nổi bật có liên quan đến một đại lượng hay hiện tượng nào đó ĐÃ được ghi lại trong dữ liệu.

2. Mô hình là phương tiện để kiểm tra một giả thuyết cụ thể hoặc đi tìm câu trả lời cho một câu hỏi nghiên cứu. Hầu như tất cả câu hỏi thường gặp đều có thể diễn đạt thành mô hình hồi quy, thí dụ : so sánh 1 đại lượng giữa 2 hay nhiều phân nhóm, xác định liệu can thiệp điều trị có hiệu quả hay không ? khảo sát vai trò của yếu tố X lên kết quả hay hiện tượng Y … vân vân.

3. Ứng dụng thứ 3 của mô hình hồi quy là để tiên lượng, hoặc để phân loại. Mục tiêu lúc này là để ước tính giá trị 1 đại lượng mà ta chưa hoặc không thể đo được, hay khả năng của một điều đã xảy ra mà ta không có khả năng nhận thức, hay xa hơn là dự đoán cho một trường hợp giả định, ở tương lai, cho một cá thể bất kì…

Nói cách khác, mô hình cho phép ta trả lời 3 câu hỏi : điều gì đã xảy ra ? Tại sao nó xảy ra ? Và điều gì sẽ xảy ra ? Trong 3 trường hợp này, ta đều phải suy diễn hay diễn giải nội dung mô hình. Nội dung suy diễn có thể mang tính chất định lượng (kich thước của hiệu ứng, khác biệt) hoặc định tính (khuynh hướng : tăng/giảm, có hay không có ý nghĩa thống kê ? độ xác tín ?).

2 Thí dụ minh họa: Hồi quy logistic

Để minh họa cho cả 3 trường hợp, Nhi chọn một bài toán Hồi quy Logistic với dataset Biopsy trong package MASS. Dataset này bao gồm 9 chỉ số về đặc tính giải phẫu bệnh lý của mô sinh thiết ở bệnh nhân có khối u Vú, và 1 biến định tính nhị phân là phân loại Lành tính hay ác tính của khối u. Nếu đặt phân loại Lành/Ác tính như biến kết quả Y được mã hóa bằng 2 con số 0 (lành tính) và 1 (ác tính), việc dựng một mô hình hồi quy cho phép trả lời cả 2 câu hỏi :

Tôi muốn biết vai trò của mỗi đặc tính tế bào học góp phần định nghĩa tính chất U lành hay ác tính ? (bài toán thăm dò /diễn dịch)

Tôi muốn phân loại một mẫu sinh thiết bất kì dựa vào các đặc tính tế bào học (mục tiêu phân loại)

Dĩ nhiên, muốn diễn giải mô hình, ta phải dựng mô hình trước đã… Nhi từng nói với các bạn là mô hình của gamlss có cấu trúc rất tinh tế và gamlss là một bảo bối lợi hại, nhưng không có nghĩa là nó chỉ dành cho những vấn đề phức tạp. Khi bạn đã có nội lực, rành rẽ chiêu thức và định rõ mục tiêu thì bất cứ thanh kiếm nào trong tay bạn cũng có thể trở thành vũ khí lợi hại cả. Như vậy không có gì khác biệt khi bạn triển khai mô hình logistic với glm, MASS, brms, gam, nlme, lme4, caret và gamlss. Câu hỏi của bài hôm nay do đó không phải là Gamlss có làm được hồi quy logistic hay không ? mà là : nó có tương thích với những quy trình, hàm thông dụng để khai thác kết quả từ mô hình logistic này không ?

library(tidyverse)

library(gamlss)

nC<-detectCores()

my_theme <- function(base_size =8, base_family = "sans"){
  theme_bw(base_size = base_size, base_family = base_family) +
    theme(
      panel.grid.major = element_line(color = "gray"),
      panel.grid.minor = element_blank(),
      panel.background = element_rect(fill = NA),
      strip.background = element_rect(fill = "#9013af", color = NA, size =0.5),
      strip.text = element_text(face = "bold", size = 8, color = "white"),
      legend.position = "bottom",
      legend.justification = "center",
      legend.background = element_blank(),
      legend.margin = margin(0.5,0.5,0.5,0.5)
    )
}

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

3 Thăm dò số liệu

df=read.csv("https://raw.github.com/vincentarelbundock/Rdatasets/master/csv/MASS/biopsy.csv")

df=read.csv("http://vincentarelbundock.github.io/Rdatasets/csv/MASS/biopsy.csv")%>%as_tibble()%>%.[,c(3:12)]%>%na.omit()

names(df)=c("clumpthickness",
            "SizeUniformity",
            "ShapeUniformity",
            "Margin_adhesion",
            "EpiCellSize",
            "Barenuclei",
            "BlandChromatin",
            "NormalNucleoli",
            "Mitoses",
            "Class"
)

df2=df%>%mutate(.,Class=as.integer(.$Class)-1L)

df%>%gather(clumpthickness:Mitoses,key="Features",value="Score")%>%
  ggplot(aes(x=Score,
             fill=Class))+
  geom_histogram(alpha=0.8,aes(y=..density..),color="black",bins=10)+
  my_theme()+
  facet_wrap(~Features,ncol=2,scales = "free")+
  scale_fill_manual(values=c("#ffe314","#ff1499"))

df%>%gather(clumpthickness:Mitoses,key="Features",value="Score")%>%
  ggplot(aes(x=Class,y=Score,
             fill=Class,
             color=Class))+
  geom_jitter(alpha=0.3,shape=21)+coord_flip()+
  my_theme()+
  facet_wrap(~Features,ncol=2,scales = "free")+
  scale_color_manual(values=c("#ffe314","#ff1499"))+
  scale_fill_manual(values=c("#ffe314","#ff1499"))

df%>%gather(clumpthickness:Mitoses,key="Features",value="Score")%>%
  ggplot(aes(x=Features,
             fill=factor(Score)))+coord_flip()+
  geom_bar(alpha=0.6,color="black")+
  my_theme()+
  scale_fill_manual(values=mycolors)+
  facet_wrap(~Class,ncol=1,scales = "free")

dfscale<-df[,-10]%>%as.matrix()%>%scale()%>%as_tibble()%>%mutate(.,Class=df$Class,Id=row.names(.))

library(viridis)

dfscale%>%gather(clumpthickness:Mitoses,key="Criteria",value="Score")%>%
  ggplot(aes(x=reorder(Id,-Score),
             y=reorder(Criteria,Score),
             fill=Score))+
  geom_tile(show.legend=T)+
  facet_wrap(~Class,ncol=2,shrink=T,scale="free")+
  scale_fill_viridis(option="C",begin=1,end=0)+
  my_theme()+
  theme(axis.text.y=element_blank(),axis.text.x = element_text(angle =45,hjust=1,vjust=1))+
  coord_flip()+
  scale_y_discrete("Criteria")+
  scale_x_discrete("Patient's Id")

df2%>%gather(clumpthickness:NormalNucleoli,key="Features",value="Score")%>%
  ggplot(aes(x=Score,y=Class,
                     color=Features,fill=Features))+
  geom_smooth(show.legend = F,se=T)+
  my_theme()+
  facet_wrap(~Features,ncol=3,scales = "free")+
  scale_color_manual(values=rev(mycolors))+
  scale_fill_manual(values=rev(mycolors))+
  scale_y_continuous("Probability of Malignant")

4 Quy trình dựng mô hình Logistic (theo trường phái Statistical learning)

Nhi sẽ hướng dẫn các bạn từng bước trong quy trình dựng một mô hình logistic theo trường phái Statistical learning (hay Machine learning).

Khác với quy trình cổ điển, chúng ta không sử dụng toàn bộ dataset mà chỉ dùng 1 phần (500 cases) để dựng mô hình, và để dành phần còn lại cho mục đích kiểm định mô hình. Nhi dùng package caret để cắt dữ liệu một cách ngẫu nhiên nhưng vẫn bảo toàn được tỉ lệ giữa 2 phân nhóm.

library(caret)

set.seed(123)
id=createDataPartition(y=df$Class, p=499/683,list=FALSE)
trainset=df2[id,]
testset=df2[-id,]

p1=trainset%>%as.data.frame()%>%ggplot(aes(x=as.factor(Class),fill=as.factor(Class)))+stat_count(show.legend=F,color="black",alpha=0.7)+scale_fill_manual(values=c("#ffe314","#ff1499"))+scale_x_discrete("Class")+coord_flip()+ggtitle("Train")+my_theme()

p2=testset%>%as.data.frame()%>%ggplot(aes(x=as.factor(Class),fill=as.factor(Class)))+stat_count(show.legend=F,color="black",alpha=0.7)+scale_fill_manual(values=c("#ffe314","#ff1499"))+scale_x_discrete("Class")+coord_flip()+ggtitle("Test")+my_theme()

p3=df2%>%as.data.frame()%>%ggplot(aes(x=as.factor(Class),fill=as.factor(Class)))+stat_count(show.legend=F,color="black",alpha=0.7)+scale_fill_manual(values=c("#ffe314","#ff1499"))+scale_x_discrete("Class")+coord_flip()+ggtitle("Origin")+my_theme()

library(gridExtra)

grid.arrange(p1,p2,p3,ncol=1)

5 Chọn lọc biến số bằng phương pháp BMA

Bước tiếp theo, Nhi sẽ sử dụng phương pháp Bayesian Model Averaging (BMA) để nhận diện những biến số quan trọng nhất (có xác suất cao nhất hiện diện trong một mô hình logistic tối ưu, dựa vào tiêu chí Bayes factor).

library(BMA)

bmalog=bic.glm(f=Class~.,data=trainset,glm.family ="binomial",OR=20)

plot(bmalog)

summary(bmalog)

## 
## Call:
## bic.glm.formula(f = Class ~ ., data = trainset, glm.family = "binomial",     OR = 20)
## 
## 
##   31  models were selected
##  Best  5  models (cumulative posterior probability =  0.462 ): 
## 
##                  p!=0    EV        SD       model 1     model 2   
## Intercept        100    -8.939767  1.16294     -9.0713     -9.5103
## clumpthickness   100.0   0.571349  0.15234      0.5829      0.6866
## SizeUniformity    12.7   0.044994  0.13358       .           .    
## ShapeUniformity   69.9   0.357343  0.28099      0.4824       .    
## Margin_adhesion   77.9   0.258642  0.17304      0.3152      0.3506
## EpiCellSize        4.1   0.009441  0.05611       .           .    
## Barenuclei       100.0   0.388159  0.09925      0.3466      0.4175
## BlandChromatin    70.7   0.333547  0.26174      0.4424      0.5174
## NormalNucleoli    43.1   0.113937  0.15050       .          0.3121
## Mitoses           27.2   0.174762  0.36006       .           .    
##                                                                   
## nVar                                              5           5   
## BIC                                         -2972.5060  -2971.6874
## post prob                                       0.150       0.100 
##                  model 3     model 4     model 5   
## Intercept           -8.4581     -8.4262     -9.4227
## clumpthickness       0.6041      0.5571      0.5122
## SizeUniformity        .           .           .    
## ShapeUniformity      0.6433      0.5086      0.4449
## Margin_adhesion      0.3656      0.3345      0.3169
## EpiCellSize           .           .           .    
## Barenuclei           0.3952      0.3906      0.3555
## BlandChromatin        .           .          0.4291
## NormalNucleoli        .          0.2384       .    
## Mitoses               .           .          0.6333
##                                                    
## nVar                   4           5           6   
## BIC              -2971.4752  -2970.8173  -2970.5663
## post prob            0.090       0.065       0.057

imageplot.bma(bmalog,color = c("#fc024d", "#0294fc", "white"))

6 Dựng mô hình Logistic bằng Gamlss

Dựa Vào kết quả của BMA, mô hình tối ưu nhất (số 1) sẽ bao gồm 5 biến số

Từ thông tin này, chúng ta sẽ đưa 5 biến số nói trên vào Gamlss. Do đã xác định trước về phân phối của biến kết quả (Binomial), ta không cần phải làm bước thăm dò phân phối như bài 2 đã mô tả. Cú pháp hàm gamlss cho mô hình logistic hoàn toàn giống như hàm glm().

Nhi cũng dựng sẵn một mô hình m0 chỉ chứa hằng số, mô hình này sẽ được dùng cho Likelihood ratio test.

m0=gamlss(Class~1,
          data=trainset,
          family=BI(),
          trace=F)

m1=gamlss(Class~clumpthickness+
            ShapeUniformity+
            Margin_adhesion+
            Barenuclei+
            BlandChromatin,
      data=trainset,
         family=BI(),
              trace=F)

Kết quả của LR test khẳng định tính phù hợp dữ liệu của mô hình logistic m1

LR.test(m0,m1)

##  Likelihood Ratio Test for nested GAMLSS models. 
##  (No check whether the models are nested is performed). 
##  
##        Null model: deviance= 647.4466 with  1 deg. of freedom 
##  Altenative model: deviance= 97.51044 with  6 deg. of freedom 
##  
##  LRT = 549.9362 with 5 deg. of freedom and p-value= 0

Log likelihood có thể tính như sau :

logL<-gen.likelihood(m1)

logL()

## [1] 48.75522

Đầu tiên ta sẽ vẽ một số biểu đồ để khảo sát đặc tính của residual, kết quả thăm dò cho thấy mô hình khá hợp lý :

newpar<-par(mfrow=c(2,2),mar=par("mar")+c(0,1,0,0),col.axis="blue4",
            col="blue4", col.main="blue4",col.lab="blue4",pch=16,
            cex=.45, cex.lab=1.2, cex.axis=1, cex.main=1.2)
plot(m1,par=newpar)

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.03426615 
##                        variance   =  0.9879639 
##                coef. of skewness  =  -0.1808607 
##                coef. of kurtosis  =  3.085871 
## Filliben correlation coefficient  =  0.9984326 
## ******************************************************************

Biểu đồ wormplot xác nhận sự phù hợp của mô hình

wp(m1)

rqres.plot(m1)

Hàm HLgof.test cho phép làm kiểm định Hosmer-Lemeshow

library(MKmisc)

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

HLgof.test(fit = fitted(m1), obs = trainset$Class) 

## $C
## 
##  Hosmer-Lemeshow C statistic
## 
## data:  fitted(m1) and trainset$Class
## X-squared = 5.0901, df = 8, p-value = 0.7479
## 
## 
## $H
## 
##  Hosmer-Lemeshow H statistic
## 
## data:  fitted(m1) and trainset$Class
## X-squared = 16.386, df = 8, p-value = 0.03717

7 Khai thác nội dung mô hình logistic trong gamlss

Như mọi package về hồi quy khác, chúng ta có thể dùng hàm summary để xem nội dung bên trong mô hình vừa dựng

Với phân phối binomial, chỉ có tham số Mu được mô hình hóa, với linkfunction là 1 hàm logit. Toàn bộ 5 biến số trong mô hình đều có ý nghĩa thống kê với các mức độ khác nhau. Hệ số hồi quy >0 gợi ý rằng 5 chỉ số này có tương quan thuận với nguy cơ khối U ác tính. Mô hình có Global deviance = 97.5 AIC = 109.5 và BIC = 134.8

summary(m1)

## ******************************************************************
## Family:  c("BI", "Binomial") 
## 
## Call:  gamlss(formula = Class ~ clumpthickness + ShapeUniformity +  
##     Margin_adhesion + Barenuclei + BlandChromatin,  
##     family = BI(), data = trainset, trace = F) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  logit
## Mu Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -9.07134    1.06824  -8.492 2.40e-16 ***
## clumpthickness   0.58294    0.14132   4.125 4.35e-05 ***
## ShapeUniformity  0.48237    0.15836   3.046 0.002443 ** 
## Margin_adhesion  0.31521    0.11347   2.778 0.005679 ** 
## Barenuclei       0.34662    0.09698   3.574 0.000386 ***
## BlandChromatin   0.44241    0.17665   2.504 0.012587 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  500 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  494 
##                       at cycle:  2 
##  
## Global Deviance:     97.51044 
##             AIC:     109.5104 
##             SBC:     134.7981 
## ******************************************************************

Tương tự, hàm confint sẽ xuất ra khoảng tin cậy 97.5% của hệ số hồi quy

coef(m1)%>%exp()%>%cbind()

##                            .
## (Intercept)     0.0001149129
## clumpthickness  1.7912905633
## ShapeUniformity 1.6199165164
## Margin_adhesion 1.3705523216
## Barenuclei      1.4142734968
## BlandChromatin  1.5564473389

confint(m1)

##                        2.5 %     97.5 %
## (Intercept)     -11.16504836 -6.9776244
## clumpthickness    0.30594745  0.8599252
## ShapeUniformity   0.17199225  0.7927570
## Margin_adhesion   0.09281228  0.5376153
## Barenuclei        0.15654793  0.5366840
## BlandChromatin    0.09617747  0.7886343

Tuy nhiên, ta phải chuyển đổi hệ số hồi quy của mô hình logistic thành Odds-ratio (OR) để có thể diễn giải ý nghĩa của chúng.

exp(cbind(OR = coef(m1), confint(m1)))

##                           OR        2.5 %       97.5 %
## (Intercept)     0.0001149129 1.416058e-05 0.0009325158
## clumpthickness  1.7912905633 1.357911e+00 2.3629840224
## ShapeUniformity 1.6199165164 1.187669e+00 2.2094795327
## Margin_adhesion 1.3705523216 1.097256e+00 1.7119196595
## Barenuclei      1.4142734968 1.169467e+00 1.7103260169
## BlandChromatin  1.5564473389 1.100954e+00 2.2003892700

Bây giờ thì mọi thứ đã rõ ràng hơn một chút. Ngưỡng trên và dưới của 97.5%CI của OR đều lớn hơn 1 ,cho thấy hiệu ứng của cả 5 biến số đều có ý nghĩa. Các bạn chắc đã quen thuộc với việc diễn giải OR rồi ? Gợi ý là bạn có thể phát biểu như sau : Chỉ số ShapeUniformity cứ tăng thêm 1 điểm thì khả năng U ác tính sẽ tăng trung bình 1.62 lần (1.187 tới 2.21), với giả định là các chỉ số còn lại không thay đổi.

Nhi vẽ thêm 2 biểu đồ « Marginal effect » của mô hình như sau :

probs_predicted <- predict(m1,type = "response") 

trainset%>%mutate(.,Pred=probs_predicted)%>%
  gather(c(clumpthickness,
           ShapeUniformity,
           Margin_adhesion,
           Barenuclei,
           BlandChromatin)
         ,key="Features",value="Score")%>%
  ggplot(aes(x=Score,y=Pred,
             color=Features,fill=Features))+
  geom_smooth(show.legend = F,se=T)+
  my_theme()+
  facet_wrap(~Features,ncol=3,scales = "free")+
  scale_color_manual(values=(mycolors))+
  scale_fill_manual(values=(mycolors))+
  scale_y_continuous("Probability of Malignant")

ordf=exp(cbind(OR = coef(m1), confint(m1)))%>%as_tibble()%>%mutate(.,Predictor=colnames(m1$mu.x))
names(ordf)=c("OR","LL","UL","Predictor")

ordf%>%ggplot(aes(x=reorder(Predictor,OR),color=reorder(Predictor,OR)))+
  geom_errorbar(aes(ymin=LL, ymax=UL),width=0.5,size=1)+
  geom_point(aes(y=OR),shape=21,size=5,fill="white",stroke=1.5)+
  my_theme()+coord_flip()+
  geom_hline(yintercept = 1,size=1,color="blue",linetype=2)+
  scale_color_manual(values=rev(mycolors))+scale_x_discrete("Predictors")+scale_y_continuous("Odds-Ratio")

8 Bootstrap cho mô hình logistic Gamlss

Thật thiếu sót nếu nói về suy diễn thống kê cho hệ số hồi quy mà không làm bootstrap. Sau đây Nhi sẽ thực hiện cả 3 quy trình bootstrap : Parametric, Semi-parametric và Non parametric :

Lưu ý là cho cả 3 hàm bootstrap, Nhi tính OR trực tiếp từ kết quả gốc. Mục tiêu cho cả 3 là ước tính khoảng tin cậy 95% cho OR của 5 biến số, dựa vào 1000 phiên bản mô hình logistic khác nhau dựa trên các mẫu được chọn ngẫu nhiên lặp lại

Parametric Bootstrap

Quy trình bootstrap parametric sử dụng hàm random với phân phối BI với Mu = giá trị dự báo của mô hình gốc để tạo ra 500 kết quả mô phỏng, sau đó dựng lại model dựa vào dữ liệu này.

#Parametric bootstrap

library(boot)

## 
## Attaching package: 'boot'

## The following object is masked from 'package:robustbase':
## 
##     salinity

## The following object is masked from 'package:survival':
## 
##     aml

## The following object is masked from 'package:lattice':
## 
##     melanoma

## The following object is masked from 'package:gamlss.data':
## 
##     aids

tdf=data.frame(trainset,fmu=fitted(m1))
          
paraboot=function(data,model,i){
  d=data
  d$Class<-if("original"%in%as.character(sys.call())) d$Class
  else rBI(nrow(d), mu=d$fmu)
  coef(update(model, data=d))%>%exp()
  }

sim<-boot(paraboot,data=tdf,model=m1, R=1000)

sim

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = tdf, statistic = paraboot, R = 1000, model = m1)
## 
## 
## Bootstrap Statistics :
##         original       bias    std. error
## t1* 0.0001149129 2.629046e-06 0.000136824
## t2* 1.7912905633 1.098208e-01 0.365027461
## t3* 1.6199165164 1.042150e-01 0.335162020
## t4* 1.3705523216 4.103648e-02 0.197446426
## t5* 1.4142734968 4.394960e-02 0.167335851
## t6* 1.5564473389 1.061723e-01 0.354000258

simdf=sim$t%>%as_tibble()
names(simdf)=c("Intercept",
               "clumpthickness",
               "ShapeUniformity",
               "Margin_adhesion",
               "Barenuclei",
               "BlandChromatin")

Hmisc::describe(simdf)

## simdf 
## 
##  6  Variables      1000  Observations
## ---------------------------------------------------------------------------
## Intercept 
##         n   missing  distinct      Info      Mean       Gmd       .05 
##      1000         0      1000         1 0.0001175 0.0001277 4.520e-06 
##       .10       .25       .50       .75       .90       .95 
## 9.721e-06 2.971e-05 7.043e-05 1.567e-04 2.845e-04 3.868e-04 
## 
## lowest : 5.148915e-09 1.561800e-08 2.209840e-08 4.660800e-08 8.156535e-08
## highest: 7.022593e-04 7.172752e-04 7.751278e-04 1.273804e-03 1.424281e-03
## ---------------------------------------------------------------------------
## clumpthickness 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.901   0.3702    1.473    1.541 
##      .25      .50      .75      .90      .95 
##    1.659    1.846    2.052    2.308    2.529 
## 
## lowest : 1.288253 1.301795 1.302413 1.305205 1.310806
## highest: 3.602368 3.650143 3.673579 3.949148 4.707007
## ---------------------------------------------------------------------------
## ShapeUniformity 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.724    0.358    1.305    1.359 
##      .25      .50      .75      .90      .95 
##    1.495    1.662    1.897    2.136    2.323 
## 
## lowest : 1.045541 1.113728 1.115859 1.118530 1.123820
## highest: 3.152913 3.186391 3.253372 3.267077 3.409626
## ---------------------------------------------------------------------------
## Margin_adhesion 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.412   0.2119    1.161    1.202 
##      .25      .50      .75      .90      .95 
##    1.274    1.383    1.514    1.651    1.770 
## 
## lowest : 0.9751873 0.9778242 0.9858825 1.0117606 1.0124739
## highest: 2.2759561 2.2865243 2.3131781 2.4038399 2.4216948
## ---------------------------------------------------------------------------
## Barenuclei 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.458   0.1823    1.231    1.268 
##      .25      .50      .75      .90      .95 
##    1.343    1.433    1.551    1.672    1.759 
## 
## lowest : 1.053116 1.112638 1.129052 1.140136 1.142995
## highest: 2.110598 2.138936 2.154450 2.225772 2.298108
## ---------------------------------------------------------------------------
## BlandChromatin 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.663   0.3796    1.195    1.275 
##      .25      .50      .75      .90      .95 
##    1.416    1.613    1.852    2.096    2.303 
## 
## lowest : 0.8833313 0.9570787 0.9706442 0.9860125 0.9864565
## highest: 3.0300504 3.0415275 3.1411191 3.3319302 4.2187493
## ---------------------------------------------------------------------------

simdf%>%gather(Intercept:BlandChromatin,key="Predictor",value="OR")%>%
  ggplot(aes(x=OR))+
  geom_histogram(aes(color=Predictor,fill=Predictor,OR),alpha=0.5,show.legend = F,bins=100)+
  geom_vline(xintercept = 1,color="blue",linetype=2,size=1)+
  my_theme()+
  facet_wrap(~Predictor,ncol=1,scales="free_y")+
  scale_color_manual(values=rev(mycolors))+
  scale_fill_manual(values=rev(mycolors))+
  scale_x_continuous("Bootstraped Odds-ratio")

simdf%>%gather(Intercept:BlandChromatin,key="Predictor",value="OR")%>%
  ggplot(aes(x=reorder(Predictor,OR),y=OR))+
  geom_boxplot(aes(color=reorder(Predictor,OR),fill=reorder(Predictor,OR)),alpha=0.7,show.legend = F)+
  geom_hline(yintercept = 1,color="blue",linetype=2,size=1)+
  my_theme()+coord_flip()+
  scale_color_manual(values=rev(mycolors))+
  scale_fill_manual(values=rev(mycolors))+
  scale_y_continuous("Bootstraped Odds-ratio")+
  scale_x_discrete("Predictors")

simdf%>%mutate(.,Id=c(1:1000))%>%gather(Intercept:BlandChromatin,key="Predictor",value="OR")%>%
  ggplot(aes(x=Id,y=OR))+
  geom_path(aes(color=reorder(Predictor,OR)),alpha=0.7,show.legend = F)+
  geom_hline(yintercept = 1,color="blue",linetype=2,size=1)+
  my_theme()+ facet_wrap(~Predictor,ncol=1,scales="free")+
  scale_color_manual(values=rev(mycolors))+
  scale_fill_manual(values=rev(mycolors))+
  scale_y_continuous("Bootstraped Odds-ratio")+
  scale_x_discrete("Iterations")

Semiparametric bootstrap

Quy trình bootstrap semi-parametric kết hợp cả giá trị dự báo cho Mu và Sai số để mô phỏng dữ liệu cho 1000 phiên bản mô hình

tdf=data.frame(trainset,res=resid(m1),fmu=fitted(m1))

semiparaboot <- function(data,model,i){
d <- data
d$Class <- qBI(pNO(d$res[i]),mu=d$fmu)
coef(update(model, data=d))%>%exp()
}

sim2<-boot(semiparaboot,data=tdf,model=m1,R=1000)

sim2

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = tdf, statistic = semiparaboot, R = 1000, model = m1)
## 
## 
## Bootstrap Statistics :
##         original        bias     std. error
## t1* 0.0001149129 -3.453094e-05 0.0001044969
## t2* 1.7912905633  1.791882e-01 0.4119618798
## t3* 1.6199165164  1.308356e-01 0.3884276941
## t4* 1.3705523216  5.890155e-02 0.2094159000
## t5* 1.4142734968  7.861180e-02 0.1924109250
## t6* 1.5564473389  1.167703e-01 0.3587770196

simdf=sim2$t%>%as_tibble()
names(simdf)=c("Intercept",
               "clumpthickness",
               "ShapeUniformity",
               "Margin_adhesion",
               "Barenuclei",
               "BlandChromatin")

Hmisc::describe(simdf)

## simdf 
## 
##  6  Variables      1000  Observations
## ---------------------------------------------------------------------------
## Intercept 
##         n   missing  distinct      Info      Mean       Gmd       .05 
##      1000         0      1000         1 8.038e-05  9.26e-05 2.329e-06 
##       .10       .25       .50       .75       .90       .95 
## 5.597e-06 1.593e-05 4.621e-05 1.018e-04 1.956e-04 2.664e-04 
## 
## lowest : 4.057892e-10 4.503245e-09 1.884417e-08 3.541566e-08 5.678025e-08
## highest: 7.265410e-04 7.422121e-04 7.514905e-04 8.241194e-04 9.145179e-04
## ---------------------------------------------------------------------------
## clumpthickness 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1     1.97   0.4093    1.489    1.565 
##      .25      .50      .75      .90      .95 
##    1.714    1.902    2.135    2.439    2.663 
## 
## lowest : 1.207148 1.253299 1.270413 1.273237 1.274463
## highest: 3.626823 3.897228 4.864944 5.313970 6.083641
## ---------------------------------------------------------------------------
## ShapeUniformity 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.751   0.3914    1.308    1.374 
##      .25      .50      .75      .90      .95 
##    1.496    1.673    1.912    2.211    2.432 
## 
## lowest : 1.043885 1.056044 1.125877 1.142147 1.160039
## highest: 3.733299 3.900524 3.948188 4.137698 5.147545
## ---------------------------------------------------------------------------
## Margin_adhesion 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.429   0.2253    1.155    1.203 
##      .25      .50      .75      .90      .95 
##    1.286    1.394    1.537    1.711    1.819 
## 
## lowest : 0.9676823 0.9716964 0.9873473 0.9881633 0.9955656
## highest: 2.2567667 2.2825977 2.3138345 2.5473243 2.6769643
## ---------------------------------------------------------------------------
## Barenuclei 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.493   0.2064    1.238    1.279 
##      .25      .50      .75      .90      .95 
##    1.367    1.462    1.593    1.734    1.825 
## 
## lowest : 1.040025 1.070333 1.105881 1.124974 1.141332
## highest: 2.303869 2.317018 2.368789 2.386330 2.450150
## ---------------------------------------------------------------------------
## BlandChromatin 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.673   0.3862    1.184    1.268 
##      .25      .50      .75      .90      .95 
##    1.438    1.628    1.864    2.123    2.285 
## 
## lowest : 0.8529395 0.8576573 0.9090451 0.9558478 0.9638761
## highest: 3.0184812 3.0255252 3.0688857 3.2787117 4.3402866
## ---------------------------------------------------------------------------

simdf%>%gather(Intercept:BlandChromatin,key="Predictor",value="OR")%>%
  ggplot(aes(x=OR))+
  geom_histogram(aes(color=Predictor,fill=Predictor,OR),alpha=0.5,show.legend = F,bins=100)+
  geom_vline(xintercept = 1,color="blue",linetype=2,size=1)+
  my_theme()+
  facet_wrap(~Predictor,ncol=1,scales="free_y")+
  scale_color_manual(values=rev(mycolors))+
  scale_fill_manual(values=rev(mycolors))+
  scale_x_continuous("Bootstraped Odds-ratio")

simdf%>%gather(Intercept:BlandChromatin,key="Predictor",value="OR")%>%
  ggplot(aes(x=reorder(Predictor,OR),y=OR))+
  geom_boxplot(aes(color=reorder(Predictor,OR),fill=reorder(Predictor,OR)),alpha=0.7,show.legend = F)+
  geom_hline(yintercept = 1,color="blue",linetype=2,size=1)+
  my_theme()+coord_flip()+
  scale_color_manual(values=rev(mycolors))+
  scale_fill_manual(values=rev(mycolors))+
  scale_y_continuous("Bootstraped Odds-ratio")+
  scale_x_discrete("Predictors")

simdf%>%mutate(.,Id=c(1:1000))%>%gather(Intercept:BlandChromatin,key="Predictor",value="OR")%>%
  ggplot(aes(x=Id,y=OR))+
  geom_path(aes(color=reorder(Predictor,OR)),alpha=0.7,show.legend = F)+
  geom_hline(yintercept = 1,color="blue",linetype=2,size=1)+
  my_theme()+ facet_wrap(~Predictor,ncol=1,scales="free")+
  scale_color_manual(values=rev(mycolors))+
  scale_fill_manual(values=rev(mycolors))+
  scale_y_continuous("Bootstraped Odds-ratio")+
  scale_x_discrete("Iterations")

Non parametric bootstrap

Quy trình bootstrap non parametric không sử dụng dữ liệu mô phỏng từ mô hình gốc, nhưng chọn mẫu ngẫu nhiên lặp lại từ dữ liệu nguyên thủy và với mỗi lần chọn mẫu như vậy gamlss dựng một mô hình mới. Đây là quy trình dễ hiểu hơn nhưng không đảm bảo tính ổn định về quy luật phân phối của biến kết quả.

Nó không phải lúc nào cũng hoạt động tốt, có khi mô hình không thể converge do dữ liệu không phù hợp giả định về phân phối.

nonparaboot <- function(data, model,i){
d<-data[i,]
coef(update(model,data=d))%>%exp()
}

sim3<-boot(nonparaboot,data=trainset,model=m1,R=1000)

sim3

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = trainset, statistic = nonparaboot, R = 1000, model = m1)
## 
## 
## Bootstrap Statistics :
##         original        bias     std. error
## t1* 0.0001149129 -8.461881e-06 0.0001147645
## t2* 1.7912905633  1.085037e-01 0.3135095642
## t3* 1.6199165164  1.377342e-01 0.4711545167
## t4* 1.3705523216  3.958977e-02 0.2128302002
## t5* 1.4142734968  8.396716e-02 0.3034540411
## t6* 1.5564473389  8.868889e-02 0.3817036564

simdf=sim3$t%>%as_tibble()
names(simdf)=c("Intercept",
               "clumpthickness",
               "ShapeUniformity",
               "Margin_adhesion",
               "Barenuclei",
               "BlandChromatin")

Hmisc::describe(simdf)

## simdf 
## 
##  6  Variables      1000  Observations
## ---------------------------------------------------------------------------
## Intercept 
##         n   missing  distinct      Info      Mean       Gmd       .05 
##      1000         0      1000         1 0.0001065 0.0001107 5.923e-06 
##       .10       .25       .50       .75       .90       .95 
## 1.139e-05 2.754e-05 7.108e-05 1.445e-04 2.465e-04 3.423e-04 
## 
## lowest : 2.669069e-09 3.917247e-08 9.523394e-08 2.421583e-07 3.716366e-07
## highest: 6.860885e-04 7.155737e-04 7.255835e-04 7.546797e-04 1.065762e-03
## ---------------------------------------------------------------------------
## clumpthickness 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1      1.9   0.3311    1.504    1.570 
##      .25      .50      .75      .90      .95 
##    1.692    1.854    2.051    2.268    2.464 
## 
## lowest : 1.225968 1.259917 1.305726 1.323940 1.331261
## highest: 3.042570 3.059415 3.607555 3.799907 4.406957
## ---------------------------------------------------------------------------
## ShapeUniformity 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.758   0.4546    1.255    1.323 
##      .25      .50      .75      .90      .95 
##    1.473    1.667    1.930    2.274    2.555 
## 
## lowest : 0.9619730 0.9910264 1.0323452 1.0622100 1.0803105
## highest: 3.9874989 4.0639971 4.5667426 5.1587490 7.4024994
## ---------------------------------------------------------------------------
## Margin_adhesion 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1     1.41   0.2275    1.128    1.176 
##      .25      .50      .75      .90      .95 
##    1.265    1.381    1.522    1.674    1.793 
## 
## lowest : 0.8763680 0.9122344 0.9349185 0.9944366 1.0024538
## highest: 2.1979817 2.2095690 2.3955696 2.4739375 3.1154100
## ---------------------------------------------------------------------------
## Barenuclei 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.498   0.2575    1.209    1.254 
##      .25      .50      .75      .90      .95 
##    1.336    1.450    1.582    1.752    1.916 
## 
## lowest : 1.045653 1.048473 1.086506 1.108247 1.115051
## highest: 3.498510 3.521692 4.552737 4.678712 5.366057
## ---------------------------------------------------------------------------
## BlandChromatin 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1000        0     1000        1    1.645   0.3973    1.174    1.246 
##      .25      .50      .75      .90      .95 
##    1.388    1.583    1.822    2.080    2.292 
## 
## lowest : 0.8837530 0.9503633 0.9990602 1.0050007 1.0055734
## highest: 3.1849364 3.4113796 3.4165698 3.8275135 4.9547556
## ---------------------------------------------------------------------------

simdf%>%gather(Intercept:BlandChromatin,key="Predictor",value="OR")%>%
  ggplot(aes(x=OR))+
  geom_histogram(aes(color=Predictor,fill=Predictor,OR),alpha=0.5,show.legend = F,bins=100)+
  geom_vline(xintercept = 1,color="blue",linetype=2,size=1)+
  my_theme()+
  facet_wrap(~Predictor,ncol=1,scales="free_y")+
  scale_color_manual(values=rev(mycolors))+
  scale_fill_manual(values=rev(mycolors))+
  scale_x_continuous("Bootstraped Odds-ratio")

simdf%>%gather(Intercept:BlandChromatin,key="Predictor",value="OR")%>%
  ggplot(aes(x=reorder(Predictor,OR),y=OR))+
  geom_boxplot(aes(color=reorder(Predictor,OR),fill=reorder(Predictor,OR)),alpha=0.7,show.legend = F)+
  geom_hline(yintercept = 1,color="blue",linetype=2,size=1)+
  my_theme()+coord_flip()+
  scale_color_manual(values=rev(mycolors))+
  scale_fill_manual(values=rev(mycolors))+
  scale_y_continuous("Bootstraped Odds-ratio")+
  scale_x_discrete("Predictors")

simdf%>%mutate(.,Id=c(1:1000))%>%gather(Intercept:BlandChromatin,key="Predictor",value="OR")%>%
  ggplot(aes(x=Id,y=OR))+
  geom_path(aes(color=reorder(Predictor,OR)),alpha=0.7,show.legend = F)+
  geom_hline(yintercept = 1,color="blue",linetype=2,size=1)+
  my_theme()+ facet_wrap(~Predictor,ncol=1,scales="free")+
  scale_color_manual(values=rev(mycolors))+
  scale_fill_manual(values=rev(mycolors))+
  scale_y_continuous("Bootstraped Odds-ratio")+
  scale_x_discrete("Iterations")

9 Tinh chỉnh cho mô hình logistic

Tuy OR cho phép Nhi diễn giải về ý nghĩa và vai trò của từng biến số đóng góp vào khả năng U ác tính, Nhi vẫn chưa đạt được mục tiêu thứ 2, đó là sử dụng mô hình logistic như 1 quy luật chẩn đoán trên lâm sàng. Trước khi làm điều này, Nhi cần phải cân chỉnh lại ngưỡng cắt cho xác suất dự báo (vì thứ mà mô hình này ước tính là 1 xác suất từ 0 đến 1 cho U ác tính, nhưng ta không biết phải đưa ra quyết định khi xác suất này lớn hơn 1 ngưỡng bằng bao nhiêu ? 0.5 chăng ? Hay cao hơn ? Hay thấp hơn ?)

Hình vẽ thứ nhất sau đây cho thấy một sự thật khá phiền toái, đó là có sự bất xứng giữa xác suất dự báo và quan sát thực tế. Tức là ngưỡng cắt sẽ thấp hơn 0.5 ; Nếu không được cân chỉnh, mô hình sẽ có nguy cơ cho ra kết quả âm tính giả và/hoặc âm tính giả.

library(classifierplots)


calibration_plot(trainset$Class, probs_predicted)+my_theme()

density_plot(trainset$Class, probs_predicted)+my_theme()+
  scale_fill_manual(values=c("#ffe314","#ff1499"))

## [1] "Generating score density plot"

lift_plot(trainset$Class, probs_predicted)+theme_bw()

Biểu đồ accuracy plot gợi ý ngưỡn cắt tối ưu nên là 30-40%, có thể là 35% thì khi đó tỉ lệ dự báo chính xác sẽ đạt tối ưu là 94% tới 95%

roc_plot(trainset$Class, probs_predicted)+theme_bw()

## [1] "Calculating AUC ..."
## [1] "(AUC) Sorting data ..."
## [1] "(AUC) Calculating ranks ..."
## [1] "AUC: 99.4268131868132"
## [1] "Bootstrapping ROC curves"
## [1] "Eval AUC"
## [1] "Producing ROC plot"

accuracy_plot(trainset$Class, probs_predicted)+theme_bw()

## Warning: Removed 40 rows containing missing values (geom_text).

10 Kiểm định hiệu năng của mô hình logistic như một quy luật thực tiễn

Bây giờ ta có thể ấn định ngưỡng cắt cho mô hình logistic là 0.35, nếu xác suất dự báo cao hơn ngưỡng này, kết quả sẽ được dán nhãn là U ác tính. Bằng confusion matrix, ta có thể kiểm tra hiệu năng của mô hình logistic này trên mẫu độc lập. Kết quả rất khả quan : mô hình này có thể được sử dụng như 1 quy luật chẩn đoán, và nó chỉ nhầm lẫn 1 lần duy nhất (1 trường hợp Âm tính giả). Độ chính xác hiệu chỉnh của mô hình là 0.992.

testpred=predict(m1,newdata=testset,type="response")

testset$Truth=ifelse(testset$Class==1,"Malignant","Benign")
testset$Pred=testpred
testset$Pred=ifelse(testset$Pred>=0.35,"Malignant","Benign")


caret::confusionMatrix(testset$Pred,reference=testset$Truth,mode="everything",positive="Malignant")

## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  Benign Malignant
##   Benign       119         1
##   Malignant      0        63
##                                           
##                Accuracy : 0.9945          
##                  95% CI : (0.9699, 0.9999)
##     No Information Rate : 0.6503          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.9879          
##  Mcnemar's Test P-Value : 1               
##                                           
##             Sensitivity : 0.9844          
##             Specificity : 1.0000          
##          Pos Pred Value : 1.0000          
##          Neg Pred Value : 0.9917          
##               Precision : 1.0000          
##                  Recall : 0.9844          
##                      F1 : 0.9921          
##              Prevalence : 0.3497          
##          Detection Rate : 0.3443          
##    Detection Prevalence : 0.3443          
##       Balanced Accuracy : 0.9922          
##                                           
##        'Positive' Class : Malignant       
## 

11 Kết luận

Suy diễn thống kê cho mô hình hồi quy giúp đưa ra câu trả lời cho 3 câu hỏi : thăm dò dữ liệu để phát hiện quy luật, dùng mô hình để chứng minh quy luật và giả thuyết, và dùng mô hình như 1 quy luật để tiên lượng hoặc phân loại. Với mô hình hồi quy, việc diễn giải dựa vào hệ số hồi quy của mô hình. Bootstrap là một giải pháp rất hữu ích để khẳng định về độ tin cậy của kết quả suy diễn hệ số hồi quy. Kiểm định độc lập mô hình là bắt buộc và rất quan trọng để có thể sử dụng mô hình như một quy tắc thực hành.

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