posterior_beta_binomial: posterior_beta_binomial

In alexandercoulter/DHBayes: Bayesian Hierarchical Sampler for Discrete (Binomial, Negative-Binomial, Poisson) Data

Description

posterior_beta_binomial

Usage

posterior_beta_binomial(
  df,
  ab_prior = NULL,
  dist = c("post", "pred"),
  pred_n = NULL,
  Nsamp = 1000,
  ...
)

Arguments

df	data.frame object, containing at least columns named 'x' containing non-negative integer values (number of successes), 'n' containing non-negative integer values (number of trials), and 'g' containing group labels.
ab_prior	length-2, positive numeric vector specifying prior hyperparameter values for a, b; if NULL, values fit by empirical Bayes (EMPB).
dist	string specifying what type of samples to return: either 'post' (for samples from the posterior distribution), or 'pred' (for samples from the posterior-predictive distribution).
pred_n	positive integer, of length 1 or length equal to number of unique elements in df$g: 'size' of binomial samples if generating samples from posterior-predictive distribution.
Nsamp	positive integer, number of samples to generate per group.
...	optional parameters to be passed to control EMPB convergence, in the case 'ab_prior' is NULL; see 'empb_beta_binomial' or 'empb_beta_binomial_c'.

Value

matrix of samples from the posterior (or posterior-predictive) distribution, where (named) columns are for group IDs included in df$g, and rows are samples; assuming df$x ~ binom(p_g, df$n), and p_g ~ beta(a, b), where 'p_g' denotes a group-level parameter.

Examples

# Generate example data:
set.seed(31)
a = 3
b = 9

# Number of groups:
NG = 10

# Creating group IDs:
g = replicate(NG, paste(sample(LETTERS, 10), sep="", collapse=""))

# Generating 'true' p parameters:
p = rbeta(length(g), a, b)

# Number of experiments, i.e. rows in df:
numexps = 100

# Filling df with pseudo data; note the requisite columns 'n', 'x', and 'g':
df = data.frame('n' = numeric(0), 'x' = numeric(0), 'g' = character(0))
for(k in 1:numexps){
  gk = sample(g, 1)
  nk = 5 + rpois(1, 10)
  xk = rbinom(1, nk, p[g == gk])
  df = rbind(df, data.frame('n' = nk, 'x' = xk, 'g' = gk))
}

# Generating 1000 posterior distribution samples for each group:
posterior_values = posterior_beta_binomial(df = df)
dim(posterior_values)

# Create histogram of posterior distribution samples for first group (by alphabetic order):
hist(posterior_values[, 1], main = colnames(posterior_values)[1])

alexandercoulter/DHBayes documentation built on Dec. 19, 2021, 12:29 a.m.