Regularisation mô hình hồi quy Bayes

1 Regularisation cho mô hình Bayes

Mô hình tuyến tính tổng quát (GLM) là một công cụ phổ quát và linh hoạt để suy diễn thống kê trong nghiên cứu khoa học. Khi dựng một mô hình, chúng ta đang cố gắng đơn giản hay thu gọn một hệ thống phức tạp (đa chiều) trong thế giới tự nhiên rộng lớn thành một hệ thống đơn giản vừa đủ cho phép tiên lượng giá trị biến kết quả (outcome) dựa vào một số ít những biến số (predictor) tiêu biểu nhất.

Không gian dữ liệu đơn giản này thường được nhà nghiên cứu chủ động tạo ra bằng thiết kế thí nghiệm, nhất là nghiên cứu lâm sàng. Đôi khi, người ta không hài lòng với mô hình quá đơn giản và lại tìm cách mở rộng nó thành mô hình đa biến, hoặc mixed model. Theo thời gian, tiến bộ về công nghệ khiến dữ liệu ngày càng phong phú và dễ thu thập hơn, do đó số lượng biến trong dữ liệu trở nên to lớn, thậm chí lớn hơn cả kích thước mẫu khảo sát. Lúc này, bài toán chọn lọc biến số, giản lược nội dung mô hình trở nên quan trọng. Việc xác định một mô hình tối giản không chỉ cần thiết cho mục tiêu suy diễn mà còn có ảnh hưởng tích cực cho độ chính xác vì giảm nguy cơ mô hình bị “overfit”.

Tương tự những phương pháp hiệu chỉnh (Regularisation) giá trị tham số hồi quy như ridge, lasso, elasticnet bên phái « frequentist », khi làm hồi quy Bayes chúng ta cũng mong muốn rút gọn nội dung mô hình bằng cách đẩy những tham số có hiệu ứng cực thấp trong mô hình về zero, nhưng bảo toàn giá trị những tham số có vai trò đóng góp quan trọng cho kết quả tiên lượng.

2 Prior Piironen-Vehtari (2017)

Năm 2017, hai tác giả Piironen và Vehtari đã tìm ra giải pháp cho câu hỏi trên, họ đề xuất 1 prior có hành vi giống như trên. Ta gọi nó là prior « hình móng ngựa » (horseshoe prior) vì đồ thị hàm mật độ xác suất của nó có dạng chữ U, tập trung về gần 0 hoặc 1.

Như ta biết, tham số hồi quy beta trong các mô hình tuyến tính Bayes thường được áp dụng một prior yếu như Gaussian với µ=0 và sigma rộng (thí dụ sigma=100) :

\[\beta_j \sim{} N(0, \sigma^2)\]

Prior hình móng ngựa cổ điển được xác định bằng một phân bố Gaussian có µ cũng =0, nhưng sigma được phân tích thành 2 tham số bộ phận khác là Tau và Lambda như sau :

\[\beta_j \sim{} N(0, \tau^2\lambda_j^2)\]

Tham số Tau và Lambda đều cùng được mô tả bằng phân phối Cauchy(0,1),

\[\tau \sim{} Cauchy(0,1)\]

\[\lambda_j \sim{} Cauchy(0,1)\] Trong đó j nhận giá trị từ 1,2,3… đến nX là số lượng predictor trong design matrix của mô hình

Như vậy, vai trò của tham số Tau là xác định cho một tham số scale parameter tổng quát (khoảng nằm giữa của phân bố) có khuynh hướng đẩy giá trị beta về zero qua một trọng số. Trong khi đó, Lambda đặc trưng cho phần đuôi của phân bố, cho phép bảo tồn giá trị của một vài tham số beta.

Piironen và Vehtari đã cải tiến prior móng ngựa cơ bản này thành một dạng mới, trong đó Tau vẫn được xác định bằng phân bố Cauchy(0,1) , tuy nhiên tham số Lambda được xác định phụ thuộc vào Tau và một hằng số c (gọi là slab) :

\[\beta_j \sim{} N(0, \tau^2\tilde{\lambda_j}^2)\]

\[\tau \sim{} cauchy(0,1)\]

