We introduce the BMRMM package implementing Bayesian inference for a class of Markov renewal mixed models which can characterize the stochastic dynamics of a collection of sequences, each comprising alternative instances of categorical states and associated continuous duration times, while being influenced by a set of exogenous factors as well as a ‘random’ individual. The default setting flexibly models the state transition probabilities using mixtures of Dirichlet distributions and the duration times using mixtures of gamma kernels while also allowing variable selection for both. Modeling such data using simpler Markov mixed models also remains an option, either by ignoring the duration times altogether or by replacing them with instances of an additional category obtained by discretizing them by a user-specified unit. The option is also useful when data on duration times may not be available in the first place. We demonstrate the package’s utility using two data sets.
Markov models (MMs) are widely used for modeling the transition dynamics of categorical state sequences. Classical Markov renewal models (MRMs) additionally model the state duration times, when available, where the state transitions follow Markov dynamics and the state duration times follow a continuous distribution that depends on the immediately preceding and following states (Figure 1). M(R)Ms have been widely used in different variations (Phelan 1990; Eichelsbacher and Ganesh 2002; Muliere et al. 2003; Alvarez 2005; Diaconis and Rolles 2006; Bulla and Muliere 2007; Etterson et al. 2007; Bacallado et al. 2009; Li 2009; Epifani et al. 2014; Siebert and Söding 2016; Holsclaw et al. 2017; Sesia et al. 2019). There are also some sparse works on covariate-dependent Markov models (Muenz and Rubinstein 1985; Gradner 1990; Alioum et al. 1998; Islam and Chowdhury 2006).
The existing literature however focuses very heavily on modeling single sequences. (Sarkar et al. 2018) developed a highly flexible and computationally efficient class of Bayesian Markov mixed models (BMMMs) for jointly modeling a collection of categorical sequences, each one associated with an individual as well as a set of time-invariant external covariates (e.g., the sex and genotype of the associated individual, an experimental condition under which the sequence was generated, etc.). BMMMs characterize the state transition probabilities using a convex combination of a fixed covariate-dependent component and a random individual-level component, both being Dirichlet distributed. They further allow covariate levels with similar effects to be probabilistically clustered together, allowing automatic selection of the significant covariates, and providing a sophisticated framework for analyzing data sets having the aforementioned structure.
BMMMs however do not model duration times of the states which are often additionally available in real-world applications. Recently, (Wu et al. 2023) extended BMMMs to Bayesian Markov renewal mixed models (BMRMMs), allowing for the additional analysis of continuous duration times which, depending on the application, can either be the duration for which a state persists or the duration between two consecutive states, i.e., inter-state intervals (ISIs). Specifically, they modeled the duration times using mixtures of gamma kernels with mixture probabilities being a convex combination a covariate-dependent effect and an individual-level effect, similar to BMMMs. Covariate levels with similar influences on mixture probabilities are clustered together in BMRMMs as well, allowing the selection of the significant covariates. BMRMMs thus holistically model both state transitions and continuous duration times, painting a comprehensive picture of the underlying stochastic dynamics.
In this article, we describe the R package BMRMM which implements BMMMs and BMRMMs, collectively referred to henceforth as BM(R)MMs. The BMRMM package runs posterior inference for categorical state transitions and continuous duration times, if available, via a Markov chain Monte Carlo (MCMC) algorithm, returning an object containing comprehensive inference results. The package also includes a suite of plotting functions to display the results graphically. Specifically for continuous duration times, when available, the package provides users with three different options: (i) ignore the duration times and model the state transitions alone as a BMMM; (ii) introduce an additional category by discretizing the continuous duration times by a user-specified unit, and analyze the appended state transitions as a BMMM; (iii) model the duration as a mixture of gamma kernels using a BMRMM, as proposed in Wu et al. (2023). Additionally, users can choose to turn off one or both of the fixed covariate effects and the random individual effects. Overall, the BMRMM package thus gives users a lot of flexibility in handling their data sets, providing inferences for both Bayesian Markov renewal or non-renewal models as needed.
The BMRMM package conveniently includes a synthetic foxp2
data set
that describes the laboratory study on the role of the FoxP2 gene
implicated in speech deficiencies for adult mice, which is also the
motivating application for the methodology of BMMMs and BMRMMs.
Mutations in the FoxP2 gene have long been associated with severe
deficits in vocal communication for mammals
(MacDermot et al. 2005). Mice with and without the mutation
singing under various "social contexts" have thus been studied in many
experiments
(Fujita et al. 2008; Castellucci et al. 2016; Chabout et al. 2016; Gaub et al. 2016).
[The FoxP2 data set (Chabout et al. 2016), e.g., comprises the sequences
of syllables making up the songs as well as the lengths of
inter-syllable intervals (ISIs). The data set foxp2
included in the
BMRMM package is taken from the simulation study of
(Wu et al. 2023). It is much shorter than the real FoxP2 data set but
closely mimics its other aspects and is used] in
this paper to demonstrate how to obtain detailed inferences for both
syllable transitions and ISI dynamics for a comprehensive analysis of
the vocal repertoire in mice with and without the FoxP2 mutation.
The utility of the BMRMM package goes well beyond the FoxP2 data set. As described above, the package is designed for scenarios where the data set consists of categorical state sequences associated with an individual as well as a number of observed covariates where additional data on continuous duration times may or may not be available. Data sets with such structures are frequently observed in different areas of scientific research and can potentially benefit from the BMRMM package. For example, Islam and Chowdhury (2006) analyzed the transitions of different rainfall orders in three districts of Bangladesh under three covariates, wind speed, humidity, and maximum temperature. In an education assessment study, Zhang et al. (2019) recorded sequences of writing states, characterized by keystroke logs, for 257 eighth graders of various genders, races, and socioeconomic statuses. Combescure et al. (2003) estimated the control states of 371 asthma patients with different body mass indices (BMIs) and disease severity over a four-year-long study and produced the asthma control data set, which we will use as an additional example to demonstrate the utility of our package in this paper. For such data sets, the BMRMM package provides a flexible, sophisticated, and principled way to model fixed effects of the covariates and random effects of the individuals in both the state transition dynamics and the distribution of the ISIs.
Other computer programs for Markov models with covariates include MARKOV (Marshall et al. 1995) and MKVPCI (Alioum and Commenges 2001). R packages for analyzing discrete Markov models include markovchain (Spedicato 2017) and msm (Jackson 2011). SemiMarkov (Listwon and Saint-Pierre 2015) and SMM (Barbu et al. 2018) provide functions for the simulation and estimation of traditional semi-Markov models. Some R packages provide the inference of hidden semi-Markov models, such as mhsmm (O’Connell and Højsgaard 2011) and hhsmm (Amini et al. 2022). Ferguson et al. (2012) built the msSurv package which provides a nonparametric estimation of semi-Markov models but does not consider covariates. There are also R packages implementing MRPs in specific application areas. For example, Kharrat et al. (2019) introduced the Countr package for flexible regression models based on MRPs. Pustejovsky (2021) developed the ARPobservation for simulating behavior streams based on alternating renewal processes. Other R packages for categorical data analysis include catdap (Katsura 1980) and vcd (Meyer et al. 2022). To our knowledge, there has not been an R package that implements flexible Bayesian M(R)MMs.
[In the following section, we summarize the technical details of BMRMMs. Documentation for the functions of the BMRMM package is then provided. Next, we demonstrate the usage of our package in analyzing two different data sets. The final section contains some concluding remarks.]
We briefly describe the BM(R)MM methodologies here – more details can be found in Sarkar et al. (2018) and Wu et al. (2023). Consider specifically a sequence s comprising Ts state instances and let ys,t denote the state at time t in sequence s. The states ys,t’s come from a set Y={1,2,…,d0}. Within a sequence s, there are Ts−1 duration times (state persistence times or inter-state intervals), denoted by {τs,t}s0,Tss=1,t=2, where τs,t is the duration between the (t−1)th and tth states in sequence s, and s0 represents the total number of sequences. Figure 1 presents a graphical summary of the data structure. Each sequence s is associated with p categorical covariates or factors, denoted by xs,j∈Xj={1,2,…,dj}, and an individual, denoted by is. Without loss of generality, we assume that the analyses of the transition probabilities and the duration times distributions both include all p covariates. Moreover, the analysis of duration times counts the previous state ys,t−1 as an additional (p+1)th covariate. In the BMRMM package, users have the flexibility to select particular covariates for each analysis and exclude the previous state from the analysis of duration times. In their original definition in (Pyke 1961), the duration time τs,t in an MRP was allowed to depend on both the preceding state ys,t−1 and the following state ys,t. To keep the notation simple and the methodology easy to understand for a broad audience, we however only include the preceding state ys,t−1 as a predictor of τs,t in this paper. This analysis can be easily modified to have the pair (ys,t−1,ys,t) as a predictor instead of just ys,t−1, as was actually done in (Wu et al. 2023).
Figure 1: Graphical model showing the data structure: ys,t denotes the observed state at the tth time location in the sth sequence; τs,t denotes the observed duration times (either state persistence times or inter-state intervals) between the states ys,t−1 and ys,t; each sequence s is also associated with an individual is and a set of exogenous time-invariant covariates xs,1,…,xs,p. The Markov mixed model considered in this article analyzes the state transitions ys,t in a collection of sequences; the Markov renewal mixed model additionally analyzes the duration times τs,t; both models accommodate fixed effects of the covariates xs,1,…,xs,p and random effects of the individuals is.
For a sequence s associated with individual i and covariate levels
x1,…,xp, the transition probabilities
Pr(ys,t=yt∣is=i,xs,1=x1,…,xs,p=xp,ys,t−1=yt−1)=P(i)trans,x1,…,xp(yt∣yt−1)
are defined as a convex combination of a fixed covariate effect
component λtrans,x1,…,xp(⋅∣yt−1) and a
random effect component λ(i)trans:
P(i)trans,x1,…,xp(yt∣yt−1)=π(i)trans,0(yt−1)λtrans,x1,…,xp(yt∣yt−1)+π(i)trans,1(yt−1)λ(i)trans(yt∣yt−1).
For each covariate j=1,…,p, it is possible that some covariate levels exert a similar effect on the transition dynamics. [For example, the components λtrans,x1=1,x2,…,xp(⋅∣yt−1) and λtrans,x1=2,x2,…,xp(⋅∣yt−1) would be equal if levels 1 and 2 of covariate 1 have similar influences on transition dynamics for fixed levels for covariates 2,…,p.] A clustering mechanism for covariate levels allows the fixed component λtrans,x1,…,xp(⋅∣yt−1) to be the same for all levels with a similar influence. In particular, for covariate j, we construct the partition C(j)trans={C(j)trans,hj}ktrans,jhj=1 of its levels, where ktrans,j is the number of clusters for covariate j and hj represents the cluster index. We introduce latent variables {ztrans,j,ℓ}p,djj=1,ℓ=1 that indicate the cluster index for the ℓth label of covariate j. [Two levels of the covariate j, ℓ1,ℓ2∈Xj={1,…,dj}, are clustered together if and only if ztrans,j,ℓ1=ztrans,j,ℓ2.] For the fixed effects, we then replace the covariate levels x1,…,xp’s with cluster indices h1,…,hp’s and present the fixed effect as λtrans,h1,…,hp(⋅∣yt−1).
We set Dirichlet priors for both fixed and individual effect components
and let them center around the same mean vector λtrans,0
to facilitate posterior computation. The probability vector
λtrans,0 is also given a Dirichlet prior with mean
λtrans,00 to capture the natural preferences of certain
states in Y:
λtrans,h1,…,hp(⋅∣yt−1)∼Dir{αtrans,0λtrans,0(1∣yt−1),…,αtrans,0λtrans,0(d0∣yt−1)},λ(i)trans(⋅∣yt−1)∼Dir{α(0)transλtrans,0(1∣yt−1),…,α(0)transλtrans,0(d0∣yt−1)},λtrans,0(⋅∣yt−1)∼Dir{αtrans,00λtrans,00(1),…,αtrans,00λtrans,00(d0)}.
We present the complete Bayesian hierarchical model for the transition
dynamics as
(ys,t∣ys,t−1=yt−1,is=i,ztrans,1,xs,1=h1,…,ztrans,p,xs,p=hp)∼Mult{P(i)trans,h1,…,hp(1∣yt−1),…,P(i)trans,h1,…,hp(d0∣yt−1)}, whereP(i)trans,h1,…,hp(⋅∣yt−1)=π(i)trans,0(yt−1)λtrans,h1,…,hp(⋅∣yt−1)+π(i)trans,1(yt−1)λ(i)trans(⋅∣yt−1),ztrans,j,ℓ∼Mult{μtrans,j(1),…,μtrans,j(dj)}, μtrans,j∼Dir(αtrans,j,…,αtrans,j),λtrans,h1,…,hp(⋅∣yt−1)∼Dir{αtrans,0λtrans,0(1∣yt−1),…,αtrans,0λtrans,0(d0∣yt−1)},λ(i)trans(⋅∣yt−1)∼Dir{α(0)transλtrans,0(1∣yt−1),…,α(0)transλtrans,0(d0∣yt−1)},λtrans,0(⋅∣yt−1)∼Dir{αtrans,00λtrans,00(1),…,αtrans,00λtrans,00(d0)},π(i)trans,0(yt−1)∼Beta(atrans,0,atrans,1),αtrans,0∼Ga(atrans,0,btrans,0), α(0)trans∼Ga(a(0)trans,b(0)trans).
The BMRMM package provides three options for analyzing the duration times: (i) ignore the durations altogether and only model the transition probabilities of the existing states, (ii) treat the durations as blocks of a new special category, with a discretization unit specified by users, (iii) model the durations as a continuous random variable with a flexible mixture of gamma distributions. For the first two options, we only need to apply the model described in the previous subsection. For the third option, we need to conduct a separate analysis of the duration times as described below.
Let K denote the number of gamma mixture components in the model for
the duration times. We let
P(i)dur(⋅∣x1,…,xp,ys,t−1) denote the
mixture probability vector given individual i, covariate levels
x1,…,xp, and preceding syllable ys,t−1. The preceding state
ys,t−1 can be removed from the formula if the users do not wish to
consider its influence in the inference of the duration times. The
distribution of the continuous duration times,
{τs,t}s0,Tss=1,t=2, is then modeled as
f(τs,t∣is=i,xs,1=x1,…,xs,p=xp,ys,t−1=yt−1)=K∑k=1P(i)dur(k∣x1,…,xp,yt−1)Ga(τs,t∣αk,βk),
f(τs,t∣zdur,s,t=k)∼Ga(τs,t∣αk,βk),Pr(zdur,s,t=k∣is=i,xs,1=x1,…,xs,p=xp,ys,t−1=yt−1)=P(i)dur(k∣x1,…,xp,yt−1).
Similar to the model for the transition probabilities, the mixture
probabilities are a convex combination of a fixed population-level
effect and a random individual-level effect:
P(i)dur(⋅∣x1,…,xp,yt−1)=π(i)dur,0(⋅)λdur,x1,…,xp,yt−1(⋅)+π(i)dur,1(⋅)λ(i)dur(⋅),
The mixture probability vectors are given Dirichlet priors with the mean
vector λdur,0, which itself centers around a global
vector λdur,00:
λdur,g1,…,gp+1(⋅)∼Dir{αdur,0λdur,0(1),…,αdur,0λdur,0(K)},λ(i)dur(⋅)∼Dir{α(0)durλdur,0(1),…,α(0)durλdur,0(K)},λdur,0(⋅)∼Dir{αdur,00λdur,00(1),…,αdur,00λdur,00(K)}.
We present the complete Bayesian hierarchical model for the continuous
duration times as
(τs,t∣zdur,s,t=k)∼Ga(τs,t∣αk,βk),(zdur,s,t∣is=i,zdur,1,xs,1=g1,…,zdur,p,xs,p=gp,zdur,p+1,ys,t−1=gp+1)∼Mult{P(i)dur,g1,…,gp+1(1),…,P(i)dur,g1,…,gp+1(K)}, whereP(i)dur,g1,…,gp+1(k)=π(i)dur,0(k)λdur,g1,…,gp+1(k)+π(i)dur,1(k)λ(i)dur(k),λdur,g1,…,gp+1(⋅)∼Dir{αdur,0λdur,0(1),…,αdur,0λdur,0(K)}, αdur,0∼Ga(adur,0,bdur,0),λ(i)dur(⋅)∼Dir{α(0)durλdur,0(1),…,α(0)durλdur,0(K)}, α(0)dur∼Ga(a(0)dur,b(0)dur),λdur,0(⋅)∼Dir{αdur,00λdur,00(1),…,αdur,00λdur,00(K)},π(i)dur,0(k)∼Beta(adur,0,adur,1),αk∼Ga(adur,0,bdur,0), βk∼Ga(adur,0,bdur,0).
Inference is based on samples drawn from the posterior using a [Metropolis-Hastings-within-Gibbs MCMC algorithm. Most full conditionals are available in closed form and can be directly sampled from. A Metropolis-Hastings step is however used for updating the discrete valued cluster configurations.] There is, however, no conjugate prior for gamma distributions with unknown shape parameters (Damsleth 1975). Recently, Miller (2019) designed a procedure that efficiently approximates the posterior full conditionals of gamma shape parameters under a gamma prior with another gamma density. We adopt this approximation in our MCMC algorithm.
The BMRMM package is developed to implement Bayesian Markov
(renewal) mixed models. The main function BMRMM
of the package carries
out detailed analyses of the state transitions and their duration times
(if applicable) as described in the previous section. [Moreover, the
package includes a number of supplementary functions that use the
results of the main function to produce numerical summaries,
visualizations, and diagnostics. Table 1 provides a brief
description of all functions. ]
Function | Description |
---|---|
BMRMM | Creates a BMRMM object. |
summary.BMRMM | Summary for an object of class BMRMM and create a BMRMMsummary object. |
plot.BMRMMsummary | Visualization of a BMRMMsummary object. |
hist.BMRMM | Returns histograms of duration times for a BMRMM object. |
diag.BMRMM | Provides MCMC diagnostic plots for a BMRMM object. |
model.selection.scores | Returns the LPML and WAIC scores of the mixture gamma model. |
The main function is BMRMM
which implements the inference for both the
state transition probabilities and the duration times. We summarize the
parameters in Table 3 and present the function as
follows.
BMRMM(data, num.cov, cov.labels = NULL, state.labels = NULL,
random.effect = TRUE, fixed.effect = TRUE,
trans.cov.index = 1:num.cov, duration.cov.index = 1:num.cov,
duration.distr = NULL, duration.incl.prev.state = TRUE,
simsize = 10000, burnin = simsize/2)
[The parameter data
specifies the target data set and needs to follow
a certain structure. ] The first column should list
the individual IDs is, followed by p columns for the values of
the p associated covariates xs,j, then two columns for the values
of the previous state ys,t−1, the current state ys,t, and
finally a column for duration times τs,t. The package supports
one to five categorical covariates that take on values 1,2,….
The duration times column is optional if the user would like to use BMMM
instead of BMRMM to analyze just the state transitions. This is shown in
Table 2. The users can look at the included
[simulated data set] foxp2
as an example.
Id Covariate 1 … Covariate p Previous State Current State State Durations/ISI |
: Table 2: Columns of the desired input data set.
[The number of covariates in the data set is specified by the argument
num.cov
. The argument cov.labels
is a list of vectors giving the
names of covariate levels in the covariate order that is presented in
data
while the parameter state.labels
is a vector providing the
names of the transition states. The default labels are Arabic numerals.] The [random.effect
]
parameter gives users the option to exclude the random individual
effects. If [random.effect
] is set to FALSE
,
the transition probabilities (and the mixture probabilities for duration
times, if applicable) will only consider the influence of the covariate
levels. Similarly, the [fixed.effect
] parameter
allows users to exclude the fixed population effects. [The default
values for random.effect
and fixed.effect
are both
TRUE
.] The covariate indices for the two analyses
can be specified by setting trans.cov.index
and duration.cov.index
.
[We note that indices specified by trans.cov.index
and
duration.cov.index
refer to the index of the covariate when the first
covariate is given index 1, thus different from its index in
data
.]
[Users can define duration.distr
in the following three
ways.]
[If users set duration.distr
to be NULL
, which is the default
setting, then the duration times will be ignored and not modeled at
all. ] The BMMM described will be implemented
to analyze the existing state transitions alone.
[If duration.distr
is set as list(‘mixDirichlet’, unit)
, the
duration times will be used to construct a new state ‘dur.state’
,
which will be analyzed along with the original set of states.] The additional argument unit
must be
defined and acts both as a threshold and as a block size for
duration times. For example, if the unit
is set to 5, then for
each duration value greater than 5 units, each block of 5 unit
in it will be treated as an instance of a new ’dur.state’
state.
[If there is a state transition from state ‘a’
to ‘b’
with a
duration time of 15 seconds and the unit
is specified at 5
seconds, then the updated Markov sequence will contain three
consecutive ‘dur.state’
states, i.e.,
(‘a’, ‘dur.state’, ‘dur.state’, ‘dur.state’, ‘b’)
.] Since we adopt the floor operation, a
duration time of say 17.68 seconds will also be replaced by three
consecutive instances of ’dur.state’
states in this example. The
BMMM model will then be implemented to analyze the resulting
appended state transitions.
These first two options may naturally result in loss of information and is therefore not recommended when a detailed analysis of the distribution of the duration times is warranted.
[If duration.distr
is set to be list(‘mixgamma’, shape, rate)
,
the duration times are modeled as a continuous random variable using
a flexible mixture of gamma kernels, as described for a BMRMM
model.] In this case, users can specify the
prior shape and rate parameters with the shape
and rate
arguments within the definition of duration.distr
. We note that
shape
and rate
must be numeric vectors of the same length.
By default, we consider the previous state ys,t−1 as a covariate
when we model the duration times as continuous variables, i.e.,
[duration.incl.prev.state
] is set to TRUE
.
Users can set this parameter to FALSE
if they wish to exclude the
previous state when analyzing the duration times. [The remaining
parameters simsize
and burnin
denote the total number of MCMC
iterations and the number of burn-ins, respectively.]
Argument | Explanation | Default value |
---|---|---|
data |
the data set to be used following the required format | |
num.cov |
an integer giving the number of observed covariates in data |
|
[cov.labels ] a list of vectors giv |
ing names of all covariate levels [NULL ] |
|
[state.labels ] a vector giving names |
of the states [NULL ] |
|
[random.effect ] TRUE if random indi |
vidual effects are included [TRUE ] |
|
[fixed.effect ] TRUE if fixed popul |
ation effects are included [TRUE ] |
|
trans.cov.index |
selects the covariates to analyze for transition probabilities | 1:num.cov |
duration.cov.index |
selects the covariates to analyze for duration times | 1:num.cov |
duration.distr |
specifies the distribution for duration times | NULL |
[duration.incl.prev.state ] TRUE if yt−1 a |
cts as a covariate for the analysis of duration times [TRUE ] |
|
simsize |
number of MCMC iterations | 10000 |
burnin |
number of burnins of the MCMC algorithm | simsize /2 |
The BMRMM
function returns [an object of class BMRMM
, which either
contains only results.trans
or both of results.trans
and
results.duration
if duration times follow a mixture gamma
distribution.] For the state transitions, the
posterior mean transition probability matrices for each combination of
the covariate levels and each individual are given by
results.trans$tp.exgns.post.mean
and
results.trans$tp.anmls.post.mean
, respectively. Additionally,
results.trans$clusters
stores cluster configurations for each
covariate from each MCMC iteration. As for duration times, the fields
results.duration$shape.samples
and results.duration$rate.samples
record the shape and rate parameters, for each mixture component in
every MCMC iteration, respectively. Meanwhile,
results.duration$comp.assignment
gives the assignment of the mixture
component for each data point in the last MCMC iteration. Similar to
transition probabilities, results.duration$clusters
gives the cluster
configurations of the covariates. Other elements of results.trans
and
results.duration
can be found in the detailed R function description.
[The BMRMM package provides an S3 method for summarizing results of
a BMRMM
object as follows. ]
summary.BMRMM(object, delta = 0.02, digits = 2, ...)
[The object
must be of class BMRMM
. The argument delta
is
associated with local tests for transition probabilities, which we will
explain further. The digit
parameter is an integer used for number
formatting, as in the general summary
function. The summary.BMRMM
function returns an object of class BMRMMsummary
with the following
fields. ]
trans.global
and dur.global
[These two fields give the global test results from the inference of transition probabilities and duration times. Global tests show the significance of the covariates in affecting the state transitions and duration times. ] Specifically, for each covariate, the empirical distribution of the size of the clusters in the stored MCMC iterations is calculated. The null hypothesis that a covariate is not important is equivalent to the event that all its levels are in the same cluster, or, in other words, that the cluster size for the covariate is just one.
trans.probs.mean
and trans.probs.sd
The two fields provide the mean and standard deviation for the posterior mean of each transition type under all combinations of covariate levels, respectively.
trans.local.mean.diff
and trans.local.null.test
[The BMRMMsummary
object also contains local test results for
transition probabilities. Local tests analyze the differences
between the transition probabilities associated with two different
levels of a covariate j, fixing the levels of the other
covariates.] For every pair of levels of
covariate j, trans.local.mean.diff
gives the absolute
differences in transition probabilities for each transition type in
the MCMC iterations. The local null hypothesis we test for each
transition type is that this difference is at least the
pre-specified value delta
. [Meanwhile, trans.local.null.test
gives the probability of the null hypothesis that the difference
between two covariate levels is not significant under each
transition type. ]
dur.mix.params
and dur.mix.probs
For each mixture component, dur.mix.params
provides the estimates
of the gamma shape and rate parameters from the last MCMC iteration.
For every covariate level, users can obtain the mixture
probabilities by calling the field dur.mix.probs
, which can be
further used to estimate the length of the duration times.
[The main plotting function of the package, plot.BMRMMsummary
, is an
S3 method for class BMRMMsummary
. It gives the barplots for global
tests as well as heatmaps for the posterior mean and standard deviation
for transition probabilities, local tests for transition probabilities,
mixture parameters and probabilities for duration times. The parameters
of plot.BMRMMsummary
include x
, which must be an object of class
BMRMMsummary
and type
, which is a single string representing the
field of x
that needs to be plotted. The function also takes general
plotting arguments such as xlab
, ylab
, etc. ]
plot.BMRMMsummary(x, type, xlab = NULL, ylab = NULL, main = NULL, col = NULL, ...)
[When duration times are analyzed as continuous variables using mixture
gamma distributions, the users can use the S3 method hist.BMRMM
to
generate histograms for duration times along with the estimated
posterior distribution. The parameter x
is an object of class BMRMM
.
The argument comp
gives the specific mixture component that the user
would like to investigate. When comp
is NULL
, which is the default
setting, the histogram of all observed duration times is plotted and
superimposed with the posterior mean of the fitted mixture gamma
distribution. When comp
is a specific integer, we will be looking at
the last MCMC iteration. The histogram for duration times assigned with
component comp
will be presented alongside the mixture gamma
distribution with the shape and rate parameters from the last MCMC
iteration. Users can refer to the documentation of the general hist
function to see the interpretation for the rest of the parameters.]
hist.BMRMM(x, comp = NULL, xlim = NULL, breaks = NULL, main = NULL,
col = 'gray', xlab = 'Duration times', ylab = 'Density', ...)
[Finally, users can check the MCMC diagnostics with the traceplots and
autocorrelation plots produced by the function diag.BMRMM
. The
object
parameter should be an object of class BMRMM
. For transition
probabilities, users can specify the covariate levels as well as the
state transitions they are interested in by defining cov.combs
and
transitions
, respectively. For duration times, users can define
components
, a numeric vector, to obtain the diagnostic plots for shape
and rate parameters of the specific component
kernels.]
diag.BMRMM(object, cov.combs = NULL, transitions = NULL, components = NULL)
When the duration times are modeled using mixtures of gamma
distributions, model selection can be performed on the number of mixture
components using the function model.selection.scores
.
model.selection.scores(object)
[The function takes an object
as its input, which must be an object of
class BMRMM
. It returns a list consisting of] the
log pseudo marginal likelihood (LPML) (Geisser and Eddy 1979) and the
widely applicable information criterion (WAIC) (Watanabe and Opper 2010)
scores of the model. Larger values of LPML and smaller values of WAIC
indicate better model fits. They are particularly suitable for complex
Bayesian hierarchical models as they can be easily computed from the
MCMC samples.
The FoxP2 data set records the songs sung by adult male mice of two
genotypes, wild type or FoxP2, denoted by W and F, respectively
(Chabout et al. 2016). The mice sang under three social contexts, U
(fresh female urine on a cotton tip placed inside the male’s cage), L
(an awake and behaving adult female placed inside the cage), and A (an
anesthetized female placed on the lid of the cage). Each song comprises
a sequence of syllables and continuous inter-syllable intervals (ISIs).
The data set can be used to analyze the effect of the FoxP2 gene on the
vocal syntax of mice, in turn providing insights into the effects of the
gene on human vocal communication abilities and related deficiencies.
[The real FoxP2 data set originates from the study by
(Chabout et al. 2016) and requires permission to use. (Wu et al. 2023)
simulated a data set that closely mimics the real one. For demonstrating
the BMRMM package, we included in it a shortened version of this
synthetic data set which we refer to as the foxp2
data set.] The foxp2
synthetic data set has 17391 rows
and 6 columns, which are Id, Genotype, Context, Prev_State, Cur_State,
and Transformed_ISI. The original FoxP2 data set records ISIs in
seconds. In the simulated data set foxp2
, following Wu et al. (2023),
log(1+ISI) values are used which give a better model fit.
Id | Genotype | Context | Prev_State | Cur_State | Transformed_ISI |
---|---|---|---|---|---|
1 | 2 | 2 | 3 | 3 | 0.20197711 |
1 | 2 | 2 | 3 | 3 | 0.06972753 |
1 | 2 | 2 | 3 | 3 | 0.07211320 |
1 | 2 | 2 | 3 | 3 | 0.15790932 |
1 | 2 | 2 | 3 | 3 | 0.06781471 |
1 | 2 | 2 | 3 | 3 | 0.09426236 |
If we are only interested in analyzing the transition probabilities with the covariates genotype and social contexts, we would use the main function as follows.
R> res.fp2 <- BMRMM(foxp2, num.cov = 2)
If we would like to pick specific covariates for our analyses, we can
define trans.cov.index
and duration.cov.index
accordingly. For
example, if we only want to use context for transition probabilities and
genotype for ISIs, we would run the following.
R> res.fp2 <- BMRMM(foxp2, num.cov = 2,
trans.cov.index = c(2), duration.cov.index = c(1))
If we would like to analyze the ISIs as part of the original state
sequence following a mixture Dirichlet distribution, as was done by
Sarkar et al. (2018), the ISIs are replaced by (possibly consecutive)
"silent" states by dividing them into blocks of 250 milliseconds. The
BMRMM
function can do this by setting duration.distr
as a list with
the string ’mixDirichlet’
and the argument unit
as
log(0.25+1), based on the log transformation.
R> res.fp2 <- BMRMM(foxp2, num.cov = 2,
duration.distr = list('mixDirichlet', unit = log(0.25+1)))
In the next example, we would like to analyze the ISIs as continuous variables following a mixture gamma distribution. For syllable transitions, we use both genotype and context as covariates. For the ISIs, in addition to these two, we also use the preceding syllable as a covariate.
R> res.fp2 <- BMRMM(data = foxp2, num.cov = 2, state.labels = c('d', 'm', 's', 'u'),
cov.labels = list(c('F', 'W'), c('U', 'L', 'A')),
duration.distr = list('mixgamma', shape = rep(1, 4), rate = rep(1, 4)))
In what follows, we show the results for the last function call. The
returned res.fp2
have two parts, which are named
res.fp2$results.trans
and res.fp2$results.duration
. Now we
demonstrate how we print and visualize the results.
[First, we obtain a BMRMMsummary
object, sm.fp2
, by calling the
summary.BMRMM
function on the returned results res.fp2
. The global
test results for identifying the significant covariates can be found by
calling the fields trans.global
and dur.global
. The function
plot.BMRMMsummary
is called to visualize the global tests using
barplots, as presented in Figure 2.] We recall that a covariate is significant when
its levels formed more than one cluster with very high posterior
probability (the bar heights). Figure 2 and the
printed results suggest that every covariate is significant for the ISIs
but only the social context is significant for the transition
probabilities.
R> sm.fp2 <- summary(res.fp2)
R> sm.fp2$trans.global
label_data
cluster_data Context Genotype
1 0 1
3 1 0
R> sm.fp2$dur.global
label_data
cluster_data Context Genotype prev_state
2 0.00 1.00 0.10
3 1.00 0.00 0.90
4 0.00 0.00 0.01
R> plot(sm.fp2, 'trans.global')
R> plot(sm.fp2, 'dur.global')
Figure 2: Results for the simulated foxp2 data set showing the global tests of significance of the covariates for the state transitions (left) and the ISIs (right). The bars represent the estimated posterior probabilities of the number of clusters formed by the levels of each covariate.
[The plotting function can be called to visualize the posterior transition probabilities under different combinations of the covariate levels. ] We show in Figure 3 the heatmaps for the posterior mean and standard deviation of the transition probabilities for each transition type for the following combinations of covariates: (F,A) and (W,L).
R> plot(sm.fp2, 'trans.probs.mean')
R> plot(sm.fp2, 'trans.probs.sd')
Figure 3: Results for the simulated foxp2 data set showing the posterior mean and standard deviation for each transition type for selected covariate combinations, (F,A) (top) and (W,L) (bottom).
We also perform the local test to assess the influence of genotype on
the transition probabilities by computing the absolute difference of the
transition probabilities between F and W among the thinned MCMC
samples after burn-ins, i.e.,
|Δλtrans,⋅,x2(yt∣yt−1)|=|λtrans,1,x2(yt∣yt−1)−λtrans,2,x2(yt∣yt−1)|.
The estimated posterior probability for the null hypothesis is therefore
the proportion of times
|Δλtrans,⋅,x2(yt∣yt−1)≤δ| is
observed in the MCMC samples, where x2 is the social context and
δ is the user-specific difference threshold delta
. [The
plotting function plot.BMRMMsummary
gives the plots for all local test
results if we set the type
to be ‘trans.local.mean.diff’
or
‘trans.local.null.test’
. Here, we show the results of local tests for
the covariate 1 (i.e., genotype) with delta
equaling the default value
of 0.02, and present the plots in
Figure 4.] From the figure, we
see that the posterior probabilities of the null hypotheses are
generally large for most transition types (e.g., transitions to the
syllable u) regardless of the social context, indicating that genotype
does not have a strong influence on transition probabilities with a
fixed context under these transition types.
R> plot(sm.fp2, 'trans.local.mean.diff')
R> plot(sm.fp2, 'trans.local.null.test')
Figure 4: Results for the simulated foxp2 data set showing local test results for genotypes fixing the social context, U (top), L (middle), and A (bottom). The averaged absolute difference in transition probabilities between F and W is presented on the left. The posterior probabilities of the corresponding null hypotheses are on the right.
Next, we turn our attention to the ISIs. We first check the fit of our estimated mixture gamma distribution presented in Figure 5. We then look further into the shape of each mixture component in Figure 5. From the histogram for each component, we see that components 2 and 4 represent longer ISIs while components 1 and 3 represent shorter ISIs.
R> hist(res.fp2, xlim = c(0,1))
R> for(comp in 1:4) {
hist(res.fp2, comp = comp)
}
Figure 5: Results for the simulated foxp2 data set showing the histograms of the ISIs with the estimated posterior gamma mixture density (left) and the histograms of the ISIs for each mixture component (right).
We examine the values of mixture parameters and mixture probabilities for each covariate level in the last MCMC iteration, which provides insights into the influence of the covariate on ISI lengths.
R> sm.fp2$dur.mix.params
shape.k rate.k
Comp 1 29.07 394.30
Comp 2 1.23 1.49
Comp 3 8.46 465.66
Comp 4 3.03 20.13
R> sm.fp2$dur.mix.probs
$Genotype
F W
Comp 1 0.46 0.48
Comp 2 0.19 0.13
Comp 3 0.10 0.15
Comp 4 0.25 0.24
$Context
U L A
Comp 1 0.58 0.35 0.48
Comp 2 0.14 0.16 0.18
Comp 3 0.08 0.21 0.08
Comp 4 0.20 0.28 0.25
$prev_state
d m s u
Comp 1 0.52 0.52 0.42 0.42
Comp 2 0.13 0.13 0.22 0.17
Comp 3 0.11 0.11 0.10 0.18
Comp 4 0.24 0.24 0.26 0.23
From the mixture probabilities, we see that mice with genotype F have a much higher mixture probability in component 2 compared to genotype W, which indicates mice with the FoxP2 mutation require a longer ISI between pronouncing two syllables, a reflection of vocal impairment.
Finally, we check the MCMC diagnostic plots and see if we had good mixing for the parameters. Here we focus on a specific transition type, u→m, for covariate combination {F,U} and a specific mixture component (component 2). We show these plots in Figure 6.
R> diag.BMRMM(res.fp2, cov.combs = list(c(1, 1)),
transitions = list(c(4, 2)), components = c(2))
Figure 6: Results for the simulated foxp2 data set showing the MCMC diagnostic plots, including traceplots and autocorrelation plots.
The BMRMM package is able to analyze duration times in detail which
could either be the ISIs, as seen in the synthetic foxp2
data set, or
the state persistence times, as in a traditional semi-Markov model. To
demonstrate the usage of our package in analyzing the state persistence
times, we use the asthma control data set from the ARIA (Association
pour la Recherche en Intelligence Artificielle) study of severe
asthmatic patients (Combescure et al. 2003) in France between 1997
and 2001. At each visit, a chest physician graded the asthma control
status of the patient using control scores (Juniper et al. 1999).
The data set contains the sojourn time of the control states as well as
three covariates: Asthma severity, sex, and the body mass index (BMI) of
the patients. Saint-Pierre et al. (2003) used a Markov model with piece-wise
constant intensities to model the asthma control evolution and proposed
a regression model for analyzing the effect of covariates.
Combescure et al. (2003) used the data set to assess the relationship
between asthma severity and control of asthma. Listwon and Saint-Pierre (2015)
fitted a semi-Markov model for the sojourn times using exponential and
Weibull distributions and analyzed the effect of covariates individually
due to complexity. Our BMRMM package is able to analyze the effect
of the three covariates while also incorporating random individual
effects exhibited by different patients on transition dynamics and state
duration times.
[The asthma
data set we use here is from the SemiMarkov package
(Listwon and Saint-Pierre 2015). We have renamed and reordered the columns such
that the data set fits the required format.]
Specifically, the data set has 928 rows, recording the asthma control
states of 371 patients, which is one of the following three transient
states: Optimal control (State 1), sub-optimal control (State 2), and
unacceptable control (State 3). Each state can transit to any other two
states and the state duration times are recorded. The data set also
contains three binary covariates of the asthma patients, including the
disease severity (1 if mild-moderate and 2 if severe), BMI (body mass
index, 1 if BMI <25 and 2 otherwise), and sex (1 if women and 2 if
men). We display part of the processed data in Table 5,
where Duration
is the sojourn time in Prev_State
.
Id | Severity | BMI | Sex | Prev_State | Cur_State | Duration |
---|---|---|---|---|---|---|
2 | 2 | 2 | 1 | 3 | 2 | 0.1533 |
2 | 2 | 2 | 1 | 2 | 2 | 4.1232 |
3 | 2 | 2 | 2 | 3 | 1 | 0.0958 |
3 | 2 | 2 | 2 | 1 | 3 | 0.2300 |
3 | 2 | 2 | 2 | 3 | 1 | 0.2656 |
3 | 2 | 2 | 2 | 1 | 1 | 5.4073 |
We investigate the transition dynamics and state persistence times of
the asthma
data set using the BMRMM
function. We consider K=4
mixture components for state persistence times. The choice of K is
derived from running the BMRMM model several times with different K’s
and comparing the fitness of the models using the LPML and WAIC scores.
R> res.asm <- BMRMM(data = asthma, num.cov = 3, state.labels = c(1, 2, 3),
cov.labels = list(c('Mild-Moderate','Severe'),
c('BMI<25','BMI>=25'),
c('Women','Men')),
duration.distr = list('mixgamma', shape = rep(1, 4),
rate = rep(1, 4)))
We name the returned BMRMM
object res.asm
and obtain the
BMRMMsummary
object sm.asm
by calling the summary.BMRMM
function.
As in the FoxP2 application, we first plot the global test results for
both transition probabilities and state persistence times in
Figure 7. For the transition probabilities,
only the severity of asthma is significant while for duration times only
the preceding state is significant. The BMI value and the sex are not
significant for either transition dynamics or state durations.
Figure 7: Results for the asthma data set showing the global tests of significance of the covariates for the state transitions (left) and the state persistence times (right). The bars represent the estimated posterior probabilities of the number of clusters formed by the levels of each covariate.
R> sm.asm <- summary(res.asm)
R> sm.asm$trans.global
label_data
cluster_data BMI Severity Sex
1 0.96 0.00 1.00
2 0.04 1.00 0.00
R> sm.asm$dur.global
label_data
cluster_data BMI prev_state Severity Sex
1 0.88 0.00 0.88 0.99
2 0.12 0.14 0.12 0.01
3 0.00 0.86 0.00 0.00
R> plot(sm.asm, 'trans.global')
R> plot(sm.asm, 'dur.global')
We show the posterior mean and standard deviations of the state transition probabilities for men and women with severe conditions and BMI ≥25 in Figure 8. We see that for severe patients with BMI ≥25, the transition probabilities are similar for men and women. We also take a look at the local test results for the BMI values fixing the severity of the patients in Figure 9. Though the absolute differences between the two covariate levels for BMI are small, the probabilities for the null hypotheses are also small, especially for transitions to state 1 and state 2. This suggests that even though the influence of BMI on state transitions is not significant globally, it is significant given the severity of asthma condition regardless of sex.
R> plot(sm.asm, 'trans.probs.mean')
R> plot(sm.asm, 'trans.probs.sd')
Figure 8: Results for the asthma data set showing the posterior mean and standard deviation for each transition type for selected covariate combinations, {Severe, BMI≥25,Men} (top) and {Severe, BMI≥25,Women} (bottom).
Figure 9: Results for the asthma data set showing local test results for BMI fixing asthma severity condition and sex, men (top) and women (bottom). The averaged absolute difference in transition probabilities between BMI <25 and BMI ≥25 is presented on the left. The posterior probability of the null hypothesis is on the right.
Figure 10 presents the histograms of the entire asthma data set, superimposed with the posterior mean of the mixture gamma distribution. With four components, the estimated mixture gamma fits the asthma data well. From the histogram for each component, we see from Figure 10 that component 1 and 2 represents shorter state persistence times while components 3 and 4 represent longer durations.
R> hist(res.asm, xlim = c(0,1))
R> for(comp in 1:4) {
hist(res.asm, comp = comp)
}
Figure 10: Results for the asthma data set showing the histograms of ISIs with the estimated posterior gamma mixture density (left) and histograms of ISIs for each mixture component (right).
We investigate the covariates’ influence by examining the mixture
probabilities from the last MCMC iteration. An interesting discovery is
that the mixture probabilities for both sexes, BMI levels, and severity
levels are the same, indicating that the levels of these three
covariates do not strongly influence the distributions of state
durations. This matches the global test results in
Figure 7. If the preceding state is state 1,
which is optimal control for asthma, the state duration time is longer
than that if the previous state is 2 or 3, as there is a lower weight in
component 1 and higher weight in components 3 and 4. On the other hand,
if the preceding state is 3, which is unacceptable control, then the
state duration time is much shorter, as the mixture probability in
component 1 is much higher when the preceding state is state 3. We
present some examples of the diagnostic plots for the asthma
data set
in Figure 11.
R> sm.asm$dur.mix.probs
$Severity
Mild-Moderate Severe
Comp 1 0.37 0.37
Comp 2 0.26 0.26
Comp 3 0.20 0.20
Comp 4 0.17 0.17
$BMI
BMI<25 BMI>=25
Comp 1 0.37 0.37
Comp 2 0.26 0.26
Comp 3 0.20 0.20
Comp 4 0.17 0.17
$Sex
Women Men
Comp 1 0.37 0.37
Comp 2 0.26 0.26
Comp 3 0.20 0.20
Comp 4 0.17 0.17
$prev_state
1 2 3
Comp 1 0.27 0.33 0.51
Comp 2 0.23 0.33 0.23
Comp 3 0.31 0.20 0.09
Comp 4 0.19 0.15 0.17
R> diag.BMRMM(res.fp2, cov.combs = list(c(1, 2, 1)),
transitions = list(c(1, 1)), components = c(3))
Figure 11: Results for the asthma data set showing the MCMC diagnostic plots, including traceplots and autocorrelation plots.
We presented the BMRMM package which implements both Bayesian Markov
mixed models (BMMM) for analyzing the state transitions and Bayesian
Markov renewal mixed models (BMRMM) for additionally analyzing the
duration times (being either state persistence times or inter-state
intervals) in a collection of categorical sequences, using flexible
Dirichlet and gamma mixtures, respectively. The BMRMM takes into account
fixed effects of the associated covariates as well as random effects of
the associated individuals while simultaneously selecting the
significant covariates separately for the state transitions and the
duration times. The package includes a synthetic foxp2
data set to
demonstrate the data framework and function usages. The package also
provides a series of plotting functions for visualizing the results of
the analyses, including various global and local hypotheses tests, MCMC
diagnostics, etc. We are committed to maintaining and further developing
the package in the future. [Future improvements to the package may
include more options for the distribution types of transition
probabilities and duration times beyond the currently available mixture
Dirichlet and mixture gamma distributions,
respectively.]
We thank two anonymous reviewers very much for their careful review of our work and their constructive comments that led to significant improvements to both the package and this paper. ::: :::::::::
Supplementary materials are available in addition to this article. It can be downloaded at RJ-2024-011.zip
BMRMM, markovchain, msm, SemiMarkov, SMM, mhsmm, hhsmm, Countr, ARPobservation, catdap, vcd
Distributions, Finance, MissingData, Survival, TeachingStatistics
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Wu & Sarkar, "BMRMM: An R Package for Bayesian Markov (Renewal) Mixed Models", The R Journal, 2025
BibTeX citation
@article{RJ-2024-011, author = {Wu, Yutong and Sarkar, Abhra}, title = {BMRMM: An R Package for Bayesian Markov (Renewal) Mixed Models}, journal = {The R Journal}, year = {2025}, note = {https://doi.org/10.32614/RJ-2024-011}, doi = {10.32614/RJ-2024-011}, volume = {16}, issue = {1}, issn = {2073-4859}, pages = {192-211} }