bayesPhat: Baysian estimation to binomial proportion using a conjugate...
In tmcd82070/evoab: Evidence of Absence (EoA)

Description Usage Arguments Details Value Author(s) See Also Examples

This routine assumes a beta(a,b) prior for a binomial proportion and computes posterior point and inteval estimates for the proportion given an observed number of 'successes' and number of 'trials'.

1	bayesPhat(x, n, a = 0.5, b = 0.5, conf = 0.9, ci = "eq")

`x`	The number of 'successes' (an integer vector)
`n`	The number of 'trials' (an integer vector)
`a`	A vector of shape parameters to specify the beta prior
`b`	A vector of scale parameters to specify the beta prior
`conf`	A scalar specifying the desired confidence level of the interval estimate. This scalar should be between 0.5 and 1.
`ci`	A character string specifying the type of posterior confidence interval to compute. A value of 'eq' computes posterior intervals with equal area in both tails. A value of 'hdi' computes highest density interval estimates. See Details.

Computations are elementary because the beta(a,b) prior is conjugate for the binomial. Nearly every text on Bayesian estimation shows that given values for x, n, a, and b, the posterior distribution of p is,

beta(x+a, n-x+b).

Hence, the Bayes point estimator is,

phat = (x+a) / (n+a+b),

which is the mean of the posterior. Standard error of the posterior is,

se.phat=sqrt(a*b/((a+b)^2*(a+b+1))).

If a>1 and b>1, mode of the posterior for p is,

(a-1)/(a+b-2).

Confidence intervals can be computed two ways. The default, ci='eq', puts (1-conf)/2 probability in the lower tail and (1-conf)/2 in the upper tail. That is, the lower limit is qbeta((1-conf)/2,x+a,n-x+b) and the upper limit is qbeta(1-(1-conf)/2,x+a,n-x+b). If ci='hdi', this routine finds the interval [l,u] such that pbeta(u,x+a,n-x+b) - pbeta(l,x+a,n-x+b) and dbeta(u,x+a,n-x+b) = dbeta(l,x+a,n-x+b).

The default values for a and b imply use of the Jeffery's prior. The Jeffery's prior is proportional to the root of Fisher's information and is equal to beta(0.5,0.5). A flat prior is beta(1,1).

A data frame with number of rows equal to max(length(x), length(n), length(a), length(b) containing the Baysian point and interval estimates. The data frame has the following columns:

phat : the Bayes point estimates equal to the mean of the posterior distribution.
phat.mode : if a>1 and b>1, this column contains the mode of the posterior. If either a<=1 or b<=1, the posterior is either multi-modal or its mode is infinity. In this case, phat.mode = NA.
se.phat : the standard error of the posterior distribution.
ll.phat : the lower limit of a 100*(1-(1-conf)/2) posterior confidence interval for phat. See Details for interpretation under the allowed values of ci.
ul.phat : the upper limit of a 100*(1-(1-conf)/2) posterior confidence interval for phat. See Details for interpretation under the allowed values of ci.

Trent McDonald

agrestiCoullPhat

bayesPhat(0:5, 100, 0.5, 0.5)  # Jeffery's prior
bayesPhat(0:5, 100, 1, 1) # flat prior

# When prior data is available:
x.prior <- 5
n.prior <- 100
bayesPhat(0:5, 100, x.prior+0.5, n.prior-x.prior+0.5)

# Simulation: point est bias and ci coverage
trueP <- 0.01
n <- 1000
x <- rbinom( 1000, n, trueP)
baPhat <- bayesPhat( x, n )
muBA <- mean(baPhat$phat)
covBA <- mean(baPhat$ll.phat <= trueP & trueP <= baPhat$ul.phat)
baStats <- c(mean=muBA,
             relBias = abs(muBA-trueP)/trueP),
             coverage = covBA)