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

Abstract

Data collected from wearable devices can shed light on an individual's pattern of behavioral and circadian routine. Phone use can be modeled as alternating processes, between the state of active use and the state of being idle. Markov chains and alternating recurrent event models are commonly used to model state transitions in cases such as these, and the incorporation of random effects can be used to introduce diurnal effects. While state labels can be derived prior to modeling dynamics, this approach omits informative regression covariates that can influence state memberships. We instead propose an alternating recurrent event proportional hazards (PH) regression to model the transitions between latent states. We propose an expectation–maximization algorithm for imputing latent state labels and estimating parameters. We show that our E-step simplifies to the hidden Markov model (HMM) forward–backward algorithm, allowing us to recover an HMM with logistic regression transition probabilities. In addition, we show that PH modeling of discrete-time transitions implicitly penalizes the logistic regression likelihood and results in shrinkage estimators for the relative risk. This new estimator favors an extended stay in a state and is useful for modeling diurnal rhythms. We derive asymptotic distributions for our parameter estimates and compare our approach against competing methods through simulation as well as in a digital phenotyping study that followed smartphone use in a cohort of adolescents with mood disorders.

1 Introduction

Diurnal and circadian rhythm studies often model physiological processes as periodic cycles, such as a person's active and rest cycle. Sleep and diurnal rhythm are essential components of many circadian physiological processes with a clear time-of-day effect on active and rest cycles (Lagona et al., 2014). Here, we consider this problem of estimating an individual's cycles between active and rest states in a mobile health (mHealth) setting based on wearable device or smartphone sensor data. If the true state labels are known, then a Markov chain can be used to model state transitions over time, otherwise a hidden Markov model (HMM) can be used to simultaneously perform classification and state transition estimation (Langrock et al., 2013). In addition, HMMs have been extended to incorporate time-of-day effects as periodicity or seasonality using random effects (Bartolucci & Farcomeni, 2015; Holsclaw et al., 2017; Stoner & Economou, 2020). Continuous-time hidden Markov models (CT-HMM) have also been used for state classification in similar contexts but have difficulty accounting for random effects (Bartolucci & Farcomeni, 2019; Bureau et al., 2003; Jackson et al., 2003). Most mixed effects HMM estimation procedures estimate logistic regression for discrete-time processes and do not account for time between state transitions (Maruotti & Rocci, 2012).

It is important to note that active and rest state transitions are ergodic processes and the sojourn time between these states can be modeled with a proportional hazards (PH) regression in both directions, active-to-rest and rest-to-active. These two directions of transitions can be viewed as an alternating recurrent event PH model (Król & Saint-Pierre, 2015; Shinohara et al., 2018; Wang et al., 2020). However, the ability for alternating recurrent event processes to accurately model sojourn times complements HMMs and provides many useful properties in addition to being computationally scalable. This novel model can be viewed as a latent state analog of an alternating recurrent event process (Król & Saint-Pierre, 2015; Wang et al., 2020). Mainly, if the underlying data-generating process is a continuous-time process, then PH models are a more appropriate modeling choice (Abbott, 1985; Ingram & Kleinman, 1989). If the data-generating process involves discrete-time transitions, we showed that PH modeling penalizes the probability of leaving a state.

We propose an approach that takes advantage of the strengths of both HMMs and alternating recurrent event models to jointly estimate latent states while simultaneously providing flexible modeling of sojourn times. Our expectation-maximization (EM) algorithm imputes latent active and rest state labels while modeling state transitions with an alternating recurrent event process using exponential PH regressions. Informative regression covariates, often omitted in state labeling, are incorporated into the latent state imputation using the EM algorithm. Under the EM algorithm, we show that the E-step in this case simplifies to imputations involving an HMM forward–backward algorithm, where state transition probabilities are defined as logistic or multinomial regression probabilities (Altman, 2007; Baum et al., 1970). We also show that the M-step reduces to fitting independent PH models weighted by E-step imputations allowing for a potentially large number of latent states, as well as, providing a means to obtain large-sample theory inference. Our EM approach involving PH models provides a scaleable M-step while returning multinomial regression transition probabilities commonly found in HMMs (Altman, 2007; Holsclaw et al., 2017; Maruotti & Rocci, 2012). Furthermore, we show that applying PH models to discrete-time transitions, implicitly penalizes logistic regression to shrink the transition probability matrices of HMMs toward the identity matrix and mitigates overfitting in many practical settings. As a result, the PH models favor processes with a low incidence of state transitions such as diurnal cycles, where we expect few state transitions within a 24 h period. This direction of penalization is in line with hidden semi-Markov models (HSMMs) extending the dwell time of staying in a state while providing a means to incorporate random effects (Zucchini et al., 2017).

