MCMCdynamicEI: Markov Chain Monte Carlo for Quinn's Dynamic Ecological Inference Model

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

Value

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

Arguments

r0: $(ntables \times 1)$ vector of row sums from row 0.
r1: $(ntables \times 1)$ vector of row sums from row 1.
c0: $(ntables \times 1)$ vector of column sums from column 0.
c1: $(ntables \times 1)$ vector of column sums from column 1.
burnin: The number of burn-in scans for the sampler.
mcmc: The number of mcmc scans to be saved.
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 then every verboseth iteration will be printed to the screen.
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.
W: Weight (not precision) matrix structuring the temporal dependence among elements of $\theta_{0}$ and $\theta_{1}$. The default value of 0 will construct a weight matrix that corresponds to random walk priors for $\theta_{0}$ and $\theta_{1}$. The default assumes that the tables are equally spaced throughout time and that the elements of $r0$, $r1$, $c0$, and $c1$ are temporally ordered.
a0: a0/2 is the shape parameter for the inverse-gamma prior on the $\sigma^2_0$ parameter.
b0: b0/2 is the scale parameter for the inverse-gamma prior on the $\sigma^2_0$ parameter.
a1: a1/2 is the shape parameter for the inverse-gamma prior on the $\sigma^2_1$ parameter.
b1: b1/2 is the scale parameter for the inverse-gamma prior on the $\sigma^2_1$ parameter.
...: further arguments to be passed

Details

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

	\| $Y=0$	\| $Y=1$	\|
---------	------------	------------	------------
$X=0$	\| $Y_{0t}$	\|	\| $r_{0t}$
---------	------------	------------	------------
$X=1$	\| $Y_{1t}$	\|	\| $r_{1t}$
---------	------------	------------	------------
	\| $c_{0t}$	\| $c_{1t}$	\| $N_t$

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.

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. tools:::Rd_expr_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.

Examples

Run this code


   if (FALSE) {
## 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))
   }

Run the code above in your browser using DataLab

	\| \(Y=0\)	\| \(Y=1\)	\|
---------	------------	------------	------------
\(X=0\)	\| \(Y_{0t}\)	\|	\| \(r_{0t}\)
---------	------------	------------	------------
\(X=1\)	\| \(Y_{1t}\)	\|	\| \(r_{1t}\)
---------	------------	------------	------------
	\| \(c_{0t}\)	\| \(c_{1t}\)	\| \(N_t\)