Chapter 17: Bayesian Analysis

Welcome to the online content for Chapter 17!

As always, I’ll assume that you’ve already read up to this chapter of the book and worked through the online content for the previous chapters. If not, please do that first.

As always, click the ‘Run Code’ buttons below to execute the R code. Remember to wait until they say ‘Run Code’ before you press them. And be careful to run these boxes in order if later boxes depend on you having done other things previously.

Using BayesFactor

The package BayesFactor that we need for Bayesian hypothesis testing doesn’t yet work on WebR, so to use the code below you’ll need to download and install R on your own computer, as described in the Introduction.

When you’ve done that, the code below will enable you to run an analysis using the data from the chapter.

Lines that begin with a # symbol are just comments, to explain what’s going on, and don’t actually do anything.

#Install the BayesFactor package
install.packages("BayesFactor")
#Load the library
library(BayesFactor)

#Read in the data from the chapter
people1 <- read.csv("https://colinfoster77.github.io/Stats/Chap17Hypothesistesting.csv")
people1

#You can see that the data consist of 20 people's heights.

# Conventional t test
t.test(people1$height, mu = 170, alternative = "two.sided")

#The mean height is significantly different from 170 cm, with p = .012.

#Bayesian t test
ttestBF(x = people1$height-170)

#The Bayes factor of 4.294801 matches the value given in the chapter.

Using rstan

The package rstan that we need for Bayesian analysis doesn’t yet work on WebR, so to use the code below you’ll need to download and install R on your own computer, as described in the Introduction.

When you’ve done that, the code below will enable you to run an analysis using the data from the chapter.

Lines that begin with a # symbol are just comments, to explain what’s going on, and don’t actually do anything.

#Install the rstan package
install.packages("rstan")
#Load the library
library(rstan)

#Read in the data from the chapter
people2 <- read.csv("https://colinfoster77.github.io/Stats/Chap17Heights.csv")
people2

#You can see that the data consist of another 20 people's heights.

#Put all of the heights into a vector called 'y'
y <- people2$height
#Let 'n' be the number of people, 20 (i.e. the length of the vector 'y')
n <- length(people2$height)

#Define the Stan model code
#Don't worry too much about this! It just tells R to use a Normal distribution and sets the prior values.
model.code <- '
data {
  int<lower=0> n;
  vector[n] y;
}

parameters {
  real mu;
  real<lower=0> sigma;
}

model {
  y ~ normal(mu, sigma);
  mu ~ normal(170,20);
  sigma ~ cauchy(6,50);
}
 '

# Compile the Stan model - this will take a few seconds, so be patient!
model <- stan_model(model_code = model.code)

fit = sampling(model,list(n=n,y=y),iter=1000,chains=4)

#Look at the results
print(fit)

#Get the parameters
parameters = extract(fit)

#Show the distributions of these the posterior mean and standard deviation, putting labels on the x axis to show which is which
hist(parameters$mu, main="", xlab="posterior mean")
hist(parameters$sigma, main="", xlab="posterior standard deviation")

#Print out the posterior distribution mean and standard deviation
mean(parameters$mu)
mean(parameters$sigma)