Package 'PrevMap'

Adjustment factor for the variance of the convolution of Gaussian noise

Description

This function computes the multiplicative constant used to adjust the value of sigma2 in the low-rank approximation of a Gaussian process.

Usage

adjust.sigma2(knots.dist, phi, kappa)

Arguments

knots.dist	a matrix of the distances between the observed coordinates and the spatial knots.
phi	scale parameter of the Matern covariance function.
kappa	shape parameter of the Matern covariance function.

Details

Let $U$ denote the $n$ by $m$ matrix of the distances between the $n$ observed coordinates and $m$ pre-defined spatial knots. This function computes the following quantity

$\frac{1}{n}\sum_{i=1}^n \sum_{j=1}^m K(u_{ij}; \phi, \kappa)^2,$

where $K(.; \phi, \kappa)$ is the Matern kernel (see matern.kernel) and $u_{ij}$ is the distance between the $i$-th sampled location and the $j$-th spatial knot.

Value

A value corresponding to the adjustment factor for sigma2.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

See Also

matern.kernel, pdist.

---

Plot of the autocorrelgram for posterior samples

Description

Plots the autocorrelogram for the posterior samples of the model parameters and spatial random effects.

Usage

autocor.plot(object, param, component.beta = NULL, component.S = NULL)

Arguments

object	an object of class 'Bayes.PrevMap'.
param	a character indicating for which component of the model the autocorrelation plot is required: param="beta" for the regression coefficients; param="sigma2" for the variance of the spatial random effect; param="phi" for the scale parameter of the Matern correlation function; param="tau2" for the variance of the nugget effect; param="S" for the spatial random effect.
component.beta	if param="beta", component.beta is a numeric value indicating the component of the regression coefficients; default is NULL.
component.S	if param="S", component.S can be a numeric value indicating the component of the spatial random effect, or set equal to "all" if the autocorrelgram should be plotted for all the components. Default is NULL.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

---

Bayesian estimation for the two-levels binary probit model

Description

This function performs Bayesian estimation for a geostatistical binary probit model. It also allows to specify a two-levels model so as to include individual-level and household-level (or any other unit comprising a group of individuals, e.g. village, school, compound, etc...) variables.

Usage

binary.probit.Bayes(
  formula,
  coords,
  data,
  ID.coords,
  control.prior,
  control.mcmc,
  kappa,
  low.rank = FALSE,
  knots = NULL,
  messages = TRUE
)

Arguments

formula	an object of class formula (or one that can be coerced to that class): a symbolic description of the model to be fitted.
coords	an object of class formula indicating the geographic coordinates.
data	a data frame containing the variables in the model.
ID.coords	vector of ID values for the unique set of spatial coordinates obtained from create.ID.coords. These must be provided in order to specify spatial random effects at household-level. Warning: the household coordinates must all be distinct otherwise see jitterDupCoords. Default is NULL.
control.prior	output from control.prior.
control.mcmc	output from control.mcmc.Bayes.
kappa	value for the shape parameter of the Matern covariance function.
low.rank	logical; if low.rank=TRUE a low-rank approximation is required. Default is low.rank=FALSE.
knots	if low.rank=TRUE, knots is a matrix of spatial knots used in the low-rank approximation. Default is knots=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.

Details

This function performs Bayesian estimation for the parameters of the geostatistical binary probit model. Let $i$ and $j$ denote the indices of the $i$-th household and $j$-th individual within that household. The response variable $Y_{ij}$ is a binary indicator taking value 1 if the individual has been tested positive for the disease of interest and 0 otherwise. Conditionally on a zero-mean stationary Gaussian process $S(x_{i})$, $Y_{ij}$ are mutually independent Bernoulli variables with probit link function $\Phi^{-1}(\cdot)$, i.e.

$\Phi^{-1}(p_{ij}) = d_{ij}'\beta + S(x_{i}),$

where $d_{ij}$ is a vector of covariates, both at individual- and household-level, with associated regression coefficients $\beta$. The Gaussian process $S(x)$ has isotropic Matern covariance function (see matern) with variance sigma2, scale parameter phi and shape parameter kappa.

Priors definition. Priors can be defined through the function control.prior. The hierarchical structure of the priors is the following. Let $\theta$ be the vector of the covariance parameters c(sigma2,phi); each component of $\theta$ has independent priors that can be freely defined by the user. However, in control.prior uniform and log-normal priors are also available as default priors for each of the covariance parameters. The vector of regression coefficients beta has a multivariate Gaussian prior with mean beta.mean and covariance matrix beta.covar.

Updating regression coefficents and random effects using auxiliary variables. To update $\beta$ and $S(x_{i})$, we use an auxiliary variable technique based on Rue and Held (2005). Let $V_{ij}$ denote a set of random variables that conditionally on $\beta$ and $S(x_{i})$, are mutually independent Gaussian with mean $d_{ij}'\beta + S(x_{i})$ and unit variance. Then, $Y_{ij}=1$ if $V_{ij} > 0$ and $Y_{ij}=0$ otherwise. Using this representation of the model, we use a Gibbs sampler to simulate from the full conditionals of $\beta$, $S(x_{i})$ and $V_{ij}$. See Section 4.3 of Rue and Held (2005) for more details.

Updating the covariance parameters with a Metropolis-Hastings algorithm. In the MCMC algorithm implemented in binary.probit.Bayes, the transformed parameters

$(\theta_{1}, \theta_{2})=(\log(\sigma^2)/2,\log(\sigma^2/\phi^{2 \kappa}))$

are independently updated using a Metropolis Hastings algorithm. At the $i$-th iteration, a new value is proposed for each parameter from a univariate Gaussian distrubion with variance $h_{i}^2$. This is tuned using the following adaptive scheme

$h_{i} = h_{i-1}+c_{1}i^{-c_{2}}(\alpha_{i}-0.45),$

where $\alpha_{i}$ is the acceptance rate at the $i$-th iteration, 0.45 is the optimal acceptance rate for a univariate Gaussian distribution, whilst $c_{1} > 0$ and $0 < c_{2} < 1$ are pre-defined constants. The starting values $h_{1}$ for each of the parameters $\theta_{1}$ and $\theta_{2}$ can be set using the function control.mcmc.Bayes through the arguments h.theta1, h.theta2 and h.theta3. To define values for $c_{1}$ and $c_{2}$, see the documentation of control.mcmc.Bayes.

Low-rank approximation. In the case of very large spatial data-sets, a low-rank approximation of the Gaussian spatial process $S(x)$ might be computationally beneficial. Let $(x_{1},\dots,x_{m})$ and $(t_{1},\dots,t_{m})$ denote the set of sampling locations and a grid of spatial knots covering the area of interest, respectively. Then $S(x)$ is approximated as $\sum_{i=1}^m K(\|x-t_{i}\|; \phi, \kappa)U_{i}$, where $U_{i}$ are zero-mean mutually independent Gaussian variables with variance sigma2 and $K(.;\phi, \kappa)$ is the isotropic Matern kernel (see matern.kernel). Since the resulting approximation is no longer a stationary process (but only approximately), sigma2 may take very different values from the actual variance of the Gaussian process to approximate. The function adjust.sigma2 can then be used to (approximately) explore the range for sigma2. For example if the variance of the Gaussian process is 0.5, then an approximate value for sigma2 is 0.5/const.sigma2, where const.sigma2 is the value obtained with adjust.sigma2.

Value

An object of class "Bayes.PrevMap". The function summary.Bayes.PrevMap is used to print a summary of the fitted model. The object is a list with the following components:

estimate: matrix of the posterior samples of the model parameters.

S: matrix of the posterior samples for each component of the random effect.

const.sigma2: vector of the values of the multiplicative factor used to adjust the values of sigma2 in the low-rank approximation.

y: binary observations.

D: matrix of covariarates.

coords: matrix of the observed sampling locations.

kappa: shape parameter of the Matern function.

ID.coords: set of ID values defined through the argument ID.coords.

knots: matrix of spatial knots used in the low-rank approximation.

h1: vector of values taken by the tuning parameter h.theta1 at each iteration.

h2: vector of values taken by the tuning parameter h.theta2 at each iteration.

call: the matched call.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

References

Diggle, P.J., Giorgi, E. (2019). Model-based Geostatistics for Global Public Health. CRC/Chapman & Hall.

Giorgi, E., Diggle, P.J. (2017). PrevMap: an R package for prevalence mapping. Journal of Statistical Software. 78(8), 1-29. doi: 10.18637/jss.v078.i08

Rue, H., Held, L. (2005). Gaussian Markov Random Fields: Theory and Applications. Chapman & Hall, London.