\[\tilde{\lambda_j}^2 =\frac{c^2\lambda_j^2}{c^2+\tau^2\lambda_j^2}\]

\[\lambda_j \sim{} Cauchy(0,1)\]

Với những biến số (predictor) có hiệu ứng mạnh, ta có :

\[\tau^2\lambda_j^2 \gg c^2\]

Khi đó, prior sẽ xấp xỉ :

\[N(0,c^2)\]

Khi đó beta sẽ không bị phạt nhiều, trái lại những biến số có hiệu ứng yếu (bao gồm cả trường hợp có quá ít dữ liệu) , ta có:

\[\tau^2\lambda_j^2 \ll c^2\] Lúc này gía trị beta sẽ bị đầy về gần zero theo quy luật chung của prior móng ngựa

3 Thí dụ minh họa

Trong thí dụ minh họa này, Nhi sử dụng bộ dữ liệu « Body Fat » của Roger W. Johnson và Carleton College, có thể tải từ:

https://ww2.amstat.org/publications/jse/v4n1/datasets.johnson.html

Mục tiêu của bài toán là xây dựng một mô hình hồi quy tuyến tính nhằm ước tính lượng mỡ trong cơ thể (biến Brozek Body Fat) dựa vào Tuổi, trọng lượng, chiều cao, và các chỉ số nhân trắc (bụng, ngực, cổ tay,…). Bài toán cũng có thể được đặt ra theo hướng diễn dịch với mục tiêu khảo sát mối tương quan bộ phận giữa mỗi chỉ số nhân trắc và giá trị Body_fat.

library(tidyverse)

df=read.table("https://ww2.amstat.org/publications/jse/datasets/fat.dat.txt")

df= df[,c(2, 5:7, 10:19)]

names(df)=c("brozek_fat","age","weight","height","neck","chest","abdomen","hip","thigh","knee","ankle","biceps","forearm","wrist")
df= df[-42,]

head(df)%>%knitr::kable()

brozek_fat	age	weight	height	neck	chest	abdomen	hip	thigh	knee	ankle	biceps	forearm	wrist
12.6	23	154.25	67.75	36.2	93.1	85.2	94.5	59.0	37.3	21.9	32.0	27.4	17.1
6.9	22	173.25	72.25	38.5	93.6	83.0	98.7	58.7	37.3	23.4	30.5	28.9	18.2
24.6	22	154.00	66.25	34.0	95.8	87.9	99.2	59.6	38.9	24.0	28.8	25.2	16.6
10.9	26	184.75	72.25	37.4	101.8	86.4	101.2	60.1	37.3	22.8	32.4	29.4	18.2
27.8	24	184.25	71.25	34.4	97.3	100.0	101.9	63.2	42.2	24.0	32.2	27.7	17.7
20.6	24	210.25	74.75	39.0	104.5	94.4	107.8	66.0	42.0	25.6	35.7	30.6	18.8

Có đến 13 predictors trong dữ liệu.

library(corrplot)
library(RColorBrewer)

cor_mat=as.matrix(cor(method="pearson",as.matrix(df)))

cor_mat%>%corrplot(.,order="hclust",type="lower",method="color",tl.col="black", tl.srt=45,tl.cex=0.5,col=rev(brewer.pal(n=10, name="RdBu")))

Nhi sẽ dựng mô hình hồi quy tuyến tính có kèm hiệu chỉnh (regularisation) tham số hồi quy cho mỗi predictor để phân lập những biến có hiệu ứng mạnh nhất (quan trọng nhất ) và những biến có ảnh hưởng yếu hoặc không có vai trò đóng góp đáng kể.

8 phương pháp sẽ được so sánh: Bên phái frequentist, đầu tiên là một mô hình hồi quy tuyến tính đa biến không có hiệu chỉnh gì cả, được dựng bằng hàm GLM, sau đó Nhi sẽ lần lượt thử : hồi quy Ridge, hồi quy Lasso và Elastic net bằng package glmnet (Friedman, Trevor Hastie và Rob Tibshirani, 2013). Bên phái Bayes, Nhi sẽ áp dụng phương pháp BMA với package BMS, và prior hình móng ngựa cải tiến của Piironen (2017) được thực hiện bằng 3 cách khác nhau: code thủ công bằng STAN, cũng như thông qua giao thức rstanarm và brms.

