Package {LaMa}

Type:	Package
Title:	Fast Numerical Maximum Likelihood Estimation for Latent Markov Models
Version:	3.0.0
Description:	A variety of latent Markov models, including hidden Markov models, hidden semi-Markov models, state-space models and continuous-time variants can be formulated and estimated within the same framework via directly maximising the likelihood function using the so-called forward algorithm. Applied researchers often need custom models that standard software does not easily support. Writing tailored 'R' code offers flexibility but suffers from slow estimation speeds. We address these issues by providing easy-to-use functions (written in 'C++' for speed) for common tasks like the forward algorithm. These functions can be combined into custom models in a Lego-type approach, offering up to 10-20 times faster estimation via standard numerical optimisers. To aid in building fully custom likelihood functions, several vignettes are included that show how to simulate data from and estimate all the above model classes.
URL:	https://janolefi.github.io/LaMa/
License:	MIT + file LICENSE
Encoding:	UTF-8
Imports:	Rcpp, stats, utils, methods, Matrix, splines2, mgcv, MASS, numDeriv, RTMBdist (≥ 1.1.0)
LinkingTo:	Rcpp, RcppArmadillo
Depends:	R (≥ 3.5.0), RTMB (≥ 2.0)
RoxygenNote:	7.3.3
Suggests:	knitr, rmarkdown, testthat (≥ 3.0.0), PHSMM, MSwM, scales
VignetteBuilder:	knitr
Config/testthat/edition:	3
LazyData:	true
NeedsCompilation:	yes
Packaged:	2026-10-28 13:04:00 UTC; jan-ole
Author:	Jan-Ole Fischer [aut, cre]
Maintainer:	Jan-Ole Fischer <jan-ole.fischer@mailbox.org>
Repository:	CRAN
Date/Publication:	2026-10-28 14:30:08 UTC

LaMa: Fast Numerical Maximum Likelihood Estimation for Latent Markov Models

Description

A variety of latent Markov models, including hidden Markov models, hidden semi-Markov models, state-space models and continuous-time variants can be formulated and estimated within the same framework via directly maximising the likelihood function using the so-called forward algorithm. Applied researchers often need custom models that standard software does not easily support. Writing tailored 'R' code offers flexibility but suffers from slow estimation speeds. We address these issues by providing easy-to-use functions (written in 'C++' for speed) for common tasks like the forward algorithm. These functions can be combined into custom models in a Lego-type approach, offering up to 10-20 times faster estimation via standard numerical optimisers. To aid in building fully custom likelihood functions, several vignettes are included that show how to simulate data from and estimate all the above model classes.

Author(s)

Maintainer: Jan-Ole Fischer jan-ole.fischer@mailbox.org (ORCID)

See Also

Useful links:

• https://janolefi.github.io/LaMa/

Sparsity-retaining matrix multiplication

Description

Standard matrix multiplication destroys automatic sparsity detection by RTMB which is essential for models with high-dimensional random effects. This can be mitigated by changing to "plain" with TapeConfig, but this can make AD tape construction very slow. Here, we provide a different version that retains sparsity. It may be slightly slower than the standard method when constructing the AD tape.

Usage

A %sp% B

Arguments

Value

the matrix product of A and B, which is of dimension n x m

Examples

A <- matrix(1:6, nrow = 2, ncol = 3)
B <- matrix(7:12, nrow = 3, ncol = 2)
A %sp% B

Generate a colour-blind-friendly palette

Description

Generate a colour-blind-friendly palette

Usage

LaMaColors(ncolors)

Arguments

Value

A colour palette of length ncolors

Examples

cb <- LaMaColors(3)

Sample parameters from approximate Gaussian posterior distribution

Description

Efficient Monte Carlo sampling of parameters (and REPORTed quantities) from the approximate posterior of an (R)TMB model. See sdreport for details on posterior variance-covariance in random effects models.

Usage

MCreport(
  obj,
  nSamples = 1000,
  include_random_pars = TRUE,
  report = FALSE,
  Q = NULL,
  ...
)

Arguments

Details

Caution: This has nothing to do with Bayesian posterior sampling. It is simple Monte Carlo sampling from a Gaussian distribution with mean at the MLE and covariance given by the inverse of the Hessian (for fixed effects models) or the joint precision matrix (for random effects models). This is a common approach to get an approximate idea of parameter uncertainty around the MLE, but it relies on large-sample asymptotics.

Sampling random effects from their posterior as compared to calculating marginal standard devitions like sdreport is particularly useful for multidimensional random effects (e.g. for locations x and y) where pointwise confidence intervals (e.g. along a path) based on standard deviations are not possible.

Value

A list structured like the original parameter list used in the MakeADFun call (potentially including additional REPORTed quantities). Each entry is a list with nSamples entries.

Examples

step <- trex$step[1:2000] # subsetting trex data
N <- 2                    # 2 states

# custom likelihood
nll <- function(par) {
  getAll(par)
  Gamma <- tpm(eta)
  delta <- stationary(Gamma)
  mu <- exp(log_mu); REPORT(mu)
  sigma <- exp(log_sigma); REPORT(sigma)
  allprobs <- matrix(1, length(step), N)
  for(j in 1:N) allprobs[,j] <- dgamma2(step, mu[j], sigma[j])
  -forward(delta, Gamma, allprobs)
}

# initial parameters in named list
par0 <- list(eta = rep(-1.5,2), 
             log_mu = log(c(0.3, 2.5)), 
             log_sigma = log(c(0.3, 1.5)))
    
# constructing AD object        
obj <- MakeADFun(nll, par0, silent = TRUE)

# optimising
opt <- nlminb(obj$par, obj$fn, obj$gr)

# sampling from distribution of the MLE
samples <- MCreport(obj, nSamples = 10, report = TRUE)

Calculate the index of the first observation of each track based on an ID variable

Description

Function to conveniently calculate the trackInd variable that is needed internally when fitting a model to longitudinal data with multiple tracks.

Usage

calc_trackInd(ID)

Arguments

Value

A vector of indices of the first observation of each track which can be passed to the forward and forward_g to sum likelihood contributions of each track

Examples

uniqueID = c("Animal1", "Animal2", "Animal3")
ID = rep(uniqueID, c(100, 200, 300))
trackInd = calc_trackInd(ID)

Trigonometric basis expansion

Description

Builds a design matrix of sin/cos pairs for use in models with periodic predictors. Can be used directly or inside formulas passed to modmat (where expansion is handled automatically).

Usage

cosinor(x, period = 24)

Arguments

Details

The resulting columns form the basis for linear predictors of the form

\eta_t = \beta_0 + \sum_k \Bigl( \beta_{1k} \sin\!\Bigl(\tfrac{2 \pi x_t}{\text{period}_k}\Bigr) + \beta_{2k} \cos\!\Bigl(\tfrac{2 \pi x_t}{\text{period}_k}\Bigr) \Bigr).

Value

A numeric matrix with 2 * length(period) columns named sin(2*pi*x/period) / cos(2*pi*x/period).

Examples

cosinor(1:24, period = 24)
cosinor(1:24, period = c(24, 12, 6))

## In model formulas (expand_cosinor handles the expansion):
form <- ~ x + temp * cosinor(hour, c(24, 12))
data <- data.frame(x = runif(24), temp = rnorm(24, 20), hour = 1:24)
modmat <- modmat(form, data = data)

State dwell-time distributions of periodically inhomogeneous Markov chains

Description

Computes the dwell-time distribution of a periodically inhomogeneous Markov chain for a given transition probability matrix.

Usage

ddwell(x, Gamma, time = NULL, state = NULL)

Arguments

Details

For Markov chains whose transition probabilities vary only periodically, which is achieved for example by expressing the transition probability matrix as a periodic function of the time of day using tpm_p or cosinor, the probability distribution of time spent in a state can be computed analytically. This function computes said distribution, either for a specific time point (conditioning on transitioning into the state at that time point) or for the overall distribution (conditioning on transitioning into the state at any time point).

Value

either time-varying dwell-time distribution(s) if time is specified, or overall dwell-time distribution if time is NULL. If more than one state is specified, a named list over states is returned.

References

When using this function, please cite:

Koslik, J. O., Feldmann, C. C., Mews, S., Michels, R., & Langrock, R. (2025). Inference on the state process of periodically inhomogeneous hidden Markov models for animal behavior. The Annals of Applied Statistics, 19(4), 2724–2737.

Examples

# setting parameters for trigonometric link
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma = tpm_p(beta = beta, degree = 1)

# at specific times and for specific state
ddwell(1:20, Gamma, time = 1:4, state = 1)
# results in 4x20 matrix

# or overall distribution for all states
ddwell(1:20, Gamma)
# results in list of length 2, each element is a vector of length 20

Reparametrised multivariate Gaussian distribution

Description

Density function of the multivariate Gaussian distribution reparametrised in terms of its precision matrix (inverse variance). This implementation is particularly useful for defining the joint log-likelihood with penalised splines or i.i.d. random effects that have a multivariate Gaussian distribution with fixed precision/ penalty matrix \lambda S. As S is fixed and only scaled by \lambda, it is more efficient to precompute the determinant of S (for the normalisation constant) and only scale the quadratic form by \lambda when multiple spline parameters/ random effects with different \lambda's but the same penalty matrix S are evaluated.

Usage

dgmrf2(x, mu = 0, S, lambda, logdetS = NULL, log = FALSE)

Arguments

Details

This implementation allows for automatic differentiation with RTMB.

Value

vector of density values

Examples

x = matrix(runif(30), nrow = 3)

# iid random effects
S = diag(10)
sigma = c(1, 2, 3) # random effect standard deviations
lambda = 1 / sigma^2
d = dgmrf2(x, 0, S, lambda)

# P-splines
L = diff(diag(10), diff = 2) # second-order difference matrix
S = t(L) %*% L
lambda = c(1,2,3)
d = dgmrf2(x, 0, S, lambda, log = TRUE)

Forward algorithm to calculate the HMM log-likelihood

Description

Calculates the log-likelihood of a sequence of observations under a hidden Markov model using the forward algorithm (Zucchini, MacDonald & Langrock, 2016).

Usage

forward(
  delta,
  Gamma,
  allprobs,
  trackID = NULL,
  logspace = FALSE,
  bw = NULL,
  report = TRUE,
  ad = NULL,
  tod = NULL
)

Arguments

Details

If trackID is provided, the total log-likelihood is the sum of each track's likelihood contribution. In this case, Gamma can be

• a matrix (same transition probabilities for each track),

• an array of dimension c(nStates, nStates, nTracks) (track-specific transition probability matrix), or

• an array of dimension c(nStates, nStates, nObs) for time-varying transition probabilities, or

• an array of dimension c(nStates, nStates, L) for periodic transition probabilities, if tod is provided.

Additionally, delta can be a vector (same initial distribution for each track) or a matrix of dimension c(nTracks, nStates) (track-specific initial distributions).

With delta = "stationary", the initial distribution of each track is computed from Gamma:

• the stationary distribution of Gamma (or of its slice for each track, if track-specific),

• the periodically stationary distribution at the first time point of each track, if tod is provided (see stationary_p),

• the stationary distribution of the first slice of each track, for time-varying transition probabilities, i.e. at the covariate values of the first observation.

This requires Gamma to contain transition probability matrices (not, e.g., sub-stochastic matrices), which is checked without automatic differentiation. The computed delta is reported as a matrix of dimension c(nTracks, nStates).

Note: forward should only be called once in the likelihood function, with a trackID argument for multiple tracks. Downstream functions like lama, viterbi, stateprobs or pseudo_res rely on the quantities it reports, and only those of the last call are kept.

Value

HMM log-likelihood for given data and parameters

References

Zucchini, W., MacDonald, I.L., & Langrock, R. (2016). Hidden Markov Models for Time Series: An Introduction Using R (2nd ed.). Chapman & Hall/CRC.

Mews, Koslik, & Langrock (2025). How to build your latent Markov model: The role of time and space. Statistical Modelling, 25(6), 481–507.

See Also

Other forward algorithms: forward_g(), forward_hsmm(), forward_ihsmm(), forward_p(), forward_phsmm()

Examples

## negative log likelihood function
nll = function(par, step) {
 # parameter transformations for unconstrained optimisation
 Gamma = tpm(par[1:2]) # multinomial logit link
 delta = stationary(Gamma) # stationary HMM
 mu = exp(par[3:4])
 sigma = exp(par[5:6])
 # calculate all state-dependent probabilities
 allprobs = matrix(1, length(step), 2)
 ind = which(!is.na(step))
 for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
 # simple forward algorithm to calculate log-likelihood
 -forward(delta, Gamma, allprobs)
}

## fitting an HMM to the trex data
par = c(-2,-2,            # initial tpm params (logit-scale)
        log(c(0.3, 2.5)), # initial means for step length (log-transformed)
        log(c(0.2, 1.5))) # initial sds for step length (log-transformed)
mod = nlm(nll, par, step = trex$step[1:1000])

Forward algorithm with time-varying transition probability matrix

Description

Calculates the log-likelihood of a sequence of observations under a hidden Markov model with time-varying transition probabilities using the forward algorithm (Zucchini, MacDonald & Langrock, 2016).

Usage

forward_g(
  delta,
  Gamma,
  allprobs,
  trackID = NULL,
  logspace = FALSE,
  bw = NULL,
  report = TRUE,
  ad = NULL
)

Arguments

Details

Deprecated: please use forward for all model types instead.

If trackID is provided, the total log-likelihood will be the sum of each track's likelihood contribution. In this case, Gamma must be an array of dimension c(nStates, nStates, nObs), matching the number of rows of allprobs. For each track, the transition matrix at the beginning of the track will be ignored (as there is no transition between tracks). Additionally, delta can be a vector (same initial distribution for each track) or a matrix of dimension c(nTracks, nStates) (different initial distribution for each track).

Note: When there are multiple tracks, for compatibility with downstream functions like viterbi_g, stateprobs_g or pseudo_res, forward_g should only be called once with a trackID argument.

Value

HMM log-likelihood for given data and parameters

References

Zucchini, W., MacDonald, I.L., & Langrock, R. (2016). Hidden Markov Models for Time Series: An Introduction Using R (2nd ed.). Chapman & Hall/CRC.

Mews, Koslik, & Langrock (2025). How to build your latent Markov model: The role of time and space. Statistical Modelling, 25(6), 481–507.

See Also

Other forward algorithms: forward(), forward_hsmm(), forward_ihsmm(), forward_p(), forward_phsmm()

Examples

## Simple usage
Gamma = array(c(0.9, 0.2, 0.1, 0.8), dim = c(2,2,10))
delta = c(0.5, 0.5)
allprobs = matrix(0.5, 10, 2)
forward_g(delta, Gamma, allprobs)

