Package {hmmTMB}

hmmTMB: Fit Hidden Markov Models using Template Model Builder

Description

Fitting hidden Markov models using automatic differentiation and Laplace approximation, allowing for fast inference and flexible covariate effects (including random effects and smoothing splines) on model parameters. The package is described by Michelot (2025) doi:10.18637/jss.v114.i05.

Details

The package hmmTMB is based on three main classes (i.e., types of objects): MarkovChain, Observation, and HMM. Type ?MarkovChain, ?Observation or ?HMM to find their documentation, or consult the package vignettes for detailed examples of the hmmTMB workflow.

Author(s)

Maintainer: Theo Michelot theo.michelot@dal.ca

Authors:

• Theo Michelot theo.michelot@dal.ca

• Richard Glennie

Other contributors:

• Jan-Ole Fischer (SPDE smooths) [contributor]

• Duncan Cameron-Steinke (Forecast class) [contributor]

• Pauline Rocu (schisq distribution) [contributor]

• Litty Mathew (hnbinom distribution) [contributor]

See Also

Useful links:

• https://github.com/TheoMichelot/hmmTMB

R6 class for probability distribution

Description

Contains the probability density/mass function, and the link and inverse link functions for a probability distribution.

Methods

Public methods

• Dist$new()

• Dist$name()

• Dist$pdf()

• Dist$cdf()

• Dist$rng()

• Dist$link()

• Dist$invlink()

• Dist$npar()

• Dist$parnames()

• Dist$parapprox()

• Dist$fixed()

• Dist$code()

• Dist$name_long()

• Dist$set_npar()

• Dist$set_parnames()

• Dist$set_code()

• Dist$pdf_apply()

• Dist$rng_apply()

• Dist$par_alt()

• Dist$n2w()

• Dist$w2n()

• Dist$clone()

Dist$new()

Create a Dist object

Usage

Dist$new(
  name,
  pdf,
  rng,
  cdf = NULL,
  link,
  invlink,
  npar,
  parnames,
  parapprox = NULL,
  fixed = NULL,
  name_long = name,
  par_alt = NULL
)

Arguments

Returns

A new Dist object

Dist$name()

Return name of Dist object

Usage

Dist$name()

Dist$pdf()

Return pdf of Dist object

Usage

Dist$pdf()

Dist$cdf()

Return cdf of Dist object

Usage

Dist$cdf()

Dist$rng()

Return random generator function of Dist object

Usage

Dist$rng()

Dist$link()

Return link function of Dist object

Usage

Dist$link()

Dist$invlink()

Return inverse link function of Dist object

Usage

Dist$invlink()

Dist$npar()

Return number of parameters of Dist object

Usage

Dist$npar()

Dist$parnames()

Return names of parameters

Usage

Dist$parnames()

Dist$parapprox()

Return function that approximates parameters

Usage

Dist$parapprox()

Dist$fixed()

Return which parameters are fixed

Usage

Dist$fixed()

Dist$code()

Return code of Dist object

Usage

Dist$code()

Dist$name_long()

Human-readable name of Dist object

Usage

Dist$name_long()

Dist$set_npar()

Set number of parameters this distribution has

Usage

Dist$set_npar(new_npar)

Arguments

Dist$set_parnames()

Set parameter names

Usage

Dist$set_parnames(new_parnames)

Arguments

Dist$set_code()

Set distribution code

Usage

Dist$set_code(new_code)

Arguments

Dist$pdf_apply()

Evaluate probability density/mass function

This method is used in the Dist$obs_probs() method. It is a wrapper around Dist$pdf(), which prepares the parameters and passes them to the function.

Usage

Dist$pdf_apply(x, par, log = FALSE)

Arguments

Returns

Probability density/mass function evaluated at x for parameters par

Dist$rng_apply()

Random number generator

This method is a wrapper around Dist$rng(), which prepares the parameters and passes them to the function.

Usage

Dist$rng_apply(n, par)

Arguments

Returns

Vector of n realisations of this distribution

Dist$par_alt()

Alternative parameter formatting

Usage

Dist$par_alt(par)

Arguments

Returns

Formatted parameters

Dist$n2w()

Natural to working parameter transformation

This method transforms parameters from the natural scale (i.e., their domain of definition) to the "working" or "linear predictor" scale (i.e., the real line). It is a wrapper for Dist$link().

Usage

Dist$n2w(par)

Arguments

Returns

Vector of parameters on the working scale

Dist$w2n()

Working to natural parameter transformation

This method transforms parameters from the "working" or "linear predictor" scale (i.e., the real line) to the natural scale (i.e., their domain of definition). It is a wrapper for Dist$invlink().

Usage

Dist$w2n(wpar, as_matrix = FALSE)

Arguments

Returns

List or matrix of parameters on natural scale

Dist$clone()

The objects of this class are cloneable with this method.

Usage

Dist$clone(deep = FALSE)

Arguments

R6 Class for Forecast in Hidden Markov Models

Description

Creates a forecast object from a fitted HMM model, projecting future state probabilities and unconditional observation densities for a user-defined horizon using a covariate data-frame if applicable. Forecasting occurs in the constructor; results are accessible via public methods.

Methods

Public methods

• Forecast$new()

• Forecast$obs_par_forecast()

• Forecast$tpm_forecast()

• Forecast$eval_range()

• Forecast$forecast_data()

• Forecast$update_eval_range()

• Forecast$hidden_state_forecast()

• Forecast$forecast_dists()

• Forecast$quantile()

• Forecast$mean()

• Forecast$plot()

• Forecast$clone()

Forecast$new()

Construct a forecast object.

Usage

Forecast$new(
  hmm = NULL,
  n = NULL,
  forecast_data = NULL,
  preset_eval_range = NULL,
  starting_state_distribution = "last"
)

Arguments

Details

Algorithm overview

1. Validate the arguments.

2. Build or accept the future design matrix (forecast_data).

3. Obtain time-varying transition matrices (tpm_forecast) and observation parameters (obs_par_forecast) via HMM$predict().

4. Forward-propagate the hidden-state distribution (hidden_state_forecast).

5. Marginalise over states to obtain unconditional pdfs on the evaluation grid (forecast_dists).

Returns

A Forecast object. Results are read back with the accessor methods below, each invoked with ():

hidden_state_forecast

Matrix of nstates rows by n forecast steps.

obs_par_forecast

Array of time-varying observation parameters.

tpm_forecast

Array of time-varying transition matrices.

forecast_dists

List of unconditional pdf matrices, one per observation variable.

eval_range

List of grids on which each pdf was evaluated.

forecast_data

The data.frame of future covariates used for prediction.

The evaluation grid can be changed after construction with update_eval_range(), which recomputes forecast_dists.

Examples

data(faithful)

# Create training and forecast data frames
training_fraction <- 0.8
n_training <- floor(nrow(faithful) * training_fraction)
training_df <- data.frame(
  ID = rep(1, n_training),
  waiting = faithful$waiting[1:n_training],
  eruptions = faithful$eruptions[1:n_training]
)
forecast_df <- data.frame(
  ID = rep(1, nrow(faithful) - n_training),
  waiting = faithful$waiting[(n_training + 1):nrow(faithful)],
  eruptions = faithful$eruptions[(n_training + 1):nrow(faithful)]
)
# Create 2-state model with non-linear effect of waiting on all transitions
hid_model <- MarkovChain$new(
  data = training_df,
  n_states = 2,
  formula = ~ s(waiting, k = 10, bs = "cs")
)

# Create observation model with normal distribution for duration
obs_model <- Observation$new(
  data = training_df,
  n_states = 2,
  dists = list(eruptions = "norm"),
  par = list(eruptions = list(
    mean = c(0, 0),
    sd = c(1, 1)
  ))
)
# Update model parameters to suggested
obs_model$update_par(par = obs_model$suggest_initial())
# Create HMM model
hmm <- HMM$new(
  hid = hid_model,
  obs = obs_model
)
hmm$fit(silent = TRUE)
# Create forecast object. Setting the evaluation grid explicitly is
# recommended; omitting it falls back to a default grid and warns.
forecast <- Forecast$new(
  hmm = hmm,
  forecast_data = forecast_df,
  preset_eval_range = list(eruptions = seq(1, 6, length.out = 100)),
  starting_state_distribution = "last"
)

# Unconditional predictive densities: grid points x forecast steps
dim(forecast$forecast_dists()$eruptions)

Forecast$obs_par_forecast()

Get predicted observation parameters.

Usage

Forecast$obs_par_forecast()

Forecast$tpm_forecast()

Get predicted transition matrices.

Usage

Forecast$tpm_forecast()

Forecast$eval_range()

Get the evaluation grid for each response variable.

Usage

Forecast$eval_range()

Forecast$forecast_data()

Get the forecast data used for predictions.

Usage

Forecast$forecast_data()

Forecast$update_eval_range()

Update the evaluation grid for forecast pdfs. This method reconfigures the evaluation range used to compute unconditional predictive pdfs. It validates the new eval_range, updates internal state, reinitializes forecast densities, and returns the object invisibly for chaining.

