Hồi quy Bayes bằng TensorFlow

1 Giới thiệu

2 Giới thiệu

2.1 TensorFlow và greta

Khi nhắc đến TensorFlow, một thư viện lập trình của Google, đa số chúng ta sẽ nghĩ đến ứng dụng phổ biến nhất của nó là Deep learning. Tuy nhiên tính năng « lập trình luồng dữ liệu » (dataflow programming) của TensorFlow còn có thể ứng dụng trong nhiều việc khác thú vị không kém, như Thống kê Bayes.

Trước đây, thống kê Bayes chủ yếu được thực hiện qua 3 hệ ngôn ngữ chính: BUGs, JAGs và STAN. Một số R package được tạo ra như giao thức gồm mcmcpack, rstanarm, rstan, và brms… Tuy nhiên sự phức tạp của các ngôn ngữ này khiến thống kê Bayes ít được phổ biến trong cộng đồng.

Vào đầu năm 2018, Nick Golding giới thiệu greta, một thư viện R code cho phép dựng mô hình thống kê Bayes bằng TensorFlow. Như vậy, có thể xem greta như một thư viện ngôn ngữ lập trình Bayes thế hệ mới (sau BUGs, JAGs và STAN). greta được gọi theo tên của nhà toán học nữ Grete Hermann (1901 – 1984) , cũng như ngôn ngữ STAN được đặt theo tên của một nhà toán học khác là Stanislaw Ulam.

Nhi đã nghe nói về greta từ tháng 4/2018, tuy nhiên chỉ đến khi trang bị xong 1 laptop với GPU cách đây vài ngày thì Nhi mới bắt đầu làm việc với ngôn ngữ này, và hết sức ngạc nhiên vì sự đơn giản và sức mạnh của nó. Công việc được giản lược tối đa, nhưng vẫn giữ được lối suy nghĩ hệ thống, thứ bậc như quy trình cổ điển bằng JAGs hay STAN.

2.2 Chuẩn bị

Nếu máy tính của bạn có GPU NVIDIA tương thích với CUDA :

1. Trước hết bạn install CUDA Toolkit 9.0 :

https://developer.nvidia.com/cuda-90-download-archive?target_os=Windows&target_arch=x86_64&target_version=10&target_type=exelocal

2. Sau đó bạn install Python - ANACONDA từ:

https://www.anaconda.com/download/#windows

3. Tiếp theo, bạn install TensorFlow và Tensorflow Probability

https://www.tensorflow.org/install/install_windows

https://github.com/tensorflow/probability

Ghi chú: chọn bản dành cho GPU, và làm trên Anaconda Commmand Prompt, một số package khác có thể được yêu cầu.

Bạn vào folder của CUDA, “C:FilesGPU Computing Toolkit.0” , copy file cudart64_90.dll và dán vào thư mục DLLS nơi install Anaconda3

4. Cuối cùng, trong R Studio: bạn thực hiện thao tác cuối cùng như hướng dẫn :

https://greta-dev.github.io/greta/get_started.html

Nếu máy tính của bạn không có GPU, không sao cả - vì TensorFlow chạy tốt trên cả CPU và GPU, bạn chỉ cần chọn đúng bản Tensorflow cho CPU ở bước 3, nhưng không cần đến CUDA.

3 Thí dụ minh họa

3.1 Mục tiêu và dữ liệu

Trong bài này, Nhi chọn một nghiên cứu có thiết kế ANOVA đơn biến, với mục tiêu là so sánh bề dày màng phế nang mao mạch phổi được đo bằng kỹ thuật không xâm lấn, giữa 3 phân nhóm: Emphysema (khí phế thũng), Fibrosis (Xơ phổi) và nhóm Chứng (bình thường).

library(tidyverse)

df=read.csv("https://raw.githubusercontent.com/kinokoberuji/R-Tutorials/master/aerodim.csv",sep=";")

head(df,5)

library(ggridges)

myfill3=c("#00c1ed","#8600ed","#ed003f")
mycol3=c("#0162bc","#58019b","#9b011d")

df%>%ggplot(aes(x=Thickness,y=Diagnostic,fill=Diagnostic,col=Diagnostic))+
  geom_density_ridges(alpha=0.7,size=1,scale=1,show.legend = F)+
  geom_rug(alpha=0.7,show.legend = F)+
  scale_x_continuous("Thickness (µm)")+
  theme_bw()+
  coord_flip()+
  scale_fill_manual(values=myfill3)+ scale_color_manual(values=mycol3)

df%>%
  group_by(Diagnostic)%>%
  summarise_at("Thickness",funs( Median=median,Mean=mean,SD=sd,
                         LL=quantile(.,probs=0.025),
                         UL=quantile(.,probs=0.975)))

3.2 Thống kê BAYES bằng 5 công cụ khác nhau

3.2.1 Mô hình

