A Within-Subject Normal-Mixture Model with Mixed-Effects for Analyzing Heart Rate Variability

Data on Heart Rate Variability (HRV) have been used extensively to indirectly assess the autonomic control of the heart. The distributions of HRV measures, such as the RR-interval, are not necessarily normally distributed and current methodology does not typically incorporate this characteristic. In this article, a mixed-effects modeling approach under the assumption of a two-component normal-mixture distribution for the within-subject observations has been proposed. Estimation of the parameters of the model was performed through an application of the EM algorithm, which is different from the traditional EM application for the normal-mixture methods. An application of this method was illustrated and the results from a simulation study were discussed. Differences among other methods were also reviewed.


Introduction
Heart rate is a standard cardiac measure that is non-invasive, inexpensive, and relatively simple to acquire. The variations in heart rate are of great interest to physicians and researchers because they capture valuable information from both the sympathetic and parasympathetic branches of the autonomic nervous system [1]. High variability in heart rate has been associated with signs of good adaptability, indicating a healthy individual with well-functioning autonomic control, whereas lower variability typically suggests abnormal or insufficient adaptability, indicative of some physiological malfunction [2]. One standard measure used to assess heart rate variability (HRV) is the RR-interval. The RR-interval quantifies the variations in heart rate by measuring the time, typically in milliseconds, between successive "R" peaks in a sinus rhythm, as shown in Figure 1.
A great deal of progress has been made in the field of HRV; however, it has mainly focused on measuring HRV, with little emphasis on modeling the measures of HRV [3]. The statistical measures used, such as the mean RR-interval and the standard deviation of the beat-tobeat RR-intervals, are primarily descriptive and may not characterize and utilize the true nature of the distribution of the data. The focus of this article is in developing methods to model the RR-interval lengths and other data that have similar within-subject characteristics. Several studies have reported that RR-interval data exhibit a two-component normal-mixture distribution in both healthy and diseased populations and that the usual assumption of normality is often inadequate [4][5][6][7][8][9]. These studies also provide some scientific rationale for the normalmixture distribution. For instance, the respiratory sinus arrhythmia (RSA) could explain the two-component normal-mixture distributions; the heart rate has been shown to increase and decrease during the inspiration and expiration phases of a breathing cycle, respectively [10,11].
Developing a mixed-effects model that could appropriately handle the mixture distributions of the RR-intervals was primarily motivated by three examples. The first example is a study involving a sample of adult, sedated and ventilated inpatients recorded for periods of up to 24 hours [12]. The primary aim was to characterize the effect of level of sedation on a variety of physiological and comfort measures, including RR-intervals. The second example is a sample of preterm infants recorded during several bottle feeding sessions, lasting up to 15 minutes each [13]. The goal was to predict bottle feeding readiness, in part, by understanding the effect of feeding on HRV. The third example is a study examining undergraduate college students during baseline, preparation, and delivery of several psychological stressor tasks [14]. In all of these examples, the data comprised of repeated measures of RR-intervals within each subject along with a variety of additional within and between subject characteristics. A common theme in these studies was the desire to identify the factors influencing the variations in RR-interval lengths. The methods proposed in this article focus on modeling these types of data under a within-subject normal-mixture assumption with mixed-effects modeling strategies. The method will be illustrated using the psychological stressor data [14].
In the literature, there is extensive research for mixed-effects models with the assumption of a normal-mixture distribution for the random between-subject effects. These do not capture the appropriate mixture of the normal-distributions for the within-subject effects relevant

