Package 'MCMCpack'

Description

This function creates an object of class BayesFactor from MCMCpack output.

Usage

BayesFactor(...)

is.BayesFactor(BF)

Arguments

Value

An object of class BayesFactor. A BayesFactor object has four attributes. They are: BF.mat an $M \times M$ matrix in which element $i,j$ contains the Bayes factor for model $i$ relative to model $j$; BF.log.mat an $M \times M$ matrix in which element $i,j$ contains the natural log of the Bayes factor for model $i$ relative to model $j$; BF.logmarglike an $M$ vector containing the log marginal likelihoods for models 1 through $M$; and BF.call an $M$ element list containing the calls used to fit models 1 through $M$.

See Also

MCMCregress

Examples

## Not run: 
data(birthwt)

model1 <- MCMCregress(bwt~age+lwt+as.factor(race) + smoke + ht,
                     data=birthwt, b0=c(2700, 0, 0, -500, -500,
                                        -500, -500),
                     B0=c(1e-6, .01, .01, 1.6e-5, 1.6e-5, 1.6e-5,
                          1.6e-5), c0=10, d0=4500000,
                     marginal.likelihood="Chib95", mcmc=10000)

model2 <- MCMCregress(bwt~age+lwt+as.factor(race) + smoke,
                     data=birthwt, b0=c(2700, 0, 0, -500, -500,
                                        -500),
                     B0=c(1e-6, .01, .01, 1.6e-5, 1.6e-5, 1.6e-5),
                     c0=10, d0=4500000,
                     marginal.likelihood="Chib95", mcmc=10000)

model3 <- MCMCregress(bwt~as.factor(race) + smoke + ht,
                     data=birthwt, b0=c(2700, -500, -500,
                                        -500, -500),
                     B0=c(1e-6, 1.6e-5, 1.6e-5, 1.6e-5,
                          1.6e-5), c0=10, d0=4500000,
                     marginal.likelihood="Chib95", mcmc=10000)

BF <- BayesFactor(model1, model2, model3)
print(BF)


## End(Not run)

Description

This function handles choice-specific covariates in multinomial choice models. See the example for an example of useage.

Usage

choicevar(var, varname, choicelevel)

Arguments

Value

The new variable used by the MCMCmnl() function.

See Also

MCMCmnl

Description

Density function and random generation from the Dirichlet distribution.

Usage

ddirichlet(x, alpha)

rdirichlet(n, alpha)

Arguments

Details

The Dirichlet distribution is the multidimensional generalization of the beta distribution.

Value

ddirichlet gives the density. rdirichlet returns a matrix with n rows, each containing a single Dirichlet random deviate.

Author(s)

Code is taken from Greg's Miscellaneous Functions (gregmisc). His code was based on code posted by Ben Bolker to R-News on 15 Dec 2000.

See Also

Beta

Examples

density <- ddirichlet(c(.1,.2,.7), c(1,1,1))
  draws <- rdirichlet(20, c(1,1,1) )

Description

dtomogplot is used to produce a tomography plot (see King, 1997) for a series of temporally ordered, partially observed 2 x 2 contingency tables.

Usage

dtomogplot(
  r0,
  r1,
  c0,
  c1,
  time.vec = NA,
  delay = 0,
  xlab = "fraction of r0 in c0 (p0)",
  ylab = "fraction of r1 in c0 (p1)",
  color.palette = heat.colors,
  bgcol = "black",
  ...
)

Arguments

Details

Consider the following partially observed 2 by 2 contingency table:

where $r_0$, $r_1$, $c_0$, $c_1$, and $N$ are non-negative integers that are observed. The interior cell entries are not observed. It is assumed that $Y_0|r_0 \sim \mathcal{B}inomial(r_0, p_0)$ and $Y_1|r_1 \sim \mathcal{B}inomial(r_1, p_1)$.

This function plots the bounds on the maximum likelihood estimates for (p0, p1) and color codes them by the elements of time.vec.

References

Gary King, 1997. A Solution to the Ecological Inference Problem. Princeton: Princeton University Press.

Jonathan C. Wakefield. 2004. “Ecological Inference for 2 x 2 Tables.” Journal of the Royal Statistical Society, Series A. 167(3): 385445.

Kevin Quinn. 2004. “Ecological Inference in the Presence of Temporal Dependence." In Ecological Inference: New Methodological Strategies. Gary King, Ori Rosen, and Martin A. Tanner (eds.). New York: Cambridge University Press.

See Also

MCMChierEI, MCMCdynamicEI,tomogplot

Examples

## Not run: 
## simulated data example 1
set.seed(3920)
n <- 100
r0 <- rpois(n, 2000)
r1 <- round(runif(n, 100, 4000))
p0.true <- pnorm(-1.5 + 1:n/(n/2))
p1.true <- pnorm(1.0 - 1:n/(n/4))
y0 <- rbinom(n, r0, p0.true)
y1 <- rbinom(n, r1, p1.true)
c0 <- y0 + y1
c1 <- (r0+r1) - c0

