R: Likelihood Based Parameter Estimation for Gaussian Random Fields

likfit {geoR}

R Documentation

Likelihood Based Parameter Estimation for Gaussian Random Fields

Description

Maximum likelihood (ML) or restricted maximum likelihood (REML) parameter estimation for (transformed) Gaussian random fields.

Usage

likfit(geodata, coords = geodata$coords, data = geodata$data,
       trend = "cte", ini.cov.pars, fix.nugget = FALSE, nugget = 0,
       fix.kappa = TRUE, kappa = 0.5, fix.lambda = TRUE, lambda = 1,
       fix.psiA = TRUE, psiA = 0, fix.psiR = TRUE, psiR = 1, 
       cov.model = "matern", method = "ML", components = FALSE,
       nospatial = TRUE, limits = likfit.limits(),
       print.pars = FALSE, messages.screen = TRUE, ...)

Arguments

`geodata`	a list containing elements `coords` and `data` as described next. Typically an object of the class `"geodata"` - a geoR data-set. If not provided the arguments `coords` and `data` must be provided instead.
`coords`	an n x 2 matrix where each row has the 2-D coordinates of the n data locations. By default it takes the component `coords` of the argument `geodata`, if provided.
`data`	a vector with n data values. By default it takes the component `data` of the argument `geodata`, if provided.
`trend`	specifies the mean part of the model. The options are: `"cte"` (constant mean), `"1st"` (a first degree polynomial on the coordinates), `"2nd"` (a second degree polynomial on the coordinates), or a formula of the type `~X` where `X` is a matrix with the covariates (external trend). Defaults to `"cte"`.
`ini.cov.pars`	initial values for the covariance parameters: sigma^2 (partial sill) and phi (range parameter). Typically a vector with two components. However a matrix can be used to provide several initial values. See `DETAILS` below.
`fix.nugget`	logical, indicating whether the parameter tau^2 (nugget variance) should be regarded as fixed (`fix.nugget = TRUE`) or should be estimated (`fix.nugget = FALSE`). Defaults to `FALSE`.
`nugget`	value of the nugget parameter. Regarded as a fixed value if `fix.nugget = TRUE` otherwise as the initial value for the minimization algorithm. Defaults to zero.
`fix.kappa`	logical, indicating whether the extra parameter kappa should be regarded as fixed (`fix.kappa = TRUE`) or should be estimated (`fix.kappa = FALSE`). Defaults to `TRUE`.
`kappa`	value of the extra parameter kappa. Regarded as a fixed value if `fix.kappa = TRUE` otherwise as the initial value for the minimization algorithm. Defaults to 0.5. This parameter is valid only if the covariance function is one of: `"matern"`, `"powered.exponential"`, `"cauchy"` or `"gneiting.matern"`. For more details on covariance functions see documentation for `cov.spatial`.
`fix.lambda`	logical, indicating whether the Box-Cox transformation parameter lambda should be regarded as fixed (`fix.lambda = TRUE`) or should be be estimated (`fix.lambda = FALSE`). Defaults to `TRUE`.
`lambda`	value of the Box-Cox transformation parameter lambda. Regarded as a fixed value if `fix.lambda = TRUE` otherwise as the initial value for the minimization algorithm. Defaults to 1. Two particular cases are lambda = 1 indicating no transformation and lambda = 0 indicating log-transformation.
`fix.psiA`	logical, indicating whether the anisotropy angle parameter psi_R should be regarded as fixed (`fix.psiA = TRUE`) or should be estimated (`fix.psiA = FALSE`). Defaults to `TRUE`.
`psiA`	value (in radians) for the anisotropy angle parameter psi_A. Regarded as a fixed value if `fix.psiA = TRUE` otherwise as the initial value for the minimization algorithm. Defaults to 0. See `coords.aniso` for further details on anisotropy correction.
`fix.psiR`	logical, indicating whether the anisotropy ratio parameter psi_R should be regarded as fixed (`fix.psiR = TRUE`) or should be estimated (`fix.psiR = FALSE`). Defaults to `TRUE`.
`psiR`	value, always greater than 1, for the anisotropy ratio parameter psi_R. Regarded as a fixed value if `fix.psiR = TRUE` otherwise as the initial value for the minimization algorithm. Defaults to 1. See `coords.aniso` for further details on anisotropy correction.
`cov.model`	a string specifying the model for the correlation function. For further details see documentation for `cov.spatial`. Defaults are equivalent to the exponential model.
`method`	options are `"ML"` for maximum likelihood and `"REML"` for restricted maximum likelihood. Defaults to `"ML"`.
`components`	an n x 3 data-frame with fitted values for the three model components: trend, spatial and residuals. See the section `DETAILS` below for the model specification.
`nospatial`	logical. If `TRUE` parameter estimates for the model without spatial component are included in the output.
`limits`	values defining lower and upper limits for the model parameters used in the numerical minimization. The auxiliary function `likfit.limits()` is called to set the limits.
`print.pars`	logical. If `TRUE` the parameters and the value of the negative log-likelihood (up to a constant) are printed each time the function to be minimised is called.
`messages.screen`	logical. Indicates whether status messages should be printed on the screen (or output device) while the function is running.
`...`	additional parameters to be passed to the minimization function. Typically arguments of the type `control()` which controls the behavior of the minimization algorithm. For further details see documentation for the minimization function `optim`.