## Full model fitting example
## negative log likelihood function
nll = function(par, step, Z) {
 # parameter transformations for unconstrained optimisation
 beta = matrix(par[1:6], nrow = 2)
 Gamma = tpm_g(Z, beta) # multinomial logit link for each time point
 delta = stationary(Gamma[,,1]) # stationary HMM
 mu = exp(par[7:8])
 sigma = exp(par[9:10])
 # calculate all state-dependent probabilities
 allprobs = matrix(1, length(step), 2)
 ind = which(!is.na(step))
 for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
 # simple forward algorithm to calculate log-likelihood
 -forward_g(delta, Gamma, allprobs)
}

## fitting an HMM to the trex data
par = c(-1.5,-1.5,        # initial tpm intercepts (logit-scale)
        rep(0, 4),        # initial tpm slopes
        log(c(0.3, 2.5)), # initial means for step length (log-transformed)
        log(c(0.2, 1.5))) # initial sds for step length (log-transformed)
mod = nlm(nll, par, step = trex$step[1:500], Z = cosinor(trex$tod[1:500]))

Forward algorithm for homogeneous hidden semi-Markov models

Description

Calculates the (approximate) log-likelihood of a sequence of observations under a homogeneous hidden semi-Markov model using a modified forward algorithm.

Usage

forward_hsmm(
  dm,
  omega,
  allprobs,
  trackID = NULL,
  delta = NULL,
  eps = 1e-10,
  report = TRUE
)

Arguments

Details

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs, where the state duration distribution is explicitly modelled by a distribution on the positive integers. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function is designed to be used with automatic differentiation based on the R package RTMB. It will be very slow without it!

Value

HSMM log-likelihood for given data and parameters

References

Langrock, R., & Zucchini, W. (2011). Hidden Markov models with arbitrary state dwell-time distributions. Computational Statistics & Data Analysis, 55(1), 715-724.

Koslik, J. O. (2025). Hidden semi-Markov models with inhomogeneous state dwell-time distributions. Computational Statistics & Data Analysis, 209, 108171.

See Also

Other forward algorithms: forward(), forward_g(), forward_ihsmm(), forward_p(), forward_phsmm()

Examples

# currently no examples

Forward algorithm for hidden semi-Markov models with inhomogeneous state durations and/ or conditional transition probabilities

Description

Calculates the (approximate) log-likelihood of a sequence of observations under an inhomogeneous hidden semi-Markov model using a modified forward algorithm.

Usage

forward_ihsmm(
  dm,
  omega,
  allprobs,
  trackID = NULL,
  delta = NULL,
  startInd = NULL,
  eps = 1e-10,
  report = TRUE
)

Arguments

Details

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs, where the state duration distribution is explicitly modelled by a distribution on the positive integers. This function can be used to fit HSMMs where the state-duration distribution and/ or the conditional transition probabilities vary with covariates. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function is designed to be used with automatic differentiation based on the R package RTMB. It will be very slow without it!

Value

HSMM log-likelihood for given data and parameters

References

Koslik, J. O. (2025). Hidden semi-Markov models with inhomogeneous state dwell-time distributions. Computational Statistics & Data Analysis, 209, 108171.

See Also

Other forward algorithms: forward(), forward_g(), forward_hsmm(), forward_p(), forward_phsmm()

Examples

# currently no examples

Forward algorithm with for periodically varying transition probability matrices

Description

Calculates the log-likelihood of a sequence of observations under a hidden Markov model with periodically varying transition probabilities using the forward algorithm.

Usage

forward_p(
  delta,
  Gamma,
  allprobs,
  tod,
  trackID = NULL,
  ad = NULL,
  report = TRUE,
  logspace = FALSE
)

Arguments

Details

Deprecated: please use forward with argument tod instead.

When the transition probability matrix only varies periodically (e.g. as a function of time of day), there are only L unique matrices if L is the period length (e.g. L=24 for hourly data and time-of-day variation). Thus, it is much more efficient to only calculate these L matrices and index them by a time variable (e.g. time of day or day of year) instead of calculating such a matrix for each index in the data set (which would be redundant). This function allows for that by only expecting a transition probability matrix for each time point in a period and an integer valued (1, \dots, L) time variable that maps the data index to the according time.

Value

HMM log-likelihood for given data and parameters

See Also

Other forward algorithms: forward(), forward_g(), forward_hsmm(), forward_ihsmm(), forward_phsmm()

Examples

## negative log likelihood function
nll = function(par, step, tod) {
 # parameter transformations for unconstrained optimisation
 beta = matrix(par[1:6], nrow = 2)
 Gamma = tpm_p(1:24, beta = beta) # multinomial logit link for each time point
 delta = stationary_p(Gamma, tod[1] %% 24 + 1) # stationary HMM
 mu = exp(par[7:8])
 sigma = exp(par[9:10])
 # calculate all state-dependent probabilities
 allprobs = matrix(1, length(step), 2)
 ind = which(!is.na(step))
 for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
 # simple forward algorithm to calculate log-likelihood
 -forward_p(delta, Gamma, allprobs, tod)
}

## fitting an HMM to the nessi data
par = c(-2,-2,            # initial tpm intercepts (logit-scale)
        rep(0, 4),        # initial tpm slopes
        log(c(0.3, 2.5)), # initial means for step length (log-transformed)
        log(c(0.2, 1.5))) # initial sds for step length (log-transformed)
mod = nlm(nll, par, step = trex$step[1:500], tod = trex$tod[1:500])

Forward algorithm for hidden semi-Markov models with periodically inhomogeneous state durations and/ or conditional transition probabilities

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs, where the state duration distribution is explicitly modelled by a distribution on the positive integers. This function can be used to fit HSMMs where the state-duration distribution and/ or the conditional transition probabilities vary with covariates. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function can be used to fit HSMMs where the state-duration distribution and/ or the conditional transition probabilities vary periodically. In the special case of periodic variation (as compared to arbitrary covariate influence), this version is to be preferred over forward_ihsmm because it computes the correct periodically stationary distribution and no observations are lost for the approximation.

This function is designed to be used with automatic differentiation based on the R package RTMB. It will be very slow without it!

Usage

forward_phsmm(
  dm,
  omega,
  allprobs,
  tod,
  trackID = NULL,
  delta = NULL,
  eps = 1e-10,
  report = TRUE
)

Arguments

Details

Calculates the (approximate) log-likelihood of a sequence of observations under a periodically inhomogeneous hidden semi-Markov model using a modified forward algorithm.

Value

HSMM log-likelihood for given data and parameters

References

Koslik, J. O. (2025). Hidden semi-Markov models with inhomogeneous state dwell-time distributions. Computational Statistics & Data Analysis, 209, 108171.

See Also

Other forward algorithms: forward(), forward_g(), forward_hsmm(), forward_ihsmm(), forward_p()

Examples

# currently no examples

Forward algorithm for hidden semi-Markov models with homogeneous transition probability matrix

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs that can be approximated by HMMs on an enlarged state space (of size M) and with structured transition probabilities.

Usage

forward_s(delta, Gamma, allprobs, sizes)

Arguments

Value

log-likelihood for given data and parameters

Examples

## generating data from homogeneous 2-state HSMM
mu = c(0, 6)
lambda = c(6, 12)
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# simulation
# for a 2-state HSMM the embedded chain always alternates between 1 and 2
s = rep(1:2, 100)
C = x = numeric(0)
for(t in 1:100){
  dt = rpois(1, lambda[s[t]])+1 # shifted Poisson
  C = c(C, rep(s[t], dt))
  x = c(x, rnorm(dt, mu[s[t]], 1.5)) # fixed sd 2 for both states
}

## negative log likelihood function
mllk = function(theta.star, x, sizes){
  # parameter transformations for unconstraint optimization
  omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE) # omega fixed (2-states)
  lambda = exp(theta.star[1:2]) # dwell time means
  dm = list(dpois(1:sizes[1]-1, lambda[1]), dpois(1:sizes[2]-1, lambda[2]))
  Gamma = tpm_hsmm2(omega, dm)
  delta = stationary(Gamma) # stationary
  mu = theta.star[3:4]
  sigma = exp(theta.star[5:6])
  # calculate all state-dependent probabilities
  allprobs = matrix(1, length(x), 2)
  for(j in 1:2){ allprobs[,j] = dnorm(x, mu[j], sigma[j]) }
  # return negative for minimization
  -forward_s(delta, Gamma, allprobs, sizes)
}

## fitting an HSMM to the data
theta.star = c(log(5), log(10), 1, 4, log(2), log(2))
mod = nlm(mllk, theta.star, x = x, sizes = c(20, 30), stepmax = 5)

Forward algorithm for hidden semi-Markov models with periodically varying transition probability matrices

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs that can be approximated by HMMs on an enlarged state space (of size M) and with structured transition probabilities. Recently, this inference procedure has been generalised to allow either the dwell-time distributions or the conditional transition probabilities to depend on external covariates such as the time of day. This special case is implemented here. This function allows for that, by expecting a transition probability matrix for each time point in a period, and an integer valued (1, \dots, L) time variable that maps the data index to the according time.

Usage

forward_sp(delta, Gamma, allprobs, sizes, tod)

Arguments

Value

log-likelihood for given data and parameters

Examples

## generating data from homogeneous 2-state HSMM
mu = c(0, 6)
beta = matrix(c(log(4),log(6),-0.2,0.2,-0.1,0.4), nrow=2)
# time varying mean dwell time
Lambda = exp(cbind(1, cosinor(1:24))%*%t(beta))
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# simulation
# for a 2-state HSMM the embedded chain always alternates between 1 and 2
s = rep(1:2, 100)
C = x = numeric(0)
tod = rep(1:24, 50) # time of day variable
time = 1
for(t in 1:100){
  dt = rpois(1, Lambda[tod[time], s[t]])+1 # dwell time depending on time of day
  time = time + dt
  C = c(C, rep(s[t], dt))
  x = c(x, rnorm(dt, mu[s[t]], 1.5)) # fixed sd 2 for both states
}
tod = tod[1:length(x)]

## negative log likelihood function
mllk = function(theta.star, x, sizes, tod){
  # parameter transformations for unconstraint optimization
  omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE) # omega fixed (2-states)
  mu = theta.star[1:2]
  sigma = exp(theta.star[3:4])
  beta = matrix(theta.star[5:10], nrow=2)
  # time varying mean dwell time
  Lambda = exp(cbind(1, cosinor(1:24))%*%t(beta))
  dm = list()
  for(j in 1:2){
    dm[[j]] = sapply(1:sizes[j]-1, dpois, lambda = Lambda[,j])
  }
  Gamma = tpm_phsmm2(omega, dm)
  delta = stationary_p(Gamma, tod[1])
  # calculate all state-dependent probabilities
  allprobs = matrix(1, length(x), 2)
  for(j in 1:2){ allprobs[,j] = dnorm(x, mu[j], sigma[j]) }
  # return negative for minimization
  -forward_sp(delta, Gamma, allprobs, sizes, tod)
}

## fitting an HSMM to the data
theta.star = c(1, 4, log(2), log(2), # state-dependent parameters
                 log(4), log(6), rep(0,4)) # state process parameters dm
# mod = nlm(mllk, theta.star, x = x, sizes = c(10, 15), tod = tod, stepmax = 5)

Reparametrised gamma distribution

Description

Density, distribution function, quantile function and random generation for the gamma distribution reparametrised in terms of mean and standard deviation.

Usage

dgamma2(x, mean = 1, sd = 1, log = FALSE)

pgamma2(q, mean = 1, sd = 1, lower.tail = TRUE, log.p = FALSE)

qgamma2(p, mean = 1, sd = 1, lower.tail = TRUE, log.p = FALSE)

rgamma2(n, mean = 1, sd = 1)

Arguments

Details

This implementation allows for automatic differentiation with RTMB.

Value

dgamma2 gives the density, pgamma2 gives the distribution function, qgamma2 gives the quantile function, and rgamma2 generates random deviates.

Examples

x = rgamma2(1)
d = dgamma2(x)
p = pgamma2(x)
q = qgamma2(p)

Computes generalised determinant

Description

Computes generalised determinant

Usage

gdeterminant(x, eps = NULL, log = TRUE)

Arguments

Value

generalised log-determinant of x

Build the generator matrix of a continuous-time Markov chain

Description

This function builds the infinitesimal generator matrix for a continuous-time Markov chain from an unconstrained parameter vector.

Usage

generator(
  beta,
  Z = NULL,
  Eta = NULL,
  byrow = FALSE,
  report = TRUE,
  param = NULL
)

Arguments

Details

Off-diagonal entries are calculated as \exp(\beta_i) to ensure positivity. The diagonal entries are then set to the negative row sums, which is required for generator matrices.

If a design matrix Z or a matrix of linear predictors Eta is provided, the function will automatically call generator_g to build the generator matrix based on the design matrix and coefficient matrix. In that case, the argument beta needs to be a matrix of coefficients of dimension c(nStates * (nStates-1), p+1), where the first column contains the intercepts.

Value

infinitesimal generator matrix of dimension c(nStates, nStates) or array of such matrices of dimension c(nStates, nStates, nObs) if Z or Eta is provided.

See Also

Other transition probability matrix functions: generator_g(), tpm(), tpm_ct(), tpm_emb(), tpm_emb_g(), tpm_g(), tpm_g2(), tpm_p()

Examples

# 2 states: 2 free off-diagonal elements
generator(rep(-1, 2))
# 3 states: 6 free off-diagonal elements
generator(rep(-2, 6))

Build generator matrices of a continuous-time Markov chain

Description

This function builds infinitesimal generator matrices for a continuous-time Markov chain based on a design matrix and coefficient matrix.

Usage

generator_g(Z, beta, Eta = NULL, byrow = FALSE, report = TRUE)

Arguments

Details

Deprecated: please use generator(beta, Z) for all model types instead.

Off-diagonal entries are calculated as \exp(Z \beta_i) to ensure positivity. The diagonal entries are then set to the negative row sums, which is required for generator matrices.

Value

array of infinitesimal generator matrices of dimension c(nStates, nStates, nObs)

See Also

Other transition probability matrix functions: generator(), tpm(), tpm_ct(), tpm_emb(), tpm_emb_g(), tpm_g(), tpm_g2(), tpm_p()

Examples

# 2 states: 2 free off-diagonal elements
generator(rep(-1, 2))
# 3 states: 6 free off-diagonal elements
generator(rep(-2, 6))

Build a LaMaModel from an optimised RTMB object

Description

Having fitted a latent Markov model using automatic differentiation via RTMB, this function calls RTMB's obj$report() and does some additional processing. This then yields estimated parameters on their natural scale, allows for convenient calculation of AIC and BIC, state-decoding, and residual calculation. print() gives a short overview of the model, summary() the estimates. report() is the deprecated old name of lama().

