mKrigMLE: Maximizes likelihood for the process marginal variance...

In fields: Tools for Spatial Data

Maximizes likelihood for the process marginal variance (sigma) and nugget standard deviation (tau) parameters (e.g. lambda) over a many covariance models or covariance parameter values.

Description

These function are designed to explore the likelihood surface for different covariance parameters with the option of maximizing over tau and sigma. They used the mKrig base are designed for computational efficiency.

Usage

mKrigMLEGrid(x, y, weights = rep(1, nrow(x)), Z = NULL,ZCommon=NULL,
                       mKrig.args = NULL,
                          cov.function = "stationary.cov", 
                         cov.args = NULL,
                           na.rm = TRUE, 
                         par.grid = NULL, 
                           reltol = 1e-06,
                             REML = FALSE,
                             GCV  = FALSE,
                       optim.args = NULL,
                 cov.params.start = NULL,
                          verbose = FALSE,
                            iseed = NA)
                          

mKrigMLEJoint(x, y, weights = rep(1, nrow(x)), Z = NULL,ZCommon=NULL,
                 mKrig.args = NULL,
		 na.rm = TRUE, cov.function = "stationary.cov",
                 cov.args = NULL, cov.params.start = NULL, optim.args =
                 NULL, reltol = 1e-06, parTransform = NULL, REML =
                 FALSE, GCV = FALSE, hessian = FALSE, iseed = 303,
		 verbose = FALSE)

profileCI(obj, parName, CIlevel = 0.95, REML = FALSE) 

mKrigJointTemp.fn(parameters, mKrig.args, cov.args, parTransform,
                 parNames, REML = FALSE, GCV = FALSE, verbose =
                 verbose, capture.env)
           

Arguments

capture.env	For the ML objective function the frame to save the results of the evaluation. This should be the environment of the function calling optim.
CIlevel	Confidence level.
cov.function	The name, a text string, of the covariance function.
cov.args	The arguments that would also be included in calls to the covariance function to specify the fixed part of the covariance model. This is the form of a list E.g.cov.args= list( aRange = 3.5)
cov.params.start	A list of initial starts for covariance parameters to perform likelihood maximization. The list contains the names of the parameters as well as the values. It usually makes sense to optimize over the important lambda parameter (tau^2/ sigma^2) is most spatial applications and so if lambda is omitted then the component lambda = .5 is added to this list.
hessian	If TRUE return the BFGS approximation to the hessian matrix at convergence.
iseed	Sets the random seed in finding the approximate Monte Carlo based GCV function and the effective degrees of freedom. This will not effect random number generation outside these functions.
mKrig.args	A list of additional parameters to supply to the mKrig function. E.g. mKrig.args= list( m=1) to set the regression function to be a constant function. mKrig function that are distinct from the covariance model. For example mKrig.args= list( m=1 ) will set the polynomial to be just a constant term (degree = m - 1 = 0). Use mKrig.args= list( m = 0 ) to omit a fixed model and assume the observations have an expectation of zero.
na.rm	Remove NAs from data.
optim.args	Additional arguments that would also be included in calls to the optim function in joint likelihood maximization. If NULL, this will be set to use the "BFGS-" optimization method. See optim for more details. The default value is: optim.args = list(method = "BFGS", control=list(fnscale = -1, ndeps = rep(log(1.1), length(cov.params.start)+1), abstol=1e-04, maxit=20)) Note that the first parameter is lambda and the others are the covariance parameters in the order they are given in cov.params.start. Also note that the optimization is performed on a transformed scale (based on the function parTransform ), and this should be taken into consideration when passing arguments to optim.
parameters	The parameter values for evaluate the likelihood.
par.grid	A list or data frame with components being parameters for different covariance models. A typical component is "aRange" comprising a vector of scale parameters to try. If par.grid is "NULL" then the covariance model is fixed at values that are given in ....
obj	List returnerd from mKrigMLEGrid
parName	Name of parameter to find confidence interval.
parNames	Names of the parameters to optimize over.
parTransform	A function that maps the parameters to a scale for optimization or effects the inverse map from the transformed scale into the original values. See below for more details.
reltol	Optim BFGS comvergence criterion.
REML	If TRUE use REML instead of the full log likelihood.
GCV	NOT IMPLEMENTED YET! A placeholder to implement optimization using an approximate cross-validation criterion.
verbose	If TRUE print out interesting intermediate results.
weights	Precision ( 1/variance) of each observation
x	Matrix of unique spatial locations (or in print or surface the returned mKrig object.)
y	Vector or matrix of observations at spatial locations, missing values are not allowed! Or in mKrig.coef a new vector of observations. If y is a matrix the columns are assumed to be independent observations vectors generated from the same covariance and measurment error model.
Z	Linear covariates to be included in fixed part of the model that are distinct from the default low order polynomial in x
ZCommon	Linear covariates to be included in fixed part of the model where a common parameter is estimated across all realizations. This option only makes sense for multiple realizations ( y is a matrix).

