Skip to contents

Changelog

LaMa 3.1.0

  • tpm_ct() is much faster: with AD and 2 states it uses an exact closed form instead of Matrix::expm() for each interval (about 100x faster gradients), and for a constant Q it computes the matrix exponential only once per unique time difference
  • qreml() no longer evaluates the Hessian a second time for the returned model, which saves about length(par) gradient evaluations per fit, and the final message now reads "Final sp: ..." as no fit follows it unless the last iterate was not the best one
  • MCreport() for qreml() models samples around the final estimate, and with joint = TRUE (or Q from joint_uncertainty()) accounts for the uncertainty in the penalty strengths; needs RTMBdist >= 1.2.1
  • residuals(mod) gives the pseudo-residuals, calling pseudo_res(mod = mod)
  • coef(mod) and vcov(mod) give the estimated fixed parameters and their covariance from the inverse Hessian (of the marginal likelihood with random effects, conditional on the smoothing parameters for qreml() models), so confint(mod) works too
  • qreml() no longer tries doubled steps after improving ones: each trial costs a full inner fit, and across the case studies tried it almost never reduced the number of inner fits
  • qreml(lsp_min = -15) bounds log(lambda) from below, so a penalty strength on a smooth the data cannot inform stops falling and the iteration can converge; a penalty strength ending there gives a warning, and summary() notes it: a larger basis dimension is suggested if no other penalty acts on the same coefficients, otherwise the penalty is reported as redundant (e.g. the wiggliness penalty of a tensor product with a random effect margin), which a larger basis does not help
  • qreml(): a step that changes the restricted likelihood by less than 1e-6 * (1 + |criterion|) counts as unchanged, instead of setting off up to max_halve refits at essentially the same penalty strengths at the optimum
  • qreml() also converges once the inner optimiser no longer moves from its warm start, so a penalty strength creeping towards zero along a flat restricted likelihood no longer runs on until lsp_min or maxiter
  • qreml() returns n_fits, the number of inner fits including rejected steps, shown by summary() next to the number of outer iterations
  • tpm() and generator() with a design matrix give a one-time message while taping if TapeConfig(matmul = 'plain') would record more than 1e8 multiply-adds for the linear predictor (counting only nonzero entries of Z, which RTMB skips otherwise), suggesting TapeConfig(matmul = 'atomic') as the tape may not fit into memory
  • qreml() returns edf_hist, the effective degrees of freedom of each smooth at the start of each outer iteration (column 1 is the initial fit)
  • qreml(): outer_hessian() now differentiates the approximate gradient of the restricted log-likelihood w.r.t. log(lambda) by central differences (2 inner fits per free penalty strength, step h = 0.1), instead of one-sided differences with an absolute step of 1e-4 in lambda, which was inaccurate and could give an indefinite Hessian. The fits are warm started at the final estimate with a tighter inner tolerance. The result is now on the log scale
  • qreml(): outer_hessian() no longer leaves the objective at the last perturbed penalty strength, so samples or Hessians computed from mod$obj afterwards use the fitted penalty strengths again
  • qreml() returns joint_uncertainty(), the joint precision matrix of the parameters and the log penalty strengths (the first-order correction of Wood, Pya and Saefken, 2016, for the uncertainty in the penalty strengths, in the form of TMB’s sdreport(getJointPrecision = TRUE)), which stays sparse with spHess = TRUE; with precision = FALSE the corrected covariance of the parameters
  • sdreport_outer() works with the outer Hessian w.r.t. log(lambda) and obtains the standard deviations of lambda and 1/lambda from it by the delta method

LaMa 3.0.0

CRAN release: 2026-10-28

This is a major release: models are now built around statedens() and lama(), which give print(), summary(), plot(), simulate(), pseudo_res() and scores() from the model object alone. Several functions are renamed (old names still work but are deprecated), a few are removed, and some outputs change (e.g. llk instead of ll, sp instead of lambda), so existing code may need small adjustments.

