###############################################

## Joshua M. Tebbs

## Date: 22 May 2010

## Update: 9 August 2011

## STAT 513 course notes: R Code

## Chapter 12

###############################################

 

# Example 12.2.

# Beta binomial Bayesian analysis

# Posterior construction

# Page 111-114

# theta is the binomial proportion

# max theta is 0.3 on plots because theta is close to zero

theta<-seq(0,0.3,0.001)

# Parameters in the beta(alpha,beta) prior

alpha=1

beta=19

# prior.theta is a function that constructs the prior distribution

# lgamma is the internal log(gamma) function

prior.theta<-function(theta,alpha,beta){

      exp(lgamma(alpha+beta)-lgamma(alpha)-lgamma(beta))*theta^(alpha-1)*(1-theta)^(beta-1)

      }

f0<-prior.theta(theta,alpha,beta)

# n = # of trials (IVDU subjects)

# u = # of "successes" (HIV positives)

# Three different values of u are used for illustration

n=100

u1=1

u2=5

u3=15

# post.theta is a function that constructs the posterior distribution

post.theta<-function(n,theta,u,alpha,beta){

      exp(lgamma(n+alpha+beta)-lgamma(u+alpha)-lgamma(n+beta-u))*theta^(u+alpha-1)*(1-theta)^(n+beta-u-1)

      }

f1<-post.theta(n,theta,u1,alpha,beta)

f2<-post.theta(n,theta,u2,alpha,beta)

f3<-post.theta(n,theta,u3,alpha,beta)

# Plot prior and three posteriors (using three different u's)

par(mfrow=c(2,2)) # sets up 2 by 2 matrix of plots

plot(theta,f0,type="l",xlab=expression(theta),ylab=expression(g(theta)),lty=1)

plot(theta,f1,type="l",xlab=expression(theta),ylab=expression(g(theta,1)),lty=1)

plot(theta,f2,type="l",xlab=expression(theta),ylab=expression(g(theta,5)),lty=1)

plot(theta,f3,type="l",xlab=expression(theta),ylab=expression(g(theta,15)),lty=1)

 

# Example 12.3.

# Poisson gamma Bayesian analysis

# Posterior construction

# Page 114-117

# theta is the Poisson mean

theta<-seq(0,20,0.01)

# Parameters in the gamma(alpha,beta) prior

alpha=2

beta=3

# prior.theta is a function that constructs the prior distribution

# lgamma is the internal log(gamma) function

prior.theta<-function(theta,alpha,beta){

      1/(exp(lgamma(alpha))*beta^alpha)*theta^(alpha-1)*exp(-theta/beta)

      }

f0<-prior.theta(theta,alpha,beta)

# n = sample size (# of rat mothers)

# u = sum (# of rat pups)

# Three different values of u are used for illustration

n=10

u1=32

u2=57

u3=90

# Posterior scale parameter

# Defining this just makes programming easier

beta.new<-1/(n+1/beta)

# post.p is a function that constructs the posterior distribution

post.theta<-function(n,theta,u,alpha,beta.new){

      1/(exp(lgamma(u+alpha))*beta.new^(u+alpha))*theta^(u+alpha-1)*exp(-theta/beta.new)

      }

f1<-post.theta(n,theta,u1,alpha,beta.new)

f2<-post.theta(n,theta,u2,alpha,beta.new)

f3<-post.theta(n,theta,u3,alpha,beta.new)

# Plot prior and three posteriors (using three different u's)

par(mfrow=c(2,2)) # sets up 2 by 2 matrix of plots

plot(theta,f0,type="l",xlab=expression(theta),ylab=expression(g(theta)),lty=1)

plot(theta,f1,type="l",xlab=expression(theta),ylab=expression(g(theta,32)),lty=1)

plot(theta,f2,type="l",xlab=expression(theta),ylab=expression(g(theta,57)),lty=1)

plot(theta,f3,type="l",xlab=expression(theta),ylab=expression(g(theta,90)),lty=1)

 

# Example 12.7.

# Beta binomial Bayesian analysis