Với thí dụ minh họa này, mục tiêu của chúng ta là dựng một mô hình hồi quy tuyến tính có nội dung:

\[Thickness \sim \beta_{1}(Emphysema) + \beta_{2}(Fibrosis) + \beta_{3}(Normal)\]

\[Thickness \sim N(\mu , \sigma )\]

\[\mu \sim \beta_{1}(Emphysema) + \beta_{2}(Fibrosis) + \beta_{3}(Normal)\]

Priors: những điều kiện phân phối tiền định

\[\beta \sim N(0,10)\]

\[\sigma \sim Cauchy(0,5)\]

3.2.2 JAGs

modelString = "
  model {

#Likelihood
for (i in 1:n) {

y[i]~dnorm(mu[i],tau)
mu[i] <- inprod(beta[],X[i,])
}

#Priors
for (j in 1:nX) {
beta[j] ~ dnorm(0,10)
}

tau <- 1 / (sigma * sigma)
sigma~dunif(0,100)
}
"

X = model.matrix(~Diagnostic-1, data = df)

data.list <- with(df, list(y = Thickness, X = X, nX = ncol(X), n = nrow(df)))


params <- c("beta", "sigma")
nChains = 2
burnInSteps = 1000
thinSteps = 2
numSavedSteps = 6000  
nIter = ceiling(burnInSteps + (numSavedSteps * thinSteps)/nChains)

library(R2jags)

## Loading required package: rjags

## Loading required package: coda

## Linked to JAGS 4.3.0

## Loaded modules: basemod,bugs

## 
## Attaching package: 'R2jags'

## The following object is masked from 'package:coda':
## 
##     traceplot

library(bayesplot)

## This is bayesplot version 1.5.0

## - Plotting theme set to bayesplot::theme_default()

## - Online documentation at mc-stan.org/bayesplot

mod.JAGs <- jags(data = data.list, inits = NULL, parameters.to.save = params,
                        model.file = textConnection(modelString), n.chains = nChains, n.iter = nIter,
                        n.burnin = burnInSteps, n.thin = thinSteps)

## module glm loaded

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 38
##    Unobserved stochastic nodes: 4
##    Total graph size: 206
## 
## Initializing model

print(mod.JAGs)

## Inference for Bugs model at "5", fit using jags,
##  2 chains, each with 7000 iterations (first 1000 discarded), n.thin = 2
##  n.sims = 6000 iterations saved
##          mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat
## beta[1]    0.304   0.032   0.242   0.283   0.305   0.325   0.367 1.001
## beta[2]    0.910   0.039   0.830   0.884   0.910   0.936   0.986 1.001
## beta[3]    0.586   0.030   0.526   0.565   0.586   0.607   0.645 1.001
## sigma      0.120   0.015   0.095   0.109   0.118   0.129   0.153 1.001
## deviance -55.089   3.133 -58.995 -57.392 -55.771 -53.570 -47.248 1.001
##          n.eff
## beta[1]   6000
## beta[2]   6000
## beta[3]   6000
## sigma     6000
## deviance  6000
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 4.9 and DIC = -50.2
## DIC is an estimate of expected predictive error (lower deviance is better).

library(coda)
data.mcmc = as.mcmc(mod.JAGs)

mcmc_areas(data.mcmc,pars=c("beta[1]","beta[2]","beta[3]"))

mcmc_trace(data.mcmc,  
           pars=c("beta[1]","beta[2]","beta[3]"),
           facet_args = list(ncol = 1, strip.position = "left"))

3.2.3 MCMCpack

library(MCMCpack)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)

## ## Copyright (C) 2003-2018 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park

## ##
## ## Support provided by the U.S. National Science Foundation

## ## (Grants SES-0350646 and SES-0350613)
## ##

mod.mcmcpack <- MCMCregress(Thickness~Diagnostic-1, data = df)

summary(mod.mcmcpack)

## 
## Iterations = 1001:11000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                Mean       SD  Naive SE Time-series SE
## DiagnosticE 0.30715 0.031580 3.158e-04      0.0003190
## DiagnosticF 0.92225 0.039487 3.949e-04      0.0003949
## DiagnosticN 0.59066 0.030939 3.094e-04      0.0003143
## sigma2      0.01406 0.003571 3.571e-05      0.0000386
## 
## 2. Quantiles for each variable:
## 
##                 2.5%     25%     50%     75%   97.5%
## DiagnosticE 0.246128 0.28608 0.30706 0.32816 0.36943
## DiagnosticF 0.844825 0.89586 0.92265 0.94838 1.00019
## DiagnosticN 0.529284 0.57057 0.59067 0.61118 0.65192
## sigma2      0.008771 0.01148 0.01347 0.01594 0.02246

library(bayesplot)

posterior <- as.array(mod.mcmcpack)

mcmc_areas(posterior,  
           pars = c("DiagnosticE", "DiagnosticF", "DiagnosticN"))