Higdon, D. (1998). A process-convolution approach to modeling temperatures in the North Atlantic Ocean. Environmental and Ecological Statistics 5, 173-190.

See Also

control.mcmc.Bayes, control.prior,summary.Bayes.PrevMap, matern, matern.kernel, create.ID.coords.

---

Bayesian estimation for the binomial logistic model

Description

This function performs Bayesian estimation for a geostatistical binomial logistic model.

Usage

binomial.logistic.Bayes(
  formula,
  units.m,
  coords,
  data,
  ID.coords = NULL,
  control.prior,
  control.mcmc,
  kappa,
  low.rank = FALSE,
  knots = NULL,
  messages = TRUE,
  mesh = NULL,
  SPDE = FALSE
)

Arguments

formula	an object of class formula (or one that can be coerced to that class): a symbolic description of the model to be fitted.
units.m	an object of class formula indicating the binomial denominators.
coords	an object of class formula indicating the geographic coordinates.
data	a data frame containing the variables in the model.
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.
control.prior	output from control.prior.
control.mcmc	output from control.mcmc.Bayes.
kappa	value for the shape parameter of the Matern covariance function.
low.rank	logical; if low.rank=TRUE a low-rank approximation is required. Default is low.rank=FALSE.
knots	if low.rank=TRUE, knots is a matrix of spatial knots used in the low-rank approximation. Default is knots=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.
mesh	an object obtained as result of a call to the function inla.mesh.2d.
SPDE	logical; if SPDE=TRUE the SPDE approximation for the Gaussian spatial model is used. Default is SPDE=FALSE.

Details

This function performs Bayesian estimation for the parameters of the geostatistical binomial logistic model. Conditionally on a zero-mean stationary Gaussian process $S(x)$ and mutually independent zero-mean Gaussian variables $Z$ with variance tau2, the linear predictor assumes the form

$\log(p/(1-p)) = d'\beta + S(x) + Z,$

where $d$ is a vector of covariates with associated regression coefficients $\beta$. The Gaussian process $S(x)$ has isotropic Matern covariance function (see matern) with variance sigma2, scale parameter phi and shape parameter kappa.

Priors definition. Priors can be defined through the function control.prior. The hierarchical structure of the priors is the following. Let $\theta$ be the vector of the covariance parameters c(sigma2,phi,tau2); then each component of $\theta$ has independent priors freely defined by the user. However, in control.prior uniform and log-normal priors are also available as default priors for each of the covariance parameters. To remove the nugget effect $Z$, no prior should be defined for tau2. Conditionally on sigma2, the vector of regression coefficients beta has a multivariate Gaussian prior with mean beta.mean and covariance matrix sigma2*beta.covar, while in the low-rank approximation the covariance matrix is simply beta.covar.

Updating the covariance parameters with a Metropolis-Hastings algorithm. In the MCMC algorithm implemented in binomial.logistic.Bayes, the transformed parameters

$(\theta_{1}, \theta_{2}, \theta_{3})=(\log(\sigma^2)/2,\log(\sigma^2/\phi^{2 \kappa}), \log(\tau^2))$

are independently updated using a Metropolis Hastings algorithm. At the $i$-th iteration, a new value is proposed for each from a univariate Gaussian distrubion with variance $h_{i}^2$ that is tuned using the following adaptive scheme

$h_{i} = h_{i-1}+c_{1}i^{-c_{2}}(\alpha_{i}-0.45),$

where $\alpha_{i}$ is the acceptance rate at the $i$-th iteration, 0.45 is the optimal acceptance rate for a univariate Gaussian distribution, whilst $c_{1} > 0$ and $0 < c_{2} < 1$ are pre-defined constants. The starting values $h_{1}$ for each of the parameters $\theta_{1}$, $\theta_{2}$ and $\theta_{3}$ can be set using the function control.mcmc.Bayes through the arguments h.theta1, h.theta2 and h.theta3. To define values for $c_{1}$ and $c_{2}$, see the documentation of control.mcmc.Bayes.

Hamiltonian Monte Carlo. The MCMC algorithm in binomial.logistic.Bayes uses a Hamiltonian Monte Carlo (HMC) procedure to update the random effect $T=d'\beta + S(x) + Z$; see Neal (2011) for an introduction to HMC. HMC makes use of a postion vector, say $t$, representing the random effect $T$, and a momentum vector, say $q$, of the same length of the position vector, say $n$. Hamiltonian dynamics also have a physical interpretation where the states of the system are described by the position of a puck and its momentum (its mass times its velocity). The Hamiltonian function is then defined as a function of $t$ and $q$, having the form $H(t,q) = -\log\{f(t | y, \beta, \theta)\} + q'q/2$, where $f(t | y, \beta, \theta)$ is the conditional distribution of $T$ given the data $y$, the regression parameters $\beta$ and covariance parameters $\theta$. The system of Hamiltonian equations then defines the evolution of the system in time, which can be used to implement an algorithm for simulation from the posterior distribution of $T$. In order to implement the Hamiltonian dynamic on a computer, the Hamiltonian equations must be discretised. The leapfrog method is then used for this purpose, where two tuning parameters should be defined: the stepsize $\epsilon$ and the number of steps $L$. These respectively correspond to epsilon.S.lim and L.S.lim in the control.mcmc.Bayes function. However, it is advisable to let $epsilon$ and $L$ take different random values at each iteration of the HCM algorithm so as to account for the different variances amongst the components of the posterior of $T$. This can be done in control.mcmc.Bayes by defning epsilon.S.lim and L.S.lim as vectors of two elements, each of which represents the lower and upper limit of a uniform distribution used to generate values for epsilon.S.lim and L.S.lim, respectively.

Using a two-level model to include household-level and individual-level information. When analysing data from household sruveys, some of the avilable information information might be at household-level (e.g. material of house, temperature) and some at individual-level (e.g. age, gender). In this case, the Gaussian spatial process $S(x)$ and the nugget effect $Z$ are defined at hosuehold-level in order to account for extra-binomial variation between and within households, respectively.

Low-rank approximation. In the case of very large spatial data-sets, a low-rank approximation of the Gaussian spatial process $S(x)$ might be computationally beneficial. Let $(x_{1},\dots,x_{m})$ and $(t_{1},\dots,t_{m})$ denote the set of sampling locations and a grid of spatial knots covering the area of interest, respectively. Then $S(x)$ is approximated as $\sum_{i=1}^m K(\|x-t_{i}\|; \phi, \kappa)U_{i}$, where $U_{i}$ are zero-mean mutually independent Gaussian variables with variance sigma2 and $K(.;\phi, \kappa)$ is the isotropic Matern kernel (see matern.kernel). Since the resulting approximation is no longer a stationary process (but only approximately), sigma2 may take very different values from the actual variance of the Gaussian process to approximate. The function adjust.sigma2 can then be used to (approximately) explore the range for sigma2. For example if the variance of the Gaussian process is 0.5, then an approximate value for sigma2 is 0.5/const.sigma2, where const.sigma2 is the value obtained with adjust.sigma2.

Value

An object of class "Bayes.PrevMap". The function summary.Bayes.PrevMap is used to print a summary of the fitted model. The object is a list with the following components:

estimate: matrix of the posterior samples of the model parameters.

S: matrix of the posterior samples for each component of the random effect.

const.sigma2: vector of the values of the multiplicative factor used to adjust the values of sigma2 in the low-rank approximation.

y: binomial observations.

units.m: binomial denominators.

D: matrix of covariarates.

coords: matrix of the observed sampling locations.

kappa: shape parameter of the Matern function.

ID.coords: set of ID values defined through the argument ID.coords.

knots: matrix of spatial knots used in the low-rank approximation.

h1: vector of values taken by the tuning parameter h.theta1 at each iteration.

h2: vector of values taken by the tuning parameter h.theta2 at each iteration.

h3: vector of values taken by the tuning parameter h.theta3 at each iteration.

acc.beta.S: empirical acceptance rate for the regression coefficients and random effects (only if SPDE=TRUE).

mesh: the mesh used in the SPDE approximation.

call: the matched call.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

References

Diggle, P.J., Giorgi, E. (2019). Model-based Geostatistics for Global Public Health. CRC/Chapman & Hall.

Giorgi, E., Diggle, P.J. (2017). PrevMap: an R package for prevalence mapping. Journal of Statistical Software. 78(8), 1-29. doi: 10.18637/jss.v078.i08

Neal, R. M. (2011) MCMC using Hamiltonian Dynamics, In: Handbook of Markov Chain Monte Carlo (Chapter 5), Edited by Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng Chapman & Hall / CRC Press.

Higdon, D. (1998). A process-convolution approach to modeling temperatures in the North Atlantic Ocean. Environmental and Ecological Statistics 5, 173-190.

