MLESpatialProcess: Estimates key covariance parameters for a spatial process.

Description

Maximizes the likelihood to determine the nugget variance (sigma^2), the sill ( rho) and the range (theta) for a spatial process.

Usage

MLESpatialProcess(
x, y, theta.grid = NULL, cov.function =
                 "stationary.cov", cov.args = list(Covariance =
                 "Matern", smoothness = 1), ngrid = 10, niter = 15, tol
                 = 0.01, Distance = "rdist", nstep.cv = 50, verbose =
                 FALSE, ...)
                           
MLESpatialProcess.fast(x, y, lambda.start = NULL, theta.start = NULL,
 cov.function = "stationary.cov", cov.args = list(Covariance = "Matern", smoothness = 1),
 relative.tolerance = 0.001, Distance = "rdist", verbose = FALSE, ...)

Arguments

Value

MLESpatialProcess: A list that includes components: theta.MLE, rho.MLE, sigma.MLE, lambda.MLE being the maximum likelihood estimates of these parameters. The component REML.grid is a two column matrix with the first column being the theta grid and the second column being the profiled and restricted likelihood for that value of theta. Here profile means that the likelihood has already been evaluated at the maximum over sigma and rho for this value of theta. eval.grid is a more complete "capture" of the evaluations being a 6 column matrix with the parameters theta, lambda, sigma, rho, profile likelihood and the effective degrees of freedom. This is just last row of lambda.est returned by the core function Krig MLESpatialProcess.fast here the returned value is limited because this function isbuilt around calls to mKrig. Returned value is a list with components: pars, the MLEs for theta, rho, sigma and lambda, logLikelihood,values of the log likelihood at the maximum, eval.grid, a table with the results from evaluating different combinations of parameters,

converge, convergence flag from optim (0=Successfull) and number of evaluations used to find maximum. and call, the calling arguments.

Details

MLESpatialProcess is designed to be a robust but perhaps slow function to maximize the likelihood for a Gaussian spatial process. For fixed range the likelihood is maximzed over the nugget and sill parameters using the Krig function. An outer optimization finds the maximum over the range parameter as well. See the help(Krig) for details of the restricted maximum likelihood criterion (REML). MLESpatialProcess.fast uses the optim function to maximize the likelihood computed from the mKrig function. It is more efficient in the computation as it does not find a full eigen decomposition with each new value of theta and maximizes the likelihood over theta and lambda simultaneously. Note the likelihood can be maximized analytically over the parameters of the fixed part and with the nugget (sigma) and sill (rho) reduced to the single parameter lambda= sigma^2/rho. So fixing any other covariance parameters the likelihood is maximzed numerically over lambda and theta. The differences between these two functions is due to the differences between the definition of the restricted likelihood used in Krig and the conventional likelihood used in mKrig.

Examples

Run this code

N<- 100
  set.seed(123)
  x<- matrix(runif(2*N), N,2)
  theta<- .2
  Sigma<-  Matern( rdist(x,x)/theta , smoothness=1.0)
  Sigma.5<- chol( Sigma)
  sigma<- .1
  F.true<- t( Sigma.5) %*%  rnorm(N)
  Y<-  F.true +  sigma*rnorm(N)
# find MLE for  sigma rho and theta  smoothness fixed  first
# data set
  obj<- MLESpatialProcess( x,Y)
  obj$pars
  # profile likelihood over theta
  plot(obj$eval.grid[,1], obj$eval.grid[,6], xlab="theta",
  ylab= "log Profile likelihood", type="p" )
  xline( obj$pars["theta"], col="red")
# log likelihood surface over theta and  log lambda
image.plot( obj$logLikelihoodSurface$x,
               obj$logLikelihoodSurface$y, obj$logLikelihoodSurface$z,
               xlab="theta (range)", ylab="log lambda" )
# MLE               
  points( obj$pars[1], log(obj$pars[2]), pch=16, col="magenta", cex=1.2)      

# using "fast" version  
obj.fast<- MLESpatialProcess.fast( x,Y)
obj.fast$pars 
# points where likelihood evaluated:
quilt.plot( log( obj.fast$eval.grid[,1:2] ), obj.fast$eval.grid[,7],
                   xlab="log(theta)",ylab="log(lambda)")

# parameters are slightly different due to the differences of REML and the full likelihood.
                 
# example with a covariate  
data(COmonthlyMet)
obj2<-  MLESpatialProcess( CO.loc, CO.tmean.MAM.climate)
obj3<-  MLESpatialProcess( CO.loc, CO.tmean.MAM.climate, Z= CO.elev)
ind<- !is.na( CO.tmean.MAM.climate)
obj4<-  MLESpatialProcess.fast( CO.loc[ind,], CO.tmean.MAM.climate[ind],
               Z= CO.elev[ind])
# elevation makes a difference
obj2$pars
obj3$pars
obj4$pars
 # fits for ozone data
data( ozone2) 	 
NDays<- nrow( ozone2$y)
O3MLE<- matrix( NA, nrow= NDays, ncol=4)
dimnames(O3MLE)<- list(NULL, c("theta",  "lambda", "rho",    "sigma"))
for( day in 1: NDays){
	cat( day, " ")
	O3MLE[day,]<- MLESpatialProcess( ozone2$lon.lat, ozone2$y[day,],
	 Distance="rdist.earth")$pars
}
plot( log(O3MLE[,1]), log(O3MLE[,3]))