5.6 The example (p.95).txt

#5.6 The example (p.95)
library(LearnBayes)

data(cancermortality)
attach(cancermortality)

fit=laplace(betabinexch, c(-7,6), cancermortality)
fit

#laplace's output has four components
#mode gives the posterior mode (theta hat)
#var gives the associated variance-covariance matrix V
#int is the approximation to the logarithm of the prior predictive density
#converge indicates if agorithm converged

#use mycontour AND lbinorm (log bivariate normal density) 
#to display contours of the approximate bivariate normal density

npar=list(m=fit$mode, v=fit$var)
mycontour(lbinorm, c(-8,-4.5,3,16.5), npar, xlab="logit eta", ylab="log K")

#construct 90% probability intervals 
#use the diagonal elements of the variance-covariance matrix

se=sqrt(diag(fit$var))
fit$mode-1.645*se

fit$mode+1.645*se

#90% interval estimate for logit(η) is (−7.28,−6.36)
#90% interval estimate for logK is (5.67, 9.49)

########

#another run

#Approximations based on posterior modes
#5.6 the example (p.95-97)
#Based on our contour plot, we start the Nelder-Mead method with
#initial guess(logit(eta), logK) = (-7,6)

library(LearnBayes)
fit=laplace(betabinexch, c(-7,6), cancermortality)
fit

#Results:
#$mode is posterior mode
#$var is variance-covariance matrix
#we use these to make a "new" contour plot

npar=list(m=fit$mode, v=fit$var)
mycontour(lbinorm, c(-8, -4.5, 3, 16.5), npar, xlab="logit eta", ylab = "log K")

#se = sqrt of diagonal elements of "var"
#we use these to construct the 90% probability intervals for the parameters

se=sqrt(diag(fit$var))
CI = c(fit$mode-1.645*se, fit$mode+1.645*se)
CI

#first 2 values are lower bounds
#last 2 values are upper bounds