Package 'HMMscpg'

Description

Tools for modelling incomplete geological point process records (such as volcanic eruptions or earthquakes) using hidden Markov models (HMMs) whose state-dependent observation distributions are gamma, shifted compound Poisson-gamma (SCPG), or time-dependent SCPG (TSCPG) renewal processes. The package implements four families of models: standard HMMs, inhomogeneous HMMs (IHMMs), hidden semi-Markov models (HSMMs), and inhomogeneous hidden semi-Markov models (IHSMMs), following Shahzadi et al. (2023, 2026).

Details

Consider a point process with event times

$H_0 < H_1 < \dots < H_t < \dots$

at jump points $0, 1, 2, \dots, t, \dots$. A renewal process is defined by interevent times $X_t = H_t - H_{t-1}$, $t = 1, 2, \dots$, which are independently and identically distributed, with $H_0 = 0$. Geological records of random events can be modelled using these processes (Wang et al., 2022). Missing events are common in geophysical point process records, and incompleteness varies over time. Given the unknown nature of missing events, Shahzadi et al.(2023, 2026) proposed several types of HMMs in which the unknown number of missing events between each pair of consecutively observed events is characterized by the distribution of the observed process, the underlying Markov chain, or both.

The package employs four types of HMMs, distinguished by the underlying discrete-time Markov chain (MC), which is described by an initial distribution, a transition probability matrix (TPM), and implicit or explicit state durations. The four types HMMs based on underlying hidden Markov chain are:

HMM

has a homogeneous Markov chain (MC) with a time-independent transition probability matrix.

IHMM

has an inhomogeneous Markov chain (IMC) with a time-dependent transition probability matrix.

HSMM

has a semi-Markov chain (SMC) with a zero-diagonal transition probability matrix and explicit state durations, following the Markov property at transition times only.

IHSMM

has an inhomogeneous semi-Markov chain (ISMC) with a zero-diagonal transition probability matrix and time-dependent explicit state durations, following the Markov property at transition times only.

For an observation sequence $\{X_t\}$ and hidden state sequence $\{S_t\}$, $t = 1, \dots, T$, with $m$ states, states are defined as:

State 1

denotes completeness: no missing events between each pair of consecutively observed events.

States $2, \dots, m$

represent incompleteness: the number of missing events between each consecutive observed pair.

Each HMM type is further categorized by its observed process as follows:

State 1

In State 1 the observed process follows a gamma renewal process;

In States $2, \dots, m$

the observed process may follow

(a)

a sum of independent gamma renewal processes when $N$ is fixed, or

(b)

an SCPG renewal process when $N$ is shifted-Poisson distributed, or

(c)

a TSCPG renewal process when $N$ is time-dependent.

The table below lists all types of HMMs considered in this package, according to the underlying Markov chain (MC) and the renewal process (RP) as the observed process (OP). In all models, State 1 represents no missing event between each consecutively observed pair; States $2, \dots, m$ represent incompleteness, with the number of missing events $N$ between each consecutive observed pair as given in the last column.

Here, $N$ denotes the number of missing events between each consecutive observed pair in a point process record.

Author(s)

Maintainer: Amina Shahzadi [email protected]

Authors:

• Amina Shahzadi [email protected]

• Ting Wang [email protected]

References

Shahzadi A, Wang T, Parry M, Bebbington M (2023). “Inhomogeneous hidden semi-Markov models for incompletely observed point processes.” Annals of Institute Statistical Mathematics, 75, 253-280. doi:10.1007/s10463-022-00843-5.

Shahzadi A, Wang T, Parry M, Bebbington M (2026). “Modelling time-inhomogeneous incomplete records of point processes using variants of hidden Markov models.” Advances in Data Analysis and Classification, 20, 663-691. doi:10.1007/s11634-025-00632-x.

Wang T, Bebbington M, Cronin S, Carman J (2022). “Forecasting Eruptions at Poorly Known Volcanoes Using Analogs and Multivariate Renewal Processes.” Geophysical Research Letters, 49(12), e2021GL096715. doi:10.1029/2021GL096715.

See Also

dscpg, pscpg, fit.scpg, WBsum

Description

Evaluates the density of the shifted compound Poisson-gamma (SCPG) distribution, of its time-dependent version (TSCPG), or of a gamma distribution with a fixed number of missing events.

Usage

dscpg(
  x,
  rate,
  shape,
  scale,
  rate.fix = FALSE,
  shape.fix = FALSE,
  time.scpg = FALSE,
  log = FALSE,
  WBsum.method = "direct.evaluation",
  WBsum.tol = 1e-12,
  max.direct.terms = 2000L
)

Arguments

Details

The model. In a complete record the interevent times $Y_k$ are independent gamma variables with shape $\alpha$ and scale $\beta$. If $N$ events are missing between two observed events, the observed interevent time is the sum of $N + 1$ of them,

$X_t = Y_1 + \dots + Y_{N+1}.$

Three cases are available.

1. Fixed number of missing events (rate.fix = TRUE): $N = n$, so $X_t$ is gamma with shape $(n+1)\alpha$ and scale $\beta$; $n = 0$ means that no event is missing.

2. SCPG (the default): at least one event is missing and $N$ is shifted Poisson, $\Pr(N = n) = e^{-\lambda}\lambda^{n-1}/(n-1)!$ for $n = 1, 2, \dots$, so on average $1 + \lambda$ events are missing. The density is

   $f(x_t) = \exp\left(-\lambda - \frac{x_t}{\beta}\right) \frac{x_t^{2\alpha-1}}{\beta^{2\alpha}}\, W(z_t; \alpha), \qquad z_t = \lambda\left(\frac{x_t}{\beta}\right)^{\alpha},$

   where $W(z; \alpha) = \sum_{k \ge 0} z^k/[k!\,\Gamma\{\alpha(k+2)\}]$ is the series computed by WBsum. The SCPG is the continuous part of the compound Poisson-gamma distribution of the Tweedie family (Tweedie, 1984; Smyth and Jorgensen, 2002), shifted so that at least one event is missing.

3. TSCPG (time.scpg = TRUE): the Poisson mean changes along the record as a logistic function of the last occurrence time $h_{t-1} = x_1 + \dots + x_{t-1}$ ($h_0 = 0$),

   $\lambda(t) = K + \frac{M}{1 + \exp(R\,h_{t-1} - P)}, \qquad K, M, R > 0.$

   It decreases from $K + M/(1 + e^{-P})$ at the start of the record (about $K + M$ when $P$ is large) towards $K$. It is halfway, $\lambda(t) = K + M/2$, at the inflection time $h_{t-1} = P/R$, and $R$ sets the speed of the change: 80 per cent of it takes place within $\log(9)/R \approx 2.2/R$ of $P/R$. The density of $x_t$ is the SCPG density with $\lambda(t)$ in place of $\lambda$, given the preceding observations, so the order of x matters.

Choosing WBsum.method. The cost of $W(z_t; \alpha)$ depends on $n^* = (z_t/\alpha^\alpha)^{1/(1+\alpha)}$, which is about the most probable number of missing events in the interval; the three methods are described in WBsum. The default "direct.evaluation" is accurate for every $x_t$ and fast for the values met in practice: $n^*$ is at most a few hundred in typical records, and 6.5 at the estimates of the fit.scpg example. Use "direct.asymptotic" when dscpg stops with the error "The direct Wright-series evaluation is refused for arguments this large", which happens only when $n^* > 10^{12}$, or to save time when $n^*$, computed from the largest $x_t$, exceeds about $10^8$. Such values arise when $x_t$ is very long compared with $\beta$, or when $\lambda$ is very large. They also arise during fit.scpg, whose optimizer may try extreme trial values of the shape and scale even though $n^*$ is small at the estimates. A trial value with any $n^* > 10^{12}$ is refused at once and treated as infeasible, but one with many values of $n^*$ between about $10^8$ and $10^{12}$ takes up to seconds per observation. If a fit is slow, use "direct.asymptotic": it gives the same estimates. "asymptotic.expansion" alone is useful only for checking, since it fails for small $z_t$.

The script verify_dscpg.R in the examples directory of the package reproduces the numbers stated here; run it in the console with

f <- system.file("examples", "verify_dscpg.R",
                 package = "HMMscpg")