Usage

lama(obj)

report(obj)

Arguments

Details

The number of observations (used by BIC) is the number of rows of the reported allprobs that are not all 1, i.e. missing observations are not counted. This assumes that forward is called only once in the likelihood function.

Value

A model object of class "LaMaModel" containing a list with the reported quantities from the RTMB object, along estimated parameters (par), the log-likelihood (llk), number of parameters (df) and observations (nobs), and the maximum gradient component (mgc) at the estimate. It also contains obj itself (not a copy), so that functions like MCreport or scores can be called on the model object.

Examples

data <- trex[1:200,]

# initial parameters and observations
par <- list(
  log_mu = log(c(0.3, 1)),      # initial means for step length (log-transformed)
  log_sigma = log(c(0.2, 0.7)), # initial sds for step length (log-transformed)
  eta = rep(-2, 2)              # initial t.p.m. parameters (on logit scale)
)    
dat <- list(
  step = data$step,             # hourly step lengths
  nStates = 2                   # number of hidden states
)

# likelihood function
nll <- function(par) {
  getAll(par, dat)
  Gamma <- tpm(eta)
  delta <- stationary(Gamma)
  mu <- exp(log_mu); REPORT(mu)
  sigma <- exp(log_sigma); REPORT(sigma)
  allprobs <- matrix(1, length(step), nStates)
  ind <- which(!is.na(step))
  for(j in 1:nStates) {
    allprobs[ind,j] <- dgamma2(step[ind], mu[j], sigma[j]) 
  } 
  -forward(delta, Gamma, allprobs)
}

# automatic differentiation and optimisation
obj <- MakeADFun(nll, par, silent = TRUE)
opt <- nlminb(obj$par, obj$fn, obj$gr)

### reporting ###
mod <- lama(obj)
mod          # short overview
summary(mod) # estimated state process and reported quantities

# estimated parameters
mod$par

# estimated quantities on natural scale
mod$mu
mod$sigma
mod$Gamma

# information criteria
AIC(mod)
BIC(mod)

# state decoding
states <- viterbi(mod = mod)   # global decoding
probs <- stateprobs(mod = mod) # local decoding

# residual calculation
pres <- pseudo_res(data$step, # observation sequence
                   "gamma2",  # distribution family
                   list(mean = mod$mu, sd = mod$sigma), # parameters for that family
                   mod = mod) # model object

# or from the model object alone, and plotted
plot(mod)

Extract log-likelihood from LaMaModel object

Description

Extract log-likelihood from LaMaModel object

Usage

## S3 method for class 'LaMaModel'
logLik(object, ...)

Arguments

Value

An object of class "logLik"

Build a standardised P-Spline design matrix and the associated P-Spline penalty matrix

Description

This function builds the B-spline design matrix for a given data vector. Importantly, the B-spline basis functions are normalised such that the integral of each basis function is 1, hence this basis can be used for spline-based density estimation, when the basis functions are weighted by non-negative weights summing to one.

Usage

make_matrices_dens(
  x,
  k,
  type = "real",
  degree = 3,
  knots = NULL,
  diff_order = 2,
  pow = 0.5,
  npoints = 10000
)

Arguments

Value

list containing the design matrix Z, the penalty matrix S, the prediction design matrix Z_predict, the prediction grid xseq, and details for the basis expansion.

The penalty matrix leaves out the last coefficient, which is fixed at zero for identifiability.

Examples

set.seed(1)
# real-valued
x <- rnorm(100)
modmat <- make_matrices_dens(x, k = 20)
# positive-continuouos
x <- rgamma2(100, mean = 5, sd = 2)
modmat <- make_matrices_dens(x, k = 20, type = "positive")
# circular
x <- rvm(100, mu = 0, kappa = 2)
modmat <- make_matrices_dens(x, k = 20, type = "circular")
# bounded in an interval
x <- rbeta(100, 1, 2)
modmat <- make_matrices_dens(x, k = 20)

AD-compatible minimum and maximum functions

Description

These functions compute the parallel minimum/ maximum of two vector-valued inputs and are compatible with automatic differentiation using RTMB.

Usage

min2(x, y)

max2(x, y)

Arguments

Value

min2 returns the parallel minimum and max2 the parallel maximum of x and y

Examples

x <- c(1, 4, 8, 2)
y <- c(2, 5, 3, 7)
min2(x, y)
max2(x, y)

Smooth approximations to max(x, 0) and min(x, 0)

Description

Smooth approximations to max(x, 0) and min(x, 0)

Usage

max0_smooth(x, rho = 20)

min0_smooth(x, rho = 20)

Arguments

Value

the approximate maximum or minimum of x and 0

Examples

x <- seq(-1, 1, by = 0.1)
min0_smooth(x)
max0_smooth(x)

Build the design and the penalty matrix for models involving penalised splines based on a formula and a data set

Description

Build the design and the penalty matrix for models involving penalised splines based on a formula and a data set

Usage

modmat(formula, data, knots = NULL)

make_matrices(formula, data, knots = NULL)

Arguments

Details

make_matrices() is the former name of modmat() and still works.

Value

a list of class LaMa_matrices containing:

See Also

predict.LaMa_matrices for prediction design matrix construction based on the model matrices object created by this function.

Examples

data = data.frame(x = runif(100), 
                  y = runif(100),
                  g = factor(rep(1:10, each = 10)))

# unvariate thin plate regression spline
modmat = modmat(~ s(x), data)
# univariate P-spline
modmat = modmat(~ s(x, bs = "ps"), data)
# adding random intercept
modmat = modmat(~ s(g, bs = "re") + s(x, bs = "ps"), data)
# tensorproduct of x and y
modmat = modmat(~ s(x) + s(y) + ti(x,y), data)
# multiple formulas at once
modmat = modmat(list(mu ~ s(x) + y, sigma ~ s(g, bs = "re")), data = data)

Optimise an RTMB object from several starting values

Description

Latent Markov models often have several local optima of the likelihood. multistart() optimises obj from several starting values: the first is obj$par (the initial parameters), the others are random perturbations of it on the working scale. Afterwards, lama(obj) builds the model from the best fit, as RTMB keeps the parameters with the lowest negative log-likelihood it evaluated.

Usage

multistart(
  obj,
  nstarts = 20,
  sd = 1,
  starts = NULL,
  scale = c("fixed", "unit_info"),
  stop = c("bayes", "hits", "none"),
  hits = 3,
  sd_max = 2,
  tol = 1e-06,
  seed = NULL,
  optimizer = stats::nlminb,
  silent = FALSE,
  ...
)

Arguments

Details

Stopping. Two runs reach the same optimum if their negative log-likelihoods differ by less than tol (relative), so permuted states (label switching) count as the same optimum. Only converged runs (convergence code 0 and maximum gradient component below 0.01) count for stopping. In any case, at most nstarts starts are run.

• stop = "bayes" uses the Bayesian stopping rule of Boender and Rinnooy Kan (1987): after k runs with w distinct optima, the expected number of optima is w(k-1)/(k-w-2), and multistart() stops when this is below w + 0.5, i.e. probably no optimum is missed. The rule assumes that the starts are independent draws from one distribution, so only the random starts count, not the first run from obj$par, and it is not used for given starts. It needs at least 8 converged random starts.

• stop = "hits" stops once the best optimum was reached hits times.

Caution. Stopping does not mean that the global optimum was found. "No optimum missed" only refers to the distribution of the starts: optima far away from it, or with a small basin of attraction, can be missed. For example, with a very small sd, all starts converge to the same optimum and the rule stops quickly, although other starts might find a better one. So choose sd wide enough and still check the fit.

Scatter.

• With scale = "fixed", the starts are obj$par plus normal noise with standard deviation sd. As most working parameters are on log or logit scale, sd = 1 is a relative change, so one value usually fits all parameters. Parameters on identity scale (e.g. means of normal distributions) may need their own sd, given as a named list like list(eta = 1, mu = 5).

• With scale = "unit_info", the starts are drawn around the first fit from a normal distribution with covariance n H^{-1}, where H is the Hessian of the negative log-likelihood and n the number of observations. This is the information of a single observation (unit information), which does not shrink with n and takes correlations of the parameters into account. Standard deviations are capped at sd_max, for weakly identified parameters. If the Hessian is not positive definite (or allprobs is not reported, so n is unknown), the fixed sd is used instead.

Value

output of optimizer for the best run, with two more elements: runs, a data frame with the negative log-likelihood nll, convergence code and maximum gradient component mgc of each run, and starts, the matrix of starting values used.

References

Boender, C. G. E., & Rinnooy Kan, A. H. G. (1987). Bayesian stopping rules for multistart global optimization methods. Mathematical Programming, 37(1), 59-80.

Examples

nll <- function(par) {
  getAll(par, dat)
  Gamma <- tpm(eta)
  mu <- exp(log_mu); REPORT(mu)
  allprobs <- cbind(dgamma2(step, mu[1], exp(log_sigma[1])),
                   dgamma2(step, mu[2], exp(log_sigma[2])))
  -forward(stationary(Gamma), Gamma, allprobs)
}
dat <- list(step = trex$step[1:500])
par <- list(eta = c(-2, -2), log_mu = log(c(0.3, 2.5)), log_sigma = log(c(0.2, 1.5)))
obj <- MakeADFun(nll, par, silent = TRUE)

opt <- multistart(obj, seed = 1)
opt$runs
mod <- lama(obj) # best fit

# starts scaled by the first fit
opt <- multistart(obj, scale = "unit_info", seed = 1)

Loch Ness Monster Acceleration Data

Description

A small group of researchers managed to put an accelerometer on the Loch Ness Monster and collected data for a few days. Now we have a data set of the overall dynamic body acceleration (ODBA) of the creature.

Usage

nessi

Format

A data frame with 5.000 rows and 3 variables:

ODBA

overall dynamci body acceleration

logODBA

logarithm of overall dynamic body acceleration

state

hidden state variable

Source

Generated for example purposes.

Computes penalty based on quadratic form

Description

This function computes quadratic penalties of the form

0.5 \sum_{i} \lambda_i b_i^T S_i b_i,

with smoothing parameters \lambda_i, coefficient vectors b_i, and fixed penalty matrices S_i.

It is intended to be used inside the penalised negative log-likelihood function when fitting models with penalised splines or simple random effects via quasi restricted maximum likelihood (qREML) with the qreml function. For qreml to work, the likelihood function needs to be compatible with the RTMB R package to enable automatic differentiation.

Usage

penalty(re_coef, S, sp, lambda)

Arguments

Details

Caution: The formatting of re_coef needs to match the structure of the parameter list in your penalised negative log-likelihood function, i.e. you cannot have two random effect vectors of different names (different list elements in the parameter list), combine them into a matrix inside your likelihood and pass the matrix to penalty. If these are seperate random effects, each with its own name, they need to be passed as a list to penalty. Moreover, the ordering of re_coef needs to match the character vector random specified in qreml.

Value

returns the penalty value and reports to qreml.

References

Koslik, J. O. (2024). Efficient smoothness selection for nonparametric Markov-switching models via quasi restricted maximum likelihood. arXiv preprint arXiv:2411.11498.

See Also

qreml for the qREML algorithm

Examples

# Example with a single random effect
re = rep(0, 5)
S = diag(5)
sp = 1
penalty(re, S, sp)

# Example with two random effects, 
# where one element contains two random effects of similar structure
re = list(matrix(0, 2, 5), rep(0, 4))
S = list(diag(5), diag(4))
sp = c(1,1,2) # length = total number of random effects
penalty(re, S, sp)

# Full model-fitting example

data = trex[1:1000,] # subset

# initial parameter list
par = list(logmu = log(c(0.3, 2.5)), # step mean
           logsigma = log(c(0.3, 1.5)), # step sd
           beta0 = c(-2,-2), # state process intercept
           betaspline = matrix(rep(0, 18), nrow = 2)) # state process spline coefs
          
# data object with initial penalty strength sp
dat = list(step = data$step, # step length
           tod = data$tod, # time of day covariate
           N = 2, # number of states
           sp = rep(10,2)) # initial penalty strength

# building model matrices
modmat = modmat(~ s(tod, bs = "cp"), 
                       data = data.frame(tod = 1:24), 
                       knots = list(tod = c(0,24))) # wrapping points
dat$Z = modmat$Z # spline design matrix
dat$S = modmat$S # penalty matrix

# penalised negative log-likelihood function
pnll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm_g(Z, cbind(beta0, betaspline)) # transition probabilities
  delta = stationary_p(Gamma, t = tod[1] %% 24 + 1) # initial distribution
  mu = exp(logmu) # step mean
  sigma = exp(logsigma) # step sd
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step)) # only for non-NA obs.
  for(j in 1:N) allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])
  -forward_g(delta, Gamma[,,tod], allprobs) +
      penalty(betaspline, S, sp) # this does all the penalization work
}

# model fitting
mod = qreml(pnll, par, dat, random = "betaspline")

Computes generalised quadratic-form penalties

Description

This function computes a quadratic penalty of the form

0.5 \sum_{i} \lambda_i b^T S_i b,

with smoothing parameters \lambda_i, coefficient vector b, and fixed penalty matrices S_i. This generalises the penalty by allowing subsets of the coefficient vector b to be penalised multiple times with different smoothing parameters, which is necessary for tensor products, functional random effects or adaptive smoothing.

It is intended to be used inside the penalised negative log-likelihood function when fitting models with penalised splines or simple random effects via quasi restricted maximum likelihood (qREML) with the qreml function. For qreml to work, the likelihood function needs to be compatible with the RTMB R package to enable automatic differentiation.

Usage

penalty2(re_coef, S, sp, lambda)

Arguments

Details

Caution: The formatting of re_coef needs to match the structure of the parameter list in your penalised negative log-likelihood function, i.e. you cannot have two random effect vectors of different names (different list elements in the parameter list), combine them into a matrix inside your likelihood and pass the matrix to penalty. If these are seperate random effects, each with its own name, they need to be passed as a list to penalty. Moreover, the ordering of re_coef needs to match the character vector random specified in qreml.

Value

returns the penalty value and reports to qreml.

See Also

qreml for the qREML algorithm

Examples

# Example with a single random effect
re = rep(0, 5)
S = diag(5)
sp = 1
penalty(re, S, sp)

# Example with two random effects, 
# where one element contains two random effects of similar structure
re = list(matrix(0, 2, 5), rep(0, 4))
S = list(diag(5), diag(4))
sp = c(1,1,2) # length = total number of random effects
penalty(re, S, sp)

# Full model-fitting example

data = trex[1:1000,] # subset

# initial parameter list
par = list(logmu = log(c(0.3, 2.5)), # step mean
           logsigma = log(c(0.3, 1.5)), # step sd
           beta0 = c(-2,-2), # state process intercept
           betaspline = matrix(rep(0, 18), nrow = 2)) # state process spline coefs
          
