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 (2022) arXiv:2211.14139.

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:

• Richard Glennie [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()

Method 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

Method name()

Return name of Dist object

Usage

Dist$name()

Method pdf()

Return pdf of Dist object

Usage

Dist$pdf()

Method cdf()

Return cdf of Dist object

Usage

Dist$cdf()

Method rng()

Return random generator function of Dist object

Usage

Dist$rng()

Method link()

Return link function of Dist object

Usage

Dist$link()

Method invlink()

Return inverse link function of Dist object

Usage

Dist$invlink()

Method npar()

Return number of parameters of Dist object

Usage

Dist$npar()

Method parnames()

Return names of parameters

Usage

Dist$parnames()

Method parapprox()

Return function that approximates parameters

Usage

Dist$parapprox()

Method fixed()

Return which parameters are fixed

Usage

Dist$fixed()

Method code()

Return code of Dist object

Usage

Dist$code()

Method name_long()

Human-readable name of Dist object

Usage

Dist$name_long()

Method set_npar()

Set number of parameters this distribution has

Usage

Dist$set_npar(new_npar)

Arguments

Method set_parnames()

Set parameter names

Usage

Dist$set_parnames(new_parnames)

Arguments

Method set_code()

Set distribution code

Usage

Dist$set_code(new_code)

Arguments

Method 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

Method 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

Method par_alt()

Alternative parameter formatting

Usage

Dist$par_alt(par)

Arguments

Returns

Formatted parameters

Method 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

Method 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

Method clone()

The objects of this class are cloneable with this method.

Usage

Dist$clone(deep = FALSE)

Arguments

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$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$AIC_marginal()

• HMM$AIC_conditional()

• HMM$print_obspar()

• HMM$print_tpm()

• HMM$formulation()

• HMM$print()

• HMM$clone()

Method new()

Create new HMM object

Usage

HMM$new(obs = NULL, hid = NULL, file = NULL, init = NULL, fixpar = 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)

Method obs()

Observation object for this model

Usage

HMM$obs()

Method hid()

MarkovChain object for this model

Usage

HMM$hid()

Method out()

Output of optimiser after model fitting

Usage

HMM$out()

Method 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()

Method 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()

Method tmb_rep()

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

Usage

HMM$tmb_rep()

Method states()

Vector of estimated states, after viterbi has been run

Usage

HMM$states()

Method coeff_fe()

Coefficients for fixed effect parameters

Usage

HMM$coeff_fe()

Method coeff_re()

Coefficients for random effect parameters

Usage

HMM$coeff_re()

Method 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()

Method fixpar()

Fixed parameters

Usage

HMM$fixpar(all = FALSE)

Arguments

Method coeff_array()

Array of working parameters

Usage

HMM$coeff_array()

Method lambda()

Smoothness parameters

Usage

HMM$lambda()

Method update_par()

Update parameters stored inside model object

Usage

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

Arguments

Method 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.

Method 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

Method set_priors()

Set priors for coefficients

Usage

HMM$set_priors(new_priors = NULL)

Arguments

Method priors()

Extract stored priors

Usage

HMM$priors()

Method iters()

Iterations from stan MCMC fit

Usage

HMM$iters(type = "response")

Arguments

Returns

see output of as.matrix in stan

Method out_stan()

fitted stan object from MCMC fit

Usage

HMM$out_stan()

Returns

the stanfit object

Method llk()

Log-likelihood at current parameters

Usage

HMM$llk()

Returns

Log-likelihood

Method 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)

Method 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

Method 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

Method 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

Method 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)

Method 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()

Method 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

Method 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)

Method viterbi()

Viterbi algorithm

Usage

HMM$viterbi()

Returns

Most likely state sequence

Method 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

Method 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

Method 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

Method 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

Method 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)

Method 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)

Method 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.

Method 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

Method 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

Method 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

Method 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

Method 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

Method 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

Method 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

Method print_obspar()

Print observation parameters at t = 1

Usage

HMM$print_obspar()

Method print_tpm()

Print observation parameters at t = 1

Usage

HMM$print_tpm()

Method formulation()

Print model formulation

Usage

HMM$formulation()

Method print()

Print HMM object

Usage

HMM$print()

Method 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()

Method 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)

Method formula()

Formula of MarkovChain model

Usage

MarkovChain$formula()

Method formulas()

List of formulas for MarkovChain model

Usage

MarkovChain$formulas()

Method tpm()

Get transition probability matrices

Usage

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

Arguments

Returns

Array with one slice for each transition probability matrix

Method ref()

Indices of reference elements in transition probability matrix

Usage

MarkovChain$ref()

Method ref_mat()

Matrix of reference elements in transition probability matrix

Usage

MarkovChain$ref_mat()

Method ref_delta0()

Indices of reference elements in initial distribution

Usage

MarkovChain$ref_delta0()

Method coeff_fe()

Current parameter estimates (fixed effects)

Usage

MarkovChain$coeff_fe()

Method 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.

Method 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)

Method stationary()

Use stationary distribution as initial distribution?

Usage

MarkovChain$stationary()

Method fixpar()