source(f, echo = TRUE)

Value

A numeric vector of densities, or log densities when log = TRUE, with the same length and names as x.

References

Tweedie MCK (1984). “An index which distinguishes between some important exponential families.” In Statistics: Applications and New Directions. Proceedings of the Indian Statistical Institute Golden Jubilee International Conference, 579–604.

Smyth GK, Jørgensen B (2002). “Fitting Tweedie's compound Poisson model to insurance claims data: dispersion modelling.” ASTIN Bulletin, 32(1), 143-157.

Smyth GK, Jorgensen B (2002). “Fitting Tweedie's compound Poisson model to insurance claims data: Dispersion modelling.” ASTIN Bulletin, 32(1), 143-157. doi:10.2143/AST.32.1.1020.

Wright EM (1935). “The asymptotic expansion of the generalized hypergeometric function.” Journal of the London Mathematical Society, 10(4), 286–293.

Shahzadi A, Wang T, Parry M, Bebbington M (2026). “Modelling time-inhomogeneous incomplete records of point processes using variants of hidden Markov models.” Advances in Data Analysis and Classification, 20, 663-691. doi:10.1007/s11634-025-00632-x.

Examples

# ===============================
# Example of an SCPG distribution
# ===============================
rate <- 3
shape <- 2
scale <- 0.5
rate.fix <- FALSE
shape.fix <- FALSE

# Simulating from an SCPG distribution

x_scpg <- rscpg(nsim=200, rate, shape, scale, seed=900)
print(x_scpg)
# PDF of an SCPG distribution

pdf_scpg <- dscpg(x_scpg, rate, shape, scale, time.scpg=FALSE,
  WBsum.method="direct.evaluation")
hist(x_scpg, probability = TRUE)
# sorted, so that the density is drawn as a curve
o <- order(x_scpg)
lines(x_scpg[o], pdf_scpg[o], col="blue")
# ===============================
# Example of a TSCPG distribution
# ===============================

# rate[1]=K, rate[2]=M, rate[3]=R, rate[4]=P: lambda(t) falls from about K+M
# to K,
# with inflection time P/R = 300
t_rate <- c(1, 2, 0.02, 6)
shape <- 2
scale <- 0.5
rate.fix <- FALSE
shape.fix <- FALSE

# Simulating from a TSCPG distribution

x_tscpg <- rscpg(nsim=200, t_rate, shape, scale, time.scpg=TRUE, seed=700)
print(x_tscpg)
# PDF of a TSCPG distribution

pdf_tscpg <- dscpg(x_tscpg, t_rate, shape, scale, time.scpg=TRUE,
  WBsum.method="direct.evaluation")
# Each x_t has its own density, at lambda(t) of its occurrence time, so the
# densities are shown as points rather than as one curve
hist(x_tscpg, probability = TRUE)
points(x_tscpg, pdf_tscpg, col="blue", pch=20)

Description

Estimates the parameters of an HMM created by hmm.scpg, of class hmm.gamma, hmm.scpg or hmm.tscpg, by maximum likelihood.

Usage

fit.hmm.scpg(
  hmmobj,
  method = "BFGS",
  maxit = 1000,
  reltol = sqrt(.Machine$double.eps),
  negative = TRUE
)

Arguments

Details

Likelihood. The likelihood sums over all hidden state sequences,

$L = \sum_{s_1, \dots, s_T} \pi_{s_1} f_{s_1}(x_1) \prod_{t=2}^{T} a_{s_{t-1}s_t} f_{s_t}(x_t),$

and is computed by the forward algorithm (Baum et al., 1970),

$\alpha_1(j) = \pi_j f_j(x_1), \qquad \alpha_t(j) = \{\sum_i \alpha_{t-1}(i)a_{ij}\} f_j(x_t), \qquad L = \sum_j \alpha_T(j),$

in $O(Tm^2)$ operations. The forward variables are rescaled to sum to one at each step, on the log scale, and $\log L$ is the sum of the logarithms of the scale factors (Zucchini et al., 2016).

Estimation. $\log L$ is maximized with optim rather than by the EM algorithm (Dempster et al., 1977), since the SCPG densities give no closed-form M-step. The parameters are mapped to unconstrained working parameters: the initial distribution by logits relative to State $m$, each row of the transition matrix by logits relative to its diagonal element $a_{ii}$, and the observation parameters as in fit.scpg, with the same bounds on the shape, the Poisson means and, for the TSCPG, $R_j h_{T-1}$, where $h_{T-1} = x_1 + \dots + x_{T-1}$. With stationary = TRUE the initial distribution is computed from the transition matrix and not estimated. The optimizer may stop at a local maximum, so several starting values should be tried and the convergence code checked. The densities are computed with scpg.par$WBsum.method; the default suits almost all records, and if a fit is slow "direct.asymptotic" gives the same estimates faster (see dscpg).

Inflection diagnostics (TSCPG). In State $j$ the Poisson mean $\lambda_j(t)$ falls from about $K_j + M_j$ to $K_j$, and the central 80 per cent of this change takes place within $2.2/R_j$ of the inflection time $P_j/R_j$. The record observes $\lambda_j(t)$ only at $h_0 = 0, h_1, \dots, h_{T-1}$, and only through the observations that the model attributes to State $j$. The data can therefore locate $P_j/R_j$ only if the change falls inside $[0, h_{T-1}]$ and enough observations of State $j$ fall within it; otherwise the likelihood is nearly flat in $P_j/R_j$. The fitted object contains the component inflection (see Value), and a warning is given when, for some State $j$,

• less than half of the change of $\lambda_j(t)$ takes place on $[0, h_{T-1}]$ (change below 0.5). This covers every $P_j/R_j$ outside $[0, h_{T-1}]$, and a $\lambda_j(t)$ that is almost constant over the record; or

• the expected number of observations of State $j$ within the central 80 per cent of the change is below 10.

These limits come from simulations of the TSCPG with known states (see fit.scpg); since the states are hidden here, a fit that meets them may still locate $P_j/R_j$ poorly. After a warning, report $P_j/R_j$ as not identified and compare the model with SCPG states (time.scpg = FALSE) by AIC. The script verify_fit_hmm_scpg.R in the examples directory of the package reproduces the statements of this help page; run it in the console with

f <- system.file("examples", "verify_fit_hmm_scpg.R",
                 package = "HMMscpg")
source(f, echo = TRUE)

Model selection. The likelihood increases with each state added, so the number of states and the observation process are chosen by $\mathrm{AIC} = -2\log L + 2\,\mathrm{nparam}$ (Akaike, 1974), together with the decoded path (vip.hmm.scpg) and the residual process (sp.hmm.scpg), as recommended by Pohle et al. (2017).

Value

An object of the same class as hmmobj: a list with

x

NULL; the data are not kept, and must be attached as x before vip.hmm.scpg or sp.hmm.scpg is used.

mc.par, scpg.par

the estimated parameters, in the form of hmm.scpg, with the settings unchanged.

logLik

the maximized log-likelihood.

nparam

the number of estimated parameters.

convergence, message, counts

as returned by optim; convergence = 0 indicates successful convergence.

objective.error

NULL, or the last reason why the log-likelihood could not be evaluated at a trial value.

inflection

present only when time.scpg = TRUE: a data frame with the inflection diagnostics of Details, one row per State $j = 2, \dots, m$:

component, label

"lambda" and "state j".

rate, intercept

the estimates of $R_j$ and $P_j$.

inflection

the inflection time $P_j/R_j$.

first, last

$P_j$ and $P_j - R_j h_{T-1}$: the position of $P_j/R_j$ relative to 0 and to $h_{T-1}$, in units of $1/R_j$. It lies inside $[0, h_{T-1}]$ when first $\ge 0 \ge$ last.

change

the fraction, between 0 and 1, of the change of $\lambda_j(t)$ that takes place on $[0, h_{T-1}]$, $|\sigma(P_j) - \sigma(P_j - R_j h_{T-1})|$ with $\sigma(u) = 1/(1 + e^{-u})$; the first limit fails when it is below 0.5.

n

the expected number of observations in State $j$, $\sum_t \Pr(S_t = j \mid x_1, \dots, x_T)$.

in.band

