Continuous-Time Hidden Markov Factor Model for Mobile Health Data: Application to Adverse Posttraumatic Neuropsychiatric Sequelae

Authors: lin  ge
xinming  an
donglin  zeng
samuel  mclean
ronald  kessler
rui  song
Venue: IEEE_org_EMBS_BHI_2024_Conference
Published: 2024-06-14 03:43:27

Abstract: Adverse posttraumatic neuropsychiatric sequelae (APNS) are common among veterans and millions of Americans after traumatic exposures, resulting in substantial health and financial burdens for trauma survivors, their families, and society. Despite numerous studies conducted on APNS over the past decades, there has been limited progress in understanding the underlying neurobiological mechanisms due to several unique challenges. One of these challenges is the reliance on subjective self-report measures to assess APNS, which can easily result in measurement errors and biases (e.g., recall bias). To mitigate this issue, in this paper, we investigate the potential of leveraging objective longitudinal mobile device data to identify homogeneous APNS states and study the dynamic transitions among them and potential risk factors after trauma exposure. To handle the unique challenges posed by longitudinal mobile device data, we developed a continuous-time hidden Markov factor model and designed a Stabilized Expectation-Maximization algorithm for parameter estimation. Simulation studies were conducted to evaluate the performance of parameter estimation and model selection. Finally, to demonstrate the practical utility of the method, we applied it to mobile device data collected from the Advancing Understanding of RecOvery afteR traumA (AURORA) study. A Python implementation of the proposed method is available at https://anonymous.4open.science/r/CTHMFM.