We apply this approach to estimate active–rest diurnal cycles in a sample of patients with affective disorders using their passively collected smartphone sensor data, namely through the accelerometer, screen on/off data of patient smartphones and time-of-day random intercepts. This quantification of routine can be correlated with a myriad of relevant clinical outcomes, as the regularity of diurnal rhythms plays an important role in psychopathology with past studies having shown associations between irregular rhythms and adverse health outcomes (Monk et al., 1991). In addition, we fit a population level HMM to study effects of individual specific covariates on state transitions.

2 Data and Methods

Our data consist of individuals, with each individual i having a sequence of covariates to be modeled with separate HMMs. The HMM of the ith individual has sequences of active or rest states ⁠, where hourly time-stamps are increasing in j. In our example, we denote active and rest states as and ⁠, respectively, with an outline of our HMM in Figure 1. Within each sequence, we define event times for state transitions as ⁠, which follow from an exponential distribution and can be fitted with a PH regression. Linking multiple exponential event time processes results in a recurrent event model. Furthermore, recurrent exponential PH models are analogous with a non-homogeneous Poisson processes which retains independent increments, allowing us to chain multiple transitions together.

The covariates used in the exponential PH regression are mean acceleration magnitudes (Euclidean norm) from the preceding hour evaluated at transitions and are outlined in Figure 1. We denote the intercept and covariates as ⁠. We make an ergodic state transition assumption, where states will inevitably communicate with each other, that is, the active-to-rest and rest-to-active transitions will eventually occur. This allows the survival function to capture the likelihood contribution of when a state transition did not occur, meaning the transition will occur at some future time. Because we do not have the true state labels ⁠, we must rely on state-dependent observations. We use screen-on counts for each time-stamp as observations from state-dependent distributions ⁠, where ⁠, are state-specific parameters and we expect for the rest state.

2.1 Alternating Recurrent Event Proportional Hazards Model

In our two states setting, rates of transition from state s to the other state are defined as ⁠. For example, 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, that is, repeated measurements. Mixed effects or frailties can be used to account for the recurrent nature of the data (McGilchrist & Aisbett, 1991; Wang et al., 2020). Modifying the standard exponential PH model with a shared log-normal frailties or normal random intercepts, state transition hazards become ⁠, where ⁠. Here, are 24 hour-of-day indicators, one-hot vectors designed to toggle the appropriate random intercepts within ⁠.

An HMM for individual i is given by the complete data likelihood

where are the true state labels. The PH likelihoods for state transitions are given as

where and are derived from the exponential distribution. Note that is the likelihood of an alternating event process (Król & Saint-Pierre, 2015; Wang et al., 2020). We denote indicators for state 1-to-2 transitions ⁠, as and 2-to-1 transitions as ⁠. We interpret failure to transition out of state 1, as censoring, denoted as and we similarly define ⁠. The screen-on count state conditional Poisson likelihood is given as where state memberships are denoted as indicators ⁠. Since the true labels are unknown, ⁠, ⁠, and are latent variables and Equation (1) becomes a mixture model.

2.2 Expectation–Maximization Algorithm for Proportional Hazards Regression and Hidden Markov Model Parameters

The log-likelihood of Equation (1) is linear functions of latent variables ⁠, ⁠, and ⁠, allowing us to use an EM algorithm our optimization (Dempster et al., 1977). Through the EM algorithm, indicators ⁠, ⁠, and 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 ⁠, which belong to a four-dimensional probability simplex, that is, values are non-negative and ⁠. Poisson mixture model (PMM) weights belong to a two-dimensional probability simplex. Our EM algorithm iteratively estimates the weights and using the forward–backward algorithm of Baum et al. (1970) and 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.