Abstract
Data on Heart Rate Variability (HRV) have been used extensively to indirectly assess the autonomic control of the heart. The distributions of HRV measures, such as the RR-interval, are not necessarily normally distributed and current methodology does not typically incorporate this characteristic. In this article, a mixed-effects modeling approach under the assumption of a two-component normal-mixture distribution for the within-subject observations has been proposed. Estimation of the parameters of the model was performed through an application of the EM algorithm, which is different from the traditional EM application for the normal-mixture methods. An application of this method was illustrated and the results from a simulation study were discussed. Differences among other methods were also reviewed. to the examples discussed above. For instance, Xu and Hedeker [15] and Verbeke and Molenberghs [16] propose a mixed-effects model, where the random-effects are assumed to be normal-mixture while the residuals are assumed to be normally distributed. Moreover, in both of these approaches, the mixture components of the normal-mixture distributions are assumed to have the same variance-covariance matrices. In this article, a new formulation is proposed. This method is specific to within-subject normal-mixture distributions and allows for the mixture components to each have distinct variance-covariance matrices. When the mixture component variance-covariance matrices are assumed to be equal, the joint distributions in the new formulation is equivalent to the previous formulations [15,16], however, the posterior distributions are entirely different. Consequently, hypothesis tests regarding the random-effects are also different. As will be shown, this formulation proposed in this article also lends itself to a straightforward application of the EM algorithm, which eliminates the need for numerical integrations for obtaining the marginal likelihood. This approach also differs from Ng et al. [17] who consider normalmixture for the marginal distributions in mixed-effects models and apply the EM algorithm. This difference will be further discussed in a subsequent section.

Within-subjects normal-mixture mixed-effects model
For simplicity, the model described here is for a two-component normal-mixture model with mixed-effects (MXME): where Y is an N×1 vector of observed responses, β is the q×1 vector of random-effects, Z is the N×q observed design matrix corresponding to the random-effects, α 1 and α 2 are the p×1 vectors of fixed-effects parameters for the two components, X is the known N×p observed design matrix corresponding to p-1 fixed-effects and an intercept, and ε is the N×1 vector of residuals. Here, N is the total number of observations from s subjects. The vector of random-effects, β, is assumed to follow a q-variate normal distribution with mean vector 0 and variance-covariance matrix G. The error vector, ε, is assumed to follow a joint normal-mixture density with two N×1 mean vectors, 0 and 0, two N×N variance-covariance matrices where ( ) ; , represents the n j dimensional multivariate normal density for k = 1, 2. Here, X j is an n j ×p dimensional design matrix with identical rows and kj Σ is the corresponding n j ×n j variancecovariance matrix for the k th component (k = 1, 2).

Estimation of the model parameters
For estimating the parameters, the marginal likelihood is given by the density of Y. This could be derived in the usual way by obtaining the joint density of Y and β, and integrating it with respect to β. That is, However, this does not lead to a closed form expression and involves the product of indefinite integrals. To obtain the maximum likelihood estimates (MLEs), (3) has to be maximized with respect to the various parameters. This requires the MLEs to be obtained iteratively through a numerical algorithm that includes numerical integrations. Alternatively, an application of the EM algorithm eliminates the need for a numerical approach.
The maximum likelihood estimation of parameters and the use of the EM algorithm for the normal-mixture models have been extensively studied in the literature [17][18][19]. The EM algorithm is adapted in the current context in such a way that the existing software for fitting mixed-effects models (such as PROC MIXED or PROC GLIMMIX in SAS) could be used. This makes our approach easy to implement. In normal-mixture models, there is a natural choice of the "incomplete" data. These are defined by an indicator variable, c ij , which designate each observation into one of the two components of the mixture normal with a certain probability. That is, th th 1, if the observation for the subject is from the first component 0, otherwise In the E-step, the missing data is represented by the probability of the i th observation from the j th subject falling into the first component. That is, ; , To implement the M-step, most of the current methods (e.g., [17]) divide the data into components of the mixture normal by strictly assigning an observation to the component that has the maximum expected probability. In the approach proposed here, rather than separating the data into distinct components, the probabilities are assigned to each observation weighting the observation appropriately into each component. That is, weights, w ijk , are assigned as In the M-step, the likelihoods for the weighted observations are then used, along with the distribution of the random-effects, to obtain the overall MLEs for the joint MXME using standard mixed-effects , across all j and k, the mixed-effects model is fitted to the weighted data, * The weighted residuals, * ε , follow an N-dimensional multivariate normal distribution, with mean 0 and variance-covariance matrix Here, one could use any standard software (e.g., PROC MIXED or PROC GLIMMIX in SAS) to fit the models. However, for each component an estimate of the mixed-effects covariance matrix, G, is produced. A consistent estimator for G is based on the weighted average of the predicted random effects. That is, where ( )