See Also

control.mcmc.Bayes, control.prior,summary.Bayes.PrevMap, matern, matern.kernel, create.ID.coords.

---

Monte Carlo Maximum Likelihood estimation for the binomial logistic model

Description

This function performs Monte Carlo maximum likelihood (MCML) estimation for the geostatistical binomial logistic model.

Usage

binomial.logistic.MCML(
  formula,
  units.m,
  coords,
  times = NULL,
  data,
  ID.coords = NULL,
  par0,
  control.mcmc,
  kappa,
  kappa.t = NULL,
  sst.model = NULL,
  fixed.rel.nugget = NULL,
  start.cov.pars,
  method = "BFGS",
  low.rank = FALSE,
  SPDE = FALSE,
  knots = NULL,
  mesh = NULL,
  messages = TRUE,
  plot.correlogram = TRUE
)

Arguments

formula	an object of class formula (or one that can be coerced to that class): a symbolic description of the model to be fitted.
units.m	an object of class formula indicating the binomial denominators in the data.
coords	an object of class formula indicating the spatial coordinates in the data.
times	an object of class formula indicating the times in the data, used in the spatio-temporal model.
data	a data frame containing the variables in the model.
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.
par0	parameters of the importance sampling distribution: these should be given in the following order c(beta,sigma2,phi,tau2), where beta are the regression coefficients, sigma2 is the variance of the Gaussian process, phi is the scale parameter of the spatial correlation and tau2 is the variance of the nugget effect (if included in the model).
control.mcmc	output from control.mcmc.MCML.
kappa	fixed value for the shape parameter of the Matern covariance function.
kappa.t	fixed value for the shape parameter of the Matern covariance function in the separable double-Matern spatio-temporal model.
sst.model	a character value that specifies the spatio-temporal correlation function. sst.model="DM" separable double-Matern. sst.model="GN1" separable correlation functions. Temporal correlation: $f(x) = 1/(1+x/\psi)$; Spatial correaltion: Matern function. Deafault is sst.model=NULL, which is used when a purely spatial model is fitted.
fixed.rel.nugget	fixed value for the relative variance of the nugget effect; fixed.rel.nugget=NULL if this should be included in the estimation. Default is fixed.rel.nugget=NULL.
start.cov.pars	a vector of length two with elements corresponding to the starting values of phi and the relative variance of the nugget effect nu2, respectively, that are used in the optimization algorithm. If nu2 is fixed through fixed.rel.nugget, then start.cov.pars represents the starting value for phi only.
method	method of optimization. If method="BFGS" then the maxBFGS function is used; otherwise method="nlminb" to use the nlminb function. Default is method="BFGS".
low.rank	logical; if low.rank=TRUE a low-rank approximation of the Gaussian spatial process is used when fitting the model. Default is low.rank=FALSE.
SPDE	logical; if SPDE=TRUE the SPDE approximation for the Gaussian spatial model is used. Default is SPDE=FALSE.
knots	if low.rank=TRUE, knots is a matrix of spatial knots that are used in the low-rank approximation. Default is knots=NULL.
mesh	an object obtained as result of a call to the function inla.mesh.2d.
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 samples of the random effect is displayed after completion of conditional simulation. Default is plot.correlogram=TRUE.

Details

This function performs parameter estimation for a geostatistical binomial logistic model. Conditionally on a zero-mean stationary Gaussian process $S(x)$ and mutually independent zero-mean Gaussian variables $Z$ with variance tau2, the observations y are generated from a binomial distribution with probability $p$ and binomial denominators units.m. A canonical logistic link is used, thus the linear predictor assumes the form

$\log(p/(1-p)) = d'\beta + S(x) + Z,$

where $d$ is a vector of covariates with associated regression coefficients $\beta$. The Gaussian process $S(x)$ has isotropic Matern covariance function (see matern) with variance sigma2, scale parameter phi and shape parameter kappa. In the binomial.logistic.MCML function, the shape parameter is treated as fixed. The relative variance of the nugget effect, nu2=tau2/sigma2, can also be fixed through the argument fixed.rel.nugget; if fixed.rel.nugget=NULL, then the relative variance of the nugget effect is also included in the estimation.

Monte Carlo Maximum likelihood. The Monte Carlo maximum likelihood method uses conditional simulation from the distribution of the random effect $T(x) = d(x)'\beta+S(x)+Z$ given the data y, in order to approximate the high-dimensiional intractable integral given by the likelihood function. The resulting approximation of the likelihood is then maximized by a numerical optimization algorithm which uses analytic epression for computation of the gradient vector and Hessian matrix. The functions used for numerical optimization are maxBFGS (method="BFGS"), from the maxLik package, and nlminb (method="nlminb").

Using a two-level model to include household-level and individual-level information. When analysing data from household sruveys, some of the avilable information information might be at household-level (e.g. material of house, temperature) and some at individual-level (e.g. age, gender). In this case, the Gaussian spatial process $S(x)$ and the nugget effect $Z$ are defined at hosuehold-level in order to account for extra-binomial variation between and within households, respectively.

Low-rank approximation. In the case of very large spatial data-sets, a low-rank approximation of the Gaussian spatial process $S(x)$ might be computationally beneficial. Let $(x_{1},\dots,x_{m})$ and $(t_{1},\dots,t_{m})$ denote the set of sampling locations and a grid of spatial knots covering the area of interest, respectively. Then $S(x)$ is approximated as $\sum_{i=1}^m K(\|x-t_{i}\|; \phi, \kappa)U_{i}$, where $U_{i}$ are zero-mean mutually independent Gaussian variables with variance sigma2 and $K(.;\phi, \kappa)$ is the isotropic Matern kernel (see matern.kernel). Since the resulting approximation is no longer a stationary process (but only approximately), the parameter sigma2 is then multiplied by a factor constant.sigma2 so as to obtain a value that is closer to the actual variance of $S(x)$.

Value

An object of class "PrevMap". The function summary.PrevMap is used to print a summary of the fitted model. The object is a list with the following components:

estimate: estimates of the model parameters; use the function coef.PrevMap to obtain estimates of covariance parameters on the original scale.

covariance: covariance matrix of the MCML estimates.

log.lik: maximum value of the log-likelihood.

y: binomial observations.

units.m: binomial denominators.

D: matrix of covariates.

coords: matrix of the observed sampling locations.

method: method of optimization used.

ID.coords: set of ID values defined through the argument ID.coords.

kappa: fixed value of the shape parameter of the Matern function.

kappa.t: fixed value for the shape parameter of the Matern covariance function in the separable double-Matern spatio-temporal model.

knots: matrix of the spatial knots used in the low-rank approximation.

mesh: the mesh used in the SPDE approximation.

const.sigma2: adjustment factor for sigma2 in the low-rank approximation.

h: vector of the values of the tuning parameter at each iteration of the Langevin-Hastings MCMC algorithm; see Laplace.sampling, or Laplace.sampling.lr if a low-rank approximation is used.

samples: matrix of the random effects samples from the importance sampling distribution used to approximate the likelihood function.

fixed.rel.nugget: fixed value for the relative variance of the nugget effect.

call: the matched call.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

References

Diggle, P.J., Giorgi, E. (2019). Model-based Geostatistics for Global Public Health. CRC/Chapman & Hall.

Giorgi, E., Diggle, P.J. (2017). PrevMap: an R package for prevalence mapping. Journal of Statistical Software. 78(8), 1-29. doi: 10.18637/jss.v078.i08

Christensen, O. F. (2004). Monte carlo maximum likelihood in model-based geostatistics. Journal of Computational and Graphical Statistics 13, 702-718.

Higdon, D. (1998). A process-convolution approach to modeling temperatures in the North Atlantic Ocean. Environmental and Ecological Statistics 5, 173-190.

See Also

Laplace.sampling, Laplace.sampling.lr, summary.PrevMap, coef.PrevMap, matern, matern.kernel, control.mcmc.MCML, create.ID.coords.

---

Extract model coefficients

Description

coef extracts parameters estimates from models fitted with the functions linear.model.MLE and binomial.logistic.MCML.

Usage

## S3 method for class 'PrevMap'
coef(object, ...)

Arguments

object	an object of class "PrevMap".
...	other arguments.

Value

coefficients extracted from the model object object.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

---

Extract model coefficients from geostatistical linear model with preferentially sampled locations

Description

coef extracts parameters estimates from models fitted with the functions lm.ps.MCML.

Usage

## S3 method for class 'PrevMap.ps'
coef(object, ...)

Arguments

object	an object of class "PrevMap.ps".
...	other arguments.

Value

a list of coefficients extracted from the model in object.

Author(s)

Emanuele Giorgi [email protected]

---

Spatially continuous sampling

