empb_norm_negbinomial: empb_norm_negbinomial

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

Description

empb_norm_negbinomial

Usage

empb_norm_negbinomial(
  df,
  lambda = 1,
  MLEeta = 0.001,
  EMPBeta = 0.001,
  tol = 0.1,
  maxIter = 10000,
  jitter = FALSE
)

Arguments

df	data.frame object, containing at least columns 'x' containing non-negative integer values, and 'g' containing group labels.
lambda	non-negative numeric regularization parameter.
MLEeta	positive numeric dampening parameter for MLE fitting in group-level normal likelihood approximations (see 'mle_negbinomial' function).
EMPBeta	positive numeric dampening parameter for gradient descent algorithm, for this function.
tol	non-negative numeric tolerance parameter for exiting optimization algorithm.
maxIter	positive integer setting maximum number of iterations for optimization algorithm.
jitter	Boolean, to add small-magnitude noise to initial precision matrix estimate to ensure invertibility; default FALSE.

Value

list object containing empirical Bayes (EMPB) estimates of mu, Sigma hyperparameters, assuming df$x ~ nbinom(r_g, p_g), and (r_g, p_g) ~ MVN(mu, Sigma), where 'p_g' and 'r_g' denote group-level parameters.

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 empirical Bayes (EMPB) solutions for mu and Sigma:
muS_fit = empb_norm_negbinomial(df = df)

# Extract r_empb, p_empb; compare fitted values to r and (a / (a + b)):
cbind(c(r, a / (a + b)), c(exp(muS_fit$mu[1]), 1 / (1 + exp(-1 * muS_fit$mu[2]))))

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