pois.krige.bayes is a function for Bayesian geostatistical
inference in the Poisson spatial model.
It can be used for an analysis with fixed values of the parameters. However, the
function pois.krige may be preferable in this case.
The Bayesian algorithm is using a discretized version of the prior
distribution for the parameter \(\phi\). This means that the prior for \(\phi\) is always a proper prior.
For simulating from the posterior distribution of \(S\) given \(y\),
we use a Langevin-Hastings type algorithm.
This algorithm is a Metropolis-Hastings algorithm, where the
proposal distribution uses gradient information from the posterior. The algorithm is described below.
For shortness of presentation, we present only the MCMC algorithm for
the case where \(\beta\) follows a uniform prior and the link
function is the canonical log-link.
When \(\beta\) follows a uniform prior and the prior for \(\sigma^2\) is a scaled inverse-\(\chi^2\)
distribution,
the marginalised improper density of \(S\) is
$$
f_I(s) \propto |D^T V^{-1}D|^{-1/2}|V|^{-1/2}\{n_{\sigma}S^2_{\sigma}
+ s^T (V^{-1}-V^{-1}D(D^T V^{-1}D)^{-1}D^T V^{-1})s \}^{-(n-p+n_{\sigma})/2},
$$
where \(V\) is the correlation matrix of \(S\). The uniform prior for \(\sigma^2\) corresponds
to \(S^2_{\sigma}=0\) and
\(n_{\sigma}=-2\), and the reciprocal prior for \(\sigma^2\) corresponds to \(S^2_{\sigma}=0\)
and \(n_{\sigma}=0\).
We use the reparametrisation \(S = Q\Gamma\), where \(Q\) is the Cholesky factorisation of \(V\) so that \(V=QQ^T\).
Posterior simulations of \(S\) are obtained by transforming
MCMC simulations from the conditional distribution of \(\Gamma\) given
\(Y=y\).
The log posterior density of
\(\Gamma\) given \(Y=y\) is
$$
\log f(\gamma|y) = const(y) - \frac{1}{2} \gamma^T (I_n - V^{-1/2}D(D^T V^{-1}D)^{-1}D^T V^{-1/2})\gamma
+\sum_{i=1}^n y_i s_i - \exp(s_i),
$$
where
\((s_1,\ldots,s_n)^T= Q \gamma\).
For the truncated Langevin-Hastings update we use a truncated form of the gradient (truncating by \(H_i\)) of the log target density,
$$
\nabla(\gamma )^{trunc}= - (I_n - Q^{-1}D(D^T V^{-1}D)^{-1}D^T(Q^{-1})^T)\gamma + Q^T\left\{y_i-\exp(s_i)\wedge H_i \right\}_{i=1}^n .
$$
The proposal \(\gamma'\) follows a multivariate normal distribution with mean vector
\(\xi(\gamma)=\gamma+(h/2)\nabla(\gamma)^{trunc}\) and covariance matrix \(hI\),
where \(h>0\) is a user-specified ``proposal variance'' (called
S.scale; see mcmc.control).
When phi.prior is not "fixed", we update the parameter \(\phi\) by a random walk Metropolis step.
Here mcmc.input$phi.scale (see mcmc.control) is the
proposal variance, which needs to be sufficiently large so that the algorithm easily can move between the discrete values in prior$phi.discrete
(see prior.glm.control).
CONTROL FUNCTIONS
The function call includes auxiliary control functions which allows
the user to specify and/or change the specification of 1) model
components
(using model.glm.control), 2) prior
distributions (using prior.glm.control), 3) options for the
MCMC algorithm (using mcmc.control), and 4) options for the
output (using output.glm.control).
Default values are available in most of the cases.
The arguments for the control functions are described in their
respective help files.
In the prediction part of the function we want to predict
\(g_{\lambda}^{-1}(S^*)\) at locations of
interest, where \(g_{\lambda}^{-1}\) is the inverse
Box-Cox transformation.
For the prediction part of the algorithm, we use the median of the
predictive distribution as the predictor and 1/4 of the length of the 95 percent predictive interval as a measure of the prediction
uncertainty. Below we describe the procedure for calculating these quantities.
First we perform a Bayesian Gaussian prediction with the given priors on \(\beta\)
and \(\sigma^2\) on each of the simulated \(S\)-``datasets'' from the
posterior distribution (and in case \(\phi\) is not fixed, for each sampled \(\phi\) value).
This Gaussian prediction is done by calling an internal
function which
is an extension of krige.bayes
allowing for more than one ``data set''.
For calculating the probability below a threshold for the predictive distribution given the data, we first calculate this
probability
for each of the simulated \(S\)-``datasets''.
This is done using the fact that the predictive distribution
for each of the simulated \(S\)-``datasets'' is a multivariate \(t\)-distribution.
Afterwards the probability below a threshold is calculated by taking the empirical mean of these conditional probabilities.
Now the median and the 2.5 percent and 97.5 percent quantiles can be calculated by an iterative procedure, where first a guess
of the value is made, and second this guess is checked by calculating the probability of being less than this value. In case
the guess is too low, it is adjusted upwards, and vise verse.