MCMCdynamicEI | R Documentation |
MCMCdynamicEI is used to fit Quinn's dynamic ecological inference model for partially observed 2 x 2 contingency tables.
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,
...
)
r0 |
|
r1 |
|
c0 |
|
c1 |
|
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
|
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
|
W |
Weight (not precision) matrix structuring the
temporal dependence among elements of |
a0 |
|
b0 |
|
a1 |
|
b1 |
|
... |
further arguments to be passed |
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.
An mcmc object that contains the sample from the posterior distribution. This object can be summarized by functions provided by the coda package.
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. \Sexpr[results=rd]{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.
MCMChierEI
,
plot.mcmc
,summary.mcmc
## 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.