# data object with initial penalty strength sp
dat = list(step = data$step, # step length
           tod = data$tod, # time of day covariate
           N = 2, # number of states
           sp = rep(10,2)) # initial penalty strength

# building model matrices
modmat = modmat(~ s(tod, bs = "cp"), 
                       data = data.frame(tod = 1:24), 
                       knots = list(tod = c(0,24))) # wrapping points
dat$Z = modmat$Z # spline design matrix
dat$S = modmat$S # penalty matrix

# penalised negative log-likelihood function
pnll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm_g(Z, cbind(beta0, betaspline)) # transition probabilities
  delta = stationary_p(Gamma, t = tod[1] %% 24 + 1) # initial distribution
  mu = exp(logmu) # step mean
  sigma = exp(logsigma) # step sd
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step)) # only for non-NA obs.
  for(j in 1:N) allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])
  -forward_g(delta, Gamma[,,tod], allprobs) +
      penalty(betaspline, S, sp) # this does all the penalisation work
}

# model fitting
mod = qreml(pnll, par, dat, random = "betaspline")

Penalty approximation of unimodality constraints for univariates smooths

Description

Penalty approximation of unimodality constraints for univariates smooths

Usage

penalty_uni(coef, m, kappa = 1000, concave = TRUE, rho = 20)

Arguments

Value

a numeric value of the penalty for the given coefficients

Examples

## coefficient vector
coef <- c(1, 2, 3, 2, 1)
# mode at position 3
penalty_uni(coef, m = 3) # basically zero
#' # mode at position 2
penalty_uni(coef, m = 2) # large positive penalty

## coefficient matrix
coef <- rbind(coef, coef)
m <- c(1, 4)
penalty_uni(coef, m)

Plot a LaMaModel object

Description

Plots the pseudo-residuals of a model from lama or qreml, computed by pseudo_res from the model object alone, with one row of plots per data stream, or the estimated state-dependent densities of each data stream built with statedens.

Usage

## S3 method for class 'LaMaModel'
plot(x, type = c("residuals", "statedens"), stream = NULL, ...)

Arguments

Details

For type = "statedens", each data stream is shown as a histogram with the state-dependent densities (probabilities for discrete distributions) weighted by the state probabilities, and their sum, the marginal density, as a dashed line. The weighted density of state j is the mean over the observations of p_{tj} f_j(x; \theta_{tj}), which allows for parameters that depend on covariates, several tracks and inhomogeneous state processes. This needs allprobs to be the product of matrices from statedens, which is checked.

Further arguments for type = "statedens":

weights

state probabilities p_{tj}: "stateprobs" (default) for the local state probabilities from stateprobs, "viterbi" for the states decoded by viterbi, "stationary" for the stationary distribution (averaged over tracks, or over the cycle for periodic models), or a vector of length nStates

col

colours of the states, defaults to LaMaColors

lwd

line width

n

number of grid points for continuous distributions

legend

position of the legend, e.g. "topright" (default), or FALSE for none

...

arguments passed to plot() of the histogram, e.g. breaks, main, xlab, xlim, ylim or border

Several streams are shown side by side, unless par(mfrow) has been set.

Value

for type = "residuals" the pseudo-residuals, for type = "statedens" a list with the grid x, the weighted state-dependent densities density (one column per state), the state weights and the observations data (or a list of these for several streams), invisibly

Examples

dat <- list(step = trex$step[1:500], angle = trex$angle[1:500])
nll <- function(par) {
  getAll(par, dat)
  Gamma <- tpm(eta)
  allprobs <- statedens(step, dgamma2, exp(log_mu), exp(log_sigma)) *
    statedens(angle, dvm, 0, exp(log_kappa))
  -forward(stationary(Gamma), Gamma, allprobs)
}
par <- list(eta = c(-2, -2), log_mu = log(c(0.3, 2.5)), log_sigma = log(c(0.2, 1.5)),
            log_kappa = log(c(0.2, 1.5)))
obj <- MakeADFun(nll, par, silent = TRUE)
opt <- nlminb(obj$par, obj$fn, obj$gr)
mod <- lama(obj)

plot(mod) # pseudo-residuals
plot(mod, type = "statedens") # state-dependent densities of both streams
plot(mod, type = "statedens", stream = "step", weights = "viterbi", 
     breaks = 40, xlab = "step length")

Plot pseudo-residuals

Description

Plot pseudo-residuals computed by pseudo_res. The QQ plot shows a pointwise 95% band (in gray) for the sample quantiles under a correct model.

Usage

## S3 method for class 'LaMaResiduals'
plot(
  x,
  col = "darkblue",
  lwd = 1.5,
  main = NULL,
  breaks = "Sturges",
  axis.lab = list(qq = c("Theoretical quantiles", "Sample quantiles"), hist =
    c("Pseudo-residuals", "Density"), acf = c("Lag", "ACF")),
  ...
)

Arguments

Value

NULL, plots the pseudo-residuals in a 3-panel layout

Examples

## pseudo-residuals for the trex data
step = trex$step[1:1000]
angle = trex$angle[1:1000]

nll = function(par){
  getAll(par)
  Gamma = tpm(logitGamma)
  delta = stationary(Gamma)
  mu = exp(logMu); REPORT(mu)
  sigma = exp(logSigma); REPORT(sigma)
  kappa = exp(logKappa); REPORT(kappa)
  allprobs = matrix(1, length(step), 2)
  ind = which(!is.na(step) & !is.na(angle))
  for(j in 1:2) {
    allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j]) * 
         dvm(angle[ind], 0, kappa[j])
  }
  -forward(delta, Gamma, allprobs)
}

par = list(logitGamma = c(-2,-2), 
           logMu = log(c(0.3, 2.5)), 
           logSigma = log(c(0.3, 1.5)),
           logKappa = log(c(0.1, 1.7)))
           
obj = MakeADFun(nll, par, silent = TRUE)
opt = nlminb(obj$par, obj$fn, obj$gr)

mod = lama(obj)

pres_step = pseudo_res(step,      # observations
                       "gamma2",  # family that is used
                       list(mean = mod$mu, sd = mod$sigma), # the family's parameters
                       mod = mod) # model object
pres_angle = pseudo_res(angle,
                        "vm",
                        list(mu = 0, kappa = mod$kappa),
                        mod = mod)
                        
# separate plots 
plot(pres_step)
plot(pres_angle)

# together
par(mfrow = c(2,3))
plot(pres_step, main = c("", "Step Length", ""), 
     axis.lab = list(hist = c("Step residuals", "Density")))
plot(pres_angle, main = c("", "Turning Angle", ""),
     axis.lab = list(hist = c("Angle residuals", "Density")))

# or from the model object alone, one row per stream
plot(pseudo_res(mod))

Build the prediction design matrix based on new data and model_matrices object created by modmat

Description

Deprecated. Use predict on the object returned by modmat instead, which does the same thing: predict(modmat, newdata) in place of pred_matrix(modmat, newdata).

Usage

pred_matrix(model_matrices, newdata, what = NULL, exclude = NULL)

Arguments

Value

prediction design matrix for newdata with the same basis as used for model_matrices

See Also

predict, which replaces this function

Examples

modmat = modmat(~ s(x), data.frame(x = 1:10))
# deprecated:
# Z_p = pred_matrix(modmat, data.frame(x = 1:10 - 0.5))
# use instead:
Z_p = predict(modmat, data.frame(x = 1:10 - 0.5))

Build the prediction design matrix based on new data and model_matrices object created by modmat

Description

Build the prediction design matrix based on new data and model_matrices object created by modmat

Usage

## S3 method for class 'LaMa_matrices'
predict(object, newdata, what = NULL, exclude = NULL, ...)

Arguments

Value

prediction design matrix for newdata with the same basis as used for object

See Also

modmat for creating objects of class LaMa_matrices which can be used for prediction by this function.

Examples

# single formula
modmat = modmat(~ s(x), data.frame(x = 1:10))
Z_p = predict(modmat, data.frame(x = 1:10 - 0.5))
# with multiple formulas
modmat = modmat(list(mu ~ s(x), sigma ~ s(x, bs = "ps")), data = data.frame(x = 1:10))
Z_p = predict(modmat, data.frame(x = 1:10 - 0.5), what = "mu")
# nested formula list
form = list(stream1 = list(mu ~ s(x), sigma ~ s(x, bs = "ps")))
modmat = modmat(form, data = data.frame(x = 1:10))
Z_p = predict(modmat, data.frame(x = 1:10 - 0.5), what = c("stream1", "mu"))

Estimated smooth densities

Description

Evaluates the smooth densities built with smoothdens for given spline coefficients, e.g. to plot the estimated state-dependent densities.

Usage

## S3 method for class 'smoothdens'
predict(object, coef, x = NULL, ...)

Arguments

Value

list with the points x and the matrix density of densities, one row per point and one column per state

Examples

x = c(rnorm(100, 0, 1), rnorm(100, 5, 1))
dens = smoothdens(x, par = list(mean = c(0, 5), sd = c(1, 1)))
# densities at the initial coefficients
pred = predict(dens, coef = dens$coef)
matplot(pred$x, pred$density, type = "l", lty = 1)

Print a LaMaModel object

Description

Prints a short overview of a model from lama or qreml: the state process, observations, parameters, log-likelihood, AIC and BIC, and the names of the reported quantities. Use summary() for the estimates.

Usage

## S3 method for class 'LaMaModel'
print(x, ...)

Arguments

Value

x, invisibly

Print a LaMa_matrices object

Description

Print a LaMa_matrices object

Usage

## S3 method for class 'LaMa_matrices'
print(x, ...)

Arguments

Value

x, invisibly. Called for its side effect of printing.

See Also

modmat

Process and standardise formulas for the state process of hidden Markov models

Description

Process and standardise formulas for the state process of hidden Markov models

Usage

process_hid_formulas(formulas, nStates, ref = NULL)

Arguments

Value

named list of formulas of length nStates * (nStates - 1), where each formula corresponds to a transition from state i to state j, excluding transitions from i to ref[i].

Examples

# single formula for all non-reference category elements
formulas = process_hid_formulas(~ s(x), nStates = 3)
# now a list of length 6 with names tr.ij, not including reference categories

# different reference categories
formulas = process_hid_formulas(~ s(x), nStates = 3, ref = c(1,1,1))

# different formulas for different entries (and only for 2 of 6)
formulas = list(tr.12 ~ s(x), tr.23 ~ s(y))
formulas = process_hid_formulas(formulas, nStates = 3, ref = c(1,1,1))
# also a list of length 6, remaining entries filled with tr.ij ~ 1

# matrix input with reference categories
formulas = matrix(c(".", "~ s(x)", "~ s(y)",
                    "~ g", ".", "~ I(x^2)",
                    "~ y", "~ 1", "."), 
                    nrow = 3, byrow = TRUE)
# dots define reference categories
formulas = process_hid_formulas(formulas, nStates = 3)

Calculate pseudo-residuals

Description

For HMMs, pseudo-residuals are used to assess the overall goodness-of-fit of the model. These are based on the cumulative distribution function (CDF)

F_{X_t}(x_t) = F(x_t \mid x_1, \dots, x_{t-1})

and can be used to quantify whether an observation is extreme relative to its model-implied distribution.

This function calculates such residuals via the probability integral transform, based on the local state probabilities obtained by stateprobs and the respective parametric family.

Usage

pseudo_res(
  obs,
  dist,
  par,
  stateprobs = NULL,
  mod = NULL,
  normal = TRUE,
  discrete = NULL,
  randomise = TRUE
)

Arguments

Details

Pseudo-residuals are based on the probability integral transform. If the cumulative distribution function (CDF) F of a random variable X is known, then F(X) is uniformly distributed on [0,1]. Applying the standard normal quantile function \Phi^{-1} then gives a residual that is approximately standard normal under the model.

For discrete-valued observations, the CDF is a step function, so the residual is not uniquely defined. One can compute a "lower" and "upper" pseudo-residual corresponding to F(X-1) and F(X), respectively. If randomise = TRUE, a value is drawn uniformly between the lower and upper bounds. In this case, the resulting pseudo-residuals are again approximately standard normal after applying the standard normal quantile function.

Model object only: pseudo_res(mod) needs no observations, family or parameters. It evaluates the likelihood function stored in the model again, with each density (e.g. dgamma2()) replaced by its CDF (pgamma2()) with the same arguments, so that the likelihood function itself builds the CDFs in place of allprobs. This works for loops over states or time, cbind(), covariate-dependent parameters, densities on the log scale, mixtures within a state and several data streams (multiplied, or summed on the log scale), giving one set of residuals per stream, named after the data variable. Before, it checks that the densities enter allprobs as a product or mixture, so that swapping in CDFs is valid. If not, e.g. for densities transformed by hand (truncation, Jacobians, zero inflation via indicators) or a latent variable integrated out inside allprobs, it stops, and pseudo_res(obs, dist, par, mod = mod) needs to be used instead. Densities in your own helper functions called by the likelihood function are found as well, and qreml models work like models from lama. The likelihood function needs to call forward once (use trackID for several tracks), and it is evaluated with the data as they are when pseudo_res() is called.

Value

vector of pseudo residuals of class LaMaResiduals or list containing lower, upper, and mean if discrete residuals are not randomised. For pseudo_res(mod) with several data streams, a named list of these, one per stream, also of class LaMaResiduals.

See Also

plot.LaMaResiduals for plotting pseudo-residuals.

Examples

## continuous-valued observations
obs = rnorm(100)
stateprobs = matrix(0.5, nrow = 100, ncol = 2)
par = list(mean = c(1,2), sd = c(1,1))
pres = pseudo_res(obs, "norm", par, stateprobs)

## discrete-valued observations
obs = rpois(100, lambda = 1)
par = list(lambda = c(1,2))
pres = pseudo_res(obs, "pois", par, stateprobs)

## custom CDF function
obs = rnbinom(100, size = 1, prob = 0.5)
par = list(size = c(0.5, 2), prob = c(0.4, 0.6))
pres = pseudo_res(obs, pnbinom, par, stateprobs, 
                  discrete = TRUE)
# if discrete CDF function is passed, 'discrete' needs to be set to TRUE

## no CDF function available, only density (artificial example)
obs = rnorm(100)
par = list(mean = c(1,2), sd = c(1,1))
# construct CDF using numerical integration
cdf = function(x, mean, sd) integrate(dnorm, -Inf, x, mean = mean, sd = sd)$value
pres = pseudo_res(obs, cdf, par, stateprobs)


### Full model fit example ###
step = trex$step[1:500]

