HMM {hmmTMB}R Documentation

R6 class for hidden Markov model

Description

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

Methods

Public methods

Method new()

Create new HMM object

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

Observation object, created with Observation$new(). This contains the formulation for the observation model.

hid

MarkovChain object, created with MarkovChain$new(). This contains the formulation for the state process model.

file

Path to specification file for HMM. If this argument is used, then obs and hid are unnecessary.

init

HMM object, used to initialise the parameters for this model. If init is passed, then all parameters that are included in init and in the present model are copied. This may be useful when fitting increasingly complex models: start from a simple model, then pass it as init to create a more complex model, and so on.

fixpar

Named list, with optional elements: 'hid', 'obs', 'delta0', 'lambda_obs', and 'lambda_hid'. Each element is a named vector of parameters in coeff_fe that should either be fixed (if the corresponding element is set to NA) or estimated to a common value (using integers or factor levels).don See examples in the vignettes, and check the TMB documentation to understand the inner workings (argument map of TMB::MakeADFun()).

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
all

Logical. If FALSE, only user-specified fixed parameters are returned, but not parameters that are fixed by definition (e.g., size of binomial distribution).

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
par_list

List with elements for coeff_fe_obs, coeff_fe_hid, coeff_re_obs, coeff_re_hid, log_delta0, log_lambda_hid, and log_lambda_obs

iter

Optional argument to update model parameters based on MCMC iterations (if using rstan). Either the index of the iteration to use, or "mean" if the posterior mean should be used.

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
t

returns parameters at time t, default is t = 1

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
new_priors

is a list of matrices for optionally coeff_fe_obs, coeff_fe_hid, log_lambda_obs log_lambda_hid each matrix has two rows (first row = mean, second row = sd) specifying parameters for Normal priors

Method priors()

Extract stored priors

Usage
HMM$priors()

Method iters()

Iterations from stan MCMC fit

Usage
HMM$iters(type = "response")
Arguments
type

Either "response" for parameters on the response (natural) scale, or "raw" for parameters on the linear predictor scale.

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
silent

Logical. If TRUE, all tracing outputs are hidden (default).

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

Arguments passed to tmbstan

silent

Logical. If FALSE, all tracing outputs are shown (default).

Method fit()

Model fitting

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

After the model has been fitted, the output of optimx() 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
silent

Logical. If FALSE, all tracing outputs are shown (default).

...

Other arguments to optimx which is used to optimise likelihood, see ?optimx

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
nsamp

Number of samples to produce

full

If TRUE and model fit by fit_stan then parameter estimates are sampled from the posterior samples before simulating each sequence

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
n_post

Number of posterior samples

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
n_post

Number of posterior samples

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
fn

Function which takes a vector of linear predictors as input and produces either a scalar or vector output

n_post

Number of posterior simulations

comp

Either "obs" for observation model linear predictor, or "hid" for hidden model linear predictor

...

Arguments passed to fn

level

Confidence interval level if required (e.g., 0.95 for 95 confidence intervals). Default is 0, i.e., confidence intervals are not returned.

return_post

Logical indicating whether to return the posterior samples. If FALSE (default), only mean estimates and confidence intervals are returned

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

Which estimates to predict? Options include transition probability matrices "tpm", stationary distributions "delta", or observation distribution parameters "obspar"

t

Time points to predict at

newdata

New dataframe to use for prediction

n_post

If greater than zero then n_post posterior samples are produced, and used to create confidence intervals.

level

Level of the confidence intervals, e.g. CI = 0.95 will produce 95% confidence intervals (default)

return_post

Logical. If TRUE, a list of posterior samples is returned.

...

Other arguments to the respective functions for hid$tpm, hid$delta, obs$par

Returns

Named array of predictions and confidence intervals, if requested

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
level

Level of confidence intervals. Defaults to 0.95, i.e., 95% confidence intervals.

Returns

List of matrices with three columns: mle (maximum likelihood estimate), lcl (lower confidence limit), and ucl (upper confidence limit). 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
n

Number of time steps to simulate

data

Optional data frame including covariates

silent

if TRUE then no messages are printed

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
check_fn

Goodness-of-fit function which accepts "data" as input and returns a statistic (either a vector or a single number) to be compared between observed data and simulations.

nsims

Number of simulations to perform

full

If model fitted with 'fit_stan', then full = TRUE will sample from posterior for each simulation

silent

Logical. If FALSE, simulation progress is shown. (Default: TRUE)

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

Name of the variable to plot.

var2

Optional name of a second variable, for 2-d plot.

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

Name of data variable

Returns

Plot of distribution 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
what

Name of model component to plot: should be one of "tpm" (transition probabilities), "delta" (stationary state probabilities), or "obspar" (state-dependent observation parameters)

var

Name of covariate to plot on x-axis

covs

Optional named list for values of covariates (other than 'var') that should be used in the plot (or dataframe with single row). If this is not specified, the mean value is used for numeric variables, and the first level for factor variables.

i

If plotting tpm then rows of tpm; if plotting delta then indices of states to plot; if plotting obspar then full names of parameters to plot (e.g., obsvar.mean)

j

If plotting tpm then columnss of tpm to plot; if plotting delta then this is ignored,; if plotting obspar then indices of states to plot

n_grid

Number of points in grid over x-axis (default: 50)

n_post

Number of posterior simulations to use when computing confidence intervals; default: 1000. See predict function for more detail.

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
deep

Whether to make a deep clone.

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)
[Package hmmTMB version 1.0.2 Index]