Summary of the EM algorithm for the two-component MXME
For the first (t = 1) iteration: • First M-step: The MLEs for the component means, ( ) 1 k α , and component variances, Σ , are obtained by fitting the two mixed-effects models separately. An estimate for the overall random-effects variance-covariance is obtained by substituting the predicted random-effects into (5). The updated estimate of the mixture proportion from the j th subject is estimated by (1) (1) 1 For the (t+1) th iteration: Λ α obtained in the t th M-step. Then determine the component weight matrices, . Compute the two sets of weighted observations, *( 1) , for k = 1, 2.
• M-step: Same as the M-step described above using the (t+1) th iteration predicted values of the random-effects.
The process is iterated until the desired level of convergence is achieved.
The model fitting could proceed in the usual way using likelihood ratio tests (for nested models) or using information criteria (such as the AIC, BIC, etc.). Whether or not a mixture model is necessary could be examined only by comparing the information criteria because the normal model is not nested within the mixture model. In this article the methods are described for a two-component mixture, but the extensions to three or more components is straight forward. However, the number of components has to be fixed (assumed to be known). In application, when the number of components is unknown, theory based research and preliminary exploratory analysis of the data (within each subject) should be assessed for the best choice for the fixed number of components.
Comparison to other methods using mixture normality in mixed-effects models As pointed out in the introduction, in order to capture the withinsubject mixture normality, the formulation in equation (1) is necessary. After some algebra, under the formulation in (1), it can be shown that the posterior mean for the jth individual is: The between-subject normal-mixture model proposed by Verbeke and Molenberghs [16], (equation (1)), assumes the random-effects to follow a normal-mixture distribution. Our formulation differs from this in three aspects. First, while the joint distributions (the integrand in equation (3)) have similar configuration in the two formulations, the posterior mean is different from the equation above. The posterior means in Verbeke and Molenberghs [16] and the one under model (1) (shown above) would be equivalent only when the normalmixture distribution is the special case of the normal distribution. In other words, if and only if,  [16]). Second, the between-subjects normal-mixture model does not lend itself to the application of the EM algorithm. Therefore computation of the marginal likelihoods using numerical methods becomes inevitable. Finally, the previous formulations assume the normal-mixture components to have the same variance covariance matrices. The formulation proposed in (1) does not require this assumption, which provides greater flexibility.
Ng et al. [17] also proposed a mixed-effects model with a normalmixture distribution and applied the EM algorithm for obtaining the MLEs of the model parameters. There are two main differences between their approach and the one proposed here. Their model formulation is theoretically different from our formulation, as well as the two mentioned above. By assuming the marginal distribution to be a normal-mixture, the formulation neither fits into the normal-mixture for the between-subjects effects nor for the within-subject observations. In [15] and [20] or our formulations, the marginal distributions do not have a closed form solution. Therefore, Ng et al.'s model [17] is not directly applicable to the example presented here. Moreover, they assume independence for the residuals, which is somewhat restrictive. The second difference is with respect to the application of the EM algorithm. Their approach is to estimate the probabilities of an observation belonging to the various components and assign them to the component that has the largest probability. This could lead to an ambiguity. For instance, if the two component probabilities for a given observation are 0.51 and 0.49 assigning that observation to either component seems inappropriate. However, their approach will assign the observation strictly to the component with probability 0.51. In our method this ambiguity is avoided, because instead of allocating the observation "outright" to one of the components, we use the probabilities as weights, thereby using all of the information available.