4 Mô hình GLM không áp dụng regularisation

# GLM

fit_glm=glm(data=df,brozek_fat~.)

summary(fit_glm)

## 
## Call:
## glm(formula = brozek_fat ~ ., data = df)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -10.2537   -2.5759   -0.0991    2.8863    9.3718  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.80071   20.60403  -0.864   0.3885    
## age           0.05684    0.03003   1.893   0.0596 .  
## weight       -0.08595    0.05747  -1.496   0.1361    
## height       -0.03734    0.16571  -0.225   0.8219    
## neck         -0.43027    0.21897  -1.965   0.0506 .  
## chest        -0.01844    0.09575  -0.193   0.8474    
## abdomen       0.89013    0.08378  10.624   <2e-16 ***
## hip          -0.19589    0.13605  -1.440   0.1512    
## thigh         0.23640    0.13596   1.739   0.0834 .  
## knee         -0.02129    0.22994  -0.093   0.9263    
## ankle         0.16723    0.20643   0.810   0.4187    
## biceps        0.15655    0.15999   0.978   0.3288    
## forearm       0.42942    0.18491   2.322   0.0211 *  
## wrist        -1.47440    0.49664  -2.969   0.0033 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 15.96847)
## 
##     Null deviance: 14915.5  on 250  degrees of freedom
## Residual deviance:  3784.5  on 237  degrees of freedom
## AIC: 1423.3
## 
## Number of Fisher Scoring iterations: 2

5 Bộ ba Lasso-Ridge-Elasticnet

Đầu tiên là mô hình LASSO

# Lasso

library(glmnet)

lasso = cv.glmnet(x=as.matrix(df[,-1]),y=df$brozek_fat,alpha=1,nfolds=5)

coef(lasso)

## 14 x 1 sparse Matrix of class "dgCMatrix"
##                       1
## (Intercept) -9.09239847
## age          0.02400386
## weight       .         
## height      -0.24340856
## neck         .         
## chest        .         
## abdomen      0.57021332
## hip          .         
## thigh        .         
## knee         .         
## ankle        .         
## biceps       .         
## forearm      .         
## wrist       -0.47897402

Sau đó là mô hình Ridge

ridge = cv.glmnet(x=as.matrix(df[,-1]),y=df$brozek_fat,alpha=0,nfolds=5)

coef(ridge)

## 14 x 1 sparse Matrix of class "dgCMatrix"
##                       1
## (Intercept) -6.15352703
## age          0.09357610
## weight       0.01731019
## height      -0.37028608
## neck        -0.12669803
## chest        0.14384149
## abdomen      0.27866578
## hip          0.08675680
## thigh        0.13744500
## knee         0.10353008
## ankle       -0.07357488
## biceps       0.07180562
## forearm      0.16700123
## wrist       -0.98190998

Cuối cùng là Elastic net

enet= cv.glmnet(x=as.matrix(df[,-1]),y=df$brozek_fat,nfolds=5)

coef(enet)

## 14 x 1 sparse Matrix of class "dgCMatrix"
##                       1
## (Intercept) -6.53342854
## age          0.02997906
## weight       .         
## height      -0.24850948
## neck         .         
## chest        .         
## abdomen      0.58853559
## hip          .         
## thigh        .         
## knee         .         
## ankle        .         
## biceps       .         
## forearm      .         
## wrist       -0.70732370

Như ta thấy, mô hình Lasso và elasticnet có khuynh hướng triệt tiêu hẳn những biến số được cho là yếu khỏi mô hình, trong khi ridge regression chỉ giảm gía trị tham số hồi quy cho những biến này xuống thấp nhất có thể.

6 Phương pháp BMA

Trong khi đó, phương pháp Bayesian Model averaging kiểm tra tất cả những tổ hợp có thể giữa 13 biến và xác định những mô hình tối ưu nhất.

library(BMS)

bmsmod= bms(df, burn = 1000, 
            iter = 10000, 
            g = "EBL",
            mprior = "uniform",
            mcmc="bdn", 
            user.int = F)

