MCMCdynamicEI {MCMCpack} | 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 |
(ntables * 1) vector of row sums from row 0. |
r1 |
(ntables * 1) vector of row sums from row 1. |
c0 |
(ntables * 1) vector of column sums from column 0. |
c1 |
(ntables * 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 verbose th 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 theta0 and theta1. The default value of 0 will construct a weight matrix that corresponds to random walk priors for theta0 and theta1. 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 |
Consider the following partially observed 2 by 2 contingency table for
unit t where t=1,...,ntables:
| Y=0 | | Y=1 | | | |
- - - - - | - - - - - | - - - - - | - - - - - |
X=0 | | Y0[t] | | | | r0[t] |
- - - - - | - - - - - | - - - - - | - - - - - |
X=1 | | Y1[t] | | | | r1[t] |
- - - - - | - - - - - | - - - - - | - - - - - |
| c0[t] | | c1[t] | | N[t] |
Where r0-t, r1[t], c0[t], c1[t], and N[t] are non-negative integers that are observed. The interior cell entries are not observed. It is assumed that Y0[t]|r0[t] ~ Binomial(r0[t], p0[t]) and Y1[t]|r1[t] ~ Binomial(r1[t],p1[t]). Let theta0[t] = log(p0[t]/(1-p0[t])), and theta1[t] = log(p1[t]/(1-p1[t])).
The following prior distributions are assumed:
p(theta0|sigma^2_0) propto sigma^(-ntables)_0 exp(-1/(2*sigma^2_0) theta0' * P * theta0)
and
p(theta1|sigma^2_1) propto sigma^(-ntables)_1 exp(-1/(2*sigma^2_1) theta1' * P * theta1)
where P[t,s] = -W[t,s] for t not equal to s and P[t,t] = sum(W[t,]). The theta0[t] is assumed to be a priori independent of theta1[t] for all t. In addition, the following hyperpriors are assumed: σ^2_0 ~ InvGamma(a0/2, b0/2), and σ^2_1 ~ InvGamma(a1/2, b1/2).
Inference centers on p0, p1, 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.
Jonathan C. Wakefield. 2003. ``Ecological inference for 2x2 tables." Read before the Royal Statistical Society, on November 12th, 2003.
Radford Neal. 2003. ``Slice Sampling" (with discussion). Annals of Statistics, 31: 705-767.
Andrew D. Martin, Kevin M. Quinn, and Daniel Pemstein. 2004. Scythe Statistical Library 1.0. http://scythe.wustl.edu.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002. Output Analysis and Diagnostics for MCMC (CODA). http://www-fis.iarc.fr/coda/.
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)