Details

This function estimate the parameters of the Gaussian random field model, specified here by:

Y(x) = mu(x) + S(x) + e

where

x defines a spatial location. Typically Euclidean coordinates on a plane.
Y is the variable been observed.
mu(x) = X %*% beta is the mean component of the model (trend).
S(x) is a stationary Gaussian process with variance sigma^2 (partial sill) and a correlation function parametrized by phi (the range parameter). Possible extra parameters for the correlation function are the smoothness parameter kappa and the anisotropy parameters phi_R and phi_A (anisotropy ratio and angle, respectively).
e is the error term with variance parameter tau^2 (nugget variance).

The additional parameter lambda allows the Box-Cox transformation. If used Y(x) above is replaced by g(Y(x)) such that

g(Y(x)) = ((Y^lambda(x)) - 1)/lambda .

Two cases of particular interest are lambda = 1 indicating no transformation and lambda = 0 indicating log-transformation.

Parameter estimation is performed numerically using the R function optim to minimize the negative log-likelihood computed by negloglik.GRF.

Lower and upper limits for parameter values can be individually specified using the function likfit.limits(). For example, including the following in the function call:
limits = likfit.limits(phi=c(0, 10), lambda=c(-2.5, 2.5)),
will change the limits for the parameters phi and lambda. Default values are used if the argument limits is not provided.

If the fix.lambda = FALSE and nospatial = FALSE the Box-Cox parameter for the model without the spatial component is obtained numerically, with log-likelihood computed by the function boxcox.ns.

Multiple initial values can be specified providing a n x 2 matrix for the argument ini.cov.pars and/or providing a vector for the values of the remaining model parameters. In this case the log-likelihood is computed for all combinations of model parameters. The set with results in the maximum value of the log-likelihood is then used to start the minimisation algorithm.

Value

An object of the classes "likGRF" and "variomodel".
The function summary.likGRF is used to print a summary of the fitted model.
The object is a list with the following components:

`cov.model`	a string with the name of the correlation function.
`nugget`	value of the nugget parameter tau^2. This is an estimate if `fix.nugget = FALSE` otherwise, a fixed value.
`cov.pars`	a vector with the estimates of the parameters sigma^2 and phi, respectively.
`kappa`	value of the smoothness parameter. Valid only if the correlation function is one of: `"matern"`, `"powered.exponential"`, `"cauchy"` or `"gneiting.matern"`.
`beta`	estimate of mean parameter beta. This can be a scalar or vector depending on the trend (covariates) specified in the model.
`beta.var`	estimated variance (or covariance matrix) for the mean parameter beta.
`lambda`	values of the Box-Cox transformation parameter. A fixed value if `fix.lambda = TRUE` otherwise the estimate value.
`aniso.pars`	fixed values or estimates of the anisotropy parameters, according to the function call.
`method`	estimation method used, `"ML"` (maximum likelihood) or `"REML"` (restricted maximum likelihood).
`loglik`	the value of the maximized likelihood.
`npars`	number of estimated parameters.
`AIC`	value of the Akaike information criteria.
`BIC`	value of the Bayesian information criteria.
`parameters.summary`	a data-frame with all model parameters, their status (estimated or fixed) and values.
`info.minimisation`	results returned by the minimisation function.
`max.dist`	maximum distance between 2 data points. This information relevant for other functions which use outputs from `likfit`.
`trend.matrix`	the trend (covariates) matrix X.
`log.jacobian`	numerical value of the logarithm of the Jacobian of the transformation.
`nospatial`	estimates for the model without the spatial component.
`call`	the function call.

Author(s)

Paulo Justiniano Ribeiro Jr. Paulo.Ribeiro@est.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

Further information about geoR can be found at:
http://www.maths.lancs.ac.uk/~ribeiro/geoR.html.

Examples

if(is.R()) data(s100)
ml <- likfit(s100, ini=c(0.5, 0.5), fix.nug = TRUE)
ml
summary(ml)
reml <- likfit(s100, ini=c(0.5, 0.5), fix.nug = TRUE, met = "REML")
summary(reml)
plot(variog(s100))
lines(ml)
lines(reml, lty = 2)