Hồi Qui Logistic Bayes Với Package Rethinking

1 Học Thống kê Bayes với Statistic Rethinking

Nếu thống kê theo trường phái Tần Số (Frequentist) đã khó hiểu thì trường phái Bayes còn trừu tượng hơn một mức vì các thông số không còn được trình bày dưới dạng các con số cụ thể (frequent: tần số) mà được khảo sát trên rất nhiều mẫu mô phỏng để xây dựng nên phân phối xác suất của các thông số để rồi lí luận và đánh giá xác xuất xảy ra của các biến số outcome.

Từ khi đọc được cuốn sách Statistic Rethinking của tác Giả Richard McElreath tôi thấy nhiều khái niệm được trình bày rõ hơn. Cùng với cuốn sách này, tác giả còn sử dụng Package Rethinking để phân tích dữ liệu, suy luận thống kê theo Bayes. Rethinking cũng giống như rstanarm, nó viết code cho stan chạy phân tích để đơn giản hóa các câu lệnh cho người dùng bớt nhức đầu.

Để có sự so sánh với các Package khác như glm và rstanarm mà Bs Nhi đã viết trong bài trước, tôi sử dụng lại dataset này để phân tích với Rethinking và so sánh với hàm glm, rstanarm.

2 Chuẩn bị dữ liệu

Đó là bộ dữ liệu Biopsy trong package MASS mà Bs Nhi đã sử dụng. Bộ dữ liệu này có 683 row và 10 collumn. Các biến số của nó là clumpthickness, SizeUniformity, ShapeUniformity, Margin_adhesion, EpiCellSize, Barenuclei, BlandChromatin, NormalNucleoli, Mitoses và Class. Đó là các biến số mô tả những đặc tính mô học, tế bào học của mẫu bệnh phẩm tuyến giáp và Class là kết luận Malignant (ác tính) hay Benign (lành tính) về các tổ chức bệnh phẩm.

library(dplyr)
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)
head(df)

## # A tibble: 6 x 10
##   clumpthickness SizeUniformity ShapeUniformity Margin_adhesion
##            <int>          <int>           <int>           <int>
## 1              5              1               1               1
## 2              5              4               4               5
## 3              3              1               1               1
## 4              6              8               8               1
## 5              4              1               1               3
## 6              8             10              10               8
## # ... with 6 more variables: EpiCellSize <int>, Barenuclei <int>,
## #   BlandChromatin <int>, NormalNucleoli <int>, Mitoses <int>,
## #   Class <fctr>

3 Phân chia dataset thành hai nhóm nhỏ để train và test model

Chúng ta phân chia dataset trên thành 2 dataset nhỏ:

trainset với 500 observations, để training model; testset với 183 observations, để kiểm định model trong bộ dữ liệu hoàn toàn mới

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

## [1] 500  10

dim(testset)

## [1] 183  10

4 Dựng mô hình glm và đánh giá khả năng classification của mô hình này

Chúng ta dựng mô hình với glm để có cơ sở so sánh với model mới.

fml="Class~clumpthickness+ ShapeUniformity+ Margin_adhesion+ Barenuclei+ BlandChromatin"
fitglm=glm(fml, data = trainset, family = binomial(link = "logit"))
summary(fitglm)

## 
## Call:
## glm(formula = fml, family = binomial(link = "logit"), data = trainset)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -3.04325  -0.14356  -0.06977   0.03671   2.34910  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -9.07134    1.06807  -8.493  < 2e-16 ***
## clumpthickness   0.58294    0.14131   4.125  3.7e-05 ***
## ShapeUniformity  0.48237    0.15835   3.046 0.002317 ** 
## Margin_adhesion  0.31521    0.11347   2.778 0.005469 ** 
## Barenuclei       0.34662    0.09697   3.574 0.000351 ***
## BlandChromatin   0.44241    0.17664   2.505 0.012259 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 647.45  on 499  degrees of freedom
## Residual deviance:  97.51  on 494  degrees of freedom
## AIC: 109.51
## 
## Number of Fisher Scoring iterations: 7

# Classification table
pred_glm<-predict(fitglm, testset,type="response")
pred_class <- ifelse(pred_glm > 0.5, 1, 0)
table(testset$Class, pred_class)

##    pred_class
##       0   1
##   0 119   0
##   1   4  60

Mô hình GLM đã phân loại đúng 119 + 60 = 179/183 quan sát, tỉ lệ chính xác là 0.9781.

5 Dựng mô hình với Rethinking

