If the installation directories for the packages are not the default locations for R packages, type:
where "PATH_TO_geoR" and "PATH_TO_geoRglm" are the paths to the directories where geoR and geoRglm are installed, respectively. If the packages are correctly loaded the following messages will be displayed:
----------------------------------------------------
geoR - a package for geostatistical analysis in R
geoR version 1.3-10 (2003-03-10) is now loaded
----------------------------------------------------
-----------------------------------------------------------
geoRglm - a package for generalised linear spatial models
geoRglm version 0.6-2 (2003-03-24) is now loaded
-----------------------------------------------------------
Typically, data are stored as an object (a list) of class "geodata" (see the geoR introductory session for more details on this). For the data sets considered here, the object will sometimes include a vector units.m consisting of observation times (for the Poisson distribution) or numbers N in bi(N,p) (for the binomial distribution).
We use the data sets p50 and rongelap included in the geoRglm distribution for the examples presented in this document. These data sets can be loaded by typing:
Helpfiles are available for geoRglm. For getting help on the function pois.krige, just type:
Here we describe conditional simulation using MCMC and spatial prediction in the Poisson-log normal model, when covariance parameters are fixed. Full Bayesian methods are also implemented and will be presented in Section 3.
The nugget effect parameter (microscale variation) in the underlying Gaussian field can be set to a fixed value.
The same applies for the smoothness and anisotropy parameters. Options for taking covariates (trends) into
account are also included.
Conditional simulation and prediction with fixed covariance parameters in the Poisson-log normal model can be performed with options for either fixed beta (OK) or flat prior on beta (SK). The function uses a Langevin-Hastings MCMC algorithm for simulating from the conditional distribution.
An example where all parameters are fixed is shown below (for illustration purposes, some parameter values are just taken).
First we need to tune the algorithm by scaling the proposal variance so that acceptance rate is approximately 60 percent (optimal acceptance rate for Langevin-Hastings algorithm). This is done by trial and error.
After a few tryouts we decide to use S.scale = 0.5. We also need to study how well the chain is mixing.
Here the functions in the coda package would be useful for assessing the convergence and inspect the mixing of the MCMC algorithm. For a small demonstration of a few CODA functions used on the output above, see here.
To reduce the autocorrelation of the samples we decide to subsample every 10 iterations (default); when working with real data sets we may need to make a more enxtensive subsampling, say, storing only every 1000 iterations.
Now we make (minimal mean square error) prediction of the intensity at the two locations (0.5, 0.5) and (1, 0.4).
The output is a list including the predicted values (test2$predict), the prediction variances
(test2$krige.var) and the estimated Monte Carlo standard errors on the predicted values (test2$mcmc.error).
Please consider printing out the predicted values and the associated Monte Carlo standard errors:
By specifying sim.predict = TRUE, simulations are drawn from the predictive intensity at the two prediction locations
(test2$simulations).
These simulations are plotted below.
The way to specify that beta should follow a uniform prior would be:
Bayesian analysis for the Poisson-log normal model and the binomial-logit model is implemented by the functions
pois.krige.bayes and binom.krige.bayes, respectively.
Model parameters can be treated as fixed or random.
As an example consider first a model without nugget and including uncertainty in the beta and sigmasq parameters (mean and variance of the random effects S, respectively).
A Bayesian analysis is made by typing commands like:
Now chose S.scale (Acc-rate=0.60 is preferable).
After having adjusted the parameters for the MCMC algorithm and checking the output we run an analysis.
The output is a list which contains the five arguments posterior, predictive, model, prior and mcmc.input.
Below we show the simulations from the posterior distribution of the signal at a few data locations.
Now we consider an example with a random correlation scale parameter phi and a positive nugget for the random effects S. The program is using a discretised prior fpr phi, where the discretisation is given by
the argument phi.discrete). The argument (tausq.rel = 0.05 gives the relative nugget for S, i.e. the relative
microscale variation).
Acc-rate=0.60 , acc-rate-phi = 0.25-0.30 are preferable.
After having adjusted the parameters for the MCMC algorithm and checking the output we run an analysis.
WARNING: RUNNING THE NEXT COMMAND CAN BE TIME-CONSUMING
Below we show the posterior distribution of the two covariance parameters and the beta parameter.
Construct similar commands using binom.krige.bayes on the data set b50 yourself
(you load the data set by typing data(b50)).
Generate a simulation from a spatial binomial random field.
output = list(threshold = 10, quantile = c(0.05,0.99)))
The posterior contains information
on the posterior distribution of the parameters, and the conditional simulations of the signal g^{-1}(S) at
the data locations.
The predictive contains information on the predictions, where
predictive$median is the predicted signal and predictive$uncertainty is the associated uncertainty.
The threshold = 10 argument gives probabilities of the predictive distribution of the signal being less
than 10 (test5$predictive$probability).
The quantiles = c(0.05,0.99) gives the 0.05 and 0.99 quantiles of the predictive distribution of the signal
(test5$predictive$quantiles).
Exercise
Exercise