fit.mesa.model: 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 loglike.

Usage

fit.mesa.model(x, mesa.data.model, type = "p", h = 0.001,
    diff.type = 1, lower = -15, upper = 15, hessian.all = FALSE,
    control = list(trace = 3, maxit = 1000))

Arguments

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

mesa.data.model

Data structure holding observations, and information regarding which geographic and spatio-temporal covariates to use when fitting the model. See create.data.model and

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 loglike.grad and gen.gradient for details.

lower, upper

Bounds on the variables, see optim.

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.

control

A list of control parameters for the optimisation. See optim for details.

Value

Returns a list containing:
res.bestA list containing the best optimisation result; elements are described below. Selection of the best result is described in details above.
res.allA 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.
messageA text string with information regarding best value and number of converged optimisations.
Most of the elements in res.best (and in res.all[[i]]) are obtained from optim. The following is a brief description: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

encoding

latin1

Details

The function estimates parameters for the spatio-temporal model using the full likelihood formulation, profile likelihood, or restricted maximum likelihood (REML). In principal full likelihood and profile likelihood should give the same results, corresponding to the maximum likelihood estimate, with the full likelihood approach being slower.

The starting point(s) for the optimisation can either contain both regression parameters and log-covariances parameters for a total of loglike.dim(mesa.data.model)$nparam parameters or only contain the log-covariances covariances parameters i.e. loglike.dim(mesa.data.model)$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 cond.expectation(x, mesa.data.model, only.pars=TRUE).

If x is a single integer then that many starting points will be created as vectors with the constant values for the log-covariance parameters being in the sequence seq(-5, 5, length.out=x); if needed the corresponding regression parameters are computed using GLS.

The estimation uses the L-BFGS-B method in optim to maximise loglike. Gradient calculations are done by loglike.grad; when specifying the type of finite differences to use note that central differences (diff.type=0) will drastically increase estimation time, often with little or no benefit to the optimisation.

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.

Examples

Run this code

##load a model object
data(mesa.data.model)

##examine the model
printMesaDataNbrObs(mesa.data.model)

##covariates
mesa.data.model$LUR.list
mesa.data.model$ST.Ind

##Important dimensions of the model
dim <- loglike.dim(mesa.data.model)
print(dim)

##Set up initial parameter values for optimization
x.init <- cbind(rep(2,dim$nparam.cov), 
                c(rep(c(1,-3),dim$m+1),-3))
##and add names to the initial values
rownames(x.init) <- loglike.var.names(mesa.data.model,
                                      all=FALSE)
print(x.init)

##estimate parameters
##This may take a while...
par.est <- fit.mesa.model(x.init, mesa.data.model, type="p",
      hessian.all=TRUE, control=list(trace=3,maxit=1000))
##Let's load precomputed results instead.
data(mesa.data.res)
par.est <- mesa.data.res$par.est

##Optimisation status message
par.est$message

##compare the estimated parameters for the two starting points
cbind(par.est$res.all[[1]]$par.all,
      par.est$res.all[[2]]$par.all)
##and values of the likelihood
cbind(par.est$res.all[[1]]$value,
      par.est$res.all[[2]]$value)

##extract the estimated parameters
x <- par.est$res.best$par.all

##and approximate uncertainties from the hessian
x.sd <- sqrt(diag(-solve(par.est$res.best$hessian.all)))
names(x.sd) <- names(x)

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

Run the code above in your browser using DataLab