coef(bmsmod, std.coefs = T, order.by.pip = T, include.constant = T)

##                PIP    Post Mean    Post SD Cond.Pos.Sign Idx
## abdomen     1.0000  1.255317146 0.09604764    1.00000000   6
## weight      0.9146 -0.434024060 0.18300505    0.00000000   2
## wrist       0.7624 -0.118947512 0.08372544    0.00000000  13
## forearm     0.5851  0.066038463 0.06619393    1.00000000  12
## neck        0.3310 -0.041071144 0.07142256    0.00000000   4
## thigh       0.2303  0.026584495 0.06242852    0.99956578   8
## biceps      0.2165  0.021031372 0.04944985    1.00000000  11
## age         0.1940  0.012537607 0.03462482    0.97010309   1
## hip         0.1910 -0.034141666 0.09343202    0.00000000   7
## height      0.1115 -0.004677901 0.02397510    0.07533632   3
## knee        0.0861  0.003214869 0.02276074    0.94308943   9
## chest       0.0834 -0.001058647 0.03027382    0.51318945   5
## ankle       0.0737  0.002293030 0.01514326    0.96200814  10
## (Intercept) 1.0000 -3.696741238         NA            NA   0

image(bmsmod,col=c("#199ac1","#e30c37"))

summary(bmsmod)

## Mean no. regressors               Draws             Burnins 
##            "4.7796"             "10000"              "1000" 
##                Time  No. models visited      Modelspace 2^K 
##    "0.8600001 secs"              "3360"              "8192" 
##           % visited         % Topmodels            Corr PMP 
##                "41"                "98"            "0.9854" 
##            No. Obs.         Model Prior             g-Prior 
##               "251"     "uniform / 6.5"               "EBL" 
##     Shrinkage-Stats 
##         "Av=0.9929"

7 Code mô hình thủ công với STAN

Trong bài báo gốc, 2 tác giả đã cung cấp 2 bộ STAN code viết sẵn cho mô hình hồi quy tuyến tính Gaussian (như thí dụ này) và logistic sử dụng phân bố bernoulli , Nhi chỉ việc dùng lại code của họ.

# Stan code

stan_code="
// data block begins here
data {
int <lower=0> n;                // number  of  instances
int <lower=0> d;               // number  of  predictors
vector[n] y;                   //numeric outcome
matrix[n,d] x;                 //  Design matrix 
real <lower=0>  scale_icept;   // prior scale parameter for  the  intercept
real <lower=0>  scale_global; // scale  for  the half-t prior  for  tau
real <lower=1>  nu_global;     // degrees  of  freedom  for  the half -t priors  for  tau
real <lower=1> nu_local;      // degrees  of  freedom  for  the half -t priors  for  lambdas
real <lower=0>  slab_scale;    // slab  scale  for  the  regularized  horseshoe
real <lower=0> slab_df;       // slab  degrees  of  freedom  for the  regularized  horseshoe
}   // data block ends here

// Parameter blocks begins here

parameters {
real  logsigma;
real  beta0;
vector[d] z;
real <lower=0>  aux1_global;
real <lower=0>  aux2_global;
vector <lower =0>[d] aux1_local;
vector <lower =0>[d] aux2_local;
real <lower=0> caux;
}  

transformed  parameters {
real <lower=0> sigma;                 // noise  scale parameter
real <lower=0> tau;                   // global  shrinkage  parameter
vector <lower =0>[d] lambda;          // local  shrinkage  parameter
vector <lower =0>[d] lambda_tilde;   // ’truncated ’ local  shrinkage  parameter
real <lower=0> c;                      // slab  scale
vector[d] beta;                       // regression  coefficients
vector[n] f;                           // latent  function  values
sigma = exp(logsigma );
lambda = aux1_local  .* sqrt(aux2_local );
tau = aux1_global * sqrt(aux2_global) * scale_global*sigma;
c = slab_scale * sqrt(caux);
lambda_tilde = sqrt( c^2 * square(lambda) ./ (c^2 + tau^2* square(lambda )) );
beta = z .*  lambda_tilde*tau;
f = beta0 + x*beta;
}
// Parameter blocks ends here

