Algorithm A REGULARITYSTATISTIC FOR MEDICAL DATA ANALYSIS Steven M. Pincus, PhD, Igor M. Gladstone, MD,* and Richard A. Ehrenkranz, MD

Pincus SM, Gladstone IM, Ehrenkranz RA. A regularity statistic for medical data analysis. J Clin Monit 1991;7:335-345 ABSTRACT. A new statistic has been developed to quantify the

amount o f regularity in data. This statistic, ApEn (approximate entropy), appears to have potential application throughout medicine, notably in electrocardiogram and related heart rate data analyses and in the analysis of endocrine hormone release pulsatility. The focus of this article is ApEn. We commence with a simple example of what we are trying to discern. We then discuss exact regularity statistics and practical difficulties o f using them in data analysis. The mathematic formula development for ApEn condudes the Solution section. We next discuss the two key input requirements, followed by an account o f a pilot study successfully applying ApEn to neonatal heart rate analysis. We conclude with the important topic of ApEn as a relative (not absolute) measure, potential applications, and some caveats about appropriate usage of ApEn. Appendix A provides example ApEn and entropy computations to develop intuition about these measures. Appendix B contains a Fortran program for computing ApEn. This article can be read from at least three viewpoints. The practitioner who wishes to use a "black box" to measure regularity should concentrate on the exact formula, choices for the two input variables, potential applications, and caveats about appropriate usage. The physician who wishes to apply ApEn to heart rate analysis should particularly note the pilot study discussion. The more mathematically inclined reader will benefit from discussions of the relative (comparative) property of ApEn and from Appendix A. KEY WORDS. Heart: monitoring.

THE PROBLEM

From the Department of Pediatrics, Yale University School of Medicine, New Haven, CT. *Current address: Department of Pediatrics, National Naval Medical Center, Bethesda, MD 20814. I. M. G. is an active-duty Naval officer. The opinions expressed herein are those of the authors and are not to be construed as reflecting the views of the Navy Department, the Naval Service at large, or the Department of Defense. Address correspondence to Dr Pincus, 990 Moose Hill Road, Guilford, CT 06437.

In the past 20 years, there has been extensive study o f the ability o f electrocardiograms (ECGs), and the derived beat-to-beat differences (BBDs) in heart rate, to identify high-risk patients. T h e physician has been trained to recognize a variety o f abnormalities in the data. In B B D analysis, insight has been gained b y l o o k ing at variability. Recently, research involving E C G data has used such techniques as Fourier spectral analysis [1] and the identification and changing behavior o f strange attractors during ventricular fibrillation [2]. There is an increasing recognition that n e w perspectives and analytic tools to evaluate E C G and B B D data m a y p r o vide n e w clinical information. O n e such n e w perspective is the quantification o f the a m o u n t o f regularity in (heart rate) time-series data. W e recently evaluated heartbeat data in a pilot study in an attempt to discern differences between healthy and sick infants, based o n a suitable quantification o f regularity in the B B D data. W e hypothesized that the heartbeats o f the healthier infants w o u l d have greater randomness Copyright © 1991 by Little, Brown and Company

335

336 Journal of Clinical Monitoring Vol 7 No 4 October 1991

and, hence, appear less regular than those o f the sicker group. We also hypothesized that the amount o f randomness in BBD data would increase with patient recovery. These hypotheses are suggested by cardiology data showing loss o f spectral reserve in Fourier analyses o f BBDs o f ill hearts [3]. T o quantify regularity (and randomness), we used the Kolmogorov-Sinai (K-S) entropy as our measure. An initial, straightforward application o f a formula for K-S entropy, posed by Grassberger and Procaccia [4] and modified by Takens [5], yielded intuitively incorrect results. Closer inspection of the formula showed that the low-magnitude noise present in the data greatly affected the calculation results. It also became apparent that to attempt to achieve convergence o f the entropy measure, rapidly growing time demands would be placed on the computational resources. The problem was to determine a suitable formula to quantify the amount o f regularity in BBD data and thereby provide a new perspective to data analysis. The general problem was to develop a single, universal statistic to quantify the amount o f regularity in time-series data as a single number. This involved a different philosophy than do algorithms that search for particular pattern features in data. Representative o f these latter algorithms are the pulse detection algorithms central to endocrine hormone analysis [6,7] that identify the number o f peaks in pulsatile data and their locations. Before presenting a detailed discussion o f regularity, we contrast two time series o f data to illustrate what we are trying to measure. In both time series, the data represent the beat-to-beat heart rate, in beats per minute, at equally spaced time intervals. In series 1, the first few terms are 90, 70, 90, 70, 90, 70, 90, 70, 90, 70, 90, 70, 90, 70, 90, 70 . . . . This long-term series alternates rates o f 90 and 70 beats/ min. In series 2, the data are 90, 70, 70, 90, 90, 90, 70, 70, 90, 90, 70, 90, 70, 70, 90, 70 . . . . Each term o f this series has either a value of 90 or 70, randomly chosen, each with probability 1/2. M o m e n t statistics, such as mean and variance, will not distinguish between these two series. The mean of both these series is 80; the standard deviation is 10. N o r will rank order statistics, such as the median (any num-

ber between 70 and 90 in this case), distinguish between these series. Yet, series 1 is "perfectly regular"; knowing that one term has the value 90 allows one to predict with certainty that the next term will have the value 70. Series 2 is randomly valued; knowing that one term has the value 90 gives no insight into whether the next term will have value 90 or 70. The regularity statistics under discussion, entropy and ApEn, distinguish between these series, with entropy and ApEn having a value o f 0 for series 1 and a value o f In 2 (natural log) for series 2. The mathematic approach that will be described will not only allow one to distinguish between such obviously different series, but also to detect subtler differences in regularity o f data. It is this property that promises potential significance in identifying subtle heart arrhythmias, abnormal endocrine hormone secretion patterns, and myriad other abnormalities. THE SOLUTION We modified the exact K-S entropy formula to provide meaningful data analysis in reasonable time. The resultant statistic, which we call ApEn (approximate entropy), yielded excellent results, clearly discriminating between two groups in a pilot study. More generally, ApEn seems w o r t h y as a robust, generally computable statistic. Before discussing the ApEn algorithm, we present a brief historic perspective and note practical difficulties in the use o f exact regularity statistics. The historic development o f mathematics to quantify regularity has centered around various types o f entropy measures. Entropy, in a different context, has been an integral part o f the modern quantitative development o f thermodynamics, statistical mechanics, and information theory. Although, intuitively, the entropy quantifications in physics address the issues o f radomness and regularity, the formulas themselves involve integrals and derivatives o f k n o w n functions, such as work, temperature, and energy [8]. In modern probability theory, entropy is explicitly defined, given a probability distribution (measure) for elements o f a set [9]. This definition coincides with intuition in that systems having more random probability distributions have greater entropy. Nonetheless, these approaches to entropy definition are not directly applicable to time-series data analysis. K-S entropy (see ref. 10 for an excellent, albeit highly technical, description and discussion) generalizes the probabilist's definition o f entropy, in a theoretical setting, and paves the way to entropy formulas for timeseries data, as discussed below. There has been particularly keen interest in the development o f these formulas