Usage

Forecast$update_eval_range(eval_range)

Arguments

Returns

The Forecast object (invisibly), allowing method chaining.

Forecast$hidden_state_forecast()

Get the forward state probabilities.

Usage

Forecast$hidden_state_forecast()

Forecast$forecast_dists()

Get the unconditional predictive pdfs.

Usage

Forecast$forecast_dists()

Forecast$quantile()

Quantiles of forecast distribution

Usage

Forecast$quantile(var, probs)

Arguments

Returns

Data frame with one column for each quantile probability and one row for each time step in the forecast

Forecast$mean()

Mean of forecast distribution

Usage

Forecast$mean(var)

Arguments

Returns

Numeric vector of means, with one element for each time step in the forecast

Forecast$plot()

Plot forecast

Usage

Forecast$plot(var, show_data = FALSE)

Arguments

Details

The plot shows the mean of the forecast distribution over time (red line), as well as bands of decreasing opacity for the 50 90

Returns

A ggplot object

Forecast$clone()

The objects of this class are cloneable with this method.

Usage

Forecast$clone(deep = FALSE)

Arguments

See Also

HMM for the fitted model that is forecast from, and MarkovChain and Observation for the hidden state and observation components it is built from.

Examples

## ------------------------------------------------
## Method `Forecast$new()`
## ------------------------------------------------


data(faithful)

# Create training and forecast data frames
training_fraction <- 0.8
n_training <- floor(nrow(faithful) * training_fraction)
training_df <- data.frame(
  ID = rep(1, n_training),
  waiting = faithful$waiting[1:n_training],
  eruptions = faithful$eruptions[1:n_training]
)
forecast_df <- data.frame(
  ID = rep(1, nrow(faithful) - n_training),
  waiting = faithful$waiting[(n_training + 1):nrow(faithful)],
  eruptions = faithful$eruptions[(n_training + 1):nrow(faithful)]
)
# Create 2-state model with non-linear effect of waiting on all transitions
hid_model <- MarkovChain$new(
  data = training_df,
  n_states = 2,
  formula = ~ s(waiting, k = 10, bs = "cs")
)

# Create observation model with normal distribution for duration
obs_model <- Observation$new(
  data = training_df,
  n_states = 2,
  dists = list(eruptions = "norm"),
  par = list(eruptions = list(
    mean = c(0, 0),
    sd = c(1, 1)
  ))
)
# Update model parameters to suggested
obs_model$update_par(par = obs_model$suggest_initial())
# Create HMM model
hmm <- HMM$new(
  hid = hid_model,
  obs = obs_model
)
hmm$fit(silent = TRUE)
# Create forecast object. Setting the evaluation grid explicitly is
# recommended; omitting it falls back to a default grid and warns.
forecast <- Forecast$new(
  hmm = hmm,
  forecast_data = forecast_df,
  preset_eval_range = list(eruptions = seq(1, 6, length.out = 100)),
  starting_state_distribution = "last"
)

# Unconditional predictive densities: grid points x forecast steps
dim(forecast$forecast_dists()$eruptions)

R6 class for hidden Markov model

Description

Encapsulates the observation and hidden state models for a hidden Markov model.

Methods

Public methods

• HMM$new()

• HMM$obs()

• HMM$hid()

• HMM$out()

• HMM$tmb_obj()

• HMM$tmb_obj_joint()

• HMM$tmb_rep()

• HMM$states()

• HMM$coeff_fe()

• HMM$coeff_re()

• HMM$coeff_list()

• HMM$fixpar()

• HMM$coeff_array()

• HMM$lambda()

• HMM$bw()

• HMM$update_bw()

• HMM$check_bw()

• HMM$update_par()

• HMM$sd_re()

• HMM$par()

• HMM$set_priors()

• HMM$priors()

• HMM$iters()

• HMM$out_stan()

• HMM$llk()

• HMM$edf()

• HMM$suggest_initial()

• HMM$setup()

• HMM$fit_stan()

• HMM$fit()

• HMM$mle()

• HMM$forward_backward()

• HMM$pseudores()

• HMM$viterbi()

• HMM$sample_states()

• HMM$state_probs()

• HMM$post_coeff()

• HMM$post_linpred()

• HMM$post_fn()

• HMM$predict()

• HMM$confint()

• HMM$simulate()

• HMM$check()

• HMM$plot_ts()

• HMM$plot_dist()

• HMM$plot()

• HMM$plot_2d()

• HMM$AIC_marginal()

• HMM$AIC_conditional()

• HMM$print_obspar()

• HMM$print_tpm()

• HMM$formulation()

• HMM$print()

• HMM$clone()

HMM$new()

Create new HMM object

Usage

HMM$new(
  obs = NULL,
  hid = NULL,
  file = NULL,
  init = NULL,
  fixpar = NULL,
  bw = NULL
)

Arguments

Returns

A new HMM object

Examples

# Load data set (included with R)
data(nottem)
data <- data.frame(temp = as.vector(t(nottem)))

# Create hidden state and observation models
hid <- MarkovChain$new(data = data, n_states = 2)
par0 <- list(temp = list(mean = c(40, 60), sd = c(5, 5)))
obs <- Observation$new(data = data, n_states = 2,
                       dists = list(temp = "norm"),
                       par = par0)

# Create HMM
hmm <- HMM$new(hid = hid, obs = obs)

HMM$obs()

Observation object for this model

Usage

HMM$obs()

HMM$hid()

MarkovChain object for this model

Usage

HMM$hid()

HMM$out()

Output of optimiser after model fitting

Usage

HMM$out()

HMM$tmb_obj()

Model object created by TMB. This is the output of the TMB function MakeADFun, and it is a list including elements

fn

Objective function

gr

Gradient function of fn

par

Vector of initial parameters on working scale

Usage

HMM$tmb_obj()

HMM$tmb_obj_joint()

Model object created by TMB for the joint likelihood of the fixed and random effects. This is the output of the TMB function MakeADFun, and it is a list including elements

fn

Objective function

gr

Gradient function of fn

par

Vector of initial parameters on working scale

Usage

HMM$tmb_obj_joint()

HMM$tmb_rep()

Output of the TMB function sdreport, which includes estimates and standard errors for all model parameters.

Usage

HMM$tmb_rep()

HMM$states()

Vector of estimated states, after viterbi has been run

Usage

HMM$states()

HMM$coeff_fe()

Coefficients for fixed effect parameters

Usage

HMM$coeff_fe()

HMM$coeff_re()

Coefficients for random effect parameters

Usage

HMM$coeff_re()

HMM$coeff_list()

List of all model coefficients

These are the parameters estimated by the model, including fixed and random effect parameters for the observation parameters and the transition probabilities, (transformed) initial probabilities, and smoothness parameters.

Usage

HMM$coeff_list()

HMM$fixpar()

Fixed parameters

Usage

HMM$fixpar(all = FALSE)

Arguments

HMM$coeff_array()

Array of working parameters

Usage

HMM$coeff_array()

HMM$lambda()

Smoothness parameters

Usage

HMM$lambda()

HMM$bw()

Bandwidth of the forward algorithm (0 = exact)

Usage

HMM$bw()

HMM$update_bw()

Update bandwidth of the forward algorithm

The exact forward algorithm accumulates log-likelihood contributions that are each conditional on the entire process history, so the Hessian with respect to latent variables spread over time is dense. The banded algorithm of Fischer (2026) splits the time series into blocks of length bw and truncates that conditioning at the previous block, which makes the Hessian banded. It costs about twice as much per evaluation, and its error decays geometrically in bw, because an ergodic Markov chain forgets its initial condition exponentially fast.

A model containing a Gaussian field is banded by default, with bw = 15, because the dense Hessian of the exact algorithm defeats the sparsity such a field exists to provide. That default is a starting point, not an answer: the bandwidth needed depends on how fast the chain forgets its initial condition, so a persistent chain needs more. Every other model is exact by default. Use HMM$check_bw() to choose a bandwidth before fitting, and afterwards check that raising it leaves the estimates unchanged.

Usage

HMM$update_bw(bw = NULL)

Arguments

HMM$check_bw()

Profile the log-likelihood against the bandwidth

Evaluates the joint log-likelihood at the model's current parameter values for a range of bandwidths. A sensible bandwidth is the smallest one beyond which the value no longer changes appreciably. This is evaluated at fixed parameters, so it is cheap and does not require the model to be re-fitted.

Usage

HMM$check_bw(bws = seq(5, 40, by = 5), silent = TRUE)

Arguments

Returns

Data frame with columns bw and llk, including a row with bw = Inf for the exact forward algorithm.

HMM$update_par()

Update parameters stored inside model object

Usage

HMM$update_par(par_list = NULL, iter = NULL)

Arguments

HMM$sd_re()

Standard deviation of smooth terms (or random effects)

