estimate.STmodel: Estimation of the Spatio-Temporal Model
In SpatioTemporal: Spatio-Temporal Model Estimation

Description Usage Arguments Details Value Author(s) See Also Examples

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.

## S3 method for class '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, ...)

`object`	`STmodel` object for which to estimate parameters.
`x`	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.

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.

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.

Johan Lindstrom

Other STmodel methods: MCMC.STmodel, c.STmodel, createSTmodel, estimateCV.STmodel, plot.STdata, predict.STmodel, print.STmodel, print.summary.STmodel, qqnorm.predCVSTmodel, scatterPlot.predCVSTmodel, simulate.STmodel, summary.STmodel

Other estimateSTmodel methods: coef.estimateSTmodel, print.estimateSTmodel

##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)

## End(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")

## End(Not run)