Description

Draws a sample of spatial locations within a spatially continuous polygonal sampling region.

Usage

continuous.sample(poly, n, delta, k = 0, rho = NULL)

Arguments

poly	boundary of a polygon.
n	number of events.
delta	minimum permissible distance between any two events in preliminary sample.
k	number of locations in preliminary sample to be replaced by near neighbours of other preliminary sample locations in final sample (must be between 0 and n/2)
rho	maximum distance between close pairs of locations in final sample.

Details

To draw a sample of size n from a spatially continuous region $A$, with the property that the distance between any two sampled locations is at least delta, the following algorithm is used.

• Step 1. Set $i = 1$ and generate a point $x_{1}$ uniformly distributed on $A$.

• Step 2. Increase $i$ by 1, generate a point $x_{i}$ uniformly distributed on $A$ and calculate the minimum, $d_{\min}$, of the distances from $x_{i}$ to all $x_{j}: j < i$.

• Step 3. If $d_{\min} \ge \delta$, increase $i$ by 1 and return to step 2 if $i \le n$, otherwise stop;

• Step 4. If $d_{\min} < \delta$, return to step 2 without increasing $i$.

Sampling close pairs of points. For some purposes, it is desirable that a spatial sampling scheme include pairs of closely spaced points. In this case, the above algorithm requires the following additional steps to be taken. Let k be the required number of close pairs. Choose a value rho such that a close pair of points will be a pair of points separated by a distance of at most rho.

• Step 5. Set $j = 1$ and draw a random sample of size 2 from the integers $1,2,\ldots,n$, say $(i_{1}; i_{2})$;

• Step 6. Replace $x_{i_{1}}$ by $x_{i_{2}} + u$ , where $u$ is uniformly distributed on the disc with centre $x_{i_{2}}$ and radius rho, increase $i$ by 1 and return to step 5 if $i \le k$, otherwise stop.

Value

A matrix of dimension n by 2 containing event locations.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

---

Contour plot of a predicted surface

Description

plot.pred.PrevMap displays contours of predictions obtained from spatial.pred.linear.MLE, spatial.pred.linear.Bayes,spatial.pred.binomial.MCML and spatial.pred.binomial.Bayes.

Usage

## S3 method for class 'pred.PrevMap'
contour(x, type = NULL, summary = "predictions", ...)

Arguments

x	an object of class "pred.PrevMap".
type	a character indicating the type of prediction to display: 'prevalence', 'odds', 'logit' or 'probit'.
summary	character indicating which summary to display: 'predictions','quantiles', 'standard.errors' or 'exceedance.prob'; default is 'predictions'. If summary="exceedance.prob", the argument type is ignored.
...	further arguments passed to contour.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

---

Control settings for the MCMC algorithm used for Bayesian inference

Description

This function defines the different tuning parameter that are used in the MCMC algorithm for Bayesian inference.

Usage

control.mcmc.Bayes(
  n.sim,
  burnin,
  thin,
  h.theta1 = 0.01,
  h.theta2 = 0.01,
  h.theta3 = 0.01,
  L.S.lim = NULL,
  epsilon.S.lim = NULL,
  start.beta = "prior mean",
  start.sigma2 = "prior mean",
  start.phi = "prior mean",
  start.S = "prior mean",
  start.nugget = "prior mean",
  c1.h.theta1 = 0.01,
  c2.h.theta1 = 1e-04,
  c1.h.theta2 = 0.01,
  c2.h.theta2 = 1e-04,
  c1.h.theta3 = 0.01,
  c2.h.theta3 = 1e-04,
  linear.model = FALSE,
  binary = FALSE
)

Arguments

n.sim	total number of simulations.
burnin	initial number of samples to be discarded.
thin	value used to retain only evey thin-th sampled value.
h.theta1	starting value of the tuning parameter of the proposal distribution for $\theta_{1} = \log(\sigma^2)/2$. See 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
h.theta2	starting value of the tuning parameter of the proposal distribution for $\theta_{2} = \log(\sigma^2/\phi^{2 \kappa})$. See 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
h.theta3	starting value of the tuning parameter of the proposal distribution for $\theta_{3} = \log(\tau^2)$. See 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
L.S.lim	an atomic value or a vector of length 2 that is used to define the number of steps used at each iteration in the Hamiltonian Monte Carlo algorithm to update the spatial random effect; if a single value is provided than the number of steps is kept fixed, otherwise if a vector of length 2 is provided the number of steps is simulated at each iteration as floor(runif(1,L.S.lim[1],L.S.lim[2]+1)).
epsilon.S.lim	an atomic value or a vector of length 2 that is used to define the stepsize used at each iteration in the Hamiltonian Monte Carlo algorithm to update the spatial random effect; if a single value is provided than the stepsize is kept fixed, otherwise if a vector of length 2 is provided the stepsize is simulated at each iteration as runif(1,epsilon.S.lim[1],epsilon.S.lim[2]).
start.beta	starting value for the regression coefficients beta.
start.sigma2	starting value for sigma2.
start.phi	starting value for phi.
start.S	starting value for the spatial random effect.
start.nugget	starting value for the variance of the nugget effect; default is NULL if the nugget effect is not present.
c1.h.theta1	value of $c_{1}$ used to adaptively tune the variance of the Gaussian proposal for the transformed parameter log(sigma2)/2; see 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
c2.h.theta1	value of $c_{2}$ used to adaptively tune the variance of the Gaussian proposal for the transformed parameter log(sigma2)/2; see 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
c1.h.theta2	value of $c_{1}$ used to adaptively tune the variance of the Gaussian proposal for the transformed parameter log(sigma2.curr/(phi.curr^(2*kappa))); see 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
c2.h.theta2	value of $c_{2}$ used to adaptively tune the variance of the Gaussian proposal for the transformed parameter log(sigma2.curr/(phi.curr^(2*kappa))); see 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
c1.h.theta3	value of $c_{1}$ used to adaptively tune the variance of the Gaussian proposal for the transformed parameter log(tau2); see 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
c2.h.theta3	value of $c_{2}$ used to adaptively tune the variance of the Gaussian proposal for the transformed parameter log(tau2); see 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
linear.model	logical; if linear.model=TRUE, the control parameters are set for the geostatistical linear model. Default is linear.model=FALSE.
binary	logical; if binary=TRUE, the control parameters are set the binary geostatistical model. Default is binary=FALSE.

Value

an object of class "mcmc.Bayes.PrevMap".

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

---

Control settings for the MCMC algorithm used for Bayesian inference using SPDE

Description

This function defines the different tuning parameter that are used in the MCMC algorithm for Bayesian inference using a SPDE approximation for the spatial Gaussian process.

Usage

control.mcmc.Bayes.SPDE(
  n.sim,
  burnin,
  thin,
  h.theta1 = 0.01,
  h.theta2 = 0.01,
  start.beta = "prior mean",
  start.sigma2 = "prior mean",
  start.phi = "prior mean",
  start.S = "prior mean",
  n.iter = 1,
  h = 1,
  c1.h.theta1 = 0.01,
  c2.h.theta1 = 1e-04,
  c1.h.theta2 = 0.01,
  c2.h.theta2 = 1e-04
)

Arguments

n.sim	total number of simulations.
burnin	initial number of samples to be discarded.
thin	value used to retain only evey thin-th sampled value.
h.theta1	starting value of the tuning parameter of the proposal distribution for $\theta_{1} = \log(\sigma^2)/2$. See 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
h.theta2	starting value of the tuning parameter of the proposal distribution for $\theta_{2} = \log(\sigma^2/\phi^{2 \kappa})$. See 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
start.beta	starting value for the regression coefficients beta. If not provided the prior mean is used.
start.sigma2	starting value for sigma2. If not provided the prior mean is used.
start.phi	starting value for phi. If not provided the prior mean is used.
start.S	starting value for the spatial random effect. If not provided the prior mean is used.
n.iter	number of iteration of the Newton-Raphson procedure used to compute the mean and coviariance matrix of the Gaussian proposal in the MCMC; defaut is n.iter=1.
h	tuning parameter for the covariance matrix of the Gaussian proposal. Default is h=1.
c1.h.theta1	value of $c_{1}$ used to adaptively tune the variance of the Gaussian proposal for the transformed parameter log(sigma2)/2; see 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
c2.h.theta1	value of $c_{2}$ used to adaptively tune the variance of the Gaussian proposal for the transformed parameter log(sigma2)/2; see 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
c1.h.theta2	value of $c_{1}$ used to adaptively tune the variance of the Gaussian proposal for the transformed parameter log(sigma2.curr/(phi.curr^(2*kappa))); see 'Details' in binomial.logistic.Bayes or linear.model.Bayes.
c2.h.theta2	value of $c_{2}$ used to adaptively tune the variance of the Gaussian proposal for the transformed parameter log(sigma2.curr/(phi.curr^(2*kappa))); see 'Details' in binomial.logistic.Bayes or linear.model.Bayes.

