Bayesian

The Bayesian Approach

Prior and Posterior Distributions

Example: Beta-Binomial

Prior on \(\pi\)

Mean of posterior distribution

Example: Belief in Hell

qbeta(c(0.025,0.975),814.5,328.5) # quantiles, beta hyperparameters

[1] 0.6860285 0.7384637

library(proportion) #can find posterior intervals directly
ciBAx(814,1142,0.05,0.5, 0.5) #y, n, error probability, beta hyperparameters

    x     LBAQx     UBAQx     LBAHx     UBAHx
1 814 0.6860285 0.7384637 0.6862811 0.7387072

ciAllx(814,1142,0.05)

      method   x LowerLimit UpperLimit LowerAbb UpperAbb ZWI
1       Wald 814  0.6865425  0.7390267       NO       NO  NO
2    ArcSine 814  0.6861994  0.7386542       NO       NO  NO
3 Likelihood 814  0.6861089  0.7385187       NO       NO  NO
4      Score 814  0.6858635  0.7382790       NO       NO  NO
5 Logit-Wald 814  0.6858446  0.7382959       NO       NO  NO
6     Wald-T 814  0.6865302  0.7390390       NO       NO  NO

Bayesian Posterior Interval Comparing Proportions

library(PropCIs)
diffci.bayes(11,11,0, 1, 1.0,1.0,1.0,1.0,0.95,nsim=1000000)

[1] 0.06424402 0.94965628

# arguments y1,n1, y2, n2, alpha1, beta1, alpha2, beta2, confidence, no.simulations

Highest Posterior Density (HPD) Posterior Intervals

Credible interval - Wikipedia

Recall the clinical trial example given above. We next obtain the HPD interval by intensive simulation:

library(HDInterval)

pi1<-rbeta(10000000,12.0,1.0) # posterior for 11 successes, 0 failures, uniform prior
pi2<-rbeta(10000000,1.0,2.0) # posterior for 0 success, 1 failure, uniform prior
hdi(pi1-pi2,credMass=0.95)

    lower     upper 
0.1255657 0.9813333 
attr(,"credMass")
[1] 0.95

hist(pi1-pi2) # posterior (not shown) is highly skewed to left

plot(density(pi1-pi2)) # density approximation

Bayesian Inference for Means

Example: Bayesian Analysis for Anorexia Therapy

(see Section 4.8.1 for more details on the inverse Gamma prior)

Anorexia<-read.table("http://stat4ds.rwth-aachen.de/data/Anorexia.dat", header=TRUE)
change<-Anorexia$after-Anorexia$before

library(MCMCpack)

Loading required package: coda

Loading required package: MASS

##
## Markov Chain Monte Carlo Package (MCMCpack)

## Copyright (C) 2003-2023 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park

##
## Support provided by the U.S. National Science Foundation

## (Grants SES-0350646 and SES-0350613)
##

fit<-MCMCregress(change[Anorexia$therapy=="cb"]~1,mcmc=5000000,
b0=0,B0=10^{-15},c0=10^{-15},d0=10^{-15})
# mean has normal prior dist. with mean b0=0, variance 1/B0 (std.dev. > 31 million)
# variance has inverse gamma prior distribution (c0/2=shape, d0/2=scale)
summary(fit)

Iterations = 1001:5001000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5e+06 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

              Mean     SD  Naive SE Time-series SE
(Intercept)  3.008  1.408 0.0006297      0.0006297
sigma2      57.531 16.626 0.0074352      0.0077093

2. Quantiles for each variable:

             2.5%   25%    50%    75%  97.5%
(Intercept)  0.23  2.08  3.007  3.936  5.787
sigma2      33.63 45.85 54.719 66.013 97.749

An improper prior

y<-Anorexia$after[Anorexia$therapy=="cb"]-Anorexia$before[Anorexia$therapy=="cb"]
n=length(y)
S=sum((y-mean(y))^2) # this is (n-1)s^2
rsigma2<-S/rchisq(1000000,n-1) # million random variables from posterior of sigma^2
mu<-rnorm(1000000,mean=mean(y),sd=sqrt(rsigma2)/sqrt(n)) # random normal means
cbind(n,mean(mu),sd(mu)) # mean and standard deviation of posterior distribution

      n                  
[1,] 29 3.004924 1.405974

quantile(mu,c(0.025,0.975)) # 95% posterior interval for population mean

     2.5%     97.5% 
0.2297127 5.7794269 

y1<-Anorexia$after[Anorexia$therapy=="cb"]-Anorexia$before[Anorexia$therapy=="cb"]
y2<-Anorexia$after[Anorexia$therapy=="c"]-Anorexia$before[Anorexia$therapy=="c"]
n1<-length(y1);n2<-length(y2)
S=sum((y1-mean(y1))^2)+ sum((y2 -mean(y2))^2) # this is [(n1-1)+(n2-1)]s^2
rsigma2<-S/rchisq(1000000,n1+n2 - 2) # random from posterior of sigma^2
mu1<-rnorm(1000000,mean=mean(y1),sd=sqrt(rsigma2)/sqrt(n1)) # random normal means
mu2<-rnorm(1000000,mean=mean(y2),sd=sqrt(rsigma2)/sqrt(n2))
cbind(n1,n2,mean(mu1-mu2),sd(mu1-mu2))

     n1 n2                  
[1,] 29 26 3.457536 2.104257

quantile(mu1-mu2,c(0.025, 0.975)) # 95% posterior interval for mu1-mu2

      2.5%      97.5% 
-0.6778715  7.6017025 

(see Section 4.8.3 for details on the improper prior)

Why Maximum Likelihood and Bayes Estimators Perform Well? - Both Enjoy Good Large Sample Properties!

(see Section 4.9 for mathematical details)

Exercises

(see Dirichlet distribution - Wikipedia)