fitGLS_opt: Fit a PARTS GLS model, with maximum likelihood spatial parameters

Description

Fit a PARTS GLS model, with maximum likelihood spatial parameters

Usage

fitGLS_opt(
  formula,
  data = NULL,
  coords,
  distm_FUN = "distm_scaled",
  covar_FUN = "covar_exp",
  start = c(range = 0.01, nugget = 0),
  fixed = c(),
  opt.only = FALSE,
  formula0 = NULL,
  save.xx = FALSE,
  save.invchol = FALSE,
  no.F = TRUE,
  trans = list(),
  backtrans = list(),
  debug = TRUE,
  ncores = NA,
  ...
)

Value

If opt.only = TRUE, fitGLS_opt returns the output from stats::optim(): see it's documentation for more details.

Otherwise, a list with two elements is returned:

opt: output from optim, as above

GLS

a "remoteGLS" object. See fitGLS for more details.

Arguments

formula: a model formula, passed to fitGLS
data: an optional data frame environment in which to search for variables given by formula; passed to fitGLS
coords: a numeric coordinate matrix or data frame, with two columns and rows corresponding to each pixel
distm_FUN: a function to calculate a distance matrix from coords
covar_FUN: a function to estimate distance-based covariances
start: a named vector of starting values for each parameter to be estimated; names must match the names of arguments in covar_FUN or "nugget"
fixed: an optional named vector of fixed parameter values; names must match the names of arguments in covar_FUN or "nugget"
opt.only: logical: if TRUE, execution will halt after estimating the parameters; a final GLS will not be fit with the estimated parameters
formula0, save.xx, save.invchol, no.F: arguments passed to fitGLS for final GLS output
trans: optional list of functions for transforming the values in start or fixed in order to constrain the parameter space within optim
backtrans: optional list of functions for back-transforming parameters to their correct scale (for use with trans)
debug: logical: debug mode (for use with trans and backtrans)
ncores: an optional integer indicating how many CPU threads to use for calculations.
...: additional arguments passed to stats::optim()

Details

Estimate spatial parameters, via maximum likelihood, from data rather than from time series residuals; Fit a GLS with these specifications.

fitGLS_opt fits a GLS by estimating spatial parameters from data. fitCor, combined with fitGLS(nugget = NA), gives better estimates of spatial parameters, but time-series residuals may not be available in all cases. In these cases, spatial parameters can be estimated from distances among points and a response vector. Mathematical optimization of the log likelihood of different GLS models are computed by calling optim() on fitGLS.

Distances are calculated with distm_FUN and a covariance matrix is calculated from these distances with covar_FUN. Arguments to to covar_FUN, except distances, are given by start and fixed. Parameters specified in start will be be estimated while those given by fixed will remain constant throughout fitting. Parameter names in start and fixed should exactly match the names of arguments in covar_FUN and should not overlap (though, fixed takes precedence).

In addition to arguments of covar_FUN a "nugget" component can also be occur in start or fixed. If "nugget" does not occur in either vector, the GLS are fit with nugget = 0. A zero nugget also allows much faster computation, through recycling the common inverse cholesky matrix in each GLS computation. A non-zero nugget requires inversion of a different matrix at each iteration, which can be substantially slower.

If opt.only = FALSE, the estimated parameters are used to fit the final maximum likelihood GLS solution with fitGLS() and arguments formula0, save.xx, save.invchol, and no.F.

Some parameter combinations may not produce valid covariance matrices. During the optimization step messages about non-positive definitive V may result on some iterations. These warnings are produced by fitGLS and NA log-likelihoods are returned in those cases.

Note that fitGLS_opt fits multiple GLS models, which requires inverting a large matrix for each one (unless a fixed 0 nugget is used). This process is very computationally intensive and may take a long time to finish depending upon your machine and the size of the data.

Sometimes optim can have a difficult time finding a reasonable solution and without any constraits on parameter space (with certain algorithms), results may even be nonsensical. To combat this, fitGLS_opt has the arguments trans and backtrans which allow you to transform (and back-transform) parameters to a different scale. For example, you may want to force the 'range' parameter between 0 and 1. The logit function can do just that, as its limits are -Inf and Inf as x approaches 0 and 1, respectively. So, we can set trans to the logit function: trans = list(range = function(x)log(x/(1-x))). Then we need to set backtrans to the inverse logit function to return a parameter value between 0 and 1: backtrans = list(range = function(x)1/(1+exp(-x))). This will force the optimizer to only search for the range parameter in the space from 0 to 1. Any other constraint function can be used for trans provided that there is a matching back-transformation.

Examples

Run this code

# \donttest{
## read data
data(ndvi_AK10000)
df = ndvi_AK10000[seq_len(200), ] # first 200 rows

## estimate nugget and range (very slow)
fitGLS_opt(formula = CLS_coef ~ 0 + land, data = df,
            coords = df[, c("lng", "lat")], start = c(range = .1, nugget = 0),
            opt.only = TRUE)

## estimate range only, fixed nugget at 0, and fit full GLS (slow)
fitGLS_opt(formula = CLS_coef ~ 0 + land, data = df,
             coords = df[, c("lng", "lat")],
             start = c(range = .1), fixed = c("nugget" = 0),
             method = "Brent", lower = 0, upper = 1)

## constrain nugget to 0 and 1
logit <- function(p) {log(p / (1 - p))}
inv_logit <- function(l) {1 / (1 + exp(-l))}

fitGLS_opt(formula = CLS_coef ~ 0 + land, data = df,
           coords = df[, c("lng", "lat")],
           start = c(range = .1, nugget = 1e-10),
           trans = list(nugget = logit), backtrans = list(nugget = inv_logit),
           opt.only = TRUE)
# }