the part of n from the $x_t$ whose $h_{t-1}$ lies within the central 80 per cent of the change, $|P_j - R_j h_{t-1}| < \log 9$.

rare

always NA here.

flag

the limits that are not met, or "".

References

Baum LE, Petrie T, Soules G, Weiss N (1970). “A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains.” The Annals of Mathematical Statistics, 41(1), 164-171. doi:10.1214/aoms/1177697196.

Zucchini W, MacDonald IL, Langrock R (2016). Hidden Markov Models for Time Series: An Introduction Using R, 2 edition. CRC Press, Boca Raton.

Dempster AP, Laird NM, Rubin DB (1977). “Maximum likelihood from incomplete data via the EM algorithm.” Journal of the Royal Statistical Society, Series B, 39(1), 1-22. doi:10.1111/j.2517-6161.1977.tb01600.x.

Akaike H (1974). “A new look at the statistical model identification.” IEEE Transactions on Automatic Control, 19(6), 716-723. doi:10.1109/TAC.1974.1100705.

Pohle J, Langrock R, van Beest FM, Schmidt NM (2017). “Selecting the number of states in hidden Markov models: pragmatic solutions illustrated using animal movement.” Journal of Agricultural, Biological and Environmental Statistics, 22(3), 270-293. doi:10.1007/s13253-017-0283-8.

Examples

# ==========================================================
# An example of a 2-state HMM with a TSCPG renewal process
# ==========================================================

# State 1 is the state of completeness (gamma interevent times). In the
# incomplete State 2 the number of missing events N between two observed
# events is shifted Poisson with parameter lambda_2(t), where
#   lambda_2(t) = K + M / (1 + exp(R * h_{t-1} - P))
# falls from about 5.3 at the start of the record towards K = 0.5, around
# the inflection time P/R = 300 of the last occurrence time h_{t-1}. The
# hidden chain is homogeneous, so only lambda_2(t) changes along the record.

m <- 2
tpm <- matrix(c(0.9, 0.1,
                0.2, 0.8), nrow = 2, byrow = TRUE)
mc.par <- list(initd = NULL, tpm = tpm, stationary = TRUE)
rate <- matrix(NA, nrow = 4, ncol = m - 1)
rate[1, ] <- 0.5    # K
rate[2, ] <- 5      # M
rate[3, ] <- 0.01   # R
rate[4, ] <- 3      # P: inflection time P/R = 300
shape <- 1
scale <- 1
scpg.par <- list(rate = rate, shape = shape, scale = scale,
                 rate.fix = FALSE, shape.fix = TRUE, time.scpg = TRUE,
                 WBsum.method = 'direct.evaluation')

# Defining an object for 2-state HMM with TSCPG distribution

obj_hmm2 <- hmm.scpg(x = NULL, mc.par, scpg.par)
class(obj_hmm2)
# Simulating a 2-state HMM with TSCPG distribution

sim_hmm2 <- sim.hmm.scpg(obj_hmm2, nsim = 200, seed = 1)
str(sim_hmm2[c("x", "states")])

# The record covers [0, h_T]; the TSCPG rate is evaluated on [0, h_{T-1}]
c(h_T1 = sum(sim_hmm2$x[-200]), h_T = sum(sim_hmm2$x))
# Fitting a 2-state HMM with TSCPG distributions
# (the parameters of obj_hmm2 are the starting values)

fit_hmm2 <- fit.hmm.scpg(sim_hmm2)
# transition matrix (true rows 0.9, 0.1 and 0.2, 0.8)
fit_hmm2$mc.par$tpm
fit_hmm2$scpg.par$rate[, 1]  # K, M, R, P of lambda_2(t)
fit_hmm2$scpg.par$scale      # true 1
fit_hmm2$logLik
# Inflection diagnostics of the fitted TSCPG rate lambda_2(t), whose true
# inflection time is 300. The change lies inside the record (first > 0 >
# last), almost all of it takes place on [0, h_{T-1}] (change >= 0.5), and
# about 35 expected observations of State 2 lie within its central 80 per
# cent (in.band >= 10), so no warning is given and the inflection time P/R,
# estimated at about 300, is identified.
fit_hmm2$inflection
# ========================================================
# An example of a 3-state HMM with a TSCPG renewal process
# ========================================================

# Fitting this 3-state model takes about 30 seconds on a laptop (the
# 2-state model above takes about 4 seconds), so the examples that fit
# it are wrapped in donttest.

m <- 3
initd <- c(1/3, 1/3, 1/3)
tpm <- matrix(c(0.2, 0.65, 0.15, 0.50, 0.30, 0.20, 0.25, 0.50, 0.25),
  nrow = 3, byrow = TRUE)
mc.par <- list(initd=NULL, tpm=tpm, stationary=TRUE)
rate <- matrix(NA, nrow = 4, ncol = m-1)
rate[1, ] <- c(1, 1) # K
rate[2, ] <- c(1.5, 2.5) # M
rate[3, ] <- c(0.02, 0.015) # R
rate[4, ] <- c(3, 3) # P: inflection times P/R = 150 and 200
shape <- 1
scale <- 1.7
scpg.par <- list(rate=rate, shape=shape, scale=scale,
            rate.fix=FALSE, shape.fix=TRUE, time.scpg=TRUE,
            WBsum.method='direct.evaluation')

# Defining an object for 3-state HMM with TSCPG distribution
obj_hmm2_m3 <- hmm.scpg(x = NULL, mc.par, scpg.par)
# Simulating a 3-state HMM with TSCPG distribution

sim_hmm2_m3 <- sim.hmm.scpg(obj_hmm2_m3, nsim=100, seed=9)
# Fitting a 3-state HMM with TSCPG distributions
# The fit takes about 30 seconds, so it is wrapped in donttest.

fit_hmm2_m3 <- fit.hmm.scpg(sim_hmm2_m3)
print(fit_hmm2_m3)

# Inflection diagnostics of the fitted TSCPG rates lambda_2(t) and
# lambda_3(t) (true inflection times 150 and 200). With 100 observations
# both are flagged here, and the fit warns: lambda_2(t) is fitted as almost
# constant on the record (its inflection time, 413, lies just after
# h_{T-1} = 402, and only 0.2 per cent of its change takes place on
# [0, h_{T-1}]), and lambda_3(t) as an abrupt step at about 99 with no
# expected observation of State 3 within it. Neither inflection time is
# identified; the homogeneous model (time.scpg = FALSE) may be compared by
# AIC.

fit_hmm2_m3$inflection

Description

Estimates the parameters of an HSMM with class hsmm.gamma, hsmm.scpg and hsmm.tscpg as typically returned by hsmm.scpg using direct nonlinear optimization.

Usage

fit.hsmm.scpg(
  hsmmobj,
  method = "BFGS",
  maxit = 1000,
  reltol = sqrt(.Machine$double.eps),
  negative = TRUE
)

Arguments

Details

Likelihood. The likelihood of a right-censored HSMM is evaluated by the forward algorithm of Guedon (2003). Let $u_t(j, d) = \prod_{l=t-d+1}^{t} f_j(x_l)$ denote the joint state-dependent observation probability of a visit to State $j$ of length $d$ ending at time $t$, and let

$\alpha_{t+1}^{*}(j) = \Pr(X_1, \dots, X_t, S_{t+1} = j, S_t \ne j \mid \theta) = \sum_{i \ne j} \alpha_t(i) a_{ij}$

denote the probability that a visit to State $j$ commences at time $t+1$, with $\alpha_1^{*}(j) = \pi_j$. The forward variable $\alpha_t(j) = \Pr(X_1, \dots, X_t, S_t = j, S_{t+1} \ne j \mid \theta)$, the $t$th observation being the last of a visit to State $j$, satisfies

$\alpha_t(j) = \sum_{d=1}^{t} \alpha_{t-d+1}^{*}(j)\, p_j(d)\, u_t(j, d), \quad t = 1, \dots, T-1.$

The sojourn time in the last visited state is right-censored at $T$, only its minimum being known, so the sojourn mass function is replaced by the survival function $S_j(d) = \sum_{k \ge d} p_j(k)$,

$\alpha_T(j) = \sum_{d=1}^{T} \alpha_{T-d+1}^{*}(j)\, S_j(d)\, u_T(j, d),$