Value

an object of class "mcmc.Bayes.PrevMap".

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

---

Control settings for the MCMC algorithm used for classical inference on a binomial logistic model

Description

This function defines the options for the MCMC algorithm used in the Monte Carlo maximum likelihood method.

Usage

control.mcmc.MCML(n.sim, burnin, thin = 1, h = NULL, c1.h = 0.01, c2.h = 1e-04)

Arguments

n.sim	number of simulations.
burnin	length of the burn-in period.
thin	only every thin iterations, a sample is stored; default is thin=1.
h	tuning parameter of the proposal distribution used in the Langevin-Hastings MCMC algorithm (see Laplace.sampling and Laplace.sampling.lr); default is h=NULL and then set internally as $1.65/n^(1/6)$, where $n$ is the dimension of the random effect.
c1.h	value of $c_{1}$ used in the adaptive scheme for h; default is c1.h=0.01. See also 'Details' in binomial.logistic.MCML
c2.h	value of $c_{2}$ used in the adaptive scheme for h; default is c1.h=0.01. See also 'Details' in binomial.logistic.MCML

Value

A list with processed arguments to be passed to the main function.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

Examples

control.mcmc <- control.mcmc.MCML(n.sim=1000,burnin=100,thin=1,h=0.05)
str(control.mcmc)

---

Priors specification

Description

This function is used to define priors for the model parameters of a Bayesian geostatistical model.

Usage

control.prior(
  beta.mean,
  beta.covar,
  log.prior.sigma2 = NULL,
  log.prior.phi = NULL,
  log.prior.nugget = NULL,
  uniform.sigma2 = NULL,
  log.normal.sigma2 = NULL,
  uniform.phi = NULL,
  log.normal.phi = NULL,
  uniform.nugget = NULL,
  log.normal.nugget = NULL
)

Arguments

beta.mean	mean vector of the Gaussian prior for the regression coefficients.
beta.covar	covariance matrix of the Gaussian prior for the regression coefficients.
log.prior.sigma2	a function corresponding to the log-density of the prior distribution for the variance sigma2 of the Gaussian process. Warning: if a low-rank approximation is used, then sigma2 corresponds to the variance of the iid zero-mean Gaussian variables. Default is NULL.
log.prior.phi	a function corresponding to the log-density of the prior distribution for the scale parameter of the Matern correlation function; default is NULL.
log.prior.nugget	optional: a function corresponding to the log-density of the prior distribution for the variance of the nugget effect; default is NULL with no nugget incorporated in the model; default is NULL.
uniform.sigma2	a vector of length two, corresponding to the lower and upper limit of the uniform prior on sigma2. Default is NULL.
log.normal.sigma2	a vector of length two, corresponding to the mean and standard deviation of the distribution on the log scale for the log-normal prior on sigma2. Default is NULL.
uniform.phi	a vector of length two, corresponding to the lower and upper limit of the uniform prior on phi. Default is NULL.
log.normal.phi	a vector of length two, corresponding to the mean and standard deviation of the distribution on the log scale for the log-normal prior on phi. Default is NULL.
uniform.nugget	a vector of length two, corresponding to the lower and upper limit of the uniform prior on tau2. Default is NULL.
log.normal.nugget	a vector of length two, corresponding to the mean and standard deviation of the distribution on the log scale for the log-normal prior on tau2. Default is NULL.

Value

a list corresponding the prior distributions for each model parameter.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

See Also

See "Priors definition" in the Details section of the binomial.logistic.Bayes function.

---

Auxliary function for controlling profile log-likelihood in the linear Gaussian model

Description

Auxiliary function used by loglik.linear.model. This function defines whether the profile-loglikelihood should be computed or evaluation of the likelihood is required by keeping the other parameters fixed.

Usage

control.profile(
  phi = NULL,
  rel.nugget = NULL,
  fixed.beta = NULL,
  fixed.sigma2 = NULL,
  fixed.phi = NULL,
  fixed.rel.nugget = NULL
)

Arguments

phi	a vector of the different values that should be used in the likelihood evalutation for the scale parameter phi, or NULL if a single value is provided either as first argument in start.par (for profile likelihood maximization) or as fixed value in fixed.phi; default is NULL.
rel.nugget	a vector of the different values that should be used in the likelihood evalutation for the relative variance of the nugget effect nu2, or NULL if a single value is provided either in start.par (for profile likelihood maximization) or as fixed value in fixed.nu2; default is NULL.
fixed.beta	a vector for the fixed values of the regression coefficients beta, or NULL if profile log-likelihood is to be performed; default is NULL.
fixed.sigma2	value for the fixed variance of the Gaussian process sigma2, or NULL if profile log-likelihood is to be performed; default is NULL.
fixed.phi	value for the fixed scale parameter phi in the Matern function, or NULL if profile log-likelihood is to be performed; default is NULL.
fixed.rel.nugget	value for the fixed relative variance of the nugget effect; fixed.rel.nugget=NULL if profile log-likelihood is to be performed; default is NULL.

Value

A list with components named as the arguments.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

See Also

loglik.linear.model

---

ID spatial coordinates

Description

Creates ID values for the unique set of coordinates.

Usage

create.ID.coords(data, coords)

Arguments

data	a data frame containing the spatial coordinates.
coords	an object of class formula indicating the geographic coordinates.

Value

a vector of integers indicating the corresponding rows in data for each distinct coordinate obtained with the unique function.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

Examples

x1 <- runif(5)
x2 <- runif(5)
data <- data.frame(x1=rep(x1,each=3),x2=rep(x2,each=3))
ID.coords <- create.ID.coords(data,coords=~x1+x2)
data[,c("x1","x2")]==unique(data[,c("x1","x2")])[ID.coords,]

---

Simulated binomial data-set over the unit square

Description

This binomial data-set was simulated by generating a zero-mean Gaussian process over a 30 by 30 grid covering the unit square. The parameters used in the simulation are sigma2=1, phi=0.15 and kappa=2. The nugget effect was not included, hence tau2=0. The variables are as follows:

• y binomial observations.

• units.m binomial denominators.

• x1 horizontal coordinates.

• x2 vertical coordinates.

• S simulated values of the Gaussian process.

Usage

data(data_sim)

Format

A data frame with 900 rows and 5 variables

---

Density plot for posterior samples

Description

Plots the autocorrelogram for the posterior samples of the model parameters and spatial random effects.

Usage

dens.plot(
  object,
  param,
  component.beta = NULL,
  component.S = NULL,
  hist = TRUE,
  ...
)

Arguments

object	an object of class 'Bayes.PrevMap'.
param	a character indicating for which component of the model the density plot is required: param="beta" for the regression coefficients; param="sigma2" for the variance of the spatial random effect; param="phi" for the scale parameter of the Matern correlation function; param="tau2" for the variance of the nugget effect; param="S" for the spatial random effect.
component.beta	if param="beta", component.beta is a numeric value indicating the component of the regression coefficients; default is NULL.
component.S	if param="S", component.S can be a numeric value indicating the component of the spatial random effect. Default is NULL.
hist	logical; if TRUE a histrogram is added to density plot.
...	additional parameters to pass to density.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

---

Spatially discrete sampling

Description

Draws a sub-sample from a set of units spatially located irregularly over some defined geographical region by imposing a minimum distance between any two sampled units.

Usage

discrete.sample(xy.all, n, delta, k = 0)

Arguments

xy.all	set of locations from which the sample will be drawn.
n	size of required sample.
delta	minimum distance between any two locations in preliminary sample.
k	number of locations in preliminary sample to be replaced by nearest neighbours of other preliminary sample locations in final sample (must be between 0 and n/2).

Details

To draw a sample of size n from a population of spatial locations $X_{i} : i = 1,\ldots,N$, with the property that the distance between any two sampled locations is at least delta, the function implements the following algorithm.

• Step 1. Draw an initial sample of size n completely at random and call this $x_{i} : i = 1,\dots, n$.

• Step 2. Set $i = 1$ and calculate the minimum, $d_{\min}$, of the distances from $x_{i}$ to all other $x_{j}$ in the initial sample.

• Step 3. If $d_{\min} \ge \delta$, increase $i$ by 1 and return to step 2 if $i \le n$, otherwise stop.

• Step 4. If $d_{\min} < \delta$, draw an integer $j$ at random from $1, 2,\ldots,N$, set $x_{i} = X_{j}$ and return to step 3.