model {
// half -t priors  for  lambdas  and tau , and  inverse -gamma  for c^2
z~normal(0, 1);

aux1_local∼normal(0, 1);

aux2_local∼inv_gamma (0.5* nu_local , 0.5* nu_local );

aux1_global∼normal(0, 1);

aux2_global∼inv_gamma (0.5* nu_global , 0.5* nu_global );

caux~inv_gamma (0.5* slab_df , 0.5* slab_df );

beta0~normal(0,scale_icept );

y~normal(f, sigma);

}"

Sau đó tạo dữ liệu cho mô hình STAN này, list dữ liệu đầu vào gồm có:

1. Design matrix Xmat (bao gồm Intercept) và vector Response

2. Khai báo các giá trị tham số: p0= số biến dự kiến có hiệu ứng mạnh, thí dụ 4; d= số biến trong design matrix = 14; n= cỡ mẫu; giá trị scale_global được ước tính bằng công thức ; các tham số cho prior khác như nu_global=1, nu_local=1, scale_icept=100 (cho beta intercept), slab_scale=2 (hằng số c), độ tự do cho slab =

Response=df$brozek_fat

Xmat <- model.matrix(~.,df[,-1])

p0 = 4
d = ncol(Xmat)
n = nrow(df)

scale_global = p0/(ncol(Xmat)-p0-1)/sqrt(n)

data_list <- list(y = Response, 
                  d = d,
                  n = n,
                  x=Xmat,
                  scale_icept = 100,
                  scale_global = scale_global,
                  nu_global = 1,
                  nu_local = 1, 
                  slab_scale = 2,
                  slab_df=4)

Sau đó ta compile code mô hình vào sampler qua hàm stan

library(rstan)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

set.seed(123)

fit_stan= stan(data = data_list, 
          model_code = stan_code , 
          chains = 1,
          iter = 2500, 
          warmup = 500, 
          thin = 2,
          control=list(adapt_delta=0.95,max_treedepth=20)
)

## In file included from C:/Users/Admin/Documents/R/win-library/3.5/BH/include/boost/config.hpp:39:0,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/BH/include/boost/math/tools/config.hpp:13,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/var.hpp:7,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:12,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file2e20582d312a.cpp:8:
## C:/Users/Admin/Documents/R/win-library/3.5/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
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
## 
## SAMPLING FOR MODEL '8399c2e5d8b36da20e9c1268d0d2033c' 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 / 2500 [  0%]  (Warmup)
## Iteration:  250 / 2500 [ 10%]  (Warmup)
## Iteration:  500 / 2500 [ 20%]  (Warmup)
## Iteration:  501 / 2500 [ 20%]  (Sampling)
## Iteration:  750 / 2500 [ 30%]  (Sampling)
## Iteration: 1000 / 2500 [ 40%]  (Sampling)
## Iteration: 1250 / 2500 [ 50%]  (Sampling)
## Iteration: 1500 / 2500 [ 60%]  (Sampling)
## Iteration: 1750 / 2500 [ 70%]  (Sampling)
## Iteration: 2000 / 2500 [ 80%]  (Sampling)
## Iteration: 2250 / 2500 [ 90%]  (Sampling)
## Iteration: 2500 / 2500 [100%]  (Sampling)
## 
##  Elapsed Time: 59.462 seconds (Warm-up)
##                162.751 seconds (Sampling)
##                222.213 seconds (Total)

library(broom)

tidyMCMC(fit_stan, conf.int = TRUE, conf.method = "HPDinterval", pars = "beta")->beta_stan

beta_stan$term=colnames(Xmat)

beta_stan