and $L(\theta) = \sum_{j=1}^{m}\alpha_T(j)$. Treating the final sojourn as complete would bias the estimated expected sojourn times downwards.

Because the forward variable is indexed by state and elapsed sojourn jointly, the recursion is of order $O(T^2m^2)$ rather than the $O(Tm^2)$ of an HMM, a visit being admissible at any length up to $t$. The observation products are accumulated by the recursion $u_t(j, d) = u_{t-1}(j, d-1) f_j(x_t)$ of Mitchell and Jamieson (1993). Scaling is embedded in the recursion at each time step and the logarithms of the scale factors accumulated.

Estimation. The likelihood is maximized by direct numerical minimization of the negative log-likelihood using optim over working parameters defined on the real line. Each row of the zero-diagonal transition probability matrix is represented by multinomial logits over its $m-1$ transition probabilities, which for $m = 2$ leaves the embedded chain fully determined and contributing no parameters; the sojourn parameters are represented by logarithms, or by scaled logits where dur.lb and dur.ub are supplied; and the observation parameters are transformed as in fit.scpg. The forward recursion is evaluated by compiled code. The fitted object is returned on the natural scale and is suitable for vip.hsmm.scpg and sp.hsmm.scpg.

Inflection diagnostics (TSCPG). These apply when the observation process is a TSCPG (time.scpg = TRUE); otherwise the model has no time-dependent component and inflection is NULL. In each state $j = 2, \dots, m$ the rate $\lambda_j(t) = K_j + M_j/\{1 + \exp(R_j h_{t-1} - P_j)\}$ changes around the inflection time $h_{t-1} = P_j/R_j$. Over the record the last occurrence time $h_{t-1}$ runs from $h_0 = 0$ ($t = 1$) to $h_{T-1} = x_1 + \dots + x_{T-1}$ ($t = T$), so $\lambda_j(t)$ is observed only on $[0, h_{T-1}]$, and $P_j/R_j$ is identified only when this interval covers the change and enough observations of state $j$ fall within it. Writing $u_t = P_j - R_j h_{t-1}$, the sensitivity of $\lambda_j(t)$ to $P_j/R_j$ is at most $R_j e^{-|u_t|}$, so observations more than a few multiples of $1/R_j$ from the inflection time carry almost no information about it; when $P_j/R_j$ lies $\delta$ multiples of $1/R_j$ outside $[0, h_{T-1}]$, the profile information decreases like $e^{-4\delta}$, by a factor of about 45 per multiple when the states are known. The state being hidden, only the observations attributed to state $j$ inform the curve, and less than they would if the states were known.

The fitted object contains the component inflection (see Value). Its limits are that at least half of the change of $\lambda_j(t)$ takes place on $[0, h_{T-1}]$ (change $\ge 0.5$), which places $P_j/R_j$ inside this interval, and that in.band is at least 10. They were derived from profile likelihoods of $P/R$ in simulations of 250 TSCPG observations (see fit.scpg). With hidden states they are necessary rather than sufficient.

A warning is issued when a limit is not met. The inflection time should then be reported as not identified rather than by its estimate, and the model with SCPG states (time.scpg = FALSE) compared by AIC. When a fitted $\lambda_j(t)$ is nearly constant on $[0, h_{T-1}]$, $P_j/R_j$ is not identified under the hypothesis of no change and the chi-squared calibration of likelihood-ratio tests and profile intervals fails.

The script verify_fit_hsmm_scpg.R in the examples directory of the package reproduces the numbers stated on this help page; run it in the console with

f <- system.file("examples", "verify_fit_hsmm_scpg.R",
                 package = "HMMscpg")
source(f, echo = TRUE)

Value

The returned object provides estimated HSMM with additional information about the function evaluation and gradient counts at the estimated minimum of negative log-likelihood function, convergence code and minimized likelihood function value as returned by optim and nparam is the number of parameters to be estimated. With TSCPG renewal processes it also contains inflection, a data frame with the inflection diagnostics of Details, one row per state $j = 2, \dots, m$:

component, label

"lambda" and "state j".

rate, intercept

$R_j$ and $P_j$.

inflection

$P_j/R_j$.

first, last

$u_1 = P_j$ and $u_T = P_j - R_j h_{T-1}$, the values of $u_t$ for the first and last interevent times, at $h_0 = 0$ and $h_{T-1}$. The inflection time lies in $(0, h_{T-1})$ when first $> 0 >$ last; otherwise -first is its distance before 0, or last its distance after $h_{T-1}$, in units of $1/R_j$.

change

the part of the change of $\lambda_j(t)$, between 0 and 1, that takes place on $[0, h_{T-1}]$, $|\sigma(u_1) - \sigma(u_T)|$ with $\sigma(u) = 1/(1 + e^{-u})$.

n

the expected number of observations in state $j$, $\sum_{t=1}^{T} \Pr(S_t = j \mid x)$, from the smoothed state probabilities.

in.band

the part of n from the observations $x_t$ whose $h_{t-1}$ lies in the central 10 to 90 per cent of the change, where $|u_t| < \log 9$.

rare

NA; used by fit.ihmm.scpg.

flag

the limits of Details that are not met, or "".

References

Guédon Y (2003). “Estimating hidden semi-Markov chains from discrete sequences.” Journal of Computational and Graphical Statistics, 12(3), 604-639. doi:10.1198/1061860032030.

Mitchell CD, Jamieson LH (1993). “Modeling duration in a hidden Markov model with the exponential family.” In Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing, 331-334. doi:10.1109/ICASSP.1993.319304.

Examples

# ==========================================================
# An example of a 2-state HSMM with an SCPG renewal process
# ==========================================================

# State 1 is the state of completeness (gamma interevent times with shape
# alpha); in the incomplete State 2 the number N of events missing between
# two observed events is shifted Poisson with parameter lambda_2, so the
# observed events are N + 1 such gamma interevent times apart. The sojourns
# are shifted Poisson with
# parameters dur (expected sojourns 1 + dur = 5 and 4 observations). With
# m = 2 the zero-diagonal transition matrix is fixed: the chain alternates
# between the two states.

m <- 2
tpm <- matrix(c(0, 1,
                1, 0), nrow = 2, byrow = TRUE)
smc.par <- list(initd = NULL, tpm = tpm, dur = c(4, 3),
                dur.lb = NA, dur.ub = NA, stationary = TRUE)
rate <- matrix(3, nrow = 1, ncol = m - 1)   # lambda_2
shape <- 1.5                                # common shape alpha
scale <- 1
scpg.par <- list(rate = rate, shape = shape, scale = scale,
                 rate.fix = FALSE, shape.fix = FALSE, time.scpg = FALSE,
                 WBsum.method = 'direct.evaluation')

# Defining an object of 2-state HSMM with SCPG distribution

obj_hsmm1 <- hsmm.scpg(x = NULL, smc.par, scpg.par)
class(obj_hsmm1)
# Simulating a 2-state HSMM with SCPG distribution

sim_hsmm1 <- sim.hsmm.scpg(obj_hsmm1, nsim = 150, seed = 1)
str(sim_hsmm1[c("x", "states")])

# Number of visits to each state in the simulated path
table(rle(sim_hsmm1$states)$values)
# Fitting a 2-state HSMM with SCPG distributions
# (the parameters of obj_hsmm1 are the starting values)

fit_hsmm1 <- fit.hsmm.scpg(sim_hsmm1)
fit_hsmm1$smc.par$dur      # sojourn parameters mu_1, mu_2 (true 4, 3)
fit_hsmm1$scpg.par$rate    # lambda_2 (true 3)
fit_hsmm1$scpg.par$shape   # alpha (true 1.5)
fit_hsmm1$scpg.par$scale   # beta (true 1)
fit_hsmm1$logLik
# =========================================================
# An example of a 3-state HSMM with an SCPG renewal process
# =========================================================

# Fitting this 3-state model takes about 15 seconds on a laptop (the
# 2-state model above takes about 1 second), so the examples that fit
# it are wrapped in donttest.

m <- 3
tpm <- matrix(c(0, 0.6, 0.4,
                0.5, 0, 0.5,
                0.3, 0.7, 0), nrow = 3, byrow = TRUE)