Main paper: 
 i. introduction: 
 Adverse posttraumatic neuropsychiatric sequelae (APNS) (e.g., pain, depression, and PTSD) are frequently observed in civilians and military veterans who have experienced traumatic events. These APNS increase the risk of chronic illnesses, including cancer and heart disease, and substantially contribute to drug abuse, suicide, and disability. Moreover, APNS impose enduring psychosocial and financial burdens not only on individuals with the disorder but also on their families, communities, and society as a whole. However, little progress has been made in advancing APNS research over the past few decades due to unique challenges. First, APNS have been evaluated through subjective self-reported measures, which lack objective reliability. Second, the heterogeneity among patients, as recognized in traditional classification and diagnoses, complicates the study of APNS. Lastly, these APNS disorders are often studied and treated independently, despite their frequent co-occurrence [1] . These obstacles hinder the identification of objective markers, the advancement in understanding the neurobiological mechanisms of APNS, and the development of effective preventative/treatment strategies.Identifying homogeneous states and exploring the dynamic prognosis of APNS in the immediate aftermath of trauma exposure holds promise for enhancing our understanding of APNS and identifying effective intervention options and appropriate timing at the individual level. Regrettably, due to the lack of appropriate data and effective statistical method, no large-scale studies have been conducted to investigate the onset, dynamic transitions (such as recovery and relapse), and associated risk factors of APNS. To help address the challenges, the National Institutes of Mental Health, joined by the US Army Medical Research and Material Command, developed the Advancing Understanding of RecOvery afteR traumA (AURORA) study [1] . This study gathered biobehavioral data from a large cohort of trauma survivors (n = 2,997) across the United States over a year, including self-reported surveys, web-based neurocognitive tests, digital phenotyping data (i.e., from wrist wearables and smartphones), psychophysical tests, neuroimaging assessments, and genomics data. In contrast to previous studies relying on self-report surveys or neuroimages [2] , our work aims to address the challenges of APNS by utilizing the objective digital phenotyping data that tracks individuals' behaviors, moods, and health statuses in real-time, real-life environments. Specifically, we developed a Hidden Markov Factor Model (HMFM) to analyze the mobile device data, allowing us to simultaneously identify homogeneous subtypes, investigate subtype-specific structures, and model individual progression with associated risk factors.Hidden Markov Models (HMMs) have been widely used in various fields [3] . However, mobile device data presents two unique challenges that standard HMMs cannot handle, including the interdependent variables with unknown interrelationship structures and unevenly spaced measurements.Mobile device sensor data, such as accelerometer and photoplethysmography (PPG) from smartwatches, involves highly intensive time series data. Instead of being directly used, these data are typically pre-processed to generate technical summaries that represent various characteristics of each time series variable, which are often highly correlated. As the number of features increases, the parameters in the covariance matrix expand exponentially, making a fully free covariance matrix impractical. Consequently, when implementing HMMs, features are usually assumed to be independent given the latent state membership, despite their correlation, potentially introducing estimation bias.To appropriately model the association between features, factor analysis models (FMs) [4] provide an efficient and parsimonious approach and have been incorporated into HMMs in various ways. For example, the factor analyzed hidden Markov model [5] combines an FM with a discrete-time HMM (DTHMM), which assumes evenly spaced measurements. It has been extensively used in a variety of real-world applications, including speech recognition [6] , environmental protection [7] , and seizure detection [8] . Similarly, [9] introduced the regime-switching factor model to handle high-dimensional financial market data. However, they all assume homogeneous transition probability matrices, limiting their ability to account for the heterogeneity of transition probabilities over time and among different subjects and explore risk factors of state transition. To simultaneously capture the interrelationships among observed features and account for the variability of transition probabilities, a joint framework incorporating HMM, FM, and a feature-based transition model was recently proposed [10] , [11] . However, it is not directly applicable to mobile device data. First, it employs a confirmatory factor model (CFM) with pre-specified structures for the factor loading matrices, which are often unknown a priori. Therefore, an exploratory factor model (EFM) is needed to explore the interrelationships among all observed features. Second, their framework assumes ordered states, which is inappropriate for our use case.Another challenge posed by mobile device data is the irregular spacing of measurements. For example, activity and heart rate variability (HRV) data were collected only when the participants wore the watches, resulting in non-uniformly spaced observations and significant variation in sampling schedules between individuals. While the aforementioned methods are all based on DTHMM, neglecting the impact of time gaps between consecutive observations on transition rates, continuoustime discrete-state HMM (CTHMM) was developed to handle irregularly spaced measurements [12] . CTHMM and its extensions that incorporate covariates to characterize transition rates are widely used in medical research [13] , [14] . However, none of them address the interrelationships among features.In this paper, to simultaneously address the two challenges and examine heterogeneous transition patterns, we propose to use Continuous-Time HMFM (CT-HMFM), integrating CTHMM, EFM, and a feature-based transition model. Our contributions are as follows: First, we examine the utility of data collected in an open environment from consumergrade mobile devices for mental health research. This contrasts with most existing studies on data collected in controlled lab environments. Second, we propose CT-HMFM to address the unique challenges introduced by mobile device data and depict the non-homogeneous state transition processes of multiple individuals. Simulation studies using synthetic data demonstrate exceptional parameter estimation and model selection performance. Finally, we analyze HRV and activity data from the AURORA study, followed by interpretations and discussions of biological findings that highlight the immense potential of mobile health data and our proposed method for mental health research.II. AURORA DATASET In the AURORA study [1] , data were collected from multiple sources, with our analysis focusing on the accelerometer and PPG data collected by the Verily smartwatches. Given the known associations with APNS [15] , [16] , the accelerometer data were used to quantify physical activity features, while the PPG data were pre-processed to extract heart rate variability (HRV) metrics (for details on the raw data and preprocessing steps, refer to [17] ). Briefly, activity features are extracted during a 24-hour window to evaluate daily activity patterns. After converting accelerometer data to activity counts [18] , meanAcc is the average activity counts. Amplitude, a cosinor rhythmometry feature, is computed to capture circadian rhythms [19] . Using the Cole-Kripke algorithm [20] , accelerometry epochs are classified into wake or sleep states, and the SWCK quantifies transition rates between wake and sleep. Additionally, the average activity during the five least active hours (L5) is calculated from raw accelerometer data [21] , representing nighttime activity.HRV features were derived from PPG data by first calculating and denoising beat-to-beat (BB) interval [22] time series to obtain normal-to-normal (NN) intervals, which were then analyzed using a 5-minute sliding window. Selected HRV features for this study include average heart rate (NNmean), skewness (NNskew), and standard deviation (SDNN) of NN intervals, indicating heart rate variability and rapid changes. Lfhf is the low-frequency to high-frequency power ratio, with extremes suggesting parasympathetic or sympathetic dominance [22] . DC serves as a mortality risk indicator in cardiac conditions [23] . The lower the DC index, the higher the mortality risk. Additionally, SD1SD2 and ApEn are used to quantify the unpredictability and regularity of successive heartbeats (R-R interval). To match the activity data, daily statistical summaries of each HRV feature are used, including mean, minimum, maximum, interquartile range, and variance.In this study, we focused on survivors of motor vehicle collision trauma. Since the data collection depends on the participants' wearing of the devices, missing data is a common issue. To assure data quality, we include only records that have complete activity data and a positive wake percentage. For HRV data, ideally, an individual can have 2,880 records per day. We retain only those days where at least 30% of these records are available to ensure that our daily summary statistics are representative. The final dataset consists of daily summaries of activity and HRV features from 258 patients, with each providing at least 50 days of records.Motivated by the structures of the processed AURORA datasets, we consider data in the form of repeated measurements of p features over T i occasions for each individual i of N subjects. The proposed models are in the framework of HMM. Let w it be the latent state of individual i at occasion t, taking value from the finite discrete set {1, • • • , J}.W i = (w i1 , • • • , w iTi )is the state sequence over T i repeated measurements. Let P it be the J × J transition probability matrix for individual i at occasion t, t = {2, • • • , T i }, of which the (k, j) entry is P it,kj = P (w it = j|w i,t-1 = k), and P it,kk = 1j:j̸ =k P it,kj . At t = 1, we assume that the initial state follows a multinomial distribution withprobabilities π = (π 1 , • • • , π J ) ′ , such that J i=1 π i = 1.The objective of the HMM is to delineate latent Markov processes given observations by estimating the transition probability matrix P and the initial state distribution. Distinct from traditional HMMs, our model further incorporates two additional components to meet our objectives. The first component, detailed in Section III-A, is a state-specific measurement model using EFM to uncover interrelationships among variables. The second component, outlined in Section III-B, is a transition model (TM) to capture heterogeneous transition patterns.Let y it denote a vector of the observed value of p features for subject i at time t. z it is a K dimensional vector of latent scores assumed to be independent of w it and following a standard multivariate normal distribution. While we assume that K is constant across states, our model can easily be extended to accommodate varying K j . For each individual i,Y i = (y i1 , • • • , y iTi ) is a p × T i matrix containing all measurements and Z i = (z i1 , • • • , z iT ) is a K × T imatrix containing all latent features. The first component of our model is an FM, with the primary goal of identifying the interrelationship structures between observed response variables and the underlying constructions of latent variables. For individual i at time t, given w it = j, the FM assumes that:EQUATIONwhere µ j is a p × 1 vector of state-specific expected mean response, Λ j is a p × K state-specific factor loading matrix, Ψ is a p × p diagonal covariance matrix for the error term e it with positive nonconstant diagonal entries. Alternatively, (1) can be expressed as[y it |w it = j]Taking into account the effects of the time interval, the continuous-time Markov process relies not on the transition probability matrix P , but on a transition intensity matrix Q [24] . This matrix Q is the limit of P as the time interval approaches zero. Suppose that δ it is the number of prespecified time units between t th and (t -1) th observation, then the transition intensity for subject i from state j to state k at time t isq jk = lim δit→0 P (w it = k|w i,t-1 = j) δ it > 0, j ̸ = k,and q jj = -k̸ =j q jk . The corresponding transition probability matrix P (δ it ) can be calculated as the matrix exponential of δ it * Q. The time intervals are assumed to be independent.To investigate the impact of covariates on transition rates, the transition intensity matrix can be modeled through a loglinear model [25] , such that log(q jk |x it ) = x ′ it B jk , where x it is a d × 1 vector of vector of covariates for individual i at time t, and B jk is a state-specific d × 1 vector of fixed effects coefficients. The B jk intends to quantify the effect of covariates on the probability of transitioning from state j to a different state k to provide an understanding of how covariates influence transition patterns and investigate the potential risk factors. Calculating the exponential of a matrix can be challenging. For computational efficiency, we approximate the exp(Q) using the (I + Q/a) a for some sufficiently large a [26] .IV. EXPECTATION-MAXIMIZATION ALGORITHM Let λ = ({µ j } J j=1 , {Λ j } J j=1 , Ψ, {B kj } J k,j=1 , π).Given the sequence of latent states W i and the latent factor scores Z i for each i, a joint probability distribution of the observations and all latent variables, L ci (λ), can be constructed as follows:P (w i1 ) × Ti t=2 P (w it |w i,t-1 , x it ) × Ti t=1 P (y it |w it , z it )P (z it ).(2) By the independence property of Y i , W i , and Z i across i, the complete likelihood function (L c ) for the whole sample can be obtained by taking the product of (2) over i.Our goal is to estimate λ by maximizing the likelihood function L c , or its logarithm l c . Since both W i and Z i are unobserved, the expectation-maximization (EM) algorithm is commonly used to identify the maximum likelihood estimator. As the name suggests, the EM algorithm finds a local maximum of the marginal likelihood by iteratively applying the expectation and maximization steps discussed below.The E-step gets the expectation of l c given observations, with respect to the current conditional distribution of unobserved variables and the current parameter estimates λ v . Denote the target expectation (i.e., E λ v [l c (λ)|Y , X]) as Ω(λ, λ v ). While an explicit form of the probability density function of z it exists, the calculation of conditional state probabilities can be computationally heavy. Therefore, we utilize a scaled version of the forward-backward algorithm (FBA) [27] to get the conditional state probabilities efficiently. Specifically, we first define the forward probability α ij (t) as P (w it = j|y i1 , • • • , y it ). Denote P j (y it ) the probability density function of y it given w it = j and c i (t) the conditional probability of observation y it given all past observations. Note that we omit the dependence of probabilities on δ it in notations for brevity. For each individual i and state j, using a recursion scheme, the forward probabilities at t = 1, • • • , T i will be calculated as:EQUATION;α ij (t) = P j (y it )[ J k=1 α ik (t -1)P itkj ] c i (t),where c i (t) = J j=1 P j (y it )[ J k=1 α ik (t -1)P itkj ]. Then, we define the backward probability β ij (t) asP (y i,t+1 ,••• ,y i,T i |wit=j,λ) ci(t+1). Similarly, we define a recursion form to update the backward probabilities at t = T i , • • • , 1:β ij (T i ) = 1; β ij (t) = J k=1 P i,t+1,jk P k (y i,t+1 )β ik (t + 1) c i (t + 1) .After that, in the smoothing step, denote ϵ v ikj (t) as P (w i,t = j, w i,t-1 = k|Y i , λ v ) and γ v ij (t) as P (w it = j|Y i , λ v ). The target conditional state probabilities are functions of the forward probability and backward probability as follows:γ v ij (t) = α ij (t)β ij (t); ϵ v ikj (t) = α v ik (t -1)P itkj P j (y it )β v ij (t) c i (t).Then, the Ω(λ, λ v ) can be written asconstant+h(π)+h({B kj } J k,j=1 )- 1 2 h(Ψ, {Λ j } J j=1 , {µ j } J j=1 ),where h(π) depends on the initial state distribution, h({B kj } J k,j=1 ) depends on the probability transition matrix, and h(Ψ, {Λ j } J j=1 , {µ j } J j=1 ) is a function of parameters Ψ, Λ j , and µ j . Explicit forms are provided in Appendix A-A.Within each M-step, since h(Ψ, {Λ j }, {µ j }), h(π), and h({B kj }) do not share parameters, we maximize each of them separately. The estimator of π, Λ j , µ j , and Ψ can be directly derived by setting h(π) = 0 and h(Ψ, {Λ j }, {µ j }) = 0 (see Appendix A-B for details). For {B kj } J k,j=1 , a one-step Fisher scoring (FS) [28] is implemented.Let θ be a vector of all transition model parameters, such that θ = vec({B ′ kj }, k ̸ = j). Recalling the first derivative of matrix exponential [14] and using Theorem 1 in [29] ,∂ ∂θ u exp(A(θ u )) = exp( A(θ u ) Ã(θ u ) 0 A(θ u ) ) 0:J,J:2J , where Ã(θ u ) = ( Ãij (θ u )) = ( ∂Aij (θu) ∂θu). Denote ∂P kj (δit) ∂θu as the (k, j) entry of the first derivative of P (δ it ) with respect to θ u (i.e., the u th entry of θ). Having the first derivative of P (δ it ) = exp(δ it * Q) with respect to each component of θ calculated accordingly, the FS can be directly implemented to update θ to forbid the calculation of the second derivative of matrix exponential. Specifically, denote S * be the score function,S * u (θ) = N i=1 Ti t=2 J j=1 J k=1 ϵ v ikj (t) P kj (δ it ) ∂P kj (δ it ) ∂θ u .Let M * be the negative Fisher information matrix. Its (u, v) entry M * uv is in the form of:M * uv (θ) = N i=1 Ti t=2 J j=1 J k=1 γ v ik (t -1) P kj (δ it ) ∂P kj (δ it ) ∂θ u ∂P kj (δ it ) ∂θ v .After getting both the score function and the Fisher information matrix, parameters θ can be updated asθ v+1 = θ v + M * (θ v ) -1 S * (θ v ).To ensure stability, we control the learning rate by updatingθ v+1 = θ v + {M * (θ v ) + S * (θ v ) T S * (θ v )} -1 S * (θ v ) in practice.Note that the algorithm requires the specification of (K, J), which are typically unknown in practice. In this study, we propose to determine (K, J) using information criteria, the efficacy of which is evaluated in Section V-C.V. SIMULATION STUDY This section assesses the proposed methods through simulation studies using synthetic data resembling the AURORA dataset. We generate data with N=200, p=23, d=3, J = 3, and K = 3. Each individual's number of observations, T i , is uniformly sampled from [50, 100]. From this, we randomly select T i time points from {1, • • • , 100} to create sequences of δ it . Each individual's initial state is drawn from a multinomial distribution with probabilities π = ( 13 , 1 3 , 1 3 ). Latent state trajectories are then generated based on transition probabilities P it (δ it ), given individuals' features. Observation vectors y it for an individual i in state j at time t are then sampled from a normal distribution with mean µ j and covariance Λ j Λ ′ j + Ψ, where Ψ = I. In the following, subsection V-A assesses model reliability by comparing empirical parameter estimates against true values; subsection V-B compares the effectiveness of our method against baseline methods; and subsection V-C explores the performance of information criteria in model selection.To validate the estimation procedure, we implement the EM algorithm with true J and K. Parameter initialization involves first fitting Gaussian Mixture Models to estimate groups, followed by EFM for each group. Guided by the insights from a pilot study, we set the maximum number of iterations for each replication at 100. The reliability and precision of the proposed methods are then evaluated from two perspectives: i) the accuracy of each individual parameter estimate and ii) the misclassification rate (C mis ), which quantifies the proportion of estimated states that diverge from the actual states.The accuracy of parameters π, µ, Λ, and Ψ is assessed by calculating the average absolute difference (AAD) between the estimates and their true values, defined asAAD(o) = r i=1 | ôi-oi| r, o i is an individual entry in matrix o and r is the total number of free parameters. The mean of AADs (standard errors in the parentheses) aggregated over 100 random seeds are presented in Table I . These mean AADs for all parameter matrices are sufficiently close to zero with small standard errors, indicating effective parameter recovery. In Table II , we report the mean bias (standard error) for each parameter in the transition model, where the biases are all close to zero. Moreover, we present the mean (standard error) of C mis in Table I . On average, misclassification rates are only .24% (0.0005), highlighting the exceptional accuracy of the proposed EM algorithm in estimating latent states.Intuitively, factors such as N, T i , sizes of J and K, variance Ψ, differences in µ j and Λ j between states, and frequency of state transitions, can affect the performance of parameter estimation. Additional simulations in Appendix B reveals that estimation accuracy for µ, Λ, Ψ, and B, as well as classification accuracy, improve when (i) common variances decrease, (ii) differences in µ j and Λ j between states increase, (iii) J decreases, or (iv) sample size (N ) or panel length (T i ) increases. Increasing K or using a B that induces infrequent transitions slightly affects most parameters but enhances the precision of transition probability estimates, thereby reducing misclassification rates. Estimation of π improves solely with increases in N or state-to-state differences in µ j .This section compares the performance of the proposed methods and the baseline approaches in correctly identifying latent states. Three benchmark methods are under our consideration: i) TM+independent HMM, which assumes independence among observed features given the states; ii) CFM+TM+HMM, which addresses interrelationships but inaccurately pre-specifies the latent factor structure by setting certain loading matrix entries to zero; and iii) EFM+HMM, which assumes a homogeneous transition probabilities matrix for all subjects. We first repeat the data generation process of Simulation 1. Then, we consider three additional scenarios by adjusting the state-to-state differences in µ j to be closer (µ: medium diff), increasing the similarity of the Λ j at different states (Λ: medium diff), and increasing the significance of the covariance matrix Ψ (Ψ = 2 × I), respectively.As depicted in Figure 1 , our proposed methods (HMFM) consistently outperform the benchmark methodologies in all settings. Regardless of sample size, our methods consistently achieve the lowest misclassification rate, nearly approximating zero, thereby emphasizing the importance of each component in our proposed models. Specifically, the comparison with TM+independent HMM shows the importance of accounting for the interrelationship between observed features; the comparison with CFM+TM+HMM reveals the risk of incorrectly specifying the interrelationship structure; and the comparison with EFM+HMM demonstrates the inadequacy of assuming homogeneous transition probabilities.Information criteria such as the Akaike information criteria (AIC) and the Bayesian information criteria (BIC) have been widely used in model selection [10] , [30] . Within this simulation study, we investigate whether the AIC or BIC is reliable for determining J and K simultaneously. We repeat the data generation process of Simulation 1, but implement the proposed methods with a different set of (J, K) for each replicate when fitting the generated data. Let J = {2, 3, 4} and K = {2, 3, 4}. We consider all possible combinations of J and K, yielding a total of nine fitted candidate models for each replicate. For 100 replications, both AIC and BIC consistently recommend the model with accurate J and K. Therefore, we believe that the sample size and the number of observations per individual in the processed AURORA data will yield reliable information criteria-based model selection results and, consequently, reliable parameter estimation. VI. ANALYSIS OF THE AURORA DATA Given the irregular data collection from mobile devices, we applied the CT-HMFM to smartwatch data from the AURORA study. Considering J = 1, 2, • • • , 6 and K = 1, 2, • • • , 9, we evaluate 54 candidate models. For each candidate model, the EM algorithm is implemented with multiple random seeds, and the seed yielding the highest estimated likelihood is selected. Finally, model comparison using AIC and BIC led to the selection of a model with J = 3 and K = 8. Subsequent subsections will discuss the interpretation of parameter estimates and biological insights from three perspectives: i) the interpretation of the three estimated states, ii) symptom cooccurrence patterns, and iii) the influence of demographic factors on transition probabilities.To investigate the biological differences between states, we first focus on the selected features. Figure 2 depicts scaled sample means for each feature across states, with a 99% confidence interval (CI). Further pairwise Tukey tests indicate significant differences between states for nearly all features, except for amplitude, SWCK, L5, and NNskew.q3 between Fig. 2 . Relative sample mean for features in each estimated state, with 99% CI error bars that are too small to distinguish. Fig. 3 . Sample mean for each symptom in each estimated state. The error bars represent the 95% CI. 0 is the least severity, 1 is the greatest. state 1 and state 2. Specifically, features related to average heart rate (NNmean), heart rate variability (SDNN), and heart deceleration capacity (dc) vary significantly across the three latent states, decreasing sequentially from state 1 to state 3. According to previous research, lower heart rate variability and deceleration capacity are associated with a higher mortality rate [23] , [31] , suggesting that states 1 through 3 represent decreasing levels of health, with state 1 being the healthiest. Additionally, activity levels are similar and higher in states 1 and 2 compared to state 3, indicating better overall health for participants in the first two states. Among features related to heart rate unpredictability (lfhf, ApEn, and SD1SD2), state 3 demonstrates significantly higher values for SD1SD2-related features but lower values for lfhf-related features compared to states 1 and 2, suggesting a different interpretation of the estimated states than our previous interpretation. However, it is important to note that the relationship between these features and the psychological or physiological state is neither straightforward nor unique [32] .To confirm the validity of the three states, we further compare their differences regarding self-reported symptoms from a flash survey. Based on the RDoC framework, ten latent constructs associated with APNS were developed using survey items selected by domain experts: Pain, Loss, Sleep Discontinuity, Nightmare, Somatic Symptoms, Mental Fatigue, Avoidance, Re-experience, and Anxious. Retaining only observations for each individual whose estimated states are known on the same day they submitted survey responses, we summarized the flash survey data with means and 95% CIs in Figure 3 . Overall, state 1 exhibits the lowest severity level for all ten symptoms, while state 3 has the highest severity level. Tukey tests reveal no significant differences between states 1 and 2 in hyperarousal, re-experience, anxiety, and somatic symptoms, but both are significantly different from state 3 in these constructs. While the differences in nightmare and sleep discontinuity between states 3 and 2 are not significant, they are statistically more severe than in state 1. For mental fatigue and depression, only the difference between state 1 and state 3 is statistically significant. In summary, both the survey data and the AURORA data support our interpretation of the three latent states. State 1 is the healthiest, while state 3 indicates having the most severe APNS symptoms.When studying the co-occurring of symptoms within each hidden state, we limit our attention to observations collected during the first week. For each state, the correlations between all ten symptoms are calculated. In state 1 (relative health state), there is a high degree of correlation (.782) between hyperarousal and anxiety, suggesting that patients in state 1 experiencing severe hyperarousal symptoms are also likely to suffer from severe anxiety. State 2 shows no highly correlated symptoms. In state 3 (the state with more severe disorders), symptoms such as depression, hyperarousal, anxiety, and reexperience are more likely to co-occur with pairwise correlation ranging from .716 to .91.This section investigates the heterogeneity of 1-day transition probabilities among subjects, focusing on transitions with a time interval δ it = 1. We estimated the transition probabilities for males and females within the sample age range, as depicted in Figure 4 . Lines with circles illustrate the probability of remaining in the same state, lines with stars indicate transitions to a more severe state, and lines with 'x' reflect the chance of improvement in psychological conditions.Overall, both males and females have a tendency to remain in their current state, with infrequent state transitions, aligning with most literature. For males, the probability of staying at states 3 and 2 increases with age, while it decreases for state 1. Moreover, while the likelihood of psychological deterioration increases with age, the chance of improvement decreases. Specifically, while the probability of transitioning from the most severe state (state 3) to the healthiest state (state 1) approaches zero as age increases, the likelihood of the reverse transition increases as age decreases, with direct transitions between state 1 and state 3 being particularly rare. The female group exhibits a similar trend to the male group, but with a higher likelihood of remaining in the most severe state (state 3) compared to males. In summary, our analysis of the AURORA data suggests that older patients are more likely to transition to more severe psychological states. Moreover, achieving psychological improvement becomes increasingly challenging as one ages.This paper investigates the unique challenges of analyzing longitudinal mobile health data, including interdependent variables with unknown interrelationship structures, heterogeneous transition probabilities, and irregular measurements. To address these issues, we propose a HMM-based model, the CT-HMFM, for multivariate longitudinal data collected irregularly. Furthermore, the performance of the corresponding Stabilized Expectation-Maximization algorithm for maximum likelihood estimation is supported by extensive simulation studies. Finally, we analyzed the AURORA data and drew biological findings comparable with previous research, implying that the mobile health data sourced from consumergrade devices, together with the proposed methods, have the immense potential to facilitate mental health diagnostics and understand the dynamic transition mechanism.The proposed methods can be extended in several ways. First, most entries in the estimated factor loading matrix are close to zero, indicating sparse factor loading matrices in real analysis. Although various methods (e.g., factor rotations and setting factor loadings below specific cutoffs to 0) are frequently used to simplify interpretation, the choice of these methods is subjective. The sparse exploratory factor loading analysis [33] , [34] provides an automated approach to set the loading entries of redundant variables to 0, thereby enhancing the interpretability of loading matrices without reliance on subjective factors. Therefore, incorporating sparse regularization into the factor loading matrix is an important extension of our current work worth studying. Second, a large number of baseline covariates are typically available in real data. However, we have no prior knowledge about the significance of each covariate in determining the transition probability. Hence, integrating regularization into the transition model to assist with variable selection can be extremely useful. Third, the exceptional diversity of mental health makes it challenging to satisfy the key conditional independence assumption of the HMM, leading to potential violations and hence model bias [35] . This issue is exacerbated by the likely presence of autocorrelation among observations collected from the same subject. Therefore, adding a random effect to the current model to account for the inter-patient heterogeneity is a natural extension [10] . Finally, previous HRV-related studies are often conducted in well-controlled laboratory environments. Thus, all existing HRV feature extraction tools rely on restingstate heart rate data. However, heart rate data collected in open environments will inevitably contain additional noise. For example, it is reasonable to expect that HRV features corresponding to different activity states (e.g., exercising and resting) would differ significantly. Therefore, recognizing the lack of tools to extract HRV features corresponding to different activity states, we believe it would be advantageous to develop a preprocessing pipeline to concurrently process heart rate and activity data to derive appropriate HRV features.Denote Λ j = (Λ j , µ j ) ∈ R p×(K+1) and z it = (z T it , 1) T ∈ R (K+1) . Each of the three parts has an explicit form: EQUATION+ γ v ij (t)y ′ it Ψ -1 y it -2γ v ij (t)y ′ it Ψ -1 Λ j E λ v ( z it |y it , w it ) + γ v ij (t)tr( Λ j ′ Ψ -1 Λ j E λ v ( z it z it ′|y it , w it )).By setting the first derivative of h(π) to zero, the parameters related to the initial state distribution are estimated as: (1) . Similarly, the parameters used to characterize the conditional distribution of y it given w it are estimated by setting the first derivative of h(Ψ, {Λ j } J j=1 , {µ j } J j=1 ) equal to 0, with π new j = N i=1 γ v ij (1)/ N i=1 J k=1 γ v ikThis section presents additional simulation results that investigate the impact of various factors on estimation performance. Using the settings for simulation 1 in Section V-A as the baseline, we conducted eight additional sets of simulations, each varying one component while maintaining rest components as the baseline setup. These components include: i) sample size (N), ii) number of measurements for each individual (T i ), iii) J, iv) K, v) size of common variance Ψ, vi) state-to-state difference in µ j , vii) state-tostate difference in Λ j , and viii) transition frequency. For common variance Ψ, we evaluated scenarios with .1I, .5I, and 1I (baseline). For µ j and Λ j , we adjusted state-to-state differences in two additional settings. For the test evaluating the effect of transition frequency, a frequent transition is defined as the probability of remaining in the same state being less than 0.70. The B in the baseline setting corresponds to infrequent transition. For tests evaluating the effects of J, K,i.i.d. ∼ N (µ j , Λ j Λ ′ j + Ψ).It is crucial to emphasize that, unlike CFM with pre-specified structures of factor loading matrices, our approach imposes no assumptions on Λ j . Therefore, the structure of Λ j is entirely data-driven, making the first component of (1) an EFM.B. Transition ModelGiven a state sequence W i , standard assumptions of HMM assume that 1) given a state w it , observations y it are independent, and 2) given a state w it and subjects' contextual features, the state at the subsequent occasion w i,t+1 is unrelated to any information from previous occasions. 
 
 Appendix:  