mcmc_trace(posterior,  
           pars = c("DiagnosticE", "DiagnosticF", "DiagnosticN"),
           facet_args = list(ncol = 1, strip.position = "left"))

3.2.4 STAN

# STAN

model.text <-  "

data { int<lower=1> n;   // total number of observations 
vector[n] Y;      // response variable 
int<lower=1> nX;  // number of effects 
matrix[n, nX] X;   // model matrix 
} 

parameters { 
vector[nX] beta;  // population-level effects 
real<lower=0> sigma;  // residual SD 
} 
transformed parameters { 
} 
model { 
vector[n] mu; 
mu = X * beta ; 
// prior specifications 
beta ~ normal(0, 10); 
sigma ~ cauchy(0, 5); 

// likelihood contribution 
Y ~ normal(mu, sigma); 
} 

generated quantities {
vector[n] log_lik;
for (i in 1:n) {
log_lik[i] = normal_lpdf(Y[i] | X[i] * beta , sigma);
} 
}
"

X = model.matrix(~Diagnostic-1, data = df)
data.list <- with(df, list(Y = Thickness, X = X, nX = ncol(X), n = nrow(df)))

nChains = 2
burnInSteps = 1000
thinSteps = 2
numSavedSteps =  5000 
nIter = ceiling(burnInSteps + (numSavedSteps * thinSteps)/nChains)

library(rstan)

## Loading required package: StanHeaders

## rstan (Version 2.17.3, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

## 
## Attaching package: 'rstan'

## The following object is masked from 'package:R2jags':
## 
##     traceplot

## The following object is masked from 'package:coda':
## 
##     traceplot

## The following object is masked from 'package:tidyr':
## 
##     extract

library(broom)

mod.STAN<- stan(data = data.list, model_code = model.text, chains = nChains,
                       iter = nIter, warmup = burnInSteps, thin = thinSteps)

## In file included from C:/Users/bacsi/Documents/R/win-library/3.5/BH/include/boost/config.hpp:39:0,
##                  from C:/Users/bacsi/Documents/R/win-library/3.5/BH/include/boost/math/tools/config.hpp:13,
##                  from C:/Users/bacsi/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/var.hpp:7,
##                  from C:/Users/bacsi/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
##                  from C:/Users/bacsi/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:12,
##                  from C:/Users/bacsi/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/bacsi/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/bacsi/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file5083ac6584e.cpp:8:
## C:/Users/bacsi/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
## In file included from C:/Users/bacsi/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:44:0,
##                  from C:/Users/bacsi/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/bacsi/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/bacsi/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file5083ac6584e.cpp:8:
## C:/Users/bacsi/Documents/R/win-library/3.5/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 '5c7a181aa62fc47b0f5448f141d8137d' 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 / 6000 [  0%]  (Warmup)
## Iteration:  600 / 6000 [ 10%]  (Warmup)
## Iteration: 1001 / 6000 [ 16%]  (Sampling)
## Iteration: 1600 / 6000 [ 26%]  (Sampling)
## Iteration: 2200 / 6000 [ 36%]  (Sampling)
## Iteration: 2800 / 6000 [ 46%]  (Sampling)
## Iteration: 3400 / 6000 [ 56%]  (Sampling)
## Iteration: 4000 / 6000 [ 66%]  (Sampling)
## Iteration: 4600 / 6000 [ 76%]  (Sampling)
## Iteration: 5200 / 6000 [ 86%]  (Sampling)
## Iteration: 5800 / 6000 [ 96%]  (Sampling)
## Iteration: 6000 / 6000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.031 seconds (Warm-up)
##                0.125 seconds (Sampling)
##                0.156 seconds (Total)
## 
## 
## SAMPLING FOR MODEL '5c7a181aa62fc47b0f5448f141d8137d' 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 / 6000 [  0%]  (Warmup)
## Iteration:  600 / 6000 [ 10%]  (Warmup)
## Iteration: 1001 / 6000 [ 16%]  (Sampling)
## Iteration: 1600 / 6000 [ 26%]  (Sampling)
## Iteration: 2200 / 6000 [ 36%]  (Sampling)
## Iteration: 2800 / 6000 [ 46%]  (Sampling)
## Iteration: 3400 / 6000 [ 56%]  (Sampling)
## Iteration: 4000 / 6000 [ 66%]  (Sampling)
## Iteration: 4600 / 6000 [ 76%]  (Sampling)
## Iteration: 5200 / 6000 [ 86%]  (Sampling)
## Iteration: 5800 / 6000 [ 96%]  (Sampling)
## Iteration: 6000 / 6000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.031 seconds (Warm-up)
##                0.169 seconds (Sampling)
##                0.2 seconds (Total)

mod.STAN%>%summary(pars=c("beta[1]","beta[2]","beta[3]"))