2.2.1 E-step

In the E-step, we derive the expectation of conditional on model parameters and observed data in the ith HMM denoted by and ⁠, as

⁠, and

are forward and backward probabilities of an HMM. Vectors are of initial-state distribution probabilities for the HMM. The transition probability matrix ⁠, derived by normalizing the alternating recurrent event exponential PH models are

where ⁠. Under PH modeling of multiple-state transitions, EM updates simplify to multinomial transitions (see Web Appendix A). The weights from can more generally be written as ⁠. The E-step involves imputing through an HMM, using a forward–backward algorithm, where transition probabilities are standard logistic functions. We denote forward and backward probabilities vectors and ⁠, where and are 2 × 1 vectors. The state-dependent distributions are contained in the 2 × 2 diagonal matrix ⁠. Note that ⁠, and are used to normalize E-step probabilities. The E-step update for iteration ⁠, simplifies to calculating probabilities as

such that the sum of all elements is equal to one, ⁠. Operations ⊗ and ⊙ are Kronecker and Hadamard products, respectively. The update for the PMM weights are ⁠, such that ⁠.

2.2.2 M-step

The M-step update for involves solving the mixture model ⁠, where solutions are known for most distributions and in the Poisson setting. From ⁠, we have updates ⁠. The update for involves fitting two frailty models for which can be accomplished by recognizing can be factored into indicators and log-likelihood weights or case weights. The weights can be interpreted as the probability that a specific state transition, for example, if ⁠, then an active-to-rest transition likely occurred. We duplicate each row of data to be both a transition event and a censored outcome and then weight the rows by ⁠, and ⁠, respectively. A data augmentation example is outlined in Web Table 1.

While there are numerous survival packages in R, we need a package that can incorporate weights, parametric PH models and normally distributed random intercepts (Therneau & Lumley, 2015). The R package tramME can be used to fit our weighted exponential parametrization of shared log-normal frailty models (Tamási & Hothorn, 2021). The package tramME uses a transformation model approach combined with an efficient implementation of the Laplace approximation to fit the shared log-normal frailty models (Hothorn et al., 2018; Hothorn, 2020; Kristensen et al., 2016). By updating we also update our transition rates, and which are used to calculate ⁠.

The M-step updates for the initial-state distribution are ⁠, such that ⁠. Finally, we iteratively calculate the E-step: ⁠, and M-step: until convergence to obtain the MLEs. Estimating shared population parameters can be done by taking the product of all individual likelihoods and applying EM.

2.3 Comparison of Relative Risks from Proportional Hazards and Logistic Regression

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; Callas et al., 1998; Ingram & Kleinman, 1989; Thompson Jr, 1977). 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 and with ⁠, where is the baseline hazard from the exponential distribution. When the event times are discrete, can be arbitrarily assigned as the event time and ⁠. Under this setting, and the Taylor expansion of the survival function at results in the power series where when the incidence rate ⁠. The respective power series of is given as where we set when the incidence rates are low. Note that in the case of low incidence rates, and writing as a function of or ⁠, we get ⁠. 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 & 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 toward zero for the relative risk under PH regression when compared to the logistic regression due to the inequality between the remainder terms and (Callas et al., 1998; Thompson Jr, 1977).

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 ⁠, with ⁠. Here, for logistic regression and 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 ⁠, which simply follows from the inequality for ⁠. The convex optimization problems of PH and logistic regression are to minimize loss functions: and ⁠. There is a positive convex penalty difference between the logistic and PH optimization problems, and with the penalty as the difference of and ⁠. It is straightforward to see that is convex with hessian and Ω as a diagonal matrix with elements ⁠. Under a simple parameterization, such as the intercept only model, it follows that the penalty favors a low incidence rate; this relationship is further discussed in Web Appendix E. The function is relatively flat for and dominated by for ⁠. However, in practice this will shrink coefficients to the solution of which tends to be near when considering that each is a different linear combination of ⁠. The penalty 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 an HMM with three 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 three-state HMM as:

