# # Chapter 5 -- Bayesian Computation With R # Beta-Binomial Model for Overdispersion # # Remove all objects just to be safe # rm(list=ls(all=TRUE)) # library(LearnBayes) # # Set up log function for problem -- In this case the log of # a ratio of two Beta Distributions -- here the y are the data -- # in this case cancer deaths and n are the number people at risk. # # y ~ Beta-Binomial(eta, 1/K) # logf <- function(y, n, K, eta) { lbeta(K*eta + y, K*(1-eta)+n-y) - lbeta(K*eta, K*(1-eta)) } # # This function computes the log of the Posterior Density -- A vague # prior is multiplied by the Beta-Binomial producing a Beta-Binomial # form. The function is passed the parameters -- eta and K -- of # distribution along with the data -- y and n. # betabinexch0 <- function(theta, data) { eta <- theta[1] K <- theta[2] y <- data[,1] n <- data[,2] N <- length(y) val <- sum(logf(y, n, K, eta)) val <- val - 2*log(1+K) - log(eta) - log(1 - eta) return(val) } # # This function computes the log of the Posterior Density -- A vague # prior is multiplied by the Beta-Binomial producing a Beta-Binomial # form. The function is passed the parameters -- Theta_1 and Theta_2 which # are Theta_1 = logit(eta) and Theta_2 = log(K) -- of # distribution along with the data -- y and n. Note that you have # to compute the Jacobian of the transformation which ends up being # exp(Theta_1 + Theta_2)/(1 + exp(Theta_1))^2. When you take the log # of the prior and the log of the Jacobian you get some cancellations # which produce the next to the last line below. # betabinexch <- function(theta, data) { eta <- exp(theta[1])/(1 + exp(theta[1])) K <- exp(theta[2]) y <- data[,1] n <- data[,2] N <- length(y) val <- sum(logf(y, n, K, eta)) # # Note that this takes into account the coefficient of the posterior PLUS # the Jacobian!!! You have to take the logs of the coefficient with # the transformation and the logs of the Jacobian!! Then you get this: # val <- val + theta[2] - 2*log(1 + exp(theta[2])) return(val) } # # cancermortality is part of the LearnBayes Library # mycontour is part of the LearnBayes Library # data(cancermortality) windows() mycontour(betabinexch0,c(.0001, .003, 1, 20000),cancermortality,main="",xlab="",ylab="",col="blue",font=2) # Main title mtext("Contour Plot of parameters eta and K\nin the Beta-Binomial Model Problem",side=3,line=1.00,cex=1.2,font=2) # x-axis title mtext("eta",side=1,line=2.75,font=2,cex=1.2) # y-axis title mtext("K",side=2,line=2.75,font=2,cex=1.2) # windows() mycontour(betabinexch,c(-8, -4.5, 3, 16.5),cancermortality,main="",xlab="",ylab="",col="red",font=2) # Main title mtext("Contour Plot of Transformed parameters eta and K\nin the Beta-Binomial Model Problem",side=3,line=1.00,cex=1.2,font=2) # x-axis title mtext("logit eta",side=1,line=2.75,font=2,cex=1.2) # y-axis title mtext("log K",side=2,line=2.75,font=2,cex=1.2)