This function transforms the smoothness parameter of each smooth term into a standard deviation, given by SD = 1/sqrt(lambda). It is particularly helpful to get the standard deviations of independent normal random effects.

Usage

HMM$sd_re()

Returns

List of standard deviations for observation model and hidden state model.

HMM$par()

Model parameters

Usage

HMM$par(t = 1)

Arguments

Returns

A list with elements:

obspar

Parameters of observation model

tpm

Transition probability matrix of hidden state model

HMM$set_priors()

Set priors for coefficients

Usage

HMM$set_priors(new_priors = NULL)

Arguments

HMM$priors()

Extract stored priors

Usage

HMM$priors()

HMM$iters()

Iterations from stan MCMC fit

Usage

HMM$iters(type = "response")

Arguments

Returns

see output of as.matrix in stan

HMM$out_stan()

fitted stan object from MCMC fit

Usage

HMM$out_stan()

Returns

the stanfit object

HMM$llk()

Log-likelihood at current parameters

Usage

HMM$llk()

Returns

Log-likelihood

HMM$edf()

Effective degrees of freedom

Usage

HMM$edf()

Returns

Number of effective degrees of freedom (accounting for flexibility in non-parametric terms implied by smoothing)

HMM$suggest_initial()

Suggest initial parameter values

Uses K-means clustering to split the data into naive "states", and estimates observation parameters within each of these states. This is meant to help with selecting initial parameter values before model fitting, but users should still think about the values carefully, and try multiple set of initial parameter values to ensure convergence to the global maximum of the likelihood function.

Usage

HMM$suggest_initial()

Returns

List of initial parameters

HMM$setup()

TMB setup

This creates an attribute tmb_obj, which can be used to evaluate the negative log-likelihood function.

Usage

HMM$setup(silent = TRUE)

Arguments

HMM$fit_stan()

Fit model using tmbstan

Consult documentation of the tmbstan package for more information. After this method has been called, the Stan output can be accessed using the method out_stan(). This Stan output can for example be visualised using functions from the rstan package. The parameters stored in this HMM object are automatically updated to the mean posterior estimate, although this can be changed using update_par().

Usage

HMM$fit_stan(..., silent = FALSE)

Arguments

HMM$fit()

Model fitting

The negative log-likelihood of the model is minimised using the function nlminb(). TMB uses the Laplace approximation to integrate the random effects out of the likelihood.

After the model has been fitted, the output of nlminb() can be accessed using the method out(). The estimated parameters can be accessed using the methods par() (for the HMM parameters, possibly dependent on covariates), predict() (for uncertainty quantification and prediction of the HMM parameters for new covariate values), coeff_fe() (for estimated fixed effect coefficients on the linear predictor scale), and coeff_re() (for estimated random effect coefficients on the linear predictor scale).

Usage

HMM$fit(silent = FALSE, ...)

Arguments

Examples

# Load data set (included with R)
data(nottem)
data <- data.frame(temp = as.vector(t(nottem)))

# Create hidden state and observation models
hid <- MarkovChain$new(data = data, n_states = 2)
par0 <- list(temp = list(mean = c(40, 60), sd = c(5, 5)))
obs <- Observation$new(data = data, n_states = 2,
                       dists = list(temp = "norm"),
                       par = par0)

# Create HMM
hmm <- HMM$new(hid = hid, obs = obs)

# Fit HMM
hmm$fit(silent = TRUE)

HMM$mle()

Get maximum likelihood estimates once model fitted

Usage

HMM$mle()

Returns

list of maximum likelihood estimates as described as input for the function update_par()

HMM$forward_backward()

Forward-backward algorithm

The forward probability for time step t and state j is the joint pdf/pmf of observations up to time t and of being in state j at time t, p(Z[1], Z[2], ..., Z[t], S[t] = j). The backward probability for time t and state j is the conditional pdf/pmf of observations between time t + 1 and n, given state j at time t, p(Z[t+1], Z[t+2], ..., Z[n] | S[t] = j). This function returns their logarithm, for use in other methods state_probs, and sample_states.

Usage

HMM$forward_backward()

Returns

Log-forward and log-backward probabilities

HMM$pseudores()

Pseudo-residuals

Compute pseudo-residuals for the fitted model. If the fitted model is the "true" model, the pseudo-residuals follow a standard normal distribution. Deviations from normality suggest lack of fit.

Usage

HMM$pseudores()

Returns

List (of length the number of variables), where each element is a vector of pseudo-residuals (of length the number of data points)

HMM$viterbi()

Viterbi algorithm

Usage

HMM$viterbi()

Returns

Most likely state sequence

HMM$sample_states()

Sample posterior state sequences using forward-filtering backward-sampling

The forward-filtering backward-sampling algorithm returns a sequence of states, similarly to the Viterbi algorithm, but it generates it from the posterior distribution of state sequences, i.e., accounting for uncertainty in the state classification. Multiple generated sequences will therefore generally not be the same.

Usage

HMM$sample_states(nsamp = 1, full = FALSE)

Arguments

Returns

Matrix where each column is a different sample of state sequences, and each row is a time of observation

HMM$state_probs()

Compute posterior probability of being in each state

Usage

HMM$state_probs()

Returns

matrix with a row for each observation and a column for each state

HMM$post_coeff()

Posterior sampling for model coefficients

Usage

HMM$post_coeff(n_post)

Arguments

Returns

Matrix with one column for each coefficient and one row for each posterior draw

HMM$post_linpred()

Posterior sampling for linear predictor

Usage

HMM$post_linpred(n_post)

Arguments

Returns

List with elements obs and hid, where each is a matrix with one column for each predictor and one row for each posterior draw

HMM$post_fn()

Create posterior simulations of a function of a model component

Usage

HMM$post_fn(fn, n_post, comp = NULL, ..., level = 0, return_post = FALSE)

Arguments

Returns

A list with elements:

post

If return_post = TRUE, this is a vector (for scalar outputs of fn) or matrix (for vector outputs) with a column for each simulation

mean

Mean over posterior samples

lcl

Lower confidence interval bound (if level !=0)

ucl

Upper confidence interval bound (if level !=0)

HMM$predict()

Predict estimates from a fitted model

By default, this returns point estimates of the HMM parameters for a new data frame of covariates. See the argument 'n_post' to also get confidence intervals.

Usage

HMM$predict(
  what,
  t = 1,
  newdata = NULL,
  n_post = 0,
  level = 0.95,
  return_post = FALSE,
  as_list = TRUE
)

Arguments

Returns

Maximum likelihood estimates (mle) of predictions, and confidence limits (lcl and ucl) if requested. The format of the output depends on whether confidence intervals are required (specified through n_post), and on the argument as_list.

Examples

# Load data set (included with R)
data(nottem)
data <- data.frame(temp = as.vector(t(nottem)))

# Create hidden state and observation models
hid <- MarkovChain$new(data = data, n_states = 2)
par0 <- list(temp = list(mean = c(40, 60), sd = c(5, 5)))
obs <- Observation$new(data = data, n_states = 2,
                       dists = list(temp = "norm"),
                       par = par0)

# Create HMM
hmm <- HMM$new(hid = hid, obs = obs)

# Fit HMM
hmm$fit(silent = TRUE)

# Get transition probability matrix with confidence intervals
hmm$predict(what = "tpm", n_post = 1000)

HMM$confint()

Confidence intervals for working parameters

This function computes standard errors for all fixed effect model parameters based on the diagonal of the inverse of the Hessian matrix, and then derives Wald-type confidence intervals.

Usage

HMM$confint(level = 0.95)

Arguments

Returns

List of matrices with three columns: mle (maximum likelihood estimate), lcl (lower confidence limit), ucl (upper confidence limit), and se (standard error). One such matrix is produced for the working parameters of the observation model, the working parameters of the hidden state model, the smoothness parameters of the observation model, and the smoothness parameters of the hidden state model.

HMM$simulate()

Simulate from hidden Markov model

Usage

HMM$simulate(n, data = NULL, silent = FALSE)

Arguments

Returns

Data frame including columns of data (if provided), and simulated data variables

HMM$check()

Compute goodness-of-fit statistics using simulation

Many time series are simulated from the fitted model, and the statistic(s) of interest are calculated for each. A histogram of those values can for example be used to compare to the observed value of the statistic. An observation far in the tails of the distribution of simulated statistics suggests lack of fit.

Usage

HMM$check(check_fn, nsims = 100, full = FALSE, silent = FALSE)

Arguments

Returns

List with elements:

obs_stat:

Vector of values of goodness-of-fit statistics for the observed data

stats:

Matrix of values of goodness-of-fit statistics for the simulated data sets (one row for each statistic, and one column for each simulation)

plot:

ggplot object

HMM$plot_ts()

Time series plot coloured by states

Creates a plot of the data coloured by the most likely state sequence, as estimated by the Viterbi algorithm. If one variable name is passed as input, it is plotted against time. If two variables are passed, they are plotted against each other.

Usage

HMM$plot_ts(var, var2 = NULL, line = TRUE)

Arguments

Returns

A ggplot object

