Combining Mixed Effects Hidden Markov Models with Latent Alternating Recurrent Event Processes to Model Diurnal Active-Rest Cycles

In our two state setting, rates of transition from state $s$ to the other state are defined as $\lambda_{s}(t_{ij})=\exp\left(\mathbf{x}^{\top}(t_{ij})\boldsymbol{\beta}_{s}\right)$. For example, $\lambda_{1}(t_{ij})$ denotes the rate of transition from state 1-to-2 (active-to-rest). Alternating recurrent event PH models often need to account for the longitudinal nature of the data, i.e., repeated measurements. Mixed effects or frailties can be used to account for the recurrent nature of the data (Wang et al., 2020; McGilchrist and Aisbett, 1991). Modifying the standard exponential PH model with a shared log-normal frailties or normal random intercepts, state transition hazards become $\lambda_{s}(t_{ij})=\exp\left(\eta_{s}(t_{ij})\right)=\exp\left(\mathbf{x}^{\top}(t_{ij})\boldsymbol{\beta}_{s}+\mathbf{z}^{\top}(t_{ij})\mathbf{b}_{s}\right)$, where $\mathbf{b}_{s}\sim\mathrm{N}(\mathbf{0}_{24},\sigma^{2}_{s}\mathbf{I}_{24})$. Here $\mathbf{z}(t_{ij})$ are 24 hour-of-day indicators, one-hot vectors designed to toggle the appropriate random intercepts within $\{\mathbf{b}_{1},\mathbf{b}_{2}\}$.

A HMM for individual $i$ is given by the complete data likelihood

where $\mathbf{A}_{i}$ are the true the state labels. The PH likelihoods for state transitions are given as

where $f(\Delta(t_{ij})|\lambda_{s}(t_{ij}))=\lambda_{s}(t_{ij})\exp\left(-\lambda_{s}(t_{ij})\Delta(t_{ij})\right)$ and $S(\Delta(t_{ij})|\lambda_{s}(t_{ij}))=\exp\left(-\lambda_{s}(t_{ij})\Delta(t_{ij})\right)$ are derived from the exponential distribution. Note that $\prod_{s=1}^{2}L(\boldsymbol{\beta}_{s},\sigma^{2}_{s},\mathbf{b}_{s}|\mathbf{A}_{i})f(\mathbf{b}_{s}|\sigma^{2}_{s})$ is the likelihood of an alternating event process (Król and Saint-Pierre, 2015; Wang et al., 2020). We denote indicators for state 1-to-2 transitions $d_{1}(t_{ij})=\mathbb{I}\left[A(t_{i(j-1)})=1,A(t_{ij})=2\right]$, as $\mathbf{d}_{1}=\left\{d_{1}(t_{11}),d_{1}(t_{12}),\dots,d_{1}(t_{ij}),\dots\right\}$ and 2-to-1 transitions as $\mathbf{d}_{2}$. We interpret failure to transition out of state 1, $c_{1}(t_{ij})=\mathbb{I}\left[A(t_{i(j-1)})=1,A(t_{ij})=1\right]$ as censoring, denoted as $\mathbf{c}_{1}=\left\{c_{1}(t_{11}),c_{1}(t_{12}),\dots,c_{1}(t_{ij}),\dots\right\}$ and we similarly define $\mathbf{c}_{2}$. The screen-on count state conditional Poisson likelihood is given as $\begin{array}[]{c}L\left(\mu_{1},\mu_{2}|\mathbf{A}_{i}\right)=\prod^{n_{i}}_{j=0}\prod^{2}_{s=1}f\left(y(t_{ij})|\mu_{s}\right)^{u_{s}(t_{ij})}\end{array}$ where state memberships $\mathbf{u}$, are denoted as indicators $u_{s}(t_{ij})=\mathbb{I}\left[A(t_{ij})=s\right]$. Since the true labels are unknown, $\mathbf{d}_{s}$, $\mathbf{c}_{s}$, and $\mathbf{u}$ are latent variables and (1) becomes a mixture model.