Reviews: review: title: title: Innovative improvement on HMM model applied on PNAS
paper_summary: paper_summary: In this paper, the authors proposed using a new method of HMM, named CTHMFM, for analysis of Adverse posttraumatic neuropsychiatric sequelae (APNS). The model is improved upon the base HMM with the addition of the State-Specific Measurement Model (which models the correlations among variables) and the a transition model (TM) to capture heterogeneous transition patterns. This CTHMFM is careful described and verified using simulated tasks, and then subsequently applied on the AURORA dataset. The resulting interpretable model is used to generate clinical insights.
strength_weakness: strength_weakness: 1. Lack of a code repo for the model, which could have other strong use cases.
2. Lack of a direct classification verification for PNAS.
questions: questions: I am confused about why the authors used simulation studies to prove the model's usefulness in classification but did not conduct the same classification task on the AURORA dataset. I may have misunderstood something, but the authors' clarification would greatly help.
title: title: Review of Continuous-Time Hidden Markov Factor Model for Mobile Health Data: Application to Adverse Posttraumatic Neuropsychiatric Sequelae
paper_summary: paper_summary: The paper presents a Continuous-Time Hidden Markov Factor Model (CT-HMFM) for analyzing longitudinal mobile health data. This model addresses the challenges of interdependent variables and irregularly spaced measurements. It combines Continuous-Time Hidden Markov Models (CTHMM) with Exploratory Factor Analysis to identify latent states and capture transitions between these states as a function of covariates in a continuous-time framework. The paper also describes an Expectation-Maximization (EM) algorithm for parameter estimation. Extensive simulation studies validate the accuracy and performance of the model, showing that it outperforms benchmark models. The CT-HMFM is applied to data from the AURORA study, focusing on heart rate variability (HRV) and activity data. It provides insights into the dynamic transitions among different health states in trauma survivors and the influence of demographic factors on these transitions. The model is presented as a robust tool for studying Adverse Posttraumatic Neuropsychiatric Sequelae (APNS) and has the potential for broader applications in mental health research.
strength_weakness: strength_weakness: The model's complexity, especially integrating CTHMM with Exploratory Factor Analysis, may make it challenging to apply in real-world settings outside of academic research. Practitioners who could benefit from this model might find it challenging to implement due to the advanced mathematical and computational requirements. The model relies on several assumptions, such as the independence of observations given the latent states and the structure of the factor model. If these assumptions are violated in practice, the model's performance could be compromised. The paper mentions the inclusion of covariates in the transition model, but the exploration of which covariates are most influential or how they impact the results is limited.
questions: questions: 1. Can you please clarify if there are specific tools or software packages that you have developed or recommended for implementing the CT-HMFM model? Considering its complexity, how do practitioners or researchers apply this model in real-world settings?