Details

The observational model follows the same as that described in the Krig function and thus the two primary covariance parameters for a stationary model are the nugget standard deviation (tau) and the marginal variance of the process (sigma). It is useful to reparametrize as sigma and lambda = tau^2/sigma. The likelihood can be maximized analytically over sigma and the parameters in the fixed part of the model, this estimate of sigma can be substituted back into the likelihood to give a expression that is just a function of lambda and the remaining covariance parameters. This operation is called concentrating the likelhood by maximizing over a subset of parameters

For these kind of computations there has to be some device to identify parameters that are fixed and those that are optimized. For mKrigMLEGrid and mKrigMLEJoint the list cov.args should have the fixed parameters. For example this is how to fix a lambda value in the model. The list cov.params.start should be list with all parameters to optimize. The values for each component are use as the starting values. This is how the optim function works.

These functions may compute the effective degrees of freedomn ( see mKrig.trace ) using the random tace method and so need to generate some random normals. The iseed arguement can be used to set the seed for this with the default being the seed 303. Note that the random number generation internal to these functions is coded so that it does not effect the random number stream outside these function calls.

For mKrigMLEJoint the default transformation of the parameters is set up for a log/exp transformation:

 parTransform <- function(ptemp, inv = FALSE) {
            if (!inv) {
                log(ptemp)
            }
            else {
                exp(ptemp)
            }
        }

Value

mKrigMLEGrid returns a list with the components:

summary	A matrix with each row giving the results of evaluating the likelihood for each covariance model.
par.grid	The par.grid argument used. A matrix where rows are the combination of parameters considered.
call	The calling arguments to this function.

mKrigMLEJoint returns a list with components:

summary	A vector giving the MLEs and the log likelihood at the maximum
lnLike.eval	A table containing information on all likelihood evaluations performed in the maximization process.
optimResults	The list returned from the optim function. Note that the parameters may be transformed values.
par.MLE	The maximum likelihood estimates.
parTransform	The transformation of the parameters used in the optimziation.

Author(s)

Douglas W. Nychka, John Paige

References

https://github.com/dnychka/fieldsRPackage

See Also

mKrig Krig stationary.cov optim

Examples

## Not run: 
#perform joint likelihood maximization over lambda and aRange. 
# NOTE: optim can get a bad answer with poor initial starts.
  data(ozone2)
  s<- ozone2$lon.lat
  z<- ozone2$y[16,]
  gridList<- list( aRange = seq( .4,1.0,length.out=20),
                 lambda = 10**seq( -1.5,0,length.out=20)
                 )
  par.grid<- make.surface.grid( gridList)
  out<- mKrigMLEGrid( s,z, par.grid=par.grid,
                          cov.args= list(smoothness=1.0,
                                     Covariance="Matern" )
                          )   
  outP<- as.surface( par.grid, out$summary[,"lnProfileLike.FULL"])
  image.plot( outP$x, log10(outP$y),outP$z,
               xlab="aRange", ylab="log10 lambda")
               

## End(Not run)
 
 ## Not run: 
  N<- 50
  set.seed(123)
  x<- matrix(runif(2*N), N,2)
  aRange<- .2
  Sigma<-  Matern( rdist(x,x)/aRange , smoothness=1.0)
  Sigma.5<- chol( Sigma)
  tau<- .1
  #  250 independent spatial data sets but a common covariance function 
  #    -- there is little overhead in
  #        MLE across independent realizations and a good test of code validity.
  M<-250
  F.true<- t( Sigma.5) %*% matrix( rnorm(N*M), N,M)
  Y<-  F.true +  tau* matrix( rnorm(N*M), N,M)