library(rstan)
library(rethinking)
m10.20 <- map( 
  alist(
    Class ~ dbinom(1, p) ,
    logit(p) <- a + b1*clumpthickness +
      b2*ShapeUniformity +
      b3*Margin_adhesion +
      b4*Barenuclei +
      b5*BlandChromatin ,
    a ~ dnorm(0,1) ,
    b1 ~ dnorm(0,1),
    b2 ~ dnorm(0,1),
    b3 ~ dnorm(0,1),
    b4 ~ dnorm(0,1),
    b5 ~ dnorm(0,1)
  ) ,
  data=data.frame(trainset) )
m10.20.stan <- map2stan(m10.20, data = data.frame(trainset), chains=2,iter=1e4 , warmup=1000)

## In file included from C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
##                  from C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file504c5a4c753d.cpp:8:
## C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
##  #  define BOOST_NO_CXX11_RVALUE_REFERENCES
##  ^
## <command-line>:0:0: note: this is the location of the previous definition
## In file included from C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/multi_array/base.hpp:28:0,
##                  from C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/multi_array.hpp:21,
##                  from C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/numeric/odeint/util/multi_array_adaption.hpp:29,
##                  from C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/numeric/odeint.hpp:61,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:13,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/arr.hpp:38,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/mat.hpp:298,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:12,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file504c5a4c753d.cpp:8:
## C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp: In static member function 'static void boost::multi_array_concepts::detail::idgen_helper<N>::call(Array&, const IdxGen&, Call_Type)':
## C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:42:43: warning: typedef 'index_range' locally defined but not used [-Wunused-local-typedefs]
##        typedef typename Array::index_range index_range;
##                                            ^
## C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:43:37: warning: typedef 'index' locally defined but not used [-Wunused-local-typedefs]
##        typedef typename Array::index index;
##                                      ^
## C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp: In static member function 'static void boost::multi_array_concepts::detail::idgen_helper<0u>::call(Array&, const IdxGen&, Call_Type)':
## C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:53:43: warning: typedef 'index_range' locally defined but not used [-Wunused-local-typedefs]
##        typedef typename Array::index_range index_range;
##                                            ^
## C:/Users/Home/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:54:37: warning: typedef 'index' locally defined but not used [-Wunused-local-typedefs]
##        typedef typename Array::index index;
##                                      ^
## In file included from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:42:0,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file504c5a4c753d.cpp:8:
## C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp: At global scope:
## C:/Users/Home/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function]
##      static void set_zero_all_adjoints() {
##                  ^
## 
## SAMPLING FOR MODEL 'Class ~ dbinom(1, p)' NOW (CHAIN 1).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 10000 [  0%]  (Warmup)
## Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Iteration: 1001 / 10000 [ 10%]  (Sampling)
## Iteration: 2000 / 10000 [ 20%]  (Sampling)
## Iteration: 3000 / 10000 [ 30%]  (Sampling)
## Iteration: 4000 / 10000 [ 40%]  (Sampling)
## Iteration: 5000 / 10000 [ 50%]  (Sampling)
## Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Iteration: 10000 / 10000 [100%]  (Sampling)
## 
##  Elapsed Time: 4.321 seconds (Warm-up)
##                36.579 seconds (Sampling)
##                40.9 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'Class ~ dbinom(1, p)' NOW (CHAIN 2).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 10000 [  0%]  (Warmup)
## Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Iteration: 1001 / 10000 [ 10%]  (Sampling)
## Iteration: 2000 / 10000 [ 20%]  (Sampling)
## Iteration: 3000 / 10000 [ 30%]  (Sampling)
## Iteration: 4000 / 10000 [ 40%]  (Sampling)
## Iteration: 5000 / 10000 [ 50%]  (Sampling)
## Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Iteration: 10000 / 10000 [100%]  (Sampling)
## 
##  Elapsed Time: 4.144 seconds (Warm-up)
##                45.001 seconds (Sampling)
##                49.145 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'Class ~ dbinom(1, p)' NOW (CHAIN 1).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## Iteration: 1 / 1 [100%]  (Sampling)
## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.001 seconds (Sampling)
##                0.001 seconds (Total)
## 
## [ 1800 / 18000 ]
[ 3600 / 18000 ]
[ 5400 / 18000 ]
[ 7200 / 18000 ]
[ 9000 / 18000 ]
[ 10800 / 18000 ]
[ 12600 / 18000 ]
[ 14400 / 18000 ]
[ 16200 / 18000 ]
[ 18000 / 18000 ]