HMM$plot_dist()

Plot observation distributions weighted by frequency in Viterbi

This is a wrapper around Observation$plot_dist, where the distribution for each state is weighted by the proportion of time spent in that state (according to the Viterbi state sequence).

Usage

HMM$plot_dist(var = NULL)

Arguments

Returns

Plot of distributions with data histogram

HMM$plot()

Plot a model component

Usage

HMM$plot(
  what,
  var = NULL,
  covs = NULL,
  i = NULL,
  j = NULL,
  n_grid = 50,
  n_post = 1000
)

Arguments

Returns

A ggplot object

HMM$plot_2d()

Plot a model component over two covariates

The two-dimensional counterpart of HMM$plot(), for a model containing a bivariate smooth such as s(x, y) or a Gaussian field. HMM$plot() varies one covariate and holds the rest fixed, which for a bivariate term shows a single slice through the surface; this varies both and draws the surface itself.

Nothing here is specific to any basis. The surface is built from HMM$predict() on a lattice of covariate values, so it works for any smoother hmmTMB accepts on two covariates – a thin plate spline, bs = "gp", an SPDE field – and for two covariates that enter the model separately, where it shows their combined effect.

Usage

HMM$plot_2d(
  what,
  var,
  var2,
  covs = NULL,
  i = NULL,
  j = NULL,
  n_grid = 50,
  n_post = 0,
  level = 0.95,
  show = "mle",
  too_far = 0.1,
  contour = TRUE,
  padding = 0
)

Arguments

Returns

A ggplot object

HMM$AIC_marginal()

Marginal Akaike Information Criterion

The marginal AIC is for example defined by Wood (2017), as AIC = - 2L + 2k where L is the maximum marginal log-likelihood (of fixed effects), and k is the number of degrees of freedom of the fixed effect component of the model

Usage

HMM$AIC_marginal()

Returns

Marginal AIC

HMM$AIC_conditional()

Conditional Akaike Information Criterion

The conditional AIC is for example defined by Wood (2017), as AIC = - 2L + 2k where L is the maximum joint log-likelihood (of fixed and random effects), and k is the number of effective degrees of freedom of the model (accounting for flexibility in non-parametric terms implied by smoothing)

Usage

HMM$AIC_conditional()

Returns

Conditional AIC

HMM$print_obspar()

Print observation parameters at t = 1

Usage

HMM$print_obspar()

HMM$print_tpm()

Print observation parameters at t = 1

Usage

HMM$print_tpm()

HMM$formulation()

Print model formulation

Usage

HMM$formulation()

HMM$print()

Print HMM object

Usage

HMM$print()

HMM$clone()

The objects of this class are cloneable with this method.

Usage

HMM$clone(deep = FALSE)

Arguments

Examples

## ------------------------------------------------
## Method `HMM$new()`
## ------------------------------------------------

# Load data set (included with R)
data(nottem)
data <- data.frame(temp = as.vector(t(nottem)))

# Create hidden state and observation models
hid <- MarkovChain$new(data = data, n_states = 2)
par0 <- list(temp = list(mean = c(40, 60), sd = c(5, 5)))
obs <- Observation$new(data = data, n_states = 2,
                       dists = list(temp = "norm"),
                       par = par0)

# Create HMM
hmm <- HMM$new(hid = hid, obs = obs)

## ------------------------------------------------
## Method `HMM$fit()`
## ------------------------------------------------

# Load data set (included with R)
data(nottem)
data <- data.frame(temp = as.vector(t(nottem)))

# Create hidden state and observation models
hid <- MarkovChain$new(data = data, n_states = 2)
par0 <- list(temp = list(mean = c(40, 60), sd = c(5, 5)))
obs <- Observation$new(data = data, n_states = 2,
                       dists = list(temp = "norm"),
                       par = par0)

# Create HMM
hmm <- HMM$new(hid = hid, obs = obs)

# Fit HMM
hmm$fit(silent = TRUE)

## ------------------------------------------------
## Method `HMM$predict()`
## ------------------------------------------------

# Load data set (included with R)
data(nottem)
data <- data.frame(temp = as.vector(t(nottem)))

# Create hidden state and observation models
hid <- MarkovChain$new(data = data, n_states = 2)
par0 <- list(temp = list(mean = c(40, 60), sd = c(5, 5)))
obs <- Observation$new(data = data, n_states = 2,
                       dists = list(temp = "norm"),
                       par = par0)

# Create HMM
hmm <- HMM$new(hid = hid, obs = obs)

# Fit HMM
hmm$fit(silent = TRUE)

# Get transition probability matrix with confidence intervals
hmm$predict(what = "tpm", n_post = 1000)

R6 class for HMM hidden process model

Description

Contains the parameters and model formulas for the hidden process model.

Methods

Public methods

• MarkovChain$new()

• MarkovChain$formula()

• MarkovChain$formulas()

• MarkovChain$tpm()

• MarkovChain$ref()

• MarkovChain$ref_mat()

• MarkovChain$ref_delta0()

• MarkovChain$coeff_fe()

• MarkovChain$delta()

• MarkovChain$delta0()

• MarkovChain$stationary()

• MarkovChain$fixpar()

• MarkovChain$coeff_re()

• MarkovChain$X_fe()

• MarkovChain$X_re()

• MarkovChain$lambda()

• MarkovChain$sd_re()

• MarkovChain$nstates()

• MarkovChain$terms()

• MarkovChain$unique_ID()

• MarkovChain$initial_state()

• MarkovChain$empty()

• MarkovChain$gam_args()

• MarkovChain$update_tpm()

• MarkovChain$update_coeff_fe()

• MarkovChain$update_coeff_re()

• MarkovChain$update_X_fe()

• MarkovChain$update_X_re()

• MarkovChain$update_delta0()

• MarkovChain$update_lambda()

• MarkovChain$update_fixpar()

• MarkovChain$make_mat()

• MarkovChain$make_mat_grid()

• MarkovChain$tpm2par()

• MarkovChain$par2tpm()

• MarkovChain$linpred()

• MarkovChain$simulate()

• MarkovChain$formulation()

• MarkovChain$print()

• MarkovChain$clone()

MarkovChain$new()

Create new MarkovChain object

Usage

MarkovChain$new(
  data = NULL,
  formula = NULL,
  n_states,
  tpm = NULL,
  initial_state = "estimated",
  fixpar = NULL,
  ref = 1:n_states,
  gam_args = NULL
)

Arguments

Returns

A new MarkovChain object

Examples

# Load data set from MSwM package
data(energy, package = "MSwM")

# Create 2-state covariate-free model and initialise transition 
# probability matrix
hid <- MarkovChain$new(data = energy, n_states = 2,
                       tpm = matrix(c(0.8, 0.3, 0.2, 0.7), 2, 2))

# Create 2-state model with non-linear effect of Oil on all transition 
# probabilities
hid <- MarkovChain$new(data = energy, n_states = 2,
                       formula = ~ s(Oil, k = 5, bs = "cs"))

# Create 2-state model with quadratic effect of Oil on Pr(1 > 2)
structure <- matrix(c(".", "~poly(Oil, 2)",
                      "~1", "."),
                    ncol = 2, byrow = TRUE)
hid <- MarkovChain$new(data = energy, n_states = 2,
                       formula = structure)

MarkovChain$formula()

Formula of MarkovChain model

Usage

MarkovChain$formula()

MarkovChain$formulas()

List of formulas for MarkovChain model

Usage

MarkovChain$formulas()

MarkovChain$tpm()

Get transition probability matrices

Usage

MarkovChain$tpm(t = 1, linpred = NULL)

Arguments

Returns

Array with one slice for each transition probability matrix

MarkovChain$ref()

Indices of reference elements in transition probability matrix

Usage

MarkovChain$ref()

MarkovChain$ref_mat()

Matrix of reference elements in transition probability matrix

Usage

MarkovChain$ref_mat()

MarkovChain$ref_delta0()

Indices of reference elements in initial distribution

Usage

MarkovChain$ref_delta0()

MarkovChain$coeff_fe()

Current parameter estimates (fixed effects)

Usage

MarkovChain$coeff_fe()

MarkovChain$delta()

Stationary distribution

Usage

MarkovChain$delta(t = NULL, linpred = NULL)

Arguments

Returns

Matrix of stationary distributions. Each row corresponds to a row of the design matrices, and each column corresponds to a state.

MarkovChain$delta0()

Initial distribution

Usage

MarkovChain$delta0(log = FALSE, as_matrix = TRUE)

Arguments

Returns

Matrix with one row for each time series ID, and one column for each state. For each ID, the i-th element of the corresponding row is the probability Pr(S[1] = i)

MarkovChain$stationary()

Use stationary distribution as initial distribution?

Usage

MarkovChain$stationary()

MarkovChain$fixpar()

Fixed parameters

Usage

MarkovChain$fixpar(all = FALSE)

Arguments

MarkovChain$coeff_re()

Current parameter estimates (random effects)

Usage

MarkovChain$coeff_re()