## $summary
##              mean      se_mean         sd      2.5%       25%       50%
## beta[1] 0.3078440 0.0004637065 0.03192986 0.2455366 0.2871185 0.3067391
## beta[2] 0.9230479 0.0005859699 0.04020586 0.8437336 0.8963110 0.9233861
## beta[3] 0.5906047 0.0004501829 0.03032828 0.5297824 0.5706931 0.5903982
##               75%     97.5%    n_eff      Rhat
## beta[1] 0.3282565 0.3737702 4741.411 1.0004421
## beta[2] 0.9505584 1.0012133 4707.910 0.9998969
## beta[3] 0.6104716 0.6502504 4538.554 1.0003613
## 
## $c_summary
## , , chains = chain:1
## 
##          stats
## parameter      mean         sd      2.5%       25%       50%       75%
##   beta[1] 0.3089488 0.03189869 0.2477589 0.2882630 0.3076675 0.3294044
##   beta[2] 0.9228189 0.04043035 0.8446953 0.8956350 0.9226338 0.9506710
##   beta[3] 0.5899114 0.03052021 0.5293108 0.5699021 0.5893263 0.6093811
##          stats
## parameter     97.5%
##   beta[1] 0.3738110
##   beta[2] 1.0010030
##   beta[3] 0.6497666
## 
## , , chains = chain:2
## 
##          stats
## parameter      mean         sd      2.5%       25%       50%       75%
##   beta[1] 0.3067391 0.03192915 0.2449166 0.2852530 0.3063527 0.3269351
##   beta[2] 0.9232770 0.03998689 0.8427820 0.8973373 0.9242943 0.9498387
##   beta[3] 0.5912979 0.03012527 0.5298942 0.5716788 0.5914112 0.6111687
##          stats
## parameter     97.5%
##   beta[1] 0.3731337
##   beta[2] 1.0013581
##   beta[3] 0.6504075

tidyMCMC(mod.STAN,pars=c("beta[1]","beta[2]","beta[3]"),
         conf.int = TRUE, 
         conf.method = "HPDinterval")

data.mcmc2 = as.array(mod.STAN)

mcmc_areas(data.mcmc2,pars=c("beta[1]","beta[2]","beta[3]"))

mcmc_trace(data.mcmc2, 
           pars=c("beta[1]","beta[2]","beta[3]"),
           facet_args = list(ncol = 1, strip.position = "left"))

3.2.5 Rstanarm

# STANARM

library(rstanarm)

## Loading required package: Rcpp

## rstanarm (Version 2.17.4, packaged: 2018-04-13 01:51:52 UTC)

## - Do not expect the default priors to remain the same in future rstanarm versions.

## Thus, R scripts should specify priors explicitly, even if they are just the defaults.

## - For execution on a local, multicore CPU with excess RAM we recommend calling

## options(mc.cores = parallel::detectCores())

## - Plotting theme set to bayesplot::theme_default().

library(broom)

mod.rstanarm = stan_glm(Thickness ~ Diagnostic - 1, 
                             data = df, 
                             iter = 7000, 
                             warmup = 1000,
                             chains = 2, thin = 2, 
                             refresh = 0, 
                             prior = normal(0, 10), 
                             prior_aux = cauchy(0, 5))

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

summary(mod.rstanarm)

## 
## Model Info:
## 
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      Thickness ~ Diagnostic - 1
##  algorithm:    sampling
##  priors:       see help('prior_summary')
##  sample:       6000 (posterior sample size)
##  observations: 38
##  predictors:   3
## 
## Estimates:
##                 mean   sd   2.5%   25%   50%   75%   97.5%
## DiagnosticE    0.3    0.0  0.2    0.3   0.3   0.3   0.4   
## DiagnosticF    0.9    0.0  0.8    0.9   0.9   0.9   1.0   
## DiagnosticN    0.6    0.0  0.5    0.6   0.6   0.6   0.7   
## sigma          0.1    0.0  0.1    0.1   0.1   0.1   0.2   
## mean_PPD       0.6    0.0  0.5    0.5   0.6   0.6   0.6   
## log-posterior 21.9    1.5 18.1   21.2  22.2  23.0  23.7   
## 
## Diagnostics:
##               mcse Rhat n_eff
## DiagnosticE   0.0  1.0  5649 
## DiagnosticF   0.0  1.0  5384 
## DiagnosticN   0.0  1.0  5364 
## sigma         0.0  1.0  5544 
## mean_PPD      0.0  1.0  5793 
## log-posterior 0.0  1.0  4234 
## 
## 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).

tidyMCMC(mod.rstanarm, 
         conf.int = TRUE, 
         conf.method = "HPDinterval")

plot(mod.rstanarm)

library(coda)

data.mcmc3 = as.array(mod.rstanarm)

mcmc_areas(data.mcmc3,  
           pars = c("DiagnosticE", "DiagnosticF", "DiagnosticN"))

mcmc_trace(data.mcmc3, 
           pars = c("DiagnosticE", "DiagnosticF", "DiagnosticN"),
           facet_args = list(ncol = 1, strip.position = "left"))