# Posterior point estimation

# Page 124-126

# Parameters in the beta(alpha,beta) prior

alpha=1

beta=19

# n = # of trials (IVDU subjects)

# u = # of "successes" (HIV positives)

# Three different values of u are used for illustration

n=100

u1=1

u2=5

u3=15

## Posterior mean

p.mean1<-(u1+alpha)/(n+alpha+beta)

p.mean2<-(u2+alpha)/(n+alpha+beta)

p.mean3<-(u3+alpha)/(n+alpha+beta)

## Posterior mode

p.mode1<-(u1+alpha-1)/(n+alpha+beta-2)

p.mode2<-(u2+alpha-1)/(n+alpha+beta-2)

p.mode3<-(u3+alpha-1)/(n+alpha+beta-2)

## Posterior median

p.med1<-qbeta(0.5,u1+alpha,n+beta-u1)

p.med2<-qbeta(0.5,u2+alpha,n+beta-u2)

p.med3<-qbeta(0.5,u3+alpha,n+beta-u3)

results.1<-data.frame(p.mean1,p.mode1,p.med1,u1/n)

results.2<-data.frame(p.mean2,p.mode2,p.med2,u2/n)

results.3<-data.frame(p.mean3,p.mode3,p.med3,u3/n)

results.1

results.2

results.3

 

# Example 12.8.

# Beta binomial Bayesian analysis

# Credible intervals

# Page 129-131

# Parameters in the beta(alpha,beta) prior

alpha=1

beta=19

# n = # of trials (IVDU subjects)

# u = # of "successes" (HIV positives)

# Three different values of u are used for illustration

n=100

u1=1

u2=5

u3=15

# Lower endpoint of 95 percent ET credible interval

p.L.1<-qbeta(0.025,u1+alpha,n+beta-u1)

p.L.2<-qbeta(0.025,u2+alpha,n+beta-u2)

p.L.3<-qbeta(0.025,u3+alpha,n+beta-u3)

# Upper endpoint of 95 percent ET credible interval

p.U.1<-qbeta(0.975,u1+alpha,n+beta-u1)

p.U.2<-qbeta(0.975,u2+alpha,n+beta-u2)

p.U.3<-qbeta(0.975,u3+alpha,n+beta-u3)

# Lower endpoint of 95 percent Wald interval

wald.L.1<-(u1/n)-qnorm(0.975)*sqrt((u1/n)*(1-u1/n)/n)

wald.L.2<-(u2/n)-qnorm(0.975)*sqrt((u2/n)*(1-u2/n)/n)

wald.L.3<-(u3/n)-qnorm(0.975)*sqrt((u3/n)*(1-u3/n)/n)

# Upper endpoint of 95 percent ET interval

wald.U.1<-(u1/n)+qnorm(0.975)*sqrt((u1/n)*(1-u1/n)/n)

wald.U.2<-(u2/n)+qnorm(0.975)*sqrt((u2/n)*(1-u2/n)/n)

wald.U.3<-(u3/n)+qnorm(0.975)*sqrt((u3/n)*(1-u3/n)/n)

# Print out results

results.1<-data.frame(p.L.1,p.U.1,wald.L.1,wald.U.1)

results.2<-data.frame(p.L.2,p.U.2,wald.L.2,wald.U.2)

results.3<-data.frame(p.L.3,p.U.3,wald.L.3,wald.U.3)

results.1

results.2

results.3

 

# Example 12.9.

# Beta binomial Bayesian analysis

# Hypothesis testing

# Page 132

# Parameters in the beta(alpha,beta) prior

alpha=1

beta=19

# n = # of trials (IVDU subjects)

# u = # of "successes" (HIV positives)

# Three different values of u are used for illustration

n=100

u1=1

u2=5

u3=15

# pr(H_0 is true), computed with respect to g(theta|u)

prob.1<-pbeta(0.05,u1+alpha,n+beta-u1)

prob.2<-pbeta(0.05,u2+alpha,n+beta-u2)

prob.3<-pbeta(0.05,u3+alpha,n+beta-u3)