Stats Week 5 lecture 2
lecture_notes_week5_2.pdf

In R

We are using the following packages for Bayesian inference

RStan

RStan is an implementation of the Stan software used for Bayesian Analysis into R

Stan gives tools which builds (compiles) a C++ program into your computer each time you run it. Unlike with most regression packages which just load the functions and let you run it. But with Stan the program needs to be built every time you generate a new model. This usually takes a couple minutes so be patient.

BRMS

BRMS implements "Bayesian Regression Models using Stan" - used for our regression analysis.

When you call a brms function, it starts by building the Stan program, then implements an algorithm to do the Markov chain Monte Carlo sampling.

Example

After first loading rstan and brms, we implement a Bayesian regression model with brm() on some data. After compiling the Stan program, it begins the Markov chain Monte Carlo sampling to get a posterior distribution sample for each model.

Sampling process:

1. Warmup

• It picks a first guess for the parameter models
• With a set of rules, it picks a random set for each function, then starts moving around
• Sometimes it gets stuck - sometimes it moves fast, sometimes slower
• It eventually approaches the point of highest probability distribution for correlation of the two variables

2. Sampling

• After it has got close to the highest likelihood zone, it starts to take random samples from the posterior distribution
• If you plot the lines between these samples in R, you can see that it zig-zags randomly around the model space
  • This shows that the sampels are independent and random which is good

This is the process of one "chain" in the Markov Chain. By default, it samples 4 chains to make sure they converge to the same place. This is a check to make sure none of the chain is falling into some local mean (there are also other complex processes to prevent this).

summary(brm)

We will compare the outputs of this summary to the linear model summary values in R:

Start

• The beginning reminds us of the call, formulas, datas, family, link functions etc

Draws

• This tells us about the number of chains used and their parameters:
  • chains: Number of Markov chains used, default is 4
  • iter: Number of total iterations including warmup used
  • warmup: The number of iterations that were used as part of the warmup. These warmup values are not actually used in model inference.
  • thin: Thinning rate - I don't really understand this, but setting this to a higher number I guess thins out the data and saves computation time.

Inference results

Coefficients

• Lists each of the coefficients, gives an estimate and some other parameters:
  • Est. Error:
  • Upper and lower 95% Confidence interval (l-95% CI and u-95% CI)
    • These are pretty much the same as what's given in the linear model with confint()
    • This is sort of a special case in this case, because we used a Gaussian observation model and linear regression in Bayesian space, so it's the same as the linear model reslts
  • Rhat
    • A measure of how well the sampling happened (how much the results can be trusted)
    • A measure of sampling variance ratios between chains to variance within chains
    • Tells us about convergence in the model
    • In our example from class it was exactly 1.00
    • The threshold is still being determined, but if it's above 1.05 then there is probably some trouble and you maybe need more samples
  • Bulk and Tail ESS (Effective Sample Size)
    • Because we have some correlation between our variables, two points that are close to each other contribute less to the model than two points that are far away
    • If we get a result of 3000 for Tail ESS out of 4000 samples, it means that our results are effectively the same as if we used 3000 samples
    • Bulk ESS is a bit complicated and we shouldn't really worry too much about it
    • The larger the autocorrelation between samples is, the lower your ESS values will be

Sigma

• An estimate of the variance (?)

Prior summary

Using the prior_summary() function we can find the prior distribution for each of our parameters.

We don't need to worry too much about this - but it tells us what the default priors used in our function is, and you might want to tweak these if you know what you're doing (not required for this course)

Tweaking the model

Some things we might want to tweak if we are running it again:

• setting the prior probabilities
• choosing the observation model
• setting number of iterations
• Setting a seed used as part of the random sampling

Assessing the model

Assessing the posterior samples

We need to determine if the Markov chain Monte Carlo (MCMC) algorithm has converged well. We need to:

• Visually inspect the sample chains and posterior distributions
• Examine the Rhat statistic - if it's larger than 1.05, it means that there are some convergence problems with MCMC

Visual inspection

• First plot histograms to make sure there is no trends other than being normal
• MCMC plot type trace tells us the probability of each chain ... (?)
  • Everything is good if the plot looks totally random - if there is a clear trend there is a problem
• MCMC plot type acf_bar tells us about autocorrelation in the data for each parameter and chain
  • Autocorrelation does not make the data useless, just reduces the accuracy
  • Not a huge problem with a large MCMC sample size
  • Not exactly sure how to actually interpret this plot

Assessing the model fit

• Plotting the posterior probability using pp_check() - this is like plotting the residuals vs fitted
  • Compares observed data to data simulated from the model
  • If there is a big difference there is a problem
  • Can be done with a histogram, boxplots or scatterplot
• Leave-One-Out cross-validation (loo())
  • Assesses the predictive skill for data that was not used in trainig the model
  • Used to compare models, like AIC
• Prediction and conditional effects
  • ?
• Extract posterior samples for model parameters

Applying to Poisson observation model

• Next up, in the lecture notes, this is applied to the Poisson Observation Model
• Starting on page 20