##           term      estimate  std.error    conf.low   conf.high
## 1  (Intercept)  0.0162140658 0.50068044 -1.00990049  1.01002596
## 2          age  0.0210251496 0.02703985 -0.02665005  0.07517061
## 3       weight -0.1047821207 0.04352636 -0.18899519 -0.02011387
## 4       height -0.0356375656 0.09617571 -0.27757721  0.12485590
## 5         neck -0.1938621839 0.21749998 -0.63916929  0.14148045
## 6        chest  0.0009930989 0.05783711 -0.11754580  0.13076210
## 7      abdomen  0.8785487422 0.07256873  0.74642848  1.01479645
## 8          hip -0.0739597853 0.10788929 -0.31094738  0.10679842
## 9        thigh  0.0929407358 0.10748415 -0.06355845  0.32501330
## 10        knee  0.0113087731 0.11149836 -0.22134925  0.27797802
## 11       ankle  0.0285924087 0.11578093 -0.18632014  0.32528990
## 12      biceps  0.0812814654 0.11278566 -0.10052852  0.33205229
## 13     forearm  0.2006673886 0.18594765 -0.06872802  0.57754006
## 14       wrist -0.8128474580 0.60129260 -1.82980265  0.12618305

8 Package rstanarm

Ngoài cách viết code thủ công, package rstanarm là một giao thức tiện dụng hơn để dùng sampler STAN. rstanarm hỗ trợ prior Piironen thông qua hàm hs()

# rstanarm

library(rstanarm)

n = nrow(df)
nX = 13
p0 = 4
global_scale = p0/(nX - p0)/sqrt(n)

fit_stanarm = stan_glm(brozek_fat~ ., 
                                  data = df, 
                                  iter=2500,
                                  warmup=500, 
                                  chains = 1, 
                                  thin = 2, 
                                  refresh = 0,
                                  prior_intercept = normal(0,100), 
                                  prior = hs(df = 4, 
                                             global_df = 1, 
                                             global_scale = global_scale),
                                  prior_aux = cauchy(0,2),
                                  control=list(adapt_delta=0.95,max_treedepth=20)
                                  )

## 
## Gradient evaluation took 0.01 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 100 seconds.
## Adjust your expectations accordingly!
## 
## 
## 
##  Elapsed Time: 7.25 seconds (Warm-up)
##                25.164 seconds (Sampling)
##                32.414 seconds (Total)

tidyMCMC(fit_stanarm$stanfit,conf.int = TRUE, 
         conf.method = "HPDinterval")->beta_stanarm

beta_stanarm

##             term      estimate   std.error      conf.low     conf.high
## 1    (Intercept) -2.837949e+01 15.47142563  -55.19909207    4.40828527
## 2            age  2.383024e-02  0.02673085   -0.02040000    0.08058551
## 3         weight -1.081275e-01  0.04660253   -0.19166770   -0.01360330
## 4         height -4.134557e-02  0.10966803   -0.30118243    0.15836417
## 5           neck -2.614049e-01  0.20917925   -0.65407528    0.07724954
## 6          chest  1.726801e-03  0.06538599   -0.12216647    0.15482328
## 7        abdomen  8.804351e-01  0.07658661    0.73167595    1.02762084
## 8            hip -9.061495e-02  0.11646706   -0.33143812    0.10055583
## 9          thigh  1.290939e-01  0.12043056   -0.05098560    0.38748269
## 10          knee  8.938019e-03  0.12439847   -0.25556659    0.27856085
## 11         ankle  2.388494e-02  0.13009575   -0.24625557    0.28512641
## 12        biceps  9.202022e-02  0.12984876   -0.09846373    0.39373817
## 13       forearm  2.183674e-01  0.17955693   -0.07128410    0.56673512
## 14         wrist -6.578419e-01  0.52219401   -1.66303209    0.10603066
## 15         sigma  4.053744e+00  0.19268328    3.67157192    4.42385861
## 16      mean_PPD  1.888053e+01  0.36285075   18.15397978   19.57771935
## 17 log-posterior -7.805212e+02  5.86496781 -791.74405554 -769.61170182

9 Package brms

Cuối cùng, brms là 1 giao thức khác của sampler STAN, brms mạnh hơn rstanarm vì hỗ trợ nhiều dạng mô hình hơn, tuy nhiên cú pháp khai báo prior Piironen trong brms phức tạp hơn rstanarm một chút. Cho mô hình Gaussian, brms sử dụng 3 class prior là Intercept, beta (b) và sigma

# brms

nX = 13
p0 = 4
par_ratio = p0/(nX - p0)

library(brms)