with for ⁠, and for ⁠. 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: ⁠. The likelihood contribution of staying in the same state is the survival function ⁠. Continuing our example, calculating the E-step for the transition from state 1 to 2 is given as multinomial probability ⁠, where we leave the derivation details to the Web Appendix A and B. Note that in the three-state HMM, E-step imputations are constrained to a nine-dimensional simplex. The M-step for the three-state HMM parameters from likelihood (2), can be estimated with two independent exponential PH models. This EM procedure can be extended to an arbitrary number of states. The M-step fits independent weighted PH models rather than a cumbersome multinomial regression.

As noted in the previous work, when the event times are heterogenous, PH regression is the correct model for estimating relative risk (Abbott, 1985; Callas et al., 1998; Ingram & Kleinman, 1989; Thompson Jr, 1977). 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 the response mechanism. The penalty 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 ⁠, toward 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 lower probabilities of transitioning out of a state (Forney, 1973; Viterbi, 1967). 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.

2.4 Competing Methods

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 F. In addition, we define a two-step estimator by calculating PMM maximum a posteriori (MAP) estimates and fitting PH regression using MAP labels (Web Appendix F). All competing methods are initialized using PMM MAP estimates for state labels and EM is repeated until ⁠.

3 Asymptotic Results

The EM algorithm converges and obtains the MLEs ⁠, with the final likelihood evaluation being equivalent to weighted exponential frailty models. As a result, we can generalize existing large-sample theory asymptotic inference for mixed effect models for our weighted setting, using weights derived in Section 2.2.1.

 

4 Simulation Studies

For the alternating event process, we use a 24 h period sine function as our time-varying covariate, and independently draw censoring times from a uniform distribution, ⁠. For the PH models, we sequentially increment by and simulate covariates ⁠. We simulate the complete data likelihood with multiple individuals but shared ⁠, by drawing from:

and with ⁠, ⁠, and ⁠. Pseudocode for generating the alternating survival data can be found in Algorithm 1 of Web Appendix G. Similarly, for a discrete-time alternating event process, we simulate transition events as and increment by ⁠. Pseudocode for generating the discrete alternating event data can be found in Algorithm 2 of Web Appendix G. Once true state labels are obtained, we simulate state-dependent observations from Poisson distributions to complete the simulation of the HMMs.

We simulated 500 replicates using four different sets of parameters and used three methods for generating the data. For each set of parameters, we simulated data using Algorithm 1 with maximum censoring times set to and to evaluate performance under heterogeneous and grouped event times. We also used Algorithm 2 in order study models fitted to data simulated from a discrete-time process. We have a total of 12 different cases outlined in Table 1. In summary, Cases 1.1–1.3 looked at low incidence rate of state transition and a large distributional difference between and ⁠; Cases 2.1–2.3 looked at high incidence rate and large distributional difference; Cases 3.1–3.3 looked at low incidence rate and small distributional difference; and Cases 4.1–4.3 looked at high incidence rate and small distributional difference. For each case, we fit models using the PMM, DT-HMM, CT-HMM, and PH-HMM approaches. The E-step update can be computed quickly in parallel where each ith HMM is processed independently. We present mean accuracy of recovering the true state labels using the MAP and empirical standard error (SE), over 500 replicates in Table 1. We also present mean parameter estimates, empirical SEs and mean square errors (MSE) in Tables 2 and 3 and coverage results in Web Appendix H.

Cases 1.1–1.3 represent data that is often associated with HMMs. These non-stationary processes often involve an extended stay in a state after a transition with clear separation in their state-dependent distributions (Holsclaw et al., 2017). We found that PH-HMM performed well in Cases 1.1–1.3 settings, due to the PH and logistic likelihoods being similar for low incidence settings as described in Section 2.3. In Cases 1.1–1.3, a region where there is little difference in PH and logistic likelihoods, with a plot of the penalty difference given in Web Figure 1. When there are heterogeneous event times (Case 1.1), PH-HMM results in more accurate coefficient estimates by specifying the correct parametric form, while DT-HMM does not achieve the same accuracy. PMM, as a two-step estimation, is less accurate due to not incorporating covariate information throughout the estimation process. Under common data generative settings of low incidence rates and for both heterogeneous and discrete times, PH-HMM is a highly reliable and robust parameterization for fitting HMMs, while logistic regression DT-HMMs does not have comparable utility. In addition, PH regression, being a derivative of the Poisson regression link function, is more numerically stable in the rare outcome settings of Cases 1.1–1.3.

