Hồi quy logistic Bayes

Sử dụng giao thức STAN

1 Giới thiệu: BAYES

Hiện nay, trường phái Bayes ngày càng được phổ biến như một giải pháp thay thế cho những phương pháp thống kê cổ điển dựa theo trường phái frequentist. Thí dụ mô hình hồi quy Bayes thay thế cho phương pháp ước lượng Maximum likelihood. Trường phái Bayes có một vài ưu điểm như:

1. Cho phép hòa hợp giữa yếu tố tiền định (thí dụ giả thuyết, giả định, thông tin trong y văn) và bằng chứng từ dữ liệu quan sát được trong nghiên cứu, thí nghiệm

2. Diễn giải kết quả linh hoạt hơn rất nhiều so với Null hypothesis test và p value, không áp đặt nhưng đưa ra chứng cứ, khả năng tin cậy, và cho phép cập nhật kết quả khi có thêm chứng cứ, dữ liệu mới.

3. Giản dị và đồng nhất : Những phân tích được diễn tả bằng mô hình, bằng xác suất điều kiện, chứ không còn bị phân tán thành từng mảnh nhỏ là các loại kiểm định .

4. Vì phương pháp Bayes yêu cầu người học phải giải phẫu mô hình hết sức chi tiết, và hiểu lý thuyết xác suất, để có thể viết code, nên các bạn thực hành Bayes thường có kiến thức vững vàng về thống kê học hơn so với các bạn chỉ dùng phần mềm thống kê theo trường phái Frequentist.

Tiếc rằng, việc ứng dụng rộng rãi phương pháp Bayes vẫn còn gặp nhiều khó khăn, nhất là tại Việt Nam. Chủ yếu là khó khăn về kỹ thuật, vì người thực hành phải có khả năng lập trình, viết code cho các Sampler để đạt được phân phối hậu nghiệm. Các ngôn ngữ và Sampler như BUGs, JAGS hay STAN do đó không mấy phổ biến.

Tuy nhiên, rất may mắn là có một số giao thức đã được phát triển trên nền tảng của JAGs và STAN trong R. Một trong những giao thức đó là rstanarm, Nhi sẽ thực hiện một số tutorial hướng dẫn cách dùng giao thức này để thực hiện các bài toán hồi quy phổ biến. Bài đầu tiên sẽ là Hồi quy Logistic.

2 STAN và rstanarm package

STAN là một ngôn ngữ lập trình thống kê xác suất,và đã được phát triển cho R, Python và Matlab. STAN cho phép người dùng dựng những mô hình thông qua xác suất có điều kiện (xác định phân phối hậu định của tham số hồi quy theta khi có dữ liệu, hàm likelihood và giả thuyết tiền định), hay còn gọi là trường phái Bayes.

Khi bạn dựng một mô hình trong STAN, nó phải được giải phẫu đến từng chi tiết nhỏ nhất (trường phái Bayes bắt buộc người học phải hiểu chi tiết algorithm, bản chất biến số, các hàm mật độ xác suất). Khi bạn đã có mô hình, thì có thể giải quyết hầu hết những bài toán thống kê, như so sánh, tương quan, hồi quy, tiên lượng… Từ chuỗi MCMC (phân phối hậu định), bạn có thể suy diễn thống kê theo nhiều cách: dùng Bayes factor, ROPE, null hypothesis testing. Mô hình Bayes có thể áp dụng để diễn giải hay tiên lượng đều được, thậm chí nó cho phép ướng lượng phân phối hậu định cho giá trị dự báo, hồi quy phân vị… một cách tự nhiên.

Khi bạn dùng STAN, đầu tiên 1 đoạn code của mô hình sẽ được phiên mã thành một chương trình (executable program) ngôn ngữ C++ (compilation), sau đó nó mới được thi hành. Cơ chế tạo ra chuỗi Markov Monte Carlo của STAN hoàn toàn khác biệt với JAGs và BUGs là 2 sampler có trước. STAN sử dụng phương pháp Hamiltonian Monte Carlo sampling dựa vào lý thuyết của Metropolis-Hastings (1953, 1970), mô phỏng chuyển động bất định của một động tử.

Bạn có thể tiếp cận STAN tại đây:

https://github.com/stan-dev

