MCMCnegbinChange: Markov Chain Monte Carlo for Negative Binomial Regression Changepoint Model

Description

This function generates a sample from the posterior distribution of a Negative Binomial regression model with multiple changepoints. For the changepoints, the sampler uses the Markov Chain Monte Carlo method of Chib (1998). The user supplies data and priors, and a sample from the posterior distribution is returned as an mcmc object, which can be subsequently analyzed with functions provided in the coda package.

Usage

MCMCnegbinChange(
  formula,
  data = parent.frame(),
  m = 1,
  fixed.m = TRUE,
  b0 = 0,
  B0 = 1,
  a = NULL,
  b = NULL,
  e = 2,
  f = 2,
  g = 10,
  burnin = 1000,
  mcmc = 1000,
  thin = 1,
  verbose = 0,
  seed = NA,
  beta.start = NA,
  P.start = NA,
  rho.start = NA,
  rho.step,
  nu.start = NA,
  marginal.likelihood = c("none", "Chib95"),
  ...
)

Value

An mcmc object that contains the posterior sample. This object can be summarized by functions provided by the coda package.

Arguments

formula: Model formula.
data: Data frame.
m: The number of changepoints.
fixed.m: A logical indicator for whether or not the number of changepoints in the sampler should be exactly equal to m or if that is simply an upper bound. Setting fixed.m to FALSE is equivalent to assuming a weak-limit approximation to a Dirichlet process mixture.
b0: The prior mean of $\beta$. This can either be a scalar or a column vector with dimension equal to the number of betas. If this takes a scalar value, then that value will serve as the prior mean for all of the betas.
B0: The prior precision of $\beta$. This can either be a scalar or a square matrix with dimensions equal to the number of betas. If this takes a scalar value, then that value times an identity matrix serves as the prior precision of beta. Default value of 0 is equivalent to an improper uniform prior for beta.
a: $a$ is the shape1 beta prior for transition probabilities. By default, the expected duration is computed and corresponding a and b values are assigned. The expected duration is the sample period divided by the number of states.
b: $b$ is the shape2 beta prior for transition probabilities. By default, the expected duration is computed and corresponding a and b values are assigned. The expected duration is the sample period divided by the number of states.
e: The hyperprior for the distribution $\rho$ See details.
f: The hyperprior for the distribution $\rho$. See details.
g: The hyperprior for the distribution $\rho$. See details.
burnin: The number of burn-in iterations for the sampler.
mcmc: The number of Metropolis iterations for the sampler.
thin: The thinning interval used in the simulation. The number of mcmc iterations must be divisible by this value.
verbose: A switch which determines whether or not the progress of the sampler is printed to the screen. If verbose is greater than 0 the iteration number, the current beta vector, and the Metropolis acceptance rate are printed to the screen every verboseth iteration.
seed: The seed for the random number generator. If NA, the Mersenne Twister generator is used with default seed 12345; if an integer is passed it is used to seed the Mersenne twister. The user can also pass a list of length two to use the L'Ecuyer random number generator, which is suitable for parallel computation. The first element of the list is the L'Ecuyer seed, which is a vector of length six or NA (if NA a default seed of rep(12345,6) is used). The second element of list is a positive substream number. See the MCMCpack specification for more details.
beta.start: The starting value for the $\beta$ vector. This can either be a scalar or a column vector with dimension equal to the number of betas. If this takes a scalar value, then that value will serve as the starting value for all of the betas. The default value of NA will use the maximum likelihood estimate of $\beta$ as the starting value for all regimes.
P.start: The starting values for the transition matrix. A user should provide a square matrix with dimension equal to the number of states. By default, draws from the Beta(0.9, 0.1) are used to construct a proper transition matrix for each raw except the last raw.
rho.start: The starting value for the $\rho$ variable. This can either be a scalar or a column vector with dimension equal to the number of regimes. If the value is scalar, it will be used for all regimes. The default value is a vector of ones.
rho.step: Tuning parameter for the slice sampling approach to sampling $rho$. Determines the size of the step-out used to find the correct slice to draw from. Lower values are more accurate, but will take longer (up to a fixed searching limit). Default is 0.1.
nu.start: The starting values for the random effect, $\nu$. The default value is a vector of ones.
marginal.likelihood: How should the marginal likelihood be calculated? Options are: none in which case the marginal likelihood will not be calculated or Laplace in which case the Laplace approximation (see Kass and Raftery, 1995) is used.
...: further arguments to be passed.