Tương tự với GLM, trong Rethinking ta vẫn chọn Class làm biến số outcome, và các biến còn lại là những predictors. Cách viết code của mô hình này được giải thích rất chi tiết trong cuốn sách Statistic Rethinking.

Class ~ dbinom(1, p): Biến số outcome Class phân phối binomial vì nhân hai giá trị 0 hoặc 1 với xác suất p cho mỗi quan sát (observation).

Nếu p rất nhỏ, thì Class có thể có phân phối Poisson, lúc đó chúng ta viết: Class ~ dpois(lambda)

Tuy nhiên, thử mô hình này tôi thấy variance của nó (WAIC) lớn hơn nhiều so với mô hình với phân phối binomial. Điều này cũng phù hợp với dataset mà chúng ta đang phân tích, có p của biến số Class (tỉ lệ Malignant) không quá nhỏ.

Stan code của mô hình trên được cho biết với hàm stancode() như sau:

stancode(m10.20.stan)

## data{
##     int<lower=1> N;
##     int Class[N];
##     int clumpthickness[N];
##     int ShapeUniformity[N];
##     int Margin_adhesion[N];
##     int Barenuclei[N];
##     int BlandChromatin[N];
## }
## parameters{
##     real a;
##     real b1;
##     real b2;
##     real b3;
##     real b4;
##     real b5;
## }
## model{
##     vector[N] p;
##     b5 ~ normal( 0 , 1 );
##     b4 ~ normal( 0 , 1 );
##     b3 ~ normal( 0 , 1 );
##     b2 ~ normal( 0 , 1 );
##     b1 ~ normal( 0 , 1 );
##     a ~ normal( 0 , 1 );
##     for ( i in 1:N ) {
##         p[i] = a + b1 * clumpthickness[i] + b2 * ShapeUniformity[i] + b3 * Margin_adhesion[i] +      b4 * Barenuclei[i] + b5 * BlandChromatin[i];
##     }
##     Class ~ binomial_logit( 1 , p );
## }
## generated quantities{
##     vector[N] p;
##     real dev;
##     dev = 0;
##     for ( i in 1:N ) {
##         p[i] = a + b1 * clumpthickness[i] + b2 * ShapeUniformity[i] + b3 * Margin_adhesion[i] +      b4 * Barenuclei[i] + b5 * BlandChromatin[i];
##     }
##     dev = dev + (-2)*binomial_logit_lpmf( Class | 1 , p );
## }

6 Khai thác mô hình trên với trainset

show(m10.20.stan)

## map2stan model fit
## 18000 samples from 2 chains
## 
## Formula:
## Class ~ dbinom(1, p)
## logit(p) <- a + b1 * clumpthickness + b2 * ShapeUniformity + 
##     b3 * Margin_adhesion + b4 * Barenuclei + b5 * BlandChromatin
## a ~ dnorm(0, 1)
## b1 ~ dnorm(0, 1)
## b2 ~ dnorm(0, 1)
## b3 ~ dnorm(0, 1)
## b4 ~ dnorm(0, 1)
## b5 ~ dnorm(0, 1)
## 
## Log-likelihood at expected values: -55.58 
## Deviance: 111.15 
## DIC: 122.25 
## Effective number of parameters (pD): 5.55 
## 
## WAIC (SE): 121.04 (16.3)
## pWAIC: 4.08

precis(m10.20.stan)

##     Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a  -5.94   0.48      -6.68      -5.13 12555    1
## b1  0.27   0.09       0.13       0.42 12436    1
## b2  0.47   0.12       0.27       0.66 15087    1
## b3  0.22   0.09       0.07       0.36 14792    1
## b4  0.33   0.08       0.20       0.46 16481    1
## b5  0.21   0.12       0.01       0.41 14331    1

precise(m10.20.stan) và show(), tương tự hàm summary() trong GLM, cho ta bảng kết quả tóm tắt về phân phối của các thông số (hệ số của các predictors trong mô hình trên).

n_eff là number of observations that is effective. Rhat là R^

Hai chỉ số trên cho thấy các thông số đều có ý nghĩa trong mô hình. Điều này tương tự như giá trị p trong thống kê frequentist.

Để hiểu ý nghĩa của những chỉ số này, các bạn có thể đọc sách Statistic Rethinking để có được những sự giãi thích chi tiết.

Một trong những điều quan tâm là tương quan giữa các thông số của các biến số predictors trong mô hình. Nếu chúng tương quan mạnh với nhau thì mô hình trên có độ chính xác không cao vì variance của mô hình lớn hơn.