smc.par <- list(initd=NULL, tpm=tpm, dur=c(3, 5, 2),
            dur.lb=NA, dur.ub=NA, stationary=TRUE)
rate <- matrix(c(8, 5), nrow = 1, ncol = m-1)
shape <- c(2, 1, 1.7)
scale <- 1
scpg.par <- list(rate=rate, shape=shape, scale=scale,
            rate.fix=FALSE, shape.fix=FALSE, time.scpg=FALSE,
            WBsum.method='direct.evaluation')

# Defining an object of 3-state HSMM with SCPG distribution

obj_hsmm1_m3 <- hsmm.scpg(x = NULL, smc.par, scpg.par)
# Simulating a 3-state HSMM with SCPG distribution

sim_hsmm1_m3 <- sim.hsmm.scpg(obj_hsmm1_m3, nsim = 80, seed = 1)
print(sim_hsmm1_m3)
# Fitting a 3-state HSMM with SCPG distributions
# The fit takes about 15 seconds, so it is wrapped in donttest.

fit_hsmm1_m3 <- fit.hsmm.scpg(sim_hsmm1_m3)
print(fit_hsmm1_m3)

# ===========================================================
# An example of a 2-state HSMM with a TSCPG renewal process
# ===========================================================

# State 1 is the state of completeness (gamma interevent times). In the
# incomplete State 2 the number of missing events N between two observed
# events is shifted Poisson with parameter lambda_2(t), where
#   lambda_2(t) = K + M / (1 + exp(R * h_{t-1} - P))
# falls from about 5.8 at the start of the record towards K = 0.5, around
# the inflection time P/R = 400 of the last occurrence time h_{t-1}. The
# sojourns are shifted Poisson with parameters dur (expected sojourns
# 1 + dur = 5 and 4 observations). With m = 2 the zero-diagonal transition
# matrix is fixed: the chain alternates between the two states.

m <- 2
tpm <- matrix(c(0, 1,
                1, 0), nrow = 2, byrow = TRUE)
smc.par <- list(initd = NULL, tpm = tpm, dur = c(4, 3),
                dur.lb = NA, dur.ub = NA, stationary = TRUE)
rate <- matrix(NA, nrow = 4, ncol = m - 1)
rate[1, ] <- 0.5    # K
rate[2, ] <- 6      # M
rate[3, ] <- 0.005  # R
rate[4, ] <- 2      # P: inflection time P/R = 400
shape <- 1
scale <- 1
scpg.par <- list(rate = rate, shape = shape, scale = scale,
                 rate.fix = FALSE, shape.fix = TRUE, time.scpg = TRUE,
                 WBsum.method = 'direct.evaluation')

# Defining an object of 2-state HSMM with TSCPG distribution

obj_hsmm2 <- hsmm.scpg(x = NULL, smc.par, scpg.par)
class(obj_hsmm2)
# Simulating a 2-state HSMM with a TSCPG distribution

sim_hsmm2 <- sim.hsmm.scpg(obj_hsmm2, nsim = 300, seed = 3)
str(sim_hsmm2[c("x", "states")])

# The record covers [0, h_T]; the TSCPG rate is evaluated on [0, h_{T-1}]
c(h_T1 = sum(sim_hsmm2$x[-300]), h_T = sum(sim_hsmm2$x))
# The fit takes several seconds, so the code below is wrapped in donttest.

# Fitting a 2-state HSMM with a TSCPG distribution
# (the parameters of obj_hsmm2 are the starting values)

fit_hsmm2 <- fit.hsmm.scpg(sim_hsmm2)
fit_hsmm2$smc.par$dur        # sojourn parameters mu_1, mu_2 (true 4, 3)
fit_hsmm2$scpg.par$rate      # K, M, R, P of lambda_2(t)
fit_hsmm2$scpg.par$scale     # true 1
fit_hsmm2$logLik

# The fit takes several seconds, so the code below is wrapped in donttest.

# Inflection diagnostics of the fitted TSCPG rate lambda_2(t), whose true
# inflection time is 400. The change lies inside the record (first > 0 >
# last), most of it takes place on [0, h_{T-1}] (change >= 0.5), and the
# expected number of observations of State 2 within its central 80 per cent
# (in.band, from the smoothed state probabilities) is well above 10, so no
# warning is given and the inflection time P/R, estimated at about 412, is
# identified.
fit_hsmm2$inflection

# ===========================================================
# An example of a 3-state HSMM with a TSCPG renewal process
# ===========================================================

# Fitting this 3-state model takes about 2 minutes on a laptop (the
# 2-state model above takes about 7 seconds), so the examples that fit
# it are wrapped in donttest.

# State 1 is the state of completeness (gamma interevent times). In the
# incomplete States 2 and 3 the number of missing events N between two
# observed events is shifted Poisson with parameter lambda_j(t), where
#   lambda_j(t) = K_j + M_j / (1 + exp(R_j * h_{t-1} - P_j))
# falls around the inflection time P_j/R_j of the last occurrence time
# h_{t-1}: lambda_2(t) from about 3.4 towards 0.5 around 600, and
# lambda_3(t) from about 10 towards 4 around 900. The sojourns are shifted
# Poisson with parameters dur (expected sojourns 1 + dur = 5, 4 and 4
# observations), and the zero-diagonal transition matrix gives the next
# state of each visit. The shape 2 makes the interevent times of the three
# states less overlapping than with shape 1.

m <- 3
tpm <- matrix(c(0, 0.6, 0.4,
                0.5, 0, 0.5,
                0.6, 0.4, 0), nrow = 3, byrow = TRUE)
smc.par <- list(initd = NULL, tpm = tpm, dur = c(4, 3, 3),
                dur.lb = NA, dur.ub = NA, stationary = TRUE)
rate <- matrix(NA, nrow = 4, ncol = m - 1)
rate[1, ] <- c(0.5, 4)       # K
rate[2, ] <- c(3, 6)         # M
rate[3, ] <- c(0.006, 0.006) # R
rate[4, ] <- c(3.6, 5.4)     # P: inflection times P/R = 600 and 900
shape <- 2
scale <- 0.5
scpg.par <- list(rate = rate, shape = shape, scale = scale,
                 rate.fix = FALSE, shape.fix = TRUE, time.scpg = TRUE,
                 WBsum.method = 'direct.evaluation')

# Defining an object of 3-state HSMM with TSCPG distribution

obj_hsmm2_m3 <- hsmm.scpg(x = NULL, smc.par, scpg.par)
class(obj_hsmm2_m3)
# Simulating a 3-state HSMM with a TSCPG distribution

sim_hsmm2_m3 <- sim.hsmm.scpg(obj_hsmm2_m3, nsim = 400, seed = 3)
str(sim_hsmm2_m3[c("x", "states")])

# The record covers [0, h_T]; the TSCPG rates are evaluated on [0, h_{T-1}]
c(h_T1 = sum(sim_hsmm2_m3$x[-400]), h_T = sum(sim_hsmm2_m3$x))

# Number of visits to each state in the simulated path
table(rle(sim_hsmm2_m3$states)$values)
# Fitting a 3-state HSMM with a TSCPG distribution
# (the parameters of obj_hsmm2_m3 are the starting values). The fit
# takes about 2 minutes, so it is wrapped in donttest.

fit_hsmm2_m3 <- fit.hsmm.scpg(sim_hsmm2_m3)
fit_hsmm2_m3$smc.par$tpm        # transition matrix of the embedded chain
# sojourn parameters mu_1, mu_2, mu_3 (true 4, 3, 3)
fit_hsmm2_m3$smc.par$dur
fit_hsmm2_m3$scpg.par$rate      # K, M, R, P of lambda_2(t) and lambda_3(t)
fit_hsmm2_m3$scpg.par$scale     # true 0.5
fit_hsmm2_m3$logLik

# Inflection diagnostics of the fitted TSCPG rates lambda_2(t) and
# lambda_3(t), whose true inflection times are 600 and 900. For each rate
# the change lies inside the record (first > 0 > last), most of it takes
# place on [0, h_{T-1}] (change >= 0.5), and the expected number of
# observations of the state within its central 80 per cent (in.band, from
# the smoothed state probabilities) is above 10, so no warning is given and
# both inflection times P_j/R_j are identified.