# find MLE for lambda with grid of ranges 
# and smoothness fixed in Matern                     
 par.grid<- list( aRange= seq( .1,.35,,8))
  obj1b<- mKrigMLEGrid( x,Y,
     cov.args = list(Covariance="Matern", smoothness=1.0), 
     cov.params.start=list( lambda = .5),
        par.grid = par.grid
                    )
  obj1b$summary # take a look
# profile over aRange
  plot( par.grid$aRange, obj1b$summary[,"lnProfileLike.FULL"],
    type="b", log="x")
 
## End(Not run)
  ## Not run: 
# m=0 is a simple switch to indicate _no_ fixed spatial drift
# (the default and highly recommended  is linear drift, m=2). 
# However, m=0 results in MLEs that are less biased, being the correct model
# -- in fact it nails it !

  obj1a<- mKrigMLEJoint(x,Y, 
                    cov.args=list(Covariance="Matern", smoothness=1.0), 
                    cov.params.start=list(aRange =.5, lambda = .5),
                     mKrig.args= list( m=0))
 
 test.for.zero( obj1a$summary["tau"], tau, tol=.007)
 test.for.zero( obj1a$summary["aRange"], aRange, tol=.015)
 

## End(Not run) 


##########################################################################
# A bootstrap example
# Here is a example of a more efficient (but less robust) bootstrap using 
# mKrigMLEJoint and tuned starting values
##########################################################################
## Not run: 
data( ozone2)
obj<- spatialProcess( ozone2$lon.lat,ozone2$y[16,] )

######### boot strap 
  set.seed(123)
  M<- 25
# create M indepedent copies of the observation vector
  ySynthetic<- simSpatialData( obj, M)
  bootSummary<- NULL
  
 aRangeMLE<- obj$summary["aRange"]
 lambdaMLE<- obj$summary["lambda"]
 
  for(  k in 1:M){
  cat( k, " " )
# here the MLEs are found using the easy top level level wrapper
# see mKrigMLEJoint for a more efficient strategy
  out <- mKrigMLEJoint(obj$x, ySynthetic[,k],
                 weights = obj$weights,
              mKrig.args = obj$mKrig.args,
                 cov.function = obj$cov.function.name,
                cov.args = obj$cov.args, 
        cov.params.start = list( aRange = aRangeMLE,
                                lambda = lambdaMLE)
                      )
  newSummary<- out$summary
  bootSummary<- rbind( bootSummary, newSummary)
  }
  
  cat(  " ", fill=TRUE )
  
  obj$summary
  stats( bootSummary)
  

## End(Not run)
## Not run: 
#perform joint likelihood maximization over lambda, aRange, and smoothness.  
#note: finding smoothness is not a robust optimiztion 
#      can get a bad answer with poor initial guesses.
obj2<- mKrigMLEJoint(x,Y, 
                      cov.args=list(Covariance="Matern"), 
                      cov.params.start=list( aRange = .18,
                                         smoothness = 1.1,
                                             lambda = .08),
                       )

#look at lnLikelihood  evaluations
obj2$summary
#compare to REML
obj3<- mKrigMLEJoint(x,Y, 
                      cov.args=list(Covariance="Matern"), 
                      cov.params.start=list(aRange = .18, 
                                       smoothness = 1.1,
                                           lambda = .08),
                       , REML=TRUE)
obj3$summary                      

## End(Not run)
## Not run: 
#look at lnLikelihood  evaluations

# check convergence of MLE to true fit with no fixed part
# 
obj4<- mKrigMLEJoint(x,Y, 
                      mKrig.args= list( m=0),
                      cov.args=list(Covariance="Matern", smoothness=1),
                      cov.params.start=list(aRange=.2, lambda=.1),
                       REML=TRUE)
#look at lnLikelihood  evaluations
obj4$summary
# nails it!

## End(Not run)

fields documentation built on Aug. 18, 2023, 1:06 a.m.