Hàm pairs(), cho ta biết tương quan giữa các thông số trong mô hình trên.

pairs(m10.20.stan)

Biểu đồ trên cho ta thấy các thông số không có tình trạng tương quan mạnh với nhau.

Phân phối của các thông số của mô hình sau khi mô phỏng với 20.000 mẫu qua 2 chuỗi MCMC (iter=1e4)

# Plots of posterior distributions
post <- extract.samples( m10.20.stan )
par(mfrow=c(3,2))
plot(post$a, type="l", col= "pink")
plot(post$b1, type="l",col= "yellow")
plot(post$b2, type="l", col= "blue")
plot(post$b3, type="l", col= "green")
plot(post$b4, type="l",col= "orange")
plot(post$b5, type="l", col= "red")

Biểu đồ tương đối đều, không có biến động bất thường, nói lên phân phối xác suất hậu nghiệm của các thông số gần với phân phối chuẩn, variance khá đồng nhất. Đó là điều chúng ta mong muốn đối với một mô hình trong một bộ dataset được phân tích.

7 OR của mô hình

Hàm precis(m10.20.stan) cho ta bảng tóm tắt về phân phôi xác suất của các thông số từ a, b1 đến b5 tương ứng với intercept và các predictors trong mô hình trên.

Vì Logistic xây dựng hàm

logit(p) = a + b1* clumpthickness + b2ShapeUniformity + b3Margin_adhesion + b4Barenuclei + b5BlandChromatin

nên OR = exp() của các thông số (e mũ logit(p))

Từ đó tôi xây dựng bộ dữ liệu của OR thủ công như sau:

Mean <- c(-5.942, 0.272, 0.47, 0.218, 0.33, 0.21)
lower <- c(-6.6778, 0.1278, 0.2668, 0.0668, 0.1986, 0.01)
upper <- c(-5.131, 0.422, 0.661, 0.362, 0.457,0.41)
OR <- exp(Mean)
lowerOR <- exp(lower)
upperOR <- exp(upper)
ordf <- data.frame(OR, lowerOR, upperOR)
ordf$predictors  <- c("a-intercept", "b1-Clumthickness", "b2-ShapeUniformity", 
                      "b3-MarginAdhesion", "b4-Barenuclei", "b5-BlandChromatin")
head(ordf)

##            OR     lowerOR     upperOR         predictors
## 1 0.002626771 0.001258544 0.005910647        a-intercept
## 2 1.312587001 1.136325715 1.525008525   b1-Clumthickness
## 3 1.599994193 1.305779264 1.936728094 b2-ShapeUniformity
## 4 1.243587068 1.069081640 1.436198942  b3-MarginAdhesion
## 5 1.390968128 1.219693991 1.579328884      b4-Barenuclei
## 6 1.233678060 1.010050167 1.506817785  b5-BlandChromatin

Biểu đồ trình bày OR của các biến số predictors

library(ggplot2)
ordf%>%ggplot(aes(x=reorder(predictors,OR),color=reorder(predictors,OR)))+
  geom_errorbar(aes(ymin=lowerOR, ymax=upperOR),width=0.5,size=1,show.legend = F)+
  geom_point(aes(y=OR),shape=21,size=5,fill="white",stroke=1.5,show.legend = F)+
  ggtitle("Marginalized effect for ORs") +
  coord_flip()+
  geom_hline(yintercept = 1,size=1,color="blue",linetype=2)+
  scale_x_discrete("Predictors")+
  scale_y_continuous("Odds-Ratio")

OR và khoảng 5.5%-94.5%, theo mặc định của hàm precis(). Tại sao không là 5% - 95% ? Câu trả lời đơn giản là không muốn chúng ta mãi ám ảnh vào thống kê frequentist với khoảng tin cậy 95% truyền thống.

OR nói lên ý nghĩa của các biến số predictors và tầm quan trọng tương đối giữa chúng với nhau, được biểu thị qua giá trị của mỗi OR tương ứng với mỗi predictors.

Thứ tự của mỗi OR tương ứng khác nhau so với khi phân tích với package rstanarm() trong bài viết trước của Bs Nhi.

Câu trả lời thỏa đáng cho sự khác biệt này tôi vẫn chưa có được. Có thể các prior khác nhau trong mô hình trên so với mô hình từ rstanarm() mà Bs Nhi xây dựng, đã tạo ra sự khác biệt này.

8 Đánh giá khả năng phân loại của mô hình trong testset.

