Journal of Electrocardiology
Vol. 25 Supplement
Autoregressive Spectral Models of Heart Rate Variability Practical Issues
Robert L. Burr, MSEE, PhD, and Marie J. Cowan, PhD, RN, FAAN
Abstract: Autoregressive time series model-based spectral estimates of heart period sequences can provide a parsimonious and visually attractive representation of the dynamics of interbeat intervals. While a corollary to Weld’s decomposition theorem implies that the discrete Fourier periodogram spectral estimate and the autoregressive spectral estimate converge asymptotically, there are practical differences between the two approaches when applied to short blocks of data. Autoregressive spectra can achieve good frequency resolution and excellent statistical stability on short segments of heart period data of sinus origin. However, the order of the autoregressive model (number of free parameters to be estimated) must be explicitly chosen, a decision that influences the trade-off of frequency resolution with statistical stability. Akaike’s Information Criterion (AIC), an information-theoretic rule for picking the optimum order, is sensitive to the aggregate amount of data in the analysis. Thus, the best model order for estimating the spectrum of a 4-minute segment of data will generally be lower than the best order for estimating an hourly spectrum based on averaging I5 4minute spectra. A major advantage of the autoregressive model approach to spectral analysis is the ease with which it can be extended to handle messy data frequently seen in heart rate variability studies. A number of autoregressivebased robust-resistant techniques are available for the analysis of heart period sequences that contain a high volume of nonsinus and other unusual beats intervals. A theoretically satisfying framework is also available for spectral analysis of unevenly sampled data and missing data. Key words: autoregressive, heart rate variability, heart period.
In the last decade there have been several hundred papers applying spectral analysis procedures in heart rate variability studies. At least 30 recent studies have used autoregressive (AR) model-based spectral
methods instead of the popular Fourier-based methods. Both approaches have proven their utility in clinical and research applications in cardiology, but many attitudes still prevail about the relative merits of Fourier-based methods and AR spectral analysis techniques. This paper discusses practical application of AR techniques in heart rate variability studies, and characterizes useful extensions of the AR model to robust-resistant spectral analysis, missing data, and unevenly sampled time series problems.
From the University of Washington, Seattle, Washington. Reprint requests: Robert L. Burr, MSEE, PhD, University of Washington, School of Nursing, Mail Stop SM-27, Seattle, WA 98195.
224
Autoregressive
Spectral Models of Heart Rate Variability
l
Burr and Cowan
225
The Heart Rate Variability Spectrum A spectrum is a kind of histogram that measures the relative strength of rhythm patterns over a graded continuum of periods or lengths. As depicted in Figure 1, fast rhythms map into peaks on the right of the spectrum, while slowly changing rhythms correspond to peaks on the left. The height of a peak in a spectrum is determined by the magnitude of a simple rhythmic process or by the additive mixture of a number of smaller rhythmic processes with similar time periods. A narrow peak usually indicates a simple and extremely regular periodicity, while a broad peak can occur either as a mixture of roughly similar rhythms in a band or as a single complex rhythm that is waxing and waning in periodicity over time. Engineers describe the spectrum as an analysis, or decomposition, of the energy of a signal. Indeed, the frequently used name “power spectrum” is derived from the physical analog. For those of us with more experience in the statistical rather than physical sciences, a better metaphor might be that the spectrum is an analysis of the variance of a time series, an isolation of variation into constitutent rhythms. Spectral analyses of R-R interval sequences allow the fractionation of the effects of sympathetic and parasympathetic influences on heart rhythm. lP3 These analyses are particularly effective at identifying and quantifying respiratory sinus arrhythmia. Midfrequency peaks in the spectrum usually reflect parasympathetically mediated respiratory sinus arrhythmia, and the low-frequency peaks indicate sympa-
Fig. 1.
A spectrum is a figure that encodes information about the amplitude and characteristic frequency of component rhythms in a signal. The location of a peak on the x axis indicates its frequency in Hertz (Hz = c/s), the inverse of the rhythm’s fundamental period in seconds. The amplitude of the rhythm determines the height of the spectral peak, usually plotted either in natural units (m?/Hz) or in logarithmic decibel units (dB = 10 * loglo {spectra}).
Fig. 2.
A Fourier frequency domain model presumes the signal is a linear combination of a very large population of
tuned oscillators. Equation: SC,,)= xl%, a,f,+,,) thetic and reciprocal sympathetic-parasympathetic reflex activity. The systemic reflexes mediated by the sympathetic nervous system appear strongly responsive only at frequencies below approximately 0.1 cycles per second (rhythm patterns 10 seconds or longer), while the parasympathetic nervous system can mediate much faster rhythmic processes. Therefore variance in the spectrum above 0.2 cycles per second (rhythms patterns 5 seconds in length or shorter) probably has been almost entirely mediated by vagal processes. The Fourier model (Fig. 2) assumes that the R-R interval sequence we observe is the weighted summation of a very large number of independent sinusoidal oscillators having different periods. The spectrum then is simply a graph of the weight or strength of the amplitude of these theoretical latent oscillators. The Fourier model is a useful mathematical artifice in many applications like this and a metaphor for understanding the time scale of parts of the signal. The analysis problem becomes one of trying to statistically estimate the amplitude of all these latent variables from the observed R-R data sequence. While this model may be conceptually convenient for some purposes, it is not a compelling physical model for the generation of the heart period sequences we observe. That is, none of us believe that there really are millions of little pure oscillators vibrating inside us that summate into the heart period signals we observe. From a philosophically scientific perspective, this is not a parsimonious model. There are a very large number of parameters to estimate with respect to the length of the data sequence. In fact, for a Fourier spectrum of a sequence n points long, there are n/2 parameters that must be estimated. This is equivalent to performing statistical regression using 50 independent variables, but fitting the regression model with only 100 cases. There are only two effective statistical degrees of freedom corresponding to each point in the raw power spectrum. This affects the statistical stability of the spectral estimate. Statistical stability
226
Journal of Electrocardiology
W, -/
Vol. 25 Supplement
Filter System h;-2wed time series)
Fig. 3. Most of the classical time domain models of a sta-
tionary time series assume that the observed signal is the
outcome of an unknown filter system operating on an unknown random uncorrelated (“white noise”) input. of the Fourier spectrum can be increased by smoothing the spectrum with a moving average filter in the frequency domain, but only at the expense of frequency resolution. Time domain models come at the problem from what initially appears to be a completely different direction. The heart period sequence we observe is assumed to be the outcome of an unknown filter system acting on an unknown white noise sequence (Fig. 3). Autoregressive spectra are derived from simple linear rules for predicting the current sample in the RR interval sequence from a set of samples in the recent past of the sequence itself (Fig. 4). The term “autoregressive” as applied to this kind of model suggests the action of prediction based on information in the local history of the signal itself. The parameters of the regression equation are usually solved using simple least-squares fitting algorithms. (As a practical matter, you can fit AR models quite adequately using the simple regression programs in your favorite statistics package as long as you take the formal hypothesis testing apparatus with a grain of salt.) The emphasis on linear processing of the local history makes the AR model somewhat compatible with physical differential equation models, which physiologists and biophysicists might use to describe credible generation mechanisms. The AR model can be contrasted with the moving average time series model, which predicts the current sample as filtered white noise (Fig. 5). The moving
s,., Fig, 4.
S”_5 S”_4 Sn-3 S”-2 S"-1
S”
l
’
l
l
An autoregressive (AR) time domain model predicts
the current sample from a finite number of samples (in this example, 4) of the recent past history of the signal itself. For each explicit determination of the regression weights of the AR model there corresponds a unique spectrum. Equation:
S, = x?= 1 bisn-i +
W”
s,,
S”_lj S”_4 Sn-3
s,., S”.i s,
l
l
l
l
Fig. 5. A pure moving
average time domain model describes the current sample of the sequence as a linear combination of an unknown white noise sequence. In this example, three random numbers are combined with fixed regression weights into an estimate of the current sample. For each choice of the regression weights of the moving average model there corresponds to unique spectrum. Equation:
S, = xf==, aiw,_i
average model does not explicitly use any information about the local history of the signal itself to make its prediction. The AR and moving average parametric models of the R-R interval sequence are time series models that fit in the time domain. However for each choice of the parameters of these time series models, there is one unique spectrum to which it corresponds. So fitting the parameters of a time series model to the data sequence can be thought of as equivalent to choosing a spectrum. Because we go through the intermediate step of computing a time domain model in order to obtain the spectrum, some authors prefer to talk about “autoregressive model-based spectra” rather than “autoregressive spectra.”
Fourier versus Autoregressive Spectra Practitioners in many fields argue about the relative merits of the Fourier and AR approaches to spectral analysis. However random process theorists tend to think these discussions are rather puerile because of a theoretical equivalence that can be traced between AR and Fourier derived spectra. The argument is as follows: 1. There is a theorem that says that an Fourier periodogram of a sequence n points long can be exactly represented as a moving average model of order n/2.4 2. As a consequence
of Wold’s decomposition theorem, the spectra derived from the moving average and AR models arbitrarily converge closely as their model orders increase.‘f6 3. Therefore Fourier spectra and AR spectra are equivalent.
Autoregressive
Spectral Models of Heart Rate Variability
This is a theoretical result that is doubly asymptotic. For convergence, it presumes that the data sequence is as long as we need it to be and that the model orders can rise without limit. There are, in fact, practical differences between Fourier periodogram estimate and an AR estimate derived from a short sequence of data. But we should expect these differences to be due to pragmatics and accidentals, like sequence boundary effects.
Practical Fourier Spectral Analysis Some generalizations that are frequently about practical Fourier application include:
l
Burr and Cowan
227
number of blocks that are aggregated into the average. 8. Sidelobe artifact can be reduced by special time domain windowing techniques that generally taper the endpoints of the finite segment. However, this reduces the frequency resolution and further increases the width of the statistical confidence intervals. 9. Although our eye is naturally drawn to the peaks of the spectrum, and much of our theorizing is in terms of the peaks as the loci of latent rhythms, the Fourier spectrum is better defined in the valleys than in the peaks. We know more about where a rhythm isn’t than where it is.
made
1. The technique is computationally efficient, whether implemented in hardware or software. 2. It is based on a well-understood least-squares lit of sinusoids. 3. It has a history of performing poorly on short lengths of data. 4. One problem is sidelobe bias, a kind of artifact that smears information from one part of the spectrum into other parts. But sidelobe distortions, if not removable, are at least well understood by practitioners. 5. There is a well-developed statistical theory for the Fourier spectral estimate. Even though it is an asymptotic theory, statisticians seem to be very comfortable with its application to finite length segments. One can write down credible 95% confidence intervals about each point in the spectrum. That is good news. 6. The bad news is that this beautiful asymptotic theory tells us that the raw periodogram has truly horrible statistical stability. Each point in the spectrum has an approximate variability of an average of two squared numbers. Of course, the average of two numbers is not a very stable estimate of the true mean. 7. There are two ways we can achieve better statistical stability. One is by blurring the spectrum with a frequency domain smoothing window, averaging together adjacent frequency bins (eg, with Daniel1 moving average smoothers). The blurring reduces the variability in the spectrum in a direct trade-off with frequency resolution. However, some information gets unavoidably smeared into other bands. Another way to reduce the variability is to average many blocks of spectra together. In Holter-based heart rate variability studies, we generally have many blocks that can be averaged together. The variability is roughly reduced as a factor of the square root of the
Practical Autoregressive Spectral Analysis If the underlying spectrum is fundamentally simple (perhaps composed of a couple of peaks and a scree slope), then an AR spectral estimator generally does a good job. It can be argued that heart rate variability spectra fall into this category. Sidelobes aren’t an issue (although analyses of short records of pure deterministic sine waves are subject to a related boundary effect called spectral line splitting. The Marple AR algorithm’ was developed to deal with this problem.) Autoregression is a conceptually simple algorithm, understandable by most people with a background in conventional statistical regression. But it is computationally more expensive than a Fourier spectrum (actually, this depends on model order: if the model order is low, the computational complexity is not necessarily much greater). In contrast to the Fourier spectrum the AR spectrum is better known in the peaks than the valleys (except when one is analyzing pure deterministic sine waves). The Model Order Selection Problem for Parametric Spectral Estimation In order to achieve all these advantages, one has to make a firm decision about the order of the parametric model. For the AR model this means specifying how many samples of the recent past to allow into the prediction equation. For the moving average model this means specifying the span of the white noise filter. Choosing too low an order will cause the model to ignore real structural detail in the signal, and the resulting spectrum will be too simple. Choosing too high an order will allow the time series model to overlit a finite segment, and will result in a spectrum having statistically unrepresentative detail.
228
Journal of Electrocardiology
Vol. 25 Supplement
The model order selection problem is a key decision in AR spectral analysis and it naturally results in some controversy among scientists. The order can be justified on the basis of theoretical knowledge of the complexity of the signal, with more complex signals needing a higher order. The order can be specified on the basis of esthetics and the analyst’s personal trade-offs between frequency resolution and statistical stability. The “best” order from a statistical or information theoretical sense can also estimated empirically using several criteria, one of the most popular and well known is the Akaike Information Criterion (AIC).8-‘o The AIC for a generalized regression problem can be described as a function of model order k: (I)
AIC(k) = log(Var(resid(k)))
+ 2-g
There are two terms here, the one on the left always decreasing, corresponding to better and better prediction with increasing order, and the one on the right always rising, corresponding to the statistical penalty for having more parameters to estimate. This tends to produce a V-shaped curve of the AIC with order (especially when depicted on a log-log plot), with the best order corresponding to the minimum. Note that the second term is basically the ratio of the number of parameters to data points. It is very important to realize that the statistically best order is a strong function of the total amount of data in the analysis. The equation is valid both for analyses of individual blocks of data, and also for multiple blocks analyzed separately and then averaged into a single spectrum. As depicted in Figure 6 for a single sudden cardiac arrest subject, the best order for estimating the spectrum of a 4-minute segment of heart period data is going to be much less than the best order for 15 4-minute spectra averaged into an hourly spectrum, and also less than the best order for 360 4minute spectra averaged into a 24-hour grand average spectrum. For a more representative view, Figure 7 shows for a sample of sudden cardiac arrest (SCA) subjects and the scatter of the statistically best model order as a function of the total amount of data in the analyses, from a single 120-second block of heart period data to 30 120-second blocks, spanning a full hour (3,600 seconds). Thus, consideration of the analytic structure of the AIC equation and empirical experiments both challenge the use of the AIC as a tool to uncover the intrinsic “true” model order of the generating process. The statistically best model order is a function of the aggregate length of the analysis sequence. If one is computing blocks of AR spectra in order to average them together, it is better to use somewhat higher model orders than would be appropriate for
Fig. 6.
The effect of block aggregation on the best AR model order. An empirical assessment of the best AR model order for a particular segment of data can be achieved with plots of the Akaike Information Criterion (AIC) as a function of model order, plotted here for three amounts of aggregate data. When plotted on a log-log scale, the characteristic “V” shape has its minimum (the MAIC) at the model order that is the “best” compromise of residual variance and number of parameters. The solid line describes the AIC for a 4-minute block of heart rate variability data from a particular subject, the dotted line corresponding to the aggregation of 15 4-minute blocks, and the dashed line characterizes the model order trade-offs for 360 4-minute blocks. The notch of the AIC or MAIC increases in model order with the total amount of data in the analysis, primarily because of the reduction in the penalty for increasing model parameters.
“8
Fig. 7.
There is great variability in the empirical “best” order determined by the minimum Akaike Information Criterion rule within subjects, between subjects, and across clinical conditions. While demonstrating considerable variation in a sample of SCA survivors, the typical best order shows a strong dependence on the aggregate amount of data in the AR analysis of heart rate variability.
Autoregressive
Spectral
Models
computing an isolated spectrum from a single block of data. The averaging over blocks will smooth out the moderate statistical jitter that comes with higher model orders. An hourly average might aggregate 15 4-minute spectra, which approximately reduces the single block variability by a factor of four. A 24-hour
of Heart Rate Variability
l
Burr and Cowan
229
average spectrum might combine 360 4-minute spectra, reducing the nominal variability by a factor of nearly 20. Illustrating the influence of analysis assumptions and technique on resulting spectral estimates, Figure 8 shows spectra computed from the same 4-minute
r
Fig. 8. Spectra derived from a single 4-minute block of heart rate variability data using six different techniques. In the first column are AR spectra using three specific choices of the model order. The second column contrasts Fourier-based spectra with different amounts of Daniel1 frequency domain smoothing. In both cases, the fundamental trade-off between frequency resolution and statistical stability can be observed. The Akaike Information Criterion best choice is the middle panel of the first column. (A) AR (5) spectrum using yule-walker; (B) AR (12) spectrum using yule-walker; (C) AR (55) spectrum using yule-walker; (D) Smoothed periodogram: bandwidth = 0.041798,95% CI is ( - 1.37089,1.63 179) dB; (E) Smoothed periodogram: bandwidth = 0.0143836, 95% CI is (-2.21383, 2.98338) dB; (F) Raw periodogram: bandwidth = 0.00120281, 95% CI is (- 5.87589, 17.5667) dB.
230
Journal of Electrocardiology
Vol. 25 Supplement
block of heart period data. On the top are AR spectra using three model orders: 5th order on the left (too smooth), 50th on the right (a little too noisy), and the AIC best choice, 12th order in the middle. Across the bottom are Fourier spectra using different amounts of Daniel1 smoothing. On the right is the raw periodogram (the vertical bar is the confidence interval), on the left is a heavily smoothed spectrum where the confidence interval has decreased but the frequency solution has deteriorated in proportion. In the middle is an esthetic trade-off between statistical stability and bandwidth for a Fourier-based estimate having about the same information content as the AR ( 12) spectrum directly above it.
Flexible Extensions of the Autoregressive Spectral Model
A major advantage of the AR model approach to spectral analysis is the ease with which it can be extended to handle messy data frequently seen in heart rate variability studies. A number of AR-based robust-resistant techniques are available for the analysis of heart period sequences that contain a high fraction of nonsinus and other unusual beats. A theoretically satisfying framework is also available for spectral analysis of unevenly sampled data and missing data.
Robust-resistant Autoregressive Spectral Estimation Techniques
The preprocessing of Holter tapes, before computation of statistics and spectra is a significant part of the operational expenses of heart rate variability-focused research projects. The current generation of intelligent Holter analyzers are excellent at identifying nonsinus beats, but careful heart rate variability analysts are still spending many hours overreading the beat coding on each tape. This extra care during the preanalysis step is expended because most heart rate variability measures, which are usually based on standard deviations and Fourier spectral estimates, are very sensitive to outliers and miscoded nonsinus beats. A single bad beat in an analysis block can contribute a relatively large amount of incoherent variance to the resulting statistical or spectral summary. There has been considerable recent research on robust-resistant techniques for regression, some of which is directly applicable to autoregession and spectral analysis. We have been experimenting with
three robust techniques for spectral analysis of heart rate variability data sets. 1. The first is an approach that provides a maximum likelihood estimate of the parameters under the assumption that the observed sequence is drawn from a mixture distribution of “good” and “bad” data. l l-l4 If the fraction of the outlier samples is sufficiently small (approximately lo%), an iterative solution, based on weighting the contribution of each data point to the analysis, can be shown to converge. This approach is fast and well regarded by statisticians, but the nonlinear iterations can get lost if the starting values are improperly chosen. 2. Martin’s filter/cleaner, and its generalizations, are also iterative techniques for fitting the essence of the data set by learning what data points are strange and what points are central. 1s-17 Odd points are adjusted to have less of an influence on model fit, the model is reestimated, and the revised residuals are examined for odd points again. These techniques employ data filtering and smoothing heuristics, and are thus quite accessible to researchers with an engineering background. 3. Least median of squared residuals is the most robust technique of the three described. 18,19Regular least squares autoregression can be considered a process of finding a fit that produces the least mean square of residuals from that fit. If the median is allowed to replace the concept of the mean, we arrive at a fitting strategy that doesn’t have an analytic solution like least squares, but which can tolerate a very high fraction of outliers. The approach essentially fits the middle “half” of the data. This technique becomes computationally expensive for model orders over 5, and essentially impossible for orders over 10. This is slightly too low to have direct utility in the heart rate variability spectral analysis problem. What has served well in our experiments is to fit a 5th or 7th order super-robust least median of squared residuals solution, then pass its data point weights and parameter estimates on to the generalized maximum likelihood estimate program to use as robust starting values. This approach combines the robust outlier detection capacity of the least median of squared residuals procedure with the well-regarded convergence properties of the generalized maximum likelihood estimate program. While experimental evaluation of these techniques in ambulatory monitoring data sets is incomplete, preliminary results are quite promising. Figure 9 shows an example where data that had 478 ectopic beats scattered thoughout 1 hour was analyzed. The top plot is the spectrum if nonsinus beat codes are ignored, simulating analysis of raw heart period data
Autoregressive
Spectral Models of Heart Rate Variability
Fig. 9. Comparison of AR (15) spectra derived from an hour of heart rate variability data using ( 1) no interpolation of bad beats, (2) conventional interpolation through bad beats, and (3) robust autoregression fit to all beats, including nonsinus rhythm SCA subject SO: hour 05 average AR ( 15) spectrum. that has not been coded. The spectrum is high and flat due to rhe amount of incoherent impulse noise from the bad beats. The middle spectrum was created using our conventional techniques, where we interpolate through bad beats with cubic spline interpolation. It has much less power at high frequencies. The lower plot is a robust AR fit to the same data seen by the top curve, that is, processed without knowledge of beat codes. Obviously, it more closely resembles the cleaned-up spectrum. This is fairly typical in our experiments on tapes collected from sudden cardiac arrest and myocardial infarction patients, with the robust spectrum having even slightly less power at high frequencies than our conventional estimate, and much less than the spectrum computed without preprocessing. We then applied the technique to a data set of archive tapes of normal subjects, who have very little ectopy. We postulated that the conventional and robust spectra would be quite comparable. This was true most of the time, but at a rate that was too high to ignore we ran into variations of Figure 10, where there were big differences between the usual AR and the robust AR spectrum of heart period data for normal subjects. Typically the main peaks are well defined, but something was very different between them. This effect alters estimates of variance within the usual bands, in particular, reducing the apparent power of the high-frequency component as estimated by a frequency band integral. The effect was much more pronounced in healthy normals than in any of the other clinical samples in our historical data set. Segments of raw time series data that corresponded to these discrepant spectra were examined, and the weights that the robust algorithms use to
l
Burr and Cowan
231
Fig. 10. Variations between spectra derived between robust and conventional AR algorithms can occur even when all the beats in the analysis frame are in sinus rhythm. This may place limits on the comparability of robust and nonrobust spectral estimates of heart rate variability. Normal subject F36: 1,000 seconds from hour 13 AR (10) spectrum.
adjust the influence of each point were also scrutinized. In most of these sections, there were visually obvious “spike” features (what we call “icicles”), as is shown in Figure 11, where a short-term arousal temporarily had increased the rate of normal beats to approximately 100 beats/min. These arousal features seem to be characteristic of healthy normal subjects, and counts of them can even be used as a quantitative measure of heart rate variability that discriminates clinical condition. The robust AR fitting algorithm has given points in this section of data of only approximately 1/20th the influence as given to other sections. Thus, this impulse nonstationarity, which contributes a significant amount of broadband power to the conventional AR (or Fourier) spectrum, is es-
Fig. 11. Spike features in normal
sinus rhythm have a major impact on the conventional AR or Fourier heart rate variability spectrum, but are essentially ignored by robust AR methods because they occupy such a small portion of the analysis frame.
232 Journal of Electrocardiology
Vol. 25 Supplement
sentially turned off in the robust spectrum. This may be viewed as a positive attribute because the robust technique is effectively ignoring short-term nonstationary structures in the analysis window. But because these generalized techniques result in qualitatively distinct spectra and quantitatively different statistical summaries than these produced by conventional techniques, experiments on a broader range of clinical and nonclinical subjects will be necessary to determine if statistics extracted from robust spectra have more or less relevance in clinical assessment and mortality studies.
Unevenly Sampled Data Most heart rate variability time series analysts interpolate the unevenly sampled natural heart period sequence onto an evenly spaced grid. This is because commonly available time series programs prefer data samples evenly spaced in time. Although there are potential differences with various interpolation techniques, the end results in ambulatory monitoringbased heart rate variability data sets don’t seem to be greatly influenced by whether we interpolate with a 0th order hold, a linear point connector interpolation, or a cubic spline fit. Because small practical differences have been noted between heart rate variability spectra derived from using different interpolation techniques on high-resolution electrocardiographic records (R-R intervals estimated at < 1 ms) collected in our laboratory, we conjecture that the spectral noise floor in the Holter-based data sets attribute to the coarser estimate of the R-R intervals inherent in current Holter technology (probably only known to be within 10 ms) swamps the spectral contributions due to interpolation errors. We also often interpolate through a certain amount of isolated missing data and short regions where nonsinus beats were coded. But interpolation is a systematic way of faking data where we really don’t know what the right value should be, or where we need to have representative data in a particular format for our analysis programs. We know that interpolants tend to reduce local variability in the sections with missing data or nonsinus beats, but introduce another kind of variance through interpolation error. Our use of interpolation is what Jaynes2’ and Tukey12 sometimes call an “ad hockery,” an expedient approximate atheoretic processing step. Finding a good spectral analysis procedure that can operate directly on heart period sequences in their natural unevenly spaced time sampling would remove this conceptual inelegance.
Jones2 1-23 has presented a strategy (and in an appendix, a functional computer program) that develops maximum likelihood estimates of the parameters of an AR model when the data is unevenly spaced or has sections that are missing. Jones creates a statespace model of a continuous AR process. He then uses a Kalman filter to consistently update the statespace parameters as each real data point arrives in unevenly sampled continuous time. The Kalman filter apparatus essentially predicts a continuous value based on its current parameters, measures the error between what it predicted and what it observed at each point when data becomes available, and uses that error to update its parameters to improve its prediction. The Kalman filter falls through sections of missing data because it only generates model parameter correction feedback when valid data are present. Preliminary experiments in the practical application of Jones’ 2l-z3 algorithm to heart rate variability spectral analysis have shown that it is highly consistent with our conventional AR spectral estimates based on cubic spline interpolation, but has a much greater computational burden. In particular, the Kalman filter-based estimate retains the nonrobust least mean of squares character of the conventional AR estimator. The current thrust of our methods’ research program is to integrate the interpolation-free Kalman filter approach with robust-resistant AR spectral techniques.
Future Directions in Heart Rate Variability Spectral Analysis Although current spectral analysis techniques allow quantification of theoretically important slow and fast modes of systematic variability, it is clear that we have only begun to scratch the surface of the information about biophysical adaptive capacity that is inherent in the instantaneous heart period sequence. The many individual generating subsystems probably have substantial nonlinearities, and these probably combine with nontrivial nonlinear interactions. Goldberger24,25 and others have even speculated that the nonlinear dynamics of heart rate variability may be chaotic, with l/f spectral behavior and subharmonic bifurcations. Heart period sequences collected during ambulatory monitoring are also unquestionably nonstationary, at a broad range of time scales and in both mean level and spectral content. Circadian nonstationarity spanning 24 hours is quite obvious in almost every record, as is beat-by-beat
Autoregressive Spectral Models of Heart Rate Variability nonstationarity associated with psychophysiological modulation. Unfortunately, both the traditional blockwise Fourier and AR approaches are inadequate to characterize nonlinear and nonstationary processes completely. Graphic depiction in three dimensions of spectral estimates of short blocks in consecutive time gives us hints of information-rich time-varying qualities in the heart period sequence, modulations that undoubtably reflect the functioning of the central autonomic nervous system. Higher order spectra may elucidate the essential nonlinear properties, and true time-varying spectral analysis techniques better may describe the nonstationary phenomenology of heart period fluctuations. Better understanding of the biological and psychophysiological subsystems that contribute to heart rate variability would allow the testing of structured nonlinear time-varying spectral models anchored on the hardstanding of physiological theory. The frankly phenomenological status of heart rate variability measures could even benefit from empirical metamodels (such as hidden Markov models) of the variation in spectral content over time. A satisfying spectral theory of heart rate variability should also directly confront the point process event sequence character of the fundamental cardiac event. Although much work remains, extensions of the versatile AR model have some relevance to each of these complexities and will continue to play a role in our understanding of the heart period process for many years to come.
6. 7.
8. 9. 10. 11. 12. 13. 14.
15.
16.
17. 18.
19. 20.
References 1. Askelrod S, Gordon D, Ubel F et al: Power spectrum analysis of heart rate fluctuations: a quantitative probe of beat-to-beat cardiovascular control. Science 2 13: 220, 1981 2. Pomeranz B, Macauley R, Caudill M et al: Assessment of autonomic function in humans by heart rate spectral analysis. Am J Physiol 248:H15 1, 1985 3. Katona P, Jih F: Respiratory sinus arrhythmia: noninvasive measure of parasympathethic cardiac control. J Appl Physiol 54:961, 1975 4. Kay SM, Marple SL: Spectrum analysis: a modem perspective. Proc IEEE 69: 1380, 1981 5. Whitman EC: The spectral analysis of discrete time series in terms of linear regressive models. Naval Ordi-
21.
22.
23.
24. 25.
l
Burr and Cowan
233
nance Labs Report NOLTR-70-109 White Oak, MD, 1974 Wold HOA: A study of the analysis of stationary time series. Almquist and Wiksell, Uppsala, 1938 Marple SL: A new autoregressive spectral analysis algorithm. IEEE Trans Acoustics, Speech, and Signal Processing ASSP-28:441, 1980 Akaike H: Statistical predictor identification. Ann Inst Statist Math 22:203, 1970 Akaike H: Autoregressive model fitting for control. Ann Inst Statist Math 23:163, 1971 Akaike H: A new look at statistical model identification. IEEE Trans Autom Contr AC-19:716, 1974 Huber PJ: Robust estimation. Annals of Math Stat 35: 73, 1964 Mosteller F, Tukey JW: Data analysis and regression. Addison-Wesley, Reading, 1977 Huber PJ: Robust statistics. Wiley, New York, 1981 Hampel FR, Ronchetti EM, Rousseeuw PJ, Stahel WA: Robust statistics: an approach based on influence functions. Wiley, New York, 1986 Martin RD: Robust estimation of autoregressive models. p. 228. In Brillinger DR and Tiao GC (eds): Directions in time series. Institute of Math Statistics, Hayward, 1980 Martin RD: Robust methods for time series. p. 683. In Findley DF (ed): Applied time series II. Academic Press, New York, 1981 Martin RD, Thomson DJ: Robust resistant spectrum estimates. Proc IEEE 70: 1097, 1983 Rousseeuw PJ: Least median of squares regression. Journal of the American Statistical Association 79: 871, 1984 Rousseeuw PJ, Leroy AM: Robust regression and outlier detection. Wiley, New York, 1987 Jaynes ET: On the rationale of maximum-entropy methods. Proc IEEE 70:939, 1982 Jones RH: Maximum likelihood fitting of ARMA models to time series with missing observations. Technometrics 22:389, 1980 Jones RH: Fitting a continuous-time autoregression to unevenly sampled discrete data. In Findley DF (ed): Applied time series analysis II. Academic Press, New York, 1981 Jones RH: Fitting multivariate time series models to unequally spaced data. Presented at the Office of Naval Research Symposium on Time Series Analysis of Irregularly Spaced Data, Texas AGM University, College Station, February 1983 Goldberger AL, Rigney D, West B: Chaos and fractals in human physiology. Sci Am 262:43, 1990 Goldberger AL, Rigney D, Mietus J et al: Nonlinear dynamics in sudden death syndrome: heart rate oscillations and bifurcations. Experientia 44:983, 1988
Vol. 25 Supplement
Autoregressive Spectral Models of Heart Rate Variability Practical Issues
Robert L. Burr, MSEE, PhD, and Marie J. Cowan, PhD, RN, FAAN
Abstract: Autoregressive time series model-based spectral estimates of heart period sequences can provide a parsimonious and visually attractive representation of the dynamics of interbeat intervals. While a corollary to Weld’s decomposition theorem implies that the discrete Fourier periodogram spectral estimate and the autoregressive spectral estimate converge asymptotically, there are practical differences between the two approaches when applied to short blocks of data. Autoregressive spectra can achieve good frequency resolution and excellent statistical stability on short segments of heart period data of sinus origin. However, the order of the autoregressive model (number of free parameters to be estimated) must be explicitly chosen, a decision that influences the trade-off of frequency resolution with statistical stability. Akaike’s Information Criterion (AIC), an information-theoretic rule for picking the optimum order, is sensitive to the aggregate amount of data in the analysis. Thus, the best model order for estimating the spectrum of a 4-minute segment of data will generally be lower than the best order for estimating an hourly spectrum based on averaging I5 4minute spectra. A major advantage of the autoregressive model approach to spectral analysis is the ease with which it can be extended to handle messy data frequently seen in heart rate variability studies. A number of autoregressivebased robust-resistant techniques are available for the analysis of heart period sequences that contain a high volume of nonsinus and other unusual beats intervals. A theoretically satisfying framework is also available for spectral analysis of unevenly sampled data and missing data. Key words: autoregressive, heart rate variability, heart period.
In the last decade there have been several hundred papers applying spectral analysis procedures in heart rate variability studies. At least 30 recent studies have used autoregressive (AR) model-based spectral
methods instead of the popular Fourier-based methods. Both approaches have proven their utility in clinical and research applications in cardiology, but many attitudes still prevail about the relative merits of Fourier-based methods and AR spectral analysis techniques. This paper discusses practical application of AR techniques in heart rate variability studies, and characterizes useful extensions of the AR model to robust-resistant spectral analysis, missing data, and unevenly sampled time series problems.
From the University of Washington, Seattle, Washington. Reprint requests: Robert L. Burr, MSEE, PhD, University of Washington, School of Nursing, Mail Stop SM-27, Seattle, WA 98195.
224
Autoregressive
Spectral Models of Heart Rate Variability
l
Burr and Cowan
225
The Heart Rate Variability Spectrum A spectrum is a kind of histogram that measures the relative strength of rhythm patterns over a graded continuum of periods or lengths. As depicted in Figure 1, fast rhythms map into peaks on the right of the spectrum, while slowly changing rhythms correspond to peaks on the left. The height of a peak in a spectrum is determined by the magnitude of a simple rhythmic process or by the additive mixture of a number of smaller rhythmic processes with similar time periods. A narrow peak usually indicates a simple and extremely regular periodicity, while a broad peak can occur either as a mixture of roughly similar rhythms in a band or as a single complex rhythm that is waxing and waning in periodicity over time. Engineers describe the spectrum as an analysis, or decomposition, of the energy of a signal. Indeed, the frequently used name “power spectrum” is derived from the physical analog. For those of us with more experience in the statistical rather than physical sciences, a better metaphor might be that the spectrum is an analysis of the variance of a time series, an isolation of variation into constitutent rhythms. Spectral analyses of R-R interval sequences allow the fractionation of the effects of sympathetic and parasympathetic influences on heart rhythm. lP3 These analyses are particularly effective at identifying and quantifying respiratory sinus arrhythmia. Midfrequency peaks in the spectrum usually reflect parasympathetically mediated respiratory sinus arrhythmia, and the low-frequency peaks indicate sympa-
Fig. 1.
A spectrum is a figure that encodes information about the amplitude and characteristic frequency of component rhythms in a signal. The location of a peak on the x axis indicates its frequency in Hertz (Hz = c/s), the inverse of the rhythm’s fundamental period in seconds. The amplitude of the rhythm determines the height of the spectral peak, usually plotted either in natural units (m?/Hz) or in logarithmic decibel units (dB = 10 * loglo {spectra}).
Fig. 2.
A Fourier frequency domain model presumes the signal is a linear combination of a very large population of
tuned oscillators. Equation: SC,,)= xl%, a,f,+,,) thetic and reciprocal sympathetic-parasympathetic reflex activity. The systemic reflexes mediated by the sympathetic nervous system appear strongly responsive only at frequencies below approximately 0.1 cycles per second (rhythm patterns 10 seconds or longer), while the parasympathetic nervous system can mediate much faster rhythmic processes. Therefore variance in the spectrum above 0.2 cycles per second (rhythms patterns 5 seconds in length or shorter) probably has been almost entirely mediated by vagal processes. The Fourier model (Fig. 2) assumes that the R-R interval sequence we observe is the weighted summation of a very large number of independent sinusoidal oscillators having different periods. The spectrum then is simply a graph of the weight or strength of the amplitude of these theoretical latent oscillators. The Fourier model is a useful mathematical artifice in many applications like this and a metaphor for understanding the time scale of parts of the signal. The analysis problem becomes one of trying to statistically estimate the amplitude of all these latent variables from the observed R-R data sequence. While this model may be conceptually convenient for some purposes, it is not a compelling physical model for the generation of the heart period sequences we observe. That is, none of us believe that there really are millions of little pure oscillators vibrating inside us that summate into the heart period signals we observe. From a philosophically scientific perspective, this is not a parsimonious model. There are a very large number of parameters to estimate with respect to the length of the data sequence. In fact, for a Fourier spectrum of a sequence n points long, there are n/2 parameters that must be estimated. This is equivalent to performing statistical regression using 50 independent variables, but fitting the regression model with only 100 cases. There are only two effective statistical degrees of freedom corresponding to each point in the raw power spectrum. This affects the statistical stability of the spectral estimate. Statistical stability
226
Journal of Electrocardiology
W, -/
Vol. 25 Supplement
Filter System h;-2wed time series)
Fig. 3. Most of the classical time domain models of a sta-
tionary time series assume that the observed signal is the
outcome of an unknown filter system operating on an unknown random uncorrelated (“white noise”) input. of the Fourier spectrum can be increased by smoothing the spectrum with a moving average filter in the frequency domain, but only at the expense of frequency resolution. Time domain models come at the problem from what initially appears to be a completely different direction. The heart period sequence we observe is assumed to be the outcome of an unknown filter system acting on an unknown white noise sequence (Fig. 3). Autoregressive spectra are derived from simple linear rules for predicting the current sample in the RR interval sequence from a set of samples in the recent past of the sequence itself (Fig. 4). The term “autoregressive” as applied to this kind of model suggests the action of prediction based on information in the local history of the signal itself. The parameters of the regression equation are usually solved using simple least-squares fitting algorithms. (As a practical matter, you can fit AR models quite adequately using the simple regression programs in your favorite statistics package as long as you take the formal hypothesis testing apparatus with a grain of salt.) The emphasis on linear processing of the local history makes the AR model somewhat compatible with physical differential equation models, which physiologists and biophysicists might use to describe credible generation mechanisms. The AR model can be contrasted with the moving average time series model, which predicts the current sample as filtered white noise (Fig. 5). The moving
s,., Fig, 4.
S”_5 S”_4 Sn-3 S”-2 S"-1
S”
l
’
l
l
An autoregressive (AR) time domain model predicts
the current sample from a finite number of samples (in this example, 4) of the recent past history of the signal itself. For each explicit determination of the regression weights of the AR model there corresponds a unique spectrum. Equation:
S, = x?= 1 bisn-i +
W”
s,,
S”_lj S”_4 Sn-3
s,., S”.i s,
l
l
l
l
Fig. 5. A pure moving
average time domain model describes the current sample of the sequence as a linear combination of an unknown white noise sequence. In this example, three random numbers are combined with fixed regression weights into an estimate of the current sample. For each choice of the regression weights of the moving average model there corresponds to unique spectrum. Equation:
S, = xf==, aiw,_i
average model does not explicitly use any information about the local history of the signal itself to make its prediction. The AR and moving average parametric models of the R-R interval sequence are time series models that fit in the time domain. However for each choice of the parameters of these time series models, there is one unique spectrum to which it corresponds. So fitting the parameters of a time series model to the data sequence can be thought of as equivalent to choosing a spectrum. Because we go through the intermediate step of computing a time domain model in order to obtain the spectrum, some authors prefer to talk about “autoregressive model-based spectra” rather than “autoregressive spectra.”
Fourier versus Autoregressive Spectra Practitioners in many fields argue about the relative merits of the Fourier and AR approaches to spectral analysis. However random process theorists tend to think these discussions are rather puerile because of a theoretical equivalence that can be traced between AR and Fourier derived spectra. The argument is as follows: 1. There is a theorem that says that an Fourier periodogram of a sequence n points long can be exactly represented as a moving average model of order n/2.4 2. As a consequence
of Wold’s decomposition theorem, the spectra derived from the moving average and AR models arbitrarily converge closely as their model orders increase.‘f6 3. Therefore Fourier spectra and AR spectra are equivalent.
Autoregressive
Spectral Models of Heart Rate Variability
This is a theoretical result that is doubly asymptotic. For convergence, it presumes that the data sequence is as long as we need it to be and that the model orders can rise without limit. There are, in fact, practical differences between Fourier periodogram estimate and an AR estimate derived from a short sequence of data. But we should expect these differences to be due to pragmatics and accidentals, like sequence boundary effects.
Practical Fourier Spectral Analysis Some generalizations that are frequently about practical Fourier application include:
l
Burr and Cowan
227
number of blocks that are aggregated into the average. 8. Sidelobe artifact can be reduced by special time domain windowing techniques that generally taper the endpoints of the finite segment. However, this reduces the frequency resolution and further increases the width of the statistical confidence intervals. 9. Although our eye is naturally drawn to the peaks of the spectrum, and much of our theorizing is in terms of the peaks as the loci of latent rhythms, the Fourier spectrum is better defined in the valleys than in the peaks. We know more about where a rhythm isn’t than where it is.
made
1. The technique is computationally efficient, whether implemented in hardware or software. 2. It is based on a well-understood least-squares lit of sinusoids. 3. It has a history of performing poorly on short lengths of data. 4. One problem is sidelobe bias, a kind of artifact that smears information from one part of the spectrum into other parts. But sidelobe distortions, if not removable, are at least well understood by practitioners. 5. There is a well-developed statistical theory for the Fourier spectral estimate. Even though it is an asymptotic theory, statisticians seem to be very comfortable with its application to finite length segments. One can write down credible 95% confidence intervals about each point in the spectrum. That is good news. 6. The bad news is that this beautiful asymptotic theory tells us that the raw periodogram has truly horrible statistical stability. Each point in the spectrum has an approximate variability of an average of two squared numbers. Of course, the average of two numbers is not a very stable estimate of the true mean. 7. There are two ways we can achieve better statistical stability. One is by blurring the spectrum with a frequency domain smoothing window, averaging together adjacent frequency bins (eg, with Daniel1 moving average smoothers). The blurring reduces the variability in the spectrum in a direct trade-off with frequency resolution. However, some information gets unavoidably smeared into other bands. Another way to reduce the variability is to average many blocks of spectra together. In Holter-based heart rate variability studies, we generally have many blocks that can be averaged together. The variability is roughly reduced as a factor of the square root of the
Practical Autoregressive Spectral Analysis If the underlying spectrum is fundamentally simple (perhaps composed of a couple of peaks and a scree slope), then an AR spectral estimator generally does a good job. It can be argued that heart rate variability spectra fall into this category. Sidelobes aren’t an issue (although analyses of short records of pure deterministic sine waves are subject to a related boundary effect called spectral line splitting. The Marple AR algorithm’ was developed to deal with this problem.) Autoregression is a conceptually simple algorithm, understandable by most people with a background in conventional statistical regression. But it is computationally more expensive than a Fourier spectrum (actually, this depends on model order: if the model order is low, the computational complexity is not necessarily much greater). In contrast to the Fourier spectrum the AR spectrum is better known in the peaks than the valleys (except when one is analyzing pure deterministic sine waves). The Model Order Selection Problem for Parametric Spectral Estimation In order to achieve all these advantages, one has to make a firm decision about the order of the parametric model. For the AR model this means specifying how many samples of the recent past to allow into the prediction equation. For the moving average model this means specifying the span of the white noise filter. Choosing too low an order will cause the model to ignore real structural detail in the signal, and the resulting spectrum will be too simple. Choosing too high an order will allow the time series model to overlit a finite segment, and will result in a spectrum having statistically unrepresentative detail.
228
Journal of Electrocardiology
Vol. 25 Supplement
The model order selection problem is a key decision in AR spectral analysis and it naturally results in some controversy among scientists. The order can be justified on the basis of theoretical knowledge of the complexity of the signal, with more complex signals needing a higher order. The order can be specified on the basis of esthetics and the analyst’s personal trade-offs between frequency resolution and statistical stability. The “best” order from a statistical or information theoretical sense can also estimated empirically using several criteria, one of the most popular and well known is the Akaike Information Criterion (AIC).8-‘o The AIC for a generalized regression problem can be described as a function of model order k: (I)
AIC(k) = log(Var(resid(k)))
+ 2-g
There are two terms here, the one on the left always decreasing, corresponding to better and better prediction with increasing order, and the one on the right always rising, corresponding to the statistical penalty for having more parameters to estimate. This tends to produce a V-shaped curve of the AIC with order (especially when depicted on a log-log plot), with the best order corresponding to the minimum. Note that the second term is basically the ratio of the number of parameters to data points. It is very important to realize that the statistically best order is a strong function of the total amount of data in the analysis. The equation is valid both for analyses of individual blocks of data, and also for multiple blocks analyzed separately and then averaged into a single spectrum. As depicted in Figure 6 for a single sudden cardiac arrest subject, the best order for estimating the spectrum of a 4-minute segment of heart period data is going to be much less than the best order for 15 4-minute spectra averaged into an hourly spectrum, and also less than the best order for 360 4minute spectra averaged into a 24-hour grand average spectrum. For a more representative view, Figure 7 shows for a sample of sudden cardiac arrest (SCA) subjects and the scatter of the statistically best model order as a function of the total amount of data in the analyses, from a single 120-second block of heart period data to 30 120-second blocks, spanning a full hour (3,600 seconds). Thus, consideration of the analytic structure of the AIC equation and empirical experiments both challenge the use of the AIC as a tool to uncover the intrinsic “true” model order of the generating process. The statistically best model order is a function of the aggregate length of the analysis sequence. If one is computing blocks of AR spectra in order to average them together, it is better to use somewhat higher model orders than would be appropriate for
Fig. 6.
The effect of block aggregation on the best AR model order. An empirical assessment of the best AR model order for a particular segment of data can be achieved with plots of the Akaike Information Criterion (AIC) as a function of model order, plotted here for three amounts of aggregate data. When plotted on a log-log scale, the characteristic “V” shape has its minimum (the MAIC) at the model order that is the “best” compromise of residual variance and number of parameters. The solid line describes the AIC for a 4-minute block of heart rate variability data from a particular subject, the dotted line corresponding to the aggregation of 15 4-minute blocks, and the dashed line characterizes the model order trade-offs for 360 4-minute blocks. The notch of the AIC or MAIC increases in model order with the total amount of data in the analysis, primarily because of the reduction in the penalty for increasing model parameters.
“8
Fig. 7.
There is great variability in the empirical “best” order determined by the minimum Akaike Information Criterion rule within subjects, between subjects, and across clinical conditions. While demonstrating considerable variation in a sample of SCA survivors, the typical best order shows a strong dependence on the aggregate amount of data in the AR analysis of heart rate variability.
Autoregressive
Spectral
Models
computing an isolated spectrum from a single block of data. The averaging over blocks will smooth out the moderate statistical jitter that comes with higher model orders. An hourly average might aggregate 15 4-minute spectra, which approximately reduces the single block variability by a factor of four. A 24-hour
of Heart Rate Variability
l
Burr and Cowan
229
average spectrum might combine 360 4-minute spectra, reducing the nominal variability by a factor of nearly 20. Illustrating the influence of analysis assumptions and technique on resulting spectral estimates, Figure 8 shows spectra computed from the same 4-minute
r
Fig. 8. Spectra derived from a single 4-minute block of heart rate variability data using six different techniques. In the first column are AR spectra using three specific choices of the model order. The second column contrasts Fourier-based spectra with different amounts of Daniel1 frequency domain smoothing. In both cases, the fundamental trade-off between frequency resolution and statistical stability can be observed. The Akaike Information Criterion best choice is the middle panel of the first column. (A) AR (5) spectrum using yule-walker; (B) AR (12) spectrum using yule-walker; (C) AR (55) spectrum using yule-walker; (D) Smoothed periodogram: bandwidth = 0.041798,95% CI is ( - 1.37089,1.63 179) dB; (E) Smoothed periodogram: bandwidth = 0.0143836, 95% CI is (-2.21383, 2.98338) dB; (F) Raw periodogram: bandwidth = 0.00120281, 95% CI is (- 5.87589, 17.5667) dB.
230
Journal of Electrocardiology
Vol. 25 Supplement
block of heart period data. On the top are AR spectra using three model orders: 5th order on the left (too smooth), 50th on the right (a little too noisy), and the AIC best choice, 12th order in the middle. Across the bottom are Fourier spectra using different amounts of Daniel1 smoothing. On the right is the raw periodogram (the vertical bar is the confidence interval), on the left is a heavily smoothed spectrum where the confidence interval has decreased but the frequency solution has deteriorated in proportion. In the middle is an esthetic trade-off between statistical stability and bandwidth for a Fourier-based estimate having about the same information content as the AR ( 12) spectrum directly above it.
Flexible Extensions of the Autoregressive Spectral Model
A major advantage of the AR model approach to spectral analysis is the ease with which it can be extended to handle messy data frequently seen in heart rate variability studies. A number of AR-based robust-resistant techniques are available for the analysis of heart period sequences that contain a high fraction of nonsinus and other unusual beats. A theoretically satisfying framework is also available for spectral analysis of unevenly sampled data and missing data.
Robust-resistant Autoregressive Spectral Estimation Techniques
The preprocessing of Holter tapes, before computation of statistics and spectra is a significant part of the operational expenses of heart rate variability-focused research projects. The current generation of intelligent Holter analyzers are excellent at identifying nonsinus beats, but careful heart rate variability analysts are still spending many hours overreading the beat coding on each tape. This extra care during the preanalysis step is expended because most heart rate variability measures, which are usually based on standard deviations and Fourier spectral estimates, are very sensitive to outliers and miscoded nonsinus beats. A single bad beat in an analysis block can contribute a relatively large amount of incoherent variance to the resulting statistical or spectral summary. There has been considerable recent research on robust-resistant techniques for regression, some of which is directly applicable to autoregession and spectral analysis. We have been experimenting with
three robust techniques for spectral analysis of heart rate variability data sets. 1. The first is an approach that provides a maximum likelihood estimate of the parameters under the assumption that the observed sequence is drawn from a mixture distribution of “good” and “bad” data. l l-l4 If the fraction of the outlier samples is sufficiently small (approximately lo%), an iterative solution, based on weighting the contribution of each data point to the analysis, can be shown to converge. This approach is fast and well regarded by statisticians, but the nonlinear iterations can get lost if the starting values are improperly chosen. 2. Martin’s filter/cleaner, and its generalizations, are also iterative techniques for fitting the essence of the data set by learning what data points are strange and what points are central. 1s-17 Odd points are adjusted to have less of an influence on model fit, the model is reestimated, and the revised residuals are examined for odd points again. These techniques employ data filtering and smoothing heuristics, and are thus quite accessible to researchers with an engineering background. 3. Least median of squared residuals is the most robust technique of the three described. 18,19Regular least squares autoregression can be considered a process of finding a fit that produces the least mean square of residuals from that fit. If the median is allowed to replace the concept of the mean, we arrive at a fitting strategy that doesn’t have an analytic solution like least squares, but which can tolerate a very high fraction of outliers. The approach essentially fits the middle “half” of the data. This technique becomes computationally expensive for model orders over 5, and essentially impossible for orders over 10. This is slightly too low to have direct utility in the heart rate variability spectral analysis problem. What has served well in our experiments is to fit a 5th or 7th order super-robust least median of squared residuals solution, then pass its data point weights and parameter estimates on to the generalized maximum likelihood estimate program to use as robust starting values. This approach combines the robust outlier detection capacity of the least median of squared residuals procedure with the well-regarded convergence properties of the generalized maximum likelihood estimate program. While experimental evaluation of these techniques in ambulatory monitoring data sets is incomplete, preliminary results are quite promising. Figure 9 shows an example where data that had 478 ectopic beats scattered thoughout 1 hour was analyzed. The top plot is the spectrum if nonsinus beat codes are ignored, simulating analysis of raw heart period data
Autoregressive
Spectral Models of Heart Rate Variability
Fig. 9. Comparison of AR (15) spectra derived from an hour of heart rate variability data using ( 1) no interpolation of bad beats, (2) conventional interpolation through bad beats, and (3) robust autoregression fit to all beats, including nonsinus rhythm SCA subject SO: hour 05 average AR ( 15) spectrum. that has not been coded. The spectrum is high and flat due to rhe amount of incoherent impulse noise from the bad beats. The middle spectrum was created using our conventional techniques, where we interpolate through bad beats with cubic spline interpolation. It has much less power at high frequencies. The lower plot is a robust AR fit to the same data seen by the top curve, that is, processed without knowledge of beat codes. Obviously, it more closely resembles the cleaned-up spectrum. This is fairly typical in our experiments on tapes collected from sudden cardiac arrest and myocardial infarction patients, with the robust spectrum having even slightly less power at high frequencies than our conventional estimate, and much less than the spectrum computed without preprocessing. We then applied the technique to a data set of archive tapes of normal subjects, who have very little ectopy. We postulated that the conventional and robust spectra would be quite comparable. This was true most of the time, but at a rate that was too high to ignore we ran into variations of Figure 10, where there were big differences between the usual AR and the robust AR spectrum of heart period data for normal subjects. Typically the main peaks are well defined, but something was very different between them. This effect alters estimates of variance within the usual bands, in particular, reducing the apparent power of the high-frequency component as estimated by a frequency band integral. The effect was much more pronounced in healthy normals than in any of the other clinical samples in our historical data set. Segments of raw time series data that corresponded to these discrepant spectra were examined, and the weights that the robust algorithms use to
l
Burr and Cowan
231
Fig. 10. Variations between spectra derived between robust and conventional AR algorithms can occur even when all the beats in the analysis frame are in sinus rhythm. This may place limits on the comparability of robust and nonrobust spectral estimates of heart rate variability. Normal subject F36: 1,000 seconds from hour 13 AR (10) spectrum.
adjust the influence of each point were also scrutinized. In most of these sections, there were visually obvious “spike” features (what we call “icicles”), as is shown in Figure 11, where a short-term arousal temporarily had increased the rate of normal beats to approximately 100 beats/min. These arousal features seem to be characteristic of healthy normal subjects, and counts of them can even be used as a quantitative measure of heart rate variability that discriminates clinical condition. The robust AR fitting algorithm has given points in this section of data of only approximately 1/20th the influence as given to other sections. Thus, this impulse nonstationarity, which contributes a significant amount of broadband power to the conventional AR (or Fourier) spectrum, is es-
Fig. 11. Spike features in normal
sinus rhythm have a major impact on the conventional AR or Fourier heart rate variability spectrum, but are essentially ignored by robust AR methods because they occupy such a small portion of the analysis frame.
232 Journal of Electrocardiology
Vol. 25 Supplement
sentially turned off in the robust spectrum. This may be viewed as a positive attribute because the robust technique is effectively ignoring short-term nonstationary structures in the analysis window. But because these generalized techniques result in qualitatively distinct spectra and quantitatively different statistical summaries than these produced by conventional techniques, experiments on a broader range of clinical and nonclinical subjects will be necessary to determine if statistics extracted from robust spectra have more or less relevance in clinical assessment and mortality studies.
Unevenly Sampled Data Most heart rate variability time series analysts interpolate the unevenly sampled natural heart period sequence onto an evenly spaced grid. This is because commonly available time series programs prefer data samples evenly spaced in time. Although there are potential differences with various interpolation techniques, the end results in ambulatory monitoringbased heart rate variability data sets don’t seem to be greatly influenced by whether we interpolate with a 0th order hold, a linear point connector interpolation, or a cubic spline fit. Because small practical differences have been noted between heart rate variability spectra derived from using different interpolation techniques on high-resolution electrocardiographic records (R-R intervals estimated at < 1 ms) collected in our laboratory, we conjecture that the spectral noise floor in the Holter-based data sets attribute to the coarser estimate of the R-R intervals inherent in current Holter technology (probably only known to be within 10 ms) swamps the spectral contributions due to interpolation errors. We also often interpolate through a certain amount of isolated missing data and short regions where nonsinus beats were coded. But interpolation is a systematic way of faking data where we really don’t know what the right value should be, or where we need to have representative data in a particular format for our analysis programs. We know that interpolants tend to reduce local variability in the sections with missing data or nonsinus beats, but introduce another kind of variance through interpolation error. Our use of interpolation is what Jaynes2’ and Tukey12 sometimes call an “ad hockery,” an expedient approximate atheoretic processing step. Finding a good spectral analysis procedure that can operate directly on heart period sequences in their natural unevenly spaced time sampling would remove this conceptual inelegance.
Jones2 1-23 has presented a strategy (and in an appendix, a functional computer program) that develops maximum likelihood estimates of the parameters of an AR model when the data is unevenly spaced or has sections that are missing. Jones creates a statespace model of a continuous AR process. He then uses a Kalman filter to consistently update the statespace parameters as each real data point arrives in unevenly sampled continuous time. The Kalman filter apparatus essentially predicts a continuous value based on its current parameters, measures the error between what it predicted and what it observed at each point when data becomes available, and uses that error to update its parameters to improve its prediction. The Kalman filter falls through sections of missing data because it only generates model parameter correction feedback when valid data are present. Preliminary experiments in the practical application of Jones’ 2l-z3 algorithm to heart rate variability spectral analysis have shown that it is highly consistent with our conventional AR spectral estimates based on cubic spline interpolation, but has a much greater computational burden. In particular, the Kalman filter-based estimate retains the nonrobust least mean of squares character of the conventional AR estimator. The current thrust of our methods’ research program is to integrate the interpolation-free Kalman filter approach with robust-resistant AR spectral techniques.
Future Directions in Heart Rate Variability Spectral Analysis Although current spectral analysis techniques allow quantification of theoretically important slow and fast modes of systematic variability, it is clear that we have only begun to scratch the surface of the information about biophysical adaptive capacity that is inherent in the instantaneous heart period sequence. The many individual generating subsystems probably have substantial nonlinearities, and these probably combine with nontrivial nonlinear interactions. Goldberger24,25 and others have even speculated that the nonlinear dynamics of heart rate variability may be chaotic, with l/f spectral behavior and subharmonic bifurcations. Heart period sequences collected during ambulatory monitoring are also unquestionably nonstationary, at a broad range of time scales and in both mean level and spectral content. Circadian nonstationarity spanning 24 hours is quite obvious in almost every record, as is beat-by-beat
Autoregressive Spectral Models of Heart Rate Variability nonstationarity associated with psychophysiological modulation. Unfortunately, both the traditional blockwise Fourier and AR approaches are inadequate to characterize nonlinear and nonstationary processes completely. Graphic depiction in three dimensions of spectral estimates of short blocks in consecutive time gives us hints of information-rich time-varying qualities in the heart period sequence, modulations that undoubtably reflect the functioning of the central autonomic nervous system. Higher order spectra may elucidate the essential nonlinear properties, and true time-varying spectral analysis techniques better may describe the nonstationary phenomenology of heart period fluctuations. Better understanding of the biological and psychophysiological subsystems that contribute to heart rate variability would allow the testing of structured nonlinear time-varying spectral models anchored on the hardstanding of physiological theory. The frankly phenomenological status of heart rate variability measures could even benefit from empirical metamodels (such as hidden Markov models) of the variation in spectral content over time. A satisfying spectral theory of heart rate variability should also directly confront the point process event sequence character of the fundamental cardiac event. Although much work remains, extensions of the versatile AR model have some relevance to each of these complexities and will continue to play a role in our understanding of the heart period process for many years to come.
6. 7.
8. 9. 10. 11. 12. 13. 14.
15.
16.
17. 18.
19. 20.
References 1. Askelrod S, Gordon D, Ubel F et al: Power spectrum analysis of heart rate fluctuations: a quantitative probe of beat-to-beat cardiovascular control. Science 2 13: 220, 1981 2. Pomeranz B, Macauley R, Caudill M et al: Assessment of autonomic function in humans by heart rate spectral analysis. Am J Physiol 248:H15 1, 1985 3. Katona P, Jih F: Respiratory sinus arrhythmia: noninvasive measure of parasympathethic cardiac control. J Appl Physiol 54:961, 1975 4. Kay SM, Marple SL: Spectrum analysis: a modem perspective. Proc IEEE 69: 1380, 1981 5. Whitman EC: The spectral analysis of discrete time series in terms of linear regressive models. Naval Ordi-
21.
22.
23.
24. 25.
l
Burr and Cowan
233
nance Labs Report NOLTR-70-109 White Oak, MD, 1974 Wold HOA: A study of the analysis of stationary time series. Almquist and Wiksell, Uppsala, 1938 Marple SL: A new autoregressive spectral analysis algorithm. IEEE Trans Acoustics, Speech, and Signal Processing ASSP-28:441, 1980 Akaike H: Statistical predictor identification. Ann Inst Statist Math 22:203, 1970 Akaike H: Autoregressive model fitting for control. Ann Inst Statist Math 23:163, 1971 Akaike H: A new look at statistical model identification. IEEE Trans Autom Contr AC-19:716, 1974 Huber PJ: Robust estimation. Annals of Math Stat 35: 73, 1964 Mosteller F, Tukey JW: Data analysis and regression. Addison-Wesley, Reading, 1977 Huber PJ: Robust statistics. Wiley, New York, 1981 Hampel FR, Ronchetti EM, Rousseeuw PJ, Stahel WA: Robust statistics: an approach based on influence functions. Wiley, New York, 1986 Martin RD: Robust estimation of autoregressive models. p. 228. In Brillinger DR and Tiao GC (eds): Directions in time series. Institute of Math Statistics, Hayward, 1980 Martin RD: Robust methods for time series. p. 683. In Findley DF (ed): Applied time series II. Academic Press, New York, 1981 Martin RD, Thomson DJ: Robust resistant spectrum estimates. Proc IEEE 70: 1097, 1983 Rousseeuw PJ: Least median of squares regression. Journal of the American Statistical Association 79: 871, 1984 Rousseeuw PJ, Leroy AM: Robust regression and outlier detection. Wiley, New York, 1987 Jaynes ET: On the rationale of maximum-entropy methods. Proc IEEE 70:939, 1982 Jones RH: Maximum likelihood fitting of ARMA models to time series with missing observations. Technometrics 22:389, 1980 Jones RH: Fitting a continuous-time autoregression to unevenly sampled discrete data. In Findley DF (ed): Applied time series analysis II. Academic Press, New York, 1981 Jones RH: Fitting multivariate time series models to unequally spaced data. Presented at the Office of Naval Research Symposium on Time Series Analysis of Irregularly Spaced Data, Texas AGM University, College Station, February 1983 Goldberger AL, Rigney D, West B: Chaos and fractals in human physiology. Sci Am 262:43, 1990 Goldberger AL, Rigney D, Mietus J et al: Nonlinear dynamics in sudden death syndrome: heart rate oscillations and bifurcations. Experientia 44:983, 1988
