HMM: R6 class for hidden Markov model
In TheoMichelot/hmmTMB: Fit Hidden Markov Models using Template Model Builder
HMM 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
obsObservation object, created with Observation$new().
This contains the formulation for the observation model.
hidMarkovChain object, created with MarkovChain$new().
This contains the formulation for the state process model.
filePath to specification file for HMM. If this argument is
used, then obs and hid are unnecessary.
initHMM 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.
fixparNamed 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
fnObjective function
grGradient function of fn
parVector 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
fnObjective function
grGradient function of fn
parVector 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
allLogical. 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_listList with elements for coeff_fe_obs,
coeff_fe_hid, coeff_re_obs, coeff_re_hid, log_delta0,
log_lambda_hid, and log_lambda_obs
iterOptional 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
treturns parameters at time t, default is t = 1
Returns
A list with elements:
obsparParameters of observation model
tpmTransition probability matrix of hidden state model
Method set_priors()
Set priors for coefficients
Usage
HMM$set_priors(new_priors = NULL)
Arguments
new_priorsis 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
typeEither "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
silentLogical. 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
silentLogical. 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
silentLogical. 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
nsampNumber of samples to produce
fullIf 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_postNumber 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_postNumber 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
fnFunction which takes a vector of linear predictors as input
and produces either a scalar or vector output
n_postNumber of posterior simulations
compEither "obs" for observation model linear predictor, or
"hid" for hidden model linear predictor
...Arguments passed to fn
levelConfidence interval level if required (e.g., 0.95 for 95
confidence intervals). Default is 0, i.e., confidence intervals are not
returned.
return_postLogical 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
whatWhich estimates to predict? Options include
transition probability matrices "tpm",
stationary distributions "delta", or
observation distribution parameters "obspar"
tTime points to predict at
newdataNew dataframe to use for prediction
n_postIf greater than zero then n_post posterior
samples are produced, and used to create confidence intervals.
levelLevel of the confidence intervals, e.g. CI = 0.95
will produce 95% confidence intervals (default)
return_postLogical. If TRUE, a list of posterior samples
is returned.
as_listLogical. 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
levelLevel 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
nNumber of time steps to simulate
dataOptional data frame including covariates
silentif 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_fnGoodness-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.
nsimsNumber of simulations to perform
fullIf model fitted with 'fit_stan', then full = TRUE
will sample from posterior for each simulation
silentLogical. 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
varName of the variable to plot.
var2Optional name of a second variable, for 2-d plot.
lineLogical. 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
varName 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
whatName of model component to plot: should be one of "tpm"
(transition probabilities), "delta" (stationary state probabilities),
or "obspar" (state-dependent observation parameters)
varName of covariate to plot on x-axis
covsOptional 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.
iIf 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)
jIf 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_gridNumber of points in grid over x-axis (default: 50)
n_postNumber 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
deepWhether 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)
TheoMichelot/hmmTMB documentation built on Dec. 13, 2024, 11:52 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| HMM | 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
obsObservation object, created with
Observation$new(). This contains the formulation for the observation model.hidMarkovChain object, created with
MarkovChain$new(). This contains the formulation for the state process model.filePath to specification file for HMM. If this argument is used, then
obsandhidare unnecessary.initHMM object, used to initialise the parameters for this model. If
initis 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.fixparNamed 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
mapofTMB::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
fnObjective function
grGradient function of fn
parVector 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
fnObjective function
grGradient function of fn
parVector 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
allLogical. 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_listList with elements for coeff_fe_obs, coeff_fe_hid, coeff_re_obs, coeff_re_hid, log_delta0, log_lambda_hid, and log_lambda_obs
iterOptional 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
treturns parameters at time t, default is t = 1
Returns
A list with elements:
obsparParameters of observation model
tpmTransition probability matrix of hidden state model
Method set_priors()
Set priors for coefficients
Usage
HMM$set_priors(new_priors = NULL)
Arguments
new_priorsis 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
typeEither "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
silentLogical. 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
silentLogical. 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
silentLogical. 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 aseval.maxanditer.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
nsampNumber of samples to produce
fullIf 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_postNumber 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_postNumber 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
fnFunction which takes a vector of linear predictors as input and produces either a scalar or vector output
n_postNumber of posterior simulations
compEither "obs" for observation model linear predictor, or "hid" for hidden model linear predictor
...Arguments passed to fn
levelConfidence interval level if required (e.g., 0.95 for 95 confidence intervals). Default is 0, i.e., confidence intervals are not returned.
return_postLogical 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
whatWhich estimates to predict? Options include transition probability matrices "tpm", stationary distributions "delta", or observation distribution parameters "obspar"
tTime points to predict at
newdataNew dataframe to use for prediction
n_postIf greater than zero then n_post posterior samples are produced, and used to create confidence intervals.
levelLevel of the confidence intervals, e.g. CI = 0.95 will produce 95% confidence intervals (default)
return_postLogical. If TRUE, a list of posterior samples is returned.
as_listLogical. 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 ifn_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
levelLevel 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
nNumber of time steps to simulate
dataOptional data frame including covariates
silentif 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_fnGoodness-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.
nsimsNumber of simulations to perform
fullIf model fitted with 'fit_stan', then full = TRUE will sample from posterior for each simulation
silentLogical. 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
varName of the variable to plot.
var2Optional name of a second variable, for 2-d plot.
lineLogical. 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
varName 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
whatName of model component to plot: should be one of "tpm" (transition probabilities), "delta" (stationary state probabilities), or "obspar" (state-dependent observation parameters)
varName of covariate to plot on x-axis
covsOptional 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.
iIf 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)
jIf 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_gridNumber of points in grid over x-axis (default: 50)
n_postNumber of posterior simulations to use when computing confidence intervals; default: 1000. See
predictfunction 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
deepWhether 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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.