estimate.STmodel: Estimation of the Spatio-Temporal Model

Description

Estimates parameters of the spatio-temporal model using maximum-likelihood, profile maximum likelihood or restricted maximum likelihood (REML). The function uses the L-BFGS-B method in optim to maximise loglikeST.

Usage

# S3 method for STmodel
estimate(object, x, x.fixed = NULL, type = "p",
  h = 0.001, diff.type = 1, hessian.all = FALSE, lower = -15,
  upper = 15, method = "L-BFGS-B", control = list(trace = 3, maxit =
  1000), restart = 0, ...)
estimate(object, x, ...)

Arguments

object

STmodel object for which to estimate parameters.

Vector or matrix of starting point(s) for the optimisation. A vector will be treated as a single starting point. If x is a matrix the optimisation will be run using each column as a separate starting point. If x is a single integer then multiple starting points will be created as a set of constant vectors with the values of each starting point taken as seq(-5, 5, length.out=x). See details below.

x.fixed

Vector with parameter to be held fixed; parameters marked as NA will still be estimated.

type

A single character indicating the type of log-likelihood to use. Valid options are "f", "p", and "r", for full, profile or restricted maximum likelihood (REML).

h, diff.type

Step length and type of finite difference to use when computing gradients, see loglikeSTGrad.

hessian.all

If type!="f" computes hessian (and uncertainties) for both regression and log-covariance parameters, not only for log-covariance parameters. See value below.

lower, upper, method

Parameter bound and optimisation method, passed to optim.

control

A list of control parameters for the optimisation. See optim for details; setting trace=0 eliminates all ouput.

restart

Number of times to restart each optimisation if optim fails to converge; can sometimes resolve issues with L-BFGS-B line search.

...

Ignored additional arguments.

Value

estimateSTmodel object containing:

res.best

A list containing the best optimisation result; elements are described below. Selection of the best result is described in details above.

res.all

A list with all the optimisations results, each element contains (almost) the same information as res.best, e.g. res.all[[i]] contains optimisation results for the i:th starting point.

summary

A list with parameter estimates and convergence information for all starting points.

Details

The starting point(s) for the optimisation can either contain both regression parameters and log-covariances parameters for a total of loglikeSTdim(object)$nparam parameters or only contain the log-covariances covariances parameters i.e. loglikeSTdim(object)$nparam.cov parameters. If regression parameters are given but not needed (type!="f") they are dropped; if they are needed but not given they are inferred through a generalised least squares (GLS) computation, obtained by calling predict.STmodel.

If multiple starting points are used this function returns all optimisation results, along with an indication of the best result. The best result is determined by first evaluating which of the optimisations have converged. Convergence is determined by checking that the output from optim has convergence==0 and that the hessian is negative definite, i.e. all(eigen(hessian)$value < -1e-10). Among the converged optimisations the one with the highest log-likelihood value is then selected as the best result.

If none of the optimisations have converged the result with the highest log-likelihood value is selected as the best result.

Most of the elements in res.best (and in res.all[[i]]) are obtained from optim. The following is a brief description:

par: The best set of parameters found.
value: Log-likelihood value corresponding to par.
counts: The number of function/gradient calls.
convergence: 0 indicates successful convergence, see optim.
message: Additional information returned by optim.
hessian: A symmetric matrix giving the finite difference Hessian of the function par.
conv: A logical variable indicating convergence; TRUE if convergence==0 and hessian is negative definite, see details above.
par.init: The initial parameters used for this optimisation.
par.all: All parameters (both regression and log-covariance). Identical to par if type="f".
hessian.all: The hessian for all parameters (both regression and log-covariance). NOTE: Due to computational considerations hessian.all is computed only for res.best.

Examples

Run this code

# NOT RUN {
##load a model object
data(mesa.model)

##create vector of initial values
dim <- loglikeSTdim(mesa.model)
x.init <- cbind(c( rep(2, dim$nparam.cov-1), 0),
                c( rep(c(1,-3), dim$m+1), -3, 0))
rownames(x.init) <- loglikeSTnames(mesa.model, all=FALSE)

# }
# NOT RUN {
  ##estimate parameters
  est.mesa.model <- estimate(mesa.model, x.init, hessian.all=TRUE)
# }
# NOT RUN {
##time consuming estimation, load pre-computed results instead
data(est.mesa.model)

#estimation results
print(est.mesa.model)

##compare the estimated parameters for the two starting points
est.mesa.model$summary$par.all
##and values of the likelihood (and convergence info)
est.mesa.model$summary$status

##extract the estimated parameters and approximate uncertainties
x <- coef(est.mesa.model)

##compare estimated parameters
##plot the estimated parameters with uncertainties
par(mfrow=c(1,1),mar=c(13.5,2.5,.5,.5))
with(x, plot(par, ylim=range(c(par-1.96*sd, par+1.96*sd)),
                xlab="", xaxt="n"))
with(x, points(par - 1.96*sd, pch=3))
with(x, points(par + 1.96*sd, pch=3))

abline(h=0, col="grey")
##add axis labels
axis(1, 1:length(x$par), rownames(x), las=2)

# }
# NOT RUN {
  ##example using a few fixed parameters
  x.fixed <- coef(est.mesa.model)$par
  x.fixed[c(1,2,5:9)] <- NA
  est.fix <- estimate(mesa.model, x.init, x.fixed, type="p")
# }

Run the code above in your browser using DataLab