likfit: 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", realisations,
       method.lik = "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

coords

an $n \times 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. See documentation of trend.spatial for further details. 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

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

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

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

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.

realisations

optional. A vector indicating replication number for each data. For more details see as.geodata.

method.lik

options are "ML" for maximum likelihood and "REML" for restricted maximum likelihood. Defaults to "ML".

components

an $n \times 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 fu

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.modela string with the name of the correlation function.
nuggetvalue of the nugget parameter $\tau^2$. This is an estimate if fix.nugget = FALSE otherwise, a fixed value.
cov.parsa vector with the estimates of the parameters $\sigma^2$ and $\phi$, respectively.
kappavalue of the smoothness parameter. Valid only if the correlation function is one of: "matern", "powered.exponential", "cauchy" or "gneiting.matern".
betaestimate of mean parameter $\beta$. This can be a scalar or vector depending on the trend (covariates) specified in the model.
beta.varestimated variance (or covariance matrix) for the mean parameter $\beta$.
lambdavalues of the Box-Cox transformation parameter. A fixed value if fix.lambda = TRUE otherwise the estimate value.
aniso.parsfixed values or estimates of the anisotropy parameters, according to the function call.
method.likestimation method used, "ML" (maximum likelihood) or "REML" (restricted maximum likelihood).
loglikthe value of the maximized likelihood.
nparsnumber of estimated parameters.
AICvalue of the Akaike information criteria.
BICvalue of the Bayesian information criteria.
parameters.summarya data-frame with all model parameters, their status (estimated or fixed) and values.
info.minimisationresults returned by the minimisation function.
max.distmaximum distance between 2 data points. This information relevant for other functions which use outputs from likfit.
trend.matrixthe trend (covariates) matrix $X$.
log.jacobiannumerical value of the logarithm of the Jacobian of the transformation.
nospatialestimates for the model without the spatial component.
callthe function call.

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)) = \frac{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 Rfunction 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 \time 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.

References

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

Examples

Run this code

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)

<testonly>ap <- grf(50, cov.pars=c(1, .3), nug=.3)
ml <- likfit(ap, ini=c(0.5, 0.5), nug=0.2)
ml <- likfit(ap, ini=c(0.5, 0.5), fix.nug=TRUE, nug=0.2)
ml <- likfit(ap, data=exp(ap$data), ini=c(0.5, 0.5), nug=0.2, fix.lambda=FALSE)
ml <- likfit(ap, ini=c(0.5, 0.5), nug=0.2, fix.psiR = FALSE)
ml <- likfit(ap, ini=c(0.5, 0.5), nug=0.2, fix.psiA = FALSE, fix.psiR = FALSE)
ml <- likfit(ap, ini=c(0.5, 0.5), cov.model="matern", fix.kappa=TRUE, kappa=2)
ml <- likfit(ap, ini=c(0.5, 0.5), cov.model="matern", fix.kappa=FALSE)
ml <- likfit(ap, ini=expand.grid(c(.5,1,1.5),c(.1,.2,.3)))
ml <- likfit(ap, ini=expand.grid(c(.5,1),c(.1,.2)), nug=c(.2,.3))

## Multiple realizations
ap1 <- grf(30, cov.pars=c(1, .3))
ap2 <- grf(20, cov.pars=c(1, .3))
ap3 <- grf(40, cov.pars=c(1, .3))
ap <- list(coords = rbind(ap1$coords, ap2$coords, ap3$coords), data = c(ap1$data, ap2$data, ap3$data))
ap$realisations <- c(rep(1,30), rep(2,20), rep(3,40))
ap.fit <- likfit(ap, ini=c(.5, .5), reali=ap$real)</testonly>

Run the code above in your browser using DataLab