## plot data
dtomogplot(r0, r1, c0, c1, delay=0.1)

## simulated data example 2
set.seed(8722)
n <- 100
r0 <- rpois(n, 2000)
r1 <- round(runif(n, 100, 4000))
p0.true <- pnorm(-1.0 + sin(1:n/(n/4)))
p1.true <- pnorm(0.0 - 2*cos(1:n/(n/9)))
y0 <- rbinom(n, r0, p0.true)
y1 <- rbinom(n, r1, p1.true)
c0 <- y0 + y1
c1 <- (r0+r1) - c0

## plot data
dtomogplot(r0, r1, c0, c1, delay=0.1)

## End(Not run)

Description

Data on head-to-head outcomes from the 2016 UEFA European Football Championship.

Format

This dataframe contains all of the head-to-head results from Euro 2016. This includes results from both the group stage and the knock-out rounds.

dummy.rater

An artificial "dummy" rater equal to 1 for all matches. Included so that Euro2016 can be used directly with MCMCpack's models for pairwise comparisons.

team1

The home team

team2

The away team

winner

The winner of the match. NA if a draw.

Source

https://en.wikipedia.org/wiki/UEFA_Euro_2016

Description

This function generates a sample from the posterior distribution of a (sticky) HDP-HMM with a Negative Binomial outcome distribution (Fox et al, 2011). 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

HDPHMMnegbin(
  formula,
  data = parent.frame(),
  K = 10,
  b0 = 0,
  B0 = 1,
  a.theta = 50,
  b.theta = 5,
  a.alpha = 1,
  b.alpha = 0.1,
  a.gamma = 1,
  b.gamma = 0.1,
  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,
  gamma.start = 0.5,
  theta.start = 0.98,
  ak.start = 100,
  ...
)

Arguments

Details

HDPHMMnegbin simulates from the posterior distribution of a sticky HDP-HMM with a Negative Binomial outcome distribution, allowing for multiple, arbitrary changepoints in the model. The details of the model are discussed in Blackwell (2017). The implementation here is based on a weak-limit approximation, where there is a large, though finite number of regimes that can be switched between. Unlike other changepoint models in MCMCpack, the HDP-HMM approach allows for the state sequence to return to previous visited states.

The model takes the following form, where we show the fixed-limit version:

$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 an upper bound on the number of states and $\beta_m$ and $\rho_m$ are parameters when a state is $m$ at $t$.

The transition probabilities between states are assumed to follow a heirarchical Dirichlet process:

$\pi_m \sim \mathcal{D}irichlet(\alpha\delta_1, \ldots, \alpha\delta_j + \kappa, \ldots, \alpha\delta_M)$

$\delta \sim \mathcal{D}irichlet(\gamma/M, \ldots, \gamma/M)$

The $\kappa$ value here is the sticky parameter that encourages self-transitions. The sampler follows Fox et al (2011) and parameterizes these priors with $\alpha + \kappa$ and $\theta = \kappa/(\alpha + \kappa)$, with the latter representing the degree of self-transition bias. Gamma priors are assumed for $(\alpha + \kappa)$ and $\gamma$.

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

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

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 Fox et al (2011), supplemental materials.

Value

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

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. 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>

Emily B. Fox, Erik B. Sudderth, Michael I. Jordan, and Alan S. Willsky. 2011.. “A sticky HDP-HMM with application to speaker diarization.” The Annals of Applied Statistics, 5(2A), 1020-1056. <doi:10.1214/10-AOAS395>

See Also

MCMCnegbinChange, HDPHMMpoisson

Examples

## Not run: 
   n <- 150
   reg <- 3
   true.s <- gl(reg, n/reg, n)
   rho.true <- c(1.5, 0.5, 3)
   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(1 + x1 * b1.true[true.s])
   y <- rpois(n, mu)

   posterior <- HDPHMMnegbin(y ~ x1, K = 10, verbose = 1000,
                          e = 2, f = 2, g = 10,
                          a.theta = 100, b.theta = 1,
                          b0 = rep(0, 2), B0 = (1/9) * diag(2),
                          rho.step = rep(0.75, times = 10),
                          seed = list(NA, 2),
                          theta.start = 0.95, gamma.start = 10,
                          ak.start = 10)

   plotHDPChangepoint(posterior, ylab="Density", start=1)
   
## End(Not run)

Description

This function generates a sample from the posterior distribution of a (sticky) HDP-HMM with a Poisson outcome distribution (Fox et al, 2011). 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

HDPHMMpoisson(
  formula,
  data = parent.frame(),
  K = 10,
  b0 = 0,
  B0 = 1,
  a.alpha = 1,
  b.alpha = 0.1,
  a.gamma = 1,
  b.gamma = 0.1,
  a.theta = 50,
  b.theta = 5,
  burnin = 1000,
  mcmc = 1000,
  thin = 1,
  verbose = 0,
  seed = NA,
  beta.start = NA,
  P.start = NA,
  gamma.start = 0.5,
  theta.start = 0.98,
  ak.start = 100,
  ...
)

