Maximum Likelihood estimation for the geostatistical linear Gaussian model

Share:

Description

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

Usage

1
2
3
4
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'β + S(x) + Z,

where Y is the measured outcome, d is a vector of coavariates, β 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}'β + 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 ω^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},…,x_{m}) and (t_{1},…,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 ∑_{i=1}^m K(\|x-t_{i}\|; φ, κ)U_{i}, where U_{i} are zero-mean mutually independent Gaussian variables with variance sigma2 and K(.;φ, κ) 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 e.giorgi@lancaster.ac.uk

Peter J. Diggle p.diggle@lancaster.ac.uk

References

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.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
data(loaloa)
# Empirical logit transformation
loaloa$logit <- log((loaloa$NO_INF+0.5)/(loaloa$NO_EXAM-loaloa$NO_INF+0.5))
fit.MLE <- linear.model.MLE(logit ~ 1,coords=~LONGITUDE+LATITUDE,
                data=loaloa, start.cov.pars=c(0.2,0.15),
                 kappa=0.5)
summary(fit.MLE)

# Low-rank approximation
data(data_sim)
n.subset <- 200
data_subset <- data_sim[sample(1:nrow(data_sim),n.subset),]

# Logit transformation
data_subset$logit <- log(data_subset$y+0.5)/
                     (data_subset$units.m-
                      data_subset$y+0.5)
knots <- as.matrix(expand.grid(seq(-0.2,1.2,length=8),seq(-0.2,1.2,length=8)))

fit <- linear.model.MLE(formula=logit~1,coords=~x1+x2,data=data_subset,
                             kappa=2,start.cov.pars=c(0.15,0.1),low.rank=TRUE,
                             knots=knots)
summary(fit,log.cov.pars=FALSE)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.