Trong R, STAN được ứng dụng dựa vào 1 số package chính, trong đó phần lõi là Rcpp, rstan. Dựa trên phần lõi này, ta có những package giao thức như bayesplot, rstanarm, brms… cho phép người dùng có thể dựng được mô hình, làm được suy diễn Bayes mà không cần biết ngôn ngữ STAN. Các giao thức này sử dụng cú pháp quy ước trong R (glm, lm, aov, gam, nlme, survival…) để khai báo mô hình, sau đó dịch sang STAN code, compile vào C++ rồi thi hành. Kết quả xuất ra có thể tương thích với những method quen thuộc, thí dụ ggplot. Đó là lợi thế của các giao thức.

3 Hồi quy logistic theo BAYES

Một cách tổng quát, suy diễn Bayes dựa trên định lý Bayes như sau :

\[ p(\theta |y,x)\propto p(y|\theta ,x) p(\theta ,x) \] theo đó, phân phối hậu định của một tham số Theta (khi ta có trong tay dữ liệu biến kết quả Y và hằng số X là một matrix các predictors trong mô hình) tỉ lệ với tích của hàm likelihood (xác suất có điều kiện của y khi có theta và x) và phân phối tiền định (prior= giả thuyết về phân phối của theta).

Bài toán dán nhãn phân loại một bệnh nhân với hai khả năng (0 và 1) tương đương với việc ước tính xác suất có y=1 trong n lượt rút thăm ngẫu nhiên từ n phần tử. Xác suất này được mô tả bằng 1 hàm mật độ xác suất (PMF) của phân phối binomial:

\[\binom{n}{y} \pi^{y} (1 - \pi)^{n - y}\]

Với pi là xác suất tìm thấy y=1.

Cho 1 link function g ta có :

\[\pi = g^{-1}(\eta)\]

Cho mô hình logistic, g là 1 hàm logit, hay còn gọi là log-odds

\[g(x) = \ln{\left(\frac{x}{1-x}\right)}\]

Từ đây ta có mô hình hồi quy logistic :

\[\binom{n}{y}\left(\text{logit}^{-1}(\eta)\right)^y \left(1 - \text{logit}^{-1}(\eta)\right)^{n-y} = \binom{n}{y} \left(\frac{e^{\eta}}{1 + e^{\eta}}\right)^{y} \left(\frac{1}{1 + e^{\eta}}\right)^{n - y}\]

Với eta là kết quả ước lượng từ một phương trình tuyến tính có dạng:

\[\eta = \alpha + \mathbf{x}^\top \boldsymbol{\beta}\]

Trong đó alpha là intercept, X là design matrix của mô hình, beta là tham số hồi quy.

Trong trường hợp hồi quy logistic với một dataset có N quan sát, K biến số, và kết quả Y là một biến nhị phân (0,1); phân phối hậu định trong Bayes có thể được viết như sau:

\[f\left(\alpha,\boldsymbol{\beta} | \mathbf{y},\mathbf{X}\right) \propto f\left(\alpha\right) \times \prod_{k=1}^K f\left(\beta_k\right) \times \prod_{i=1}^N { g^{-1}\left(\eta_i\right)^{y_i} \left(1 - g^{-1}\left(\eta_i\right)\right)^{n_i-y_i}}\]

Với alpha và beta(k) lần lượt là intercept và các tham số hồi quy cho K biến số trong mô hình.

Để xác định phân phối hậu định cho alpha và beta, ta phảcung cấp cho Sampler một khái niệm (giả định) về phân phối của chúng (priors) là f(alpha) và f(beta). Prior này có thể là normal hoặc Student-t.

4 Dữ liệu minh họa

Để minh họa, 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)

Dataset này đã được thăm dò rất chi tiết trong 1 tutorial trước đây (Gamlss bài 3):

http://rpubs.com/lengockhanhi/294782

Chúng ta chia dữ liệu ra 2 phần: trainset và testset như thường lệ:

library(tidyverse)

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)

library(caret)

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

5 Mô hình logistic dùng GLM

Đầu tiên, Nhi dựng một mô hình logistic theo phương pháp quen thuộc trong R là hàm glm:

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

6 Dựng một mô hình logistic thủ công dùng STAN