Fixed parameters

Usage

MarkovChain$fixpar(all = FALSE)

Arguments

Method coeff_re()

Current parameter estimates (random effects)

Usage

MarkovChain$coeff_re()

Method X_fe()

Fixed effect design matrix

Usage

MarkovChain$X_fe()

Method X_re()

Random effect design matrix

Usage

MarkovChain$X_re()

Method lambda()

Smoothness parameters

Usage

MarkovChain$lambda()

Method 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.

Usage

MarkovChain$sd_re()

Method nstates()

Number of states

Usage

MarkovChain$nstates()

Method terms()

Terms of model formulas

Usage

MarkovChain$terms()

Method unique_ID()

Number of time series

Usage

MarkovChain$unique_ID()

Method initial_state()

Initial state (see constructor argument)

Usage

MarkovChain$initial_state()

Method empty()

Empty model? (for simulation only)

Usage

MarkovChain$empty()

Method gam_args()

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

Usage

MarkovChain$gam_args()

Method update_tpm()

Update transition probability matrix

Usage

MarkovChain$update_tpm(tpm)

Arguments

Method update_coeff_fe()

Update coefficients for fixed effect parameters

Usage

MarkovChain$update_coeff_fe(coeff_fe)

Arguments

Method update_coeff_re()

Update coefficients for random effect parameters

Usage

MarkovChain$update_coeff_re(coeff_re)

Arguments

Method update_X_fe()

Update design matrix for fixed effects

Usage

MarkovChain$update_X_fe(X_fe)

Arguments

Method update_X_re()

Update design matrix for random effects

Usage

MarkovChain$update_X_re(X_re)

Arguments

Method update_delta0()

Update initial distribution

Usage

MarkovChain$update_delta0(delta0)

Arguments

Method update_lambda()

Update smoothness parameters

Usage

MarkovChain$update_lambda(lambda)

Arguments

Method update_fixpar()

Update information about fixed parameters

Usage

MarkovChain$update_fixpar(fixpar)

Arguments

Method 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

Method 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.

Method 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

Method 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

Method linpred()

Linear predictor for transition probabilities

Usage

MarkovChain$linpred()

Method simulate()

Simulate from Markov chain

Usage

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

Arguments

Returns

Sequence of states of simulated chain

Method formulation()

Print model formulation

Usage

MarkovChain$formulation()

Method print()

Print MarkovChain object

Usage

MarkovChain$print()

Method 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()

Method 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)

Method data()

Data frame

Usage

Observation$data()

Method dists()

List of distributions

Usage

Observation$dists()

Method nstates()

Number of states

Usage

Observation$nstates()

Method 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))

Method 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

Method inipar()

Return initial parameter values supplied

Usage

Observation$inipar()

Method coeff_fe()

Fixed effect parameters on working scale

Usage

Observation$coeff_fe()

Method coeff_re()

Random effect parameters

Usage

Observation$coeff_re()

Method X_fe()

Fixed effect design matrix

Usage

Observation$X_fe()

Method X_re()

Random effect design matrix

Usage

Observation$X_re()

Method lambda()

Smoothness parameters

Usage

Observation$lambda()

Method 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.

Usage

Observation$sd_re()

Method formulas()

List of model formulas for observation model

Usage

Observation$formulas(raw = FALSE)

Arguments

Method terms()

Terms of model formulas

Usage

Observation$terms()

Method obs_var()

Data frame of response variables

Usage

Observation$obs_var(expand = FALSE)

Arguments

Returns

Data frame of observation variables

Method known_states()

Vector of known states

Usage

Observation$known_states(mat = TRUE)

Arguments

Method fixpar()

Fixed parameters

Usage

Observation$fixpar(all = FALSE)

Arguments

Method empty()

Empty model? (for simulation only)

Usage

Observation$empty()

Method gam_args()

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

Usage

Observation$gam_args()

Method 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

Method update_coeff_fe()

Update coefficients for fixed effect parameters

Usage

Observation$update_coeff_fe(coeff_fe)

Arguments

Method update_coeff_re()

Update random effect parameters

Usage

Observation$update_coeff_re(coeff_re)

Arguments

Method update_X_fe()

Update fixed effect design matrix

Usage

Observation$update_X_fe(X_fe)

Arguments

Method update_X_re()

Update random effect design matrix

Usage

Observation$update_X_re(X_re)

Arguments

Method update_lambda()

Update smoothness parameters

Usage

Observation$update_lambda(lambda)

Arguments

Method update_data()

Update data

Usage

Observation$update_data(data)

Arguments

Method update_fixpar()

Update information about fixed parameters

Usage

Observation$update_fixpar(fixpar)

Arguments

Method 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

Method 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.

Method 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

Method 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

Method linpred()

Compute linear predictor

Usage

Observation$linpred()

Method 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.

Method 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.

Method 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()

Method 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

Method formulation()

Print model formulation

Usage

Observation$formulation()

Method print()

Print Observation object Check constructor arguments

Usage

Observation$print()

Method 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?

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).

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.

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

Format

An object of class character of length 6.

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

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

• 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

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 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)

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()