MarkovChain$X_fe()

Fixed effect design matrix

Usage

MarkovChain$X_fe()

MarkovChain$X_re()

Random effect design matrix

Usage

MarkovChain$X_re()

MarkovChain$lambda()

Smoothness parameters

Usage

MarkovChain$lambda()

MarkovChain$sd_re()

Standard deviation of smooth terms

This function transforms the smoothness parameter of each smooth term into a standard deviation, given by SD = 1/sqrt(lambda). It is particularly helpful to get the standard deviations of independent normal random effects. A smooth with its own parameterisation rather than a single smoothness parameter – a Gaussian field, whose parameters are already a marginal standard deviation and a range – has no such transformation and returns NA; read those off lambda() directly.

Usage

MarkovChain$sd_re()

MarkovChain$nstates()

Number of states

Usage

MarkovChain$nstates()

MarkovChain$terms()

Terms of model formulas

Usage

MarkovChain$terms()

MarkovChain$unique_ID()

Number of time series

Usage

MarkovChain$unique_ID()

MarkovChain$initial_state()

Initial state (see constructor argument)

Usage

MarkovChain$initial_state()

MarkovChain$empty()

Empty model? (for simulation only)

Usage

MarkovChain$empty()

MarkovChain$gam_args()

Extra arguments for mgcv::gam (passed to make_matrices)

Usage

MarkovChain$gam_args()

MarkovChain$update_tpm()

Update transition probability matrix

Usage

MarkovChain$update_tpm(tpm)

Arguments

MarkovChain$update_coeff_fe()

Update coefficients for fixed effect parameters

Usage

MarkovChain$update_coeff_fe(coeff_fe)

Arguments

MarkovChain$update_coeff_re()

Update coefficients for random effect parameters

Usage

MarkovChain$update_coeff_re(coeff_re)

Arguments

MarkovChain$update_X_fe()

Update design matrix for fixed effects

Usage

MarkovChain$update_X_fe(X_fe)

Arguments

MarkovChain$update_X_re()

Update design matrix for random effects

Usage

MarkovChain$update_X_re(X_re)

Arguments

MarkovChain$update_delta0()

Update initial distribution

Usage

MarkovChain$update_delta0(delta0)

Arguments

MarkovChain$update_lambda()

Update smoothness parameters

Usage

MarkovChain$update_lambda(lambda)

Arguments

MarkovChain$update_fixpar()

Update information about fixed parameters

Usage

MarkovChain$update_fixpar(fixpar)

Arguments

MarkovChain$make_mat()

Make model matrices

Usage

MarkovChain$make_mat(data, new_data = NULL)

Arguments

Returns

A list with elements:

X_fe

Design matrix for fixed effects

X_re

Design matrix for random effects

S

Smoothness matrix for random effects

ncol_fe

Number of columns of X_fe for each parameter

ncol_re

Number of columns of X_re and S for each random effect

MarkovChain$make_mat_grid()

Design matrices for grid of covariates

Used in plotting functions such as HMM$plot_tpm and HMM$plot_stat_dist

Usage

MarkovChain$make_mat_grid(var, data, covs = NULL, n_grid = 1000)

Arguments

Returns

A list with the same elements as the output of make_mat, plus a data frame of covariates values.

MarkovChain$tpm2par()

Transform transition probabilities to working scale

Apply the multinomial logit link function to get the corresponding parameters on the working scale (i.e., linear predictor scale).

Usage

MarkovChain$tpm2par(tpm)

Arguments

Returns

Vector of parameters on linear predictor scale

MarkovChain$par2tpm()

Transform working parameters to transition probabilities

Apply the inverse multinomial logit link function to transform the parameters on the working scale (i.e., linear predictor scale) into the transition probabilities.

Usage

MarkovChain$par2tpm(par)

Arguments

Returns

Transition probability matrix

MarkovChain$linpred()

Linear predictor for transition probabilities

Usage

MarkovChain$linpred()

MarkovChain$simulate()

Simulate from Markov chain

Usage

MarkovChain$simulate(n, data = NULL, new_data = NULL, silent = FALSE)

Arguments

Returns

Sequence of states of simulated chain

MarkovChain$formulation()

Print model formulation

Usage

MarkovChain$formulation()

MarkovChain$print()

Print MarkovChain object

Usage

MarkovChain$print()

MarkovChain$clone()

The objects of this class are cloneable with this method.

Usage

MarkovChain$clone(deep = FALSE)

Arguments

Examples

## ------------------------------------------------
## Method `MarkovChain$new()`
## ------------------------------------------------

# Load data set from MSwM package
data(energy, package = "MSwM")

# Create 2-state covariate-free model and initialise transition 
# probability matrix
hid <- MarkovChain$new(data = energy, n_states = 2,
                       tpm = matrix(c(0.8, 0.3, 0.2, 0.7), 2, 2))

# Create 2-state model with non-linear effect of Oil on all transition 
# probabilities
hid <- MarkovChain$new(data = energy, n_states = 2,
                       formula = ~ s(Oil, k = 5, bs = "cs"))

# Create 2-state model with quadratic effect of Oil on Pr(1 > 2)
structure <- matrix(c(".", "~poly(Oil, 2)",
                      "~1", "."),
                    ncol = 2, byrow = TRUE)
hid <- MarkovChain$new(data = energy, n_states = 2,
                       formula = structure)

R6 class for HMM observation model

Description

Contains the data, distributions, parameters, and formulas for the observation model from a hidden Markov model.

Methods

Public methods

• Observation$new()

• Observation$data()

• Observation$dists()

• Observation$nstates()

• Observation$par()

• Observation$par_alt()

• Observation$inipar()

• Observation$coeff_fe()

• Observation$coeff_re()

• Observation$X_fe()

• Observation$X_re()

• Observation$lambda()

• Observation$sd_re()

• Observation$formulas()

• Observation$terms()

• Observation$obs_var()

• Observation$known_states()

• Observation$fixpar()

• Observation$empty()

• Observation$gam_args()

• Observation$update_par()

• Observation$update_coeff_fe()

• Observation$update_coeff_re()

• Observation$update_X_fe()

• Observation$update_X_re()

• Observation$update_lambda()

• Observation$update_data()

• Observation$update_fixpar()

• Observation$make_mat()

• Observation$make_newdata_grid()

• Observation$n2w()

• Observation$w2n()

• Observation$linpred()

• Observation$obs_probs()

• Observation$cdf()

• Observation$suggest_initial()

• Observation$plot_dist()

• Observation$formulation()

• Observation$print()

• Observation$clone()

Observation$new()

Create new Observation object

Usage

Observation$new(
  data = NULL,
  dists,
  formulas = NULL,
  n_states = NULL,
  par,
  fixpar = NULL,
  gam_args = NULL
)

Arguments

Returns

A new Observation object

Examples

# Load data set from MSwM package
data(energy, package = "MSwM")

# Initial observation parameters
par0 <- list(Price = list(mean = c(3, 6), sd = c(2, 2)))

# Model "energy" with normal distributions
obs <- Observation$new(data = energy, 
                       dists = list(Price = "norm"),
                       par = par0)
                       
# Model "energy" with gamma distributions
obs <- Observation$new(data = energy, 
                       dists = list(Price = "gamma2"),
                       par = par0)
                       
# Model with non-linear effect of EurDol on mean price
f <- list(Price = list(mean = ~ s(EurDol, k = 5, bs = "cs")))
obs <- Observation$new(data = energy, 
                       dists = list(Price = "norm"),
                       par = par0, 
                       formula = f)

Observation$data()

Data frame

Usage

Observation$data()

Observation$dists()

List of distributions

Usage

Observation$dists()

Observation$nstates()

Number of states

Usage

Observation$nstates()

Observation$par()

Parameters on natural scale

Usage

Observation$par(t = 1, full_names = TRUE, linpred = NULL, as_list = FALSE)

Arguments

Returns

Array of parameters with one row for each observation parameter, one column for each state, and one slice for each time step. (See as_list argument for alternative output format.)

Examples

# Load data set from MSwM package
data(energy, package = "MSwM")

# Initial observation parameters
par0 <- list(Price = list(mean = c(3, 6), sd = c(2, 2)))

# Model with linear effect of EurDol on mean price
f <- list(Price = list(mean = ~ EurDol))
obs <- Observation$new(data = energy, 
                       dists = list(Price = "norm"),
                       par = par0,
                       n_states = 2, 
                       formula = f)

# Set slope coefficients
obs$update_coeff_fe(coeff_fe = c(3, 2, 6, -2, log(2), log(2)))

# Observation parameter values for given data rows
obs$par(t = c(1, 10, 20))

Observation$par_alt()

Alternative parameter output

This function is only useful for the categorical and multivariate normal distributions, and it formats the parameters in a slightly nicer way.

Usage

Observation$par_alt(var = NULL, t = 1)

Arguments

Returns

List of distribution parameters, with one element for each state

Observation$inipar()

Return initial parameter values supplied

Usage

