mcmcsae (version 0.8.1)

f_negbinomial: Specify a negative binomial sampling distribution

Description

This function can be used in the family argument of create_sampler or generate_data to specify a negative binomial sampling distribution.

Usage

f_negbinomial(
  link = "log",
  shape.vec = ~1,
  inv.shape.prior = pr_invchisq(df = 1, scale = 1),
  control = negbinomial_control()
)

Value

A family object.

Arguments

link

the name of a link function. Currently the only allowed link function for the negative binomial sampling distribution is "log".

shape.vec

optional formula specification of unequal shape values. The negative binomial (vector) shape parameter is then equal to this vector of shape values, multiplied by the scalar shape parameter, whose prior is specified through inv.shape.prior.

inv.shape.prior

Prior on the (scalar) reciprocal shape parameter, i.e. the overdispersion parameter. Supported prior distributions are pr_fixed with a default value of 1, pr_invchisq and pr_gig. The current default is pr_invchisq(df=1, scale=1).

control

a list with computational options. These options can be specified using function negbinomial_control.

Details

The negative binomial distribution with shape r and probability p has density $$p(y|r, p) = {\Gamma(y + r)\over y!\Gamma(r)}(1-p)^r p^y$$ with mean \(\mu = E(y|r,p) = {rp\over 1-p}\) and variance \(V(y|r,p) = \mu(1 + \mu/r)\). The second term of the variance can be interpreted as overdispersion with respect to a Poisson distribution, which would correspond to the limit \(r \rightarrow \infty\). So the reciprocal shape \(1/r\) is an overdispersion parameter, which typically is inferred. It is assigned a default prior, which may be changed through argument inv.shape.prior.

The only supported link function is "log". Strictly speaking the relation between mean \(\mu\) and linear predictor \(\eta\) is $$\log\mu = \log r + \log{p\over 1-p} = \log r + \eta$$ This way the likelihood function has the same form as that of logistic binomial regression, so that a Polya-Gamma data augmentation sampling algorithm can be employed. Note that the fact that the linear predictor \(\eta\) does not include \(\log r\) effectively changes the interpretation of its intercept.

References

N. Polson, J.G. Scott and J. Windle (2013). Bayesian Inference for Logistic Models Using Polya-Gamma Latent Variables. Journal of the American Statistical Association 108(504), 1339-1349.

M. Zhou and L. Carin (2015). Negative Binomial Process Count and Mixture Modeling. IEEE Transactions on Pattern Analysis and Machine Intelligence 37(2), 307-320.

Examples

if (FALSE) {
# first generate some data based on a negative binomial structured additive regression model
n <- 1200
nT <- 50
dat <- data.frame(
  x = runif(n),
  t = factor(sample(1:nT, n, replace=TRUE), levels=1:nT)
)
# generate randow walk over time t
v <- cumsum(rnorm(nT, sd=0.1))
v <- v - mean(v)  # subtract mean to separate random walk from intercept
lp <- with(dat, 1 - 0.5*x + v[t])
shape <- 2  # inverse negative binomial dispersion parameter
dat$y <- rnbinom(n, size=shape, mu = shape * exp(lp))

# now fit a negative binomial model to this data
sampler <- create_sampler(
  y ~ reg(~ x, name="beta") + gen(factor = ~ RW1(t), name="v"),
  data=dat, family = f_negbinomial()
)
sim <- MCMCsim(sampler, store.all=TRUE, plot.trace=c("v_sigma", "negbin_shape_"))
(summ <- summary(sim))
loo(sim)
bayesplot::mcmc_recover_intervals(as.array(sim$beta), c(1, -0.5))
bayesplot::mcmc_scatter(as.array(sim$beta))
bayesplot::mcmc_recover_hist(as.array(sim$negbin_shape_), shape)
bayesplot::mcmc_recover_hist(as.array(sim$v_sigma), 0.1)
plot.ts(v); lines(summ$v[, "Mean"], col=2)
pred <- predict(sim, type="response")
summpred <- summary(pred)
# NB prediction with type="response" does not include factor shape
plot(exp(lp), summpred[, "Mean"]); abline(0, 1)
plot(dat$y, shape * summpred[, "Mean"]); abline(0, 1)
mean(dat$y); shape * mean(summpred[, "Mean"])
# posterior predictive p-values at the observation level:
ppred <- predict(sim, ppcheck=TRUE)
hist(attr(ppred, "ppp"))
# NB the peak at 1 is due to zeros in the data:
dat$y[attr(ppred, "ppp") == 1]
# as Pr(y_rep >= y) = 1 by construction for these cases
}