Description of the physiological stressor task data
In this section, the method proposed is applied to heart rate (RRinterval) data obtained from adults during physiological stressor tasks. The RR-interval data were collected from nine healthy undergraduate During the study, several cardiovascular and psychosocial measures were recorded. The ECG was consistently recorded while each subject performed one of six psychosocial stressor tasks. The psychosocial stressor tasks include five speech related tasks and one verbal mental arithmetic task, assigned in a random order to each subject. The mental arithmetic task (MA) had four minutes of baseline and four minutes of task. The "Saab" (BS) speech task had four minutes of baseline, two minutes of preparation, and two minutes of delivery. The remaining four speech tasks, "why I'm likeable" (LS), "ask for date" (AS), "describe way to school" (WS), and "describe inanimate objects in room" (IS), each had two minutes of baseline, two minutes of preparation, and two minutes of delivery. A total of 32 minutes of ECG data from each subject was available for the data analysis. Each analog ECG signal, for every minute, and from every subject, was first acquired at 1000 Hz then decimated to 500 Hz. Next, the ECG data was carefully edited and checked for any artifacts. The R-peaks were identified using waveform matching templates and then a time/amplitude criterion. Once the R-peaks had been identified, the RR-interval series was calculated using the time between successive R-peaks. For further details see Mandrekar [9]. In the original analysis [9], each subject's data were analyzed separately, mainly to illustrate the inadequacy of a normal model and to suggest that a normal-mixture model might be more appropriate. Here, we use these all of the subject's data to fit a single model that accounts for the repeated measures on each subject and simultaneously include fixed effects (age, gender, and the type of task).
Prior to fitting the model, several data cleaning steps were administered. A total of 588 RR-intervals (out of 20,416), were identified as artifacts and subsequently excluded from the analysis due to missed or spurious beats. An interval was considered an artifact if it was "large" or "small" in relation to all the other intervals for that subject during the specified task and period combination, and in relation to the intervals surrounding the artifact. On average 2.88% of the data from each subject was identified as artifacts. After excluding artifacts, the data were de-trended in order to remove the ultra-low frequency trends and to reduce some of the non-stationarity, as suggested by Litvack et al. [11]. This was accomplished by fitting first-order polynomials to the series of RR-intervals obtained from each subject, for each combination of task and period. The de-trended, artifact free data (rescaled by a factor of 10) was used in the subsequent analysis.

Estimation of the model
The EM algorithm described previously was implemented with a SAS macro using the PROC IML and PROC MIXED in SAS v.9. The PROC IML code was primarily used for the E-Step and the PROC MIXED was called repeatedly within this code for the EM iterations. In this application, the within-subject residuals were assumed to be independent. The degrees of freedom in MIXED procedure were appropriately adjusted using the DDFM statement. The specific hypotheses tested are representative of typical research questions that may be of interest: 1. Do the demographic variables have significant effects on RRintervals?
2. Is the mixture proportion different across the subjects?
3. Do task and period have significant effects on heart rate? If so, is the period effect consistent across the tasks?
4. Is there significant subject to subject variability in the heart rate?
5. Are the estimates of the within-subject variance-covariances different across the components?
The estimated mixture proportions ranged from 0.3809 to 0.5628 across the nine subjects, with standard errors (SEs) ranging from 0.0082 to 0.1022. To test if there were differences among mixture proportions across the subjects, a model with a single mixture proportion for all nine subjects was fitted and the corresponding estimate of the overall mixture proportion was ˆ0 .4796 , p-value = 0.96). This test was performed under a full fixed-effects model. For subsequent analyses, the models were fit with a single mixture proportion. Table 1 shows the likelihood ratio tests for the fixed-effects.
The estimates and standard errors for the variance components along with 95% Wald confidence intervals are given in Table 2. The first two rows correspond to the residual variances of the first and second components and the last row represents the variance for the random subject effects (i.e. the between subject variability).
The least squares (LS) means and standard errors for the combinations of task and period effects were computed at the mean age (19) for each component. The component LS means are shown in Table  3.The LS means for the normal-mixture distribution (combining both components) could be computed as , where gender is denoted by f, task is denoted by u, and period is denoted by v. Figure 2, at an average age of 19.22, are plotted separately for males (panel 3a and 3c) and females (panels 3b and 3d) across each period, baseline (B), preparation (P), and delivery (D), with a separate line for each of the five tasks.