In Cases 2.1–2.3, settings with clear separation in their state-dependent distributions and high incidence of state transitions, PH-HMM has low bias when the data-generating process is a survival model and DT-HMM yields a poor fit in this setting. For Case 2.3, when the data-generating model is a logistic regression, PH-HMM results in biased parameter estimation. In Cases 2.1–2.3 where ⁠, parameter estimation for Case 2.3 from PH-HMM fits with similar results for ⁠. This attenuation is in line with Section 2.3, penalizing ⁠, noting that for ⁠. Though this induces bias, it is useful for robust estimation of HMMs for common non-stationary processes with few incidences of state transitions. This bias is reflected in low coverage rates which can be found in Web Appendix H. Our data-generating model sets to illustrate a case with many incidences of state transitions. In practice, an extended region of the ⁠, where and is a level of certainty seldom encountered in real data. Existing HMM methods such as Viterbi algorithm and HSMMs work to attenuate transition probabilities in order to mitigate instances of erroneous state transitions (Bulla et al., 2010; Viterbi, 1967). Our PH-HMM method estimates milder transition probabilities through the implicit penalization of PH regression without changing the covariate link function to the HSMM parameterization which can be difficult to interpret.

For Cases 3.1–4.3, we present simulation studies that mirror the regression parameters used in Cases 1.1–2.3 but with small distributional differences. As we close the gap between state-dependent distributions, we observe less accurate coefficient estimation overall. In practice, settings with small distributional differences often require additional consideration when deciding model parameterization. In general, DT-HMM tends to over estimate the range of ⁠, while PH-HMM tends to be more conservative with probability values of being closer to 0. The CT-HMM estimation procedure that calculates transition probabilities using continuous-time Markov chains work well in grouped and discrete-time settings. Additional discussion of simulation studies can be found in Web Appendix I.

From our simulation studies, we found that PH-HMM excels in a number of different cases. DT-HMM and CT-HMM are sensitive to heterogeneous event times leading to inaccurate state label recovery and bias parameter estimation. While DT-HMM is more accurate than PH-HMM at state label recovery in Cases 3.1 and 4.1, their parameter estimates have higher MSE and bias. DT-HMM generally has higher MSE than PH-HMM for cases where the data are generated from PH models (Cases 1.1, 1.2, 2.1, 2.2, 3.1, 3.2, 4.1, and 4.2). The shrinkage properties outlined in Section 2.3 carry over into our simulation study. In most cases and when the data are generated from a logistic regression, DT-HMM and PH-HMM have comparable accuracy and state-dependent distribution estimates, even though PH-HMM coefficients exhibit shrinkage. This shrinkage property has many uses which extend beyond simple simulation studies which we outline next in our real-data analysis examples.

5 Application: mHealth Data

A sample of individuals, recruited via the Penn/CHOP Lifespan Brain Institute or through the Outpatient Psychiatry Clinic at the University of Pennsylvania as part of a study of affective instability in youth (Xia et al., 2022). All participants provided informed consent to all study procedures. This study was approved by the University of Pennsylvania Institutional Review Board. For each individual, roughly 3 months of data was collected using the Beiwe platform (Torous et al., 2016). Accelerometer measures (meters per second squared) for x, y, and z-axes, screen-on events for Android devices and screen-unlock events for iOS (Apple) devices were acquired through the Beiwe platform—we refer to both as “screen-on events” in this paper. In general, we recommend a minimum of 30 days of data in order to fit our individual-specific model with hour-of-day random intercepts.

