Estimates key covariance parameters for a spatial process.
Maximizes the likelihood to determine the nugget variance (sigma^2), the sill ( rho) and the range (theta) for a spatial process.
MLESpatialProcess( x, y, weights = rep(1, nrow(x)), Z = NULL,
mKrig.args = NULL,
cov.function = "stationary.cov",
cov.args = list(Covariance = "Matern",
smoothness = 1),
gridTheta = NULL,
gridN = 15,
optim.args = NULL,
na.rm = TRUE,
verbose = FALSE,
abstol = 1e-4,
REML = FALSE,
cov.params.start = NULL,
...)x |
A matrix of spatial locations with rows indexing location and columns the dimension (e.g. longitude/latitude) |
y |
The spatial observations. This can be a matrix of reoplicated spatial data. |
weights |
Precision ( 1/variance) of each observation |
Z |
Linear covariates to be included in fixed part of the
model that are distinct from the default low order
polynomial in |
mKrig.args |
A list containing other objects to pass to mKrig. |
gridN |
Number of points to use in grid search over theta. |
cov.function |
The name of the covariance function (See help on Krig for details. ) |
gridTheta |
A grid of range parameters to search over. The default is to use a range based on the quantiles of pairwise distances of the spatial locations. |
cov.args |
A list with arguments for the covariance functions. These are usually parameters and other options such as the type of distance function. |
optim.args |
Additional arguments passed to the optim function for likelihood
maximization. The default value is:
|
na.rm |
If TRUE remove missing values in y and corresponding locations in x. |
verbose |
If TRUE print out intermediate information for debugging. |
abstol |
Absolute tolerance used to judeg convergence in optim. |
REML |
If TRUE use maximize the restricted Likelihood instead of the concentrated likelihood.(Preliminary experience suggests this does not make much difference.) |
cov.params.start |
A list with each component being the name of a covariance paramater to optimize over and the value being the starting value. This is the same format used by optim. If NULL no additional parameters are optimized and their values are fixed in the cov.args list. |
... |
Additional arguments to pass to the mKrig function. |
MLESpatialProcess is designed to be a simple and easy to use function for
maximizing the likelihood for a Gaussian spatial process. For other fixed,
covariance parameters, the likelihood is maximized over the nugget and sill
parameters using the mKrig function. lambda and theta
are optimized using the mKrigMLEJoint function on a log scale. This is
a wrapper for two lower level functions to make spatialProcess readable. We
recommend using mKrigMLEJoint and mKrigMLEGrid for direct
optimzation and grid searches over parameters and consulting
those two commented functions for an exact description of how
the computations are being done. At the very basic
level the likelihood function is evaluated by calling the
mKrig function with the covariance parameters being
adjusted to test different values. So the lowest level
likelihood evaluate and what one might use at the highest
level are the same workhorse function, mKrig.
Note the likelihood can be maximized analytically over the parameters of the fixed part of the spatial model and with the nugget (sigma) and sill (rho) reduced to the single parameter lambda= sigma^2/rho. The likelihood is maximized numerically over lambda and theta if there are additional covariance parameters ( such as smoothness for the Matern) these need to be fixed and so the MLE is found for the covariance conditional on these additional parameter values. From a practical point of view it is often difficult to estimate just these three from a moderate spatial data set and the user is encourage to try different combinations of fixing covariance parameters with ML for the remaining ones.
MLESpatialProcess.fast is an older fields function also using
the optim
function to maximize the likelihood computed from the
mKrig function. It will
eventually be removed from later versions of fields but perhaps
is
still useful as a cross
check on these newer functions
MLESpatialProcess:
A list (with sublists!) that includes components:
summary A vector with the optimized parameters (aka MLEs) along with log likelihood
at the maximum and information from optim on convergence. Also included are the effecive
degrees of freedom of the model at the MLEs and the GCV function value.
MLEJoint A list with components:
summary, lnLike.eval, optimResults, pars.MLE, parTransform.
summary gives the results from the joint
optimization over
all parameters. lnLike.eval is a table
that has captured all the likelihood evaluations as optim has
searched for the maxmimum. Note that
if the method has converged many of these parameter value
choices will be close the final converged
value but often still useful to investigate the likelihood
surface around the maximum. optimResults is the list
returned by the optim function.
pars.MLE are the parameters specifically optimized over.
par.Transform the function used to
transform the parameters for a more suitable optimization.
The default is a log transformation for
lambda and theta.
MLEGrid A list with results from profiling the likelihood over theta, the range parameter.
Here the component summary is a table giving the values tried for theta and the corresponding parameters that
maximize the likelihood with this fixed value for theta. Also included in the table are some convergence
information, effective degrees of freedom, and the GCV criterion.
MLEProfileLambda Results from profiling the likelihood over lambda. Format is the same as that described for MLEGrid.
call The call made to this function.
timingTable Some timing information for the joint optimzation, as the column labeled timeOptim, and
the two profiling computations with timeGrid being for theta and timeProfile being for lambda.
MLESpatialProcess.fast has been depreciated and is included for backward compatibility.
Doug Nychka, John Paige
#
# Trying out on a simulated data set
# Generate observation locations (50 is really too small but this is
# just to make this run quickly)
n <- 50
set.seed(124)
x = matrix(runif(2*n), nrow=n)
#generate observations at the locations
trueTheta = .1
trueSigma = .01
covMat = exp( -rdist(x,x) /trueTheta )
# As rcode: y = t(covMat)%*% (rnorm(n)) + trueSigma * rnorm( n)
y <- t(covMat) %*% (rnorm(n)) + trueSigma * rnorm( n)
# Use exponential covariance estimate constant function for mean
out <- MLESpatialProcess(x, y,
smoothness = .5,
mKrig.args = list( m = 1),
)
out$summary
## Not run:
#Now a (small) grid search over the smoothness
# add replicated data fields ( i.e. independent copies drawn from same covariance model)
# to make this stable
set.seed(124)
# Y = t(covMat)
Y = t(covMat)%*% matrix(rnorm(n*200), n,200) +
trueSigma * matrix((rnorm(n*200)), n,200)
#This may take a few seconds
testSmoothness = c(.5, 1.5, 2.0)
for( nu in testSmoothness){
out = MLESpatialProcess(x, Y, cov.args = list(Covariance="Matern"),
smoothness = nu, cov.params.start= list( lambda=.5))
print( out$MLEJoint$summary)
}
## End(Not run)
# example with a covariate
## Not run:
data(COmonthlyMet)
ind<- !is.na( CO.tmean.MAM.climate)
x<- CO.loc[ind,]
y<- CO.tmean.MAM.climate[ind]
elev<- CO.elev[ind]
obj2<- MLESpatialProcess( x,y)
obj3<- MLESpatialProcess( x,y, Z=elev)
# elevation makes a difference
obj2$MLEJoint$summary
obj3$MLEJoint$summary
## End(Not run)
## Not run:
# fits for first 10 days from ozone data
data( ozone2)
NDays<- 5
O3MLE<- NULL
for( day in 1: NDays){
cat( day, " ")
ind<- !is.na(ozone2$y[day,] )
x<- ozone2$lon.lat[ind,]
y<- ozone2$y[day,ind]
print( length( y))
out<- MLESpatialProcess( x,y,
Distance="rdist.earth")$MLEJoint$summary
O3MLE<- rbind( O3MLE, out)
}
# NOTE: names from summary:
#[1] "lnProfileLike.FULL" "lambda" "theta"
#[4] "sigmaMLE" "rhoMLE" "funEval"
#[7] "gradEval" "eff.df" "GCV"
plot( log(O3MLE[,"lambda"]), log(O3MLE[,"theta"]))
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.