drimmR : R Package for Estimation, Simulation and Reliability of drifting Markov Models

library(knitr)
knitr::opts_chunk$set(
    fig.align = "center",
    fig.height = 5.5,
    fig.width = 6,
    warning = FALSE,
    collapse = TRUE,
    dev.args = list(pointsize = 10),
    out.width = "65%"
)
# For reproductibility
set.seed(42)
library(drimmR)

Introduction

In this work we focus on multi-state systems modeled by means of a particular class of non-homogeneous Markov processes introduced in @Ver08, called drifting Markov processes. Most of the estimation methods and reliability results are developed in @BaVe2018.

Note that in many mathematical models it is assumed the homogeneity with respect to time, which is inappropriate in most of the applications. But, considering general non-homogeneous processes could be unrealistic from a practical point of view. For this reason, the drifting Markov chains assume that the Markov transition matrix is a linear/polynomial function of two/several Markov transition matrices. Thus we obtain a "smooth" non-homogeneity, with sensibly less parameters than in the general case.

Few R packages have been developed to handle similar types of models of Markov, semi-Markov or hidden semi-Markov type, useful in reliability or DNA analysis. For semi-Markov models we have semiMarkov R package @krol_semimarkov_2015 that performs maximum likelihood estimation for parametric continuous-time semi-Markov processes, smm R package @BaBeCeSaVe which performs parametric and non-parametric estimation and simulation for multi-state discrete-time semi-Markov processes. Two R packages are also dedicated to hidden semi-Markov models, implementing estimation and prediction methods: the hsmm R package @bulla_hsmm_2010 and the mhsmm R package @oconnell_hidden_2011.

Note that there is no R package developed for drifting Markov models (DMMs). Thus the purpose of this paper is to present an R package that we have developed, called drimmR, which performs estimation and simulation for such models, as well as the estimation of associated reliability measures. The aim of this paper is to describe the different possibilities of this package. To summarize, the package drimmR that we present deals with different problems:

We would like to mention that a web interface called WebDRIMM has been also developed [cf. @web2018] for simulating and estimating drifting Markov models, as well as associated reliability indicators; it is available at http://bioinfo.univ-rouen.fr/WebDRIMM\ The paper is organized as follows. Section 2{reference-type="ref" reference="section_DMM"} describes the drifting Markov models used in this package, present associated reliability indicators and corresponding estimation results and techniques. Section 3{reference-type="ref" reference="section_package"} illustrates the different functions of the drimmR package and we end the paper by presenting some concluding remarks on this R package in Section 4{reference-type="ref" reference="concluding"}.

