R: Estimation of a drifting semi-Markov chain

fit_dsmm {dsmmR}

R Documentation

Estimation of a drifting semi-Markov chain

Description

Estimation of a drifting semi-Markov chain, given one sequence of states. This estimation can be parametric or non-parametric and is available for the three types of drifting semi-Markov models.

Usage

fit_dsmm(
  sequence,
  degree,
  f_is_drifting,
  p_is_drifting,
  states = NULL,
  initial_dist = "unif",
  estimation = "nonparametric",
  f_dist = NULL
)

Arguments

`sequence`	Character vector that represents a sequence of states from the state space `E`.
`degree`	Positive integer that represents the polynomial degree `d` for the drifting semi-Markov model.
`f_is_drifting`	Logical. Specifies if `f` is drifting or not.
`p_is_drifting`	Logical. Specifies if `p` is drifting or not.
`states`	Character vector that represents the state space `E`, with length equal to `s = \|E\|`. Default value is set equal to the sorted, unique states present in the given sequence.
`initial_dist`	Optional. Character that represents the method to estimate the initial distribution. `"unif"` : The initial distribution of each state is equal to `1/s` (default value). `"freq"` : The initial distribution of each state is equal to the frequency that it has in the sequence.
`estimation`	Optional. Character. Represents whether the estimation will be nonparametric or parametric. `"nonparametric"` : The estimation will be non-parametric (default value). `"parametric"` : The estimation will be parametric.
`f_dist`	Optional. It can be defined in two ways: If `estimation = "nonparametric"`, it is equal to `NULL` (default value). If `estimation = "parametric"`, it is a character array that specifies the distributions of the sojourn times, for every state transition. The list of possible values is: `["unif", "geom", "pois", "dweibull", "nbinom", NA]`. It can be defined in two ways: If `f` is not drifting, it has dimensions of `s \times s`. If `f` is drifting, it has dimensions of `s \times s \times (d+1)` (see more in Details, Parametric Estimation.) It is defined similarly to the attribute `f_dist` in dsmm_parametric.

Details

This function estimates a drifting semi-Markov model in the parametric and non-parametric case. The parametric estimation can be achieved by the following steps:

We obtain the non-parametric estimation of the sojourn time distributions.
We estimate the parameters for the distributions defined in the attribute f_dist through the probabilities that were obtained in step 1.

Three different models are possible for to be estimated for each case. A normalization technique is used in order to correct estimation errors from small sequences. We will use the exponentials (1), (2), (3) to distinguish between the drifting semi-Markov kernel \widehat{q}_{\frac{t}{n}} and the semi-Markov kernels \widehat{q}_\frac{i}{d} used in Model 1, Model 2, Model 3. More about the theory of drifting semi-Markov models in dsmmR.

Non-parametric Estimation

Model 1

When the transition matrix p of the embedded Markov chain (J_{t})_{t\in \{0,\dots,n\}} and the conditional sojourn time distribution f are both drifting, the drifting semi-Markov kernel can be estimated as:

\widehat{q}_{\frac{t}{n}}^{\ (1)}(u,v,l) = \sum_{i = 0}^{d}A_{i}(t)\ \widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l),

The semi-Markov kernels \widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l), i = 0, \dots, d, are estimated through Least Squares Estimation (LSE) and are obtained after solving the following system, \forall t \in \{0, \dots, n\}, \forall u, v \in E and \forall l \in \{1, \dots, k_{max}\}:

MJ = P,

where the matrices are written as:

M = (M_{ij})_{i,j \in \{0, \dots, d\} } = \left(\sum_{t=1}^{n}1_{u}(t)A_{i}(t)A_{j}(t)\right)_{ i,j \in \{0, \dots, d\}}
J = (J_i)_{i \in \{0, \dots, d\} } = \left(\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)\right)_{ i \in \{0, \dots, d\}}
P=(P_i)_{i\in \{0, \dots, d\} }= \left(\sum_{t=1}^{n}1_{uvl}(t)A_{i}(t)\right)_{ i \in \{0, \dots, d\}}

and we use the following indicator functions:

1_{u}(t) = 1_{ \{J_{t-1} = u \} } = 1, if at t the previous state is u, 0 otherwise.
1_{uvl}(t) = 1_{ \{J_{t-1} = u, J_{t} = v, X_{t} = l \} } = 1, if at t the previous state is u with sojourn time l and next state v, 0 otherwise.

