posterior_beta_negbinomial: posterior_beta_negbinomial

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

Description

posterior_beta_negbinomial

Usage

posterior_beta_negbinomial(
  df,
  rab_prior = NULL,
  dist = c("post", "pred"),
  Nsamp = 1000,
  ...
)

Arguments

df	data.frame object, containing at least columns named 'x' containing non-negative integer values and 'g' containing group labels.
rab_prior	length-3, positive numeric vector specifying prior hyperparameter values for r, a, and 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).
Nsamp	positive integer, number of samples to generate per group.
...	optional parameters to be passed to control EMPB convergence, in the case 'rab_prior' is NULL; see 'empb_gamma_poisson'.

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 ~ nbinom(r, p_g), and p_g ~ beta(a, b), where 'p_g' denotes a group-level parameter.

Examples

# Generate example data:
set.seed(31)
r = 4
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 'x' and 'g':
df = data.frame('x' = numeric(0), 'g' = character(0))
for(k in 1:numexps){
  gk = sample(g, 1)
  xk = rnbinom(1, r, p[g == gk])
  df = rbind(df, data.frame('x' = xk, 'g' = gk))
}

# Generating 1000 posterior distribution samples for each group:
posterior_values = posterior_beta_negbinomial(df = df, method = 'gdescent')
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.