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