The component LS means, shown in
The main results and conclusions are as follows: • A normal-mixture model with two-components with the same mixture proportion across the subjects fits the data adequately. As suggested in the Introduction, the means of the two components might represent the mean increase and decrease in the heart rate during the inspiration and expiration phased of a breathing cycle due to respiratory sinus arrhythmia, respectively. The cardiac parasympathetic nerve is minimized during inspiration, leading to shorter RR-interval lengths, while the cardiac parasympathetic nerve is maximized during expiration, leading to longer RRintervals [10].
• The heart rate measures vary significantly across the demographic variables, age and gender. While this is consistent with previous research, the approach proposed here provides additional information regarding how these demographic characteristics perform with respect to each of the two components and regarding the variability associated with each component. If the two components actually represent inspiration and expiration, the proposed model could be used to better examine how the demographic variables are associates with each phase. On average, the RR-interval length increases (i.e., heart rate decreases) about   12 and 11 milliseconds, for the first (inspiration) and second (expiration) components, respectively, with each year increase in age. Also, females had greater interval lengths than males by an average of 27 and 16 milliseconds, for the first and second components, respectively.
• The model produced two distinct variance estimates for the two components, namely 2 1 σ = 22.6 (95% CI = 22.19, 23.08) and 2 2 σ = 24.13 (95% CI = 23.66, 24.61). Examining the overlap in the CIs is indicative of the statistically significant differences in the variances of the two components (i.e., between the inspiration and expiration phases).
• The task by period interaction was also significant. However, the magnitude of the interaction for the first component (inspiration) was consistently higher than the second component, assuming a mean age of 19 and regardless of gender (Table 3, Figure 2). In general, the RR-interval lengths decreased from baseline to preparation to delivery. However, this decrease was not consistent across tasks. The IS and WS tasks exhibited quadratic trends. The IS RR-interval lengths increased from baseline to preparation but decreased at delivery, while the WS lengths decreased from baseline to preparation but increase at delivery. The lengths of all other tasks tended to decrease gradually from baseline to preparation to delivery.

Discussion
The use of the mixed-effects model is a valuable tool for modeling HRV type data. However, the current framework does not adequately address the issues of non-normality observed within-subjects. The HRV data measured by the RR-interval lengths tend to follow a twocomponent normal-mixture distribution. To accommodate this, a two-component MXME with normal random-effects seems to fit the data well. Under this framework, use of the EM solution to estimate the model parameters was proposed in this article. Although the EM solution has been applied for the estimation of parameters in normalmixture distributions, its application to the mixed-effects situation is novel. Furthermore, in contrast with the traditional methods, where the normal-mixture is assumed for the between-subject random-effects, the proposed method is specific to the within-subject non-normality. The within-subjects normal-mixture lends itself to the application of the EM algorithm and provides a much simpler posterior mean. By defining the incomplete data in terms of the weights of each observation belonging to one of the two components, the estimation could be carried out using existing mixed-effects methodology and software, as was demonstrated here with PROC IML and PROC MIXED in SAS.
The within-subject mixture models approach proposed here is more applicable than the widely studied between-subject mixture models. The model also could provide alternative clinical interpretations than the time and frequency domain approaches. In the application presented here, the two-component normal-mixture model was interpreted to be representative of the inspiration and expiration phases     of the breathing cycle. However, this is simply a conjecture that needs to be further examined and other possible interpretations should be explored, perhaps in consultation with the scientists in the HRV field. In situations where data on when inspiration and expiration phases occur is available along with the RR-interval data, one could compare the fit of the model proposed here with known phases to determine how well the model fits the RSA theory. In other types of HRV data, the components of the mixture could have other interpretations in HRV as well as outside of HRV. For example, if data are observed over a 24 hour period, the two dominating components might represent the circadian rhythm. Finally, there may be a need for fitting mixtures with more than two components and the proposed approach would easily extend for those situations.
To study the characteristics of the proposed method, a Monte-Carlo simulation study was also performed [21]. In this study, the mixture proportion, between-subject variability, number of subjects, number of observations within each subject, and the percentage of overlap in the two component distributions were allowed to vary. For each combination of these factors, 100 simulations were performed and the proposed MXME was applied. The bias and the mean squared error (MSE) of the estimates of the model parameters were computed. In order to study the characteristics of the method a (single replicate) full-factorial analysis was performed using absolute bias and MSE as the outcome variables. The results suggest that the method reproduces the mixture proportions the best and reproduces the between subject variance the least accurately, when the sample size is small. However, when the number of subjects increased from 5 subjects to 20, this estimation improves dramatically. The EM algorithm always converged to the global maxima and was found to be insensitive to the initial values.
The main advantage of the modeling is that it utilizes all of the data, thereby capturing information about within subject correlations, and allows for the researchers to include covariates in the model that may be associated with HRV (e.g. gender, age, etc.). Because this approach is shifting the paradigm from one that relies on the summary statistics such as the time and frequency domain measures (which have been well understood and utilized by the clinical researchers in the field) to modeling, the burden of interpreting the models and translating the results to the clinician rests on the statistician. Continued dialogue between the statistician and the clinician regarding the modeling approach could achieve this goal.