The log-likelihood of (1) are linear functions of latent variables $\mathbf{d}_{s}$, $\mathbf{c}_{s}$, and $\mathbf{u}$, lending our optimization approach to an EM algorithm (Dempster et al., 1977). Through the EM algorithm, indicators $\mathbf{d}_{s}$, $\mathbf{c}_{s}$, and $\mathbf{u}$ are imputed as continuous probabilities, in the process of obtaining maximum likelihood estimates (MLEs). As a result, the alternating recurrent event exponential PH model reduces to two weighted frailty models. We denote PH model weights as $\mathbf{w}(t_{ij})=\left\{{c}_{1}(t_{ij}),{d}_{1}(t_{ij}),{c}_{2}(t_{ij}),{d}_{2}(t_{ij})\right\}$, which belong to a 4-dimensional probability simplex, i.e., values are non-negative and $\|\mathbf{w}(t_{ij})\|_{1}=1$. Poisson mixture model weights $\mathbf{u}(t_{ij})=\left\{{u}_{1}(t_{ij}),u_{2}(t_{ij})\right\}$ belong to a 2-dimensional probability simplex. Our EM algorithm iteratively estimates the weights $\mathbf{w}(t_{ij})$ and $\mathbf{u}(t_{ij})$ using the forward-backward algorithm of Baum et al. (1970) and $\left\{\boldsymbol{\beta}_{1},\mathbf{b}_{1},\sigma^{2}_{1},\boldsymbol{\beta}_{2},\mathbf{b}_{2},\sigma^{2}_{2}\right\}$ using survival modeling. While the complete data likelihood is written as an alternating event processes, we see an equivalence with logistic regression transition probabilities in the E-step calculations, effectively retooling alternating event processes to fit HMMs to data that have heterogeneous event times transitions.

The estimated relative risk or coefficient estimates from PH and logistic regression have been shown to be similar under a variety of situations (Abbott, 1985; Ingram and Kleinman, 1989; Thompson Jr, 1977; Callas et al., 1998). Many useful properties of PH modeling can be leveraged in Markov chains which can improve the robustness of parameter estimation and reduce computational burden. Though many analyses look at the Cox PH case, we can adapt these approaches to the exponential PH which is a special case of the Cox PH. Following the derivations of Abbott (1985), we have $\lambda(t)=-d\log\{S(t)\}/dt$ and $S(\Delta(t))=\exp\left\{-\exp\left(\zeta_{0}(\Delta(t))+\sum_{k}\beta_{k}x_{k}\right)\right\}$ with $\zeta_{0}(\Delta(t))=\log\left\{\int_{0}^{\Delta(t)}\lambda_{0}(z)dz\right\}=\log(\Delta(t))+\beta_{0}$ where $\lambda_{0}(z)=\exp(\beta_{0})$ is the baseline hazard from the exponential distribution. When the event times are discrete, $\Delta(t)=1$ can be arbitrarily assigned as the event time and $\zeta_{0}(\Delta(t))=\beta_{0}$. Under this setting, $\log\left\{-\log\left(S(\Delta(t_{ij}))\right)\right\}=\eta_{s}(t_{ij})$ and the Taylor expansion of the survival function $S(\Delta(t_{ij})=1\mid\lambda_{s}(t_{ij}))=\exp(-\lambda_{s}(t_{ij}))$ at $\lambda_{s}(t_{ij})=0$ results in the power series $\exp(-\exp(\eta_{s}(t_{ij})))=1-\lambda_{s}(t_{ij})+\lambda_{s}(t_{ij})^{2}/2!-\lambda_{s}(t_{ij})^{3}/3!\ldots+(-1)^{n}\lambda_{s}(t_{ij})^{n}/n!\ldots=1-\lambda_{s}(t_{ij})+R_{\text{PH}}\left(\lambda_{s}(t_{ij})\right)$ where $1-P_{ij}=S(\Delta(t_{ij})=1\mid\lambda_{s}(t_{ij}))\approx 1-\lambda_{s}(t_{ij})$ when the incidence rate $\lambda_{s}(t_{ij})\approx 0$. The respective power series of $1-\text{expit}(\eta_{s}(t_{ij}))$ is given as $1-\text{expit}(\eta_{s}(t_{ij}))=1-\lambda_{s}(t_{ij})+\lambda_{s}(t_{ij})^{2}-\lambda_{s}(t_{ij})^{3}\ldots+(-1)^{n}\lambda_{s}(t_{ij})^{n}\ldots=1-\lambda_{s}(t_{ij})+R_{\text{log}}\left(\lambda_{s}(t_{ij})\right)$ where we set $\exp(-\exp(\eta_{s}(t_{ij})))=1-\text{expit}(\eta_{s}(t_{ij}))=1-P_{ij}\approx 1-\lambda_{s}(t_{ij})$ when the incidence rates are low. Note that in the case of low incidence rates, $R_{\text{PH}}\left(\lambda_{s}(t_{ij})\right)\leq R_{\text{log}}\left(\lambda_{s}(t_{ij})\right)$ and writing $\eta_{s}(t_{ij})$ as a function of $R_{\text{PH}}\left(\lambda_{s}(t_{ij})\right)$ or $R_{\text{log}}\left(\lambda_{s}(t_{ij})\right)$, we get $\log\left(P_{ij}+R_{\text{PH}}\left(\lambda_{s}(t_{ij})\right)\right)\leq\log\left(P_{ij}+R_{\text{log}}\left(\lambda_{s}(t_{ij})\right)\right)$. As noted in previous studies, PH and logistic regression estimate similar relative risk under a group event time, low incidence rate setting (Abbott, 1985; Ingram and Kleinman, 1989). In many practical settings, discrete-time applications of binary outcome models involve low incidence rates, letting exponential PH regression serve as an alternative to logistic regression. The analogous Markov chain with rare outcomes is a process with an extended stay in the current state, commonly found in diurnal biological processes. However, when event times are heterogeneous, then PH regression is the more appropriate choice for estimating relative risk. As noted in Abbott (1985), in many cases, we observed shrinkage towards zero for the relative risk under PH regression when compared to the logistic regression due to the inequality between the remainder terms $R_{\text{PH}}\left(\lambda_{s}(t_{ij})\right)$ and $R_{\text{log}}\left(\lambda_{s}(t_{ij})\right)$ (Ingram and Kleinman, 1989; Thompson Jr, 1977; Callas et al., 1998).

