Bayesian Modelling with Rstan

R Bayesian Statistics Series

An Example

# https://www.google.com/search?q=rstan+example&ei=6ZQ7YMvdJKaQr7wP84mfoA4&start=10&sa=N&ved=2ahUKEwiL9uWn0ozvAhUmyIsBHfPEB-QQ8tMDegQIBRA0&biw=1653&bih=902

#---------------------------------
#  Prepare data for stan program
#---------------------------------
N <- nrow(trees)
x <- trees$Height
y <- trees$Volume

stan_data <- list(N = N, x = x, y = y)

#-----------------------------
#          Stan code
#-----------------------------

stan_model <- "

data {
  int < lower = 1 > N; // Sample size
  vector[N] x; // Predictor
  vector[N] y; // Outcome
}

parameters {
  real alpha; // Intercept
  real beta; // Slope (regression coefficients)
  real < lower = 0 > sigma; // Error SD
}

model {
  y ~ normal(alpha + x * beta , sigma);
   
}"

#--------------------------------------------
#   Fit Bayesian Model with proper priors
#--------------------------------------------

# fit2 <- stan(model_code = stan_model, data = stan_data, iter = 1000, chains = 4)

library(rstan)

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

rstan_options(auto_write = TRUE)

system.time(
  fit <- stan(file = "bayesian_reg.stan", 
              data = stan_data, 
              warmup = 1000, 
              iter = 2000, 
              chains = 4, 
              algorithm = "HMC", 
              verbose = TRUE, 
              seed = 29)
  )

## 
## TRANSLATING MODEL 'bayesian_reg' FROM Stan CODE TO C++ CODE NOW.
## successful in parsing the Stan model 'bayesian_reg'.
## 
## CHECKING DATA AND PREPROCESSING FOR MODEL 'bayesian_reg' NOW.
## 
## COMPILING MODEL 'bayesian_reg' NOW.
## 
## STARTING SAMPLER FOR MODEL 'bayesian_reg' NOW.

##    user  system elapsed 
##    1.39    0.23   16.67

# Results: 

fit

## Inference for Stan model: bayesian_reg.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##         mean se_mean    sd    2.5%     25%    50%    75%  97.5% n_eff Rhat
## alpha -87.12    0.25 30.13 -145.71 -107.03 -87.30 -67.32 -26.44 14408 1.00
## beta    1.54    0.00  0.39    0.76    1.29   1.54   1.80   2.30 14408 1.00
## sigma  13.89    0.03  1.82   10.85   12.64  13.68  14.96  17.88  3364 1.00
## lp__  -93.81    0.05  1.20  -96.93  -94.36 -93.51 -92.93 -92.45   710 1.01
## 
## Samples were drawn using HMC(diag_e) at Sun Feb 28 22:05:27 2021.
## 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).

traceplot(fit, pars = c("alpha", "beta", "sigma"), nrow = 3)

# Compare: 

lm(y ~ x)

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##     -87.124        1.543

# More information: 

fit_ex <- rstan::extract(fit)

mean(fit_ex$alpha)

## [1] -87.12175

length(fit_ex$alpha)

## [1] 4000