Arguments

Details

HDPHMMpoisson simulates from the posterior distribution of a sticky HDP-HMM with a Poisson outcome distribution, allowing for multiple, arbitrary changepoints in the model. The details of the model are discussed in Blackwell (2017). The implementation here is based on a weak-limit approximation, where there is a large, though finite number of regimes that can be switched between. Unlike other changepoint models in MCMCpack, the HDP-HMM approach allows for the state sequence to return to previous visited states.

The model takes the following form, where we show the fixed-limit version:

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

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

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

The transition probabilities between states are assumed to follow a heirarchical Dirichlet process:

$\pi_m \sim \mathcal{D}irichlet(\alpha\delta_1, \ldots, \alpha\delta_j + \kappa, \ldots, \alpha\delta_M)$

$\delta \sim \mathcal{D}irichlet(\gamma/M, \ldots, \gamma/M)$

The $\kappa$ value here is the sticky parameter that encourages self-transitions. The sampler follows Fox et al (2011) and parameterizes these priors with $\alpha + \kappa$ and $\theta = \kappa/(\alpha + \kappa)$, with the latter representing the degree of self-transition bias. Gamma priors are assumed for $(\alpha + \kappa)$ and $\gamma$.

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

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

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 states are updated as in Fox et al (2011), supplemental materials.

Value

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

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. 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>

Emily B. Fox, Erik B. Sudderth, Michael I. Jordan, and Alan S. Willsky. 2011.. “A sticky HDP-HMM with application to speaker diarization.” The Annals of Applied Statistics, 5(2A), 1020-1056. <doi:10.1214/10-AOAS395>

See Also

MCMCpoissonChange, HDPHMMnegbin

Examples

## Not run: 
   n <- 150
   reg <- 3
   true.s <- gl(reg, n/reg, n)
   b1.true <- c(1, -2, 2)
   x1 <- runif(n, 0, 2)
   mu <- exp(1 + x1 * b1.true[true.s])
   y <- rpois(n, mu)

   posterior <- HDPHMMpoisson(y ~ x1, K = 10, verbose = 1000,
                          a.theta = 100, b.theta = 1,
                          b0 = rep(0, 2), B0 = (1/9) * diag(2),
                          seed = list(NA, 2),
                          theta.start = 0.95, gamma.start = 10,
                          ak.start = 10)

   plotHDPChangepoint(posterior, ylab="Density", start=1)
   
## End(Not run)

Description

This function generates a sample from the posterior distribution of a Hidden Semi-Markov Model with a Heirarchical Dirichlet Process and a Negative Binomial outcome distribution (Johnson and Willsky, 2013). 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

HDPHSMMnegbin(
  formula,
  data = parent.frame(),
  K = 10,
  b0 = 0,
  B0 = 1,
  a.alpha = 1,
  b.alpha = 0.1,
  a.gamma = 1,
  b.gamma = 0.1,
  a.omega,
  b.omega,
  e = 2,
  f = 2,
  g = 10,
  r = 1,
  burnin = 1000,
  mcmc = 1000,
  thin = 1,
  verbose = 0,
  seed = NA,
  beta.start = NA,
  P.start = NA,
  rho.start = NA,
  rho.step,
  nu.start = NA,
  omega.start = NA,
  gamma.start = 0.5,
  alpha.start = 100,
  ...
)

Arguments

Details

HDPHSMMnegbin simulates from the posterior distribution of a HDP-HSMM with a Negative Binomial outcome distribution, allowing for multiple, arbitrary changepoints in the model. The details of the model are discussed in Johnson & Willsky (2013). The implementation here is based on a weak-limit approximation, where there is a large, though finite number of regimes that can be switched between. Unlike other changepoint models in MCMCpack, the HDP-HSMM approach allows for the state sequence to return to previous visited states.

The model takes the following form, where we show the fixed-limit version:

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

$\mu_t = x_t ' \beta_k,\;\; k = 1, \ldots, K$

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

Where $K$ is an upper bound on the number of states and $\beta_k$ and $\rho_k$ are parameters when a state is $k$ at $t$.

In the HDP-HSMM, there is a super-state sequence that, for a given observation, is drawn from the transition distribution and then a duration is drawn from a duration distribution to determin how long that state will stay active. After that duration, a new super-state is drawn from the transition distribution, where self-transitions are disallowed. The transition probabilities between states are assumed to follow a heirarchical Dirichlet process:

$\pi_k \sim \mathcal{D}irichlet(\alpha\delta_1, \ldots , \alpha\delta_K)$

$\delta \sim \mathcal{D}irichlet(\gamma/K, \ldots, \gamma/K)$

In the algorithm itself, these $\pi$ vectors are modified to remove self-transitions as discussed above. There is a unique duration distribution for each regime with the following parameters:

