MCMChierEI {MCMCpack}R Documentation

Markov Chain Monte Carlo for Wakefield's Hierarchial Ecological Inference Model

Description

‘MCMChierEI’ is used to fit Wakefield's hierarchical ecological inference model for partially observed 2 x 2 contingency tables.

Usage

MCMChierEI(
  r0,
  r1,
  c0,
  c1,
  burnin = 5000,
  mcmc = 50000,
  thin = 1,
  verbose = 0,
  seed = NA,
  m0 = 0,
  M0 = 2.287656,
  m1 = 0,
  M1 = 2.287656,
  a0 = 0.825,
  b0 = 0.0105,
  a1 = 0.825,
  b1 = 0.0105,
  ...
)

Arguments

r0

(ntables×1)(ntables \times 1) vector of row sums from row 0.

r1

(ntables×1)(ntables \times 1) vector of row sums from row 1.

c0

(ntables×1)(ntables \times 1) vector of column sums from column 0.

c1

(ntables×1)(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.

m0

Prior mean of the μ0\mu_0 parameter.

M0

Prior variance of the μ0\mu_0 parameter.

m1

Prior mean of the μ1\mu_1 parameter.

M1

Prior variance of the μ1\mu_1 parameter.

a0

a0/2 is the shape parameter for the inverse-gamma prior on the σ02\sigma^2_0 parameter.

b0

b0/2 is the scale parameter for the inverse-gamma prior on the σ02\sigma^2_0 parameter.

a1

a1/2 is the shape parameter for the inverse-gamma prior on the σ12\sigma^2_1 parameter.

b1

b1/2 is the scale parameter for the inverse-gamma prior on the σ12\sigma^2_1 parameter.

...

further arguments to be passed

Details

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

| Y=0Y=0 | Y=1Y=1 |
--------- ------------ ------------ ------------
X=0X=0 | Y0tY_{0t} | | r0tr_{0t}
--------- ------------ ------------ ------------
X=1X=1 | Y1tY_{1t} | | r1tr_{1t}
--------- ------------ ------------ ------------
| c0tc_{0t} | c1tc_{1t} | NtN_t

Where r0tr_{0t}, r1tr_{1t}, c0tc_{0t}, c1tc_{1t}, and NtN_t are non-negative integers that are observed. The interior cell entries are not observed. It is assumed that Y0tr0tBinomial(r0t,p0t)Y_{0t}|r_{0t} \sim \mathcal{B}inomial(r_{0t}, p_{0t}) and Y1tr1tBinomial(r1t,p1t)Y_{1t}|r_{1t} \sim \mathcal{B}inomial(r_{1t}, p_{1t}). Let θ0t=log(p0t/(1p0t))\theta_{0t} = log(p_{0t}/(1-p_{0t})), and θ1t=log(p1t/(1p1t))\theta_{1t} = log(p_{1t}/(1-p_{1t})).

The following prior distributions are assumed: θ0tN(μ0,σ02)\theta_{0t} \sim \mathcal{N}(\mu_0, \sigma^2_0), θ1tN(μ1,σ12)\theta_{1t} \sim \mathcal{N}(\mu_1, \sigma^2_1). θ0t\theta_{0t} is assumed to be a priori independent of θ1t\theta_{1t} for all t. In addition, we assume the following hyperpriors: μ0N(m0,M0)\mu_0 \sim \mathcal{N}(m_0, M_0), μ1N(m1,M1)\mu_1 \sim \mathcal{N}(m_1, M_1), σ02IG(a0/2,b0/2)\sigma^2_0 \sim \mathcal{IG}(a_0/2, b_0/2), and σ12IG(a1/2,b1/2)\sigma^2_1 \sim \mathcal{IG}(a_1/2, b_1/2).

The default priors have been chosen to make the implied prior distribution for p0p_{0} and p1p_{1} approximately uniform on (0,1).

Inference centers on p0p_0, p1p_1, μ0\mu_0, μ1\mu_1, σ02\sigma^2_0, and σ12\sigma^2_1. Univariate slice sampling (Neal, 2003) along with Gibbs sampling is used to sample from the posterior distribution.

See Section 5.4 of Wakefield (2003) for discussion of the priors used here. MCMChierEI departs from the Wakefield model in that the mu0 and mu1 are here assumed to be drawn from independent normal distributions whereas Wakefield assumes they are drawn from logistic distributions.

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

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

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

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

See Also

MCMCdynamicEI, plot.mcmc,summary.mcmc

Examples


   ## Not run: 
## simulated data example
set.seed(3920)
n <- 100
r0 <- round(runif(n, 400, 1500))
r1 <- round(runif(n, 100, 4000))
p0.true <- pnorm(rnorm(n, m=0.5, s=0.25))
p1.true <- pnorm(rnorm(n, m=0.0, s=0.10))
y0 <- rbinom(n, r0, p0.true)
y1 <- rbinom(n, r1, p1.true)
c0 <- y0 + y1
c1 <- (r0+r1) - c0

## plot data
tomogplot(r0, r1, c0, c1)

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

p0meanHier <- colMeans(post)[1:n]
p1meanHier <- colMeans(post)[(n+1):(2*n)]

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


[Package MCMCpack version 1.7-0 Index]