1% Generated by roxygen2: do not edit by hand
2% Please edit documentation in R/MCMCpoissonChange.R
3\name{MCMCpoissonChange}
4\alias{MCMCpoissonChange}
5\title{Markov Chain Monte Carlo for a Poisson Regression Changepoint Model}
6\usage{
7MCMCpoissonChange(
8  formula,
9  data = parent.frame(),
10  m = 1,
11  b0 = 0,
12  B0 = 1,
13  a = NULL,
14  b = NULL,
15  c0 = NA,
16  d0 = NA,
17  lambda.mu = NA,
18  lambda.var = NA,
19  burnin = 1000,
20  mcmc = 1000,
21  thin = 1,
22  verbose = 0,
23  seed = NA,
24  beta.start = NA,
25  P.start = NA,
26  marginal.likelihood = c("none", "Chib95"),
27  ...
28)
29}
30\arguments{
31\item{formula}{Model formula.}
32
33\item{data}{Data frame.}
34
35\item{m}{The number of changepoints.}
36
37\item{b0}{The prior mean of \eqn{\beta}.  This can either be a scalar
38or a column vector with dimension equal to the number of betas. If this
39takes a scalar value, then that value will serve as the prior mean for all
40of the betas.}
41
42\item{B0}{The prior precision of \eqn{\beta}.  This can either be a
43scalar or a square matrix with dimensions equal to the number of betas.  If
44this takes a scalar value, then that value times an identity matrix serves
45as the prior precision of beta. Default value of 0 is equivalent to an
46improper uniform prior for beta.}
47
48\item{a}{\eqn{a} is the shape1 beta prior for transition probabilities.
49By default, the expected duration is computed and corresponding a and b
50values are assigned. The expected duration is the sample period divided by
51the number of states.}
52
53\item{b}{\eqn{b} is the shape2 beta prior for transition probabilities.
54By default, the expected duration is computed and corresponding a and b
55values are assigned. The expected duration is the sample period divided by
56the number of states.}
57
58\item{c0}{\eqn{c_0} is the shape parameter for Gamma prior on
59\eqn{\lambda} (the mean). When there is no covariate, this should be
60provided by users. No default value is provided.}
61
62\item{d0}{\eqn{d_0} is the scale parameter for Gamma prior on
63\eqn{\lambda} (the mean). When there is no covariate, this should be
64provided by users. No default value is provided.}
65
66\item{lambda.mu}{The mean of the Gamma prior on \eqn{\lambda}.
67\eqn{sigma.mu} and \eqn{sigma.var} allow users to
68choose the Gamma prior by choosing its mean and variance.}
69
70\item{lambda.var}{The variacne of the Gamma prior on \eqn{\lambda}.
71\eqn{sigma.mu} and \eqn{sigma.var} allow users to
72choose the Gamma prior by choosing its mean and variance.}
73
74\item{burnin}{The number of burn-in iterations for the sampler.}
75
76\item{mcmc}{The number of MCMC iterations after burn-in.}
77
78\item{thin}{The thinning interval used in the simulation.  The number of
79MCMC iterations must be divisible by this value.}
80
81\item{verbose}{A switch which determines whether or not the progress of the
82sampler is printed to the screen.  If \code{verbose} is greater than 0, the
83iteration number and the posterior density samples are printed to the screen
84every \code{verbose}th iteration.}
85
86\item{seed}{The seed for the random number generator.  If NA, current R
87system seed is used.}
88
89\item{beta.start}{The starting values for the beta vector. This can either
90be a scalar or a column vector with dimension equal to the number of betas.
91The default value of NA will use draws from the Uniform distribution with
92the same boundary with the data as the starting value. If this is a scalar,
93that value will serve as the starting value mean for all of the betas. When
94there is no covariate, the log value of means should be used.}
95
96\item{P.start}{The starting values for the transition matrix. A user should
97provide a square matrix with dimension equal to the number of states. By
98default, draws from the \code{Beta(0.9, 0.1)} are used to construct a proper
99transition matrix for each raw except the last raw.}
100
101\item{marginal.likelihood}{How should the marginal likelihood be calculated?
102Options are: \code{none} in which case the marginal likelihood will not be
103calculated, and \code{Chib95} in which case the method of Chib (1995) is
104used.}
105
106\item{...}{further arguments to be passed}
107}
108\value{
109An mcmc object that contains the posterior sample. This object can
110be summarized by functions provided by the coda package.  The object
111contains an attribute \code{prob.state} storage matrix that contains the
112probability of \eqn{state_i} for each period, and the log-marginal
113likelihood of the model (\code{logmarglike}).
114}
115\description{
116This function generates a sample from the posterior distribution of a
117Poisson regression model with multiple changepoints. The function uses the
118Markov chain Monte Carlo method of Chib (1998).  The user supplies data and
119priors, and a sample from the posterior distribution is returned as an mcmc
120object, which can be subsequently analyzed with functions provided in the
121coda package.
122}
123\details{
124\code{MCMCpoissonChange} simulates from the posterior distribution of a
125Poisson regression model with multiple changepoints using the methods of
126Chib (1998) and Fruhwirth-Schnatter and Wagner (2006).  The details of the
127model are discussed in Park (2010).
128
129The model takes the following form:
130
131\deqn{y_t \sim \mathcal{P}oisson(\mu_t)}
132
133\deqn{\mu_t = x_t ' \beta_m,\;\; m = 1, \ldots, M}
134
135Where
136\eqn{M} is the number of states and \eqn{\beta_m} is paramters
137when a state is \eqn{m} at \eqn{t}.
138
139We assume Gaussian distribution for prior of \eqn{\beta}:
140
141\deqn{\beta_m \sim \mathcal{N}(b_0,B_0^{-1}),\;\; m = 1, \ldots, M}
142
143And:
144
145\deqn{p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M} Where \eqn{M} is the number of states.
146}
147\examples{
148
149    \dontrun{
150    set.seed(11119)
151    n <- 150
152    x1 <- runif(n, 0, 0.5)
153    true.beta1 <- c(1,  1)
154    true.beta2 <- c(1,  -2)
155    true.beta3 <- c(1,  2)
156
157    ## set true two breaks at (50, 100)
158    true.s <- rep(1:3, each=n/3)
159    mu1 <- exp(1 + x1[true.s==1]*1)
160    mu2 <- exp(1 + x1[true.s==2]*-2)
161    mu3 <- exp(1 + x1[true.s==3]*2)
162
163    y <- as.ts(c(rpois(n/3, mu1), rpois(n/3, mu2), rpois(n/3, mu3)))
164    formula = y ~ x1
165
166    ## fit multiple models with a varying number of breaks
167    model0 <-  MCMCpoissonChange(formula, m=0,
168            mcmc = 1000, burnin = 1000, verbose = 500,
169            b0 = rep(0, 2), B0 = 1/5*diag(2), marginal.likelihood = "Chib95")
170    model1 <-  MCMCpoissonChange(formula, m=1,
171            mcmc = 1000, burnin = 1000, verbose = 500,
172            b0 = rep(0, 2), B0 = 1/5*diag(2), marginal.likelihood = "Chib95")
173    model2 <-  MCMCpoissonChange(formula, m=2,
174            mcmc = 1000, burnin = 1000, verbose = 500,
175            b0 = rep(0, 2), B0 = 1/5*diag(2), marginal.likelihood = "Chib95")
176    model3 <-  MCMCpoissonChange(formula, m=3,
177            mcmc = 1000, burnin = 1000, verbose = 500,
178            b0 = rep(0, 2), B0 = 1/5*diag(2), marginal.likelihood = "Chib95")
179    model4 <-  MCMCpoissonChange(formula, m=4,
180            mcmc = 1000, burnin = 1000, verbose = 500,
181            b0 = rep(0, 2), B0 = 1/5*diag(2), marginal.likelihood = "Chib95")
182    model5 <-  MCMCpoissonChange(formula, m=5,
183            mcmc = 1000, burnin = 1000, verbose = 500,
184            b0 = rep(0, 2), B0 = 1/5*diag(2), marginal.likelihood = "Chib95")
185
186    ## find the most reasonable one
187    print(BayesFactor(model0, model1, model2, model3, model4, model5))
188
189    ## draw plots using the "right" model
190    par(mfrow=c(attr(model2, "m") + 1, 1), mai=c(0.4, 0.6, 0.3, 0.05))
191    plotState(model2, legend.control = c(1, 0.6))
192    plotChangepoint(model2, verbose = TRUE, ylab="Density", start=1, overlay=TRUE)
193
194    ## No covariate case
195    model2.1 <- MCMCpoissonChange(y ~ 1, m = 2, c0 = 2, d0 = 1,
196             mcmc = 1000, burnin = 1000, verbose = 500,
197             marginal.likelihood = "Chib95")
198    print(BayesFactor(model2, model2.1))
199    }
200
201}
202\references{
203Jong Hee Park. 2010. ``Structural Change in the U.S. Presidents'
204Use of Force Abroad.'' \emph{American Journal of Political Science} 54:
205766-782. <doi:10.1111/j.1540-5907.2010.00459.x>
206
207Sylvia Fruhwirth-Schnatter and Helga Wagner 2006. ``Auxiliary Mixture
208Sampling for Parameter-driven Models of Time Series of Counts with
209Applications to State Space Modelling.'' \emph{Biometrika}. 93:827--841.
210
211Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point
212models.'' \emph{Journal of Econometrics}. 86: 221-241.
213<doi: 10.1016/S0304-4076(97)00115-2>
214
215Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.  ``MCMCpack:
216Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}.
21742(9): 1-21.  \doi{10.18637/jss.v042.i09}.
218
219Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.''
220\emph{Journal of the American Statistical Association}. 90: 1313-1321.
221<doi: 10.1080/01621459.1995.10476635>
222}
223\seealso{
224\code{\link{MCMCbinaryChange}}, \code{\link{plotState}},
225\code{\link{plotChangepoint}}
226}
227\keyword{models}
228