Observation$inipar()

Observation$coeff_fe()

Fixed effect parameters on working scale

Usage

Observation$coeff_fe()

Observation$coeff_re()

Random effect parameters

Usage

Observation$coeff_re()

Observation$X_fe()

Fixed effect design matrix

Usage

Observation$X_fe()

Observation$X_re()

Random effect design matrix

Usage

Observation$X_re()

Observation$lambda()

Smoothness parameters

Usage

Observation$lambda()

Observation$sd_re()

Standard deviation of smooth terms

This function transforms the smoothness parameter of each smooth term into a standard deviation, given by SD = 1/sqrt(lambda). It is particularly helpful to get the standard deviations of independent normal random effects. A smooth with its own parameterisation rather than a single smoothness parameter – a Gaussian field, whose parameters are already a marginal standard deviation and a range – has no such transformation and returns NA; read those off lambda() directly.

Usage

Observation$sd_re()

Observation$formulas()

List of model formulas for observation model

Usage

Observation$formulas(raw = FALSE)

Arguments

Observation$terms()

Terms of model formulas

Usage

Observation$terms()

Observation$obs_var()

Data frame of response variables

Usage

Observation$obs_var(expand = FALSE)

Arguments

Returns

Data frame of observation variables

Observation$known_states()

Vector of known states

Usage

Observation$known_states(mat = TRUE)

Arguments

Observation$fixpar()

Fixed parameters

Usage

Observation$fixpar(all = FALSE)

Arguments

Observation$empty()

Empty model? (for simulation only)

Usage

Observation$empty()

Observation$gam_args()

Extra arguments for mgcv::gam (passed to make_matrices)

Usage

Observation$gam_args()

Observation$update_par()

Update parameters

Updates the 'par' attribute to the list passed as input, and updates the intercept elements of 'coeff_fe' using the list passed as input

Usage

Observation$update_par(par)

Arguments

Observation$update_coeff_fe()

Update coefficients for fixed effect parameters

Usage

Observation$update_coeff_fe(coeff_fe)

Arguments

Observation$update_coeff_re()

Update random effect parameters

Usage

Observation$update_coeff_re(coeff_re)

Arguments

Observation$update_X_fe()

Update fixed effect design matrix

Usage

Observation$update_X_fe(X_fe)

Arguments

Observation$update_X_re()

Update random effect design matrix

Usage

Observation$update_X_re(X_re)

Arguments

Observation$update_lambda()

Update smoothness parameters

Usage

Observation$update_lambda(lambda)

Arguments

Observation$update_data()

Update data

Usage

Observation$update_data(data)

Arguments

Observation$update_fixpar()

Update information about fixed parameters

Usage

Observation$update_fixpar(fixpar)

Arguments

Observation$make_mat()

Make model matrices

Usage

Observation$make_mat(new_data = NULL)

Arguments

Returns

A list with elements:

X_fe

Design matrix for fixed effects

X_re

Design matrix for random effects

S

Smoothness matrix for random effects

ncol_fe

Number of columns of X_fe for each parameter

ncol_re

Number of columns of X_re and S for each random effect

Observation$make_newdata_grid()

Design matrices for grid of covariates

Usage

Observation$make_newdata_grid(var, covs = NULL, n_grid = 1000)

Arguments

Returns

A list with the same elements as the output of make_mat, plus a data frame of covariates values.

Observation$n2w()

Natural to working parameter transformation

This function applies the link functions of the distribution parameters, to transform parameters from their natural scale to the working scale (i.e., linear predictor scale)

Usage

Observation$n2w(par)

Arguments

Returns

Vector of parameters on working scale

Observation$w2n()

Working to natural parameter transformation

This function applies the inverse link functions of the distribution parameters, to transform parameters from the working scale (i.e., linear predictor scale) to their natural scale.

Usage

Observation$w2n(wpar)

Arguments

Returns

List of parameters on natural scale

Observation$linpred()

Compute linear predictor

Usage

Observation$linpred()

Observation$obs_probs()

Observation likelihoods

Usage

Observation$obs_probs(data = NULL)

Arguments

Returns

Matrix of likelihoods of observations, with one row for each time step, and one column for each state.

Observation$cdf()

Cumulative probabilities of observations

Usage

Observation$cdf()

Returns

List of cumulative probabilities, with one element for each observed variable. Matrix rows correspond to time steps, and columns correspond to states.

Observation$suggest_initial()

Suggest initial observation parameters

The K-means algorithm is used to define clusters of observations (supposed to approximate the HMM states). Then, for each cluster, the parapprox function of the relevant Dist object is used to obtain parameter values.

Usage

Observation$suggest_initial()

Returns

List of initial parameters for each observation variable

Examples

# Load data set from MSwM package
data(energy, package = "MSwM")

# Initial observation parameters
par0 <- list(Price = list(mean = c(3, 6), sd = c(2, 2)))

# Model "energy" with normal distributions
obs <- Observation$new(data = energy, 
                       dists = list(Price = "norm"),
                       par = par0,
                       n_states = 2)

# Print observation parameters
obs$par()

# Suggest initial parameters
par0_new <- obs$suggest_initial()
par0_new

# Update model parameters to suggested
obs$update_par(par = par0_new)
obs$par()

Observation$plot_dist()

Plot histogram of data and pdfs

Plot histogram of observations for the variable specified by the argument name, overlaid with the pdf of the specified distribution for that data stream. Helpful to select initial parameter values for model fitting, or to visualise fitted state-dependent distributions.

Usage

Observation$plot_dist(var = NULL, weights = NULL, t = 1)

Arguments

Returns

A ggplot object

Observation$formulation()

Print model formulation

Usage

Observation$formulation()

Observation$print()

Print Observation object Check constructor arguments

Usage

Observation$print()

Observation$clone()

The objects of this class are cloneable with this method.

Usage

Observation$clone(deep = FALSE)

Arguments

Examples

## ------------------------------------------------
## Method `Observation$new()`
## ------------------------------------------------

# Load data set from MSwM package
data(energy, package = "MSwM")

# Initial observation parameters
par0 <- list(Price = list(mean = c(3, 6), sd = c(2, 2)))

# Model "energy" with normal distributions
obs <- Observation$new(data = energy, 
                       dists = list(Price = "norm"),
                       par = par0)
                       
# Model "energy" with gamma distributions
obs <- Observation$new(data = energy, 
                       dists = list(Price = "gamma2"),
                       par = par0)
                       
# Model with non-linear effect of EurDol on mean price
f <- list(Price = list(mean = ~ s(EurDol, k = 5, bs = "cs")))
obs <- Observation$new(data = energy, 
                       dists = list(Price = "norm"),
                       par = par0, 
                       formula = f)

## ------------------------------------------------
## Method `Observation$par()`
## ------------------------------------------------

# Load data set from MSwM package
data(energy, package = "MSwM")

# Initial observation parameters
par0 <- list(Price = list(mean = c(3, 6), sd = c(2, 2)))

# Model with linear effect of EurDol on mean price
f <- list(Price = list(mean = ~ EurDol))
obs <- Observation$new(data = energy, 
                       dists = list(Price = "norm"),
                       par = par0,
                       n_states = 2, 
                       formula = f)

# Set slope coefficients
obs$update_coeff_fe(coeff_fe = c(3, 2, 6, -2, log(2), log(2)))

# Observation parameter values for given data rows
obs$par(t = c(1, 10, 20))

## ------------------------------------------------
## Method `Observation$suggest_initial()`
## ------------------------------------------------

# Load data set from MSwM package
data(energy, package = "MSwM")

# Initial observation parameters
par0 <- list(Price = list(mean = c(3, 6), sd = c(2, 2)))

# Model "energy" with normal distributions
obs <- Observation$new(data = energy, 
                       dists = list(Price = "norm"),
                       par = par0,
                       n_states = 2)

# Print observation parameters
obs$par()

# Suggest initial parameters
par0_new <- obs$suggest_initial()
par0_new

# Update model parameters to suggested
obs$update_par(par = par0_new)
obs$par()

Read formula with as.character without splitting

Description

Read formula with as.character without splitting

Usage

as_character_formula(x, ...)

Arguments

Details

Citation: this function was taken from the R package formula.tools: Christopher Brown (2018). formula.tools: Programmatic Utilities for Manipulating Formulas, Expressions, Calls, Assignments and Other R Objects. R package version 1.7.1. https://CRAN.R-project.org/package=formula.tools

Transforms matrix to dgTMatrix

Description

Transforms matrix to dgTMatrix

Usage

as_sparse(x)

Arguments

Value

Sparse matrix of class dgTMatrix

Create block diagonal matrix (safe version)

Description

This version of bdiag checks whether the matrices passed as arguments are NULL. This avoids errors that would arise if using bdiag directly.

Usage

bdiag_check(...)

Arguments

Value

Block diagonal matrix

Check values in vector are contiguous

Description

Check values in vector are contiguous

Usage

check_contiguous(x)

Arguments

Value

Logical: are values contiguous?