First, we analyzed the data under a discrete-time setting where we impute data for each missing hour. By collecting both screen-on events and accelerometer data, we are able to construct a missing at random (MAR) model for imputing missing accelerometer data. In our case, there is missing accelerometer data during of periods screen-on activity, that is, accelerometer data are missing over a given hour, while there are many screen-on events observed over the same period. On the other hand, there are periods of dormancy: accelerometer data are missing and there are also no screen-on events. More specifically periods of dormancy occur when accelerometer measurements are missing due to user and device-related factors such as the phone being powered off or being in the airplane mode. Periods of dormancy have greater probability of missing accelerometer features and are identified using a two-state HSMM with Bernoulli state-dependent distributions (Bulla et al., 2010). Missing mean acceleration magnitudes from dormant periods were imputed using the minimum of acceleration features (excluding outliers). While missing data assigned to the periods screen-on activity were imputed by regressing accelerometer features on over all hours where data are completely observed. For the heterogeneous event time example, we did not impute missing data and absorbed the duration of dormancy into the event times ⁠. Periods of consecutive missing acceleration magnitudes over 24 h constitutes an end of a Markov chain and the start of a new chain, where the likelihoods of multiple HMM sequences can be multiplied together for parameter estimation.

5.1 Estimating Strength of Routine in Youth with Affective Disorders

In psychiatric studies, regularity of a rhythm is defined as the association of time-of-day and state membership; the effect of hour-of-day on activity and rest state membership represents the diurnal rhythm. As an illustrative example, we fit a model with hour-of-day effects, as normally distributed random intercepts in our HMMs and for each individual we fit PH models ⁠, where are 24 h-of-day random intercepts. By fitting hour-of-day random intercepts, and are each of length 24 and only if is in the rth h of the day. Rates of transition from active-to-rest states are given as ⁠, rates of transition from rest-to-active states are given as and an example of PH-HMM outputs can be found in Figure 2. The variances of the random intercepts can be interpreted as an L2 penalty on hour-of-day effects and disappears as ⁠. We quantify the strength of diurnal effects for an individual, by looking at the variances ⁠, where large variances correspond with large hour-of-day effect sizes and greater regularity in diurnal rhythms with an example in Figure 2.

5.2 Population HMM: Differences between Operating Systems

In addition, we can fit a population model, with random intercepts being specific to each individual through estimation using the likelihood ⁠. However, for iOS devices we only have screen-unlock events, that is, entering in a passcode to unlock the phone. Android devices have screen-on events which occur when the phone screen turns on such as when receiving a message; the phone does not need to be unlocked for the screen to be turned on. Screen-unlock events are less frequent and a subset of screen-on events, causing the counts ⁠, to be lower for iOS devices. The relationship between acceleration, and ⁠, may be experience effect modification due to the operating system (OS). We can test interaction between OS and acceleration, while controlling for the interaction with user sex in our regression and other individual effects with random intercepts. Android devices and males serve as baseline in this analysis. For the active-to-rest model: and rest-to-active model: ⁠, ⁠, where ⁠, are 41 individual specific random intercepts and are individual specific indicators. We fit our competing methods and test interaction for the discrete and heterogeneous event time settings, where estimates can be found in Table 4. Sensitivity analysis for our discrete-time model with MAR missing data imputation can be found in Web Appendix J.

We found that there was no significant interaction between OS and acceleration in the discrete-time active-to-rest model but did find significant interaction in the heterogeneous event time active-to-rest model. In addition, we found that there was significant interaction between OS and acceleration in the discrete-time rest-to-active model but did not find significant interaction in the heterogeneous event time rest-to-active model. During the active state, we know that counts, ⁠, are lower for iOS devices. In order to compensate for lower screen-on counts in iOS devices, iOS rest-to-active transitions require a higher magnitude of acceleration to achieve the same transition rate as an Android device, hence the negative sign of β₂₃. During the rest state, we know that iOS devices have zero-inflated counts, where excess zeros are due to not being able to record a screen-on event. In the active-to-rest model, we may achieve the same transition rate while having higher acceleration in iOS devices than Android devices, hence the positive sign of β₁₃. Our results suggest that the magnitude of effect acceleration has on state transition depends on OS, but further investigation is needed.