nll = function(par){
  getAll(par)
  Gamma = tpm(logitGamma)
  delta = stationary(Gamma)
  mu = exp(logMu); REPORT(mu)
  sigma = exp(logSigma); REPORT(sigma)
  allprobs = matrix(1, length(step), 2)
  ind = which(!is.na(step))
  for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
  -forward(delta, Gamma, allprobs)
}

par = list(logitGamma = c(-2,-2), 
           logMu = log(c(0.3, 2.5)), 
           logSigma = log(c(0.3, 1.5)))
           
obj = MakeADFun(nll, par, silent = TRUE)
opt = nlminb(obj$par, obj$fn, obj$gr)

mod = lama(obj)

pres = pseudo_res(step,      # observation sequence
                  "gamma2",  # parametric family that was used
                  list(mean = mod$mu, sd = mod$sigma), # parameters for that family
                  mod = mod) # model object

# or, from the model object alone
pres = pseudo_res(mod)

plot(pres)

Automatic smoothness selection for arbitrary statistical models with penalised splines, simple random effects, or tensorproducts

Description

Efficiently fits models involving quadratic penalties of the form

\sum_{i} \lambda_i b^T S_i b,

by combining automatic differentiation via RTMB with a custom implementation of the extended Fellner-Schall update.

Users only need to supply a custom penalised log-likelihood function that calls penalty() or penalty2().

Usage

qreml(
  pnll,
  par,
  dat,
  random,
  map = NULL,
  silent = 1,
  spname = NULL,
  tol_edf = 0.01,
  maxiter = 100,
  alpha = 0.1,
  tol = 0.01,
  smoothing = 1,
  method = "BFGS",
  lsp_max = 15,
  step_small = 0.05,
  max_halve = 6,
  control = list(),
  spHess = FALSE,
  joint_unc = FALSE,
  saveall = FALSE
)

Arguments

Details

Step size control follows the extended Fellner-Schall implementation in mgcv while convergence is judged on the effective degrees of freedom (differing from mgcv).

The criterion. The penalty strengths are chosen to maximise the Laplace-approximate restricted log-likelihood

V(\lambda) = \ell_p(\hat{b}; \lambda) + \tfrac{1}{2} \log |S_\lambda|_+ - \tfrac{1}{2} \log |J|,

where \ell_p is the penalised log-likelihood at its mode \hat{b}, S_\lambda = \sum_i \lambda_i S_i, |\cdot|_+ denotes the product of the positive eigenvalues, and J is the Hessian of the penalised negative log-likelihood at \hat{b}. Every outer iteration refits the inner problem to obtain \hat{b} and J, so one call is a sequence of complete model fits, each warm started at the previous one.

The update. The penalty strengths are updated multiplicatively, \lambda_i \leftarrow \lambda_i r_i, by the extended Fellner-Schall ratio

r_i = \frac{\mathrm{tr}(S_\lambda^- S_i) - \mathrm{tr}(J^{-1} S_i)}{\hat{b}^T S_i \hat{b}},

which equals one exactly where \partial V / \partial \lambda_i = 0 and otherwise moves \lambda_i in the direction of the gradient, without ever making it negative. For a smooth carrying a single penalty the first trace is \mathrm{rank}(S_i) / \lambda_i, obtained once from an eigendecomposition instead of at every iteration. Only the diagonal blocks of J^{-1} enter, and these come from the Cholesky factor by triangular solves rather than from a full inverse.

Step control. The step is taken in \log \lambda and scaled by a multiplier that persists across iterations: it doubles while the criterion keeps improving and is halved, up to max_halve times, while it does not. A step that still worsens the criterion after that many halvings is accepted, which is the only way the iteration can move downhill. In addition no penalty strength may fall by more than a factor alpha in one iteration, and \log \lambda is capped at lsp_max. The fit returned is the best iterate seen, not necessarily the last one.

Convergence. The iteration stops once no smooth's effective degrees of freedom

\mathrm{edf}_i = K_i - \mathrm{tr}(J^{-1} \lambda_i S_i)

has moved by more than tol_edf over the last four iterations, provided the step in \log \lambda is also small; see tol_edf. Being eigenvalue sums, the effective degrees of freedom must lie in [0, K_i], and they leave that range only if the data Hessian is indefinite, in which case the traces driving the update are not the quantities the method assumes. That is reported rather than repaired, through the returned edf_valid and a warning, and the criterion then falls back to tol on the restricted log-likelihood.

Numerics. J is the exact automatic differentiation Hessian, which matters because the finite differencing error of optimHess is of the same order as the criterion differences the convergence test has to resolve. It is factorised by a pivoted Cholesky decomposition after scaling to a unit diagonal; should it not be positive definite, a ridge is added and reported through the returned hessian_repaired.

Uncertainty. The returned outer_hessian() gives the Hessian of the negative restricted log-likelihood in the penalty strengths, which sdreport_outer turns into standard errors. Note that these describe the curvature of V only, and treat the approximation to it as exact.

Value

model object of class "qremlModel", inheriting from "LaMaModel" (see lama), carrying the fitted quantities, the smoothness selection diagnostics shown by summary.qremlModel, and outer_hessian(), the outer Hessian that sdreport_outer turns into standard errors

References

Wood, S. N., & Fasiolo, M. (2017). A generalized Fellner-Schall method for smoothing parameter optimization with application to Tweedie location, scale and shape models. Biometrics, 73(4), 1071-1081.

Koslik, J. O. (2024). Efficient smoothness selection for nonparametric Markov-switching models via quasi restricted maximum likelihood. arXiv preprint arXiv:2411.11498.

See Also

penalty and penalty2 to compute the penalty inside the likelihood function, and qreml_old for the original algorithm

Examples

data = trex[1:2000,] # subset of the data

# initial parameter list
par = list(logmu = log(c(0.3, 2.5)), # step mean
           logsigma = log(c(0.3, 1.5)), # step sd
           beta0 = c(-2,-2), # state process intercept
           beta_spline = matrix(rep(0, 18), nrow = 2)) # state process spline coefs
          
# data object with initial penalty strength sp
dat = list(step = data$step, # step length
           tod = data$tod, # time of day covariate
           N = 2, # number of states
           sp = rep(20,2)) # initial penalty strength

# building model matrices
modmat = modmat(~ s(tod, bs = "cp"), 
                       data = data.frame(tod = 1:24), 
                       knots = list(tod = c(0,24))) # wrapping points
dat$Z = modmat$Z # spline design matrix
dat$S = modmat$S # penalty matrix

# penalised negative log-likelihood function
pnll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm(cbind(beta0, beta_spline), Z) # transition probabilities
  delta = stationary_p(Gamma, t = tod[1] %% 24 + 1) # initial distribution
  mu = exp(logmu) # step mean
  sigma = exp(logsigma) # step sd
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step)) # only for non-NA obs.
  for(j in 1:N) allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])
  -forward(delta, Gamma[,,tod], allprobs) +
      penalty(beta_spline, S, sp) # this does all the penalization work
}

# model fitting
mod = qreml(pnll, par, dat, random = "beta_spline", silent = 2)

Quasi restricted maximum likelihood (qREML) algorithm, previous version

Description

This algorithm can be used very flexibly to fit statistical models that involve penalised splines or simple i.i.d. random effects, i.e. that have penalties of the form

0.5 \sum_{i} \lambda_i b_i^T S_i b_i,

with smoothing parameters \lambda_i, coefficient vectors b_i, and fixed penalty matrices S_i.

The qREML algorithm is typically much faster than REML or marginal ML using the full Laplace approximation method, but may be slightly less accurate regarding the estimation of the penalty strength parameters.

qreml_old is the previous implementation, kept for reference. Use qreml, which estimates the same quantities with a better behaved outer iteration.

Under the hood, qreml_old uses the R package RTMB for automatic differentiation in the inner optimisation. The user has to specify the penalised negative log-likelihood function pnll structured as dictated by RTMB and use the penalty function to compute the quadratic-form penalty inside the likelihood.

Usage

qreml_old(
  pnll,
  par,
  dat,
  random,
  map = NULL,
  silent = 1,
  psname = "lambda",
  alpha = 0.3,
  smoothing = 1,
  maxiter = 100,
  tol = 1e-04,
  max_halve = 5,
  method = "BFGS",
  control = list(),
  conv_crit = "relchange",
  spHess = FALSE,
  joint_unc = FALSE,
  saveall = FALSE
)

Arguments

Value

model object of class 'qremlModel'. This is a list containing:

References

Koslik, J. O. (2024). Efficient smoothness selection for nonparametric Markov-switching models via quasi restricted maximum likelihood. arXiv preprint arXiv:2411.11498.

See Also

penalty and penalty2 to compute the penalty inside the likelihood function

Examples

data = trex[1:1000,] # subset

# initial parameter list
par = list(logmu = log(c(0.3, 2.5)), # step mean
           logsigma = log(c(0.3, 1.5)), # step sd
           beta0 = c(-2,-2), # state process intercept
           betaspline = matrix(rep(0, 18), nrow = 2)) # state process spline coefs
          
# data object with initial penalty strength lambda
dat = list(step = data$step, # step length
           tod = data$tod, # time of day covariate
           N = 2, # number of states
           lambda = rep(10,2)) # initial penalty strength

# building model matrices
modmat = modmat(~ s(tod, bs = "cp"), 
                       data = data.frame(tod = 1:24), 
                       knots = list(tod = c(0,24))) # wrapping points
dat$Z = modmat$Z # spline design matrix
dat$S = modmat$S # penalty matrix

# penalised negative log-likelihood function
pnll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm_g(Z, cbind(beta0, betaspline), ad = TRUE) # transition probabilities
  delta = stationary_p(Gamma, t = tod[1] %% 24 + 1, ad = TRUE) # initial distribution
  mu = exp(logmu) # step mean
  sigma = exp(logsigma) # step sd
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step)) # only for non-NA obs.
  for(j in 1:N) allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])
  -forward_g(delta, Gamma[,,tod], allprobs) +
      penalty(betaspline, S, lambda) # this does all the penalization work
}

# model fitting
mod = qreml_old(pnll, par, dat, random = "betaspline", silent = 2)

Repeat initial penalty strengths for replicated smooths

Description

Expands the sp0 element of a modmat object into a sp vector for penalty or penalty2, when each smooth is used more than once, for example one spline per state or per off-diagonal element of the transition probability matrix.

Usage

rep_sp(model_matrices, N)

Arguments

Details

sp runs over replicates within each smooth, so for a tensor product the margins vary fastest. A single rep() over sp0 therefore only gives the right answer for some models, silently producing the right length but the wrong values for others. This function repeats each smooth's penalty strengths separately, which is always correct.

Value

named numeric vector of initial penalty strengths, of the length and in the order expected for sp

See Also

modmat, which computes sp0

Examples

modmat = modmat(~ s(x) + ti(x, y), data.frame(x = runif(100), y = runif(100)))
modmat$sp0         # one value per penalty
rep_sp(modmat, 2)  # ready to use as 'sp' for two replicates of each smooth

Per-observation score contributions

Description

Computes the score contribution of each observation,

u_t = \frac{\partial}{\partial \theta} \log p(x_t \mid x_1, \dots, x_{t-1}),

for a model whose negative log-likelihood function calls forward. Their column sums are the gradient of the log-likelihood. They are the building block for robust (sandwich) standard errors, e.g. cluster-robust ones via rowsum(U, ID).

Usage

scores(obj, par = NULL)

Arguments

Details

Each log-likelihood contribution in forward gets a weight w_t. The gradient of the weighted log-likelihood is linear in the weights, so its derivative w.r.t. w_t is u_t. This needs one tape of the gradient and one reverse sweep per parameter. The weights are handed to forward internally, so the likelihood function does not need to change.

The likelihood function is taken from obj and evaluated again, with the data as they are when scores() is called. Parameters fixed via map are respected, and the columns correspond to obj$par, like obj$gr() and obj$he().

If the likelihood function calls forward several times, e.g. once per track, the rows of the result follow the order of these calls. Only forward (and its deprecated _g and _p versions) is supported, and models with random effects are not.

Value

matrix of score contributions of dimension c(nObs, length(obj$par))

Examples

nll <- function(par) {
  getAll(par)
  Gamma <- tpm(eta)
  delta <- stationary(Gamma)
  allprobs <- cbind(dgamma2(step, mu[1], sigma[1]), 
                    dgamma2(step, mu[2], sigma[2]))
  -forward(delta, Gamma, allprobs)
}
step <- trex$step[1:200]
par <- list(eta = c(-2, -2), mu = c(0.3, 2.5), sigma = c(0.2, 1.5))
obj <- MakeADFun(nll, par, silent = TRUE)
opt <- nlminb(obj$par, obj$fn, obj$gr)
mod <- lama(obj)
U <- scores(mod)
colSums(U) # ~ 0 at the MLE

Monte Carlo version of sdreport

Description

After optimisation of an AD model, sdreportMC can be used to calculate samples of confidence intervals of all model parameters and REPORT()ed quantities including nonlinear functions of random effects and parameters.

Usage

sdreportMC(
  obj,
  what,
  nSamples = 1000,
  Hessian = NULL,
  CI = FALSE,
  probs = c(0.025, 0.975)
)

Arguments

Details

This function simply samples from the approximate multivariate normal distribution of the maximum likelihood estimate (MLE) of the parameters

\hat{\theta} \sim N(\theta, H^{-1}),

where H is the Hessian matrix of the negative log-likelihood function at the MLE. It then returns either the sampled parameters or REPORT()ed transformations of them. If CI = TRUE, it does not return the samples but directly computes confidence intervals.

If you are interested in several quantities, calling sdreportMC once with a vector what will generally be faster than calling it several times with single elements of what.

Value

named list corresponding to the elements of what. Each element has the structure of the corresponding quantity with an additional dimension added for the samples. For example, if a quantity is a vector, the list contains a matrix. If a quantity is a matrix, the list contains an array. If quantity is an array, the list contains an array with one extra dimension.

Examples

# fitting an HMM to the trex data and running sdreportMC
## negative log-likelihood function
nll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm(eta) # computes transition probability matrix from unconstrained eta
  delta = stationary(Gamma) # computes stationary distribution
  # exponentiating because all parameters strictly positive
  mu = exp(logmu)
  sigma = exp(logsigma)
  kappa = exp(logkappa)
  # reporting statements for sdreportMC
  REPORT(mu)
  REPORT(sigma)
  REPORT(kappa)
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step) & !is.na(angle)) # only for non-NA obs.
  for(j in 1:N){
    allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])*dvm(angle[ind],0,kappa[j])
  }
  -forward(delta, Gamma, allprobs) # simple forward algorithm
}