Check that only supported smooths are used

Description

Currently, this raises an error message if tensor product smooths are used (te, ti, t2).

Usage

check_smooths(form)

Arguments

Grid of covariates

Description

Grid of covariates

Usage

cov_grid(var, data = NULL, obj = NULL, covs = NULL, formulas, n_grid = 1000)

Arguments

Value

Data frame of covariates, with 'var' defined over a grid, and other covariates fixed to their mean (numeric) or first level (factor).

Grid of covariate values over two covariates

Description

The two-dimensional counterpart of cov_grid, used by HMM$plot_2d(). Both var and var2 are varied over a regular n_grid by n_grid lattice spanning their ranges in the data, and every other covariate is held at one value, as in cov_grid: the one given in covs, otherwise the mean of a numeric covariate or the first level of a factor.

Usage

cov_grid_2d(
  var,
  var2,
  data = NULL,
  obj = NULL,
  covs = NULL,
  formulas,
  n_grid = 40,
  padding = 0
)

Arguments

Details

The lattice is regular because that is what ggplot2::geom_raster() needs, and var varies fastest, so that the result can be reshaped into a matrix column by column.

Value

Data frame with n_grid^2 rows and one column for each covariate

Covariate names that are columns of the data

Description

all.vars() on a model formula returns every symbol in it, including objects a smooth refers to through its xt argument – the mesh of an SPDE field, say – which are not covariates and are not columns of the data. Those have to be dropped before the data are indexed by name.

Usage

cov_names_in_data(var_names, data)

Arguments

Details

A covariate that is genuinely missing from the data is dropped here too, rather than raising an error. mgcv::gam() reports it when the design matrices are built, which is a better place to say so.

Value

The subset of var_names naming columns of data

Density function of von Mises distribution

Description

Density function of von Mises distribution

Usage

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

Arguments

Value

Von Mises density

Density function of wrapped Cauchy distribution

Description

Density function of wrapped Cauchy distribution

Usage

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

Arguments

Value

Wrapped Cauchy density

Find s(, bs = "re") terms in formula

Description

This function is used to identify the variables "x" which are included as s(x, bs = "re") in the formula, in particular to check that they are factors.

Usage

find_re(form)

Arguments

Value

Vector of names of variables for which a random effect term is included in the model.

Shell gam object for building prediction matrices

Description

mgcv::predict.gam() with type = "lpmatrix" uses only the model frame, the terms objects, the factor levels and contrasts, and the smooth objects – all of which mgcv::gam(fit = FALSE) already produces. This wraps that unfitted setup in something predict.gam() accepts, so that a prediction matrix can be built without fitting anything.

Usage

gam_shell(G)

Arguments

Details

The coefficients are set to zero. predict.gam() requires the slot to exist and to have the right length, but for type = "lpmatrix" it returns the design matrix itself and never multiplies by them.

Value

An object of class "gam", only usable for predict(type = "lpmatrix")

Generalized matrix determinant

Description

Generalized determinant = product of non-zero eigenvalues (see e.g., Wood 2017). Used for (log)determinant of penalty matrices, required in log-likelihood function.

Usage

gdeterminant(x, eps = 1e-10, log = TRUE)

Arguments

Value

Generalized determinant of input matrix

hmmTMB colour palette

Description

hmmTMB colour palette

Usage

hmmTMB_cols

Multivarite inverse logit function

Description

Multivarite inverse logit function

Usage

invmlogit(x)

Arguments

Check if number of whole number

Description

Check if number of whole number

Usage

is_whole_number(x, tol = 1e-10)

Arguments

Value

TRUE if it is a whole number within tolerance

logLik function for SDE objects

Description

This function makes it possible to call generic R methods such as AIC and BIC on HMM objects. It is based on the number of degrees of freedom of the *conditional* AIC (rather than marginal AIC), i.e., including degrees of freedom from the smooth/random effect components of the model.

Usage

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

Arguments

Value

Maximum log-likelihood value for the model, with attributes df (degrees of freedom) and nobs (number of observations)

Log of sum of exponentials

Description

Log of sum of exponentials

Usage

logsumexp(x)

Arguments

Make covariance matrix from standard deviations and correlations

Description

Make covariance matrix from standard deviations and correlations

Usage

make_cov(sds, corr)

Arguments

Value

An m by m covariance matrix

Process formulas and store in nested list

Description

Process formulas and store in nested list

Usage

make_formulas(input_forms, var_names, par_names, n_states)

Arguments

Details

Formulas for the observation parameters can be different for the different states, using special functions of the form "state1", "state2", etc. This method processes the list of formulas passed by the user to extract the state-specific formulas. Missing formulas are assumed to be intercept-only ~1.

Value

Nested list of formulas, with three levels: observed variable, parameter of the observation distribution, and state.

Examples

input_forms <- list(step = list(shape = ~ state1(x1) + x2,
                                scale = ~ x1),
                    count = list(lambda = ~ state1(x1) + state2(s(x2, bs = "cs"))))

make_formulas(input_forms = input_forms, 
              var_names = names(input_forms), 
              par_names = lapply(input_forms, names), 
              n_states = 2)
              

Create model matrices

Description

Create model matrices

Usage

make_matrices(formulas, data, new_data = NULL, gam_args = NULL)

Arguments

Details

Note the two indices: ncol_re, log_det_S and gmrf have one entry per penalty, while L, sp_names, sp_gmrf and theta_start have one per smoothing parameter. The two coincide for an ordinary smooth, which has one of each, but an SPDE field has three penalties and two parameters.

Value

A list of

• X_fe Design matrix for fixed effects

• X_re Design matrix for random effects

• S Smoothness matrix

• log_det_S Vector of log-determinants of smoothness matrices, one for each penalty

• ncol_fe Number of columns of X_fe for each parameter

• ncol_re Number of columns of X_re and S for each random effect