Samples generated in this way will exhibit a more regular spatial arrangement than would a random sample of the same size. The degree of regularity achievable will be influenced by the spatial arrangement of the population $X_{i} : i = 1,\ldots,N$, the specified value of delta and the sample size n. For any given population, if n and/or delta are too large, a sample of the required size with the distance between any two sampled locations at least delta will not be achievable; the suggested solution is then to run the algorithm with a smaller value of delta.

Sampling close pairs of points. For some purposes, it is desirable that a spatial sampling scheme include pairs of closely spaced points. In this case, the above algorithm requires the following additional steps to be taken. Let k be the required number of close pairs.

• Step 5. Set $j = 1$ and draw a random sample of size 2 from the integers $1, 2,\ldots,n$, say $(i_{1}, i_{2} )$.

• Step 6. Find the integer $r$ such that the distances from $x_{i_{1}}$ to $X_{r}$ is the minimum of all $N-1$ distances from $x_{i_{1}}$ to the $X_{j}$.

• Step 7. Replace $x_{i_{2}}$ by $X_{r}$, increase $i$ by 1 and return to step 5 if $i \le k$, otherwise stop.

Value

A matrix of dimension n by 2 containing the final sampled locations.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

Examples

x<-0.015+0.03*(1:33)
xall<-rep(x,33)
yall<-c(t(matrix(xall,33,33)))
xy<-cbind(xall,yall)+matrix(-0.0075+0.015*runif(33*33*2),33*33,2)
par(pty="s",mfrow=c(1,2))
plot(xy[,1],xy[,2],pch=19,cex=0.25,xlab="Easting",ylab="Northing",
   cex.lab=1,cex.axis=1,cex.main=1)

set.seed(15892)
# Generate spatially random sample
xy.sample<-xy[sample(1:dim(xy)[1],50,replace=FALSE),]
points(xy.sample[,1],xy.sample[,2],pch=19,col="red")
points(xy[,1],xy[,2],pch=19,cex=0.25)
plot(xy[,1],xy[,2],pch=19,cex=0.25,xlab="Easting",ylab="Northing",
   cex.lab=1,cex.axis=1,cex.main=1)

set.seed(15892)
# Generate spatially regular sample
xy.sample<-discrete.sample(xy,50,0.08)
points(xy.sample[,1],xy.sample[,2],pch=19,col="red")
points(xy[,1],xy[,2],pch=19,cex=0.25)

---

Heavy metal biomonitoring in Galicia

Description

This data-set relates to two studies on lead concentration in moss samples, in micrograms per gram dry weight, collected in Galicia, norther Spain. The data are from two surveys, one conducted in October 1997 and on in July 2000. The variables are as follows:

• x x-coordinate of the spatial locations.

• y y-coordinate of the spatial locations.

• lead lead concentration.

• survey year of the survey (either 1997 or 2000).

Usage

data(galicia)

Format

A data frame with 195 rows and 4 variables

Source

Diggle, P.J., Menezes, R. and Su, T.-L. (2010). Geostatistical analysis under preferential sampling (with Discussion). Applied Statistics, 59, 191-232.

---

Boundary of Galicia

Description

This data-set contains the geographical coordinates of the boundary of the Galicia region in northern Spain.

The variables are as follows:

• x x-coordinate of the spatial locations.

• y y-coordinate of the spatial locations.

Usage

data(galicia.boundary)

Format

A data frame with 42315 rows and 2 variables

---

Maximum Likelihood estimation for generalised linear geostatistical models via the Laplace approximation

Description

This function performs the Laplace method for maximum likelihood estimation of a generalised linear geostatistical model.

Usage

glgm.LA(
  formula,
  units.m = NULL,
  coords,
  times = NULL,
  data,
  ID.coords = NULL,
  kappa,
  kappa.t = 0.5,
  fixed.rel.nugget = NULL,
  start.cov.pars,
  method = "nlminb",
  messages = TRUE,
  family,
  return.covariance = TRUE
)

Arguments

formula	an object of class formula (or one that can be coerced to that class): a symbolic description of the model to be fitted.
units.m	an object of class formula indicating the binomial denominators in the data.
coords	an object of class formula indicating the spatial coordinates in the data.
times	an object of class formula indicating the times in the data, used in the spatio-temporal model.
data	a data frame containing the variables in the model.
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.
kappa	fixed value for the shape parameter of the Matern covariance function.
kappa.t	fixed value for the shape parameter of the Matern covariance function in the separable double-Matern spatio-temporal model.
fixed.rel.nugget	fixed value for the relative variance of the nugget effect; fixed.rel.nugget=NULL if this should be included in the estimation. Default is fixed.rel.nugget=NULL.
start.cov.pars	a vector of length two with elements corresponding to the starting values of phi and the relative variance of the nugget effect nu2, respectively, that are used in the optimization algorithm. If nu2 is fixed through fixed.rel.nugget, then start.cov.pars represents the starting value for phi only.
method	method of optimization. If method="BFGS" then the maxBFGS function is used; otherwise method="nlminb" to use the nlminb function. Default is method="BFGS".
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.
family	character, indicating the conditional distribution of the outcome. This should be "Gaussian", "Binomial" or "Poisson".
return.covariance	logical; if return.covariance=TRUE then a numerical estimation of the covariance function for the model parameters is returned. Default is return.covariance=TRUE.

Details

This function performs parameter estimation for a generealized linear geostatistical model. Conditionally on a zero-mean stationary Gaussian process $S(x)$ and mutually independent zero-mean Gaussian variables $Z$ with variance tau2, the observations y are generated from a GLM with link function $g(.)$ and linear predictor

$\eta = d'\beta + S(x) + Z,$

where $d$ is a vector of covariates with associated regression coefficients $\beta$. The Gaussian process $S(x)$ has isotropic Matern covariance function (see matern) with variance sigma2, scale parameter phi and shape parameter kappa. The shape parameter is treated as fixed. The relative variance of the nugget effect, nu2=tau2/sigma2, can also be fixed through the argument fixed.rel.nugget; if fixed.rel.nugget=NULL, then the relative variance of the nugget effect is also included in the estimation.

Laplace Approximation The Laplace approximation (LA) method uses a second-order Taylor expansion of the integrand expressing the likelihood function. The resulting approximation of the likelihood is then maximized by a numerical optimization as defined through the argument method.

Using a two-level model to include household-level and individual-level information. When analysing data from household sruveys, some of the avilable information information might be at household-level (e.g. material of house, temperature) and some at individual-level (e.g. age, gender). In this case, the Gaussian spatial process $S(x)$ and the nugget effect $Z$ are defined at hosuehold-level in order to account for extra-binomial variation between and within households, respectively.

Value

An object of class "PrevMap". The function summary.PrevMap is used to print a summary of the fitted model. The object is a list with the following components:

estimate: estimates of the model parameters; use the function coef.PrevMap to obtain estimates of covariance parameters on the original scale.

covariance: covariance matrix of the MCML estimates.

log.lik: maximum value of the log-likelihood.

y: binomial observations.

units.m: binomial denominators.

D: matrix of covariates.

coords: matrix of the observed sampling locations.

times: vector of the time points used in a spatio-temporal model.

method: method of optimization used.

ID.coords: set of ID values defined through the argument ID.coords.

kappa: fixed value of the shape parameter of the Matern function.

kappa.t: fixed value for the shape parameter of the Matern covariance function in the separable double-Matern spatio-temporal model.

fixed.rel.nugget: fixed value for the relative variance of the nugget effect.

call: the matched call.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

References

Diggle, P.J., Giorgi, E. (2019). Model-based Geostatistics for Global Public Health. CRC/Chapman & Hall.

Giorgi, E., Diggle, P.J. (2017). PrevMap: an R package for prevalence mapping. Journal of Statistical Software. 78(8), 1-29. doi: 10.18637/jss.v078.i08

Christensen, O. F. (2004). Monte carlo maximum likelihood in model-based geostatistics. Journal of Computational and Graphical Statistics 13, 702-718.

Higdon, D. (1998). A process-convolution approach to modeling temperatures in the North Atlantic Ocean. Environmental and Ecological Statistics 5, 173-190.

See Also

Laplace.sampling, Laplace.sampling.lr, summary.PrevMap, coef.PrevMap, matern, matern.kernel, control.mcmc.MCML, create.ID.coords.

---

Langevin-Hastings MCMC for conditional simulation

Description

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

Usage

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

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/Poisson observations.
units.m	vector of binomial denominators, or offset if the Poisson model is used.
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.
poisson.llik	logical; if poisson.llik=TRUE a Poisson model is used or, if poisson.llik=FALSE, a binomial model is used.

Details

Binomial model. 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.$