Algorithm: Pincus et ah Regularity Statisticfor Data 337

in the last 10 years, since entropy has been shown to be a critical " s u m m a r y " statistic in nonlinear dynamical system analysis and chaos [11]. In 1983, Grassberger and Procaccia developed a formula, based on the K-S entropy, to measure the entropy o f a time series [4]; this formula, and a slight variation produced by Takens [5], have become the "standard" entropy and regularity measures for use with time-series data.

Regularity Statistics: The Need for an Approximation Statistic These formulas, however, were developed with strange attractors, central objects in chaos, in mind; such objects result from noiseless data taken from deterministic systems. The formulas are technical in appearance, represented as "triple limits." There are two paramount issues in applying these formulas to experimental data: computational demands and noise effects. COMPUTATIONALDEMANDS. W o l f e t al [12] conclude that an accurate estimation o f the characteristic exponents (one o f the three principal statistics used in chaos theory) for an attractor o f dimension d requires from 10a to 30 a data values. The same difficulty is seen with fractal dimension calculations [13]. O f crucial importance to the development here is the fact that the similarity between the algorithms to compute dimension [4] and K-S entropy implies that an accurate entropy calculation similarly requires vast amounts o f data. This data requirement yields computationally unusable algorithms for medium- and large-dimensional systems, since the associated algorithms grow exponentially in computation time with the dimensionatity o f the underlying system. There is no reason to anticipate that data typically encountered in medicine, especially from such complex systems as those that generate heart potential and brain wave pulses, are low dimensional, so this difficulty is crucial in our context. NOISE EFFECTS. Noise is an even greater problem in regularity statistics. The formulas o f Grassberger and Procaccia and Takens for the entropy are "wildly" discontinuous to system noise. Let the experimental time series u(1), u(2) . . . . correspond to scalar measurements taken uniformly in time. If v(t), v(2) . . . . is generated as a small quantity e times white noise, independent o f the sequence u(1), u(2) . . . . . then it is straightforward to verify that the entropy o f the sequence u(l) + v(1), u(2) + v(2) . . . . is infinity, independent o f the u(i)s and e. Hence, a blind application o f this regularity statistic will only evaluate system noise, not underlying system properties. See Appendix A, especially series 7, for illustrative examples o f this crucial problem.

This acute sensitivity to noise has been noted by the developers o f regularity statistics. There are three approaches to handling the noise. The first, adopted in the definition o f ApEn, is to do data comparisons on a scale larger than most o f the noise. The second approach, discussed by W o l f et al [12], is to perform a low-pass filtering o f the data before doing calculations. The third approach, discussed by Fraser [14], is based on separating noise from deterministic effects via algorithms to calculate marginal redundancy and mutual information.

Approximate Entropy Given the above, two choices present themselves. One can hope that the underlying system is very low dimensional, with negligible noise; then, the above statistics may be usable. Alternatively, one can modify an existing regularity statistic to derive a usable, robust measure for a general time series, at the expense o f some theoretic properties. We adopted the tatter philosophy, expecting the modification to be sufficient for its intended use. We next describe the algorithms for computing K-S entropy [4] and ApEn, noting the differences between them.

Both Algorithms STEP 1. Form a time series o f data u(1), u(2) . . . . u(N). These are N raw data values from measurements equally spaced in time. STEP 2. Fix m, an integer, and r, a positive real number. The value o f m represents the length o f compared runs o f data, and r specifies a filtering level (see discussion in Constraints section). STEP 3. Form a sequence o f vectors x(1), x(2) . . . . . x(N - m + 1) in R m, real m-dimensional space, defined by x(i) = [u(i) . . . . . u(i + m - 1)]. STEP 4. Use the sequence x(1), x(2) . . . . . x(N - m + 1) to construct, for each i, 1 -< i -< N - m + 1, cim(r) = (number o f x ( j ) such that d[x(i),x(j)] -< r ) / ( N - m + 1). We must define d(x(i),x(j)] for vectors x(i) and x(j). We follow the Takens modification o f the formula as given in ref. [5] by defining d[x,x*] = maxlu(a) - u*(a)l, a

where the u(a) are the m scalar components o f x. d represents the distance between the vectors x(i) and x(j), given by the m a x i m u m difference in their respective scalar components.

338 Journal of Clinical Monitoring Vol 7 No 4 October 1991

STEP 5 .

Next,

define N-m+l

qbm(r) = (N - In + 1) -1

E

lnCim(r)'

i=1