We further expand on this observed shrinkage by showing that the exponential PH regression is a penalized logistic regression. The exponential canonical form of exponential PH and logistic regression is given as $L\left(\boldsymbol{\beta}_{s}\mid\mathbf{A}_{i}\right)=\prod^{n_{i}}_{j=1}\exp\left\{d_{s}(t_{ij})\eta_{s}(t_{ij})\right\}\exp\left\{-\Psi\left(\eta_{s}(t_{ij})\right)\right\},\$ with $d_{s}(t_{ij})\in\{0,1\}$. Here, $\Psi_{\text{log}}(\eta)=\log(1+\exp(\eta))$ for logistic regression and $\Psi_{\text{PH}}(\eta)=\Delta(t)\exp(\eta)=\exp(\eta)$ for exponential PH regression under the discrete-time setting. We see that the logistic regression likelihood is uniformly bounded below by the PH likelihood due to $\exp(-\log(1+\exp(\eta)))>\exp(-\exp(\eta))$, which simply follows from the inequality $\log(1+z)<z$ for $z>0$. The convex optimization problems of PH and logistic regression are to minimize loss functions: $J_{\text{log}}(\boldsymbol{\beta}_{s})=-\log L_{\text{log}}\left(\boldsymbol{\beta}_{s}\mid\mathbf{A}_{i}\right)$ and $J_{\text{PH}}(\boldsymbol{\beta}_{s})=-\log L_{\text{PH}}\left(\boldsymbol{\beta}_{s}\mid\mathbf{A}_{i}\right)$. There is a positive convex penalty difference between the logistic and PH optimization problems, $J_{\text{PH}}(\boldsymbol{\beta}_{s})=J_{\text{log}}(\boldsymbol{\beta}_{s})+\mathcal{P}(\boldsymbol{\beta}_{s})$ and $\mathcal{P}(\boldsymbol{\beta}_{s})=\sum^{n_{i}}_{j=1}-\log(1+\exp(\eta_{s}(t_{ij})))+\exp(\eta_{s}(t_{ij}))>0$ with the penalty as the difference of $\Psi_{\text{log}}(\eta)$ and $\Psi_{\text{PH}}(\eta)$. It is straight forward to see that $\mathcal{P}(\boldsymbol{\beta}_{s})$ is convex with hessian $\nabla_{\boldsymbol{\beta}_{s}}^{2}\mathcal{P}(\boldsymbol{\beta}_{s})=\mathbf{X}^{\top}\boldsymbol{\Omega}\mathbf{X}\succeq 0$ and $\boldsymbol{\Omega}$ as a diagonal matrix with elements $\exp(\eta_{s}(t_{ij}))\left[1-\{1+\exp(\eta_{s}(t_{ij}))\}^{-2}\right]>0$. Under a simple parameterization, such as the intercept only model, it follows that the penalty $\mathcal{P}(\boldsymbol{\beta})$ favors a low incidence rate. The function $-\log(1+\exp(\eta))+\exp(\eta)$ is relatively flat for $\eta<0$ and dominated by $\exp(\eta)$ for $\eta\gg 0$. However, in practice this will shrink coefficients to the solution of $\text{argmin}_{\boldsymbol{\beta}_{s}}\{\mathcal{P}(\boldsymbol{\beta}_{s})\}$ which tends to be near $\boldsymbol{\beta}_{s}=0$ when considering that each $\eta_{s}(t_{ij})=\mathbf{x}^{\top}(t_{ij})\boldsymbol{\beta}_{s}$ is a different linear combination of $\boldsymbol{\beta}_{s}$. The penalty $\mathcal{P}(\boldsymbol{\beta}_{s})$ modifies the convex hull of the logistic regression loss function to induce shrinkage and results in a PH loss function which has many readily available software for estimation.