Poisson model. Conditionally on the random effect $S$, the data y follow a Poisson distribution with mean $m\lambda$, where $m$ is an offset set through the argument units.m. The log link function is used for the linear predictor, which assumes the form

$\log(\lambda)=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{\Sigma}^{-1/2}(S-\tilde{s})$ where $\tilde{\Sigma}$ 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}}(\alpha_{i}-0.547),$

where $c_{1} > 0$ and $0 < c_{2} < 1$ are pre-defined constants, and $\alpha_{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}))=\mu_{ij}+S_{i}$

where the random effects $S_{i}$ are now defined at household level and have mean zero. Warning: this modelling option is available only for the binomial model.

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 [email protected]

Peter J. Diggle [email protected]

See Also

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

---

Langevin-Hastings MCMC for conditional simulation (low-rank approximation)

Description

This function simulates from the conditional distribution of the random effects of binomial and Poisson models.

Usage

Laplace.sampling.lr(
  mu,
  sigma2,
  K,
  y,
  units.m,
  control.mcmc,
  messages = TRUE,
  plot.correlogram = TRUE,
  poisson.llik = FALSE
)

Arguments

mu	mean vector of the linear predictor.
sigma2	variance of the random effect.
K	random effect design matrix, or kernel matrix for the low-rank approximation.
y	vector of binomial/Poisson observations.
units.m	vector of binomial denominators, or offset if the Poisson model is used.
control.mcmc	output from control.mcmc.MCML.
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.
poisson.llik	logical; if poisson.llik=TRUE a Poisson model is used or, if poisson.llik=FALSE, a binomial model is used.

Details

Binomial model. Conditionally on $Z$, the data y follow a binomial distribution with probability $p$ and binomial denominators units.m. Let $K$ denote the random effects design matrix; a logistic link function is used, thus the linear predictor assumes the form

$\log(p/(1-p))=\mu + KZ$

where $\mu$ is the mean vector component defined through mu. Poisson model. Conditionally on $Z$, the data y follow a Poisson distribution with mean $m\lambda$, where $m$ is an offset set through the argument units.m. Let $K$ denote the random effects design matrix; a log link function is used, thus the linear predictor assumes the form

$\log(\lambda)=\mu + KZ$

where $\mu$ is the mean vector component defined through mu. The random effect $Z$ has iid components distributed as zero-mean Gaussian variables with variance sigma2.

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

$\tilde{z}_{curr}+(h/2)\nabla \log f(\tilde{Z} | y),$

where $\tilde{z}_{curr}$ is the current value for $\tilde{Z}$, $h$ is a tuning parameter and $\nabla \log f(\tilde{Z} | y)$ is the the gradient of the log-density of the distribution of $\tilde{Z}$ 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}}(\alpha_{i}-0.547),$