fit_hsmm2_m3$inflection

Description

Estimates the parameters of an IHMM with class ihmm.gamma, ihmm.scpg and ihmm.tscpg as typically returned by ihmm.scpg using direct nonlinear optimization.

Usage

fit.ihmm.scpg(
  ihmmobj,
  method = "BFGS",
  maxit = 1000,
  reltol = sqrt(.Machine$double.eps),
  negative = TRUE
)

Arguments

Details

Likelihood. The likelihood of an IHMM is obtained from that of an HMM by replacing the transition probabilities with their time-varying counterparts (Zucchini et al., 2016),

$L(\theta) = \Pr(X_{1:T} \mid H_{1:T}; \theta) = \sum_{s_1, \dots, s_T} \pi_{s_1} f_{s_1}(x_1) \prod_{t=2}^{T} a_{s_{t-1}s_t}(t) f_{s_t}(x_t),$

where $\theta = \{\pi_j, r_{ij}, w_{ij}, \theta_j\}$. With the forward variable $\alpha_t(j) = \Pr(X_1, \dots, X_t, S_t = j \mid H_{t-1}; \theta)$, the recursion becomes

$\alpha_1(j) = \pi_j f_j(x_1), \qquad \alpha_t(j) = \left\{\sum_{i=1}^{m}\alpha_{t-1}(i)a_{ij}(t)\right\} f_j(x_t), \quad t = 2, \dots, T,$

and $L(\theta) = \sum_{j=1}^{m}\alpha_T(j)$. The transition probability matrix is reconstructed from its logistic coefficients at each step, using the occurrence time accumulated from the preceding observations; the likelihood is consequently not invariant under permutation of the observations. As in the homogeneous case, the recursion is scaled at each step and the logarithms of the scale factors accumulated.

Estimation. The likelihood is maximized by direct numerical minimization of the negative log-likelihood using optim over working parameters defined on the real line. The transition structure is described by $2m(m-1)$ logistic coefficients. Each pair $(r_{ij}, w_{ij})$ enters the log-odds $w_{ij} - r_{ij} h_{t-1}$, evaluated at the last occurrence time $h_{t-1} = x_1 + \dots + x_{t-1}$, and is represented by $\log(r_{ij} h_{T-1})$ and $w_{ij}/(r_{ij} h_{T-1})$, the time $w_{ij}/r_{ij}$ as a fraction of $h_{T-1} = x_1 + \dots + x_{T-1}$, the value of $h_{t-1}$ for the last interevent time $t = T$. The estimates of $\log r_{ij}$ and $w_{ij}$ are almost perfectly correlated whenever $w_{ij}/r_{ij}$ lies in $(0, h_{T-1})$, whereas these working parameters are much less correlated and do not depend on the unit of time. The initial distribution is represented by multinomial logits and the observation parameters as in fit.scpg. The inverse transformations are applied to the minimizer, so that the fitted object is a valid ihmm.scpg object.

Inflection diagnostics. The transition probabilities change along the record through the log-odds $u_t = w_{ij} - r_{ij} h_{t-1}$, which is zero at $h_{t-1} = w_{ij}/r_{ij}$: for $m = 2$ the inflection time of $a_{i2}(t)$, and for $m \ge 3$ the time at which $a_{ij}(t) = a_{i1}(t)$. The departure at time $t$ uses $a_{ij}(t)$ at $h_{t-1}$, so over the record $a_{ij}(t)$ is observed for $h_{t-1}$ between $h_1$ ($t = 2$) and $h_{T-1} = x_1 + \dots + x_{T-1}$ ($t = T$). The information about $(r_{ij}, w_{ij})$ comes from the departures from state $i$, each weighted by $a_{ij}(t)\{1 - a_{ij}(t)\}$: where the transition probability is close to 0 or 1 the observed transitions are almost all of one kind. The likelihood is therefore nearly flat in $w_{ij}/r_{ij}$ when the rarer outcome of a departure from state $i$ (moving to $j$, or not) is expected fewer than about 10 times, and exactly flat in one direction when it is never observed. Because the log-odds is linear, the profile information about $w_{ij}/r_{ij}$ decreases more slowly than for the four-parameter curves of the TSCPG and the IHSMM: with the states known, by a factor of about 3.7 (against about 45) for each multiple $\delta$ of $1/r_{ij}$ by which it lies outside $[0, h_{T-1}]$. When the states are TSCPG (time.scpg = TRUE), the rates $\lambda_j(t)$ are checked as in fit.hmm.scpg; how the two kinds of curve affect each other is described at the end of Details.

The fitted object contains the component inflection (see Value). The limits for a transition are that at least half of the logistic change takes place on $[0, h_{T-1}]$ (change $\ge 0.5$), which places $w_{ij}/r_{ij}$ inside this interval, and that rare is at least 10; the limits for $\lambda_j(t)$ are those of fit.hmm.scpg. In simulations of a two-state IHMM with 300 observations the 95 per cent profile interval of $w_{12}/r_{12}$ was bounded when the rarer outcome of a departure from State 1 was expected 8.6 times (at the true values) and unbounded when it was expected 5.4 times or fewer, even with 293 departures.

A warning is issued when a limit is not met. The crossing time should then be reported as not identified rather than by its estimate, and a model with fewer time-dependent parameters, such as the homogeneous fit.hmm.scpg, compared by AIC. When a fitted transition probability is nearly constant on $[0, h_{T-1}]$, its crossing time is not identified under the hypothesis of no change and the chi-squared calibration of likelihood-ratio tests and profile intervals fails.

Two kinds of time dependence (class ihmm.tscpg). With time.scpg = TRUE the record can become more complete in two ways: the transition probabilities $a_{ij}(t)$ govern how often the record is in an incomplete state, and the Poisson means $\lambda_j(t)$ how many events are missed while it is there. Both change the mean interevent time, $\alpha\beta$ in State 1 and $\alpha\beta\{\lambda_j(t) + 2\}$ in State $j$, so the data tell them apart only through the shape of the distribution of $x_t$ and the persistence of the states. This affects their inflection times in two ways. First, $\lambda_j(t)$ is informed only by the observations of State $j$: if the transitions make State $j$ rare where $\lambda_j(t)$ changes, few observations inform it, which in.band shows, since it counts the expected observations of State $j$. Second, as $\lambda_j(t)$ falls, the interevent times of State $j$ (sums of at least two complete ones) come closer to those of State 1, the state probabilities become less certain, and every curve is then less precisely determined, which the counts in.band and rare do not show. In simulations of a two-state IHMM with 500 observations, whose $a_{12}(t)$ and $a_{22}(t)$ change around 304 and 600 and whose $\lambda_2(t)$ falls between the levels 5.5 and 0.5 around 200, 400 or 600 (30 records each), the estimated inflection time of $\lambda_2(t)$ and the crossing times of $a_{12}(t)$ and $a_{22}(t)$ were almost uncorrelated (Spearman correlations between $-0.01$ and 0.21), also when $\lambda_2(t)$ and $a_{22}(t)$ changed at the same time. Their precision depended on how long the two states stayed distinct: with the inflection time of $\lambda_2(t)$ at 200 rather than 600, the spread (interquartile range divided by 1.349) of the estimated inflection times was 124 rather than 49 for $\lambda_2(t)$ and 149 rather than 49 for $a_{22}(t)$, and the mean absolute error of the proportion of events recorded, $p_t = \Pr(S_t = 1 \mid x) + \sum_{j \ge 2}\Pr(S_t = j \mid x)/\{2 + \lambda_j(t)\}$, averaged over each quarter of the record, was 0.060 rather than 0.030, although in.band was about the same (medians 49 and 53). When both kinds are fitted, compare the model by AIC with those that have only one: fit.hmm.scpg with time.scpg = TRUE (only $\lambda_j(t)$) and fit.ihmm.scpg with time.scpg = FALSE (only $a_{ij}(t)$); and report $p_t$, which combines them, alongside the inflection times.

The script verify_fit_ihmm_scpg.R in the examples directory of the package reproduces the numbers stated on this help page; run it in the console with

f <- system.file("examples", "verify_fit_ihmm_scpg.R",
                 package = "HMMscpg")
source(f, echo = TRUE)

Value