where In is the natural logarithm. T h r o u g h step 5, the K-S entropy and ApEn algorithms are identical. The next step distinguishes between K-S entropy and ApEn. STEP 6 (K-S). K-S entropy = lira lira lira [(I)m(r) -- (I~m+l(r)] r - ~ m-~o N--+oo

(1)

This is an abstract formula, rather than an algorithm, since we have not indicated an explicit procedure to take the limits and have not guaranteed the existence o f a limit. Furthermore, close examination o f this formula will convince the reader that the time to compute the expression within the brackets o f equation 1 grows exponentially with m, compounding the difficulty o f approximating the limit. We can explicitly compute ApEn. STEP 6 (ApEn). We define approximate entropy by

ApEn =

(2)

qbm(r) -- O p m + i ( r )

for m and r fixed as in Step 2. Typically, we choose m = 2 or m = 3; r depends greatly on the application. This choice o f m is made to insure that the conditional probabilities defined in equation 3 (below) for fixed m and r, are reasonably estimated from the N input data points. Theoretic calculations indicate that reasonable estimates o f these probabilities, for fixed m and r chosen as discussed below, are achieved with between 10TM and 30 m points, analogous to the result noted earlier by Wolf et al [12]. These values o f m give computation times for ApEn that range from 5 seconds (VAX 8800) to 10 minutes (MacIntosh PC) for data length N = 1,000. Noise much smaller than the value r is filtered out. As an important observation that provides insight into these two statistics, on unravelling definitions, we deduce that - A p E n = ~m+l(r) - ~ ( r ) = average over i ofln

CONSTRAINTS

[conditional probability that lu(j + m) u(i + m) l --<-r, given that lu(j + k) u(i + k)l --< r for k = O, 1. . . . .

m - 11.

main close on next incremental comparisons. The cim(r)s measure within a tolerance r the regularity, or frequency, o f patterns similar to a given pattern. @m+~(r) - ~m(r) represents the average stability o f these similar patterns on incrementing. Also note that while entropy is badly compromised by steady, tiny amounts o f noise, both it and ApEn are very stable to infrequent, large outliers or numeric artifacts. This follows from equation 3, vchere an infrequent outlier is a low probability event, and hence makes a tiny contribution to these measures. The magnitude does not affect the computations, beyond the filter level. This contrasts vividly with m o m e n t statistics, which are robust to (barely affected by) steady, small noise, yet dramatically affected by large outliers. As noted above, in using ApEn, we compromise some theoretic properties to gain computational efficacy and robustness to noise. Nonetheless, ApEn and entropy share some properties. Entropy is nonnegative, and a positive value implies chaos [15]. Entropy is 0 for any eventually periodic system and unchanged under translation or scaling. Intuitively, entropy measures regularity, in that if sequence A has larger entropy than sequence B, it is less ordered, more " r a n d o m " over time. ApEn is also nonnegative, although it can be positive for a nonchaotic sequence, even for an appropriately constructed periodic sequence. The loss o f this property may be too great a price to pay in some theoretical applications. ApEn is also unchanged under translation or scaling. The regularity intuition appears to carry over to ApEn. In theoretic analyses we have performed, whenever entropy(A) --<- entropy(B) for noiseless systems, then ApEn(A) ApEn(B). Furthermore, ApEn(m,r) can be defined as a limit as N goes to infinity. For a stochastic process, K-S entropy is often infinite while ApEn(m,r) is finite; thus, ApEn(m,r) can provide useful system information to distinguish differing stochastic processes. The goal o f ApEn is to provide a widely applicable statistic for the data analyst with no expertise in entropy or chaos that will distinguish data sets in terms o f their amounts o f regularity. This goal appears to be compatible with the general clinical need to distinguish healthy from abnormal.

(3)

Heuristically, entropy and ApEn measure the (logarithmic) likelihood that runs o f patterns that are close re-

The input data for ApEn is a scalar time series with measurements taken at equally spaced times. A quite similar development, although somewhat more tedious in detail, can be made for unequally spaced measurements; also, a similar development can be made for vec-

Algorithm: Pincus et ah Regularity Statisticfor Data 339

tor time series. The time series will typically contain between 50 and 5,000 numbers. Fewer than 50 numbers will likely yield a less meaningful result, especially for m = 2 or m = 3. More than 5,000 numbers will likely cause unacceptably long computer run times for the ApEn computations. With a given time series specified, there are two input parameters, m and r, that must be set; m is the length o f compared runs and r is effectively a filter level. The output o f the ApEn calculation is a single, nonnegative real number. It must be emphasized that m and r are fixed for a given application o f ApEn. More explicitly, one should consider ApEn(m,r) as the statistic. While ApEn(m,r) converges to entropy as m tends to ~ and r tends to 0, for small m and nontrivial r, the approximation may not be accurate. ApEn(2, 0.1) might be quite different from ApEn(4, 0.01), so the question arises as to which to use. We discuss choices o f m and r below, but the most important requirement is consistency. For noiseless, theoretically described systems, such as those given by the logistic map or the Henon map (two muchstudied chaotic systems), we have found that when entropy(A) -< entropy(B), then ApEn(m,r)(A) -< ApEn(m,r)(B), and conversely. Furthermore, for both theoretically described systems and those described by experimental data, we have found that when ApEn (ml,rl)(A) ~< ApEn(m~,rl)(B), then ApEn(m2,r2)(A ) --< ApEn(mz, r2)(B) (for r larger than most o f the noise), and conversely. This latter property also generally holds for parametrized systems o f stochastic processes, in which K-S entropy is infinite. We call this ability o f ApEn to preserve order a relative property. It is the key to the utility o f ApEn. There is no reason to expect that for the same systems, ApEn(ml,rl)(A ) --< ApEn(m2,rz)(B ). Since ApEn values for a given system can vary significantly with different m and r values, we do not view ApEn as an absolute measure. The power o f ApEn is its ability to compare systems. This can be accomplished with differing m and r values, which result in differing ApEn values. T o avoid a significant contribution from noise in the ApEn calculation, one must choose r larger than most o f the noise. While exact guidelines depend on the distribution o f the noise and the nature o f the underlying system, we have had success with r at least three times an estimated mean noise amplitude. Given the probabilistic nature o f the ApEn statistic, the key is that the magnitude o f the noise should rarely approach r. Noise can be categorized as measurement noise (lack o f resolution) and dynamic noise (fluctuations in system dynamics or parameter values). In the former case, one may be able to generate good estimates for the mean

and distribution o f the noise. The latter case is generally more difficult to estimate. White noise is typically assumed, although often unwisely. In some fields, such as endocrinology, estimates o f noise, e.g., intra-assay variation, are well developed [16], which facilitates a reliable choice o f r. The other consideration in choosing r is that it be large enough for the conditioning part o f equation 3 to have a sufficient number o f sequences o f x-vectors within distance r o f most specified vectors. This is to ensure reasonable estimates o f the conditional probabilities in equation 3. We have had good results with r ranging from 0.1 to 0.25 standard deviations o f the u(i) data. For smaller r values, one usually achieves numerically unstable conditional probability estimates in equation 3, given approximately 1,000 data points, while for larger r values, too much detailed system information is lost due to filter coarseness. If the noise is o f this order, use o f ApEn may be seriously compromised. In the pilot study, the ApEn statistics with m = 2 and m = 3 effectively discriminated between the healthy and sick subgroups. For very tow dimensional underlying systems, ApEn with m = 3 and suitably small r may indeed be an accurate approximation o f entropy. Even calculations with m = 1 have shown important system distinctions in related theoretic studies. The usual issues surrounding Nyquist critical sampling frequency apply in this setting. One cannot discern regularity on smaller levels than the sampling frequency o f the u(i) data. EVALUATIONOF PERFORMANCE:CLINICALPILOTSTUDY T o test the performance o f ApEn on clinical data, we analyzed heart rate data in an attempt to discern differences between healthy and sick infants. We hypothesized that ApEn would be smaller in sick neonates than in healthy neonates. This is suggested by cardiology data showing loss o f spectral reserve in Fourier series o f BBDs o f hearts with ventricular fibrillation present [3], and by spectral analysis studies on preterm babies [17] and in eventual sudden infant death syndrome victims [1]. These reports have produced suggestive results, but lack a clear-cut statistic that summarizes the frequency spectra or underlying system structure. We also hypothesized that ApEn increases with recovery from severe illness.

Methods for Pilot Study Approval for the study was granted by the H u m a n Investigation Committee o f the Yale University School

340 Journal of Clinical Monitoring VoI 7 No 4 October 1991

o f Medicine. The sample population for this study consisted o f 24 infants being cared for in the Neonatal Special Care Unit, Y a l e - N e w Haven Hospital. The sick subgroup (n = 9) consisted o f all neonates in the unit within a 2-week period who were severely ill, as defined by the clinical team. The sick infants included cases o f severe respiratory distress, persistent pulmonary hypertension (PPH), tracheoesophageal fistula with aspiration, and heart failure. The healthy neonates (n = 15) were randomly selected on discharge from the nursery in the same time period. There were no significant differences between the sick and healthy groups in postconceptional age (32.9 +- 6.1 wk vs 35.6 + 4.7 wk) and weight (2,000 -+ 1,010 g vs 2,210 + 1,000 g). Data were also collected serially on one o f the sick neonates as he recovered from PPH; 7 data sets were recorded during the first postnatal week. In addition, data from 1 fetus with sinusoidal heart rhythm (SHR) was obtained and analyzed separately. For each subject, beat-to-beat heart rate data were recorded for 90 minutes during quiet sleep to minimize movement artifacts. R-R interval lengths were derived from these raw data, based on analog voltage measurements taken every 2 milliseconds from a standard ECG monitor (Hewlett-Packard Model 78333, Palo Alto, CA). Every 5 seconds, the R-R interval lengths were averaged and recorded by a microcomputer, in units o f beats per minute. Use o f 5-second blocks was based on spectral data [1,3,17]. A time series o f 1,000 consecutive averages formed the time series for ApEn evaluation. From the same series, the beat-to-beat variability (VAR, standard deviation) was also calculated. We used analysis values o f m = 2, r = 2.0, and N = 1,000, with these variables as defined in the Constraints section. For the stated analysis parameters, each A p E n calculation consumed approximately 90 seconds o f run time on a Macintosh Plus Personal Computer. We tested for statistically significant differences between the sick and healthy groups using a Student's t-test with u n k n o w n means and variances; results are given as mean -+ 1 SD.

Results of Pilot Study The results confirmed the hypotheses. ApEn(2, 2.0) was significantly lower for the sick group, 0.80 +- 0.31 versus 1.22 -+ 0.12 for the healthy group, p = 0.003. VAR showed a difference between the two groups, although not with the same level o f confidence (8.48 -+ 2.42 vs 10.27 -+ 3.33, p = 0.143). There was no correlation o f ApEn with either postconceptional age or weight. Three deaths occurred in the sick group and none occurred in the healthy group.

1.5 (F)

=(G)

PRELIMINARYHEALTHYRANGE J.

1,0

(E)

(D)

•[c) 0.5

(A)= i(B)

i

i__

0

20

,

i

i

i

J

i

40

60

80

100

120

HOURS

POST-DELIVERY

Fig 1. Serial study of ApEn with recovering infant (birthweight, 1,090 g; gestational age, 27 wk). (A and B) PPH, on tolazotine and dopamine, ventilated; (C) on tolazoline and dopamine, minimally ventilated, resolving PPH; (D) on dopamine only, no PPH detectible, still ventilated; (E) off dopamine; (F) on CPAP; and (G) in room air. Results o f the serial measurements o f ApEn(2, 2.0) in 1 child are illustrated in Figure 1. This study o f an infant with P P H (birthweight, 1,090 g; gestational age, 27 wk) showed a large, steady increase in ApEn with recovery. The fetus with SHR gave an ApEn(2, 2.0) value o f 0.37 and a VAR o f 7.2. He died during labor.

Discussion of Pilot Study The sick infants in this study were identified better by ApEn than by beat-to-beat variability. A preliminary healthy ApEn(2, 2.0) range o f 1.0 to 1.4 (mean -+ 2 SD) is suggested by the data. The serial study result (Fig 1) is consistent with the hypothesis that ApEn is a marker o f general (neonatal) well-being. The use o f dopamine in this child would be expected to increase sympathetic tone and decrease heart rate variability. The effects o f chronotropes, inotropes, and more general drugs on regularity and, thus, ApEn are unknown and will require further study. Figure 2 can provide an intuition for the capability o f ApEn to discern subtle features o f ECGs that might indicate lack o f well-being. The two tracings, both taken from infants in the pilot study, have a similar amount o f variability, yet the one corresponding to the healthy neonate has twice as large a value o f ApEn(2, 2.0) and so has less regularity. As another example, the fetus with SHR had very low ApEn, yet nearly normal VAR. Beat-to-beat variability can detect deviations from an average but can give no information about regularity, and by itself cannot identify SHR as an ominous tracing. This study suggests the potential o f ApEn as a statis-

Algorithm: Pincus et al: Regularity Statistic for Data

341

Table 1. ApEn Statistics for Various Run Lengths and Filter Levels: Pilot Study Data 180 t

0

50

100

150

200

250

18of

~160

ed140[ 120[,

A. SICK TRACING 0

50

ApEn= 0.69

l I A R = 10.49 ,

100 150 5-SECOND INTERVALS

200

280

Fig 2. Comparison of E C G tracings with similar VAR. The upper tracing is from a healthy neonate; the lower tracing is fiom an infant with congestive heart failure. tic that c o u l d enhance p r e d i c t i v e and d i a g n o s t i c capabilities for the ill n e w b o r n . T h e s e results are p r e l i m i n a r y ; for clinical use, h e a l t h y ranges o f A p E n m u s t still be established, b a s e d o n variables such as age, general act i v i t y level, a n d sleep state d u r i n g m e a s u r e m e n t .

ma

rb

~H~

SDHd

IXse

SDs f

Probability Valueg

1 1 1 1 1 2 2 2 2 2 3 3 3 3 3

1.0 2.0 3.0 4.0 5.0 1.0 2.0 3.0 4.0 5.0 1.0 2.0 3.0 4.0 5.0

1.945 1.540 1.255 1.045 0. 882 1.303 1.221 1.063 0.914 0.784 O.657 0.838 0.820 0.759 0.682

0.219 0.242 0.248 0.241 0. 233 0.095 0.123 0.173 0.196 0.202 O.202 0.100 0.083 0.123 0.144

1.309 0.971 0.762 0.625 0. 526 1.018 0.799 0.639 0.530 0.452 O.720 0.655 0.552 0.466 0.402

0.451 0.393 0.331 0.281 0. 244 0.173 0.308 0.285 0.254 0.218 O. 145 0.210 0.223 0.214 0.198

0.0025 0.0022 0.0018 0.0020 0. 0028 0.0221 0.0032 0.0019 0.0017 0.0022 O.3922 0.0364 0.0067 0.0030 0.0027

aRun length. bFilter level. CMean ApEn value for healthy group. dStandard deviation of ApEn values for healthy group. ~Mean ApEn value for sick group. fStandard deviation of ApEn values for sick group. gProbability that sick and healthy groups have the same ApEn mean values (t-test).

DISCUSSIONOF ApEn

Relative, N o t Absolute, Measure In the C o n s t r a i n t s section, w e n o t e d that A p E n values for a g i v e n s y s t e m can v a r y significantly w i t h different m and r values. W e h a v e also c l a i m e d that the utility o f A p E n is as a r e l a t i v e m e a s u r e ; for fixed m and r, A p E n can d i s t i n g u i s h g r o u p s o f data. W e e l a b o r a t e o n this i m p o r t a n t p o i n t b y l o o k i n g m o r e d e e p l y into the data analysis for the p i l o t s t u d y . In the results o f the p i l o t s t u d y discussed a b o v e , m = 2 and r = 2.0. W e recalculated A p E n values for each n e o n a t e in the s t u d y for a range o f m and r; w e c o n s i d e r e d all possible c o m b i n a t i o n s o f m = 1, m = 2, a n d r e = 3 w i t h r = 1.0, r = 2.0, 3.0, 4.0, a n d 5.0, y i e l d i n g a total o f 15 A p E n ( m , r ) statistics. F o r each m and r, w e also calculated s u m m a r y statistics for A p E n for the sick and h e a l t h y g r o u p s : m e a n a n d s t a n d a r d d e v i a t i o n and a p value for the p r o b a b i l i t y that these h e a l t h y and sick neonates arise f r o m a single p o o l o f infants (identical m e a n A p E n values). T h e s e s u m m a r y statistics, for each m and r, are g i v e n in T a b l e 1. F u r t h e r m o r e , for three r e p r e s e n t a t i v e pairs o f m a n d r values (m = 1 and r = 1.0, m = 2 and r = 2.0, m = 3 and r = 4.0), w e give the A p E n values for each h e a l t h y a n d sick n e o n a t e in T a b l e 2.

T h e d e g r e e to w h i c h A p E n can vary, as a function o f m a n d r, is a p p a r e n t in T a b l e 2. F o r neonate H1, A p E n ( 1 , 1 . 0 ) = 1.990, A p E n ( 2 , 2 . 0 ) = 1.268, and A p E n ( 3 , 4 . 0 ) = 0,825. O n e notes the s a m e level o f distinction in the m e a n A p E n values for these three pairs o f values: the m e a n A p E n value for m = 1 and r = 1.0 is 1.945, the value for m = 2 a n d r = 2.0 is 1.220, and the value for m = 3 a n d r = 4.0 is 0.759. N o t e , h o w ever, that the p r o b a b i l i t y values that d i s t i n g u i s h the h e a l t h y a n d sick g r o u p s , via A p E n , are all v e r y significant for these three pairs: p = 0.0025, p = 0.0032, and p = 0.0030, respectively. Indeed, f r o m T a b l e 1, one deduces that A p E n is greater for the h e a l t h y neonates than the sick neonates, at the p < 0.05 significance level, for all b u t one o f the pairs o f m and r u n d e r c o n s i d e r ation, i n d e p e n d e n t o f the disparate set o f absolute A p E n values. T h e results f r o m T a b l e 1 are even m o r e d r a m a t i c i f one s i m p l y considers the A p E n calculations w i t h m = 1 and m = 2. F o r the 10 such (m,r) pairs considered here, 9 p r o d u c e a significant A p E n d i s t i n c t i o n b e t w e e n the h e a l t h y and sick g r o u p s at t h e p < 0.005 level, w h i l e the 10th (m = 2, r = 1.0) p r o d u c e s a p r o b a b i l i t y value o f 0.0221. T h e r e is one o b v i o u s l y a n o m a l o u s result in T a b l e 1;

342 Journal of Clinical Monitoring Vol 7 No 4 October 1991

Table 2. ApEn Valuesfor Three (m,r) Pairs: Pilot Study Subject~ m = 1,r = 1.0 m = 2,r = 2.0

m = 3,r = 4.0

H1 H2 H3 H4 H5 H6 H7 H8 H9 H10 Hll H12 H13 H14 H15 S1 $2 $3 $4 $5 $6 $7 $8 $9

0.825 0.791 0.910 0.646 0.502 0.941 0.579 0.726 0.751 0.642 0.836 0.773 0.827 0.764 0.879 0.187 0.314 0.561 0.193 0.332 0.765 0.593 0.558 0.690

1.990 2.001 2.172 1.776 1.508 2.185 1.658 1.856 1.902 1.675 2.125 2.278 2.078 1.897 2.068 0.897 1.061 1.596 0.563 t.255 1.857 1.061 1.667 1.821

1.268 1.296 1.378 1.122 0.935 1.403 1.076 1.158 1.243 1.114 1.285 1.253 1.270 1.201 1.303 0.457 0.585 1.016 0.322 0.688 1.186 0.731 1.028 1.149

~H from healthy group; S from sick group. for r = 1.0 and m = 3, the mean ApEn for the healthy group (0.657) is less than the mean value for the sick group (0.720). The probable cause o f this result is that the conditional probabilities calculated in equation 3 are poorly estimated. For comparisons with this choice o f m and r, we are requiring similar blocks o f data o f lengths 3 and 4 that are less than 1 beat/rain apart at each corresponding point in the block. The input o f 1,000 data points is likely not producing sufficient n u m bers o f similar sequences to ensure reasonably accurate estimates o f the desired probabilities. We anticipate that with 30,000 data points (which would require 40 hours o f stationary data and would thus be impractical), the result for this choice o f m and r would be similar to the other results in Table 1. These results reiterate the need for a statistic that can quantify a measure o f regularity, in a sound manner, for 1,000 time-series data values. Recall that for fixed rn, we noted that theoretic calculations afforded reasonable estimates o f the conditional probabilities in equation 3, and hence o f ApEn, for between 10 m and 30 m points. We r e c o m m e n d , as a preliminary choice, that m = 2 be used in the A p E n statistic definition for 1,000

points, especially when accompanied by a relatively small input r value. It thus becomes clear that to compare A p E n values, to establish normal ranges, and so forth, it is imperative to report the m and r values chosen. T o determine definitive normal A p E n values for heart rate applications, a controlled resting state, such as quiet sleep, should also likely be used for the measurement. Another important practical point is that each o f the above A p E n calculations, for 1,000 input data points, consumed about 90 seconds o f Macintosh Plus Personal C o m p u t e r time. It seems essential, for an algorithm to have wide utility, that an implemented computer prog r a m could be run in "real time" on a personal c o m puter or on a minicomputer. M a n y o f the theoretic chaos papers report numeric calculations that could only be performed on a supercomputer that devoted a significant amount o f time to the calculations; such hardware requirements would be unrealizable in most practical environments.

Potential Applications Potential applications o f A p E n are myriad. We list several medical possibilities; potential nonmedical applications o f A p E n abound, but are separate f r o m the focus o f this report. A p E n can be computed for E C G and respiratory data, perhaps providing insights into adult sudden death, sudden infant death syndrome, subtle arrhythmias, and fetal distress. ApEn m a y also be used for interpreting electroencephalographic data in neurologic and neuropharmacologic investigations. It could be used as a figure o f merit to evaluate appropriate dosage levels o f (single--opiate) anesthetics. A p E n could be used within endocrinology to characterize normal and abnormal pulsatility o f hormones [18]. For example, it m a y detect abnormal pulsatile secretion patterns for insulin before the onset of diabetes. It could also be used to test the efficacy o f estrogen replacement therapy in perimenopausal women, m i n imizing " h o t flashes." CAVEATSAND CONCLUSIONS The principal contribution o f this paper is the description o f a new, versatile, readily usable statistic, ApEn, to quantify the a m o u n t o f regularity in time-series data. The potential uses o f A p E n to provide new insights, in both medical and nonmedical settings, are rife. One such application has been discussed: the use o f ApEn on (neonatal) heart rate data. Unfortunately, potential misuses o f A p E n are also rife. We conclude b y pointing out some pitfalls to be avoided.

Algorithm: Pincus et al: Regularity Statisticfor Data 343

First, if system noise is very large (consistent signalto-noise ratio o f less than 3), the validity o f ApEn calculations may be seriously compromised. Also, for specified m (run length) in the ApEn calculation, it is important to have at least 10m, perferably 30 m data points to analyze. Otherwise, the conditional probability estimates that form the backbone o f the ApEn calculation may not be accurate. Furthermore, as emphasized in the discussion o f relative versus absolute measure, the user must fix m and r for the application o f ApEn. ApEn is a regularity, not a magnitude, statistic. While it affords a new view to traditional data analysis, it does not replace m o m e n t statistics, such as the mean and standard deviation. The series 1,0,1,0,1,0,1,0 . . . . is different from 100,2,100,2,100,2,100,2 . . . . . yet both have entropy and ApEn values o f 0. As such, we recommend use o f ApEn in conjunction with other statistics, not as a sole indicator o f system characteristics. We anticipate that the primary use o f ApEn will be to uncover abnormalities in long-term data that are not otherwise apparent. We do not anticipate that it will be so useful in discerning acute developments that are associated with relatively infrequent data measurements. This again follows from the probabilistic nature o f equation 3 and the fact that infrequent or few aberrant data will not dramatically alter ApEn. A priori, there is no reason to claim that for all systems, higher ApEn values are either good or bad. A reasonable goal is to establish normal ApEn ranges for any system under study. Each physical system will probably have a unique "healthy" or normal range, even for set m and r. Finally, ApEn decrease may correlate with standard deviation decrease (although exact entropy will not correlate in such a manner). A statistic such as A p E n / standard deviation may ultimately provide a better, normalized measure o f " p u r e " regularity change in situations in which test and control groups have significantly different standard deviations. Presented in part at the Society for Pediatric Research, Washington, DC, May 1989. APPENDIXA

Entropy and A p E n Example Calculations In this Appendix, we discuss several entropy and ApEn calculations to provide greater intuition for these measures. In calculations with series 1 through 4, we illustrate that greater entropy corresponds to the intuition o f more randomness, tn each o f these series, we consider

sequences o f integers, always either 0 or 1. The appearance o f an asterisk (*) in the series means that the integer chosen for that event has been randomly generated, with probability 1/2 o f either 0 or 1. Series 1 is given by 1,0,1,0,1,0,1,0,1,0,1,0,1,0, 1,0,1,0, . . . Alternating terms are set to 1 and 0. Series 2 is given by *,0,1,%0,1,*,0,1,*,0,1,%0, 1 , * , 0 , 1 . . . The second term in every group o f three is set to 0, the third is set to 1, and the first term in each group o f three is randomly 0 or 1. Series 3 is given by *,*,0,*,*,1,*,*,0,*,*,1,*,*,0, *,*, 1 . . . The third term in every group o f three is set, alternating 0 and 1 values; the first and second terms in each group o f three are randomly chosen as 0 or 1. Series 4 is given by *, ~, *, *, ~, ~, ~, *,*, *, *,*,*,*, *, *, *, * . . . All terms in this series are randomly chosen as 0 or 1, with probability 1/2 o f either. The K-S entropy o f series 1 is 0, o f series 2 is 1/3 (ln 2), o f series 3 is 2/3 (In 2), and o f series 4 is In 2. These calculations are relatively straightforward from the definition o f entropy as given by equation 1 in the main text. These four series are progressively more random, with the first series perfectly ordered and deterministic, the second series random for one third o f its terms, the third series random for two thirds o f its terms, and the fourth series entirely randomly generated. In calculations with series 5, we illustrate the point that to discern the full nature o f a pattern may often require larger values o f m in the ApEn statistic. Series 5 is given by 1,2,1,3,1,2,1,3,1,2,1,3,1,2,1,3, 1 , 2 , 1 , 3 . . . For r < 1.0 and m = 1, ApEn for this series is approximately 1/2(ln2), with convergence as the n u m b e r o f terms in this sequence goes to infinity. For r < 1.0 and m -> 2, ApEn for this series equals zero (and hence the entropy o f this series also equals zero). The periodicity o f this series is discerned only by comparisons o f runs o f at least two consecutive points. In calculations with series 6, we illustrate that entropy quantifies order and randomness, not necessarily magnitude. Series 6 is given by * , * , * , * , * , * , * , * , * , * , * , * , * , * , * , *,*,* . . . . All terms in series 6 are randomly chosen as 0.002 or - 0.001, with probability 1/2 o f either. The K-S entropy for series 6 is In 2; the calculation is virtually identical to that for series 4; merely consider r < 0.001 in equation 1 in the main text and note the natural correspondence o f the two series. The important conclusion here is that we could replace 0.002 by 0.00004 and - 0 . 0 0 1 by 0.0000005 and deduce the same entropy value, so long as the two values occur randomly with probability 1/2 each. Entropy is similarly unaltered by scale reduction or enlargement applied uniformly to all terms.

344 Journal of Clinical Monitoring Vol 7 No 4 October 1991

In calculations with series 7, we illustrate the crucial property that the K-S entropy formula is discontinuous to the effects o f even very small noise. Series 7 is given by 1 + * , * , 1 + * , * , 1 + * , * , 1 + * , *,1 + * , * , t + * , % 1 + *,*,1 + * , * , 1 + % * . . . All * terms in series 7 are randomly chosen as 0.002 or - 0 . 0 0 i , with probability 1/2 o f either. Note that series 7 is the sum o f the two series, 1 and 6. Also note that series 7 has values very close to those for series 1. A particular realization o f series 7 choices for the random *s) could commence: 1.002, 0.002, 0.999, 0.002, 0.999, -0.001, . . In actual data, this could represent 1,0,1,0,1,0 . . . . disturbed by some tiny system noise. The analyst who graphs series 1 and 7 notes that the respective timeseries plots o f these series are virtually coincident. Yet, the K-S entropy o f series 7 = In 2, just as for series 4 and series 6; it is not close to the entropy value o f 0 for series 1. Unlike m o m e n t statistics (e.g., mean and standard deviation), where small changes in the sequences result in small changes in the statistic, entropy is not continuous. The tiniest system noise can dramatically alter the entropy. This consideration is central to the raison d'etre for ApEn. ApEn is not discontinuous in this fashion, given that the level, r, is chosen larger than most o f the system noise. For series 7, for 0.995 > r > 0.005, and for any m, ApEn = 0. APPENDIX B

C C C

200 C C C C C

C C C C 900

901

902

PROGRAM TO COMPUTE APPROXIMATE E N T R O P Y OF T I M E SERIES D I M E N S I O N U(5000), IC(5000), ID(5000)

C C C C

c C C C C C

310

311

U S E R I N P U T , U N I T 9 IS I N T E R A C T I V E UNIT WRITE(9,900) F O R M A T (1X, ' E N T E R M, R U N L E N G T H S T O BE C O M P A R E D ' , IX) READ(9,*)M WRITE(9,901) F O R M A T (1X, ' E N T E R R, T H E FILTER LEVEL', 1X) READ(9,*)R WRITE(9,902) F O R M A T (1X, ' E N T E R N, N U M B E R OF INP U T D A T A P O I N T S ' , IX) READ (9, *)N

O P E N ( U N I T = 15, FILE = ' I N P U T ' ) O P E N ( U N I T = 16, FILE = ' O U T P U T ' ) R E W I N D 15 R E W I N D 16 D O 2 0 0 I = 1,N READ(15,*)U(I) CONTINUE C O M M E N C E ApEn C A L C U L A T I O N : INITIALIZE " V E C T O R C L O S E N E S S " COUNTERS NSUM = N - M D O 300 I = 1 , N S U M ID(I) = 0 IC(I) = 0

Fortran Program f o r A p E n

C C

O P E N t / O FILES, READ IN TIME SERIES

320 301 300 C C C C C C C C C C

C A L C U L A T E , VIA N E S T E D D O L O O P , N U M B E R OF V E C T O R S T H A T (i) ARE C L O S E T O G E T H E R , T H E

IC(i)'s, A N D (ii) R E M A I N C L O S E T O G E T H E R UPON INCREMENTING, THE ID(I)'S (NEXT COORDINATE COMPARISON) D O 301J = 1 , N S U M K=I CONTINUE DIFF = ABS(U(I + K - 1 ) - U(J + K - 1 ) ) IF (DIFF. GT. R) G O T O 320 K=K+I IF (K.GT.M) G O T O 311 G O T O 310 IC(I) = IC(I) + 1 DIFF = ABS(U(I + M) - U(J + M)) IF (DIFF. GT. R) G O T O 320 ID(I) = ID(I) + 1 CONTINUE CONTINUE CONTINUE C O M P U T E ApEn F R O M T H E ID(I)'S A N D IC(I)'S (i) T A K E A N A V E R A G E OF LOGS OF C O N D I T I O N A L PROBABILITIES OF " S T A Y I N G CLOSE, G I V E N A L R E A D Y CLOSE," AS M E A S U R E D BY ID(.) A N D IC(.) (ii) THIS GIVES PHI (M) (R) - PHI (M + 1) (R), E Q U A T I O N 2 IN T E X T

Algorithm: Pincus et al: Regularity Statistic for Data 345

400

401

A P E N = 0.0 D O 400 I = 1 , N S U M RATIO = FLOAT(ID(I))/FLOAT(IC(I)) APEN = APEN + (ALOG(RATIO)) CONTINUE APEN = - 1.0*APEN/FLOAT(NSUM) WRITE(16,401)M, R , A P E N F O R M A T ( 1 X , ' A p E n ( ' , I I , ' , ' , F 6 . 3 , ' ) = ', 1X,

F8.4, ~X) STOP END

GLOSSARY CHAOS A rapidly developing science devoted to understanding certain nonlinear dynamic phenomena and their occurrences in a wide variety o f physical and mathematic contexts. The interdisciplinary effort offers a mathematic description o f behavior that was previously classified as r a n d o m or irregular. CONmTIONAL WO~ABIUTY The probability o f the occurrence o f an event, conditioned on (given) the occurrence o f another prescribed event. DISCONTINUOUS (TO SYSTEM NOISE) A function is continuous if its value for prescribed input is the same as its limiting value for very tiny changes to the input, as the size o f the changes goes to 0. A function is discontinuous (e. g., entropy) if it is not continuous. ENTROVY, K-S ENTROVY Entropy is a concept that addresses system randomness and predictability, with greater entropy associated with m o r e randomness and less system order. There are m a n y different entropy formulations, in physics, information theory, and other branches o f mathematics, both o f a metric (distance) and a probabilistic nature. Most o f the entropy formulas contain a logarithmic weighting. Confusingly, not all entropy definitions can be related to each other. K-S entropy is a particular entropy measure, developed b y K o l m o g o r o v and Sinai, that allows one to classify dynamic systems b y rates o f information generation. NOISE The difference between the "true" and measured values o f a system, which can be due to lack of resolution, fluctuations in system dynamics or parameter values (dynamic noise), or measurement tool imprecision. REGULAR~TY The tendency that patterns within data recur in the same manner throughout the data. SeECTl~AL RESERVE A term coined by A. Goldberger to describe a property o f a Fourier povcer spectrum. Greater spectral reserve corresponds to a wider, more broadband spectrum, while a loss o f spectral reserve corresponds to narrowing o f a broadband spectrum.

STRANGE ATTRACTOR A geometric object, defining the "equilibrium" behavior o f a deterministic nonlinear system that is bounded, yet not periodic or multiply periodic. These objects can be low dimensional, and typically have a self-similar, or fractal, structure.

REFERENCES 1. Kluge KA, Harper RM, Schechtman VL, et aL Spectral analysis assessment of respiratory sinus arrhythmia in normal infants and infants who subsequently died of sudden infant death syndrome. Pediatr Res 1988;24:677-682 2. Goldberger AL, Bhargava V, West BJ, Mandell AJ. On a mechanism of cardiac instability: the fractal hypothesis. Biophys J 1985;48:525-528 3. Goldberger AL, West BJ. Fractats in physiology and medicine. Yale J Biol Med 1987;60:421-435 4. Grassberger P, Procaccia I. Estimation of the Kolmogorov entropy from a chaotic signal. Phys Rev A 1983; 28:2591-2593 5. Takens F. Invariants related to dimension and entropy. In: Atas do 13. Rio deJaneiro: Col. Brasiliero de Matematicas, 1983 6. Santen RJ, Bardin CW, Episodic luteinizing hormone secretion in man. Pulse analysis, clinical interpretation, physiologic mechanisms. J Clin Invest 1973;52:2617-2628 7. Van Cauter E. Quantitative methods for the analysis of circadian and episodic hormone fluctuations. In: Van Cauter E, Copinschi G, eds. Human pituitary hormones: circadian and episodic variations. The Hague: Martinus Nijhoff, 1981:1-25 8. Feynman RP. The Feynman lectures on physics, vol. 1. Reading, MA: Addison-Wesley, 1963:44.10-44.13 9. Billingsley P. Ergodic theory and information. New" York: Wiley, 1965:60-94 10. Eckmann JP, Ruelle D. Ergodic theory of chaos and strange attractors. Rev Mod Phys 1985;57:617-656 11. Crutchfield JP, Packard NH. Symbolic dynamics of onedimensional maps: entropies, finite precursor, and noise. Int J Theor Phys 1982;21:433-465 12. Wolf A, SwiftJB, Swinney HL, VastanoJA. Determining Lyapunov exponents from a time-series. Physica 1985; 16D:285-317 13. Greenside HS, Wolf A, Swift J, Pignataro T. The impracticality of a box-counting algorithm for calculating the dimensionality of strange attractors. Phys Rev A 1982; 25:3453-3456 14. Fraser AM. Information and entropy in strange attractors. IEEE Trans IT 1989;35(2):245-262 15. Ruelle D. An inequality for the entropy of differentiable maps. Bol Soc Bras Mat 1978;9:83 16. Lambalk CB, de Koning J, van Kessel H, et al. Calculation of the intra-assay variation per assay and its relevance to luteinizing hormone pulse detection. In: Crowley WF, Hofler JG, eds. The episodic secretion of hormones. New York: Wiley, 1987:67-73 17. Aarimaa T, Oja R, Antila K, Valimaki I. Interaction of heart rate and respiration in newborn babies. Pediatr Res 1988;24:745-750 18. Merriam GR. Methods of characterizing episodic hormone secretion. In: Crowley WF, Hofler JG, eds. The episodic secretion of hormones. New York: Wiley, 1987:47-65