3.2.6 brms

library(brms)

## Loading 'brms' package (version 2.4.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## Run theme_set(theme_default()) to use the default bayesplot theme.

## 
## Attaching package: 'brms'

## The following objects are masked from 'package:rstanarm':
## 
##     exponential, kfold, lasso, ngrps

library(coda)

mod.BRMS = brm(Thickness ~ Diagnostic -1, 
               data = df, 
               iter = 7000, 
               warmup = 1000,
                    chains = 2, 
                    thin = 2, 
               refresh = 0, 
               prior = c(prior(normal(0, 100), class = "b"), 
                        prior(cauchy(0,5), class = "sigma")
                        )
               )

## Compiling the C++ model

## Start sampling

## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## 
##  Elapsed Time: 0.038 seconds (Warm-up)
##                0.224 seconds (Sampling)
##                0.262 seconds (Total)
## 
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## 
##  Elapsed Time: 0.032 seconds (Warm-up)
##                0.156 seconds (Sampling)
##                0.188 seconds (Total)

summary(mod.BRMS)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Thickness ~ Diagnostic - 1 
##    Data: df (Number of observations: 38) 
## Samples: 2 chains, each with iter = 7000; warmup = 1000; thin = 2;
##          total post-warmup samples = 6000
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## DiagnosticE     0.31      0.03     0.25     0.37       6000 1.00
## DiagnosticF     0.92      0.04     0.84     1.00       5540 1.00
## DiagnosticN     0.59      0.03     0.53     0.65       5865 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     0.12      0.02     0.09     0.15       5403 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

data.mcmc4 = as.mcmc(mod.BRMS)

mcmc_areas(data.mcmc4,  
           pars = c("b_DiagnosticE", "b_DiagnosticF", "b_DiagnosticN"))

mcmc_trace(data.mcmc4, 
           pars = c("b_DiagnosticE", "b_DiagnosticF", "b_DiagnosticN"),
           facet_args = list(ncol = 1, strip.position = "left"))

3.3 Dựng mô hình ANOVA Bayes với greta và TensorFlow

Tiến trình dựng mô hình bằng GRETA bắt đầu bằng việc chuẩn bị (khai báo) 2 bộ phận dữ liệu là biến kết quả (Thickness) và design matrix (ma trận gồm 3 biến giả - dummy variable – yếu tố phân nhóm : Diagnostic = E, F và N). Điểm đặc biệt, đó là chúng ta sẽ dùng TensorFlow, nên vector Thickness và design matrix cần được chuyển đổi thành các tensor (hay array). Việc chuyển đổi này được làm với hàm as_data của greta :

Hàm as_data chuyển design matrix thành 2D tensor, và vector Thickness thành 1D tensor

Ghi chú: bước chuyển dạng dữ liệu có thể không cần thiết, bạn hoàn toàn có thể tính toán với R data object, chúng sẽ được tự động chuyển thành tensor trong khi tính toán, tuy nhiên việc dung hàm as_data( ) có nhiều ích lợi, vì nó kiểm tra các bất thường về dữ liệu (NA, Inf…) và tập thói quen làm việc có trình tự, thứ bậc. Lưu ý: GRETA phân biệt <- và = , toán tử <- để tạo ra các thành phần xác định (deterministic), còn = cho các giá trị bất định (stochastic). Toán tử = tương đương với ~ trong STAN

library(greta)

## 
## Attaching package: 'greta'

## The following objects are masked from 'package:brms':
## 
##     bernoulli, categorical, exponential, lognormal, mixture,
##     student, weibull

## The following objects are masked from 'package:rstanarm':
## 
##     cauchy, dirichlet, exponential, laplace, normal

## The following object is masked from 'package:coda':
## 
##     mcmc

## 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

# Design matrix
design_mat<-model.matrix(~Diagnostic-1, data=df)

head(design_mat,5)

##   DiagnosticE DiagnosticF DiagnosticN
## 1           1           0           0
## 2           1           0           0
## 3           1           0           0
## 4           1           0           0
## 5           1           0           0

Bước tiếp theo trong quy trình, ta sẽ thiết lập các tham số trong mô hình (chú ý: GRETA gọi chúng là Variable, nhưng Nhi sẽ dùng thuật ngữ tham số (parameter), để nhất quán với các ngôn ngữ Bayes khác). Điểm khác biệt giữa GRETA và các ngôn ngữ như STAN, JAGs… đó là trong khi bạn viết code, các object (biến, tensor, tham số ) sẽ xuất hiện ngay lập tức ở thời gian thực, chứ không chờ đến khi toàn bộ code được đóng gói (compile). Ở đây, 2 tham số beta và sigma sẽ được tạo ra dưới dạng object rỗng không (chưa có dữ liệu bên trong), đặc biệt GRETA không yêu cầu bạn phải khai báo chi tiết đặc tính các object này như STAN.

