Method Arguments
The estimate method accepts the following arguments:
model$estimate(
data,
method = "ekf",
ode.solver = "euler",
ode.timestep = diff(data$t),
loss = "quadratic",
loss_c = NULL,
control = list(trace=1,iter.max=1e5,eval.max=1e5),
use.hessian = FALSE,
laplace.residuals = FALSE,
unconstrained.optim = FALSE,
estimate.initial.state = FALSE,
silent = FALSE
)method
The method argument decides which estimation/filtering
algorithm is used. The following methods are available:
method='lkf': The Linear Kalman Filter.method='ekf': The Extended Kalman Filter (default).method='ukf': The Unscented Kalman Filter.method='laplace': The Laplace Approximation.
There is a difference between the Kalman filtering methods, and the Laplace method in the way they approach likelihood computations, and thus what information the filters produce:
- The Kalman filters produce prior and posterior state estimates. The prior estimates are conditioned on observations from the starting time \(t_0\) up til the previous time-point \(t_{k-1}\). The posterior estimates are conditioned also on the observations available at the “current” time-point \(t_{k}\). These are denoted respectively by
\[ \text{Prior:} \quad \mathrm{E}\left[ x_{t_k} \mid y_{t_{k-1}}, y_{t_{k-2}},...,y_{t_0} \right] \]
\[ \text{Posterior:} \quad \mathrm{E}\left[ x_{t_k} \mid y_{t_{k}}, y_{t_{k-1}},...,y_{t_0} \right] \]
The Laplace filter produces smoothed state estimates. This is a state estimate based observations at all time-points \(t_{0}\) to \(t_{N}\) where \(N\) is the last index in the time-series. These are denoted by
\[ \text{Smoothed:} \quad \mathrm{E}\left[ x_{t_k} \mid y_{N}, y_{N-1},...,y_{t_{k}},...,y_{t_1}, y_{t_0} \right] \]
- The likelihood contributions from the Kalman filters are based on
the prior estimates, which are one-step-ahead predictions. This
gives rise to independent one-step-ahead residuals, ideal for
residual analysis and goodness-of-fit model validation. The Laplace
filter does not produce such residuals inherently, and must instead
compute these on the side. Residual calculations are disabled by default
and determined with the
laplace.residualsargument. The computations are very costly and slow. The user is referred to the documentation inTMB::oneStepPredictfor further information.
ode.solver
This argument is used for the Kalman filtering methods only to determine the ODE integrator used to solving the moment differential equations i.e.
\[ \dfrac{d\mathrm{E}\left[x_t\right]}{dt} = f(\mathrm{E}\left[x_t\right]) \]
\[ \dfrac{d\mathrm{V}\left[x_t\right]}{dt} = \dfrac{df}{dx}\bigg\vert_{(\mathrm{E}\left[x_t\right])}\mathrm{V}\left[x_t\right] + \mathrm{V}\left[x_t\right]\left[\dfrac{df}{dx}\bigg\vert_{(\mathrm{E}\left[x_t\right])}\right]^{T} + g(\mathrm{E}\left[x_t\right])g(\mathrm{E}\left[x_t\right])^{T} \]
The ctsmTMB package implements the Explicit Forward-Euler euler and the Explicit 4th Order Runge-Kutta rk4 methods.
ode.timestep
This argument has two different implications depending on whether
Kalman filtering or Laplace filtering is carried out. The method accepts
either a single scalar value to be used as the global time-step, or a
vector of length diff(data$t) specifying all individual
time-steps. The input values are interpolated linearly between
time-points if more than one step is taken.
Kalman filters: In this case the argument controls the time-step used for the euler and rk4 ODE solvers.
Laplace filter: In this case the argument controls the number of added intermediate time-points between observations each of which represents an additional state (random effect) at that particular time-point.
If a provided time-step \(\Delta t_{i}\) does not divide the correspond time-difference \(t_{i+1}-t_{i}\) in the data then it is rounded down such that it does. Consider the following example where the time-difference between two observations in the data is 3 seconds, but a time-step of 0.7. This produces a non-integer number of steps i.e.:
\[ N_{i} = \dfrac{t_{i+1}-t_{i}}{\Delta t_i} = \dfrac{3}{0.7} = 4.28... \]
Thus the number of steps taken is rounded up to \(N^{*}_i = \left\lceil N_i \right\rceil = \left\lceil 4.28... \right\rceil = 5\) and the corrected time-step then becomes
\[ \Delta t^{*}_{i} = \dfrac{3}{N^{*}_{i}} = \dfrac{3}{5} = 0.6 \]
Note: The exception to this rule is when the remainder is less than \(\epsilon = 10^{-3}\) i.e. if for instance \(N_i = 4.0001\) then the time-step is accepted, and the number of steps rounded down to \(N^{*}_{i}= \left\lfloor N_i \right\rfloor = 4\).
loss / loss_c
The following losses are currently available:
loss='quadratic'(default)loss='huber'loss='tukey'
This argument only affects the Kalman filtering methods, and is used to regularize the likelihood contributions, removing the influence of large outliers.
The ith likelihood contributions is given by:
\[ -\log L_{i}(\theta) \propto f(r_i) \]
where \(r_i = \sqrt{e_{i}^{T} \Sigma_{i}^{-1} e_{i}}\) is a normalized residual, and where \(e_{i}\) is the ith residual vector, and \(\Sigma_{i}^{-1}\) the ith residual precision matrix.
The loss argument changes the function \(f\) as follows:
- If
loss='quadratic'then \(f\) is quadratic in the residuals i.e.
\[ f(r) = r^2 \]
and the likelihood contributions are exactly those from a Gaussian.
- If
loss='huber'then \(f\) is the Huber’s \(\psi\) function given by
\[ \psi_{c}(r) = \left\{ \begin{array}{l} r^2 & \text{for} \,\, r \leq c \\ c(2r-c) & \text{otherwise} \end{array} \right\} \]
which is quadratic/linear in the residuals below/above the threshold determined by \(c\).
- If
loss='tukey'then \(f\) is Tukey’s byweight function given by
\[ l_{c}(r) = \left\{ \begin{array}{l} r^2 & \text{for} \,\, r \leq c \\ c^2 & \text{otherwise} \end{array} \right\} \]
which is quadratic/constant in the residuals below/above the threshold determined by \(c\).
In practice a smooth approximation to both Huber and Tukey are implemented in practice using the construction
\[ \tilde{\psi}_{c}(r) = r^2 (1-\sigma_{c}(r)) + c^2 \sigma_{c}(r) \]
\[ \tilde{l}_{c}(r) = r^2 (1-\sigma_{c}(r)) + c(2r-c) \sigma_{c}(r) \]
where \(\sigma(r)\) is the sigmoid function
\[ \sigma_{c}(r) = \dfrac{1}{1+\exp(-5(r-c))} \]
The plot below shows the actual and implemented loss functions (almost indistinguishable) for \(c=5\). The threshold value is marked by with dashed line for the line \(x = c\).
Loss Functions
The loss_c argument is used to determine the value of
\(c\). The default values are chosen
based on the fact that under the assumed (multivariate) normal
distribution of the residuals, the squared (normalized) residuals
follows a \(\chi^{2}_{m}\) distribution
with degrees of freedom equal to the number of elements of \(e_{i}\) i.e:
\[ r_{i}^2 = e_{i}^{T} \Sigma_{i}^{-1} e_{i} \sim \chi^{2}_{m} \]
It is therefore reasonable to choose the threshold level (value of \(c\)) which determines whether or not \(r_{i}\) is an outlier, as the level at which the null-hypothesis
\[ H_{0}: r_{i}^2 \sim \chi^{2}_{m} \]
is rejected for some critical value \(1 - \alpha\). Choosing the significance level \(\alpha = 0.05\) the appropriate \(c\) threshold value becomes
where \(m\) is the number of observation equations.
Note: The significance level will be higher than the chosen if there are missing observations in some indices of \(i\) for systems with multiple observation equations.
use_hessian
This argument is a boolean which determines whether or not the
likelihood hessian constructed by automatic differentiation from
TMB is used during the optimization procedure by
providing it to the hessian argument of
stats::nlminb. The default is
use.hessian=FALSE. The argument has no effect if
method=laplace due to restrictions in
TMB/RTMB.
One advantage of using the hessian is that the optimizer is more robust against parameters being on different scales, and the minimum is typically found in fewer (but more expensive) iterations. It is usually faster just to use the gradient.
laplace.residuals
This boolean controls whether or not model residuals are calculated
with TMB::oneStepPredict when the Laplace filtering is used
(method=laplace). This takes a considerable amount of time
- typically much longer than the estimation itself.
unconstrainted.optim
This boolean allows for quick unconstrained estimation removing the
parameter boundaries specified by setParameter. This may be
useful sometimes, to quickly check whether estimation issues occur
because of ill boundaries.
estimate.initial.state
This boolean determines whether or not the initial-time state
distribution (mean and covariance) should be estimated or not. By
default estimate.initial.state=FALSE and it is not
estimated but simply taken as the values provided by
setInitialState. When
estimate.initial.state=TRUE the mean and covariance are
estimated as the stationary solution of the moment differential
equations using the input values are the initial time-point.
The stationary solution is obtained first by solving
\[ \dfrac{d\mathrm{E}\left[x_{\infty}\right]}{dt} = f(\mathrm{E}\left[x_\infty\right]) = 0 \]
for \(\mathrm{E}\left[x_\infty\right]\) using Newton’s method. This stationary mean is then used to solve
control
This argument is a list which controls various settings of the
stats::nlminb optimizer. See the documentation on
?stats::nlminb for more information.
The default is list(trace=1,iter.max=1e5,eval.max=1e5)
which prints the iteration steps, and increases the default number of
iterations and function calls allowed before the optimization procedure
terminates.
Note:: The user should remember that disabling
tracing by passing control = list(trace=0) will remove the
ìter.maxand eval.max arguments, so they should
be provided as well if needed.
silent
This boolean argument controls whether or not various information messages are printed by ctsmTMB during model building, compilation and estimation.