For the discrete-time setting, we estimated active-state distribution and rest-state distribution ⁠. For the heterogeneous event time setting, we estimated and ⁠. Screen-on counts separate into stark clusters, where μ₁ and μ₂ are similar to the large distributional difference from our simulation study. We found that rate of transition from active-to-rest states is negatively associated with acceleration; rate of transition for rest-to-active states is positively associated with acceleration and is statistically significant (⁠⁠) for all competing methods. These HMM parameter estimates related to screen-on counts and acceleration align with common intuition. For the PMM method, modeling the MAP estimates first and then combining the estimates to obtain a population level model, does not account for acceleration when imputing state labels and resulted in a poor fit. CT-HMM and PH-HMM estimates are comparable and aligned with the findings from the simulation study. In addition to large distributional differences, diurnal active–rest cycles are expected to be low-incidence rate processes, where we anticipate a few transition between states in a 24 h period. We observed that the magnitude of the relative risk estimates for CT-HMM and PH-HMM is comparable to each other but less than the DT-HMM estimates. This difference becomes more noticeable for the heterogeneous event time setting, which DT-HMM is not equipped to handle. For many parameters, DT-HMM relative risk estimates are three times that of PH-HMM estimates, while estimates for state-dependent distribution parameters are similar. PH-HMM allows us to achieve comparable estimation of state-dependent distributions while leveraging shrinkage to avoid overfitting coefficients.

6 Discussion

For a latent state setting, we proposed a method for estimating alternating recurrent event exponential PH model with shared log-normal frailties using the EM algorithm. Our E-step imputations involves a discrete-time HMM using logistic or multinomial regression transition probabilities with normally distributed random intercepts. The HMM obtained during the E-step of our EM algorithm is an alternative method for estimating mixed HMM which are typically obtained by estimating logistic or multinomial regressions (Altman, 2007; Maruotti & Rocci, 2012). The M-step conveniently involves fitting several independent PH models rather than multinomial regression with many states. In addition, we showed that the PH model applied to the discrete-time setting is a penalized logistic regression which shrinks transition probability matrices toward the identity matrix. Our framework can accommodate random intercepts to account for longitudinal data, such as data collected from the same individual or hour-of-day periodic effects. We derived asymptotic distributions for the PH regression coefficients and random intercepts, where coefficients have a hazard ratio interpretation akin to the Cox PH model. Our PH-HMM approach is a flexible method for modeling complex mHealth datasets, where heterogeneous event times can be incorporated into the PH regression while accounting for latent states. If the underlying data are a discrete-time process, then PH modeling offers penalization to mitigate overfitting, otherwise PH models are more appropriate for heterogeneous event time processes. Detail coverage rate results are presented in Web Appendix H and show that DT-HMM performs favorably in settings with frequent state changes, despite the PH-HMM excelling in numerous other situations.

Though we focused on alternating two-state models, our proposed framework can be generalized to arbitrary number of states. However, we caution that the number of model parameters grows quadratically with the number of states. In this case, model selection criteria such as Akaike information criterion/Bayesian information criterion can be used to inform the number of states that best fit the data. Furthermore, we provide R code for our analyses in Web Appendix K, which can be easily extended to consider model selection in addition to alternative parameter initialization.

As for our MAR assumption in our imputation model, we demonstrate via a sensitivity analysis that this assumption is robust even in the presence of moderate MNAR-induced systematic bias. Web Appendix J shows robust estimates for regression coefficients in the presence of biased imputation. Overall we present a flexible regression procedure that can accommodate different degrees of missing data, parameterization of random effects, and the use of other statistical structures such as semiparametric regression and multilevel models. However, parameter initialization, computational complexity, and interpretability should be considered when parameterizing complicated models such as HMMs (Maruotti & Punzo, 2021). Finally, a strong appeal of our approach is that complicated statistical structures can be incorporated into independent PH regressions which then ultimately simplify to multinomial transition probabilities during the E-step. Indeed, our application to quantify routine activity data is but one example of a framework that can flexibly adapt to a variety of different data types and corresponding modeling choices.

Data Availability Statement

The data that support the findings in this paper are available on request from the corresponding author (Benny Ren [email protected]). The data are not publicly available due to privacy or ethical restrictions.

Acknowledgments

This work was supported by grants from the NIMH: R01MH116884. Support for data collection was provided by the AE Foundation: R01MH113550 and R01MH107703. The authors wish to thank the Lifespan Informatics & Neuroimaging Center for providing the data

References