# R6 class for hidden Markov model Encapsulates the observation and hidden state models for a hidden Markov model. ## Examples ```r ## ------------------------------------------------ ## 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) ``` ## Methods ### Public methods * [list("HMM$new()")](#method-HMM-new) * [list("HMM$obs()")](#method-HMM-obs) * [list("HMM$hid()")](#method-HMM-hid) * [list("HMM$out()")](#method-HMM-out) * [list("HMM$tmb_obj()")](#method-HMM-tmb_obj) * [list("HMM$tmb_obj_joint()")](#method-HMM-tmb_obj_joint) * [list("HMM$tmb_rep()")](#method-HMM-tmb_rep) * [list("HMM$states()")](#method-HMM-states) * [list("HMM$coeff_fe()")](#method-HMM-coeff_fe) * [list("HMM$coeff_re()")](#method-HMM-coeff_re) * [list("HMM$coeff_list()")](#method-HMM-coeff_list) * [list("HMM$fixpar()")](#method-HMM-fixpar) * [list("HMM$coeff_array()")](#method-HMM-coeff_array) * [list("HMM$lambda()")](#method-HMM-lambda) * [list("HMM$update_par()")](#method-HMM-update_par) * [list("HMM$sd_re()")](#method-HMM-sd_re) * [list("HMM$par()")](#method-HMM-par) * [list("HMM$set_priors()")](#method-HMM-set_priors) * [list("HMM$priors()")](#method-HMM-priors) * [list("HMM$iters()")](#method-HMM-iters) * [list("HMM$out_stan()")](#method-HMM-out_stan) * [list("HMM$llk()")](#method-HMM-llk) * [list("HMM$edf()")](#method-HMM-edf) * [list("HMM$suggest_initial()")](#method-HMM-suggest_initial) * [list("HMM$setup()")](#method-HMM-setup) * [list("HMM$fit_stan()")](#method-HMM-fit_stan) * [list("HMM$fit()")](#method-HMM-fit) * [list("HMM$mle()")](#method-HMM-mle) * [list("HMM$forward_backward()")](#method-HMM-forward_backward) * [list("HMM$pseudores()")](#method-HMM-pseudores) * [list("HMM$viterbi()")](#method-HMM-viterbi) * [list("HMM$sample_states()")](#method-HMM-sample_states) * [list("HMM$state_probs()")](#method-HMM-state_probs) * [list("HMM$post_coeff()")](#method-HMM-post_coeff) * [list("HMM$post_linpred()")](#method-HMM-post_linpred) * [list("HMM$post_fn()")](#method-HMM-post_fn) * [list("HMM$predict()")](#method-HMM-predict) * [list("HMM$confint()")](#method-HMM-confint) * [list("HMM$simulate()")](#method-HMM-simulate) * [list("HMM$check()")](#method-HMM-check) * [list("HMM$plot_ts()")](#method-HMM-plot_ts) * [list("HMM$plot_dist()")](#method-HMM-plot_dist) * [list("HMM$plot()")](#method-HMM-plot) * [list("HMM$AIC_marginal()")](#method-HMM-AIC_marginal) * [list("HMM$AIC_conditional()")](#method-HMM-AIC_conditional) * [list("HMM$print_obspar()")](#method-HMM-print_obspar) * [list("HMM$print_tpm()")](#method-HMM-print_tpm) * [list("HMM$formulation()")](#method-HMM-formulation) * [list("HMM$print()")](#method-HMM-print) * [list("HMM$clone()")](#method-HMM-clone) ### 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 named list of matrices with optional elements coeff_fe_obs, coeff_fe_hid, log_lambda_obs, andlog_lambda_hid. Each matrix has two columns (first col = mean, second col = 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 `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 - **`silent`**: Logical. If FALSE, all tracing outputs are shown (default). - **`...`**: Other arguments to nlminb which is used to optimise the likelihood. This currently only supports the additional argument `control`, which is a list of control parameters such as `eval.max` and `iter.max` (see `?nlminb`) #### 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, as_list = TRUE ) ``` #### 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. - **`as_list`**: Logical. If confidence intervals are required for the transition probabilities or observation parameters, this argument determines whether the MLE, lower confidence limit and upper confidence limit are returned as separate elements in a list (if TRUE; default), or whether they are combined into a single array (if FALSE). Ignored if `what = "delta"` or if `n_post = 0`. - **`...`**: Other arguments to the respective functions for hid$tpm, hid$delta, obs$par #### 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 - **`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), 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 - **`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, line = TRUE) ``` #### Arguments - **`var`**: Name of the variable to plot. - **`var2`**: Optional name of a second variable, for 2-d plot. - **`line`**: Logical. If TRUE (default), lines are drawn between successive data points. Can be set to FALSE if another geom is needed (e.g., geom_point). #### 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 - **`var`**: Name of data variable. If NULL, a list of plots are returned (one for each observation variable) #### 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 - **`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.

HMM function