## initial parameter list
par = list(
 logmu = log(c(0.3, 1)),       # initial means for step length (log-transformed)
  logsigma = log(c(0.2, 0.7)), # initial sds for step length (log-transformed)
  logkappa = log(c(0.2, 0.7)), # initial concentration for turning angle (log-transformed)
  eta = rep(-2, 2)             # initial t.p.m. parameters (on logit scale)
)   
## data and hyperparameters
dat = list(
  step = trex$step[1:500],   # hourly step lengths
  angle = trex$angle[1:500], # hourly turning angles
  N = 2
)

## creating AD function
obj = MakeADFun(nll, par, silent = TRUE) # creating the objective function

## optimising
opt = nlminb(obj$par, obj$fn, obj$gr) # optimization

## running sdreportMC
# `mu` has report statement, `delta` is automatically reported by `forward()`
sdrMC = sdreportMC(obj, 
                   what = c("mu", "delta"), 
                   nSamples = 50)
dim(sdrMC$delta)
# now a matrix with 50 samples (rows)

Report uncertainty of the estimated smoothing parameters or variances

Description

Computes standard deviations for the smoothing parameters of a model object returned by qreml using the delta method.

Usage

sdreport_outer(mod, invert = FALSE)

Arguments

Details

The computations are based on the approximate gradient of the restricted log likelihood. The outer Hessian is computed by finite differencing of this gradient. If the inverse smoothing parameters are requested, the standard deviations are transformed to the variances using the delta method.

Value

list containing report matrix summarising parameters and standard deviations as well as the outer Hessian matrix.

Examples

## no examples

Simulate from a fitted model

Description

Simulates states and observations from a model fitted with forward and statedens. By default, this is done on the observed design: the tracks and their lengths, the initial distribution delta, the transition probabilities Gamma and the parameters of the state-dependent densities are those reported by the model, i.e. conditional on the observed covariates. Simulating for new covariates is not possible; instead, rebuild the likelihood function for the new data and call lama on the resulting object.

Usage

## S3 method for class 'LaMaModel'
simulate(
  object,
  nsim = 1,
  seed = NULL,
  rdens = NULL,
  random = FALSE,
  fill_na = FALSE,
  ...
)

Arguments

Details

The states are drawn as in the decoding functions, so that all types of Gamma that forward accepts (homogeneous, track-specific, time-varying, periodic) work, as does a continuous-time model with one transition matrix per time interval (the simulation is then conditional on the observed time points). Each data stream is drawn from the random number generator belonging to its density, e.g. rgamma2 for dgamma2.

Missing observations of the data stay missing (unless fill_na = TRUE), except for a stream that is missing completely: such a placeholder (e.g. step = rep(NA, n)) lets you simulate from a likelihood function without data, see the examples.

Models with random effects are by default simulated conditional on the estimated random effects. With random = TRUE, RTMB draws the random effects first, which requires that the density of each random effect is evaluated in the likelihood function before the random effect is used, either by u %~% dnorm(...) or by nll <- nll - sum(dnorm(u, ..., log = TRUE)) (or another density). The states and data are then simulated given these draws, so that every simulation has new random effects. This also requires that allprobs is the product of the data streams built by statedens, with each data stream passed once.

Value

data frame with the columns trackID (index of the track, if used), state and one column per data stream, for nsim = 1, else a list of nsim such data frames.

Examples

nll <- function(par) {
  getAll(par, dat)
  Gamma <- tpm(eta)
  allprobs <- statedens(step, dgamma2, exp(log_mu), exp(log_sigma)) *
    statedens(angle, dvm, 0, exp(log_kappa))
  -forward(stationary(Gamma), Gamma, allprobs)
}
dat <- list(step = trex$step[1:500], angle = trex$angle[1:500])
par <- list(eta = c(-2, -2), log_mu = log(c(0.3, 2.5)), log_sigma = log(c(0.2, 1.5)),
            log_kappa = log(c(0.2, 1.5)))
obj <- MakeADFun(nll, par, silent = TRUE)
opt <- nlminb(obj$par, obj$fn, obj$gr)
mod <- lama(obj)

sim <- simulate(mod, seed = 1)
head(sim)

# simulate without data: all-NA streams as placeholders, parameters of your choice, then refit
dat$step <- dat$angle <- rep(NA, 500)
obj0 <- MakeADFun(nll, par, silent = TRUE)
sim <- simulate(lama(obj0), seed = 1)
dat$step <- sim$step; dat$angle <- sim$angle
obj1 <- MakeADFun(nll, par, silent = TRUE)
opt1 <- nlminb(obj1$par, obj1$fn, obj1$gr)

# random effects: RTMB draws them first, then the states and data are simulated given them
nll_re <- function(par) {
  getAll(par, dat)
  u %~% dnorm(0, exp(log_su)) # track-level random effect on the mean step length in state 1
  mu <- exp(cbind(log_mu[1] + u[id], log_mu[2]))
  Gamma <- tpm(eta)
  .nll <- .nll - forward(stationary(Gamma), Gamma, statedens(step, dgamma2, mu, exp(log_sigma)), id)
  .nll
}
dat <- list(step = trex$step[1:500], id = rep(1:5, each = 100))
par_re <- list(eta = c(-2, -2), log_mu = log(c(0.3, 2.5)), log_sigma = log(c(0.2, 1.5)),
               log_su = log(0.3), u = rep(0, 5))
obj_re <- MakeADFun(nll_re, par_re, random = "u", silent = TRUE)
sim_re <- simulate(lama(obj_re), random = TRUE, seed = 1)

Skew normal distribution

Description

Density, distribution function, quantile function and random generation for the skew normal distribution.

Usage

dskewnorm(x, xi = 0, omega = 1, alpha = 0, log = FALSE)

pskewnorm(q, xi = 0, omega = 1, alpha = 0, lower.tail = TRUE, log.p = FALSE)

qskewnorm(p, xi = 0, omega = 1, alpha = 0, lower.tail = TRUE, log.p = FALSE)

rskewnorm(n, xi = 0, omega = 1, alpha = 0)

Arguments

Details

This implementation of dskewnorm allows for automatic differentiation with RTMB while the other functions are imported from the sn package.

Value

dskewnorm gives the density, pskewnorm gives the distribution function, qskewnorm gives the quantile function, and rskewnorm generates random deviates.

Examples

x = rskewnorm(1)
d = dskewnorm(x)
p = pskewnorm(x)
q = qskewnorm(p)

Build the design and penalty matrices for smooth density estimation

Description

smooth_dens_construct() is deprecated, use smoothdens instead. It returns the same objects, but ordered by component first and data stream second.

Usage

smooth_dens_construct(
  data,
  par,
  type = "real",
  k = 25,
  knots = NULL,
  degree = 3,
  diff_order = 2
)

Arguments

Value

a nested list containing the design matrices Z, the penalty matrices S, the initial coefficients coef the prediction design matrices Z_predict, the prediction grids xseq, and details for the basis expansion for each data stream.

Examples

data = data.frame(x = rnorm(100))
# deprecated:
# SmoothDens = smooth_dens_construct(data, par = list(x = list(mean = 0, sd = 1)))
# use instead:
dens = smoothdens(data$x, par = list(mean = 0, sd = 1))

Build the design and penalty matrices for smooth density estimation

Description

This high-level function can be used to prepare objects needed to estimate mixture models of smooth densities using P-Splines.

Usage

smoothdens(
  x,
  par,
  type = "real",
  k = 25,
  knots = NULL,
  degree = 3,
  diff_order = 2,
  pow = 0.5,
  npoints = 10000
)

Arguments

Details

Under the hood, make_matrices_dens is used for the actual construction of the design and penalty matrices.

You can provide one or multiple data streams of different types (real, positive, circular) and specify initial means and standard deviations/ concentrations for each data stream. This information is then converted into suitable spline coefficients. smoothdens then constructs the design and penalty matrices for standardised B-splines basis functions (integrating to one) for each data stream. For types "real" and "circular" the knots are placed equidistant in the range of the data, for type "positive" the knots are placed using polynomial spacing.

The last spline coefficient of each state is fixed at zero for identifiability, hence coef has k - 1 columns. density(beta) takes care of this: it maps the coefficient matrix to the allprobs matrix, with rows of missing observations equal to one.

Value

For one data stream, an object of class "smoothdens": a list containing the design matrix Z, the penalty matrix S, the initial coefficients coef, the prediction design matrix Z_predict, the prediction grid xseq, details for the basis expansion basis, and the function density(beta). For several data streams, a named list of such objects.

See Also

predict for the estimated densities

Examples

## one data stream, mixture of two normals
x = c(rnorm(100, 0, 1), rnorm(100, 5, 1))
dens = smoothdens(x, par = list(mean = c(0, 5), sd = c(1, 1)))

Z = dens$Z # design matrix
S = dens$S # penalty matrix
beta = dens$coef # initial coefficients
allprobs = dens$density(beta) # use this in the likelihood function

## 3 data streams
data = data.frame(x1 = rnorm(100, mean = 0, sd = 1),
                  x2 = rgamma2(100, mean = 5, sd = 3),
                  x3 = rvm(100, mu = 0, kappa = 2))

par = list(x1 = list(mean = 0, sd = 1),
           x2 = list(mean = 5, sd = 3),
           x3 = list(mean = 0, kappa = 2))

dens = smoothdens(data, par, type = c("real", "positive", "circular"))

# objects for x1
Z1 = dens$x1$Z
allprobs = dens$x1$density(dens$x1$coef) *
  dens$x2$density(dens$x2$coef) *
  dens$x3$density(dens$x3$coef)

State-dependent densities of one data stream

Description

Builds the matrix of state-dependent densities (or probabilities) of one data stream, as used in the allprobs matrix of forward. Several data streams that are conditionally independent given the state are combined by multiplying their matrices.

Usage

statedens(x, dens, ..., nStates = NULL, log = FALSE)

Arguments

Details

Each parameter in ... can be

• a single value, used for all states,

• a vector of length nStates, one value per state, or

• a matrix of dimension c(nObs, nStates), one value per observation and state, e.g. from a regression on covariates.

Missing observations contribute nothing to the likelihood, so their rows are 1 (0 with log = TRUE). As each data stream gets its own matrix, a missing value in one stream does not remove the observations of the others.

With log = TRUE, the log-densities are returned, to be added over data streams and used with forward(..., logspace = TRUE). They are computed by dens(..., log = TRUE) if dens has an argument log, otherwise as log(dens(...)) (e.g. for densities from zero_inflate).

In obj$simulate() of RTMB, x is registered as data and simulated by forward, see simulate.LaMaModel.

Writing the loop over states by hand is equally valid and gives the same result. statedens() is just a shortcut.

Value

matrix of dimension c(nObs, nStates)

Examples

step <- trex$step[1:200]
angle <- trex$angle[1:200]
mu <- c(0.3, 2.5); sigma <- c(0.2, 1.5); kappa <- c(0.2, 1.5)

# one data stream
allprobs <- statedens(step, dgamma2, mu, sigma)

# the same by hand
allprobs2 <- matrix(1, length(step), 2)
ind <- which(!is.na(step))
for(j in 1:2) allprobs2[ind, j] <- dgamma2(step[ind], mu[j], sigma[j])
all.equal(allprobs, allprobs2)

# two data streams, mean of the turning angle the same in both states
allprobs <- statedens(step, dgamma2, mu, sigma) * statedens(angle, dvm, 0, kappa)

# on the log scale, for forward(..., logspace = TRUE)
logprobs <- statedens(step, dgamma2, mu, sigma, log = TRUE) +
  statedens(angle, dvm, 0, kappa, log = TRUE)

Calculate conditional local state probabilities in HMMs

Description

Computes

\Pr(\text{State}_t = j \mid X_1, ..., X_T)

for a given HMM.

Usage

stateprobs(
  delta,
  Gamma,
  allprobs,
  trackID = NULL,
  mod = NULL,
  forecast = FALSE,
  tod = NULL
)

Arguments

Value

matrix of conditional state probabilities of dimension c(nObs, nStates)

See Also

Other decoding functions: stateprobs_g(), stateprobs_p(), viterbi(), viterbi_g(), viterbi_p()

Examples

Gamma = tpm(c(-1,-2))
delta = stationary(Gamma)
allprobs = matrix(runif(10), nrow = 10, ncol = 2)

probs = stateprobs(delta, Gamma, allprobs)

Calculate conditional local state probabilities for inhomogeneous HMMs

Description

Computes

\Pr(\text{State}_t = j \mid X_1, ..., X_T)

for inhomogeneous HMMs

Usage

stateprobs_g(
  delta,
  Gamma,
  allprobs,
  trackID = NULL,
  mod = NULL,
  forecast = FALSE
)

Arguments

Details

Deprecated: please use stateprobs for all model types instead.

Value

matrix of conditional state probabilities of dimension c(nObs, nStates)

See Also

Other decoding functions: stateprobs(), stateprobs_p(), viterbi(), viterbi_g(), viterbi_p()

Examples

Gamma = tpm_g(runif(10), matrix(c(-1,-1,1,-2), nrow = 2, byrow = TRUE))
delta = c(0.5, 0.5)
allprobs = matrix(runif(20), nrow = 10, ncol = 2)

probs = stateprobs_g(delta, Gamma[,,-1], allprobs)

Calculate conditional local state probabilities for periodically inhomogeneous HMMs

Description

Computes

\Pr(S_t = j \mid X_1, ..., X_T)

for periodically inhomogeneous HMMs

Usage

stateprobs_p(
  delta,
  Gamma,
  allprobs,
  tod,
  trackID = NULL,
  mod = NULL,
  forecast = FALSE
)

Arguments

Details

Deprecated: please use stateprobs with argument tod instead.

Value

matrix of conditional state probabilities of dimension c(n,N)

See Also

Other decoding functions: stateprobs(), stateprobs_g(), viterbi(), viterbi_g(), viterbi_p()

Examples

L = 24
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma = tpm_p(1:L, L, beta, degree = 1)
delta = stationary_p(Gamma, 2) # for tod[1] = 1
allprobs = matrix(runif(200), nrow = 100, ncol = 2)
tod = rep(1:24, 5)[1:100]

probs = stateprobs_p(delta, Gamma, allprobs, tod)

Compute the stationary distribution of a homogeneous Markov chain

Description

A homogeneous, finite state Markov chain that is irreducible and aperiodic converges to a unique stationary distribution, here called \delta. As it is stationary, this distribution satisfies

\delta \Gamma = \delta,

subject to \sum_{j=1}^N \delta_j = 1, where \Gamma is the transition probability matrix. This function solves the linear system of equations above.

Usage

stationary(Gamma)

Arguments

Value

either a single stationary distribution of the Markov chain (vector of length N) or a matrix of stationary distributions of dimension c(nTracks,N) with one stationary distribution in each row

See Also

tpm to create a transition probabilty matrix using the multinomial logistic link (softmax)

Other stationary distribution functions: stationary_ct(), stationary_p()

Examples

