Introduction to Bayesian models

Murray Logan

12/02/2013

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

$P(D|H)$
can a null ever actually be true

Frequentist vs Bayesian

	Frequentist	Bayesian
Observed data	One possible	Fixed, true
Parameters	Fixed, true	Random, distribution
Inferences	Data	Parameters
Probability	Long-run frequency P(D\|H)	Degree of belief P(H\|D)

Bayesian

Bayes' rule

\begin{align*} 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{align*}

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

Bayesian

Advantages

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

Bayesian

Disadvantages

subjectivity?
intractable

$$P(H|D) = \frac{P(D|H) \times P(H)}{P(D)}$$ $$P(D) - \mathsf{probability~of~data~from~all~possible~hypotheses}$$

MCMC sampling

Marchov Chain Monte Carlo sampling

draw samples proportional to likelihood

two parameters $\alpha$ and $\beta$
infinitely vague priors - posterior likelihood only
likelihood multivariate normal

MCMC sampling

Marchov Chain Monte Carlo sampling

initial values

two parameters $\alpha$ and $\beta$
infinitely vague priors - posterior likelihood only
likelihood multivariate normal

MCMC sampling

Marchov Chain Monte Carlo sampling

ratio of new and old likelihoods

two parameters $\alpha$ and $\beta$
infinitely vague priors - posterior likelihood only
likelihood multivariate normal

MCMC sampling

Marchov Chain Monte Carlo sampling

chain of samples

two parameters $\alpha$ and $\beta$
infinitely vague priors - posterior likelihood only

MCMC sampling

Marchov Chain Monte Carlo sampling

1000 samples

two parameters $\alpha$ and $\beta$
infinitely vague priors - posterior likelihood only
likelihood multivariate normal

MCMC sampling

Marchov Chain Monte Carlo sampling

10,00 samples

two parameters $\alpha$ and $\beta$
infinitely vague priors - posterior likelihood only
likelihood multivariate normal

MCMC sampling

Marchov Chain Monte Carlo sampling

Aim: samples reflect posterior
the sharper the multidimensional 'features' - more samples
chain should have traversed entire posterior
inital location should not influence

Trace plots

Autocorrelation

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

Autocorrelation

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

Autocorrelation

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

Plots of distributions

Freq vs Bayes


n: `10` Slope: `-0.0792` t: `-2.3399` p: `0.0474`	n: `10` Slope: `-7.8107` t: `-2.1922` p: `0.0597`	n: `100` Slope: `-10.1985` t: `-7.4121` p: `0`

JAGS

$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}$)
- $\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

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)
+  }
+ " 
> writeLines(modelString,con="BUGSscripts/model1.txt")

JAGS

Create a data list

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

$y
[1]  1.374  3.684  4.164  8.095  8.330
[6]  8.680 11.487

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

$n
[1] 7

JAGS

Define chain parameters

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

JAGS

Perform MCMC sampling

> library(R2jags)
> data.r2jags <- jags(data=data.list, 
+ 	  inits=NULL, 
+ 	  parameters.to.save=params, 
+ 	  model.file="BUGSscripts/model1.txt",
+ 	  n.chains=3,
+ 	  n.iter=nIter, 
+ 	  n.burnin=burnInSteps,
+       n.thin=thinSteps
+ 	  )

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 38

Initializing model

JAGS

Evaluate the model - traceplot
> plot(as.mcmc(data.r2jags))

JAGS

Evaluate the model - Raftery diagnostics

> raftery.diag(as.mcmc(data.r2jags))

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                     
          Burn-in  Total Lower bound 
          (M)      (N)   (Nmin)      
 beta0    3        4586  3746        
 beta1    4        7366  3746        
 deviance 3        4054  3746        
 sigma    4        4910  3746        
           
 Dependence
 factor (I)
 1.22      
 1.97      
 1.08      
 1.31      


[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                     
          Burn-in  Total Lower bound 
          (M)      (N)   (Nmin)      
 beta0    9        13827 3746        
 beta1    6        8180  3746        
 deviance 3        4008  3746        
 sigma    4        4773  3746        
           
 Dependence
 factor (I)
 3.69      
 2.18      
 1.07      
 1.27      


[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                     
          Burn-in  Total Lower bound 
          (M)      (N)   (Nmin)      
 beta0    6        8440  3746        
 beta1    6        8192  3746        
 deviance 2        3919  3746        
 sigma    4        5081  3746        
           
 Dependence
 factor (I)
 2.25      
 2.19      
 1.05      
 1.36

JAGS

Evaluate the model - Autocorrelation diagnostics

> autocorr.diag(as.mcmc(data.r2jags))

           beta0     beta1 deviance
Lag 0   1.000000  1.000000  1.00000
Lag 1  -0.008001 -0.001116  0.66897
Lag 5   0.012437  0.002244  0.26645
Lag 10  0.003184  0.005771  0.11655
Lag 50  0.001979  0.004254  0.00254
          sigma
Lag 0  1.000000
Lag 1  0.793155
Lag 5  0.390997
Lag 10 0.203675
Lag 50 0.001013

JAGS

Examine the model parameters

> print(data.r2jags)

Inference for Bugs model at "BUGSscripts/model1.txt", fit using jags,
 3 chains, each with 16667 iterations (first 2000 discarded)
 n.sims = 44001 iterations saved
         mu.vect sd.vect   2.5%    25%
beta0      1.770   1.045 -0.302  1.238
beta1      1.590   0.287  1.023  1.440
sigma      1.348   0.721  0.633  0.915
deviance  21.650   4.111 17.057 18.667
            50%    75%  97.5%  Rhat
beta0     1.773  2.311  3.811 1.001
beta1     1.590  1.739  2.163 1.001
sigma     1.161  1.546  3.226 1.001
deviance 20.539 23.467 32.494 1.001
         n.eff
beta0    38000
beta1    44000
sigma    44000
deviance 44000

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 = 8.5 and DIC = 30.1
DIC is an estimate of expected predictive error (lower deviance is better).