2. How sensitive is the CT-HMFM model to the key assumptions, such as the independence of observations given latent states? Have you conducted any sensitivity analyses or explored extensions that might relax some of these assumptions?

3. Can you provide more details on the role of covariates in the transition model? How did you select the covariates used in your analysis, and what insights did you gain from them? Are there any covariates you found particularly influential or surprising in their impact?
title: title: Continuous-Time Hidden Markov Factor Model for Mobile Health Data: Application to Adverse Posttraumatic Neuropsychiatric Sequelae
paper_summary: paper_summary: This paper develops a Continuous-Time Hidden Markov Factor Model (CT-HMFM) and a Stabilized Expectation-Maximization algorithm for parameter estimation to study Adverse Posttraumatic Neuropsychiatric Sequelae (APNS), a mental health disease common in veterans and trauma survivors, by leveraging objective mobile device data to improve understanding of these conditions. Traditional research methods rely heavily on subjective self-reports, leading to measurement errors and biases. The authors propose a novel approach using longitudinal data from mobile devices in AURORA study to identify and track homogeneous APNS states and their transitions over time.
strength_weakness: strength_weakness: The dataset section is not described clearly. Only the number of subjects is included in the description without the wearables data description.
questions: questions: Why do the patients are included specifically on Day 67? Is 50 days a special threshold for APNS patients?

Could you please provide more details about the AURORA dataset? What is the sampling frequency of signals? Does 258 subjects wear the same type of device?