5.4 A Beta-Binomial Model for Overdispersion

5.4 A Beta-Binomial Model for Overdispersion
library('LearnBayes')
betabinexch0=function(theta, data)
{
	eta = theta[1]
	K = theta[2]
	y = data[, 1]
	n = data[, 2]
	N = length(y)
	logf = function(y, n, K, eta) 
		lbeta(K * eta + y, K * (1 -
		eta) + n - y) - lbeta(K * eta, K * (1 - eta))
	val = sum(logf(y, n, K, eta))
	val = val - 2 * log(1 + K) - log(eta) - log(1 - eta)
	return(val)
}

#log likelihood is 
#lbeta(K * eta + y, K * (1 -eta) + n - y) - lbeta(K * eta, K * (1 - eta))
#sum of log likelihood is "val"
#see math in page 90-91.
#note here that logf is a function, within a function
#lbeta returs the natural log of the beta function

data(cancermortality)
attach(cancermortality)
mycontour(betabinexch0,c(.0001,.003,1,20000),cancermortality, 
	xlab="eta",ylab="K")


#following guideline in section 5.3 - transformed parameter
#the log posterior density of transformed parameters is programmed
#using the betabinexch function

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)
	logf = function(y, n, K, eta) lbeta(K * eta + y, K * (1 -
	eta) + n - y) - lbeta(K * eta, K * (1 - eta))
	val = sum(logf(y, n, K, eta))
	val = val + theta[2] - 2 * log(1 + exp(theta[2]))
	return(val)
}

mycontour(betabinexch, c(-8,-4.5,3,16.5), cancermortality, 
	xlab="logit eta", ylab="log K")

#now we see that ALTHOUGH the density has an unusual shape, 
#the strong skewness has been reduced and the distribution 
#is more amenable to the computational methods to be described 
#in the following chapters.

#(testing changes)