The returned object provides estimated HMM with additional information about the function evaluation and gradient counts at the estimated minimum of negative log-likelihood function, convergence code and minimized likelihood function value as returned by optim and nparam is the number of parameters to be estimated. It also contains inflection, a data frame with the inflection diagnostics of Details, with one row for each transition $i \to j$, $j = 2, \dots, m$, and, for TSCPG states, one row for each $\lambda_j(t)$. For the transitions:

component, label

"a" and "i -> j".

rate, intercept

$r_{ij}$ and $w_{ij}$.

inflection

$w_{ij}/r_{ij}$.

first, last

$w_{ij}$ and $u_T = w_{ij} - r_{ij} h_{T-1}$, the log-odds at $h_{t-1} = 0$ and at $h_{T-1}$. The crossing lies in $(0, h_{T-1})$ when first $> 0 >$ last; otherwise -first is its distance before 0, or last its distance after $h_{T-1}$, in units of $1/r_{ij}$.

change

the part of the logistic change, between 0 and 1, that takes place on $[0, h_{T-1}]$.

n

the expected number of departures from state $i$, $\sum_{t=2}^{T}\Pr(S_{t-1} = i \mid x)$, from the smoothed state probabilities.

in.band

the part of n from departures whose log-odds $u_t$ lies in $(-\log 9, \log 9)$.

rare

the expected number of the rarer outcome of a departure from state $i$, $\sum_{t=2}^{T}\Pr(S_{t-1} = i \mid x)\min\{a_{ij}(t), 1 - a_{ij}(t)\}$.

flag

the limits of Details that are not met, or "".

References

Zucchini W, MacDonald IL, Langrock R (2016). Hidden Markov Models for Time Series: An Introduction Using R, 2 edition. CRC Press, Boca Raton.

Examples

# ==========================================================
# An example of a 2-state IHMM with an SCPG renewal process
# ==========================================================

# The transition probabilities a_12(t) and a_22(t) into the incomplete
# State 2 decrease along the record: the log-odds of moving from state i to
# State 2 rather than to State 1 is w_i2 - r_i2 * h_{t-1}, which is zero at
# the crossing time w_i2 / r_i2 (250 for a_12(t) and 500 for a_22(t)). In
# State 2 the number N of events missing between two observed events is
# shifted Poisson with parameter lambda_2. The chain starts mostly in State 2
# (initd = (0.2, 0.8)),
# since the early part of a record is usually the least complete.

m <- 2
initd <- c(0.2, 0.8)
ttpm <- array(NA, dim = c(m, m-1, 2))
ttpm[, , 1] <- c(0.004, 0.006)   # rates of change r_12, r_22
ttpm[, , 2] <- c(1, 3)           # intercepts w_12, w_22
imc.par <- list(initd=initd, ttpm=ttpm, stationary=FALSE)
rate <- matrix(3, nrow = 1, ncol = m-1)   # lambda_2
shape <- 1.5
scale <- 1
scpg.par <- list(rate=rate, shape=shape, scale=scale,
            rate.fix=FALSE, shape.fix=FALSE, time.scpg=FALSE,
            WBsum.method='direct.evaluation')

# Defining an object of 2-state IHMM with SCPG distribution

obj_ihmm1 <- ihmm.scpg(x = NULL, imc.par, scpg.par)
class(obj_ihmm1)
# Simulating a 2-state IHMM with SCPG distribution

sim_ihmm1 <- sim.ihmm.scpg(obj_ihmm1, nsim = 200, seed = 1)
str(sim_ihmm1[c("x", "states")])

# Observations in each state: State 1 becomes more frequent as the
# transitions into State 2 become less likely
table(state = sim_ihmm1$states,
      half = ifelse(seq_along(sim_ihmm1$states) <= 100, "first", "second"))
# The fit takes about 6 s, so the code below is wrapped in donttest.

# Fitting a 2-state IHMM with SCPG distributions
# (the parameters of obj_ihmm1 are the starting values)

fit_ihmm1 <- fit.ihmm.scpg(sim_ihmm1)
fit_ihmm1$imc.par$ttpm[, 1, ]   # r_i2 (first column) and w_i2 (second)
fit_ihmm1$scpg.par$rate         # lambda_2 (true 3)
fit_ihmm1$scpg.par$shape        # true 1.5
fit_ihmm1$scpg.par$scale        # true 1
fit_ihmm1$logLik

# The fit takes about 6 s, so the code below is wrapped in donttest.

# Inflection diagnostics of the fitted transition probabilities a_12(t) and
# a_22(t) (true crossing times 250 and 500; labelled "i -> j"). Each curve
# is checked against the limits change >= 0.5 and rare >= 10 expected
# departures of the rarer kind (moving to State 2, or not). Here both are
# met (rare about 19 and 14), so no warning is given; the crossing times
# are estimated at about 442 and 486.
fit_ihmm1$inflection

# =========================================================
# An example of a 3-state IHMM with an SCPG renewal process
# =========================================================

# Fitting this 3-state model takes about 1.5 minutes on a laptop (the
# 2-state model above takes about 6 seconds), so the examples that fit
# it are wrapped in donttest.

m <- 3
initd <- c(0.45, 0.25, 0.30)
ttpm <- array(NA, dim = c(m, m-1, 2))
# rates of change r_ij (rows i = 1, 2, 3; columns j = 2, 3)
ttpm[, , 1] <- matrix(c(0.0008, 0.0010, 0.0006,
                        0.0006, 0.0008, 0.0010), nrow = m, ncol = m-1)
# intercepts w_ij; a_ij(t) = a_i1(t) at h = w_ij/r_ij, between 2000 and 3000
ttpm[, , 2] <- ttpm[, , 1] * matrix(c(2500, 2000, 3000,
                                      2000, 3000, 2500), nrow = m, ncol = m-1)
stationary <- FALSE
imc.par <- list(initd=initd, ttpm=ttpm, stationary=stationary)
rate <- matrix(c(7, 2), nrow = 1, ncol = m-1)
shape <- 1.7
scale <- 3.5
scpg.par <- list(rate=rate, shape=shape, scale=scale,
            rate.fix=FALSE, shape.fix=FALSE, time.scpg=FALSE,
            WBsum.method='direct.evaluation')

# Defining an object of 3-state IHMM with SCPG distribution

obj_ihmm1_m3 <- ihmm.scpg(x = NULL, imc.par, scpg.par)
# Simulating a 3-state IHMM with SCPG distribution

sim_ihmm1_m3 <- sim.ihmm.scpg(obj_ihmm1_m3, nsim = 200, seed = 2)
print(sim_ihmm1_m3)
# Fitting a 3-state IHMM with SCPG distributions
# The fit takes about 1.5 minutes, so it is wrapped in donttest.

fit_ihmm1_m3 <- fit.ihmm.scpg(sim_ihmm1_m3)
print(fit_ihmm1_m3)

# Inflection diagnostics of the fitted transition probabilities (true
# crossing times between 2000 and 3000). The rates r_ij of this example are
# small (0.0006 to 0.001), so each change is wide compared with the record
# (h_{T-1} = 4733), and with 200 observations only a_12(t) meets the
# limits here: the crossing times of a_13(t), a_22(t), a_32(t) and a_33(t)
# lie outside the record, with less than half of the change on
# [0, h_{T-1}], and a_13(t), a_22(t) and a_23(t) have fewer than 10
# expected departures of the rarer kind. The fit warns, and the flagged
# crossing times should be reported as not identified.

fit_ihmm1_m3$inflection

# ==========================================================
# An example of a 2-state IHMM with a TSCPG renewal process
# ==========================================================

# Two kinds of time dependence describe how completeness improves along the
# record (see the Details of ?fit.ihmm.scpg):
#  - the transition probabilities a_12(t) and a_22(t) into the incomplete
#    State 2 decrease: the log-odds of moving from state i to State 2
#    rather than to State 1 is w_i2 - r_i2 * h_{t-1}, which is zero at the
#    crossing time w_i2 / r_i2 (about 304 for a_12(t) and 600 for a_22(t));
#  - the Poisson mean lambda_2(t) of the TSCPG falls from about 5.4 at the
#    start of the record towards K = 0.5, around the inflection time
#    P/R = 400.
# The chain starts mostly in State 2 (initd = (0.2, 0.8)), since the early
# part of a record is usually the least complete.

