spagmix (version 0.3-5)

lgcpmix: Generate a spatial log-Gaussian Cox process intensity

Description

Generate a realisation of a (possibly inhomogeneous) log-Gaussian Cox process (LGCP) spatial intensity function with an identifiable mean structure.

Usage

lgcpmix(lambda, covmodel = "exp", covpars = NULL)

Arguments

lambda

A pixel image giving the deterministic spatial intensity as the mean structure of the process. The generated Gaussian field will match the dimensions, resolution and domain of this object.

covmodel

A character string giving the short name of a spatial covariance model available in the RandomFields package (Schlather et al., 2015). The code searches for a function of this name after adding the prefix "RM", piggybacking on the internals of rLGCP. See `Details'.

covpars

A named list of values for the parameters required by the chosen covmodel. This will normally include var and scale. Both variance and scale parameters default to 1 if this argument is left NULL.

Value

A pixel image giving the generated intensity function, comprised of the product of lambda (fixed, and unchanging in repeated calls to this function) and the exponentiated Gaussian field (with expected value 1, this is stochastic and varies in repeated calls).

Details

This function allows the user to generate a spatial intensity function \(\Gamma\) of the form

\(\Gamma(x) = \lambda(x)\exp[Y(x)]\)

for \(x \in W\), where \(\lambda(x)\) (passed to lambda) is the deterministic spatial intensity over the spatial domain \(W\), and \(Y(x)\) is a Gaussian random field on \(W\). This Gaussian field is defined with a particular spatial covariance function (specified in covmodel) with variance and scale parameters \(\sigma^2\) and \(\phi\) respectively, as well as any additionally required parameter values (all specified in covpars).

The mean parameter \(\mu\) of the Gaussian field \(Y\) is internally fixed at \(-\sigma^2/2\); negative half the variance. This is for identifiability of the mean structure, forcing \(E[Y(x)] = 1\) for all \(x \in W\) (see theoretical properties in M<U+00F8>ller et al., 1998). This means the deterministic intensity function \(\lambda(x)\) is solely responsible for describing fixed heterogeneity in spatial intensity over \(W\), with the randomly generated Gaussian field left to describe residual stochastic spatial correlation. This presents a highly flexible class of model, even with stationarity and isotropy of the Gaussian field itself, and is intuitively sensible in a variety of applications. See Diggle et al. (2005) and Davies & Hazelton (2013) for example.

As such, the pixel image supplied to lambda as \(\lambda(x)\) must be non-negatively-valued and yield a finite integral. The choice of covariance model and correspondingly required parameters as well as actual simulation of the Gaussian field is deferred to functionality in the RandomFields package; see RMmodel for possible choices. For example, requesting covmodel = "exp" (default) will search for the RandomFields function RMexp and imposes an exponential covariance structure on the generated field whereby \(Cov(u) = \sigma^2\exp(-u/\phi)\) for the Euclidean distance between any two spatial locations \(u\).

To generate a subsequent dataset, use e.g. rpoispp or rpoispoly.

References

Davies, T.M. and Hazelton, M.L. (2013), Assessing minimum contrast parameter estimation for spatial and spatiotemporal log-Gaussian Cox processes, Statistica Neerlandica, 67(4) 355--389.

Diggle, P.J., Rowlingson, B. and Su, T. (2005), Point process methodology for on-line spatio-temporal disease surveillance, Environmetrics, 16 423--434.

M<U+00F8>ller, J., Syversveen, A.R. and Waagepetersen, R.P. (1998), Log-Gaussian Cox processes, Scandinavian Journal of Statistics, 25(3) 451--482.

Schlather, M., Malinowski, A., Menck, P.J., Oesting, M. and Strokorb, K. (2015) Analysis, Simulation and Prediction of Multivariate Random Fields with Package RandomFields, Journal of Statistical Software, 63(8) 1--25.

See Also

rLGCP, rpoispp, rpoispoly

Examples

# NOT RUN {
## Homogeneous example ##

# Create constant intensity image integrating to 500

homog <- as.im(as.mask(toywin))
homog <- homog/integral(homog)*500


# Corresponding LGCP realisations using exponential covariance structure
par(mfrow=c(2,2),mar=rep(1.5,4))
for(i in 1:4){
  temp <- lgcpmix(homog,covmod="exp",covpars=list(var=1,scale=0.2))
  plot(temp,main=paste("Realisation",i),log=TRUE)
}



## Inhomogeneous examples ##

# Create deterministic trend

mn <- cbind(c(0.25,0.8),c(0.31,0.82),c(0.43,0.64),c(0.63,0.62),c(0.21,0.26))
v1 <- matrix(c(0.0023,-0.0009,-0.0009,0.002),2)
v2 <- matrix(c(0.0016,0.0015,0.0015,0.004),2)
v3 <- matrix(c(0.0007,0.0004,0.0004,0.0011),2)
v4 <- matrix(c(0.0023,-0.0041,-0.0041,0.0099),2)
v5 <- matrix(c(0.0013,0.0011,0.0011,0.0014),2)
vr <- array(NA,dim=c(2,2,5))
for(i in 1:5) vr[,,i] <- get(paste("v",i,sep=""))
intens <- sgmix(mean=mn,vcv=vr,window=toywin,p0=0.1,int=500)


# Two realisations (identical calls to function), exponential covariance structure

r1exp <- lgcpmix(lambda=intens,covmodel="exp",covpars=list(var=2,scale=0.05))
r2exp <- lgcpmix(lambda=intens,covmodel="exp",covpars=list(var=2,scale=0.05))


# Two more realisations, Matern covariance with smoothness 1

r1mat <- lgcpmix(lambda=intens,covmodel="matern",covpars=list(var=2,scale=0.05,nu=1))
r2mat <- lgcpmix(lambda=intens,covmodel="matern",covpars=list(var=2,scale=0.05,nu=1))


# Plot everything, including 'intens' alone (no correlation)

par(mar=rep(2,4))
layout(matrix(c(1,2,4,1,3,5),3))
plot(intens,main="intens alone",log=TRUE)
plot(r1exp,main="realisation 1\nexponential covar",log=TRUE)
plot(r2exp,main="realisation 2\nexponential covar",log=TRUE)
plot(r1mat,main="realisation 1\nMatern covar",log=TRUE)
plot(r2mat,main="realisation 2\nMatern covar",log=TRUE)


# Plot example datasets
dint <- rpoispoly(intens,w=toywin)
d1exp <- rpoispoly(r1exp,w=toywin)
d2exp <- rpoispoly(r2exp,w=toywin)
d1mat <- rpoispoly(r1mat,w=toywin)
d2mat <- rpoispoly(r2mat,w=toywin)

par(mar=rep(2,4))
layout(matrix(c(1,2,4,1,3,5),3))
plot(dint,main="intens alone",log=TRUE)
plot(d1exp,main="realisation 1\nexponential covar",log=TRUE)
plot(d2exp,main="realisation 2\nexponential covar",log=TRUE)
plot(d1mat,main="realisation 1\nMatern covar",log=TRUE)
plot(d2mat,main="realisation 2\nMatern covar",log=TRUE)
# }