where $c_{1} > 0$ and $0 < c_{2} < 1$ are pre-defined constants, and $\alpha_{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.

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 [email protected]

Peter J. Diggle [email protected]

See Also

control.mcmc.MCML.

---

Independence sampler for conditional simulation of a Gaussian process using SPDE

Description

This function simulates from the conditional distribution of a Gaussian process given binomial y. The Guassian process is also approximated using SPDE.

Usage

Laplace.sampling.SPDE(
  mu,
  sigma2,
  phi,
  kappa,
  y,
  units.m,
  coords,
  mesh,
  control.mcmc,
  messages = TRUE,
  plot.correlogram = TRUE,
  poisson.llik
)

Arguments

mu	mean vector of the Gaussian process to approximate.
sigma2	variance of the Gaussian process to approximate.
phi	scale parameter of the Matern function for the Gaussian process to approximate.
kappa	smothness parameter of the Matern function for the Gaussian process to approximate.
y	vector of binomial observations.
units.m	vector of binomial denominators.
coords	matrix of two columns corresponding to the spatial coordinates.
mesh	mesh object set through inla.mesh.2d.
control.mcmc	control parameters of the Independence sampler set through control.mcmc.MCML.
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.
poisson.llik	logical: if poisson.llik=TRUE then conditional conditional distribution of the data is Poisson; poisson.llik=FALSE then conditional conditional distribution of the data is Binomial.

Details

Binomial model. 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.

Value

A list with the following components

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

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

See Also

control.mcmc.MCML.

---

Heavy metal biomonitoring in Galicia

Description

This data-set contains counts of reported onchocerciasis (or riverblindess) cases from 91 villages sampled across Liberia. The variables are as follows:

• lat latitude of the of sampled villages.

• long longitude of the sampled villages.

• ntest number of tested people for the presence nodules.

• npos number of people that tested positive for the presence of nodules.

• elevation the elevation in meters of the sampled village.

• log_elevation the log-transformed elevation in meters of the sampled village.

Usage

data(liberia)

Format

A data frame with 90 rows and 6 variables

Source

Zouré, H. G. M., Noma, M., Tekle, Afework, H., Amazigo, U. V., Diggle, P. J., Giorgi, E., and Remme, J. H. F. (2014). The Geographic Distribution of Onchocerciasis in the 20 Participating Countries of the African Programme for Onchocerciasis Control: (2) Pre-Control Endemicity Levels and Estimated Number Infected. Parasites & Vectors, 7, 326.

---

Bayesian estimation for the geostatistical linear Gaussian model

Description

This function performs Bayesian estimation for the geostatistical linear Gaussian model.

Usage

linear.model.Bayes(
  formula,
  coords,
  data,
  kappa,
  control.mcmc,
  control.prior,
  low.rank = FALSE,
  knots = NULL,
  messages = TRUE
)

Arguments

formula	an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted.
coords	an object of class formula indicating the geographic coordinates.
data	a data frame containing the variables in the model.
kappa	shape parameter of the Matern covariance function.
control.mcmc	output from control.mcmc.Bayes.
control.prior	output from control.prior.
low.rank	logical; if low.rank=TRUE a low-rank approximation is fitted.
knots	if low.rank=TRUE, knots is a matrix of spatial knots used in the low-rank approximation. Default is knots=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.

Details

This function performs Bayesian estimation for the geostatistical linear Gaussian model, specified as

$Y = d'\beta + S(x) + Z,$

where $Y$ is the measured outcome, $d$ is a vector of coavariates, $\beta$ is a vector of regression coefficients, $S(x)$ is a stationary Gaussian spatial process and $Z$ are independent zero-mean Gaussian variables with variance tau2. More specifically, $S(x)$ has an isotropic Matern covariance function with variance sigma2, scale parameter phi and shape parameter kappa. The shape parameter kappa is treated as fixed.

Priors definition. Priors can be defined through the function control.prior. The hierarchical structure of the priors is the following. Let $\theta$ be the vector of the covariance parameters $(\sigma^2,\phi,\tau^2)$; then each component of $\theta$ can have independent priors freely defined by the user. However, uniform and log-normal priors are also available as default priors for each of the covariance parameters. To remove the nugget effect $Z$, no prior should be defined for tau2. Conditionally on sigma2, the vector of regression coefficients beta has a multivariate Gaussian prior with mean beta.mean and covariance matrix sigma2*beta.covar, while in the low-rank approximation the covariance matrix is simply beta.covar.

Updating the covariance parameters using a Metropolis-Hastings algorithm. In the MCMC algorithm implemented in linear.model.Bayes, the transformed parameters

$(\theta_{1}, \theta_{2}, \theta_{3})=(\log(\sigma^2)/2,\log(\sigma^2/\phi^{2 \kappa}), \log(\tau^2))$

are independently updated using a Metropolis Hastings algorithm. At the $i$-th iteration, a new value is proposed for each from a univariate Gaussian distrubion with variance, say $h_{i}^2$, tuned according the following adaptive scheme

$h_{i} = h_{i-1}+c_{1}i^{-c_{2}}(\alpha_{i}-0.45),$

where $\alpha_{i}$ is the acceptance rate at the $i$-th iteration (0.45 is the optimal acceptance rate for a univariate Gaussian distribution) whilst $c_{1} > 0$ and $0 < c_{2} < 1$ are pre-defined constants. The starting values $h_{1}$ for each of the parameters $\theta_{1}$, $\theta_{2}$ and $\theta_{3}$ can be set using the function control.mcmc.Bayes through the arguments h.theta1, h.theta2 and h.theta3. To define values for $c_{1}$ and $c_{2}$, see the documentation of control.mcmc.Bayes.

Low-rank approximation. In the case of very large spatial data-sets, a low-rank approximation of the Gaussian spatial process $S(x)$ might be computationally beneficial. Let $(x_{1},\dots,x_{m})$ and $(t_{1},\dots,t_{m})$ denote the set of sampling locations and a grid of spatial knots covering the area of interest, respectively. Then $S(x)$ is approximated as $\sum_{i=1}^m K(\|x-t_{i}\|; \phi, \kappa)U_{i}$, where $U_{i}$ are zero-mean mutually independent Gaussian variables with variance sigma2 and $K(.;\phi, \kappa)$ is the isotropic Matern kernel (see matern.kernel). Since the resulting approximation is no longer a stationary process (but only approximately), sigma2 may take very different values from the actual variance of the Gaussian process to approximate. The function adjust.sigma2 can then be used to (approximately) explore the range for sigma2. For example if the variance of the Gaussian process is 0.5, then an approximate value for sigma2 is 0.5/const.sigma2, where const.sigma2 is the value obtained with adjust.sigma2.

Value

An object of class "Bayes.PrevMap". The function summary.Bayes.PrevMap is used to print a summary of the fitted model. The object is a list with the following components:

estimate: matrix of the posterior samples for each of the model parameters.

S: matrix of the posterior samplesfor each component of the random effect. This is only returned for the low-rank approximation.

y: response variable.

D: matrix of covariarates.

coords: matrix of the observed sampling locations.

kappa: vaues of the shape parameter of the Matern function.

knots: matrix of spatial knots used in the low-rank approximation.

const.sigma2: vector of the values of the multiplicative factor used to adjust the sigma2 in the low-rank approximation.

h1: vector of values taken by the tuning parameter h.theta1 at each iteration.

h2: vector of values taken by the tuning parameter h.theta2 at each iteration.

h3: vector of values taken by the tuning parameter h.theta3 at each iteration.

call: the matched call.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

References

Diggle, P.J., Giorgi, E. (2019). Model-based Geostatistics for Global Public Health. CRC/Chapman & Hall.

Giorgi, E., Diggle, P.J. (2017). PrevMap: an R package for prevalence mapping. Journal of Statistical Software. 78(8), 1-29. doi: 10.18637/jss.v078.i08

Higdon, D. (1998). A process-convolution approach to modeling temperatures in the North Atlantic Ocean. Environmental and Ecological Statistics 5, 173-190.

See Also

control.prior, control.mcmc.Bayes, shape.matern, summary.Bayes.PrevMap, autocor.plot, trace.plot, dens.plot, matern, matern.kernel, adjust.sigma2.

---

Maximum Likelihood estimation for the geostatistical linear Gaussian model

Description

This function performs maximum likelihood estimation for the geostatistical linear Gaussian Model.

Usage

linear.model.MLE(
  formula,
  coords = NULL,
  data,
  ID.coords = NULL,
  kappa,
  fixed.rel.nugget = NULL,
  start.cov.pars,
  method = "BFGS",
  low.rank = FALSE,
  knots = NULL,
  messages = TRUE,
  profile.llik = FALSE,
  SPDE = FALSE,
  mesh = NULL,
  SPDE.analytic.hessian = FALSE
)

Arguments

formula	an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted.
coords	an object of class formula indicating the geographic coordinates.
data	a data frame containing the variables in the model.
ID.coords	vector of ID values for the unique set of spatial coordinates obtained from create.ID.coords. These must be provided in order to define a geostatistical model where locations have multiple observations. Default is ID.coords=NULL. See the Details section for more information.
kappa	shape parameter of the Matern covariance function.
fixed.rel.nugget	fixed value for the relative variance of the nugget effect; default is fixed.rel.nugget=NULL if this should be included in the estimation.
start.cov.pars	if ID.coords=NULL, a vector of length two with elements corresponding to the starting values of phi and the relative variance of the nugget effect nu2, respectively, that are used in the optimization algorithm; if ID.coords is provided, a third starting value for the relative variance of the individual unexplained variation nu2.star = omega2/sigma2 must be provided. If nu2 is fixed through fixed.rel.nugget, then start.cov.pars represents the starting value for phi only, if ID.coords=NULL, or for phi and nu2.star, otherwise.
method	method of optimization. If method="BFGS" then the maxBFGS function is used; otherwise method="nlminb" to use the nlminb function. Default is method="BFGS".
low.rank	logical; if low.rank=TRUE a low-rank approximation of the Gaussian spatial process is used when fitting the model. Default is low.rank=FALSE.
knots	if low.rank=TRUE, knots is a matrix of spatial knots that are used in the low-rank approximation. Default is knots=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.
profile.llik	logical; if profile.llik=TRUE the maximization of the profile likelihood is carried out. If profile.llik=FALSE the full-likelihood is used. Default is profile.llik=FALSE.
SPDE	logical; if SPDE=TRUE the SPDE approximation for the Gaussian spatial model is used. Default is SPDE=FALSE.
mesh	an object obtained as result of a call to the function inla.mesh.2d.
SPDE.analytic.hessian	logical; if SPDE.analytic.hessian=TRUE computation of the hessian matrix using the SPDE approximation is carried out using analytical expressions, otherwise a numerical approximation is used. Defauls is SPDE.analytic.hessian=FALSE.

Details

This function estimates the parameters of a geostatistical linear Gaussian model, specified as

$Y = d'\beta + S(x) + Z,$

where $Y$ is the measured outcome, $d$ is a vector of coavariates, $\beta$ is a vector of regression coefficients, $S(x)$ is a stationary Gaussian spatial process and $Z$ are independent zero-mean Gaussian variables with variance tau2. More specifically, $S(x)$ has an isotropic Matern covariance function with variance sigma2, scale parameter phi and shape parameter kappa. In the estimation, the shape parameter kappa is treated as fixed. The relative variance of the nugget effect, nu2=tau2/sigma2, can be fixed though the argument fixed.rel.nugget; if fixed.rel.nugget=NULL, then the variance of the nugget effect is also included in the estimation.

Locations with multiple observations. If multiple observations are available at any of the sampled locations the above model is modified as follows. Let $Y_{ij}$ denote the random variable associated to the measured outcome for the j-th individual at location $x_{i}$. The linear geostatistical model assumes the form

$Y_{ij} = d_{ij}'\beta + S(x_{i}) + Z{i} + U_{ij},$

where $S(x_{i})$ and $Z_{i}$ are specified as mentioned above, and $U_{ij}$ are i.i.d. zer0-mean Gaussian variable with variance $\omega^2$. his model can be fitted by specifing a vector of ID for the unique set locations thourgh the argument ID.coords (see also create.ID.coords).

Low-rank approximation. In the case of very large spatial data-sets, a low-rank approximation of the Gaussian spatial process $S(x)$ can be computationally beneficial. Let $(x_{1},\dots,x_{m})$ and $(t_{1},\dots,t_{m})$ denote the set of sampling locations and a grid of spatial knots covering the area of interest, respectively. Then $S(x)$ is approximated as $\sum_{i=1}^m K(\|x-t_{i}\|; \phi, \kappa)U_{i}$, where $U_{i}$ are zero-mean mutually independent Gaussian variables with variance sigma2 and $K(.;\phi, \kappa)$ is the isotropic Matern kernel (see matern.kernel). Since the resulting approximation is no longer a stationary process, the parameter sigma2 is adjusted by a factorconstant.sigma2. See adjust.sigma2 for more details on the the computation of the adjustment factor constant.sigma2 in the low-rank approximation.

Value

An object of class "PrevMap". The function summary.PrevMap is used to print a summary of the fitted model. The object is a list with the following components:

estimate: estimates of the model parameters; use the function coef.PrevMap to obtain estimates of covariance parameters on the original scale.

covariance: covariance matrix of the ML estimates.

log.lik: maximum value of the log-likelihood.

y: response variable.

D: matrix of covariates.

coords: matrix of the observed sampling locations.

ID.coords: set of ID values defined through the argument ID.coords.

method: method of optimization used.

kappa: fixed value of the shape parameter of the Matern function.

knots: matrix of the spatial knots used in the low-rank approximation.

const.sigma2: adjustment factor for sigma2 in the low-rank approximation.

fixed.rel.nugget: fixed value for the relative variance of the nugget effect.

mesh: the mesh used in the SPDE approximation.

call: the matched call.

Author(s)

Emanuele Giorgi [email protected]

Peter J. Diggle [email protected]

References

Diggle, P.J., Giorgi, E. (2019). Model-based Geostatistics for Global Public Health. CRC/Chapman & Hall.

Giorgi, E., Diggle, P.J. (2017). PrevMap: an R package for prevalence mapping. Journal of Statistical Software. 78(8), 1-29. doi: 10.18637/jss.v078.i08

Higdon, D. (1998). A process-convolution approach to modeling temperatures in the North Atlantic Ocean. Environmental and Ecological Statistics 5, 173-190.

See Also

shape.matern, summary.PrevMap, coef.PrevMap, matern, matern.kernel, maxBFGS, nlminb.

---