Quy trình sau đây tương đối phức tạp, nên nếu các bạn không có hứng thú có thể bỏ qua. Nhi muốn minh họa cách code 1 mô hình trong STAN. Đầu tiên chúng ta phải khai báo về data, bao gồm tên biến số. Sau đó, ta phải khai báo về biến số thật chi tiết, kiểu biến số (integer), các khoảng chặn, danh sách tham số hồi quy, chuyển dạng thành xác suất, cuối cùng là nội dung mô hình. Sau đó ta cho STAN thi hành việc lấy mẫu cho chuỗi MCMC với 2000 lượt bao gồm 500 lượt warm-u.

Kết quả cuối cùng là phân phối hậu định của 6 tham số hồi quy beta:

library(rstan)

## Create Stan data
data=list(N= nrow(trainset),
            p        = 6,
            Class    = trainset$Class,
            clmptkns=trainset$clumpthickness,
            sunif=trainset$ShapeUniformity,
            mad=trainset$Margin_adhesion,
            bn=trainset$Barenuclei,
            bchro=trainset$BlandChromatin
          )

stancode="data {
 // Khai bao bien so 

// Khai bao co mau (so nguyen)
int<lower=0> N;

// Khai bao so luong bien so

int<lower=0> p;

// Khai bao ten bien so va ban chat (so nguyen)

int <lower=0,upper=1> Class[N];
int<lower=0>  clmptkns[N];
int<lower=0>  sunif[N];
int<lower=0> mad[N];
int<lower=0>  bn[N];
int<lower=0> bchro[N];
}

parameters {
// Khai bao tham so hoi quy beta
real beta[p];
}

transformed parameters  {
   // Hoan chuyen sang xac suat
real<lower=0> odds[N];
real<lower=0, upper=1> prob[N];

for (i in 1:N) {
odds[i] = exp(beta[1] + beta[2]*clmptkns[i] + beta[3]*sunif[i] + beta[4]*mad[i] + beta[5]*bn[i] + beta[6]*bchro[i]);
prob[i] = odds[i] / (odds[i] + 1);
}
}

model {
// Prior cua mo hinh (flat)

// Thanh phan Likelihood cua suy dien Bayesian
Class ~ bernoulli(prob);
}"

mod=stan(model_code = stancode, 
         data = data,
        chains = 2, iter = 2000, warmup = 500,
        cores=parallel::detectCores())

## In file included from C:/Users/Admin/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file2ac052fc7dd9.cpp:8:
## C:/Users/Admin/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
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"

print(mod,pars = c("beta"))

## Inference for Stan model: 09413ad80138a7d73f1d3cf363aaa40d.
## 2 chains, each with iter=2000; warmup=500; thin=1; 
## post-warmup draws per chain=1500, total post-warmup draws=3000.
## 
##          mean se_mean   sd   2.5%    25%   50%   75% 97.5% n_eff Rhat
## beta[1] -9.70    0.03 1.18 -12.14 -10.49 -9.65 -8.87 -7.56  1529    1
## beta[2]  0.63    0.00 0.15   0.35   0.52  0.62  0.73  0.94  1646    1
## beta[3]  0.52    0.00 0.16   0.22   0.41  0.52  0.63  0.86  2442    1
## beta[4]  0.34    0.00 0.12   0.11   0.26  0.34  0.42  0.59  3000    1
## beta[5]  0.37    0.00 0.10   0.17   0.30  0.37  0.44  0.57  2447    1
## beta[6]  0.47    0.00 0.18   0.14   0.35  0.47  0.59  0.83  2311    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Sep 06 09:22:13 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

7 Mô hình logistic sử dụng rstanarm package

Bây giờ, chúng ta sẽ đi theo con đường đơn giản hơn rất nhiều đó là dùng giao thức rstanarm.

Cú pháp của hàm stan_glm hoàn toàn giống như hàm glm(), chỉ khác là ta có thêm 1 số arguments về: giả thuyết tiền định prior cho intercept và các tham số hồi quy, quy trình lấy mẫu cho chuỗi MCMC (ở đây Nhi dùng 2 chuỗi, mỗi chuỗi 1500 lượt, bao gồm 500 lượt warming-up, chạy song song trên 4 core). Bạn nhớ xác định giá trị seed để tái lập kết quả.

Ta có thể thí dụ: giả thuyết tiền định (piror) là Student-t với df=5, trung bình = 0 và scale=2.5;

Mô hình sẽ được rstan compile thành C++ program và chạy với tốc độ converge rất nhanh chỉ trong vài phút:

library(rstanarm)

t_prior=student_t(df = 5, location = 0, scale = 2.5)

fml="Class~clumpthickness+
  ShapeUniformity+
Margin_adhesion+
Barenuclei+
BlandChromatin"

fittprior=stan_glm(
  fml,
  data = trainset, 
  family = binomial(link = "logit"), 
  prior = t_prior,
  prior_intercept = t_prior,
  algorithm="sampling",
  chains = 2, iter = 1500,warmup = 500,
  cores=parallel::detectCores(),
  seed=123
)

8 Khai thác mô hình logistic Bayes

Đầu tiên ta dùng hàm summary: Nội dung mô hình rất gần với kết quả của phương pháp REML

Sử dụng hàm exp(), ta tính được OR

summary(fittprior,digits=5,probs=c(0.025,0.975))

## 
## Model Info:
## 
##  function:  stan_glm
##  family:    binomial [logit]
##  formula:   "Class~clumpthickness+\n ShapeUniformity+\nMargin_adhesion+\nBarenuclei+\nBlandChromatin"
##  algorithm: sampling
##  priors:    see help('prior_summary')
##  sample:    2000 (posterior sample size)
##  num obs:   500
## 
## Estimates:
##                   mean      sd        2.5%      97.5%  
## (Intercept)      -9.35274   1.07046 -11.58034  -7.39699
## clumpthickness    0.59474   0.13949   0.33598   0.88405
## ShapeUniformity   0.51011   0.16115   0.21630   0.85396
## Margin_adhesion   0.32270   0.11846   0.09652   0.56287
## Barenuclei        0.36470   0.09297   0.18309   0.55348
## BlandChromatin    0.46691   0.16797   0.15067   0.81076
## mean_PPD          0.35065   0.00999   0.33000   0.37000
## log-posterior   -58.89030   1.63128 -62.69783 -56.54771
## 
## Diagnostics:
##                 mcse    Rhat    n_eff
## (Intercept)     0.02697 0.99942 1575 
## clumpthickness  0.00364 1.00030 1472 
## ShapeUniformity 0.00392 0.99948 1693 
## Margin_adhesion 0.00293 0.99995 1633 
## Barenuclei      0.00227 1.00025 1674 
## BlandChromatin  0.00406 0.99968 1713 
## mean_PPD        0.00022 1.00023 2000 
## log-posterior   0.05397 1.00064  914 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

posterior_interval(fittprior)%>%exp()

##                           5%          95%
## (Intercept)     1.396778e-05 0.0004701105
## clumpthickness  1.448342e+00 2.3156126545
## ShapeUniformity 1.295290e+00 2.1984733924
## Margin_adhesion 1.141195e+00 1.6927596223
## Barenuclei      1.230699e+00 1.6758497047
## BlandChromatin  1.215031e+00 2.1168243484

Ta vẽ được 1 biểu đồ marginalized effect cho OR:

ordf=exp(cbind(OR = coef(fittprior), posterior_interval(fittprior)))%>%
  as_tibble()%>%
  mutate(.,Predictor=fittprior$covmat%>%rownames)

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,show.legend = F)+
  geom_point(aes(y=OR),shape=21,size=5,fill="white",stroke=1.5,show.legend = F)+
  theme_bw()+
  coord_flip()+
  geom_hline(yintercept = 1,size=1,color="blue",linetype=2)+
  scale_x_discrete("Predictors")+
  scale_y_continuous("Odds-Ratio")

Hàm as.data.frame cho phép trích xuất trực tiếp 6 chuỗi MCMC từ mô hình Bayes:

Bây giờ ta có thể vẽ được phân phối hậu định của odds-ratio cho 6 tham số hồi quy

posterior=as.data.frame(fittprior)

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

(Intercept)	clumpthickness	ShapeUniformity	Margin_adhesion	Barenuclei	BlandChromatin
-7.975938	0.4612697	0.5978315	0.0613562	0.4350187	0.3733282
-7.582937	0.3877390	0.5682261	0.0235291	0.4186240	0.4031363
-9.319668	0.5947666	0.5949226	0.2925004	0.4407699	0.3757688
-9.526398	0.6752297	0.2945156	0.3742914	0.3731516	0.4677688
-8.085083	0.4918975	0.5138523	0.2328746	0.3753085	0.3650310
-11.037517	0.7963245	0.5443347	0.4349765	0.4546275	0.3139365