$D_k \sim \mathcal{N}egBin(r, \omega_k)$

$\omega_k \sim \mathcal{B}eta(a_{\omega,k}, b_{\omega, k})$

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

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

The overdispersion parameters have a prior with the following form:

$f(\rho_k|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_t$ are updated their (conjugate) full conditional, which is also Gamma. The states and their durations are drawn as in Johnson & Willsky (2013).

Value

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

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. 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>

Matthew J. Johnson and Alan S. Willsky. 2013. “Bayesian Nonparametric Hidden Semi-Markov Models.” Journal of Machine Learning Research, 14(Feb), 673-701.

See Also

MCMCnegbinChange, HDPHMMnegbin,

Examples

## Not run: 
   n <- 150
   reg <- 3
   true.s <- gl(reg, n/reg, n)
   rho.true <- c(1.5, 0.5, 3)
   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(1 + x1 * b1.true[true.s])
   y <- rpois(n, mu)

   posterior <- HDPHSMMnegbin(y ~ x1, K = 10, verbose = 1000,
                          e = 2, f = 2, g = 10,
                          b0 = 0, B0 = 1/9,
                          a.omega = 1, b.omega = 100, r = 1,
                          rho.step = rep(0.75, times = 10),
                          seed = list(NA, 2),
                          omega.start = 0.05, gamma.start = 10,
                          alpha.start = 5)

   plotHDPChangepoint(posterior, ylab="Density", start=1)
   
## End(Not run)

Description

HMMpanelFE generates a sample from the posterior distribution of the fixed-effects model with varying individual effects model discussed in Park (2011). The code works for both balanced and unbalanced panel data as long as there is no missing data in the middle of each group. This model uses a multivariate Normal prior for the fixed effects parameters and varying individual effects, an Inverse-Gamma prior on the residual error variance, and Beta prior for transition probabilities. 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

HMMpanelFE(
  subject.id,
  y,
  X,
  m,
  mcmc = 1000,
  burnin = 1000,
  thin = 1,
  verbose = 0,
  b0 = 0,
  B0 = 0.001,
  c0 = 0.001,
  d0 = 0.001,
  delta0 = 0,
  Delta0 = 0.001,
  a = NULL,
  b = NULL,
  seed = NA,
  ...
)

Arguments

Details

HMMpanelFE simulates from the fixed-effect hidden Markov pbject level:

$\varepsilon_{it} \sim \mathcal{N}(\alpha_{im}, \sigma^2_{im})$

We assume standard, semi-conjugate priors:

$\beta \sim \mathcal{N}(b_0,B_0^{-1})$

And:

$\sigma^{-2} \sim \mathcal{G}amma(c_0/2, d_0/2)$

And:

$\alpha \sim \mathcal{N}(delta_0,Delta_0^{-1})$

$\beta$, $\alpha$ and $\sigma^{-2}$ are assumed a priori independent.

And:

$p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M$

Where $M$ is the number of states.

OLS estimates are used for starting values.

Value

An mcmc object that contains the posterior sample. This object can be summarized by functions provided by the coda package. The object contains an attribute sigma storage matrix that contains time-varying residual variance, an attribute state storage matrix that contains posterior samples of hidden states, and an attribute delta storage matrix containing time-varying intercepts.

References

Jong Hee Park, 2012. “Unified Method for Dynamic and Cross-Sectional Heterogeneity: Introducing Hidden Markov Panel Models.” American Journal of Political Science.56: 1040-1054. <doi: 10.1111/j.1540-5907.2012.00590.x>

Siddhartha Chib. 1998. “Estimation and comparison of multiple change-point models.” Journal of Econometrics. 86: 221-241. <doi: 10.1016/S0304-4076(97)00115-2>

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. doi:10.18637/jss.v042.i09.

Examples

## Not run: 
  ## data generating
  set.seed(1974)
  N <- 30
  T <- 80
  NT <- N*T

  ## true parameter values
  true.beta <- c(1, 1)
  true.sigma <- 3
  x1 <- rnorm(NT)
  x2 <- runif(NT, 2, 4)

  ## group-specific breaks
  break.point = rep(T/2, N); break.sigma=c(rep(1, N));
  break.list <- rep(1, N)

  X <- as.matrix(cbind(x1, x2), NT, );
  y <- rep(NA, NT)
  id  <-  rep(1:N, each=NT/N)
  K <-  ncol(X);
  true.beta <- as.matrix(true.beta, K, 1)

  ## compute the break probability
  ruler <- c(1:T)
  W.mat <- matrix(NA, T, N)
  for (i in 1:N){
    W.mat[, i] <- pnorm((ruler-break.point[i])/break.sigma[i])
  }
  Weight <- as.vector(W.mat)

  ## draw time-varying individual effects and sample y
  j = 1
  true.sigma.alpha <- 30
  true.alpha1 <- true.alpha2 <- rep(NA, N)
  for (i in 1:N){
    Xi <- X[j:(j+T-1), ]
    true.mean <- Xi  %*% true.beta
    weight <- Weight[j:(j+T-1)]
    true.alpha1[i] <- rnorm(1, 0, true.sigma.alpha)
    true.alpha2[i] <- -1*true.alpha1[i]
    y[j:(j+T-1)] <- ((1-weight)*true.mean + (1-weight)*rnorm(T, 0, true.sigma) +
    		    (1-weight)*true.alpha1[i]) +
    		    (weight*true.mean + weight*rnorm(T, 0, true.sigma) + weight*true.alpha2[i])
    j <- j + T
  }

  ## extract the standardized residuals from the OLS with fixed-effects
  FEols <- lm(y ~ X + as.factor(id) -1 )
  resid.all <- rstandard(FEols)
  time.id <- rep(1:80, N)

  ## model fitting
  G <- 100
  BF <- testpanelSubjectBreak(subject.id=id, time.id=time.id,
         resid= resid.all, max.break=3, minimum = 10,
         mcmc=G, burnin = G, thin=1, verbose=G,
         b0=0, B0=1/100, c0=2, d0=2, Time = time.id)

  ## get the estimated break numbers
  estimated.breaks <- make.breaklist(BF, threshold=3)

  ## model fitting
  out <- HMMpanelFE(subject.id = id, y, X=X, m =  estimated.breaks,
             mcmc=G, burnin=G, thin=1, verbose=G,
             b0=0, B0=1/100, c0=2, d0=2, delta0=0, Delta0=1/100)

  ## print out the slope estimate
  ## true values are 1 and 1
  summary(out)

  ## compare them with the result from the constant fixed-effects
  summary(FEols)

## End(Not run)

Description

HMMpanelRE generates a sample from the posterior distribution of the hidden Markov random-effects model discussed in Park (2011). The code works for panel data with the same starting point. The sampling of panel parameters is based on Algorithm 2 of Chib and Carlin (1999). This model uses a multivariate Normal prior for the fixed effects parameters and varying individual effects, an Inverse-Wishart prior on the random-effects parameters, an Inverse-Gamma prior on the residual error variance, and Beta prior for transition probabilities. 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

HMMpanelRE(
  subject.id,
  time.id,
  y,
  X,
  W,
  m = 1,
  mcmc = 1000,
  burnin = 1000,
  thin = 1,
  verbose = 0,
  b0 = 0,
  B0 = 0.001,
  c0 = 0.001,
  d0 = 0.001,
  r0,
  R0,
  a = NULL,
  b = NULL,
  seed = NA,
  beta.start = NA,
  sigma2.start = NA,
  D.start = NA,
  P.start = NA,
  marginal.likelihood = c("none", "Chib95"),
  ...
)

Arguments

Details

HMMpanelRE simulates from the random-effect hidden Markov panel model introduced by Park (2011).

The model takes the following form:

$y_i = X_i \beta_m + W_i b_i + \varepsilon_i\;\; m = 1, \ldots, M$

Where each group $i$ have $k_i$ observations. Random-effects parameters are assumed to be time-varying at the system level:

$b_i \sim \mathcal{N}_q(0, D_m)$

$\varepsilon_i \sim \mathcal{N}(0, \sigma^2_m I_{k_i})$

And the errors: We assume standard, conjugate priors:

$\beta \sim \mathcal{N}_p(b0, B0)$

And:

$\sigma^{2} \sim \mathcal{IG}amma(c0/2, d0/2)$

And:

$D \sim \mathcal{IW}ishart(r0, R0)$

See Chib and Carlin (1999) for more details.

And:

$p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M$

Where $M$ is the number of states.

NOTE: We do not provide default parameters for the priors on the precision matrix for the random effects. When fitting one of these models, it is of utmost importance to choose a prior that reflects your prior beliefs about the random effects. Using the dwish and rwish functions might be useful in choosing these values.

Value

An mcmc object that contains the posterior sample. This object can be summarized by functions provided by the coda package. The object contains an attribute prob.state storage matrix that contains the probability of $state_i$ for each period, and the log-marginal likelihood of the model (logmarglike).

References

Jong Hee Park, 2012. “Unified Method for Dynamic and Cross-Sectional Heterogeneity: Introducing Hidden Markov Panel Models.” American Journal of Political Science.56: 1040-1054. <doi: 10.1111/j.1540-5907.2012.00590.x>

Siddhartha Chib. 1998. “Estimation and comparison of multiple change-point models.” Journal of Econometrics. 86: 221-241. <doi: 10.1016/S0304-4076(97)00115-2>

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. doi:10.18637/jss.v042.i09.

Examples

## Not run: 
  ## data generating
  set.seed(1977)
  Q <- 3
  true.beta1   <-  c(1, 1, 1) ; true.beta2   <-  c(-1, -1, -1)
  true.sigma2 <-  c(2, 5); true.D1 <- diag(.5, Q); true.D2 <- diag(2.5, Q)
  N=30; T=100;
  NT <- N*T
  x1 <- runif(NT, 1, 2)
  x2 <- runif(NT, 1, 2)
  X <- cbind(1, x1, x2);   W <- X;   y <- rep(NA, NT)

  ## true break numbers are one and at the center
  break.point = rep(T/2, N); break.sigma=c(rep(1, N));
  break.list <- rep(1, N)
  id  <-  rep(1:N, each=NT/N)
  K <-  ncol(X);
  ruler <- c(1:T)

  ## compute the weight for the break
  W.mat <- matrix(NA, T, N)
  for (i in 1:N){
    W.mat[, i] <- pnorm((ruler-break.point[i])/break.sigma[i])
  }
  Weight <- as.vector(W.mat)

  ## data generating by weighting two means and variances
  j = 1
  for (i in 1:N){
    Xi <- X[j:(j+T-1), ]
    Wi <- W[j:(j+T-1), ]
    true.V1 <- true.sigma2[1]*diag(T) + Wi%*%true.D1%*%t(Wi)
    true.V2 <- true.sigma2[2]*diag(T) + Wi%*%true.D2%*%t(Wi)
    true.mean1 <- Xi%*%true.beta1
    true.mean2 <- Xi%*%true.beta2
    weight <- Weight[j:(j+T-1)]
    y[j:(j+T-1)] <- (1-weight)*true.mean1 + (1-weight)*chol(true.V1)%*%rnorm(T) +
      weight*true.mean2 + weight*chol(true.V2)%*%rnorm(T)
    j <- j + T
  }
  ## model fitting
  subject.id <- c(rep(1:N, each=T))
  time.id <- c(rep(1:T, N))

  ## model fitting
  G <- 100
  b0  <- rep(0, K) ; B0  <- solve(diag(100, K))
  c0  <- 2; d0  <- 2
  r0  <- 5; R0  <- diag(c(1, 0.1, 0.1))
  subject.id <- c(rep(1:N, each=T))
  time.id <- c(rep(1:T, N))
  out1 <- HMMpanelRE(subject.id, time.id, y, X, W, m=1,
                     mcmc=G, burnin=G, thin=1, verbose=G,
                     b0=b0, B0=B0, c0=c0, d0=d0, r0=r0, R0=R0)

  ## latent state changes
  plotState(out1)

  ## print mcmc output
  summary(out1)




## End(Not run)

Description

Density function and random generation from the inverse gamma distribution.

Usage

dinvgamma(x, shape, scale = 1)

rinvgamma(n, shape, scale = 1)

Arguments

Details

An inverse gamma random variable with shape $a$ and scale $b$ has mean $\frac{b}{a-1}$ (assuming $a>1$) and variance $\frac{b^2}{(a-1)^2(a-2)}$ (assuming $a>2$).

Value

dinvgamma evaluates the density at x.

rinvgamma takes n draws from the inverse Gamma distribution. The parameterization is consistent with the Gamma Distribution in the stats package.

References

Andrew Gelman, John B. Carlin, Hal S. Stern, and Donald B. Rubin. 2004. Bayesian Data Analysis. 2nd Edition. Boca Raton: Chapman & Hall.

See Also

GammaDist

Examples

density <- dinvgamma(4.2, 1.1)
draws <- rinvgamma(10, 3.2)

Description

Density function and random generation from the Inverse Wishart distribution.

Usage

riwish(v, S)

diwish(W, v, S)

Arguments

Details

The mean of an inverse Wishart random variable with v degrees of freedom and scale matrix S is $(v-p-1)^{-1}S$.

Value

diwish evaluates the density at positive definite matrix W. riwish generates one random draw from the distribution.

Examples

density <- diwish(matrix(c(2,-.3,-.3,4),2,2), 3, matrix(c(1,.3,.3,1),2,2))
draw <- riwish(3, matrix(c(1,.3,.3,1),2,2))

Description

This function generates a vector of break numbers using the output of testpanelSubjectBreak. The function performs a pairwise comparison of models using Bayes Factors.

Usage

make.breaklist(BF, threshold = 3)

Arguments

Value

Vector fo break numbers.

References

Jong Hee Park, 2012. “Unified Method for Dynamic and Cross-Sectional Heterogeneity: Introducing Hidden Markov Panel Models.” American Journal of Political Science.56: 1040-1054. <doi: 10.1111/j.1540-5907.2012.00590.x>

Harold Jeffreys, 1961. The Theory of Probability. Oxford University Press.

See Also

testpanelSubjectBreak

Description

This function generates a sample from the posterior distribution of a binomial likelihood with a Beta prior.

Usage

MCbinomialbeta(y, n, alpha = 1, beta = 1, mc = 1000, ...)

Arguments

Details

MCbinomialbeta directly simulates from the posterior distribution. This model is designed primarily for instructional use. $\pi$ is the probability of success for each independent Bernoulli trial. We assume a conjugate Beta prior:

$\pi \sim \mathcal{B}eta(\alpha, \beta)$

$y$ is the number of successes in $n$ trials. By default, a uniform prior is used.

Value

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

See Also

plot.mcmc, summary.mcmc

Examples

## Not run: 
posterior <- MCbinomialbeta(3,12,mc=5000)
summary(posterior)
plot(posterior)
grid <- seq(0,1,0.01)
plot(grid, dbeta(grid, 1, 1), type="l", col="red", lwd=3, ylim=c(0,3.6),
  xlab="pi", ylab="density")
lines(density(posterior), col="blue", lwd=3)
legend(.75, 3.6, c("prior", "posterior"), lwd=3, col=c("red", "blue"))

## End(Not run)

Description

This function generates a sample from the posterior distribution of a binary model with multiple changepoints. The function 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

MCMCbinaryChange(
  data,
  m = 1,
  c0 = 1,
  d0 = 1,
  a = NULL,
  b = NULL,
  burnin = 10000,
  mcmc = 10000,
  thin = 1,
  verbose = 0,
  seed = NA,
  phi.start = NA,
  P.start = NA,
  marginal.likelihood = c("none", "Chib95"),
  ...
)

Arguments

Details

MCMCbinaryChange simulates from the posterior distribution of a binary model with multiple changepoints.

The model takes the following form:

$Y_t \sim \mathcal{B}ernoulli(\phi_i),\;\; i = 1, \ldots, k$

Where $k$ is the number of states.

We assume Beta priors for $\phi_{i}$ and for transition probabilities:

$\phi_i \sim \mathcal{B}eta(c_0, d_0)$

And:

$p_{mm} \sim \mathcal{B}eta{a}{b},\;\; m = 1, \ldots, k$

Where $M$ is the number of states.

Value

An mcmc object that contains the posterior sample. This object can be summarized by functions provided by the coda package. The object contains an attribute prob.state storage matrix that contains the probability of $state_i$ for each period, and the log-marginal likelihood of the model (logmarglike).

References

Jong Hee Park. 2011. “Changepoint Analysis of Binary and Ordinal Probit Models: An Application to Bank Rate Policy Under the Interwar Gold Standard." Political Analysis. 19: 188-204. <doi:10.1093/pan/mpr007>

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. doi:10.18637/jss.v042.i09.

Siddhartha Chib. 1995. “Marginal Likelihood from the Gibbs Output.” Journal of the American Statistical Association. 90: 1313-1321. <doi: 10.1080/01621459.1995.10476635>

See Also

MCMCpoissonChange,plotState, plotChangepoint

Examples

## Not run: 
    set.seed(19173)
    true.phi<- c(0.5, 0.8, 0.4)

    ## two breaks at c(80, 180)
    y1 <- rbinom(80, 1,  true.phi[1])
    y2 <- rbinom(100, 1, true.phi[2])
    y3 <- rbinom(120, 1, true.phi[3])
    y  <- as.ts(c(y1, y2, y3))

    model0 <- MCMCbinaryChange(y, m=0, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
    	      marginal.likelihood = "Chib95")
    model1 <- MCMCbinaryChange(y, m=1, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
    	      marginal.likelihood = "Chib95")
    model2 <- MCMCbinaryChange(y, m=2, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
    	      marginal.likelihood = "Chib95")
    model3 <- MCMCbinaryChange(y, m=3, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
    	      marginal.likelihood = "Chib95")
    model4 <- MCMCbinaryChange(y, m=4, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
    	      marginal.likelihood = "Chib95")
    model5 <- MCMCbinaryChange(y, m=5, c0=2, d0=2, mcmc=100, burnin=100, verbose=50,
    	      marginal.likelihood = "Chib95")

    print(BayesFactor(model0, model1, model2, model3, model4, model5))

    ## plot two plots in one screen
    par(mfrow=c(attr(model2, "m") + 1, 1), mai=c(0.4, 0.6, 0.3, 0.05))
    plotState(model2, legend.control = c(1, 0.6))
    plotChangepoint(model2, verbose = TRUE, ylab="Density", start=1, overlay=TRUE)

    
## End(Not run)

Description

MCMCdynamicEI is used to fit Quinn's dynamic ecological inference model for partially observed 2 x 2 contingency tables.

Usage

MCMCdynamicEI(
  r0,
  r1,
  c0,
  c1,
  burnin = 5000,
  mcmc = 50000,
  thin = 1,
  verbose = 0,
  seed = NA,
  W = 0,
  a0 = 0.825,
  b0 = 0.0105,
  a1 = 0.825,
  b1 = 0.0105,
  ...
)

Arguments

Details

Consider the following partially observed 2 by 2 contingency table for unit $t$ where $t=1,\ldots,ntables$:

Where $r_{0t}$, $r_{1t}$, $c_{0t}$, $c_{1t}$, and $N_t$ are non-negative integers that are observed. The interior cell entries are not observed. It is assumed that $Y_{0t}|r_{0t} \sim \mathcal{B}inomial(r_{0t}, p_{0t})$ and $Y_{1t}|r_{1t} \sim \mathcal{B}inomial(r_{1t}, p_{1t})$. Let $\theta_{0t} = log(p_{0t}/(1-p_{0t}))$, and $\theta_{1t} = log(p_{1t}/(1-p_{1t}))$.

The following prior distributions are assumed:

$p(\theta_0|\sigma^2_0) \propto \sigma_0^{-ntables} \exp \left(-\frac{1}{2\sigma^2_0} \theta'_{0} P \theta_{0}\right)$

and

$p(\theta_1|\sigma^2_1) \propto \sigma_1^{-ntables} \exp \left(-\frac{1}{2\sigma^2_1} \theta'_{1} P \theta_{1}\right)$

where $P_{ts}$ = $-W_{ts}$ for $t$ not equal to $s$ and $P_{tt}$ = $\sum_{s \ne t}W_{ts}$. The $\theta_{0t}$ is assumed to be a priori independent of $\theta_{1t}$ for all t. In addition, the following hyperpriors are assumed: $\sigma^2_0 \sim \mathcal{IG}(a_0/2, b_0/2)$, and $\sigma^2_1 \sim \mathcal{IG}(a_1/2, b_1/2)$.

Inference centers on $p_0$, $p_1$, $\sigma^2_0$, and $\sigma^2_1$. Univariate slice sampling (Neal, 2003) together with Gibbs sampling is used to sample from the posterior distribution.

Value

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

References

Kevin Quinn. 2004. “Ecological Inference in the Presence of Temporal Dependence." In Ecological Inference: New Methodological Strategies. Gary King, Ori Rosen, and Martin A. Tanner (eds.). New York: Cambridge University Press.

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. doi:10.18637/jss.v042.i09.

Radford Neal. 2003. “Slice Sampling" (with discussion). Annals of Statistics, 31: 705-767.

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/.

Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. “Output Analysis and Diagnostics for MCMC (CODA)”, R News. 6(1): 7-11. https://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf.

Jonathan C. Wakefield. 2004. “Ecological Inference for 2 x 2 Tables.” Journal of the Royal Statistical Society, Series A. 167(3): 385445.

See Also

MCMChierEI, plot.mcmc,summary.mcmc

Examples

## Not run: 
## simulated data example 1
set.seed(3920)
n <- 100
r0 <- rpois(n, 2000)
r1 <- round(runif(n, 100, 4000))
p0.true <- pnorm(-1.5 + 1:n/(n/2))
p1.true <- pnorm(1.0 - 1:n/(n/4))
y0 <- rbinom(n, r0, p0.true)
y1 <- rbinom(n, r1, p1.true)
c0 <- y0 + y1
c1 <- (r0+r1) - c0

## plot data
dtomogplot(r0, r1, c0, c1, delay=0.1)

## fit dynamic model
post1 <- MCMCdynamicEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100,
                    seed=list(NA, 1))

## fit exchangeable hierarchical model
post2 <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100,
                    seed=list(NA, 2))

p0meanDyn <- colMeans(post1)[1:n]
p1meanDyn <- colMeans(post1)[(n+1):(2*n)]
p0meanHier <- colMeans(post2)[1:n]
p1meanHier <- colMeans(post2)[(n+1):(2*n)]

## plot truth and posterior means
pairs(cbind(p0.true, p0meanDyn, p0meanHier, p1.true, p1meanDyn, p1meanHier))


## simulated data example 2
set.seed(8722)
n <- 100
r0 <- rpois(n, 2000)
r1 <- round(runif(n, 100, 4000))
p0.true <- pnorm(-1.0 + sin(1:n/(n/4)))
p1.true <- pnorm(0.0 - 2*cos(1:n/(n/9)))
y0 <- rbinom(n, r0, p0.true)
y1 <- rbinom(n, r1, p1.true)
c0 <- y0 + y1
c1 <- (r0+r1) - c0

## plot data
dtomogplot(r0, r1, c0, c1, delay=0.1)

## fit dynamic model
post1 <- MCMCdynamicEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100,
                    seed=list(NA, 1))

## fit exchangeable hierarchical model
post2 <- MCMChierEI(r0,r1,c0,c1, mcmc=40000, thin=5, verbose=100,
                    seed=list(NA, 2))

p0meanDyn <- colMeans(post1)[1:n]
p1meanDyn <- colMeans(post1)[(n+1):(2*n)]
p0meanHier <- colMeans(post2)[1:n]
p1meanHier <- colMeans(post2)[(n+1):(2*n)]

## plot truth and posterior means
pairs(cbind(p0.true, p0meanDyn, p0meanHier, p1.true, p1meanDyn, p1meanHier))
   
## End(Not run)