Đến đây chúng ta đã có được cái nhìn về khả năng của mô hình trên trong trainset, nơi training mô hình này.

Để xem mô hình này chính xác như thế nào với một dataset hoàn toàn mới, chúng ta áp dụng vào testset sau đó phân tích các chỉ số performance của nó.

Hàm link() trong package Rethinking cho ta phân phối hậu nghiệm của xác suất p của các quan sát (observations) thông qua mô hình stan ở trên.

p <- link(m10.20.stan, data = testset)

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

str(p)

##  num [1:1000, 1:183] 0.0365 0.0415 0.0454 0.0394 0.0315 ...

p_mean <- apply(p,2,mean)
p_HPDI <- apply(p,2,HPDI,prob=0.9)
pdf <- rbind(p_mean,p_HPDI)
head(p_mean)

## [1] 0.04260364 0.06080281 0.28674465 0.03805447 0.33708011 0.98954545

Đây là kết quả của hàm link(), cho ra một matrix xác suất p. 1:1000 cho biết có từ 1 đến 1000 mẫu được mô phỏng qua chuỗi MCMC, 1:183 cho biết từ 1 đến 183 quan sát (observations) của testset.

Mỗi một quan sát có một xác suất để nó được phân loại là 0 (Benign) hoặc 1 (Malignant). Một quan sát được mô phỏng lại 1000 lần, các giá trị này tạo ra một phân phối, đây là giá trị cốt lõi của Bayes.

Từ đó tính ra được xác suất p trung bình cho mỗi quan sát. Dựa trên xác suất trung bình này chúng ta tiên đoán (predict) quan sát đó là 0 (nếu p_mean < 0.5, tức là Benign) hoặc 1 (nếu p_mean >= 0.5, tức là Malignant).

Với mỗi quan sát, không chỉ có giá 1 trị p trung bình mà còn có khoảng mật độ xác suất cao (HPDI) với mặc định là 90%.

Để xem mô hình trên phân loại “Benign” hoặc “Malignant” chính xác thế nào chúng ta lập bảng với predictions và Class phân loại thực trong testset.

predictions <- ifelse(p_mean < 0.5, 0, 1)
table(predictions)

## predictions
##   0   1 
## 122  61

table(testset$Class)

## 
##   0   1 
## 119  64

table(predictions, testset$Class)

##            
## predictions   0   1
##           0 119   3
##           1   0  61

Trên bộ dữ liệu testset, mô hình rethinking đã phân loại đúng 119 + 61 = 180/183 quan sát, tỉ lệ chính xác là 0.9836, cao hơn so với mô hình glm (0.9781).

Để biết mức độ chính xác, sensitivity và specificity của mô hình trong testset ta sử dụng hàm caret để lập confusion matrix của bảng kết quả phân loại.

testset$Truth = ifelse(testset$Class==1,"Malignant","Benign")
testset$Pred = predictions
testset$Pred = ifelse(testset$Pred == 1,"Malignant","Benign")

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

## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  Benign Malignant
##   Benign       119         3
##   Malignant      0        61
##                                           
##                Accuracy : 0.9836          
##                  95% CI : (0.9528, 0.9966)
##     No Information Rate : 0.6503          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.9636          
##  Mcnemar's Test P-Value : 0.2482          
##                                           
##             Sensitivity : 0.9531          
##             Specificity : 1.0000          
##          Pos Pred Value : 1.0000          
##          Neg Pred Value : 0.9754          
##               Precision : 1.0000          
##                  Recall : 0.9531          
##                      F1 : 0.9760          
##              Prevalence : 0.3497          
##          Detection Rate : 0.3333          
##    Detection Prevalence : 0.3333          
##       Balanced Accuracy : 0.9766          
##                                           
##        'Positive' Class : Malignant       
## 

Kết quả cho thấy mức độ chính xác cao, sensitivity và specificity tương đương với kết quả từ mô hình được xây dựng với hàm rstanarm() mà Bs Nhi đã trình bày trong bài trước.

Tóm lại, với rethinking() chúng ta xây dựng mô hình phân loại logistic regression cho ra kết quả có khác biệt ít nhiều với mô hình từ hàm rstanarm() về tầm quan trọng tương đối của các biến số predictors (thứ tự quan trọng của các biến số). Tuy nhiên, performance của mô hình với rethinking() thì cũng tương tự với rstanarm().

Cám ơn các bạn đã đọc bài viết. Mong nhận được những ý kiến, phân tích để làm rõ hơn những vấn đề của chủ đề này.