m <- 2
initd <- c(0.2, 0.8)
ttpm <- array(NA, dim = c(m, m-1, 2))
ttpm[, , 1] <- c(0.0028, 0.005)   # rates of change r_12, r_22
ttpm[, , 2] <- c(0.85, 3)         # intercepts w_12, w_22
imc.par <- list(initd=initd, ttpm=ttpm, stationary=FALSE)
rate <- matrix(c(0.5, 5, 0.01, 4), nrow = 4)   # K, M, R, P: P/R = 400
shape <- 1
scale <- 1
scpg.par <- list(rate=rate, shape=shape, scale=scale,
            rate.fix=FALSE, shape.fix=TRUE, time.scpg=TRUE,
            WBsum.method='direct.evaluation')

# Defining an object of 2-state IHMM with TSCPG distribution

obj_ihmm2 <- ihmm.scpg(x = NULL, imc.par, scpg.par)
class(obj_ihmm2)
# Simulating a 2-state IHMM with a TSCPG distribution

sim_ihmm2 <- sim.ihmm.scpg(obj_ihmm2, nsim = 200, seed = 2)
str(sim_ihmm2[c("x", "states")])

# Observations in each state: State 1 becomes more frequent as the
# transitions into State 2 become less likely
table(state = sim_ihmm2$states,
      half = ifelse(seq_along(sim_ihmm2$states) <= 100, "first", "second"))
# The fit takes several seconds, so the code below is wrapped in donttest.

# Fitting a 2-state IHMM with a TSCPG distribution
# (the parameters of obj_ihmm2 are the starting values)

fit_ihmm2 <- fit.ihmm.scpg(sim_ihmm2)
fit_ihmm2$imc.par$ttpm[, 1, ]   # r_i2 (first column) and w_i2 (second)
fit_ihmm2$scpg.par$rate[, 1]    # K, M, R, P of lambda_2(t)
fit_ihmm2$scpg.par$scale        # true 1
fit_ihmm2$logLik

# The fit takes several seconds, so the code below is wrapped in donttest.

# Inflection diagnostics of the three fitted logistic curves: lambda_2(t)
# (true inflection time 400; component "lambda") and a_12(t), a_22(t)
# (true crossing times 304 and 600; component "a", labelled "i -> j").
# Each curve is checked against its own limits: change >= 0.5 for all of
# them, in.band >= 10 expected observations of State 2 for lambda_2(t), and
# rare >= 10 expected departures of the rarer kind for a_ij(t). Here every
# limit is met, so no warning is given. lambda_2(t) is located well (402),
# but a_12(t) changes slowly (1/r_12 is about 360 time units) and its
# crossing time is estimated at 542, while that of a_22(t), at 637, lies
# near the end of the record (h_{T-1} = 691), with only 55 per cent of its
# change on [0, h_{T-1}]: identified is not the same as precisely
# estimated. How the two kinds of curve affect each other is described in
# the Details of ?fit.ihmm.scpg.
fit_ihmm2$inflection

# =========================================================
# An example of a 3-state IHMM with a TSCPG renewal process
# =========================================================

# Fitting this 3-state model takes about 2.5 minutes on a laptop (the
# 2-state model above takes about 8 seconds), so the examples that fit
# it are wrapped in donttest.

# Two kinds of time dependence describe how completeness improves along the
# record (see the Details of ?fit.ihmm.scpg):
#  - the transition probabilities a_ij(t): moves into the incomplete
#    States 2 and 3 become less likely as time goes on. The log-odds of
#    moving from state i to state j rather than to State 1 is
#    w_ij - r_ij * h_{t-1}, which is zero at the crossing time w_ij / r_ij;
#  - the Poisson means lambda_j(t) of the TSCPG: in the incomplete states
#    fewer events are missed after the inflection times P_j / R_j.
# The shape 2 makes the interevent times of the three states less
# overlapping than with shape 1, so that the states are easier to tell apart.

m <- 3
initd <- c(0.4, 0.3, 0.3)
ttpm <- array(NA, dim = c(m, m-1, 2))
# rates of change r_ij (rows i = 1, 2, 3; columns j = 2, 3)
ttpm[, , 1] <- matrix(c(0.006, 0.006, 0.006,
                        0.006, 0.006, 0.006), nrow = m, ncol = m-1)
# intercepts w_ij; crossing times w_ij/r_ij between 400 and 600
ttpm[, , 2] <- ttpm[, , 1] * matrix(c(500, 600, 500,
                                      400, 500, 600), nrow = m, ncol = m-1)
imc.par <- list(initd=initd, ttpm=ttpm, stationary=FALSE)
rate <- matrix(NA, nrow = 4, ncol = m-1)
rate[1, ] <- c(0.5, 4)       # K: lambda_2(t) falls from about 3.4 towards
rate[2, ] <- c(3, 6)         # M: 0.5, lambda_3(t) from about 10 towards 4
rate[3, ] <- c(0.006, 0.006) # R
rate[4, ] <- c(3, 4.2)       # P: inflection times P/R = 500 and 700
shape <- 2
scale <- 0.5
scpg.par <- list(rate=rate, shape=shape, scale=scale,
            rate.fix=FALSE, shape.fix=TRUE, time.scpg=TRUE,
            WBsum.method='direct.evaluation')

# Defining an object of 3-state IHMM with TSCPG distribution

obj_ihmm2_m3 <- ihmm.scpg(x = NULL, imc.par, scpg.par)
class(obj_ihmm2_m3)
# Simulating a 3-state IHMM with a TSCPG distribution

sim_ihmm2_m3 <- sim.ihmm.scpg(obj_ihmm2_m3, nsim = 500, seed = 1)
str(sim_ihmm2_m3[c("x", "states")])

# Observations in each state: State 1 becomes more frequent as the
# transitions into the incomplete states become less likely
table(state = sim_ihmm2_m3$states,
      half = ifelse(seq_along(sim_ihmm2_m3$states) <= 250, "first", "second"))
# Fitting a 3-state IHMM with a TSCPG distribution
# (the parameters of obj_ihmm2_m3 are the starting values). The fit
# takes about 2.5 minutes, so it is wrapped in donttest.

fit_ihmm2_m3 <- fit.ihmm.scpg(sim_ihmm2_m3)
fit_ihmm2_m3$imc.par$ttpm[, , 1]   # r_ij (rows i, columns j = 2, 3)
fit_ihmm2_m3$imc.par$ttpm[, , 2]   # w_ij
fit_ihmm2_m3$scpg.par$rate         # K, M, R, P of lambda_2(t) and lambda_3(t)
fit_ihmm2_m3$logLik

# Inflection diagnostics of the eight fitted logistic curves: lambda_2(t)
# and lambda_3(t) (true inflection times 500 and 700; component "lambda")
# and the six transition curves (true crossing times 400 to 600; component
# "a", labelled "i -> j"). Each curve is checked against its own limits:
# change >= 0.5 for all of them, in.band >= 10 expected observations of the
# state for lambda_j(t), and rare >= 10 expected departures of the rarer
# kind for a_ij(t).
#
# Here lambda_2(t) and the six transition curves meet their limits, while
# lambda_3(t) is flagged: it is fitted as a step, with no observation of
# State 3 inside its change. This shows how the two kinds of time
# dependence interact: lambda_3(t) is informed only by the observations of
# State 3, and the transitions make State 3 rare later in the record (in
# the simulated path, 31 of the 70 observations before h = 500 but only 16
# of the 374 after h = 700, where lambda_3(t) changes), so the fit can put
# the whole change between two of them (see the Details
# of ?fit.ihmm.scpg). Its inflection time should be reported
# as not identified; a model with fewer time-dependent parameters may be
# compared by AIC.

fit_ihmm2_m3$inflection

Description

Estimates the parameters of an IHSMM with class ihsmm.gamma, ihsmm.scpg and ihsmm.tscpg as typically returned by ihsmm.scpg using direct nonlinear optimization.

Usage

fit.ihsmm.scpg(
  ihsmmobj,
  method = "BFGS",
  maxit = 1000,
  reltol = sqrt(.Machine$double.eps),
  negative = TRUE
)

Arguments