As noted in our E-step, the conditional expectations of the EM algorithms reduces to transition probability matrices that are comprised of logistic regression probabilities. In the case of a HMM with 3 or more states, the E-step reduces to transition probabilities comprised of multinomial logistic regression probabilities. However, the M-step conveniently involves fitting several independent exponential PH models. In order to illustrate this point, we define the complete data likelihood of transitioning out of state 1 in a 3 state HMM as,

with $d_{1m}(t_{ij})=\mathbb{I}\left[A(t_{i(j-1)})=1,A(t_{ij})=m\right]$ for $m\in\{1,2,3\}$, $\sum_{m=1}^{3}d_{1m}(t_{ij})=1$ and $\lambda_{1k}(t_{ij})=\exp(\eta_{1k}(t_{ij}))=\exp(\mathbf{x}^{\top}(t_{ij})\boldsymbol{\beta}_{1k})$ for $k\in\{2,3\}$. The minimum of two or more independent exponential random variables follows an exponential distribution with a new rate equal to the sum of rates, in our case: $\sum_{k=2}^{3}\lambda_{1k}(t_{ij})$. The likelihood contribution of staying in the same state is the survival function $S\left(\Delta(t_{ij})|\sum_{k=2}^{3}\lambda_{1k}(t_{ij})\right)=\prod_{k=2}^{3}S\left(\Delta(t_{ij})|\lambda_{1k}(t_{ij})\right)=\prod_{k=2}^{3}\exp(\Delta(t_{ij})\lambda_{1k}(t_{ij}))$. Continuing our example, calculating the E-step for the transition from state 1 to 2 is given as multinomial probability $\mathbb{E}\left[d_{12}(t_{ij})\right]=\lambda_{12}(t_{ij})/\left\{1+\sum_{k=2}^{3}\lambda_{1k}(t_{ij})\right\}$ where we leave the derivation details to the Web Appendix B. Note that in the 3 state HMM, E-step imputations are constrained to a 9-dimensional simplex. The M-step for the 3 state HMM parameters from likelihood (2), can be estimated with 2 independent exponential PH models. This EM procedure can be extended to an arbitrary number of states. In the M-step we fit independent weighted PH models rather than a cumbersome multinomial regression.

As noted in previous work, when the event times are heterogenous, PH regression is the correct model for estimating relative risk (Abbott, 1985; Ingram and Kleinman, 1989; Thompson Jr, 1977; Callas et al., 1998). However, PH and logistic regression estimate similar relative risk for the discrete and grouped event time setting with low incidence rates. We showed that under the discrete-time setting, exponential PH regression is an implicitly penalized logistic regression resulting in shrinkage of the relative risk estimates. This is desirable in many situations, specifically in our case where we are building HMMs, a model with a complicated error in response mechanism. The penalty $\mathcal{P}(\boldsymbol{\beta}_{s})$ slightly favors low probabilities of transitioning out of a state, which is useful to mitigate false positives. As a result, penalization shrinks the transition probability matrices $\boldsymbol{\Gamma}(t_{ij})$, towards the identity matrix and favors an extended sojourn time. Of note, post estimation HMM procedures such as the Viterbi algorithm for finding the most likely path also favors low probabilities of transitioning out of a state (Viterbi, 1967; Forney, 1973). In addition, we also invoke mixed effect modeling which numerically benefits from penalization. Next, fitting multiple independent weighted frailty models in the M-step is a more straight forward procedure than fitting an analogous weighted mixed multinomial logistic regression. The PH model is the correct model for heterogenous event times, but even in the case of discrete-time transitions, PH modeling’s implicit penalization has several useful properties for estimating HMMs.

We denote the alternating recurrent event exponential PH model outlined in Section 2.2 as PH-HMM. For our competing methods, we use the CT-HMM and discrete-time mixed effect logistic regression HMM (denoted as DT-HMM), described in our Web Appendix A. In addition, we define a two step estimator based on Poisson mixture model (PMM) maximum a posteriori (MAP) estimates in the Web Appendix A. All competing methods are initialized using PMM MAP estimates for state labels and EM is repeated until $\|\Theta^{(l+1)}-\Theta^{(l)}\|_{1}\leq 10^{-4}$.

For the alternating event process, we use a 24h period sine function as our time varying covariate, $x(t)=\sin\left({2\pi t}/{24}\right)\in[-1,1]$ and independently draw censoring times from a uniform distribution, $r_{ij}\sim U(0,h_{\text{max}})$. For the PH models we sequentially increment $t_{ij}$ by $\Delta(t_{ij})$ and simulate covariates $x(t_{ij})$. We simulate the complete data likelihood with multiple individuals but shared $\boldsymbol{\beta}_{s}$, by drawing from