Laplace.sampling: Langevin-Hastings MCMC for conditional simulation

In barryrowlingson/PrevMap: Geostatistical Modelling of Spatially Referenced Prevalence Data

Description

This function simulates from the conditional distribution of a Gaussian random effect, given binomial observations y.

Usage

Laplace.sampling(mu, Sigma, y, units.m, control.mcmc, ID.coords = NULL,
  messages = TRUE, plot.correlogram = TRUE)

Arguments

mu	mean vector of the marginal distribution of the random effect.
Sigma	covariance matrix of the marginal distribution of the random effect.
y	vector of binomial observations.
units.m	vector of binomial denominators.
control.mcmc	output from control.mcmc.MCML.
ID.coords	vector of ID values for the unique set of spatial coordinates obtained from create.ID.coords. These must be provided if, for example, spatial random effects are defined at household level but some of the covariates are at individual level. Warning: the household coordinates must all be distinct otherwise see jitterDupCoords. Default is NULL.
messages	logical; if messages=TRUE then status messages are printed on the screen (or output device) while the function is running. Default is messages=TRUE.
plot.correlogram	logical; if plot.correlogram=TRUE the autocorrelation plot of the conditional simulations is displayed.

Details

Conditionally on the random effect S, the data y follow a binomial distribution with probability p and binomial denominators units.m. The logistic link function is used for the linear predictor, which assumes the form

\log(p/(1-p))=S.

The random effect S has a multivariate Gaussian distribution with mean mu and covariance matrix Sigma.

Laplace sampling. This function generates samples from the distribution of S given the data y. Specifically a Langevin-Hastings algorithm is used to update \tilde{S} = \tilde{Σ}^{-1/2}(S-\tilde{s}) where \tilde{Σ} and \tilde{s} are the inverse of the negative Hessian and the mode of the distribution of S given y, respectively. At each iteration a new value \tilde{s}_{prop} for \tilde{S} is proposed from a multivariate Gaussian distribution with mean

\tilde{s}_{curr}+(h/2)\nabla \log f(\tilde{S} | y),

where \tilde{s}_{curr} is the current value for \tilde{S}, h is a tuning parameter and \nabla \log f(\tilde{S} | y) is the the gradient of the log-density of the distribution of \tilde{S} given y. The tuning parameter h is updated according to the following adaptive scheme: the value of h at the i-th iteration, say h_{i}, is given by

h_{i} = h_{i-1}+c_{1}i^{-c_{2}}(α_{i}-0.547),

where c_{1} > 0 and 0 < c_{2} < 1 are pre-defined constants, and α_{i} is the acceptance rate at the i-th iteration (0.547 is the optimal acceptance rate for a multivariate standard Gaussian distribution). The starting value for h, and the values for c_{1} and c_{2} can be set through the function control.mcmc.MCML.

Random effects at household-level. When the data consist of two nested levels, such as households and individuals within households, the argument ID.coords must be used to define the household IDs for each individual. Let i and j denote the i-th household and the j-th person within that household; the logistic link function then assumes the form

\log(p_{ij}/(1-p_{ij}))=μ_{ij}+S_{i}

where the random effects S_{i} are now defined at household level and have mean zero.

Value

A list with the following components

samples: a matrix, each row of which corresponds to a sample from the predictive distribution.

h: vector of the values of the tuning parameter at each iteration of the Langevin-Hastings MCMC algorithm.

Author(s)

Emanuele Giorgi [e-mail]

Peter J. Diggle [e-mail]

See Also

control.mcmc.MCML, create.ID.coords.

Examples

set.seed(1234)
data(data_sim)
n.subset <- 50
data_subset <- data_sim[sample(1:nrow(data_sim),n.subset),]
mu <- rep(0,50)
Sigma <- varcov.spatial(coords=data_subset[,c("x1","x2")],
              cov.pars=c(1,0.15),kappa=2)$varcov
control.mcmc <- control.mcmc.MCML(n.sim=1000,burnin=0,thin=1,
                           h=1.65/(n.subset^2/3))
invisible(Laplace.sampling(mu=mu,Sigma=Sigma,
                                       y=data_subset$y,units.m=data_subset$units.m,
                                       control.mcmc=control.mcmc))

barryrowlingson/PrevMap documentation built on May 11, 2019, 6:24 p.m.