Ghi chú: Thực ra, sự phân biệt Variable và Parameter cũng chính là phân biệt giữa 2 trường phái: frequentist và Bayes; và GRETA cho phép thực hiện mô hình theo cả 2 phái : Nếu bạn theo phái Frequentist, beta, intercept, sigma … được xem là các variable, tạo ra bằng hàm variable( ) và không cần phải có prior (có thể truncate), chúng đơn giản là những hộp rỗng không và nội dung sẽ định hình duy nhất dựa vào dữ liệu và algorithm hồi quy tuyến tính. Còn ngược lại, nếu bạn theo phái Bayes, thì những hộp rỗng này bắt buộc phải có prior, và được xác định bằng các chuỗi MCMC theo định lý Bayes. Variable hay Parameter có thể là 1D (sigma) hay 2D array (beta)

# Para and Priors
betas <- normal(0, 10, dim = ncol(design_mat))
sigma <- cauchy(0, 5, truncation = c(0, Inf))

betas

## greta array (variable following a normal distribution)
## 
##      [,1]
## [1,]  ?  
## [2,]  ?  
## [3,]  ?

sigma

## greta array (variable following a cauchy distribution)
## 
##      [,1]
## [1,]  ?

Một khi các variable/parameter đã có mặt, ta có thể tính toán với chúng, và từ đó tạo ra them những parameter mới, thí dụ ở đây là ước tính Mu bằng phép nhân ma trận giữa beta và design_matrix. Kết quả tạo ra cũng là hộp rỗng không, chưa có nội dung bên trong. Ghi chú: việc tính toán này bao gồm khả năng viết các hàm đặc biệt nếu bạn thích – do đó ngôn ngữ GRETA hoàn toàn có tiềm năng tương đương để tạo ra những mô hình phức tạp như STAN hay JAGs. Tuy nhiên, GRETA linh hoạt hơn rất nhiều so với STAN, vì tất cả những gì bạn viết đều là mã R tự nhiên – có nghĩa là bạn có thể gọi các hàm từ R package khác, viết vòng lặp… tùy thích, chứ không bị trói buộc vào ngôn ngữ bên ngoài như STAN và JAGs (mô hình GRETA không cần phải là 1 khối code hoàn chỉnh và nhất quán như STAN code hay JAGs code, bạn có thể phân chia thành nhiều mảnh, và tất cả đều là R code, R object.

# Mu para
mu <- design_mat %*% betas

head(mu,5)

## greta array (operation)
## 
##      [,1]
## [1,]  ?  
## [2,]  ?  
## [3,]  ?  
## [4,]  ?  
## [5,]  ?

Bước thứ 3 trong quy trình có tính chất bản lề, đó là hàm likelihood. Nội dung của bước này là thiết lập 1 hàm mật độ phân bố cho phép ước lượng giá trị của biến số ngẫu nhiên (kết quả) tùy thuộc vào các tham số của họ phân phối được chọn.

Thí dụ, ở đây là ước lượng Thickness như 1 biến ngẫu nhiên theo quy luật phân bố Gaussian với 2 tham số mu và sigma (đã được tạo ra trước).

# Likelihood 
distribution(df$Thickness) <- normal(mu, sigma)

Bước cuối cùng, ta đóng gói mô hình bằng hàm model( ), lưu ý: chuỗi MCMC sẽ được tạo ra cho mỗi tham số được liệt kê trong hàm này, thí dụ nếu bạn chỉ liệt kê beta (coef), intercept, sd, thì bạn không thể lấy mẫu MCMC cho tham số Mu đến từng cá thể; ngược lại, nếu cho thêm Mu vào model, giá trị Mu cho từng cá thể sẽ được ước lượng nhưng tiến trình lấy mẫu sẽ nặng nhọc hơn nhiều cho máy tính.

# defining the model
m <- model(betas, mu, sigma)

Khác với STAN hay JAGs, đó là Object model mà GRETA vừa tạo ra không phải chỉ là các khái niệm trừu tượng, mà có thực: bạn có thể nhìn thấy sơ đồ cấu trúc của mô hình bằng hàm plot. Điểm thú vị của phân tích Bayes với TensorFlow, đó là sơ đồ này được dung như một mạch những liên kết và hộp chứa rỗng, và những Tensor sẽ chảy qua mạch này theo cùng cơ chế như những mô hình Neural network, để cuối cùng mỗi nút / hộp rỗng (tham só) sẽ được đổ đầy các giá trị của chuỗi MCMC, tùy theo kích thước và hình dạng của chúng.

plot(m)

Như vậy việc dựng mô hình Bayes với GRETA gồm 2 công đoạn chính: Thiết kế cấu trúc mô hình, và lấy mẫu cho MCMC, ta đã làm xong giai đoạn 1, bây giờ đến phần tạo MCMC. Sampler được sử dụng là loại Hamiltonian, với tùy chỉnh chế độ lấy mẫu như STAN và JAG, tức là bạn có thể điều chỉnh số lượt khởi động, số lượt lấy mẫu chính thức, tỉ lệ rút gọn (thinning), số chuỗi MCMC.

Những gì xảy ra tiếp theo tùy thuộc vào máy tính của bạn (CPU hay GPU ? Bao nhiêu CPU ?…), đặc biệt TensorFlow không giới hạn số CPU như rstan hay brms, bạn có thể sampling với CPU cluster. Nhi sử dụng GPU và tốc độ khá nhanh.

draws <- mcmc(m, 
              chains=3,
              warmup = 2000,
              n_samples = 4000,             
              thin=2,
              pb_update = 100)

## 
## chain 1/3
## 
## chain 2/3
## 
## chain 3/3

Sau khi tiến trình lấy mẫu hoàn thành, những bước tiếp theo như khai thác kết quả, suy diễn Bayes… sẽ được làm tương tự như những gì chúng ta đã bàn trong nhiều bài trước. Tất cả dựa vào nội dung các chuỗi MCMC.

summary(draws)

## 
## Iterations = 1:2000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 2000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##              Mean      SD  Naive SE Time-series SE
## betas[1,1] 0.3075 0.03198 0.0004128      0.0005882
## betas[2,1] 0.9238 0.04076 0.0005262      0.0008908
## betas[3,1] 0.5901 0.03156 0.0004074      0.0006952
## mu[1,1]    0.3075 0.03198 0.0004128      0.0005882
## mu[2,1]    0.3075 0.03198 0.0004128      0.0005882
## mu[3,1]    0.3075 0.03198 0.0004128      0.0005882
## mu[4,1]    0.3075 0.03198 0.0004128      0.0005882
## mu[5,1]    0.3075 0.03198 0.0004128      0.0005882
## mu[6,1]    0.3075 0.03198 0.0004128      0.0005882
## mu[7,1]    0.3075 0.03198 0.0004128      0.0005882
## mu[8,1]    0.3075 0.03198 0.0004128      0.0005882
## mu[9,1]    0.3075 0.03198 0.0004128      0.0005882
## mu[10,1]   0.3075 0.03198 0.0004128      0.0005882
## mu[11,1]   0.3075 0.03198 0.0004128      0.0005882
## mu[12,1]   0.3075 0.03198 0.0004128      0.0005882
## mu[13,1]   0.3075 0.03198 0.0004128      0.0005882
## mu[14,1]   0.3075 0.03198 0.0004128      0.0005882
## mu[15,1]   0.9238 0.04076 0.0005262      0.0008908
## mu[16,1]   0.9238 0.04076 0.0005262      0.0008908
## mu[17,1]   0.9238 0.04076 0.0005262      0.0008908
## mu[18,1]   0.9238 0.04076 0.0005262      0.0008908
## mu[19,1]   0.9238 0.04076 0.0005262      0.0008908
## mu[20,1]   0.9238 0.04076 0.0005262      0.0008908
## mu[21,1]   0.9238 0.04076 0.0005262      0.0008908
## mu[22,1]   0.9238 0.04076 0.0005262      0.0008908
## mu[23,1]   0.9238 0.04076 0.0005262      0.0008908
## mu[24,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[25,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[26,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[27,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[28,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[29,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[30,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[31,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[32,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[33,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[34,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[35,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[36,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[37,1]   0.5901 0.03156 0.0004074      0.0006952
## mu[38,1]   0.5901 0.03156 0.0004074      0.0006952
## sigma      0.1195 0.01515 0.0001956      0.0003342
## 
## 2. Quantiles for each variable:
## 
##               2.5%    25%    50%    75%  97.5%
## betas[1,1] 0.24319 0.2859 0.3078 0.3285 0.3704
## betas[2,1] 0.84084 0.8977 0.9241 0.9509 1.0039
## betas[3,1] 0.52632 0.5695 0.5902 0.6112 0.6509
## mu[1,1]    0.24319 0.2859 0.3078 0.3285 0.3704
## mu[2,1]    0.24319 0.2859 0.3078 0.3285 0.3704
## mu[3,1]    0.24319 0.2859 0.3078 0.3285 0.3704
## mu[4,1]    0.24319 0.2859 0.3078 0.3285 0.3704
## mu[5,1]    0.24319 0.2859 0.3078 0.3285 0.3704
## mu[6,1]    0.24319 0.2859 0.3078 0.3285 0.3704
## mu[7,1]    0.24319 0.2859 0.3078 0.3285 0.3704
## mu[8,1]    0.24319 0.2859 0.3078 0.3285 0.3704
## mu[9,1]    0.24319 0.2859 0.3078 0.3285 0.3704
## mu[10,1]   0.24319 0.2859 0.3078 0.3285 0.3704
## mu[11,1]   0.24319 0.2859 0.3078 0.3285 0.3704
## mu[12,1]   0.24319 0.2859 0.3078 0.3285 0.3704
## mu[13,1]   0.24319 0.2859 0.3078 0.3285 0.3704
## mu[14,1]   0.24319 0.2859 0.3078 0.3285 0.3704
## mu[15,1]   0.84084 0.8977 0.9241 0.9509 1.0039
## mu[16,1]   0.84084 0.8977 0.9241 0.9509 1.0039
## mu[17,1]   0.84084 0.8977 0.9241 0.9509 1.0039
## mu[18,1]   0.84084 0.8977 0.9241 0.9509 1.0039
## mu[19,1]   0.84084 0.8977 0.9241 0.9509 1.0039
## mu[20,1]   0.84084 0.8977 0.9241 0.9509 1.0039
## mu[21,1]   0.84084 0.8977 0.9241 0.9509 1.0039
## mu[22,1]   0.84084 0.8977 0.9241 0.9509 1.0039
## mu[23,1]   0.84084 0.8977 0.9241 0.9509 1.0039
## mu[24,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[25,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[26,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[27,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[28,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[29,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[30,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[31,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[32,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[33,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[34,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[35,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[36,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[37,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## mu[38,1]   0.52632 0.5695 0.5902 0.6112 0.6509
## sigma      0.09442 0.1089 0.1179 0.1286 0.1536

c1<-draws[[1]]%>%as.data.frame()%>%mutate(Chain=1,iteration=c(1:2000))
c2<-draws[[2]]%>%as.data.frame()%>%mutate(Chain=2,iteration=c(1:2000))
c3<-draws[[3]]%>%as.data.frame()%>%mutate(Chain=3,iteration=c(1:2000))

pdf=rbind(c1,c2,c3)%>%dplyr::select(c(1:3,43,44))

colnames(pdf)<-c("Emphysema","Fibrosis","Normal","Chain","Iteration")

pdf%>%gather(Emphysema:Normal,key="Para",value="Value")%>%
  ggplot(aes(x=Value,fill=factor(Chain),col=factor(Chain)))+
  geom_density(alpha=0.3)+
  theme_bw()+
  facet_wrap(~Para,ncol=1,scales="free")

pdf%>%gather(Emphysema:Normal,key="Para",value="Value")%>%
  ggplot(aes(x=Iteration,y=Value,fill=factor(Chain),col=factor(Chain)))+
  geom_path(alpha=0.5)+
  theme_bw()+
  facet_wrap(~Para,ncol=1,scales="free")

pdf%>%gather(Emphysema:Normal,key="Para",value="Value")%>%
  group_by(Para)%>%
  summarise_at("Value",funs( Median=median,
                         LL=quantile(.,probs=0.025),
                         UL=quantile(.,probs=0.975)))

Ta có thể kiểm tra lại nội dung mô hình bằng hàm glm thông thường…

glm(Thickness~Diagnostic-1,data=df)%>%summary()

## 
## Call:
## glm(formula = Thickness ~ Diagnostic - 1, data = df)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.21078  -0.08879  -0.01242   0.08684   0.18183  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## DiagnosticE  0.30750    0.03078    9.99 8.69e-12 ***
## DiagnosticF  0.92296    0.03839   24.04  < 2e-16 ***
## DiagnosticN  0.59098    0.02974   19.87  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.01326313)
## 
##     Null deviance: 14.69352  on 38  degrees of freedom
## Residual deviance:  0.46421  on 35  degrees of freedom
## AIC: -51.551
## 
## Number of Fisher Scoring iterations: 2

4 Nhận xét:

GRETA là một ý tưởng mới để thực hiện thống kê Bayes trong R, có thể xem đây là thế hệ ngôn ngữ lập trình Bayes thứ 4 (sau BUGs, JAGs và STAN). GRETA có 3 ưu điểm quan trọng so với những phương pháp đi trước, bao gồm : 1) Sử dụng trực tiếp R code tự nhiên, không phụ thuộc vào một ngôn ngữ bên ngoài như BUGs, STAN và JAGs. Ngay cả các giao thức như mcmcPack, rstanarm, jagR, brms … cũng không làm được điều này, vì mô hình vẫn phải được đóng gói và phiên dịch thành mã JAGs hay STAN, và bị giới hạn bởi những tính năng, họ phân phối và hàm của các ngôn ngữ này. Với GRETA, ta có thể làm việc độc lập và hoàn toàn bằng R.

2. GRETA cho phép viết mô hình từng bước, và giữa các bước độc lập với nhau. Cú pháp khá đơn giản, dễ hiểu.

3. Việc sử dụng TensorFlow cho phép làm việc với những dữ liệu kích thước lớn, phức tạp, bao gồm Time series hoặc ảnh chụp, ngôn ngữ…

Có thể tiên đoán rằng GRETA sẽ thành một ngôn ngữ Bayes nổi bật không kém gì STAN trong vài năm sắp tới.