Hồi quy Bayes bằng TensorFlow

1 Giới thiệu

Câu chuyện về mô hình hồi quy logistic đã được lặp lại khá nhiều lần trong quá khứ, tuy nhiên trong bài thực hành ngắn này Nhi muốn trình bày 2 biến tấu của loại mô hình này với bản chất khác nhau, nhưng cả 2 đều liên quan đến Tensorflow.

2 Minh họa : phân biệt u vú lành tính hay ác tính

Nhi dùng lại thí dụ minh họa như trong bài Bayes hồi quy logistic lần trước 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 tiên lượng và diễn giải quan hệ giữa các biến.

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.2.1 --

## v ggplot2 3.0.0     v purrr   0.2.5
## v tibble  1.4.2     v dplyr   0.7.6
## v tidyr   0.8.1     v stringr 1.3.1
## v readr   1.1.1     v forcats 0.3.0

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

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

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

head(df)

3 Hàm glm()

Một mô hình logistic có thể được thực hiện đơn giản bằng hàm glm()

glm_mod<-glm(as.factor(Class)~.,data=df,family=stats::binomial())

summary(glm_mod)

## 
## Call:
## glm(formula = as.factor(Class) ~ ., family = stats::binomial(), 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4841  -0.1153  -0.0619   0.0222   2.4698  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -10.10394    1.17488  -8.600  < 2e-16 ***
## clumpthickness    0.53501    0.14202   3.767 0.000165 ***
## SizeUniformity   -0.00628    0.20908  -0.030 0.976039    
## ShapeUniformity   0.32271    0.23060   1.399 0.161688    
## Margin_adhesion   0.33064    0.12345   2.678 0.007400 ** 
## EpiCellSize       0.09663    0.15659   0.617 0.537159    
## Barenuclei        0.38303    0.09384   4.082 4.47e-05 ***
## BlandChromatin    0.44719    0.17138   2.609 0.009073 ** 
## NormalNucleoli    0.21303    0.11287   1.887 0.059115 .  
## Mitoses           0.53484    0.32877   1.627 0.103788    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 884.35  on 682  degrees of freedom
## Residual deviance: 102.89  on 673  degrees of freedom
## AIC: 122.89
## 
## Number of Fisher Scoring iterations: 8

library(broom)

glm_mod%>%tidy()%>%mutate(OR=exp(estimate),
                           UL=exp(estimate+1.96*std.error),
                           LL=exp(estimate-1.96*std.error))%>%
  select(estimate,OR,LL,UL,p.value)%>%round(.,4)

4 Logistic như mạng neuron 1 lớp

Tuy nhiên, chúng ta cũng có thể xem hồi quy logistic là một trường hợp đơn giản nhất của mạng Neuron (Neural network), chỉ có 1 lớp neuron duy nhất, được kích hoạt bằng 1 hàm sigmoid và nối với 1 neuron xuất kết quả.

Ta sẽ thử dựng neural network này với keras và TensorFlow backend, cấu trúc của mạng neuron gồm 1 lớp duy nhất để tiếp nhận 1 Intercept và 9 biến, kích hoạt bằng hàm sigmoid, xuất ra 1 neuron kết quả. Mô hình được huấn luyện với hàm loss là binary_crossentropy và tiêu chí kiểm định là accuracy

library(keras)

build_model <- function(){
  model_keras <- keras_model_sequential()%>% 
    layer_dense(
      units              = 1, 
      kernel_initializer = "uniform",
      activation         = "sigmoid", 
      input_shape        = 10)%>% 
    compile(
      optimizer = 'sgd',
      loss      = 'binary_crossentropy',
      metrics   = c('binary_accuracy')
    )
}

design_mat<-model.matrix(~., data=df[,-10])

head(design_mat,5)

##   (Intercept) clumpthickness SizeUniformity ShapeUniformity
## 1           1              5              1               1
## 2           1              5              4               4
## 3           1              3              1               1
## 4           1              6              8               8
## 5           1              4              1               1
##   Margin_adhesion EpiCellSize Barenuclei BlandChromatin NormalNucleoli
## 1               1           2          1              3              1
## 2               5           7         10              3              2
## 3               1           2          2              3              1
## 4               1           3          4              3              7
## 5               3           2          1              3              1
##   Mitoses
## 1       1
## 2       1
## 3       1
## 4       1
## 5       1

Với 20 batchs dữ liệu và 100 lượt train, mô hình nhanh chóng đạt trạng thái bão hòa với Accuracy = 0.99

model<-build_model()

model

## Model
## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #     
## ===========================================================================
## dense_1 (Dense)                  (None, 1)                     11          
## ===========================================================================
## Total params: 11
## Trainable params: 11
## Non-trainable params: 0
## ___________________________________________________________________________