• L Matrix with log(lambda) = L * theta, mapping the smoothing parameters of each smooth to the weights of its penalties (mgcv's convention). It is the identity unless a smooth combines several penalties through fewer parameters. See ?mgcv::smooth.construct for details.

• gmrf For each penalty, 1 if its log-determinant depends on the parameters and so has to be computed inside the likelihood, and 0 if the penalty is a fixed matrix scaled by exp(theta)

• sp_names Name of each smoothing parameter

• theta_start Starting value of each log smoothing parameter

Multivariate logit function

Description

Multivariate logit function

Usage

mlogit(x)

Arguments

Multivariate Normal inverse link function

Description

Multivariate Normal inverse link function

Usage

mvnorm_invlink(x)

Arguments

Multivariate Normal link function

Description

Multivariate Normal link function

Usage

mvnorm_link(x)

Arguments

Fill in NAs

Description

Replace NA entries in a vector by the last non-NA value. If the first entry of the vector is NA, it is replaced by the first non-NA value. If the vector passed as input doesn't contain NAs, it is returned as is.

Usage

na_fill(x)

Arguments

Value

Copy of x in which NAs have been replaced by nearest available value.

Get covariance matrix from precision matrix

Description

The covariance matrix is the inverse of the precision matrix. By default, the function solve is used for inversion. If it fails (e.g., singular system), then MASS::ginv is used instead, and returns the Moore-Penrose generalised inverse of the precision matrix.

Usage

prec_to_cov(prec_mat)

Arguments

Value

Precision matrix

Solve for positive root of quadratic ax^2 + bx + c = 0 when it exists

Description

Solve for positive root of quadratic ax^2 + bx + c = 0 when it exists

Usage

quad_pos_solve(a, b, c)

Arguments

Value

real positive root if it exists

Sample from a multivariate normal given its precision matrix

Description

Draws from N(mu, Q^-1) through a sparse Cholesky factorisation of the precision, Q = P' L L' P, as x = mu + P' L^-T z with z standard normal. Nothing is inverted, and the Cholesky factor keeps the sparsity of Q, so the cost is that of the factorisation rather than of a dense inverse followed by a dense n-by-n Cholesky of it. A dense precision is simply a sparse matrix with no zeros, and costs no more this way than the old route did.

Usage

rmvn_prec(n, mu, prec_mat)

Arguments

Details

If prec_mat is not positive definite – a fit that has not converged, or a singular Hessian – there is no Cholesky factor to be had, and this falls back on prec_to_cov, which uses a generalised inverse and warns.

Value

Matrix with one row for each sample and one column for each element of mu

Sample from von Mises distribution

Description

Sample from von Mises distribution

Usage

rvm(n, mu, kappa)

Arguments

Details

Uses basic rejection sampling, based on dvm(), which might be inefficient for large kappa. Could be improved following Best & Fisher (1979), Efficient simulation of the von Mises distribution, JRSSC, 28(2), 152-157.

Value

Vector of n samples from vm(mu, kappa)

Sample from wrapped Cauchy distribution

Description

Sample from wrapped Cauchy distribution

Usage

rwrpcauchy(n, mu, rho)

Arguments

Details

Uses basic rejection sampling, based on dwrpcauchy(), which might be inefficient for large rho.

Value

Vector of n samples from wrpcauchy(mu, rho)

Matern SPDE smooth

Description

A smooth term for a spatially (or temporally) continuous Gaussian random field, usable as s(x, y, bs = "spde", xt = list(mesh = mesh)) in any hmmTMB formula, and so on any transition probability or any parameter of any observation distribution.

Usage

## S3 method for class 'spde.smooth.spec'
smooth.construct(object, data, knots)

## S3 method for class 'spde.smooth'
Predict.matrix(object, data)

Arguments

Value

A smoothCon object of class spde.smooth, with sparse X and sparse penalties.

The mesh in one dimension

A one-dimensional field needs no mesh: s(x, bs = "spde") builds one by default, as k quadratic B-spline basis functions on evenly spaced knots, so

f <- ~ s(x, bs = "spde")          # and s(x, bs = "spde", k = 40) for a finer one

is the whole of it. The knots are spread over the range of the covariate widened by a fifth on each side, which is the one-dimensional version of the outer extension of a two-dimensional mesh: with the knots stopping at the data, the boundary condition inflates the field's variance at the two ends by a third or more.

k is a resolution, not a smoothing parameter, and this is the one place where a field does not behave like the smooths around it. The range is estimated and is what does the smoothing, so k is not a cap on wiggliness that the penalty then pulls back from: it only has to be large enough to resolve the range the data ask for, knots about a third of a range apart as in two dimensions. That is why the default is k = 15 rather than mgcv's 10 – with the widening, 15 knots put the knot spacing at about half the range the smoother starts from, which is a fifth of the spread of the covariate.

Too small a k does not simply oversmooth. Once the knots are further apart than the range the data want, the weights are near enough independent that the likelihood cannot tell one short range from another, and the range collapses towards zero while the standard deviation grows to compensate. The fitted curve can still look right when that happens, because the basis is then doing the smoothing on its own, but the two parameters are meaningless. So read lambda() before reading the curve: a range far below the knot spacing, or a standard deviation orders of magnitude above the scale of the linear predictor, means raise k and refit.

Passing your own fmesher::fm_mesh_1d() overrides all of this, and is worth doing when the covariate is unevenly spread and you want the knots placed accordingly.

Building a mesh in two dimensions

In two dimensions the mesh is yours to build and pass in, because a good one depends on the domain, the data and the range you expect. fmesher::fm_mesh_2d() takes locations, a maximum triangle edge length inside the domain and in the outer extension, and a cutoff below which nearby points are merged:

mesh <- fmesher::fm_mesh_2d(loc = cbind(data$x, data$y),
                            max.edge = c(0.1, 0.3),
                            cutoff = 0.05, offset = c(0.1, 0.3))
f <- ~ s(x, y, bs = "spde", xt = list(mesh = mesh))

A rule of thumb is max.edge no larger than a third of the range you expect, and an offset about equal to the range; without the outer extension the field's variance is inflated at the boundary.

The field is evaluated at every time step, so missing covariates are replaced by the last non-missing value as elsewhere in hmmTMB. If that is too crude for locations, interpolate them before fitting.

Parameterisation

The precision is

Q(\tau, \kappa) = \tau^2(\kappa^4 C + 2\kappa^2 G_1 + G_2)

with C, G_1, G_2 the finite element matrices from fmesher::fm_fem(). alpha = 2 throughout, so the Matern smoothness is \nu = 2 - d/2: \nu = 1 in two dimensions and \nu = 3/2 in one. The range is \rho = \sqrt{8\nu}/\kappa and the marginal variance is \sigma^2 = \Gamma(\nu) / (\Gamma(2)(4\pi)^{d/2}\kappa^{2\nu}\tau^2), which is 1/(\tau^2\kappa^2 4\pi) in two dimensions and 1/(4\tau^2\kappa^3) in one. Both constants depend on d, so the one- and two-dimensional smoothers do not share them.

The two estimated parameters are \log\sigma and \log\rho, not \log\tau and \log\kappa, and appear as sd and range in lambda(). The change is exact – a linear map of the log parameters, absorbed into the constants in front of the three matrices – and worth making: the parameters mean something on their own, and the likelihood's long ridge along constant marginal variance runs diagonally in (\tau, \kappa) but along an axis in (\sigma, \rho), which the optimiser finds much easier.

Both are estimated freely, with no prior. A single realisation identifies them only loosely, and a field whose range approaches the size of the domain is close to improper: the range then runs off and the standard deviation follows it, with the fitted surface barely changing. The fit is still usable if that happens, but the two parameters are not, and a finer mesh will not help – it is the data that are uninformative.

The field is proper for any positive \kappa, so no identifiability constraint is imposed, as in INLA. It is not centred, and its level is only weakly separated from the intercept when the range is large.

Bandwidth

A model containing this smooth is fitted with the banded forward algorithm, because the exact one makes the Hessian with respect to the field weights dense and so defeats the sparsity the SPDE representation exists to provide. See HMM$update_bw() and HMM$check_bw().

References

Lindgren, F., Rue, H. and Lindstrom, J. (2011). An explicit link between Gaussian fields and Gaussian Markov random fields. JRSS-B 73, 423-498.

Miller, D. L., Glennie, R., & Seaton, A. E. (2020). Understanding the stochastic partial differential equation approach to smoothing. Journal of Agricultural, Biological and Environmental Statistics, 25(1), 1-16.

Fischer, J.-O. (2026). Fast and scalable inference in hidden Markov models with Gaussian fields.

Examples

# A one-dimensional field over a covariate, with a default mesh
if(requireNamespace("fmesher", quietly = TRUE)) {
  d <- data.frame(x = runif(100))
  sm <- mgcv::smoothCon(mgcv::s(x, bs = "spde"), data = d)[[1]]
  c(k = ncol(sm$X), sm$theta.names)
}

A length scale for the meshed domain, for starting values

Description

The larger side of the bounding box of the mesh nodes that carry basis weight, which is the region the data occupy rather than the outer extension. A one-dimensional mesh stores its knots in 'loc' and the midpoints of its basis functions in 'mid', and only the latter are one per column of X, so 'mid' is the one to use where there is a choice.

Usage

spde_extent(mesh, X)

Arguments

Value

A positive number

Observation locations for an SPDE smooth

Description

Observation locations for an SPDE smooth

Usage

spde_loc(term, data)

Arguments

Value

Numeric vector (1d) or two-column matrix (2d) of locations

Default mesh for a one-dimensional SPDE smooth

Description

k quadratic B-spline basis functions on evenly spaced knots, over the range of the covariate widened by a fifth on each side. The widening is the one-dimensional counterpart of the outer extension of a two-dimensional mesh: fmesher's default boundary condition is Neumann, under which a field whose knots stop at the data has its variance inflated by a third or more at the two ends. A fifth of the spread of the covariate is about one range at the starting value used in smooth.construct.spde.smooth.spec(), which is the same rule of thumb used for the offset in two dimensions.

Usage

spde_mesh_1d(x, k)

Arguments

Details

A degree-2 mesh has one basis function fewer than it has knots, so k knots would not give k of them; k keeps its mgcv meaning of a basis dimension here, and k + 1 knots are laid down.

The default of 15 is above mgcv's 10, because k is a resolution here and not a cap on wiggliness: with the widening, 15 knots put the knot spacing at about half the range the smoother starts from, in the spirit of max.edge in two dimensions. Ten of them leave the knots further apart than that range, and it collapses towards zero as a result.

Value

An fm_mesh_1d with k basis functions

Strip comments marked with a hash from a character vector

Description

Strip comments marked with a hash from a character vector

Usage

strip_comments(str)

Arguments

Value

character vector with comments removed (and lines with only comments completely removed)

Update a model to a new model by changing one formula

Description

Update a model to a new model by changing one formula

Usage

## S3 method for class 'HMM'
update(object, type, i, j, change, fit = TRUE, silent = FALSE, ...)

Arguments

Examples

# Load data set from MSwM package
data(energy, package = "MSwM")

# Create hidden state and observation models
hid <- MarkovChain$new(data = energy, n_states = 2)
par0 <- list(Price = list(mean = c(3, 6), sd = c(2, 3)))
obs <- Observation$new(data = energy, n_states = 2,
                       dists = list(Price = "norm"),
                       par = par0)

# Create HMM (no covariate effects)
hmm <- HMM$new(hid = hid, obs = obs)
hmm$hid()$formula()
hmm$obs()$formulas()

# Update transition probability formulas (one at a time)
hmm <- update(hmm, type = "hid", i = 1, j = 2, 
              change = ~ . + Oil, fit = FALSE)
hmm <- update(hmm, type = "hid", i = 2, j = 1, 
              change = ~ . + Gas + Coal, fit = FALSE)
hmm$hid()$formula()

# Update observation parameter formulas (one at a time)
hmm <- update(hmm, type = "obs", i = "Price", j = "mean", 
              change = ~ . + EurDol, fit = FALSE)
hmm$obs()$formulas()