Construct objective functions with derivatives based on a compiled C++ template.

Construct objective functions with derivatives based on the users C++ template.

Usage

MakeADFun(
  data,
  parameters,
  map = list(),
  type = c("ADFun", "Fun", "ADGrad"[!intern && (!is.null(random) || !is.null(profile))]),
  random = NULL,
  profile = NULL,
  random.start = expression(last.par.best[random]),
  hessian = FALSE,
  method = "BFGS",
  inner.method = "newton",
  inner.control = list(maxit = 1000),
  MCcontrol = list(doMC = FALSE, seed = 123, n = 100),
  ADreport = FALSE,
  atomic = TRUE,
  LaplaceNonZeroGradient = FALSE,
  DLL = getUserDLL(),
  checkParameterOrder = TRUE,
  regexp = FALSE,
  silent = FALSE,
  intern = FALSE,
  integrate = NULL,
  ...
)

Arguments

data

List of data objects (vectors, matrices, arrays, factors, sparse matrices) required by the user template (order does not matter and un-used components are allowed).

parameters

List of all parameter objects required by the user template (both random and fixed effects).

map

List defining how to optionally collect and fix parameters - see details.

type

Character vector defining which operation stacks are generated from the users template - see details.

random

Character vector defining the random effect parameters. See also regexp.

profile

Parameters to profile out of the likelihood (this subset will be appended to random with Laplace approximation disabled).

random.start

Expression defining the strategy for choosing random effect initial values as function of previous function evaluations - see details.

hessian

Calculate Hessian at optimum?

method

Outer optimization method.

inner.method

Inner optimization method (see function "newton").

inner.control

List controlling inner optimization.

MCcontrol

List controlling importance sampler (turned off by default).

ADreport

Calculate derivatives of macro ADREPORT(vector) instead of objective_function return value?

atomic

Allow tape to contain atomic functions?

LaplaceNonZeroGradient

Allow Taylor expansion around non-stationary point?

DLL

Name of shared object file compiled by user (without the conventional extension, .so, .dll, ...).

checkParameterOrder

Optional check for correct parameter order.

regexp

Match random effects by regular expressions?

silent

Disable all tracing information?

intern

Do Laplace approximation on C++ side ? See details (Experimental - may change without notice)

integrate

Specify alternative integration method(s) for random effects (see details)

...

Currently unused.

Value

List with components (fn, gr, etc) suitable for calling an R optimizer, such as nlminb or optim.

Details

A call to MakeADFun will return an object that, based on the users DLL code (specified through DLL), contains functions to calculate the objective function and its gradient. The object contains the following components:

• par A default parameter.

• fn The likelihood function.

• gr The gradient function.

• report A function to report all variables reported with the REPORT() macro in the user template.

• env Environment with access to all parts of the structure.

and is thus ready for a call to an R optimizer, such as nlminb or optim. Data (data) and parameters (parameters) are directly read by the user template via the macros beginning with DATA_ and PARAMETER_. The order of the PARAMETER_ macros defines the order of parameters in the final objective function. There are no restrictions on the order of random parameters, fixed parameters or data in the template.

Note

Do not rely upon the default arguments of any of the functions in the model object obj$fn, obj$gr, obj$he, obj$report. I.e. always use the explicit form obj$fn(obj$par) rather than obj$fn().

Parameter mapping

Optionally, a simple mechanism for collecting and fixing parameters from R is available through the map argument. A map is a named list of factors with the following properties:

• names(map) is a subset of names(parameters).

• For a parameter "p" length(map$p) equals length(parameters$p).

• Parameter entries with NAs in the factor are fixed.

• Parameter entries with equal factor level are collected to a common value.

More advanced parameter mapping, such as collecting parameters between different vectors etc., must be implemented from the template.

Specifying random effects

Random effects are specified via the argument random: A component of the parameter list is marked as random if its name is matched by any of the characters of the vector random (Regular expression match is performed if regexp=TRUE). If some parameters are specified as random effects, these will be integrated out of the objective function via the Laplace approximation. In this situation the functions fn and gr automatically perform an optimization of random effects for each function evaluation. This is referred to as the 'inner optimization'. Strategies for choosing initial values of the inner optimization can be controlled via the argument random.start. The default is expression(last.par.best[random]) where last.par.best is an internal full parameter vector corresponding to the currently best likelihood. An alternative choice could be expression(last.par[random]) i.e. the random effect optimum of the most recent - not necessarily best - likelihood evaluation. Further control of the inner optimization can be obtained by the argument inner.control which is a list of control parameters for the inner optimizer newton. Depending of the inner optimization problem type the following settings are recommended:

1. Quasi-convex: smartsearch=TRUE (the default).

2. Strictly-convex: smartsearch=FALSE and maxit=20.

3. Quadratic: smartsearch=FALSE and maxit=1.

The model environment env

Technically, the user template is processed several times by inserting different types as template parameter, selected by argument type:

• "ADFun" Run through the template with AD-types and produce a stack of operations representing the objective function.

• "Fun" Run through the template with ordinary double-types.

• "ADGrad" Run through the template with nested AD-types and produce a stack of operations representing the objective function gradient.

Each of these are represented by external pointers to C++ structures available in the environment env.

Further objects in the environment env:

• validpar Function defining the valid parameter region (by default no restrictions). If an invalid parameter is inserted fn immediately return NaN.

• parList Function to get the full parameter vector of random and fixed effects in a convenient list format.

• random An index vector of random effect positions in the full parameter vector.

• last.par Full parameter of the latest likelihood evaluation.

• last.par.best Full parameter of the best likelihood evaluation.

• tracepar Trace every likelihood evaluation ?

• tracemgc Trace maximum gradient component of every gradient evaluation ?

• silent Pass 'silent=TRUE' to all try-calls ?

The argument intern

By passing intern=TRUE the entire Laplace approximation (including sparse matrix calculations) is done within the AD machinery on the C++ side. This requires the model to be compiled using the 'TMBad framework' - see compile. For any serious use of this option one should consider compiling with supernodal=TRUE - again see compile - in order to get performance comparable to R's matrix calculations. The benefit of the 'intern' LA is that it may be faster in some cases and that it provides an autodiff hessian (obj$he) wrt. the fixed effects which would otherwise not work for random effect models. Another benefit is that it gives access to fast computations with certain hessian structures that do not meet the usual sparsity requirement. A detailed list of options are found in the online doxygen documentation in the 'newton' namespace under the 'newton_config' struct. All these options can be passed from R via the `inner.control` argument. However, there are some drawbacks of running the LA on the C++ side. Notably, random effects are no longer visible in the model environment which may break assumptions on the layout of internal vectors (`par`, `last.par`, etc). In addition, model debugging becomes harder when calculations are moved to C++.

Controlling tracing

A high level of tracing information will be output by default when evaluating the objective function and gradient. This is useful while developing a model, but may eventually become annoying. Disable all tracing by passing silent=TRUE to the MakeADFun call.