In order to obtain the estimations of \widehat{p}_{\frac{i}{d}}(u,v) and \widehat{f}_{\frac{i}{d}}(u,v,l), we use the following formulas:

\widehat{p}_{\frac{i}{d}}(u,v) = \sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l),

\widehat{f}_{\frac{i}{d}}(u,v,l) = \frac{\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{ \sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.

Model 2

In this case, p is drifting and f is not drifting. Therefore, the estimated drifting semi-Markov kernel will be given by:

\widehat{q}_{\frac{t}{n}}^{\ (2)}(u,v,l) = \sum_{i=0}^{d} A_{i}(t)\ \widehat{q}_{\frac{i}{d}}^{\ (2)}(u,v,l),

\forall t \in \{0,\dots,n\}, \forall u,v\in E, \forall l\in \{1,\dots, k_{max} \}, where k_{max} is the maximum sojourn time that was observed in the sequence and A_i, i = 0, \dots, d are d + 1 polynomials with degree d (see dsmmR). In order to obtain the estimators \widehat{p} and \widehat{f}, we use the estimated semi-Markov kernels \widehat{q}_{\frac{i}{d}}^{\ (1)} from Model 1. Since p is drifting, we define the estimation \widehat{p} the same way as we did in Model 1. In total, we have the following estimations, \forall u,v \in E, \forall l \in \{1,\dots, k_{max} \}:

\widehat{p}_{\frac{i}{d}}(u,v) = \sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l),

\widehat{f}(u,v,l) = \frac{\sum_{i=0}^{d}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{ \sum_{i=0}^{d}\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.

Thus, the estimated semi-Markov kernels for Model 2, \widehat{q}_{\frac{i}{d}}^{\ (2)}(u,v,l) = \widehat{p}_{\frac{i}{d}}(u,v)\widehat{f}(u,v,l), can be written with regards to the estimated semi-Markov kernels of Model 1, \widehat{q}_{\frac{i}{d}}^{\ (1)}, as in the following:

\widehat{q}_{\frac{i}{d}}^{\ (2)}(u,v,l) = \frac{ (\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)) (\sum_{i = 0}^{d}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l))}{ \sum_{i = 0}^{d}\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.

Model 3

In this case, f is drifting and p is not drifting. Therefore, the estimated drifting semi-Markov kernel will be given by:

\widehat{q}_{\frac{t}{n}}^{\ (3)}(u,v,l) = \sum_{i=0}^{d} A_{i}(t)\ \widehat{q}_{\frac{i}{d}}^{\ (3)}(u,v,l),

\forall t \in \{0,\dots,n\}, \forall u,v\in E, \forall l\in \{1,\dots, k_{max} \}, where k_{max} is the maximum sojourn time that was observed in the sequence and A_i, i = 0, \dots, d are d + 1 polynomials with degree d (see dsmmR). In order to obtain the estimators \widehat{p} and \widehat{f}, we use the estimated semi-Markov kernels estimated semi-Markov kernels \widehat{q}_{\frac{i}{d}}^{\ (1)} from Model 1. Since f is drifting, we define the estimation \widehat{f} the same way as we did in Model 1. In total, we have the following estimations, \forall u,v \in E, \forall l \in \{1,\dots, k_{max} \}:

\widehat{p}\ (u,v) = \frac{\sum_{i=0}^{d}\sum_{l=1}^{k_{max}} \widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{d+1},

\widehat{f}_{\frac{i}{d}}(u,v,l) = \frac{\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{ \sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.

Thus, the estimated semi-Markov kernels for Model 3, \widehat{q}_{\frac{i}{d}}^{\ (3)}(u,v,l) = \widehat{p}\ (u,v)\widehat{f}_{\frac{i}{d}}(u,v,l), can be written with regards to the estimated semi-Markov kernels of Model 1, \widehat{q}_{\frac{i}{d}}^{\ (1)}, as in the following:

\widehat{q}_{\frac{i}{d}}^{\ (3)}(u,v,l) = \frac{ \widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l) \sum_{i=0}^{d}\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)} {(d+1)\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.

Parametric Estimation

In this package, the parametric estimation of the sojourn time distributions defined in the attribute f_dist is achieved as follows:

We obtain the non-parametric LSE of the sojourn time distributions f.
We estimate the parameters for the distributions defined in f_dist through the probabilities of f, estimated in previously in 1.

The available distributions for the modelling of the conditional sojourn times of the drifting semi-Markov model, defined from the argument f_dist, have their parameters estimated through the following formulas:

Geometric (p):

f(x) = p (1-p)^{x-1}, where x = 1, 2, \dots,k_{max}. We estimate the probability of success \widehat{p} as such:

\widehat{p} = \frac{1}{E(X)}
Poisson (\lambda):

f(x) = \frac{\lambda^{x-1} exp(-\lambda)}{(x-1)!}, where x = 1, 2,\dots,k_{max}. We estimate \widehat{\lambda} > 0 as such:

\widehat{\lambda} = E(X)
Negative binomial (\alpha, p):

f(x)=\frac{\Gamma(x + \alpha - 1)}{\Gamma(\alpha)(x-1)!} p^{\alpha}(1-p)^{x-1}, where x = 1, 2,\ldots,k_{max}. \Gamma is the Gamma function, p is the probability of success and \alpha \in (0, +\infty) is the parameter describing the number of successful trials, or the dispersion parameter (the shape parameter of the gamma mixing distribution). We estimate them as such:

\widehat{p} = \frac{E(X)}{Var(X)},

\widehat{\alpha} = E(X)\frac{\widehat{p}}{1 - \widehat{p}} = \frac{E(X)^2}{Var(X) - E(X)}.
Discrete Weibull of type 1 (q, \beta):

f(x)=q^{(x-1)^{\beta}}-q^{x^{\beta}}, where x= 1, 2, \dots,k_{max}, q is the first parameter with 0 < q < 1 and \beta \in (0, +\infty) the second parameter. We estimate them as such:

\widehat{q} = 1 - f(1),

\widehat{\beta} = \frac{\sum_{i = 2}^{k_{max}} \log_{i}(\log_{\widehat{q}}(\sum_{j = 1}^{i}f(j)))}{k_{max} - 1}.

Note that we require k_{max} \geq 2 for estimating \widehat{\beta}.
Uniform (n): f(x) = 1/n where x = 1, 2, \dots, n, for n a positive integer. We use a numerical method to obtain an estimator for \widehat{n} in this case.

Value

Returns an object of S3 class (dsmm_fit_nonparametric, dsmm) or (dsmm_fit_parametric, dsmm). It has the following attributes:

dist : List. Contains 2 or 3 arrays.
- If estimation = "nonparametric" we have 2 arrays:
  - p_drift or p_notdrift, corresponding to whether the defined p transition matrix is drifting or not.
  - f_drift or f_notdrift, corresponding to whether the defined f sojourn time distribution is drifting or not.
- If estimation = "parametric" we have 3 arrays:
  - p_drift or p_notdrift, corresponding to whether the defined p transition matrix is drifting or not.
  - f_drift_parametric or f_notdrift_parametric, corresponding to whether the defined f sojourn time distribution is drifting or not.
  - f_drift_parameters or f_notdrift_parameters, which are the defined f sojourn time distribution parameters, depending on whether f is drifting or not.
emc : Character vector that contains the embedded Markov chain (J_{t})_{t\in \{0,\dots,n\}} of the original sequence. It is this attribute of the object that describes the size of the model n. Last state is also included, for a total length of n+1, but it is not used for any calculation.
soj_times : Numerical vector that contains the sojourn times spent for each state in emc before the jump to the next state. Last state is also included, for a total length of n+1, but it is not used for any calculation.
initial_dist : Numerical vector that contains an estimation for the initial distribution of the realized states in sequence. It always has values between 0 and 1.
states : Character vector. Passing down from the arguments. It contains the realized states given in the argument sequence.
s : Positive integer that contains the length of the character vector given in the attribute states, which is equal to s = |E|.
degree : Positive integer. Passing down from the arguments. It contains the polynomial degree d considered for the drifting of the model.
k_max : Numerical value that contains the maximum sojourn time, which is the maximum value in soj_times, excluding the last state.
model_size : Positive integer that contains the size of the drifting semi-Markov model n, which is equal to the length of the embedded Markov chain (J_{t})_{t\in \{0,\dots,n\}}, minus the last state. It has a value of length(emc) - 1, for emc as defined above.
f_is_drifting : Logical. Passing down from the arguments. Specifies if f is drifting or not.
p_is_drifting : Logical. Passing down from the arguments. Specifies if p is drifting or not.
Model : Character. Possible values:
- "Model_1" : Both p and f are drifting.
- "Model_2" : p is drifting and f is not drifting.
- "Model_3" : f is drifting and p is not drifting.
estimation : Character. Specifies whether parametric or nonparametric estimation was used.
A_i : Numerical Matrix. Represents the polynomials A_i(t) with degree d that were used for solving the system MJ = P. Used for the methods defined for the object. Not printed when viewing the object.
J_i : Numerical Array. Represents the estimated semi-Markov kernels of the first model (\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l))_{i\in\{0,\dots,d\}} that were obtained after solving the system MJ = P. Not printed when viewing the object.

References

V. S. Barbu, N. Limnios. (2008). semi-Markov Chains and Hidden semi-Markov Models Toward Applications - Their Use in Reliability and DNA Analysis. New York: Lecture Notes in Statistics, vol. 191, Springer.

Vergne, N. (2008). Drifting Markov models with Polynomial Drift and Applications to DNA Sequences. Statistical Applications in Genetics Molecular Biology 7 (1).

Barbu V. S., Vergne, N. (2019). Reliability and survival analysis for Drifting Markov models: modelling and estimation. Methodology and Computing in Applied Probability, 21(4), 1407-1429.

T. Nakagawa and S. Osaki. (1975). The discrete Weibull distribution. IEEE Transactions on Reliability, R-24, 300-301.

Examples

# Create a random sequence
sequence <- create_sequence("DNA", len = 2000, seed = 1)
## Alternatively, we could obtain a sequence as follows:
## > data("lambda", package = "dsmmR")
## > sequence <- c(lambda)
states <- sort(unique(sequence))
degree <- 3


# ===========================================================================
# Nonparametric Estimation.
# Fitting a random sequence under distributions of unknown shape.
# ===========================================================================


# ---------------------------------------------------------------------------
# Both p and f are drifting - Model 1.
# ---------------------------------------------------------------------------

obj_model_1 <- fit_dsmm(sequence = sequence,
                        degree = degree,
                        f_is_drifting = TRUE,
                        p_is_drifting = TRUE,
                        states = states,
                        initial_dist = "freq",
                        estimation = "nonparametric", # default value
                        f_dist = NULL # default value
                        )

cat(paste0("We fitted a sequence with ", obj_model_1$Model, ",\n",
           "model size: n = ", obj_model_1$model_size, ",\n",
           "length of state space: s = ", obj_model_1$s, ",\n",
           "maximum sojourn time: k_max = ", obj_model_1$k_max, " and\n",
           "polynomial (drifting) Degree: d = ", obj_model_1$degree, ".\n"))

# Get the drifting p and f arrays.
p_drift <- obj_model_1$dist$p_drift
f_drift <- obj_model_1$dist$f_drift

cat(paste0("Dimension of p_drift: (s, s, d + 1) = (",
           paste(dim(p_drift), collapse = ", "), ").\n",
           "Dimension of f_drift: (s, s, k_max, d + 1) = (",
           paste(dim(f_drift), collapse = ", "), ").\n"))

# We can even check the embedded Markov chain and the sojourn times
# directly from the returned object, if we wish to do so.
# This is achieved through the `base::rle()` function, used on `sequence`.
model_emc <- obj_model_1$emc
model_sojourn_times <- obj_model_1$soj_times

# ---------------------------------------------------------------------------
# Fitting the sequence when p is drifting and f is not drifting - Model 2.
# ---------------------------------------------------------------------------


obj_model_2 <- fit_dsmm(sequence = sequence,
                        degree = degree,
                        f_is_drifting = FALSE,
                        p_is_drifting = TRUE)


cat(paste0("We fitted a sequence with ", obj_model_2$Model, ".\n"))


# Get the drifting p and non-drifting f arrays.
p_drift_2 <- obj_model_2$dist$p_drift
f_notdrift <- obj_model_2$dist$f_notdrift


all.equal.numeric(p_drift, p_drift_2) # p is the same as in Model 1.


cat(paste0("Dimension of f_notdrift: (s, s, k_max) = (",
           paste(dim(f_notdrift), collapse = ", "), ").\n"))



# ---------------------------------------------------------------------------
# Fitting the sequence when f is drifting and p is not drifting - Model 3.
# ---------------------------------------------------------------------------

obj_model_3 <- fit_dsmm(sequence = sequence,
                        degree = degree,
                        f_is_drifting = TRUE,
                        p_is_drifting = FALSE)

cat(paste0("We fitted a sequence with ", obj_model_3$Model, ".\n"))
# Get the drifting f and non-drifting p arrays.
p_notdrift <- obj_model_3$dist$p_notdrift
f_drift_3 <- obj_model_3$dist$f_drift
all.equal.numeric(f_drift, f_drift_3) # f is the same as in Model 1.
cat(paste0("Dimension of f_notdrift: (s, s) = (",
           paste(dim(p_notdrift), collapse = ", "), ").\n"))

# ===========================================================================
# Parametric Estimation
# Fitting a random sequence under distributions of known shape.
# ===========================================================================
### Comments
### 1.  For the parametric estimation it is recommended to use a common set
###     of distributions while only the parameters (of the sojourn times)
###     are drifting. This results in (generally) higher accuracy.
### 2.  This process is similar to that used in `dsmm_parametric()`.


s <- length(states)
# Getting the distributions for the states.
# Rows correspond to previous state `u`.
# Columns correspond to next state `v`.
f_dist_1 <- matrix(c(NA,         "unif",     "dweibull", "nbinom",
                     "pois",      NA,        "pois",     "dweibull",
                     "geom",     "pois",      NA,        "geom",
                     "dweibull", 'geom',     "pois",      NA),
                   nrow = s, ncol = s, byrow = TRUE)
f_dist <- array(f_dist_1, dim = c(s, s, degree + 1))
dim(f_dist)

# ---------------------------------------------------------------------------
# Both p and f are drifting - Model 1.
# ---------------------------------------------------------------------------

obj_fit_parametric <- fit_dsmm(sequence = sequence,
                               degree = degree,
                               f_is_drifting = TRUE,
                               p_is_drifting = TRUE,
                               states = states,
                               initial_dist = 'unif',
                               estimation = 'parametric',
                               f_dist = f_dist)
cat("The class of `obj_fit_parametric` is : (",
    paste0(class(obj_fit_parametric), collapse = ', '), ").\n")

# Estimated parameters.
f_params <- obj_fit_parametric$dist$f_drift_parameters

# The drifting sojourn time distribution parameters.
f_0 <- f_params[,,,1]
f_1.3 <- f_params[,,,2]
f_2.3 <- f_params[,,,3]
f_1 <- f_params[,,,4]

params <- paste0('q = ', round(f_params["c", "t", 1, ], 3),
                 ', beta = ', round(f_params["c", "t", 2, ], 3))
f_names <- c("f_0", paste0("f_", 1:(degree-1), "/", degree), "f_1")
all_names <- paste(f_names, ":", params)
cat("The drifting of the parameters for passing from \n",
    "`u` = 'c' to `v` = 't' under a discrete Weibull distribution is:",
    "\n", all_names[1], "\n", all_names[2],
    "\n", all_names[3], "\n", all_names[4])

# ---------------------------------------------------------------------------
# f is not drifting, only p is drifting - Model 2.
# ---------------------------------------------------------------------------

obj_fit_parametric_2 <- fit_dsmm(sequence = sequence,
                                 degree = degree,
                                 f_is_drifting = FALSE,
                                 p_is_drifting = TRUE,
                                 initial_dist = 'unif',
                                 estimation = 'parametric',
                                 f_dist = f_dist_1)

cat("The class of `obj_fit_parametric_2` is : (",
    paste0(class(obj_fit_parametric_2), collapse = ', '), ").\n")
# Estimated parameters.
f_params_2 <- obj_fit_parametric_2$dist$f_notdrift_parameters

params_2 <- paste0('q = ', round(f_params_2["c", "t", 1], 3),
                   ', beta = ', round(f_params_2["c", "t", 2], 3))

cat("Not-drifting parameters for passing from ",
    "`u` = 'c' to `v` = 't' \n under a discrete Weibull distribution are:\n",
    paste("f :", params_2))


# ===========================================================================
# `simulate()` and `get_kernel()` can be used for the two objects,
# `dsmm_fit_nonparametric` and `dsmm_fit_parametric`.
# ===========================================================================


sim_seq_nonparametric <- simulate(obj_model_1, nsim = 10)
str(sim_seq_nonparametric)


kernel_drift_parametric <- get_kernel(obj_fit_parametric, klim = 10)
str(kernel_drift_parametric)

[Package dsmmR version 1.0.5 Index]