fit_brms= brm(brozek_fat~ ., 
              data = df,
              iter = 2500, 
              warmup = 500,
              chains = 1, 
              thin = 2,
              refresh = 0,
              prior = c(prior(normal(0, 100),class = "Intercept"), 
                        prior(horseshoe(df = 4,
                                        par_ratio = par_ratio),
                              class = "b"), 
                        prior(cauchy(0, 5),class = "sigma")
                        ),
              control=list(adapt_delta=0.95,max_treedepth=20)
              )

## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## 
##  Elapsed Time: 4.166 seconds (Warm-up)
##                16.871 seconds (Sampling)
##                21.037 seconds (Total)

tidyMCMC(fit_brms$fit,conf.int = TRUE, 
         conf.method = "HPDinterval")->beta_brms

beta_brms

##           term      estimate   std.error     conf.low    conf.high
## 1  b_Intercept -2.831626e+01 15.93898643 -54.19379742  7.492338784
## 2        b_age  2.068595e-02  0.02619802  -0.02400919  0.073916499
## 3     b_weight -1.084322e-01  0.04729948  -0.19390625 -0.009047489
## 4     b_height -4.742829e-02  0.10774528  -0.28352938  0.156593600
## 5       b_neck -2.462822e-01  0.21273066  -0.67589728  0.054713620
## 6      b_chest  1.822540e-03  0.06475696  -0.12749344  0.143689846
## 7    b_abdomen  8.791064e-01  0.07475039   0.74674926  1.016718300
## 8        b_hip -8.854441e-02  0.10865496  -0.33725519  0.067052421
## 9      b_thigh  1.217147e-01  0.11135051  -0.06446909  0.338543806
## 10      b_knee  3.788469e-04  0.12919469  -0.26014343  0.291824528
## 11     b_ankle  2.604558e-02  0.12551830  -0.22191156  0.309632577
## 12    b_biceps  9.238613e-02  0.11849940  -0.09810031  0.346194821
## 13   b_forearm  2.010712e-01  0.17969870  -0.08467930  0.542849804
## 14     b_wrist -5.968681e-01  0.53614990  -1.67633265  0.103817062
## 15       sigma  4.046756e+00  0.18951528   3.70842548  4.454642806
## 16       hs_c2  1.958997e+00  2.82946496   0.18636004  5.326918895

beta_brms$term=c(colnames(Xmat),"sigma","mean_PPD")

10 So sánh giá trị tham số beta của các mô hình

Ta có thể trích xuất giá trị của 14 tham số hồi quy beta trong 8 mô hình đã dựng và so sánh chúng với nhau một cách trực quan bằng heatmap

# Coefficient

glm_dat=data_frame(Term=fit_glm$coefficients%>%names(),
                   Betas=fit_glm$coefficients%>%as.numeric(),
                   Model="GLM")

lasso_dat=data_frame(Term=coef(lasso)%>%as.matrix()%>%rownames(),
                   Betas=coef(lasso)%>%as.numeric(),
                   Model="LASSO")

enet_dat=data_frame(Term=coef(enet)%>%as.matrix()%>%rownames(),
                     Betas=coef(enet)%>%as.numeric(),
                     Model="ENET")

ridge_dat=data_frame(Term=coef(ridge)%>%as.matrix()%>%rownames(),
                    Betas=coef(ridge)%>%as.numeric(),
                    Model="RIDGE")

coef_bms=coef(bmsmod, std.coefs = T, order.by.pip = T, include.constant = T)%>%as.matrix()

bma_dat=data_frame(Term=coef_bms%>%rownames(),
                     Betas=coef_bms[,3]%>%as.numeric(),
                     Model="BMA")

brms_dat=data_frame(Term=beta_brms[1:14,]%>%.$term,
                   Betas=beta_brms[1:14,]%>%.$estimate,
                   Model="BRMS")

stanarm_dat=data_frame(Term=beta_stanarm[1:14,]%>%.$term,
                    Betas=beta_stanarm[1:14,]%>%.$estimate,
                    Model="STANARM")

stan_dat=data_frame(Term=beta_stan[1:14,]%>%.$term,
                       Betas=beta_stan[1:14,]%>%.$estimate,
                       Model="STAN")

