Workshop 5: Introduction to Bayesian models

Murray Logan

09 Apr 2016

Frequentist

P(D|H)

long-run frequency

simple analytical methods to solve roots

conclusions pertain to data, not parameters or hypotheses

compared to theoretical distribution when NULL is true

probability of obtaining observed data or MORE EXTREME data

\(P(D|H)\)

can a null ever actually be true

Frequentist

P-value
- probabulity of rejecting NULL
- NOT a measure of the magnitude of an effect or degree of significance!
- measure of whether the sample size is large enough

95% CI
- NOT about the parameter it is about the interval
- does not tell you the range of values likely to contain the true mean

Frequentist vs Bayesian

-------------------------------------------------
                Frequentist          Bayesian  
--------------  ------------         ------------
Obs. data       One possible         Fixed, true

Parameters      Fixed, true          Random, 
                                     distribution

Inferences      Data                 Parameters

Probability     Long-run frequency   Degree of belief
                $P(D|H)$             $P(H|D)$
-------------------------------------------------

Frequentist vs Bayesian

n: 10
Slope: -0.1022
t: -2.3252
p: 0.0485

n: 10
Slope: -10.2318
t: -2.2115
p: 0.0579

n: 100
Slope: -10.4713
t: -6.6457
p: 1.7101362 × 10^-9

Frequentist vs Bayesian

	Population A	Population B
Percentage change	0.46	45.46
Prob. >5% decline	0	0.86

Bayesian

Bayes rule

\[ \begin{aligned} P(H|D) &= \frac{P(D|H) \times P(H)}{P(D)}\\[1em] \mathsf{posterior\\belief\\(probability)} &= \frac{likelihood \times \mathsf{prior~probability}}{\mathsf{normalizing~constant}} \end{aligned} \]

Bayesian

Bayes rule

\[ \begin{aligned} P(H|D) &= \frac{P(D|H) \times P(H)}{P(D)}\\ \mathsf{posterior\\belief\\(probability)} &= \frac{likelihood \times \mathsf{prior~probability}}{\mathsf{normalizing~constant}} \end{aligned} \]

The normalizing constant is required for probability - turn a frequency distribution into a probability distribution

Estimation: OLS

Estimation: Likelihood

\(P(D|H)\)

Bayesian

conclusions pertain to hypotheses
computationally robust (sample size,balance,collinearity)
inferential flexibility - derive any number of inferences

Bayesian

subjectivity?
intractable \[P(H|D) = \frac{P(D|H) \times P(H)}{P(D)}\] \(P(D) \)- probability of data from all possible hypotheses

MCMC sampling

Marchov Chain Monte Carlo sampling

draw samples proportional to likelihood

 <ul>
<li>two parameters $\alpha$ and $\beta$</li>
<li>infinitely vague priors - posterior likelihood only</li>
<li>likelihood multivariate normal</li>

MCMC sampling

Marchov Chain Monte Carlo sampling

draw samples proportional to likelihood

 <ul>
<li>two parameters $\alpha$ and $\beta$</li>
<li>infinitely vague priors - posterior likelihood only</li>
<li>likelihood multivariate normal</li>

MCMC sampling

Marchov Chain Monte Carlo sampling

draw samples proportional to likelihood

MCMC sampling

Marchov Chain Monte Carlo sampling

chain of samples

MCMC sampling

Marchov Chain Monte Carlo sampling

1000 samples

MCMC sampling

Marchov Chain Monte Carlo sampling

10,000 samples

MCMC sampling

Marchov Chain Monte Carlo sampling

Aim: samples reflect posterior frequency distribution
samples used to construct posterior prob. dist.
the sharper the multidimensional “features” - more samples
chain should have traversed entire posterior
inital location should not influence

MCMC diagnostics

Trace plots

MCMC diagnostics

Autocorrelation

Summary stats on non-independent values are biased
Thinning factor = 1

MCMC diagnostics

Autocorrelation

Summary stats on non-independent values are biased
Thinning factor = 10

MCMC diagnostics

Autocorrelation

Summary stats on non-independent values are biased
Thinning factor = 10, n=10,000

MCMC diagnostics

Plot of Distributions

Native options in R

MCMCpack
MCMCglmm

JAGS/BUGS

WinBUGS - object pascal
- made Bayesian analyses ‘available’ to masses
- models mirror written definitions
- very slow

JAGS - c++
- same declarative language
- much faster

JAGS/BUGS

stand along application

file input
- model declaration
- data as a list

R2jags - interface from R

JAGS/BUGS

\(y_i = \beta_0 + \beta_1 x_i+\epsilon_i \quad \epsilon \sim{} N(0,\sigma^2)\)

\(y_i \sim{} N(\beta_0 + \beta_1 x_i, \sigma^2)\)

\(y_i \sim{} N(\beta_0 + \beta_1 x_i, \tau)\)
- \(\tau\) is precision (\(\frac{1}{\sigma^2}\))

\(y_i \sim{} N(\mu_i, \tau)\)
\(\mu_i = \beta_0 + \beta_1 x_i\)

\(\beta_0 \sim{} N(0,0.000001)\)

\(\beta_1 \sim{} N(0,0.000001)\)

\(\tau = \frac{1}{\sigma^2}\)

\(\sigma \sim{} Uniform(0,100)\)

JAGS/BUGS

Define the model

> modelString="
+ model {
+    #Likelihood
+    for (i in 1:n) {
+       y[i]~dnorm(mu[i],tau)
+       mu[i] <- beta0+beta1*x[i]
+    }
+ 
+    #Priors
+    beta0 ~ dnorm (0,1.0E-6)
+    beta1 ~ dnorm(0,1.0E-6)
+    tau <- 1 / (sigma * sigma)
+    sigma~dunif(0,100)
+  }
+ "

\(y_i \sim{} N(\mu_i, \tau)\)
\(\mu_i = \beta_0 + \beta_1 x_i\)

\(\beta_0 \sim{} N(0,0.000001)\)
\(\beta_1 \sim{} N(0,0.000001)\)
\(\tau = \frac{1}{\sigma^2}\)
\(\sigma \sim{} Uniform(0,100)\)

> writeLines(modelString,con="BUGSscripts/regression.txt")

Error in file(con, "w"): cannot open the connection

JAGS/BUGS

Create the data list

Y	X
3	0
2.5	1
6	2
5.5	3
9	4
8.6	5
12	6

JAGS/BUGS

Create the data list

> data.list <- with(DATA,
+            list(y=Y, 
+            x=X,n=nrow(DATA))
+ )
> data.list

$y
[1]  3.0  2.5  6.0  5.5  9.0  8.6 12.0

$x
[1] 0 1 2 3 4 5 6

$n
[1] 7

JAGS/BUGS

Define the chain parameters

> #params <- c("beta0","beta1","sigma")
> #burnInSteps = 2000
> #nChains = 3 
> #numSavedSteps = 50000
> #thinSteps = 1
> #nIter = ceiling((numSavedSteps * thinSteps)/nChains)

JAGS/BUGS

Perform MCMC sampling

> library(R2jags)