MCMCpack (version 1.4-4)

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

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

Value

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

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}\)
--------- ------------ ------------ ------------

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. http://www.jstatsoft.org/v42/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.

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

Run this code
# NOT RUN {
   
# }
# 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))
   
# }
# NOT RUN {
# }

Run the code above in your browser using DataLab