posterior%>%mutate(.,Iteration=as.numeric(rownames(.)))%>%
  gather(`(Intercept)`:BlandChromatin,
                        key="Predictors",
                        value="Coefficients")%>%
  ggplot(aes(x=Iteration,
             y=exp(Coefficients),
             color=Predictors))+
  geom_line(alpha=0.7,show.legend = F)+
  geom_hline(yintercept=1,size=1,linetype=2,color="red4")+
  scale_y_continuous("Odds-ratio")+
  facet_wrap(~Predictors,scales = "free",ncol=2)+
  theme_bw()

library(ggjoy)

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

posterior[,-1]%>%gather(clumpthickness:BlandChromatin,
                        key="Predictors",
                        value="Coefficients")%>%
  ggplot(aes(x=exp(Coefficients),
             y=reorder(Predictors,Coefficients),
             fill=reorder(Predictors,Coefficients)))+
  geom_joy(alpha=0.6,show.legend = F,scale=3)+
  scale_y_discrete("Predictors")+
  scale_x_continuous("Odds-ratio")+
  geom_vline(xintercept=c(1,2),size=1,linetype=2,color=c("red4","blue"))+
  theme_bw()

## Picking joint bandwidth of 0.0425

Ta có thể kiểm định kết quả của mô hình BAYES trên testset với package caret:Mô hình này hoạt động rất tốt, không thua kém mô hình dựng bằng GLM hay Gamlss

testpred=predict(fittprior,newdata=testset,type="link")%>%exp()

testset$Truth=ifelse(testset$Class==1,"Malignant","Benign")
testset$Pred=testpred
testset$Pred=ifelse(testset$Pred>=0.5,"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       
## 

Bonus: Cơ chế hoạt động của mô hình, từ góc nhìn 2 chiều giữa 2 biến số bất kì trong mô hình:

pr=function(x, y, ests) plogis(ests[1] + 
                                 ests[2] * x +
                                 ests[3]*3.164+
                                 ests[4]*3.23+
                                 ests[6] *y +
                                 ests[5]*3.564)

grid=expand.grid(clumpthickness = seq(0, 10, length.out = 50), 
                 BlandChromatin= seq(0, 10, length.out = 50))

grid$prob=with(grid, pr(clumpthickness, BlandChromatin, coef(fittprior)))

ggplot(grid, aes(x = clumpthickness, 
                 y = BlandChromatin)) + 
  geom_tile(aes(fill = prob),alpha=0.6) + 
  geom_point(data = testset,
             aes(color = factor(Class)),
             size = 3, alpha = 0.8) + 
  scale_fill_gradient(low="#7ee1ea",high="#ed283b")+
  scale_color_manual("Class", 
                     values = c("blue3", "red3"), 
                     labels = c("Neg", "Pos")
                     )+
  scale_x_continuous(breaks=c(0,1,2,3,4,5,6,7,8,9,10))+
  scale_y_continuous(breaks=c(0,1,2,3,4,5,6,7,8,9,10))+
  theme_bw()

9 Kết luận

Với những giao thức hồi quy Bayes trong R, chúng ta hoàn toàn có thể thay thế các phương pháp REML cổ điển bằng trường phái BAYES mà không cần phải biết về ngôn ngữ lập trình cho sampler như STAN.

Bản thân ngôn ngữ STAN là một phát triển rất hữu ích, nó cho phép dựng những mô hình phức tạp mà trước đó người dùng gặp trở ngại với JAGs hoặc BUGs. Có 2 giao thứccho STAN trong R đó là rstanarm và brms. Giao thức brms đã được giới thiệu vào năm 2016 trong dự án Bayes for Vietnam, và vẫn đang được phát triển.

Giao thức rstanarm đơn giản và dễ sử dụng hơn brms, nó phù hợp hơn cho những nghiên cứu thông thường, và hỗ trợ hầu hết các mô hình GLM và MGLM, bao gồm hồi quy với random effect, beta, gamma, và hồi quy cho count data (negative binomial).

Trong thời gian tới Nhi sẽ giới thiệu với các bạn thêm một vài mô hình Bayes khác.

Nếu các bạn có hứng thú tham gia reboot cho dự án Bayes for Vietnam, xin liên lạc với nhóm chúng tôi.

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