Drifting Markov models, estimation and associated reliability measures {#section_DMM}

Let us consider a random system with finite state space $E={1,\ldots, s}$, $s < \infty.$ We assume that the time evolution of a system is governed by a discrete-time stochastic process with values in $E.$ In the following we will introduce a class of processes according to @Ver08 in Section 2.1{reference-type="ref" reference="subsection_DMM"}, we will briefly present associated reliability indicators according to @BaVe2018 in Section 2.2{reference-type="ref" reference="subsection_rel"} and estimation of the parameters in Section 2.3{reference-type="ref" reference="subsection_estimation"}, according to @Ver08 @BaVe2018.

Drifting Markov models {#subsection_DMM}

Let $\boldsymbol{\Pi}0 = (\Pi_0(u, v)){u,v \in E}$ and $\boldsymbol{\Pi}1 = (\Pi_1(u, v)){u,v \in E}$ be two Markov transition matrices of order 1 over the state space $E.$

::: {#def_DriftLin .definition} Definition 1 (linear drifting Markov chain of order $1$ and of length $n$). A sequence $X_0,$ $X_1, \ldots, X_n$ with state space $E=\mathbb {1, 2, \ldots, s}$ is said to be a linear drifting Markov chain (of order 1) of length $n$ between the Markov transition matrices $\boldsymbol{\Pi}0$ and $\boldsymbol{\Pi}_1$ if the distribution of $X_t,$ $t = 1, \ldots, n,$ is defined by $$\label{eq_DriftLin1} \mathbb P(X_t=v \mid X{t-1} = u, X_{t-2}, \ldots ) = \Pi_{\frac{t}{n}}(u, v), \; u, v \in E,$$ where $$\label{eq_DriftLin2} \Pi_{\frac{t}{n}}(u, v) = \left( 1 - \frac{t}{n} \right) \Pi_0(u, v) + \frac{t}{n} \Pi_1(u, v), \; u, v \in E.$$ :::

Let us denote by $\boldsymbol{\alpha} = (\alpha(1), \ldots, \alpha(s))$ the initial distribution of the chain, that is the distribution of $X_0,$ $\alpha(u)=\mathbb P (X_0 = u)$ for any state $u\in E.$

The linear drifting Markov model of order $1$ can be generalized to polynomial drifting Markov model of order $k$ and degree $d$. Let $\boldsymbol{\Pi}{\frac{i}{d}} = (\boldsymbol{\Pi}{\frac{i}{d}}(u_1, \dots, u_k, v))_{u_1, \dots, u_k,v \in E}$ be $d$ Markov transition matrices (of order $k$) over a state space $E.$

::: {#def_DriftPol .definition} Definition 2 (polynomial drifting Markov chain of order $k$ and of length $n$). A sequence $X_0,$ $X_1,$ $\ldots,$ $X_n$ with state space $E=\mathbb {1, 2, \ldots, s}$ is said to be a polynomial drifting Markov chain of order k and of length $n$ if the distribution of $X_t,$ $t = 1, \ldots, n,$ is defined by $$\label{eq_DriftPol1} \mathbb P(X_t=v \mid X_{t-1} = u_k, X_{t-2}=u_{k-1}, \ldots ) = \boldsymbol{\Pi}{\frac{t}{n}}(u_1, \dots, u_k, v), \; u_1, \dots, u_k, v \in E,$$ where $$\label{eq_DriftPol2} \boldsymbol{\Pi}{\frac{t}{n}}(u_1, \dots, u_k, v) = \displaystyle\sum_{i=0}^d A_i(t)\boldsymbol{\Pi}{\frac{i}{d}}(u_1, \dots, u_k,v), \; u_1, \dots, u_k, v \in E,$$ with $A_i$ polynomials of degree $d$ such as, for any $i, j \in {0, 1, \ldots, d},$ $A_i\left(\frac{nj}{d}\right)=\mathbb{1}{{i=j}}.$ :::

We would like to stress that the coherence between notations implies the choice of the notation $\boldsymbol{\Pi}_{\frac{i}{d}}$ and that, in fact, $A_i$ are Lagrange polynomials; see @Ver08 for more details on these two points.

Reliability of drifting Markov models {#subsection_rel}

In order to undertake a reliability analysis of a system modeled by a DMM, let us assume that the state space of the system is partitioned into working and failure states, $E=U \cup D,$ with $U \cap D=\emptyset,$ where $U={1,\ldots,s_1}$ represents the working states and $D={s_1+1,\ldots,s}$ the failure states of the system. According to this partition of the state space we partition any matrix of vector we are working with and we denote the corresponding partitions accordingly (e.g., $\boldsymbol{\Pi}_0^{UU},$ $\boldsymbol{\Pi}_0^{DU},$ $\boldsymbol{\alpha}^{U}$ etc.).

::: proposition []{#prop_RelDrift label="prop_RelDrift"} For a linear drifting Markov chain of order $1$ $(X_t){0 \leq t \leq n},$ the reliability at time $l,$ $l \in \mathbb N,$ is given by $$\label{eq_RelDrift} R(l) = \boldsymbol{\alpha}^{U} \; \prod{t=1}^l \left( \left(1 - \frac{t}{n}\right) \boldsymbol{\Pi}0^{UU} + \frac{t}{n} \boldsymbol{\Pi}_1^{UU} \right) \; \mathbb{1}^U, \textrm{ where } \mathbb{1}^U = (\underbrace{1, \cdots, 1}{s_1})^{\top}.$$ :::

::: proposition []{#prop_AvDrift label="prop_AvDrift"} For a linear drifting Markov chain of order $1$ $(X_t){0 \leq t \leq n},$ the pointwise (or instantaneous) availability at time $l,$ $l \in \mathbb N,$ is given by $$\label{eq_AvDrift} A(l) = \boldsymbol{\alpha} \; \prod{t=1}^l \left( \left(1 - \frac{t}{n}\right) \boldsymbol{\Pi}0 + \frac{t}{n} \boldsymbol{\Pi}_1 \right) \; \mathbb{1}^{E,U}, \textrm{ where } \mathbb{1}^{E,U} = (\underbrace{1, \cdots, 1}{s_1}, \underbrace{0, \cdots, 0}_{s-s_1} )^{\top}.$$ :::

::: proposition []{#prop_MaintDrift label="prop_MaintDrift"} For a linear drifting Markov chain of order $1$ $(X_t){0 \leq t \leq n},$ the maintainability at time $l,$ $l \in \mathbb N,$ is given by $$\label{eq_MaintDrift_gen} M(l) = 1 - \boldsymbol{\alpha}^{D} \; \prod{t=1}^l \left( \left(1 - \frac{t}{n}\right) \boldsymbol{\Pi}0^{DD} + \frac{t}{n} \boldsymbol{\Pi}_1^{DD} \right) \; \mathbb{1}^D, \textrm{ where } \mathbb{1}^D = (\underbrace{1, \cdots, 1}{s-s_1})^{\top}.$$ :::

::: proposition []{#prop_BMPDrift label="prop_BMPDrift"} For a linear drifting Markov chain of order $1$ $(X_t){0 \leq t \leq n},$ the BMP-failure rate at time $l at time$l,$$l N,$is given by % \begin{equation} \label{eq_BMPDrift} \lambda(l) \left{ \begin{array}{ll} 1- \frac{\mu_0^U \ \prod{t=1}^{l}( \ (1-\frac{t}{n}) \pi_0^{UU} + (\frac{t}{n}) \pi_1^{UU}) \ \mathbb{1}^U}{\mu_0^U \ \prod_{t=1}^{l-1}( \ (1-\frac{t}{n}) \pi_0^{UU} + (\frac{t}{n}) \pi_1^{UU}) \ \mathbb{1}^U} \ \ , \ \text{si R(l-1) != 0 }\ 0 \ , \ otherwise \ \end{array} \right. \end{equation}

\end{proposition}

\begin{proposition}\label{prop_RGDrift} For a linear drifting Markov chain of order$1$$(X_t)_0t n,$ the RG-failure rate at time$l at time $l,$ $l \in \mathbb N,$ is given by $$\label{eq_RGDrift} r(l) \left{ \begin{array}{ll} -\ln \frac{\mu_0^U \ \prod_{t=1}^{l}( \ (1-\frac{t}{n}) \pi_0^{UU} + (\frac{t}{n}) \pi_1^{UU}) \ \mathbb{1}^U}{\mu_0^U \ \prod_{t=1}^{l-1}( \ (1-\frac{t}{n}) \pi_0^{UU} + (\frac{t}{n}) \pi_1^{UU}) \ \mathbb{1}^U} \ \ , \ if \ l \ge 1 \ ,\ -\ln R(0) \ , \ if \ \ l = 0 \ \end{array} \right.$$ :::

See @BaVe2018 for more details and @Bar2004b; @Bar2008b for similar questions in a discrete-time semi-Markov framework.

Estimation of drifting Markov models {#subsection_estimation}

In this section we will consider different types of data for which the estimators of the characteristics of a drifting Markov chain and of the associated reliability indicators will be derived.

One can observe one sample path, that will be denoted by ${\mathcal H}(m,n):= (X_0,X_1, \ldots,X_{m}),$ where $m$ denotes the length of the sample path and $n$ the length of the drifting Markov chain. Two cases can be considered: (a1) $m=n$ (a complete sample path); (a2) $m < n$ (an incomplete sample path).

One can also observe $H$ i.i.d. sample paths, ${\mathcal H}_i(m_i,n_i), i=1, \ldots, H.$ Four cases can be considered here: (b1) $m_i=n_i=n$ for all $i=1, \ldots, H$ (complete sample paths of drifting Markov chains of the same length) ; (b2) $n_i=n$ for all $i=1, \ldots, H$ (incomplete sample paths of drifting Markov chains of the same length); (b3) $m_i=n_i$ for all $i=1, \ldots, H$ (complete sample paths of drifting Markov chains of different lengths) ; (b4) $m_i \leq n_i$ for all $i=1, \ldots, H$ (incomplete sample paths of drifting Markov chains of different lengths).

We have developed cf. @Ver08; @BaVe2018 mean square estimators starting from data under these frameworks; we present here only the case of a linear drifting Markov chain of order $1$ under the sample framework (b1) (complete sample paths of drifting Markov chains of the same length).

::: proposition []{#prop_B1 label="prop_B1"} Under the setting (b1), starting with $H$ complete sample paths of drifting Markov chains of the same length of a linear drifting Markov chain between two Markov transition matrices (of order 1) $\boldsymbol{\Pi}_0$ and $\boldsymbol{\Pi}_1,$ for any states $u, v \in E,$ the estimators of $\Pi_0(u,v)$ and $\Pi_1(u,v)$ are given by:

$$\begin{aligned} &&\widehat{\Pi}{0; (n,H)} (u,v) = \frac{P_1(H, \boldsymbol{m}, \boldsymbol{n}) P_2(H, \boldsymbol{m}, \boldsymbol{n}) - P_3(H, \boldsymbol{m}, \boldsymbol{n})P_4(H, \boldsymbol{m}, \boldsymbol{n})} {P_5(H, \boldsymbol{m}, \boldsymbol{n}) P_1(H, \boldsymbol{m}, \boldsymbol{n}) - P_3(H, \boldsymbol{m}, \boldsymbol{n})^2}\ && \widehat{\Pi}{1; (n,H)} (u,v)=\frac{P_5(H, \boldsymbol{m}, \boldsymbol{n}) P_4(H, \boldsymbol{m}, \boldsymbol{n}) - P_3(H, \boldsymbol{m}, \boldsymbol{n}) P_2(H, \boldsymbol{m}, \boldsymbol{n})} {P_5(H, \boldsymbol{m}, \boldsymbol{n})P_1(H, \boldsymbol{m}, \boldsymbol{n}) - P_3(H, \boldsymbol{m}, \boldsymbol{n})^2},\end{aligned}$$

where we have introduced the following notation:\

::: small $\displaystyle P_1(H, \boldsymbol{m}, \boldsymbol{n}) = \sum_{t=1}^{m}\sum_{h=1}^H \mathbb{1}{{X{t-1}^h=u }} \left(\frac{t}{n}\right)^2,$ $\displaystyle P_2(H, \boldsymbol{m}, \boldsymbol{n}) = \sum_{t=1}^{m} \sum_{h=1}^H\mathbb{1}{{X{t-1}^h=u, X_t^h=v }} \left(1-\frac{t}{n}\right),$\ $\displaystyle P_3(H, \boldsymbol{m}, \boldsymbol{n}) = \sum_{t=1}^{m}\sum_{h=1}^H\mathbb{1}{{X{t-1}^h=u }} \left(1-\frac{t}{n}\right)\left(\frac{t}{n}\right),$ $\displaystyle P_4(H, \boldsymbol{m}, \boldsymbol{n}) = \sum_{t=1}^{m}\sum_{h=1}^H\mathbb{1}{{X{t-1}^h=u, X_t^h=v }}\left(\frac{t}{n}\right),$\ $\displaystyle P_5(H, \boldsymbol{m}, \boldsymbol{n}) = \displaystyle\sum_{t=1}^{m}\sum_{h=1}^H\mathbb{1}{{X{t-1}^h=u }}\left(1-\frac{t}{n}\right)^2,$ with $\boldsymbol{m}=(m_1, \ldots, m_H)$ and $\boldsymbol{n}=(n_1, \ldots, n_H).$ ::: :::

Note we can adapt the estimation procedures that we have previously obtained in order to get estimators of the drifting Markov models in the other cases; one can see @BaVe2018 for more details. Using the expression of the reliability indicators of a drifting Markov chain previously obtained and the estimators of the characteristics of a drifting Markov chain, one immediately obtains the associated plug-in estimators of the reliability metrics.

The package {#section_package}

The package is principally devoted to the simulation and estimation of drifting Markov models, as well as to the estimation of associated reliability measures and to the computation of the probabilities of a word occurrence under a given model. All the different possibilities of the package are illustrated in Figure 2{reference-type="ref" reference="fig"}.

Estimation of drifting Markov models {#subsection_pack_estim}

The estimation of DMMs is carried out by the functions fitdmm, when starting from one (several) sample path (and obtain LSE). We will describe this function in the sequel.

  1. The function fitdmm

The different arguments of this function are:

Here we have an example of an estimation of a drifting Markov model of order 1 and degree 1 using the function fitdmm, starting from a DNA sequence called lambda.

data(lambda, package = "drimmR")
states <- c("a","c","g","t")
dmm <- fitdmm(lambda, 1, 1, states, init.estim="freq", fit.method="sum")
dmm

Estimation of corresponding model characteristics {#subsection_pack_charact}

Once a Markov drifting model has been estimated/constructed as described in the previous subsection, various characteristics of the model can be computed. These are: the log-likelihood of one or several sequences, the AIC and BIC information criteria, the stationary distribution on an entire sequence or only on a part of it, the distribution of the chain on the entire sequence or only on a part of it.

The functions loglik, aic and bic have the following arguments:

They return the numerical values of the log-likelihood, the AIC, the BIC of the sequence, respectively.

data(lambda, package = "drimmR")
sequence <- c("a","g","g","t","c","g","a","t","a","a","a")
dmm <-fitdmm(lambda, 1, 1, c('a','c','g','t'), init.estim = "freq")

loglik(dmm, sequence)

aic(dmm, sequence)

bic(dmm, sequence)

The function getTransitionMatrix evaluates the transition matrix at a given position. The function has the following arguments :

data(lambda)
dmm <- fitdmm(lambda, 1, 1, c('a','c','g','t'),init.estim = "freq")
t <- 10
getTransitionMatrix(dmm,pos=t)

The function getStationaryLaw evaluates the stationary law at a given position or for all positions along the list of sequence(s).

data(lambda, package = "drimmR")
sequence <- sample(lambda, 30, replace=TRUE )
dmm <- fitdmm(sequence, 1, 1, c('a','c','g','t'), init.estim = "freq")
t <- 10
getStationaryLaw(dmm,pos=t)
getStationaryLaw(dmm,all.pos=TRUE)

Analysis of a model {#subsection_pack_analysis}

Once a Markov drifting model has been estimated/constructed as described at the beginning of this section, the obtained model can be analyzed by means of several functions.\

  1. The computation of word probabilities according to an estimated/constructed DMM is carried out by means of the functions word_probability, word_probabilities and words_probabilities.

The function word_probability computes the probability of a word at a given position; it has the following arguments:

It returns the numerical value of the probability.\

data(lambda, package = "drimmR")
dmm <- fitdmm(lambda, 1, 1, c('a','c','g','t'), init.estim = "freq")
word_probability("cgt",10,dmm,output_file=tempfile("drimmR"))

The function word_probabilities computes the probabilities of a word at given positions. It has the following arguments:

It returns a numerical vector of probabilities computed at the positions given in pos.\

data(lambda, package = "drimmR")
dmm <- fitdmm(lambda, 1, 1, c('a','c','g','t'), init.estim = "freq", fit.method="sum")
res <- word_probabilities("cgt", c(100, length(lambda)-3), dmm, output_file=tempfile("drimmR"), plot=TRUE)
head(res[[1]], n=10)
res[2]

The function words_probabilities computes the probabilities of a word at given positions; it has the following arguments:

It returns a data frame of probabilities computed for each element of words at the positions given in pos.\

data(lambda, package = "drimmR")
dmm <- fitdmm(lambda, 1, 1, c('a','c','g','t'), init.estim = "freq")
res <- words_probabilities(c("atcgattc", "taggct", "ggatcgg"), c(100, length(lambda)-8), dmm, output_file=tempfile("drimmR"), plot=TRUE)
head(res[[1]], n=10)
res[2]
  1. The function lengthWord_probabilities computes the probabilities of occurrence of the observed word of given size in a sequence at several positions. This function has the following arguments:

  2. n: Integer, the given length of the word

  3. sequence: Vector of characters representing the sequence

  4. pos: Vector of integer positions of the word along the sequence. It it is given by a start and end point

  5. x: An object of class dmm

  6. output_file A file containing the vector of probabilities

  7. plot display several figure plots of probabilities of appearance of words along the positions if TRUE

It returns a data frame of probability by position\

data(lambda, package = "drimmR")
dmm <- fitdmm(lambda, 1, 1, c('a','c','g','t'), init.estim = "freq")
res <- lengthWord_probabilities(2, lambda, c(100, length(lambda)-2), dmm, output_file=tempfile("drimmR"), plot=TRUE)
res[2]

Simulation of drifting Markov models {#subsection_pack_simul}

The simulation of drifting Markov models is carried out by the function simulation that has the following arguments:

data(lambda, package = "drimmR")
dmm <- fitdmm(lambda, 1, 1, c('a','c','g','t'), init.estim = "freq")
# simulate a sequence of length 20 000 from dmm
simulate(dmm, tempfile("drimmR"), 20000)

Reliability of drifting Markov models {#subsection_pack_reliab}

The reliability measures are computed by means of the following functions: A (availability), R (reliability,survival function), M (maintainability), errorRate (the classical failure rate, called BMP-failure rate and a more recent failure rate adapted to discrete data, called RG-failure rate). For more details on these reliability measures one can see @Bar2004b; @Bar2008b. All these functions have the following arguments:

The function errorRate enables to select the type of evaluated failure rate with an additional argument :

These functions return vectors with the values of the corresponding reliability measure.

Estimation of reliability {#subsection_pack_reliab}

data(lambda, package = "drimmR")
dmm <- fitdmm(lambda, 1, 1, c("a","c","g","t"), init.estim ="freq", fit.method="sum")
reliability(dmm, k1=1, k2=10, upstates=c("a","c"),output_file = tempfile("drimmR"), plot=TRUE)

Estimation of availability {#subsection_pack_avail}

data(lambda, package = "drimmR")
dmm <- fitdmm(lambda, 1, 1, c("a", "c", "g", "t"), init.estim = "freq", fit.method="sum")
availability(dmm, k1=1, k2=10, upstates=c("a","c"), output_file = tempfile("drimmR"), plot=TRUE)

Estimation of maintainability {#subsection_pack_maintain}

data(lambda, package = "drimmR")
dmm <- fitdmm(lambda, 1, 1, c("a","c","g","t"), init.estim="freq", fit.method="sum")
maintainability(dmm, k1=1, k2=10, upstates=c("a","c"), output_file = tempfile("drimmR"), plot=TRUE)

Estimation of failure rates {#subsection_pack_reliab}

data(lambda, package = "drimmR")
dmm <- fitdmm(lambda, 1, 1, c("a","c","g","t"), init.estim="freq", fit.method="sum")
failureRate(dmm, k1=1, k2=10, upstates=c("a","c"), failure.rate="BMP", output_file=tempfile("drimmR"), plot=TRUE)

Concluding remarks {#concluding}

To conclude, in this paper we have presented drimmR, an R package for simulation, estimation and reliability and survival analysis of drifting Markov models. These are versatile stochastic models of Markov type capable of taking into account a time non-homogeneity of a known, controlled shape. For this reason these models can represent interesting modeling alternatives to classical models (like Markov models, semi-Markov models, etc.) and can be useful for researchers, practitioners and engineers in various fields.

The package addresses several items important from practical point of view, when carrying out the estimation: we consider one or several samples, the sample paths can be complete or not, can come from models of the same length or of different lengths.

The fields of application of the DMMs can be numerous; we only want to point out three of them.

Flowchart of the drimmR package.{#fig width="70%"}



Try the drimmR package in your browser

Any scripts or data that you put into this service are public.

drimmR documentation built on Oct. 9, 2026, 5:08 p.m.