comp_dat=rbind(glm_dat,lasso_dat,bma_dat,
               ridge_dat,enet_dat,
               brms_dat,stanarm_dat,stan_dat)

comp_dat%>%filter(Term!="(Intercept)")%>%
  ggplot(aes(x=reorder(Model,Betas),
             y=reorder(Term,Betas),fill=Betas))+
  geom_tile()+
  theme_bw()+
  scale_fill_gradient2(low="#199ac1",high="#e30c37",mid=NA,midpoint = 0)+
  scale_x_discrete()

Kết quả cho thấy cả 3 mô hình STAn thủ công, brms và rstanarm là tương đương, vì thực chất chúng đều dựa vào lõi là STAN sampler, và cùng code cho hàm likelihood, cùng prior. Có thể dùng 1 trong 3 mô hình để so sánh với những mô hình còn lại.

Prior móng ngựa của Piironen gây ra một sự điều chỉnh rất nhẹ cho tham số hồi quy, tương đương với Ridge regularisation và không khác biệt rõ so với mô hình GLM không có hiệu chỉnh. Nó không triệt tiêu biến số như Lasso hay enet. Có vẻ như các tham số của prior đã không được tối ưu ? Hoặc hiệu ứng này chỉ rõ ràng hơn khi áp dụng cho dataset với rất nhiều biến, thí dụ như dataset Leukemia chẳng hạn ?

Những biến có hiệu ứng mạnh thì kết quả đồng nhất ở các mô hình, thí dụ abdomen, wrist, …

11 So sánh trực quan kết quả của 5 mô hình với giá trị thực tế

pred_brms=data_frame(Pred=predict(fit_brms,newdata=df)%>%as.matrix()%>%.[,1],
                     Model="Piironen")

pred_enet=data_frame(Pred=predict(enet,newx=as.matrix(df[,-1]))%>%as.numeric(),
                     Model="ENET")

pred_lasso=data_frame(Pred=predict(lasso,newx=as.matrix(df[,-1]))%>%as.numeric(),
                      Model="LASSO")

pred_ridge=data_frame(Pred=predict(ridge,newx=as.matrix(df[,-1]))%>%as.numeric(),
                      Model="RIDGE")

truth=data_frame(Pred=df$brozek_fat,Model="Truth")

pred_glm=data_frame(Pred=fit_glm$fitted.values,Model="GLM")

comp_pred=rbind(pred_brms,pred_enet,pred_glm,pred_lasso,pred_ridge)

comp_pred%>%mutate(Truth=rep(df$brozek_fat,5))%>%
  ggplot()+
  geom_violin(alpha=0.7,aes(x=Model,y=Truth,group=Model),fill="white")+
  geom_violin(alpha=0.5,aes(x=Model,y=Pred,fill=Model,group=Model),linetype=2)+
  geom_jitter(alpha=0.2,aes(x=Model,y=Pred,group=Model),shape=21,col="black")+
  theme_bw()+scale_y_continuous("Value")+coord_flip()

Trên chính tập dữ liệu gốc, mô hình Bayes với prior Piironen cho ra kết quả tiên lượng có phân phối đồng dạng với giá trị thực tế,độ chính xác của nó tương đương với mô hình tuyến tính cơ bản bằng hàm glm.

12 Kết luận

Trong bài này Nhi đã giới thiệu với các bạn thêm 1 phương pháp chọn lọc biến số/mô hình nữa, và lần này là cho mô hình Bayes viết bằng STAN code. Một công cụ Regularisation cho mô hình Bayes có ý nghĩa quan trọng, vì trong hoàn cảnh dữ liệu hạn chế, nó sẽ giúp tránh nguy cơ diễn giải sai phân bố hậu nghiệm.

Các bạn có thể tải bài báo gốc của Piironen và Vehtari : https://arxiv.org/pdf/1707.01694.pdf, trong phần phụ lục, tác giả có hướng dẫn 2 cách code khác nhau cho mô hình Gaussian, và chỉ thay đổi 3 dòng code thì bạn đã có thể làm mô hình logistic.

Nếu chưa quen code thủ công, bạn cũng có thể chọn cách dùng giao thức qua stanarm và brms, dễ dàng và nhanh chóng hơn.