Details

MCMCnegbinChangesimulates from the posterior distribution of a Negative Binomial regression model with multiple changepoints using the methods of Chib (1998) and Fruehwirth-Schnatter et al (2009). The details of the model are discussed in Blackwell (2017).

The model takes the following form:

$$y_t \sim \mathcal{P}oisson(\nu_t\mu_t)$$

$$\mu_t = x_t ' \beta_m,\;\; m = 1, \ldots, M$$

$$\nu_t \sim \mathcal{G}amma(\rho_m, \rho_m)$$

Where $M$ is the number of states and $\beta_m$ and $\rho_m$ are parameters when a state is $m$ at $t$.

We assume Gaussian distribution for prior of $\beta$:

$$\beta_m \sim \mathcal{N}(b_0,B_0^{-1}),\;\; m = 1, \ldots, M$$

And:

$$p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M$$ Where $M$ is the number of states.

The overdispersion parameters have a prior with the following form:

$$f(\rho_m|e,f,g) \propto \rho^{e-1}(\rho + g)^{-(e+f)}$$

The model is simulated via blocked Gibbs conditonal on the states. The $\beta$ being simulated via the auxiliary mixture sampling method of Fuerhwirth-Schanetter et al. (2009). The $\rho$ is updated via slice sampling. The $\nu_i$ are updated their (conjugate) full conditional, which is also Gamma. The states are updated as in Chib (1998)

References

Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', Journal of Statistical Software. 42(9): 1-21. tools:::Rd_expr_doi("10.18637/jss.v042.i09").

Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. Scythe Statistical Library 1.0. http://scythe.wustl.edu.s3-website-us-east-1.amazonaws.com/.

Sylvia Fruehwirth-Schnatter, Rudolf Fruehwirth, Leonhard Held, and Havard Rue. 2009. ``Improved auxiliary mixture sampling for hierarchical models of non-Gaussian data'', Statistics and Computing 19(4): 479-492. <doi:10.1007/s11222-008-9109-4>

Matthew Blackwell. 2017. ``Game Changers: Detecting Shifts in Overdispersed Count Data,'' Political Analysis 26(2), 230-239. <doi:10.1017/pan.2017.42>

Examples

Run this code


 if (FALSE) {
   n <- 150
   reg <- 3
   true.s <- gl(reg, n/reg, n)
   rho.true <- c(1.5, 0.5, 3)
   b0.true <- c(1, 3, 1)
   b1.true <- c(1, -2, 2)
   x1 <- runif(n, 0, 2)
   nu.true <- rgamma(n, rho.true[true.s], rho.true[true.s])
   mu <- nu.true * exp(b0.true[true.s] + x1 * b1.true[true.s])
   y <- rpois(n, mu)

   posterior <- MCMCnegbinChange(y ~ x1, m = 2, verbose = 1000,
                          marginal.likelihood = "Chib95",
                          e = 2, f = 2, g = 10,
                          b0 = rep(0, 2), B0 = (1/9) * diag(2),
                          rho.step = rep(0.75, times = 3),
                          seed = list(NA, 2))

   par(mfrow=c(attr(posterior, "m") + 1, 1), mai=c(0.4, 0.6, 0.3, 0.05))
   plotState(posterior, legend.control = c(1, 0.6))
   plotChangepoint(posterior, verbose = TRUE, ylab="Density",
  start=1, overlay=TRUE)


   open.ended <- MCMCnegbinChange(y ~ x1, m = 10, verbose = 1000,
                          fixed.m = FALSE, mcmc = 2000, burnin = 2000,
                          e = 2, f = 2, g = 10,
                          a = 100, b = 1,
                          b0 = rep(0, 2), B0 = (1/9) * diag(2),
                          rho.step = 0.75,
                          seed = list(NA, 2))

   plotState(open.ended, legend.control = c(1, 0.6))
   }

Run the code above in your browser using DataLab