fit_keras <- fit(
  object=model,
  x=design_mat, 
  y=df$Class,
  batch_size = 20,
  epochs= 100,
  validation_split = 0.3
)

fit_keras

## Trained on 478 samples, validated on 205 samples (batch_size=20, epochs=100)
## Final epoch (plot to see history):
##            val_loss: 0.108
## val_binary_accuracy: 0.9951
##                loss: 0.1709
##     binary_accuracy: 0.9477

Những giá trị trọng số cho lớp neuron trong mô hình:

get_weights(model)%>%.[[1]]->wdf

rownames(wdf)<-colnames(design_mat)

wdf

##                        [,1]
## (Intercept)     -1.67289507
## clumpthickness   0.08638906
## SizeUniformity   0.32727787
## ShapeUniformity  0.19069567
## Margin_adhesion  0.11432683
## EpiCellSize     -0.18441433
## Barenuclei       0.36628929
## BlandChromatin  -0.20194633
## NormalNucleoli   0.20442222
## Mitoses          0.03903888

5 Logistic theo Bayes

Tiếp theo, ta sẽ dùng công cụ “greta” đã giới thiệu ở bài trước để dựng 2 mô hình logistic theo phái Bayes, một dùng quy luật phân bố nhị thức (Binomial), một dùng họ phân phối Bernoulli.

Trong bài trước, Nhi đã giới thiệu về greta, thế hệ ngôn ngữ Bayes thứ tư và sử dụng Tensorflow. Thí dụ minh họa trong bài là một mô hình ANOVA, sử dụng design matrix làm dữ liệu đầu vào và likelihood là phân phối Gaussian. Cho bài toán mô hình logistic, chúng ta cùng đi ngược từ kết quả: Class là 1 biến nhị phân với giá trị 0/1. Để ước lượng Class bằng một mô hình hồi quy tuyến tính, ta phải tìm cách ánh xạ giá trị của biến kết quả Y trong mô hình trên khoảng 0-1 ; hai quy luật phân bố có thể áp dụng cho hàm likelihood của Y đó là binomial :

\[ Y \sim Binomial (N,p)\]

và Bernoulli

\[ Y \sim Bernoulli(p)\]

Trong đó, p lsinh ra từ hàm liên kết logit cho kết quả của một mô hình hồi quy tuyến tính. Mục tiêu của chúng ta là xác định 9 tham số beta và Intercept của mô hình.

Giả thuyết tiền định về beta và intercept, đó là chúng có phân bố chuẩn, với mu=0, và sigma = 10

5.1 Mô hình sử dụng Binomial (n,p)

Ta viết mô hình đầu tiên sử dụng likelihood là phân phối binomial(n,p)

library(greta)

## 
## Attaching package: 'greta'

## The following object is masked _by_ '.GlobalEnv':
## 
##     model

## The following objects are masked from 'package:stats':
## 
##     binomial, poisson

## The following objects are masked from 'package:base':
## 
##     %*%, backsolve, beta, colMeans, colSums, diag, eigen,
##     forwardsolve, gamma, rowMeans, rowSums, sweep, tapply

# create matrices to greta arrays
features <- as_data(df[,-10])

Y <- as_data(df$Class)

# create greta arrays for random variables
intercept <- normal(0, 10)
beta <- normal(0, 10, dim = 9L)

# logit-linear model
linear_predictor <- intercept + features %*% beta

p <- ilogit(linear_predictor)

# distribution (likelihood) over observed values
distribution(Y) <- binomial(nrow(df),p)

binom_mod <- greta::model(intercept,beta)

plot(binom_mod)

Sau đó ta sampling cho 10 chuỗi MCMC

draws <- mcmc(binom_mod,
              warmup = 1000,
              n_samples = 5000,             
              pb_update = 100)

Các chuỗi MCMC được thực hiện rất nhanh chóng, ta có thể xuất kết quả của mô hình

summary(draws)

## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                Mean      SD  Naive SE Time-series SE
## intercept -10.00825 0.21935 0.0031020      0.0038928
## beta[1,1]   0.12156 0.02870 0.0004058      0.0007045
## beta[2,1]   0.05680 0.04318 0.0006107      0.0030992
## beta[3,1]   0.02538 0.04411 0.0006239      0.0035118
## beta[4,1]   0.01947 0.02461 0.0003481      0.0014159
## beta[5,1]   0.02885 0.03119 0.0004411      0.0007063
## beta[6,1]   0.12853 0.02340 0.0003310      0.0007970
## beta[7,1]   0.03810 0.03494 0.0004942      0.0011135
## beta[8,1]   0.04888 0.02252 0.0003185      0.0008728
## beta[9,1]  -0.02137 0.02816 0.0003982      0.0006779
## 
## 2. Quantiles for each variable:
## 
##                2.5%        25%       50%       75%    97.5%
## intercept -10.44277 -10.155088 -10.00590 -9.860210 -9.59030
## beta[1,1]   0.06550   0.102452   0.12148  0.140905  0.17798
## beta[2,1]  -0.02692   0.027700   0.05702  0.085334  0.14238
## beta[3,1]  -0.06406  -0.004164   0.02385  0.056799  0.11109
## beta[4,1]  -0.02916   0.002400   0.01996  0.037015  0.06578
## beta[5,1]  -0.03158   0.008081   0.02859  0.049745  0.09144
## beta[6,1]   0.08371   0.112772   0.12851  0.143577  0.17471
## beta[7,1]  -0.03078   0.015994   0.03796  0.061085  0.10775
## beta[8,1]   0.00629   0.033133   0.04811  0.064587  0.09267
## beta[9,1]  -0.07776  -0.040327  -0.02124 -0.001916  0.03364

Cũng như phân phối hậu nghiệm của 9 Odds-ratios

post_df<-draws[[1]]%>%
  as.data.frame()%>%.[,c(1:10)]%>%
  mutate(iteration=c(1:5000))

colnames(post_df)=c("Intercept",colnames(df[,-10]),"Iteration")

post_df%>%gather(colnames(df[,-10]),key="Para",value="Value")%>%
  mutate(OR=exp(Value))%>%
  ggplot(aes(x=OR))+
  geom_density(alpha=0.5,aes(fill=Para,col=Para),show.legend = F)+
  theme_bw()+
  geom_vline(xintercept=1,col="blue4",linetype=2)+
  facet_wrap(~Para,ncol=3,scales="free")

library(ggridges)

## 
## Attaching package: 'ggridges'

## The following object is masked from 'package:ggplot2':
## 
##     scale_discrete_manual

library(viridis)

## Loading required package: viridisLite

post_df%>%gather(colnames(df[,-10]),key="Para",value="Value")%>%
  mutate(OR=exp(Value))%>%
  ggplot()+
  geom_density_ridges_gradient(aes(x=OR,y=reorder(Para,OR),fill=..x..),
                               scale=1,alpha=0.01,
                               show.legend =F)+
  scale_fill_viridis(option="C",direction = -1)+  
  labs(x="Odds-ratio",y="Betas")+
  geom_vline(xintercept = 1,col="blue",linetype=2)+
  theme_bw()

## Picking joint bandwidth of 0.00533

5.2 Mô hình sử dụng Bernoulli(p)

# create matrices to greta arrays
features <- as_data(df[,-10])

Y <- as_data(df$Class)

# create greta arrays for random variables
intercept <- normal(0, 10)
beta <- normal(0, 10, dim = 9L)

# logit-linear model
linear_predictor <- intercept + features %*% beta

prob <- ilogit(linear_predictor)

# distribution (likelihood) over observed values
distribution(Y) <- bernoulli(prob)

mod_bernoulli <- greta::model(intercept,beta)

plot(mod_bernoulli)

Ta làm tương tự như trên để có 10 chuỗi MCMC và xuất kết quả

drawsB <- mcmc(mod_bernoulli,
              warmup = 2000,
              n_samples = 5000,             
              pb_update = 100)

post_df2<-drawsB[[1]]%>%as.data.frame()%>%mutate(iteration=c(1:5000))

colnames(post_df2)=c("Intercept",colnames(df[,-10]),"Iteration")


post_df2%>%gather(colnames(df[,-10]),key="Para",value="Value")%>%
  mutate(OR=exp(Value))%>%
  ggplot(aes(x=OR))+
  geom_density(alpha=0.5,aes(fill=Para,col=Para),show.legend = F)+
  theme_bw()+
  geom_vline(xintercept=1,col="blue4",linetype=2)+
  facet_wrap(~Para,ncol=3,scales="free")

post_df2%>%gather(colnames(df[,-10]),key="Para",value="Value")%>%
  mutate(OR=exp(Value))%>%
  ggplot()+
  geom_density_ridges_gradient(aes(x=OR,y=reorder(Para,OR),fill=..x..),
                               scale=1,alpha=0.01,
                               show.legend =F)+
  scale_fill_viridis(option="D",direction = -1)+  
  labs(x="Odds-ratio",y="Betas")+
  geom_vline(xintercept = 1,col="blue",linetype=2)+
  theme_bw()

## Picking joint bandwidth of 0.0326

6 Tổng kết

Chúng ta vừa thực hiện 2 biến thể của mô hình hồi quy logistic, đó là mạng neuron 1 lớp và hồi quy Bayes. Tuy khác nhau về cơ chế, cả 2 đều dựa vào Tensorflow. Hy vọng bài viết giúp các bạn mua vui được vài phút.

Cuối tuần vui vẻ …