# single matrix
Gamma = tpm(c(rep(-2,3), rep(-3,3)))
delta = stationary(Gamma)
# multiple matrices
Gamma = array(Gamma, dim = c(3,3,10))
Delta = stationary(Gamma)

Compute the stationary distribution of a continuous-time Markov chain

Description

A well-behaved continuous-time Markov chain converges to a unique stationary distribution, here called \pi. This distribution satisfies

\pi Q = 0,

subject to \sum_{j=1}^N \pi_j = 1, where Q is the infinitesimal generator of the Markov chain. This function solves the linear system of equations above for a given generator matrix.

Usage

stationary_ct(Q)

stationary_cont(Q)

Arguments

Value

either a single stationary distribution of the continuous-time Markov chain (vector of length N) or a matrix of stationary distributions of dimension c(nTracks,N) with one stationary distribution in each row

See Also

generator to create a generator matrix

Other stationary distribution functions: stationary(), stationary_p()

Examples

# single matrix
Q = generator(c(-2,-2))
Pi = stationary_ct(Q)
# multiple matrices
Q = array(Q, dim = c(2,2,10))
Pi = stationary_ct(Q)

Periodically stationary distribution of a periodically inhomogeneous Markov chain

Description

Computes the periodically stationary distribution of a periodically inhomogeneous Markov chain.

Usage

stationary_p(Gamma, t = NULL, ad = NULL)

Arguments

Details

If the transition probability matrix of an inhomogeneous Markov chain varies only periodically (with period length L), it converges to a so-called periodically stationary distribution. This happens, because the thinned Markov chain, which has a full cycle as each time step, has homogeneous transition probability matrix

\Gamma_t = \Gamma^{(t)} \Gamma^{(t+1)} \dots \Gamma^{(t+L-1)}

for all t = 1, \dots, L. The stationary distribution for time t satifies \delta^{(t)} \Gamma_t = \delta^{(t)}. Here, Gamma[,,t] is the transition from t to t+1. forward with tod uses Gamma[,,tod[t]] for the transition into t, so there, the initial distribution is stationary_p(Gamma, tod[1] %% L + 1), and stationary_p(Gamma)[c(2:L, 1), ] has the distribution at tod = k in row k.

This function calculates said periodically stationary distribution.

Value

either the periodically stationary distribution at time t or all periodically stationary distributions.

References

When using this function, please cite:

Koslik, J. O., Feldmann, C. C., Mews, S., Michels, R., & Langrock, R. (2025). Inference on the state process of periodically inhomogeneous hidden Markov models for animal behavior. The Annals of Applied Statistics, 19(4), 2724–2737.

See Also

tpm_p and tpm_g to create multiple transition matrices based on a cyclic variable or design matrix

Other stationary distribution functions: stationary(), stationary_ct()

Examples

# setting parameters for trigonometric link
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma = tpm_p(beta = beta, degree = 1)
# periodically stationary distribution for specific time point
delta = stationary_p(Gamma, 4)

# all periodically stationary distributions
Delta = stationary_p(Gamma)

Sparse version of stationary_p

Description

This is function computes the periodically stationary distribution of a Markov chain given a list of L sparse transition probability matrices. Compatible with automatic differentiation by RTMB

Usage

stationary_p_sparse(Gamma, t = NULL)

Arguments

Value

either the periodically stationary distribution at time t or all periodically stationary distributions.

Examples

## periodic HSMM example (here the approximating tpm is sparse)
N = 2 # number of states
L = 24 # cycle length
# time-varying mean dwell times
Z = cosinor(1:L, period = L) # trigonometric basis functions design matrix
beta = matrix(c(2, 2, 0.1, -0.1, -0.2, 0.2), nrow = 2)
Lambda = exp(cbind(1, Z) %*% t(beta))
sizes = c(20, 20) # approximating chain with 40 states
# state dwell-time distributions
dm = lapply(1:N, function(i) sapply(1:sizes[i]-1, dpois, lambda = Lambda[,i]))
omega = matrix(c(0,1,1,0), nrow = N, byrow = TRUE) # embedded t.p.m.

# calculating extended-state-space t.p.m.s
Gamma = tpm_phsmm(omega, dm)
# Periodically stationary distribution for specific time point
delta = stationary_p_sparse(Gamma, 4)

# All periodically stationary distributions
Delta = stationary_p_sparse(Gamma)

Sparse version of stationary

Description

This is function computes the stationary distribution of a Markov chain with a given sparse transition probability matrix. Compatible with automatic differentiation by RTMB

Usage

stationary_sparse(Gamma)

Arguments

Value

stationary distribution of the Markov chain with the given transition probability matrix

Examples

## HSMM example (here the approximating tpm is sparse)
# building the t.p.m. of the embedded Markov chain
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# defining state aggregate sizes
sizes = c(20, 30)
# defining state dwell-time distributions
lambda = c(5, 11)
dm = list(dpois(1:sizes[1]-1, lambda[1]), dpois(1:sizes[2]-1, lambda[2]))
# calculating extended-state-space t.p.m.
Gamma = tpm_hsmm(omega, dm)
delta = stationary_sparse(Gamma)

Summary of a LaMaModel object

Description

Summarises a model from lama or qreml. Besides the overview of print(), this contains

• for a homogeneous state process (or one per track) with at most 10 states: the transition probability matrix (or generator matrix Q in continuous time, if reported, e.g. by tpm_ct), its stationary distribution and mean dwell times, averaged over tracks if there is one per track,

• for a periodic state process: the stationary distribution averaged over the cycle, i.e. the long-run proportion of time spent in each state,

• the quantities reported by the likelihood function via REPORT(): numeric vectors of length equal to the number of states are shown as one table with a column per state, other small objects separately, and large ones only by name,

• log-likelihood, AIC, BIC and the maximum gradient component at the estimate.

Usage

## S3 method for class 'LaMaModel'
summary(object, ...)

## S3 method for class 'summary.LaMaModel'
print(x, digits = max(3, getOption("digits") - 3), ...)

Arguments

Value

summary() returns an object of class "summary.LaMaModel", a list containing the quantities shown, e.g. Gamma, stationary, dwell and bystate. print() returns x invisibly.

Summary of a qremlModel object

Description

Summarises a model from qreml: the summary of summary.LaMaModel, plus the effective degrees of freedom and smoothing parameters of the smooth terms and diagnostics of the smoothness selection.

Usage

## S3 method for class 'qremlModel'
summary(object, ...)

## S3 method for class 'summary.qremlModel'
print(x, digits = max(3, getOption("digits") - 3), ...)

Arguments

Value

summary() returns an object of class "summary.qremlModel", extending "summary.LaMaModel" by edf, lambda and selection. print() returns x invisibly.

Build the transition probability matrix from unconstrained parameter vector

Description

Markov chains are parametrised in terms of a transition probability matrix \Gamma, for which each row contains a conditional probability distribution of the next state given the current state. Hence, each row has entries between 0 and 1 that need to sum to one.

For numerical optimisation, we parameterise in terms of unconstrained parameters, thus this function computes said matrix from an unconstrained parameter vector via the inverse multinomial logistic link (also known as softmax) applied to each row.

Usage

tpm(
  beta,
  Z = NULL,
  Eta = NULL,
  byrow = FALSE,
  ref = NULL,
  ad = NULL,
  report = TRUE,
  param = NULL
)

Arguments

Value

Transition probability matrix of dimension c(nStates, nStates) or array of such matrices of dimension c(nStates, nStates, nObs) if Z or Eta is provided.

See Also

Other transition probability matrix functions: generator(), generator_g(), tpm_ct(), tpm_emb(), tpm_emb_g(), tpm_g(), tpm_g2(), tpm_p()

Examples

## homogeneous Markov chain
# 2 states: 2 = 2*(2-1) free off-diagonal elements
par <- rep(-2, 2)
Gamma <- tpm(par)
# 3 states: 6 = 3*(3-1) free off-diagonal elements
par <- rep(-3, 6)
Gamma <- tpm(par)
# 4 states: 12 = 4*(4-1) free off-diagonal elements
par <- rep(-4, 12)
Gamma <- tpm(par)

## inhomogeneous Markov chain
# t.p.m. depends on covariates
z1 <- runif(100); z2 <- runif(100) # 2 covariates
Z <- cbind(1, z1, z2) # design matrix
beta0 <- c(-2, -2); beta1 = c(1, -2); beta2 = c(2, -1) # coefficients for intercept and covariates
beta <- cbind(beta0, beta1, beta2) # coefficient matrix; with intercepts!
Gamma <- tpm(beta, Z) # array with 100 slices

Calculate continuous time transition probabilities

Description

A continuous-time Markov chain is described by an infinitesimal generator matrix Q. When observing data at time points t_1, \dots, t_n the transition probabilites between t_i and t_{i+1} are caluclated as

\Gamma(\Delta t_i) = \exp(Q \Delta t_i),

where \exp() is the matrix exponential. The mapping \Gamma(\Delta t) is also called the Markov semigroup. This function calculates all transition matrices based on a given generator and time differences.

Usage

tpm_ct(Q, timediff, rates = NULL, ad = NULL, report = TRUE)

tpm_cont(Q, timediff, rates = NULL, ad = NULL, report = TRUE)

Arguments

Value

array of continuous-time transition matrices of dimension c(nStates, nStates, nObs-1)

See Also

Other transition probability matrix functions: generator(), generator_g(), tpm(), tpm_emb(), tpm_emb_g(), tpm_g(), tpm_g2(), tpm_p()

Examples

# building a Q matrix for a 3-state cont.-time Markov chain
Q = generator(rep(-2, 6))

# draw random time differences
timediff = rexp(100, 10)

# compute all transition matrices
Gamma = tpm_cont(Q, timediff)

Build the embedded transition probability matrix of an HSMM from unconstrained parameter vector

Description

Hidden semi-Markov models are defined in terms of state durations and an embedded transition probability matrix that contains the conditional transition probabilities given that the current state is left. This matrix necessarily has diagonal entries all equal to zero as self-transitions are impossible.

This function builds such an embedded/ conditional transition probability matrix from an unconstrained parameter vector. For each row of the matrix, the inverse multinomial logistic link is applied.

For a matrix of dimension c(nStates, nStates), the number of free off-diagonal elements is nStates * (nStates-2), hence also the length of param. This means, for 2 states, the function needs to be called without any arguments, for 3-states with a vector of length 3, for 4 states with a vector of length 8, etc.

Compatible with automatic differentiation by RTMB

Usage

tpm_emb(param = NULL)

Arguments

Value

embedded/ conditional transition probability matrix of dimension c(nStates, nStates)

See Also

Other transition probability matrix functions: generator(), generator_g(), tpm(), tpm_ct(), tpm_emb_g(), tpm_g(), tpm_g2(), tpm_p()

Examples

# 2 states: no free off-diagonal elements
omega = tpm_emb()

# 3 states: 3 free off-diagonal elements
param = rep(0, 3)
omega = tpm_emb(param)

# 4 states: 8 free off-diagonal elements
param = rep(0, 8)
omega = tpm_emb(param)

Build all embedded transition probability matrices of an inhomogeneous HSMM

Description

Hidden semi-Markov models are defined in terms of state durations and an embedded transition probability matrix that contains the conditional transition probabilities given that the current state is left. This matrix necessarily has diagonal entries all equal to zero as self-transitions are impossible. We can allow this matrix to vary with covariates, which is the purpose of this function.

It builds all embedded/ conditional transition probability matrices based on a design and parameter matrix. For each row of the matrix, the inverse multinomial logistic link is applied.

For a matrix of dimension c(nStates, nStates), the number of free off-diagonal elements is nStates * (nStates-2) which determines the number of rows of the parameter matrix.

Compatible with automatic differentiation by RTMB

Usage

tpm_emb_g(Z, beta, report = TRUE)

Arguments

Value

array of embedded/ conditional transition probability matrices of dimension c(nStates, nStates, nObs)

See Also

Other transition probability matrix functions: generator(), generator_g(), tpm(), tpm_ct(), tpm_emb(), tpm_g(), tpm_g2(), tpm_p()

Examples

## parameter matrix for 3-state HSMM
beta = matrix(c(rep(0, 3), -0.2, 0.2, 0.1), nrow = 3)
# no intercept
Z = rnorm(100)
omega = tpm_emb_g(Z, beta)
# intercept
Z = cbind(1, Z)
omega = tpm_emb_g(Z, beta)

Build all transition probability matrices of an inhomogeneous HMM

Description

In an HMM, we often model the influence of covariates on the state process by linking them to the transition probabiltiy matrix. Most commonly, this is done by specifying a linear predictor

\eta_{ij}^{(t)} = \beta^{(ij)}_0 + \beta^{(ij)}_1 z_{t1} + \dots + \beta^{(ij)}_p z_{tp}

for each off-diagonal element (i \neq j) of the transition probability matrix and then applying the inverse multinomial logistic link (also known as softmax) to each row. This function efficiently calculates all transition probabilty matrices for a given design matrix Z and parameter matrix beta.

Usage

tpm_g(Z, beta, Eta = NULL, byrow = FALSE, ref = NULL, ad = NULL, report = TRUE)

Arguments

Details

Deprecated: please use tpm(beta, Z) for all model types instead.

Value

array of transition probability matrices of dimension c(nStates, nStates, nObs)

See Also

Other transition probability matrix functions: generator(), generator_g(), tpm(), tpm_ct(), tpm_emb(), tpm_emb_g(), tpm_g2(), tpm_p()

Examples

## inhomogeneous Markov chain
# t.p.m. depends on covariates
z1 <- runif(100); z2 <- runif(100) # 2 covariates
Z <- cbind(1, z1, z2) # design matrix
beta0 <- c(-2, -2); beta1 = c(1, -2); beta2 = c(2, -1) # coefficients for intercept and covariates
beta <- cbind(beta0, beta1, beta2) # coefficient matrix; with intercepts!
#' Gamma <- tpm(beta, Z) # array with 100 slices

Build all transition probability matrices of an inhomogeneous HMM

Description

In an HMM, we often model the influence of covariates on the state process by linking them to the transition probabiltiy matrix. Most commonly, this is done by specifying linear predictors

\eta_{ij}^{(t)} = \beta^{(ij)}_0 + \beta^{(ij)}_1 z_{t1} + \dots + \beta^{(ij)}_p z_{tp}

for each off-diagonal element (i \neq j) of the transition probability matrix and then applying the inverse multinomial logistic link (also known as softmax) to each row. This function efficiently calculates all transition probabilty matrices for a given design matrix Z and parameter matrix beta.

Usage

tpm_g2(Z, beta, byrow = FALSE, ad = NULL, report = TRUE, ref = NULL)

Arguments

Value

array of transition probability matrices of dimension c(N,N,n)

See Also