New features

  • statedens() builds the matrix of state-dependent densities of one data stream for allprobs, with parameters given per state, for all states or per observation and state; missing observations are handled for each stream separately, and log = TRUE gives log-densities for forward(..., logspace = TRUE); pseudo_res(mod) works with it

  • forward() without AD gives an error for negative allprobs, suggesting logspace = TRUE

  • forward(delta = "stationary") computes the initial distribution from Gamma: the stationary distribution (per track if track-specific, of the first slice of each track if time-varying) or, with tod, the periodically stationary distribution at the first time point of each track

  • plot(mod, type = "statedens") shows the state-dependent densities of each data stream built with statedens() on a histogram, weighted by the local state probabilities (or decoded states, the stationary distribution or given weights), which also works with covariates, several tracks and inhomogeneous models; the weighted densities are returned invisibly

  • plot(mod, stream = ...) selects data streams, for residuals and state-dependent densities

  • plot(mod, type = "statedens", xlim = ...) evaluates the density curves on the plotted range, not only the range of the data

  • simulate(mod) simulates states and observations of models built with statedens() on the observed design: same tracks and lengths, with delta, Gamma and the parameters as reported (conditional on the covariates, for continuous-time models on the observed time points); missing observations stay missing unless fill_na = TRUE, and rdens gives random number generators for densities without a matching r-function; a data stream that is NA throughout is a placeholder, so a likelihood function can be used to simulate data and then be fitted to them

  • simulate(mod, random = TRUE) and obj$simulate() of RTMB work with forward() and statedens(): statedens() registers its data like OBS(x), so random effects (with their density evaluated first, e.g. u %~% dnorm(...)) are drawn first, then the states and data given them; forward() reports the simulated states as states

  • print() and summary() for models from lama() and qreml(): print() gives a short overview, summary() the estimated state process (transition probability or generator matrix, stationary distribution and mean dwell times, averaged over tracks if track-specific) and the reported quantities, with those of length nStates in one table; summary.qremlModel() adds the smooth terms and returns an object now

  • qreml() models have class c("qremlModel", "LaMaModel"), so all LaMaModel methods apply to them

  • forward() reports type "track-specific" for one Gamma per track (before "homogeneous")

  • multistart(obj) optimises an RTMB object from several starting values (obj$par and perturbations of it, or of the first fit scaled by the unit information with scale = "unit_info"), stops by a Bayesian stopping rule or once the best optimum is reached often enough, and returns the best fit with a table of all runs; lama(obj) then builds the model from the best fit

  • modmat() returns sp0, initial penalty strengths from mgcv::initial.sp() (times 100, capped at 1e6), one per penalty matrix and named as S; rep_sp() repeats them as lambda for smooths used more than once, e.g. one per state

  • The penalty strengths are now called sp: qreml() takes them from dat$sp if present, otherwise from dat$lambda as before, so lambda is free for other things, e.g. a Poisson mean; spname still overrides this. penalty() and penalty2() call their argument sp, lambda still works

  • print() for objects from modmat() shows the formulas and the initial penalty strengths

  • smoothdens() builds the objects for smooth (P-spline) densities from a data vector, or one object per stream from a data frame or list; density(beta) in the returned object gives the allprobs matrix and predict() the estimated densities, e.g. on the grid xseq

Changes

Bug fixes

  • Examples and vignettes: the periodically stationary initial distribution matches forward()’s tod, where Gamma[,,tod[t]] is the transition into t, i.e. stationary_p(Gamma, tod[1] %% L + 1)

  • lama(): the log-likelihood was evaluated at the last parameters the optimiser tried, not the best

  • summary.qremlModel() failed for reported lists or character vectors, and if Gamma was not reported but delta was

  • pseudo_res(): missing observations give NA in all variants (before 0 with normal = FALSE, -Inf for non-randomised discrete residuals)

  • pseudo_res(): parameters given per observation (matrix or vector of length nObs) work with missing observations

  • forward(): with track-specific delta or Gamma, the AD version matched them to the wrong tracks if the track IDs were not sorted (e.g. "ID1", ..., "ID10")

  • forward(): tracks of length one failed with AD

  • logspace = TRUE is now exact (delta is no longer modified) and also stable without AD

  • one-state models are correct with logspace = TRUE and report their quantities

  • forward_p() works if tod does not contain all values 1, ..., L

  • forward_g() works without AD for two observations

  • stateprobs_g(forecast = TRUE) with trackID used the transition matrix of the previous time point, which also affected pseudo_res() for inhomogeneous models with tracks

  • stateprobs() passes forecast on for time-varying Gamma

  • viterbi() and stateprobs() work for tracks of length one

  • stateprobs(forecast = TRUE): the first row of each track is delta, not the filtered probabilities given x_1, so pseudo_res() is correct for the first observation

  • modmat(): smooths with several penalties but not one per margin (e.g. bs = "fs", "ad" or t2(full = TRUE)) kept only the first, so later smooths got the wrong penalty matrices

  • modmat(): a nested formula list failed if one data stream had no smooth terms

  • modmat(): tensor product margins over several covariates are named e.g. "x,y" (two such margins failed); unpenalised smooths (fx = TRUE) give an informative error

  • predict() on objects from modmat() ignored exclude

  • smooth_dens_construct(): coefficients of one state are a matrix (before a vector, so cbind(coef, 0) was wrong); circular data with custom knots failed; angles outside [-pi, pi) are wrapped; par accepts kappa next to concentration

  • make_matrices_dens(): the basis is normalised on a midpoint grid (before a fixed 1e-3 offset, which gave negative weights for data on a small scale), which changes the weights very slightly

Minor changes

LaMa 2.1.1 (2026-05-13)

CRAN release: 2026-05-13

Major changes

  • vignettes 1-3 have been updated to work with automatic differentiation

  • tpm(), forward(), stateprobs(), viterbi() can now be used with inhomogeneous models: call their inhomogeneous counterparts internally

  • function MCreport() has been added for sampling

    • from the distribution of the MLE (fixed effects only models)
    • from the joint post. distribution of parameters and random effects (random effects models)

Minor changes