Other transition probability matrix functions: generator(), generator_g(), tpm(), tpm_ct(), tpm_emb(), tpm_emb_g(), tpm_g(), tpm_p()

Examples

Z = matrix(runif(200), ncol = 2)
beta = matrix(c(-1, 1, 2, -2, 1, -2), nrow = 2, byrow = TRUE)
Gamma = tpm_g(Z, beta)

Builds the transition probability matrix of an HSMM-approximating HMM

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs, where the state duration distribution is explicitly modelled. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function computes the transition matrix to approximate a given HSMM by an HMM with a larger state space.

Usage

tpm_hsmm(omega, dm, Fm = NULL, sparse = TRUE, eps = 1e-10)

Arguments

Value

extended-state-space transition probability matrix of the approximating HMM

Examples

# building the t.p.m. of the embedded Markov chain
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# defining state aggregate sizes
sizes = c(20, 30)
# defining state dwell-time distributions
lambda = c(5, 11)
dm = list(dpois(1:sizes[1]-1, lambda[1]), dpois(1:sizes[2]-1, lambda[2]))
# calculating extended-state-space t.p.m.
Gamma = tpm_hsmm(omega, dm)

Build the transition probability matrix of an HSMM-approximating HMM

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities. This function computes the transition matrix of an HSMM.

Usage

tpm_hsmm2(omega, dm, eps = 1e-10)

Arguments

Value

extended-state-space transition probability matrix of the approximating HMM

Examples

# building the t.p.m. of the embedded Markov chain
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# defining state aggregate sizes
sizes = c(20, 30)
# defining state dwell-time distributions
lambda = c(5, 11)
dm = list(dpois(1:sizes[1]-1, lambda[1]), dpois(1:sizes[2]-1, lambda[2]))
# calculating extended-state-space t.p.m.
Gamma = tpm_hsmm(omega, dm)

Builds all transition probability matrices of an inhomogeneous-HSMM-approximating HMM

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function computes the transition matrices of a periodically inhomogeneos HSMMs.

Usage

tpm_ihsmm(omega, dm, eps = 1e-10)

Arguments

Value

list of dimension length n - max(sapply(dm, ncol)), containing sparse extended-state-space transition probability matrices for each time point (except the first max(sapply(dm, ncol)) - 1).

Examples

N = 2
# time-varying mean dwell times
n = 100
z = runif(n)
beta = matrix(c(2, 2, 0.1, -0.1), nrow = 2)
Lambda = exp(cbind(1, z) %*% t(beta))
sizes = c(15, 15) # approximating chain with 30 states
# state dwell-time distributions
dm = lapply(1:N, function(i) sapply(1:sizes[i]-1, dpois, lambda = Lambda[,i]))

## homogeneous conditional transition probabilites
# diagonal elements are zero, rowsums are one
omega = matrix(c(0,1,1,0), nrow = N, byrow = TRUE)

# calculating extended-state-space t.p.m.s
Gamma = tpm_ihsmm(omega, dm)

## inhomogeneous conditional transition probabilites
# omega can be an array
omega = array(omega, dim = c(N,N,n))

# calculating extended-state-space t.p.m.s
Gamma = tpm_ihsmm(omega, dm)

Build all transition probability matrices of a periodically inhomogeneous HMM

Description

Given a periodically varying variable such as time of day or day of year and the associated cycle length, this function calculates the transition probability matrices by applying the inverse multinomial logistic link (also known as softmax) to linear predictors of the form

\eta^{(t)}_{ij} = \beta_0^{(ij)} + \sum_{k=1}^K \bigl( \beta_{1k}^{(ij)} \sin(\frac{2 \pi k t}{L}) + \beta_{2k}^{(ij)} \cos(\frac{2 \pi k t}{L}) \bigr)

for the off-diagonal elements (i \neq j) of the transition probability matrix. This is relevant for modeling e.g. diurnal variation and the flexibility can be increased by adding smaller frequencies (i.e. increasing K).

Usage

tpm_p(
  tod = 1:24,
  L = 24,
  beta,
  degree = 1,
  Z = NULL,
  byrow = FALSE,
  ad = NULL,
  report = TRUE
)

Arguments

Details

Note that using this function inside the negative log-likelihood function is convenient, but it performs the basis expansion into sine and cosine terms each time it is called. As these do not change during the optimisation, using tpm_g with a pre-calculated (by trigBasisExp) design matrix would be more efficient.

Value

array of transition probability matrices of dimension c(nStates, nStates, length(tod))

See Also

Other transition probability matrix functions: generator(), generator_g(), tpm(), tpm_ct(), tpm_emb(), tpm_emb_g(), tpm_g(), tpm_g2()

Examples

# hourly data 
tod = seq(1, 24, by = 1)
L = 24
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma = tpm_p(tod, L, beta, degree = 1)

# half-hourly data
## integer tod sequence
tod = seq(1, 48, by = 1)
L = 48
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma1 = tpm_p(tod, L, beta, degree = 1)

## equivalent specification
tod = seq(0.5, 24, by = 0.5)
L = 24
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma2 = tpm_p(tod, L, beta, degree = 1)

all(Gamma1 == Gamma2) # same result

Builds all transition probability matrices of an periodic-HSMM-approximating HMM

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function computes the transition matrices of a periodically inhomogeneos HSMMs.

Usage

tpm_phsmm(omega, dm, eps = 1e-10)

Arguments

Value

list of dimension length L, containing sparse extended-state-space transition probability matrices of the approximating HMM for each time point of the cycle.

Examples

N = 2 # number of states
L = 24 # cycle length
# time-varying mean dwell times
Z = cosinor(1:L, period = L) # trigonometric basis functions design matrix
beta = matrix(c(2, 2, 0.1, -0.1, -0.2, 0.2), nrow = 2)
Lambda = exp(cbind(1, Z) %*% t(beta))
sizes = c(20, 20) # approximating chain with 40 states
# state dwell-time distributions
dm = lapply(1:N, function(i) sapply(1:sizes[i]-1, dpois, lambda = Lambda[,i]))

## homogeneous conditional transition probabilites
# diagonal elements are zero, rowsums are one
omega = matrix(c(0,1,1,0), nrow = N, byrow = TRUE)

# calculating extended-state-space t.p.m.s
Gamma = tpm_phsmm(omega, dm)

## inhomogeneous conditional transition probabilites
# omega can be an array
omega = array(omega, dim = c(N,N,L))

# calculating extended-state-space t.p.m.s
Gamma = tpm_phsmm(omega, dm)

Build all transition probability matrices of an periodic-HSMM-approximating HMM

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities. This function computes the transition matrices of a periodically inhomogeneos HSMMs.

Usage

tpm_phsmm2(omega, dm, eps = 1e-10)

Arguments

Value

array of dimension c(N,N,L), containing the extended-state-space transition probability matrices of the approximating HMM for each time point of the cycle.

Examples

N = 3
L = 24
# time-varying mean dwell times
Lambda = exp(matrix(rnorm(L*N, 2, 0.5), nrow = L))
sizes = c(25, 25, 25) # approximating chain with 75 states
# state dwell-time distributions
dm = list()
for(i in 1:3){
  dmi = matrix(nrow = L, ncol = sizes[i])
  for(t in 1:L){
    dmi[t,] = dpois(1:sizes[i]-1, Lambda[t,i])
  }
  dm[[i]] = dmi
}

## homogeneous conditional transition probabilites
# diagonal elements are zero, rowsums are one
omega = matrix(c(0,0.5,0.5,0.2,0,0.8,0.7,0.3,0), nrow = N, byrow = TRUE)

# calculating extended-state-space t.p.m.s
Gamma = tpm_phsmm(omega, dm)

## inhomogeneous conditional transition probabilites
# omega can be an array
omega = array(rep(omega,L), dim = c(N,N,L))
omega[1,,4] = c(0, 0.2, 0.8) # small change for inhomogeneity

# calculating extended-state-space t.p.m.s
Gamma = tpm_phsmm(omega, dm)

Compute the transition probability matrix of a thinned periodically inhomogeneous Markov chain.

Description

If the transition probability matrix of an inhomogeneous Markov chain varies only periodically (with period length L), it converges to a so-called periodically stationary distribution. This happens, because the thinned Markov chain, which has a full cycle as each time step, has homogeneous transition probability matrix

\Gamma_t = \Gamma^{(t)} \Gamma^{(t+1)} \dots \Gamma^{(t+L-1)}

for all t = 1, \dots, L. This function calculates the matrix above efficiently as a preliminery step to calculating the periodically stationary distribution.

Usage

tpm_thinned(Gamma, t)

Arguments

Value

thinned transition probabilty matrix of dimension c(nStates, nStates)

Examples

# setting parameters for trigonometric link
beta = matrix(c(-1, -2, 2, -1, 2, -4), nrow = 2, byrow = TRUE)
# calculating periodically varying t.p.m. array (of length 24 here)
Gamma = tpm_p(beta = beta)
# calculating t.p.m. of thinned Markov chain
tpm_thinned(Gamma, 4)

T-Rex Movement Data

Description

Hourly step lengths and turning angles of a Tyrannosaurus rex, living 66 million years ago.

Usage

trex

Format

A data frame with 10.000 rows and 4 variables:

tod

time of day variable ranging from 1 to 24

step

hourly step lengths in kilometres

angle

hourly turning angles in radians

state

hidden state variable

Source

Generated for example purposes.

Compute the design matrix for a trigonometric basis expansion

Description

Given a periodically varying variable such as time of day or day of year and the associated cycle length, this function performs a basis expansion to efficiently calculate a linear predictor of the form

\eta^{(t)} = \beta_0 + \sum_{k=1}^K \bigl( \beta_{1k} \sin(\frac{2 \pi k t}{L}) + \beta_{2k} \cos(\frac{2 \pi k t}{L}) \bigr).

This is relevant for modeling e.g. diurnal variation and the flexibility can be increased by adding smaller frequencies (i.e. increasing K).

Usage

trigBasisExp(tod, L = 24, degree = 1)

Arguments

Value

design matrix (without intercept column), ordered as sin1, cos1, sin2, cos2, ...

Examples

# no examples

Viterbi algorithm for state decoding in HMMs

Description

The Viterbi algorithm decodes the most probable state sequence of an HMM.

Usage

viterbi(delta, Gamma, allprobs, trackID = NULL, mod = NULL, tod = NULL)

Arguments

Value

vector of decoded states of length nObs

See Also

Other decoding functions: stateprobs(), stateprobs_g(), stateprobs_p(), viterbi_g(), viterbi_p()

Examples

delta = c(0.5, 0.5)
Gamma = matrix(c(0.9, 0.1, 0.2, 0.8), nrow = 2, byrow = TRUE)
allprobs = matrix(runif(200), nrow = 100, ncol = 2)
states = viterbi(delta, Gamma, allprobs)

Viterbi algorithm for state decoding in inhomogeneous HMMs

Description

The Viterbi algorithm decodes the most probable state sequence of an HMM.

Usage

viterbi_g(delta, Gamma, allprobs, trackID = NULL, mod = NULL)

Arguments

Details

Deprecated: please use viterbi for all model types instead.

Value

vector of decoded states of length nObs

See Also

Other decoding functions: stateprobs(), stateprobs_g(), stateprobs_p(), viterbi(), viterbi_p()

Examples

delta = c(0.5, 0.5)
Gamma = tpm_g(runif(10), matrix(c(-2,-2,1,-1), nrow = 2))
allprobs = matrix(runif(20), nrow = 10, ncol = 2)
states = viterbi_g(delta, Gamma[,,-1], allprobs)

Viterbi algorithm for state decoding in periodically inhomogeneous HMMs

Description

The Viterbi algorithm allows one to decode the most probable state sequence of an HMM.

Usage

viterbi_p(delta, Gamma, allprobs, tod, trackID = NULL, mod = NULL)

Arguments

Details

Deprecated: please use viterbi with argument tod instead.

Value

vector of decoded states of length n

See Also

Other decoding functions: stateprobs(), stateprobs_g(), stateprobs_p(), viterbi(), viterbi_g()

Examples

delta = c(0.5, 0.5)
beta = matrix(c(-2, 1, -1,
                -2, -1, 1), nrow = 2, byrow = TRUE)
Gamma = tpm_p(1:24, 24, beta)

tod = rep(1:24, 5)
n = length(tod)

allprobs = matrix(runif(2*n), nrow = n, ncol = 2)
states = viterbi_p(delta, Gamma, allprobs, tod)

Von Mises distribution

Description

Density, distribution function and random generation for the von Mises distribution.

Usage

dvm(x, mu = 0, kappa = 1, log = FALSE)

pvm(q, mu = 0, kappa = 1, from = NULL, tol = 1e-20)

rvm(n, mu = 0, kappa = 1, wrap = TRUE)

Arguments

Details

This implementation of dvm allows for automatic differentiation with RTMB. rvm and pvm are simply wrappers of the corresponding functions from circular.

Value

dvm gives the density, pvm gives the distribution function, and rvm generates random deviates.

Examples

set.seed(1)
x = rvm(10, 0, 1)
d = dvm(x, 0, 1)
p = pvm(x, 0, 1)

Wrapped Cauchy distribution

Description

Density and random generation for the wrapped Cauchy distribution.

Usage

dwrpcauchy(x, mu = 0, rho, log = FALSE)

rwrpcauchy(n, mu = 0, rho, wrap = TRUE)

Arguments

Details

This implementation of dwrpcauchy allows for automatic differentiation with RTMB. rwrpcauchy is simply a wrapper for rwrappedcauchy imported from circular.

Value

dwrpcauchy gives the density and rwrpcauchy generates random deviates.

Examples

set.seed(1)
x = rwrpcauchy(10, 0, 0.3)
d = dwrpcauchy(x, 0, 0.3)

Zero-inflated density constructer

Description

Constructs a zero-inflated density function from a given probability density function

Usage

zero_inflate(dist, discrete = NULL)

Arguments

Details

The definition of zero-inflation is different for discrete and continuous distributions. For discrete distributions with p.m.f. f and zero-inflation probability p, we have

\Pr(X = 0) = p + (1 - p) \cdot f(0),

and

\Pr(X = x) = (1 - p) \cdot f(x), \quad x > 0.

For continuous distributions with p.d.f. f, we have

f_{\text{zinfl}}(x) = p \cdot \delta_0(x) + (1 - p) \cdot f(x),

where \delta_0 is the Dirac delta function at zero.

Value

zero-inflated density function with first argument x, second argument zeroprob, and additional arguments ... that will be passed to dist.

Examples

dzinorm <- zero_inflate(dnorm)
dzinorm(c(NA, 0, 2), 0.5, mean = 1, sd = 1)

zipois <- zero_inflate(dpois)
zipois(c(NA, 0, 1), 0.5, 1)