Abstract
We describe an analysis of dynamic behavior apparent in timesseries recordings of infant breathing during sleep. Three principal techniques were used: estimation of correlation dimension, surrogate data analysis, and reduced linear (autoregressive) modeling (RARM). Correlation dimension can be used to quantify the complexity of time series and has been applied to a variety of physiological and biological measurements. However, the methods most commonly used to estimate correlation dimension suffer from some technical problems that can produce misleading results if not correctly applied. We used a new technique of estimating correlation dimension that has fewer problems. We tested the significance of dimension estimates by comparing estimates with artificial data sets (surrogate data). On the basis of the analysis, we conclude that the dynamics of infant breathing during quiet sleep can best be described as a nonlinear dynamic system with largescale, lowdimensional and smallscale, highdimensional behavior; more specifically, a noisedriven nonlinear system with a twodimensional periodic orbit. Using our RARM technique, we identified the second period as cyclic amplitude modulation of the same period as periodic breathing. We conclude that our data are consistent with respiration being chaotic.
 control of breathing
 periodic breathing
 dynamical systems
 chaos
 surrogate data techniques
this paper describes and summarizes a study of infant breathing by using dataanalysis techniques derived from dynamic systems theory (DST). Such techniques have been useful for examining other complex physiological rhythms such as heart rate (7, 46, 47, 63,67), electroencephalogram (2, 24, 26, 29, 36, 37), parathyroid hormone secretion (31), and opticokinetic nystagmus (45) and can distinguish variations that are random from those that are deterministic (“chaotic”).
Any model that describes the control of respiration must account for two important behaviors, rhythm and pattern generation. There are two distinct approaches to the problem of rhythm: the first requires the existence of discrete “pacemaker” cells with intrinsic activity that drives other respiratory neurons. The output of various respiratory centers or pools of motoneurons is then organized by a pattern generator. A second approach implies that networks of cells with oscillatory behavior interact in a complex way to produce respiratory rhythms, which are either further organized by a pattern generator or might be selforganizing (11). The first model is, in effect, a linearly coupled model with pattern generation dependent on feedback mechanisms from the various components of the respiratory system, whereas the second is more likely to be complex. Advances in neurobiology have allowed recordings to be made from individual neurons and groups of neurons in the brain. Using these techniques, various studies have demonstrated that the concept of discrete respiratory centers made up of neurons with specific functions defined by the nature of a particular “center” is obsolete (11). Whereas there is organization of neurons into functional networks or pools, these are not necessarily anatomically discrete. Also, there are conflicting data in regard to the presence of a specific pattern generator. Given the complexity of the connections between the various groups of oscillating, respiratoryrelated neurons and the capacity for interactions between simple oscillating systems to produce complex behavior, we have argued that information about the organization of respiratory control can be determined by using DST. In essence, the argument that there is a simple pattern generator that coordinates the output from various respiratory centers is unnecessary if the output from interacting networks is dynamic and selforganizing. To examine such a hypothesis, respiration must first be adequately shown to be chaotic, then examined further to describe its dynamic structure.
Most studies of the dynamic behavior of biological systems have used fractal dimension estimation to establish that a system’s behavior is chaotic or to classify distinct types of behavior by their complexity. Recent studies have suggested that respiration in humans is chaotic. If that is the case, then techniques derived from DST should allow the dynamic structure of respiratory behavior to be better described, thus improving our understanding of the control of breathing. However, these earlier studies have important limitations. Most studies have used the Grassberger and Procaccia algorithm (14, 15) for estimating fractal dimension, which is simple and easy to implement. Unfortunately, it is now recognized (19, 33) that this algorithm has some technical problems that can lead to misinterpretations of data. The most serious problems occur with small data sets or when the system incorporates a substantial noise component. The studies described below each suffer these limitations to a greater or lesser extent and none examined respiratory patterns in quiet sleep when breathing is most regular.
Donaldson (9) used DST to examine respiratory behavior in humans. He studied resting adults and, by estimating Lyapunov exponents, concluded that resting respiration is chaotic. However, this relatively simple approach could not distinguish a nonlinear dynamic system from linearly filtered noise.
Pilgram et al. (30) studied the dimensional complexity in the respiratory traces of infants during rapideyemovement (REM) sleep. They concluded that breathing during REM sleep is chaotic but did not study patterns of breathing during quiet sleep. In reaching their conclusion, Pilgram et al. also employed the method of surrogate data. Dimensionestimation algorithm used by Pilgram et al. failed to give reliable estimates for their noisedriven surrogates, and so the authors inferred that the control of respiration is possibly deterministic.
In another set of studies, Sammon and colleagues (3842) examined dynamic respiratory behavior in rats. From their observations they concluded that, in anesthetized vagotomized rats, the respiratory system behaves as an oscillator with a single degree of freedom. With the vagus intact, however, respiratory behavior was more complex, exhibiting loworder chaos, which, the authors speculated, was due to feedback from various types of pulmonary afferent activity.
In this paper, we test for cyclic modulation in the amplitude of respiration. Cyclic fluctuation in the amplitude of respiration indicates the presence of additional complexity in the respiratory system and is consistent with the results of our correlation dimension calculations. Other authors (see, e.g., Ref. 64) have identified cyclic behavior in breath sizes. Waggener and colleagues (64) identify cyclic variation in the breath sizes of resting adults at simulated altitude.
The study reported here employs a new algorithm to determine fractal dimension that is technically more complex but is in practice more reliable, more robust under the restrictions of finite data, and less prone to misinterpretation. Estimates of fractal dimension were used as an aid to identify the dynamic system that produced the data we have measured.
From dimension estimations we conclude that the dynamics of breathing during quiet sleep are consistent with a largescale, lowdimensional system with a substantial smallscale, highdimensional component, i.e., a periodic orbit with a few (perhaps 2 or 3) degrees of freedom supplemented by smaller, more complex deterministic fluctuations.
The nature of the lowdimensional system was investigated further by constructing surrogate data, which enabled us to test whether the dynamics were consistent with linearly filtered noise or with a nonlinear dynamic system. Our results show clearly that, in almost all cases, the dynamics are best described as a lowdimensional, nonlinear dynamic system being driven by a highdimensional noise source. In those cases when such a model is inconsistent with the data, the measured data have strong indications of nonstationarity; that is, the breathing patterns changed during the recording (for example, a sudden switch to periodic breathing occurred).
To describe a complex dynamic system as chaotic may not be helpful when the components that produce the structure cannot be described. Understanding these components is more important than simply classifying the system as chaotic. However, before the mechanism can be understood, it is important to determine basic properties of the system: Is it dynamic? Is it nonlinear? Is it chaotic?
We, therefore, investigated the structure of the dynamics we described in greater detail. We used a new form of linear modeling to expose a possible source of our multidimensional periodic orbit. By applying a reduced autoregressive modeling (RARM) algorithm to some data from sleeping infants we observe cyclic amplitude modulation (CAM) indicative of substantial structure in breathtobreath variability. These calculations lead us to conclude that respiration in quiet sleep is most probably controlled by multidimensional oscillator supplemented by smallscale, complex (chaotic) behavior. Moreover, our results are inconsistent with a noisedriven, onedimensional model of eupnea (or any other model that does not include substantial structure in breathtobreath variability).
After a brief introduction to the new dimension estimation algorithm, we describe the experimental methodology, including a description of our surrogate datageneration methods and our RARM algorithm. Finally, we discuss the dimension calculations and the results of the hypothesis testing by using the surrogate data sets.
METHODS
Using standard noninvasive inductiveplethysmography techniques, we obtained a measurement proportional to the crosssectional area of the chest or abdomen, which is a gauge of the lung volume. The present study collected measurements of the crosssectional area of the abdomen of infants during natural sleep. The study was approved by the Princess Margaret Hospital Ethics Committee.
Subjects
Ten healthy infants were studied at 2 mo of age in the sleep laboratory at Princess Margaret Hospital. Data from 15 infants (including the 10 studied at 2 mo) between 1 and 6 mo of age were recorded and are used in our calculations to detect CAM.
Data Collection
The unfiltered analog signal from a NIMS Respitrace+ [NonInvasive Monitoring systems (NIMS); trading through Sensormedics, Yorba Linda, CA] inductance plethysmograph was passed through a directcurrent amplifier and 12bit analogtodigital converter (sampling at 50 Hz). The digital data were recorded in ASCII format by using Labdat and Anadat software packages (RHTInfodat, Montreal, Quebec, Canada) installed on an IBM 286 microcomputer. This ASCII data were then transferred to Unix workstations at the University of Western Australia for analysis with the use of Matlab and C programs.
The only practical limitation on the length of time for which data could be collected is the period during which the infant remains asleep and still. The crosssectional area of the lung varies with the position of the infant. However, in this study, we were interested only in the variation due to the breathing, and so we have been careful to avoid artifacts due to changes in position or band slippage. We have made observations of up to 2 h that are free from significant movement artifacts, although, typically, observations are in the range of 5–30 min. All 27 observations used to calculate dimension are between 240 and 360 s; those used to identify CAM are between 400 and 1,400 s.
In contrast to the study by Pilgram and colleagues (30) that examined breathing in REM sleep, we have studied infants in quiet sleep. From measurements of electroencephalogram, electromyogram, and electrooculogram, sleep stage was determined by using standard polysomnographic criteria (6). During quiet sleep, breathing often appears relatively regular. The possibly chaotic features of most interest are the small variations from this regular periodic behavior. Because we wished to observe such fine detail, we did not filter signals. The analog output of the respiratory plethysmograph (operating in its directcurrent mode) has no builtin filtering. Filtering methods, such as linear filters and singularvalue decomposition methods (30), can remove some features that we wished to observe. Furthermore, filtering (even to avoid aliasing) has been shown in some cases to lead to erroneous identification of chaos (25, 28, 34).
The abdominal signal is not necessarily proportional to lung volume. Takens’ embedding theorem (56) (and, therefore, the methods of this paper; see appendix , Embedding) only requires a diffeomorphism (a smooth function) of a measurement of the system. Moreover, present methods are not capable of dealing well with multichannel data, and, therefore, use of both rib and abdominal signals to approximate actual lung volume is difficult. Of the available measurements, we found that the abdominal cross section was the easiest to measure experimentally.
The 27 observations used to calculate dimension were selected based on sleep state (quiet stage 3–4 sleep) and then on the basis of sufficient stationarity and a minimum of 4 min in length. From each of these observations, 240 s of stationary data (the 240 s that had the most stationary moving average) were used to calculate dimension. Data for RARM calculations were selected based on being at least 400 s in length and taken during uninterrupted quiet sleep.
DATA ANALYSIS
In this study, we employed three main analysis methods: correlation dimension estimation, surrogate data analysis, and RARM. This section provides a description of these methods. The mathematical detail is deferred to appendixes a–d. First, we discuss generalized dimension and correlation dimension estimation, and then we provide an overview of surrogatedata techniques and describe RARM procedures.
Glossary
An autoregressive process of order n
C _{N}(ε)
Correlation function
d _{c}, d _{c}(ε_{0})
Correlation dimension
d _{e}, n
Embedding dimension
I(X )
An indicator function
MDL
Minimum description length
N
Length of a data set
P _{j},P _{j}(⋅)
An orthogonal projection
Z ^{+}
Set of positive integers
ε
Length scale
ε_{0}
Observation scale
τ
Embedding lag
∥⋅∥
Euclidean norm—the usual distance function
Dimension Estimation
We are accustomed to thinking of realworld objects as one, two, or three dimensional. However, there exist complex mathematical objects, called fractals, that have noninteger dimension, a socalled “fractal dimension.” Many realworld phenomena, in particular chaotic dynamic systems, can be observed to have properties of a fractal, including a noninteger dimension. A meaningful definition of fractal dimension comes from a generalization, or extension, of wellknown properties of integer dimension objects. For a more detailed discussion of generalized fractal dimension and estimation of correlation dimension d _{c}, see appendix.
Most applications of correlation dimension to physiological sciences have utilized the Grassberger and Procaccia algorithm (14,15). However, this paper employs a new algorithm, which is technically more complex but is in practice more reliable and less prone to misinterpretation. Unlike previous estimation methods, this new algorithm recognizes that the dimension of an object (its structural complexity) may vary depending on how closely you examine it. Hence, the value of the estimate of correlation dimension may change with scale. It therefore offers a more informative and accurate estimate of dimension.
Computing correlation dimension d _{c} as a function of scale d _{c}(ε_{0}) can tell us much more about the structure of an object; for example, it can indicate the presence of large scale “periodic” motion and simultaneously detect smaller scale, higher dimensional chaotic motion. Quoting a single number as the correlation dimension of a data set ignores much of this information; in many respects, it produces an “average dimension.” Plots of dimension as a function of scale are particularly important when complex physiological behavior is studied because they yield far more information than a single estimate at a fixed scale.
The estimation algorithm.
The estimation algorithm used for the calculations in this paper is described in detail by Judd (19, 20), and an alternative treatment may be found in (for example) Ref. 17. The accuracy of this algorithm has been independently verified by Galka and colleagues (13). One important advantage of the new method is that it is possible to calculate error bars for dimension estimates. The confidence intervals on the dimension that the algorithm provides are dependent on the length of the time series.
For each time series, the dimension was calculated for timedelay embedding (see appendix , Embedding) in two, three, four, five, seven, and nine dimensions, a far greater range than necessary, but one which encompasses suitable values of embedding dimension suggested by false nearest neighbor methods (appendix, Embedding dimension).
Hence, for each data set, our dimension estimation methods produced a graph with many lines on it. Each line on the graph is the dimension estimate for the same data set with a different embedding dimension. These lines are a plot of the change in correlation dimension (vertical axis) with scale (horizontal axis). Scale is calculated as the logarithm of “viewing scale,” so moving to the right on a plot indicates increasing scale. The righthand end of the plots is the estimate of dimension at the largest scale (the most obvious features), whereas the lefthand end is the dimension estimate at the smallest scale (the finest details).
Phenomenological Models
Throughout this paper, the term “model” is used in the sense of a phenomenological model and not in the sense of a mechanistic model. In this section, we briefly highlight this distinction.
Surrogatedata techniques involve (implicitly or explicitly) building a model from data, as do RARM procedures. It is important to emphasize that the types of models we are describing are phenomenological rather than mechanistic. That is, these models are not described as a set of equations that have been motivated by physical laws or physiological insight. Rather, the models we describe here are fitted to experimental data to predict future behavior of the system or to simulate the system on a computer and produce an output indistinguishable from observations of the actual system.
Usually, a phenomenological model is selected from a broad abstract class of models by fitting data using some statistical criteria, that is, to find the most probable model within the class. Nonlinear models and reduced autoregressive models are particular examples of broad abstract classes. Mechanistic models include, for example, a model of the lungs as bellows or a physical model of bloodgas concentrations. In a mechanistic model, parameters and variables represent physical quantities, for example, the elasticity of the lung or blood oxygen concentration. In a phenomenological model, this is not the case. However, phenomenological models can provide clear and useful information about a system, such as the periodicities of the system.
Surrogate Data
Estimating the dimension of the data set gave valuable information about the geometric structure of these data, but dimension estimation alone is not enough to give a sure indication of the presence of lowdimensional chaos or even nonlinear dynamics. Any experimentally obtained data will include some observational noise and, when added to a deterministic linear process, can produce dimension estimates not dissimilar to the results of our calculations.
To determine whether our results indicate the presence of anything more complicated than a noisy linear system, we employed the surrogatedata methods described by Theiler et al. (59). The principle of surrogate data is the following. One first assumes that the data come from some specific class of dynamic system, possibly fitting a parametric (phenomenological) model to the data. One then generates surrogate data from this hypothetical system and calculates various statistics of the surrogate data and original data. The surrogate data give the expected distribution of statistical values, and one can check that the original has a typical value. If the original data have atypical statistics, then we reject the hypothesis that the system that generated the original data is of the assumed class. One always begins with simple assumptions and progresses to more sophisticated models if the data are inconsistent with the surrogate data.
We expect our data to be most consistent with some type of nonlinear dynamic system. Before considering this type of surrogate, it is necessary to determine whether a simpler description of the data would be sufficient. To do this, we compared our data with surrogates generated by the traditional (linear) methods (see appendix, Traditional Surrogate Data). Many studies in biological sciences have employed these traditional surrogate methods (2, 31, 42, 63, 67). These methods determine whether experimental data are significantly different from specific (broad) categories of linear systems. In addition to these linear surrogate tests, we applied a new, more complicated, nonlinear surrogate test (50, 53). This method was used to determine whether the data are distinguishable from these generated by a broad class of nonlinear models (see appendix, NoiseDriven Nonlinear System Surrogates, andappendix ).
Reduced Autoregressive Models
In this section, we briefly describe the rationale of the final mathematical tool we utilize: RARM. We applied these methods to time series consisting only of the size of successive breaths (the difference between peak and trough values). From this we deduced the period of CAM present in the original time series. One could also examine the length of successive breaths to determine the presence of cyclic durational modulation. However, in this paper, we focus on CAM.
Our new method extends the traditional autoregressive (AR) model of order n [AR(n)], which predicts the next value in a time series as a weighted average of the last n values. We consider instead a RARM where any past values may be used to predict the upcoming value but only those that are important are used. To determine which past values are important, we employed Rissanen’s minimum description length (MDL) criterion (35) using a modeling procedure described by Judd and Mees (21, 22); our implementation of these methods has been presented elsewhere (55) and is discussed elsewhere (48, 54). A computer software implementation of these techniques is available online (49) or from the first author. A brief description of the mathematical methods is presented in appendix. For now, let us assume that RARM can produce a model consisting only of those previous values that are useful in predicting future values. This is not necessarily a particularly good model in terms of prediction—it is only an approximation to the breathtobreath dynamics that we utilize to extract important information. Using this information, we deduce the period of quasiperiodic behavior in the time series from the temporal separation of the previous values. Hence, RARM can identify the period of CAM in much the same way as autocorrelation may, except our methods prove to be more sensitive and more discriminatory. Our method involves automated application of statistical procedures to time series and can be applied even when amplitude modulation is not apparent to the eye.
By reducing the original data to a breathtobreath time series, we effectively rescale the time axes so that each breath is of equal length. However, CAM on a breathtobreath basis does not (necessarily) suppress timedependent dynamics. A hypothesized cyclic variation in breath duration may be related to a cyclic variation in breath amplitude (and, hence, evident in the breath amplitude time series). These could essentially be two separate observations of the same periodic behavior. The duration of a single breath is also far more difficult to measure accurately because of (relatively) long flat peaks and troughs.
RESULTS
We first present our results from applying our dimension estimation algorithm. Next, we describe the results of our surrogatedata and RARM calculations.
Dimension Estimation
The results of the calculations ofd _{c}(ε_{0}), as shown in Fig.1, can be summarized as follows. All calculations fall into two broad categories. Most of the estimates ofd _{c}(ε_{0}) produced curves that increase, more or less linearly, with decreasing scale log ε_{0}, but some showed an initial decrease in dimension before increasing with decreasing scale (Fig. 1, subjects 1 and4). For any particular data set, it was generally found that the graph of d _{c}(ε_{0}) was shifted to higher dimensions as the embedding dimension was increased, although the shape of the graph varied little with changes in embedding dimension. In nearly all cases, the dimension estimates at the largest scale lay between two and three.
The more or less linear increase in dimension with decreasing scale ε_{0}, and the shift to higher dimensions as the embedding dimension is increased, are both indications that the system, or measurements, have a substantial component of highdimensional dynamics, or noise, at smalltomoderate scales. The increase of dimension with decreasing scale is an obvious effect of highdimensional dynamics or noise. The shifting to higher dimensions with increasing embedding dimension occurs because in a higher dimensional embedding the points “move away” from their neighbors and tend to become equidistant from each other, which in effect amplifies, or propagates, the smallscale, highdimensional properties to large scales. (This effect is related to the counterintuitive fact that spheres in higher dimensions have most of their volume close to their surfaces rather than near their centers, as is the case in two and three dimensions.)
Some of the dimension estimates, particularly in two and three dimensions, produced curves that linearly increased for largelength scales but appeared to level off as length scale decreased. For most of the estimates we have computed, this is the case when the data are embedded in two dimensions. Furthermore, for these embeddings in twodimensional space, the correlation dimension estimate seemed to approach two. This indicates that as we look closer at the data (that is, at a smaller length scale), the data appear to fill up all of our embedding space. For many of the dimension estimates (Fig. 1,subjects 7 and 9), the embedding in three dimensions also leveled at values slightly less than three. This behavior can be attributed to an attractor with correlation dimension of ∼2.8–2.9. However, it is probably more likely that this too is simply due to the data “filling up” the threedimensional space. This is consistent with the results of our false nearest neighbor calculations, which suggested that three or fourdimensional space would be required to successfully embed the data.
There is one particular estimate that appeared to behave quite differently from all the others. Some of the curves of the estimates for subject 8 appeared to increase, decrease, and then increase again. This could indicate that, as we look more closely at the structure, there is some length scale for which the embedding structure seems to be relatively high in dimension, whereas by looking at an even smaller length scale the behavior has significantly lower dimension. These observations are supported by what we can observe directly from the data. This time series includes an episode of periodic breathing—increasing the complexity of the largescale behavior (see Fig. 2). Similarly, some of the data sets for subject 2 include large sighs, causing the dimension estimate to increase at large scales (see Fig.3).
Finally, the remainder of the estimates (e.g., Fig. 1, subjects 1, 2, 4, 6, 7, 8, and10) behaved in yet another manner. These estimates were approximately constant for a small range of large length scales and gradually increased over small length scales. The estimates at large length scales were generally about two to three, indicating that the largescale behavior is slightly above two dimensional. The increase in dimension estimate for smaller length scales can again be attributed to either noise or highdimensional dynamics. However, the scale of “smallscale structure” in the dimension estimates is at a larger scale than the instrumentation noise level. Typically, the smallest scale is log_{e}(ε_{0}) ≈ −2.5, a scale of ∼5% of the attractor (e ^{−3} ≈ 0.049787 ≈ 0.05). The digitized signal will typically use at least 10 bits of the analogtodigital converter (2^{−10} = 1/1,024 < 0.001); other sources of instrumental error are certainly at levels <5%.
The approximately twodimensional behavior is probably due to the inspirationexpiration cycle along with the breathtobreath variation within that cycle. This is easily visualized as the orbit of a point around the surface of a torus. A dimension estimate of two could indicate that the attractor was any twodimensional surface; the embedded data, however, have an approximately toroidal shape. In this motion, there are two characteristic cycles, first, the motion around the center of the torus and, second, a twisting motion around the surface. Our estimates slightly over two indicate that this behavior is complicated further by some other roughness over the surface of the torus. The shape of such an attractor would very closely resemble the textured surface of a doughnut.
SurrogateData Analysis
Dimension estimation gave information about the shape of the dynamic system we are studying. In an attempt to classify this system, we applied surrogatedata techniques. First, we compared breathing dynamics to linear systems. Following this, we compared the breathing dynamics with nonlinear dynamic systems by fitting a type of nonlinear model to the data.
Linear surrogates.
By comparing the value of dimension obtained from our data and surrogates consistent with each of these three null hypotheses, we were able to reject all three null hypotheses (see Fig.4 for an example of such a calculation). These results are summarized in Table1, which shows the number of SDs between the values of d _{c}(ε_{0}) for data and surrogate, for the value of log (ε_{0}), which gave the greatest difference. This is calculated over the range −2.5 ≤ log (ε_{0}) ≤ −0.5, and ford _{e} = 3, 4, 5. Data are from infants at 2 mo of age. The symbol n/a indicates that none of the surrogates produced convergent dimension estimate at any value of ε_{0}. For each data set and each hypothesis test, there are three pairs of numbers. These three pairs of numbers are the results ford _{e} = 3, 4, 5, respectively. The first number is the number of SDs by which the mean value of dimension for the surrogates exceeded that for the data. The second number (in parentheses) is the value of log (ε_{0}) for which this occurred.
The work of Pilgram et al. (30) with respiratory traces during REM sleep produced similar observations for a different physiological phenomena. By rejecting these null hypotheses, we may make two important observations. First, the data are not a (monotonic) transformation of linearly filtered noise. Second, correlation dimension alone is sufficient to distinguish between our data and data consistent with these hypotheses.
These results, however comforting, are not particularly surprising. Our data are regular and periodic, and the surrogates’ are not (see, e.g., Fig. 5).
Nonlinear surrogates.
For each set of data, we have calculated its correlation dimension. Using a modeling algorithm described elsewhere (21, 22, 51), we constructed a radialbasis model of the data. From this model, we constructed 30 surrogate data sets. The surrogates were embedded in two, three, and four dimensions, by using the same embedding strategy as the true data set. We then calculated the correlation dimension curve for each set of surrogate data.
The results (see Fig. 6) of these calculations fell into two very distinct categories. For many of the data sets, the surrogates very closely resembled the true dimension estimate, whereas for some others the data and the surrogates appeared to be very different. On a closer examination of the time series, it appears that the model failed to produce accurate surrogates only when the data set was significantly nonstationary. Although no data set used in these calculations had an obvious drift or changed sleep state, nonstationarity occurred with sudden change in respiratory behavior (see, e.g., Fig. 2). Hence, when the data were sufficiently stationary (as was the case with 24 of our 27 data sets), the modeling algorithm produced surrogate data that were indistinguishable (according to the method of surrogate data, with respect to correlation dimension) from the true data. Furthermore, the models exhibited a toroidal attractor, as suggested by the correlation dimension estimates.
CAM Detected by Using RARM
In this section, we present some preliminary results of the detection of CAM in the respiratory traces of infants in quiet sleep. The data used for these calculations are different from those for which the correlation dimension was calculated. The data requirements of this algorithm are moderately large (typically 10 min of recording), calculation of correlation dimension and radialbasis models for such data sets proved prohibitive. Moreover, the two types of models are entirely distinct: RARM is more robust to nonstationarity, whereas radialbasis models are better at capturing qualitative (and many quantitative) features of respiration.
Table 2 outlines the results of our calculations applied to 14 data sets from 14 infants. Data used for these calculations were recorded during quiet sleep. Subjects 1–10 are the same subjects as used for correlation dimension estimates. Data for subjects 1–6 were recorded during the same study as those used for dimension estimates. Data for subjects 7–12 were recorded at 4 mo of age, for subjects 13–14 at 6 mo. Respiratory rate is the average respiratory rate over the duration of the recording. Note that, although there was some variability in both the respiratory rate and the period expressed as a number of breaths, the period in seconds was relatively constant. In most cases, this period also falls within the range of periodic breathing. In subject 11, periodic breathing with cycle times from 13.5 to 15.5 s occurred during the same study.
DISCUSSION
This study has confirmed that apparently regular breathing during quiet sleep is probably chaotic. This conclusion should be qualified. As Rapp (33) observed, to conclude that a phenomenon is chaotic is both difficult and often irrelevant. In real data sets, noise contamination will always increase the dimensional complexity of the data, and almost any experimental data will exhibit noninteger correlation dimension. Identification of apparently chaotic behavior is, however, a good first step in dynamic analysis. Our data are the most convincing evidence that respiratory variability in humans is deterministic and not random, due to noise. Furthermore, we have extended our observations and analyses to describe the dynamic structure of the system in greater detail. We have demonstrated the presence of a twodimensional periodic system: one dimension corresponds to the regular breathbybreath cycle of respiration, the other to a regular periodic modulation of breath amplitude (CAM), with cycle length approximating that of periodic breathing.
By applying a new dimension estimation algorithm, we have produced reliable, accurate, and informative dimension estimates. This method produces error bars which, for moderately long time series (such as ours), are reassuringly small. This, however, must not be taken to indicate that we were able to estimate the dimension of the dynamic system responsible for this behavior with similar confidence. As with any experiment, the accuracy of the results is limited by physical constraints such as observational error and noise interference. Even if the data can be collected in such a way that they accurately represented the volume of the subject’s lung, there may be other phenomena contributing to that measurement, which are independent of the dynamic system we are attempting to observe. Hence, the tightness of the error bars we have calculated was due more to the quantity of the data than to lack of experimental error.
Our dimension estimate results indicate that, on a large scale, there is lowdimensional behavior, whereas the smallscale behavior is often dominated by very highdimensional dynamics or noise (that is, extremely highdimensional dynamics). Even though false nearest neighbor techniques suggest that we were embedding in high enough dimensions, there was still some smallscale behavior that filled the embedding space. The scale at which the embedding space is filled by the dynamics could indicate the level of experimental noise.
The most conclusive estimates from this study indicated that the structure of the attractor is likely to be similar to a torus with smallscale, very highdimensional dynamics. Hence, at large length scales, the structure looked like the surface of a torus, whereas at smaller length scales dimension increased. This indicates that the attractor appears to be a torus with a very rough surface. The most important conclusion from these data is that this twodimensional periodic system indicates two levels of periodicity. Hence, in addition to the periodic inspirationexpiration motion, it is likely that there was some cyclic breathtobreath variation in breath depth.
By applying the method of surrogate data we demonstrated that the correlation dimension is related to the data from which we estimate it in a nontrivial way. The surrogates produced by algorithms 0, 1, and 2 are clearly inadequate. It is apparent that they should fail, and this was confirmed by our results. We have constructed our own surrogates using a nonlinear modeling process and compared surrogates and data to test the accuracy of the model. For the majority of our data sets, we found that the data and nonlinear surrogates were indistinguishable according to correlation dimension. For those data sets that were distinguishable from their surrogates, we found that there were several possible reasons for this. Usually, if the data were nonstationary, the model simply failed to produce surrogates that were close enough to the data. The model is stationary and periodic, whereas the data are not.
Occasionally, with nonstationary data, the model failed to produce even periodic surrogates. If this was the case, then the model had a stable fixed point. In these cases, the dimension estimates of data and surrogate were obviously different, and a better model is required. The fact that this modeling algorithm failed in cases where the data were not stationary is not particularly surprising—both modeling and dimension estimation algorithms require stationarity. Perhaps, with improved modeling techniques, similar results could be obtained in these cases.
Even if both data and surrogate were stationary, the dimension estimates of the surrogates could still be different from those of the data. In all these cases, however, this has been found to be a problem with the level of dynamic noise introduced to the model to generate the surrogates. By changing the noise level, the dimension would also change, effectively moving the dimension estimate vertically. Because, in these cases, the shapes of the dimension estimate curves were approximately the same, by altering the noise level it was possible to produce surrogate estimates that were indistinguishable from the data. In all cases, however, the dynamic noise was substantially less than the model’s root mean square prediction error. The root mean square prediction error is the noise level predicted by the modeling algorithm. However, this is the total noise and includes both dynamic and observational noise.
Dynamic noise and observational noise have a different effect on the correlation dimension estimates. Observational noise will increase the value of correlation dimension at length scales less than and equal to the noise level. It appears from our calculations that an increase in the level of dynamic noise increases the correlation dimension estimate equally across all length scales, effectively producing a vertical shift in the estimate. An increase in dynamic noise will certainly have a greater affect on d _{c}(ε_{0}) for larger ε_{0} than would a similar increase in observational noise. Assuming one has correctly identified the underlying deterministic dynamics, it may be possible to “tune” the level of dynamic noise so that surrogates and data have approximately the same value of correlation dimension estimates at moderatetolarge length scales and then alter the level of observational noise to tune the dimension estimates at small length scales. Hence, the level of noise required to be unable to reject the surrogate data is an indication of the relative proportion of dynamic and observational noise in the system. That is, we can distinguish between random behavior within the system (dynamic noise) and experimental error (observational noise).
Our calculations have demonstrated that the respiration of infants in quiet sleep can be modeled well with radialbasis methods. We conclude that the behavior is similar to a noisedriven periodic orbit. Furthermore, the radialbasis surrogates confirm that the periodic orbit associated with this behavior has at least two degrees of freedom, as suggested by our initial dimension estimates. Finally, our RARM calculations to detect CAM support these results.
The observation of CAM is intriguing. We serendipitously recorded periodic breathing from one infant. The cycle time of CAM (13.3–15.6 s) in the same infant corresponded almost exactly with that of the observed periodic breathing (13.5–15.5 s) as demonstrated in Fig. 2. The relationship to periodic breathing needs further investigation, but we believe that these two behaviors with identical cycle lengths (CAM and periodic breathing) are likely to be related and determined by similar factors, whatever they might be. These data support the hypothesis that oscillatory activity responsible for periodic breathing is ubiquitously present but masked during apparently regular breathing by the normal regular stimulation from respiratory motoneurons. Periodic breathing occurs when this normal regular drive is decreased [i.e., in infants when core body temperature is raised (18)]. The adoption of one particular physiological state, i.e., regular tonic respiration with CAM or periodic breathing is likely to be dependent on the environmental conditions and maturity of respiratory control as well as the presence of any pathological conditions.
In conclusion, this study has addressed the limitations of previous studies that have examined whether respiration is chaotic. We investigated children in quiet sleep when breathing appears most regular. Correlation dimension estimates are consistent with a chaotic system. Furthermore, unlike in most previous studies, we used surrogate data analyses to test whether the apparently chaotic behavior was due to linearly filtered noise. We found this unlikely and concluded that the simplest system consistent with our data is a noisedriven, nonlinear, radialbasis model. Our data are the most convincing evidence that respiratory variability in infants is deterministic and not random, due to noise. A recent study has demonstrated reduced variability of respiratory movements in infants, who subsequently died of sudden infant death syndrome (SIDS) [Schechtman et al. (43)]. This observation was retrospective but suggests that, because the variability that we have observed during quiet breathing is deterministic, then further study using DST could allow early identification of infants at risk of SIDS from simple measurements of respiratory patterns.
We speculate that since CAM is an important contributor to the complexity observed during quiet breathing, further studies might demonstrate distinct patterns of CAM in infants with respiratory control problems; for example, absence of CAM might explain the reduction in variability observed by Schechtman et al. (43) in infants who died of SIDS.
We believe that the techniques employed by us in this study could be used in the future for studying the complex dynamic behavior of respiration and its control in health and disease.
Acknowledgments
For part of this project, M. Small was funded by the National SIDS Council of Australia.
Footnotes

Address for reprint requests: M. Small, Dept of Physics, HeriotWatt University, Riccarton, Edinburgh, EH14 4AS, UK. (Email:watchman{at}maths.uwa.edu.au).

Received 16 October 1997; accepted in final form 16 September 1998.

↵1 To within some arbitrary (possibly the machine) precision.
 Copyright © 1999 the American Physiological Society
DIMENSION ESTIMATION
In this appendix, we discuss the mathematical aspects of dimension estimation. Next, we present the new algorithm. Dimension estimation is logically, and operationally, divided into two steps: embedding the observed time series in a mathematical space and estimation of dimension from the embedded structure of the data. We begin with a brief description of the embedding process; see references for more details.
Embedding
Let x
_{1}, x
_{2},x
_{3}, … , x
_{N}represent a time series of N measurements of a system recorded with a fixed sampling rate. Construct an ndimensional vector time series v
_{t}, the timedelay embedding of the scalar time series, by
The Takens embedding theorem (27, 56) guarantees that the structure of the cloud of points (see Fig. 7) defined by a timedelay embedding is diffeomorphic (equivalent) to that of the measured system in most circumstances—actually Takens’ theorem says much more than this. The correlation dimension is invariant under diffeomorphism (smooth change of coordinates) and, hence, it is possible to analyze embedded data to discover a system’s fractal dimension.
An embedding depends on two parameters, the lag τ and the embedding dimension n. For an embedding to be suitable for successful estimation of dimension, one must choose suitable values of these parameters.
Embedding dimension.
A fundamental topological result (56) guarantees that anndimensional manifold can be successfully embedded in at mostR ^{2n + 1}, but more recent results (8) show that estimation of dimension requires that an embedding dimension need only be greater or equal to n. In practice, one could guess a suitable value for n by successively embedding in higher dimensions and looking for consistency of results; this is the method that is generally employed. However, other methods, such as the false nearest neighbor technique (10, 57), are now available to suggest the value of n. For the data used in our studies, the false nearest neighbor technique suggested that an embedding dimension of three or four would be sufficient. Using this knowledge, we repeated the dimension estimation procedure described below for increasing values of n, starting at n = 2.
Embedding lag.
There are two main methods (33) for choosing an appropriate value of the lag τ: the first zero of the autocorrelation function (4, 5) and the first minimum of the mutual information (1, 12, 23). The rationale of both of them, however, is to choose the lag so that the coordinate components of v _{t} are reasonably uncorrelated while still being “close” to one another. When the data exhibit strong periodicity, as is the case with respiratory patterns, a value of τ that is onequarter of the length of the average breath generally gives a good embedding. This lag is approximately the same as the time of the first zero of the autocorrelation function. Coordinates produced by this method are within a few breaths of each other (even in relatively highdimensional embeddings) while being spread out as much as possible over a single breath. Moreover, for embedding in three or four dimensions (as suggested by false nearest neighbor techniques), the data are spread out over onehalf to threequarters of a breath. This means that the coordinates of a single point in the three or fourdimensional vector time series v _{t} represent most of the information for an entire breath. This choice of lag is extremely easy to calculate and, for the data sets that we consider, it also seems to give much more reliable results than the mutual information criterion.
Dimension
Once we have embedded the data properly, we wish to measure the complexity of the “cloud” of pointsv _{t}. The measure we use in this paper is the correlation dimension.
To define the correlation dimension in a meaningful way, we generalize the concept of integer dimension to fractal objects with noninteger dimension. In dimensions of one, two, three, or more, it is easily established, and intuitively obvious, that a measure of volume V (e.g., length, area, volume, and hypervolume) varies as
The method most often employed to estimate the correlation dimension is the GrassbergerProcaccia algorithm. In this method, one calculates the correlation function and plots log C _{N}(ε) against log ε. The gradient of this graph in the limit as ε → 0 should approach the correlation dimension. Unfortunately, when a finite amount of data is used, the graph will jump about irregularly for small values of ε. To avoid this, one looks instead at the behavior of this graph for moderately small ε. A typical correlation integral plot will contain a “scaling region” over which the slope of log C _{N}(ε) remains relatively constant. A common way to examine the slope in the scaling region is to numerically differentiate the plot of log ε against logC _{N}(ε). This ought to produce a function that is constant over the scaling region, and its value on this region should be the correlation dimension (see Fig.8).
Unfortunately, as Judd (19) points out, there are several problems with this procedure. The most obvious of these is that the choice of the scaling region is entirely subjective (Fig. 8). For many data sets, a slight change in the region used can lead to substantially different results. Judd demonstrates that for many objects, including many fractals, a better description of C
_{N}(ε) is that for ε less than some ε_{0}
The GrassbergerProcaccia method assumes thatC _{N}(ε) ∝ ε^{dc}, but this new method allows for the presence of a further polynomial term that takes into account variations of the slope within and outside of a scaling region. This new method dispenses with the need for a scaling region and substitutes a singlescale parameter ε_{0}, which has an interesting benefit. For many natural objects, the dimension is not the same at all length scales. If one observes a large river stone, its surface at its largest length scale is very nearly two dimensional, but at smaller length scales, one can discern the details of grains that add to the complexity and increase the dimension. Consequently, it is natural to consider dimension d _{c} as a function of ε_{0} and writed _{c}(ε_{0}). For an alternative treatment of this algorithm, see, for example, Ref. 17.
By allowing our dimension to be a function of scale, we produce estimates that are both more accurate and more informative. We avoid some of the approximation necessary to define correlation dimension as a single number and we can extract more detailed information about the changes in dimension with scale. Finally, the data requirements of this algorithm are modest enough to be satisfied by our data.
SURROGATEDATA ANALYSIS
This appendix briefly outlines the surrogatedata methods utilized in this article. A more comprehensive description of these methods may be found in Ref. 59, and some caveats and recent extensions to these methods are outlined in Ref. 60. In Ref. 53, Small and Judd describe the application of correlation dimension as a test statistic for linear and nonlinear hypothesis testing. First, we discuss our choice of test statistic and then the rationale of linear and nonlinear surrogatedata methods.
Test Statistics
To compare the data to their surrogates, a suitable test statistic must be selected. A useful statistic must measure a nontrivial invariant of a dynamic system that is independent of the way surrogates are generated. Because we are interested in correlation dimension, which is independent of most surrogategeneration methods, it is natural to use this as our test statistic. Correlation dimension, as we have defined it, is a function of ε_{0}. There are several obvious ways to compare these curves. On many occasions, it is sufficient to compare the value of dimension for some fixed values of ε_{0}, and this is the method we used. Other possibilities include the mean value of the dimension estimate or the slope of the line of best fit. More sophisticated methods include statistical tests such as the χ^{2} test or the KolmogorovSmirnov statistic applied to the distribution of interpoint distances to determine whether the distributions are the same.
Traditional Surrogate Data
Different types of surrogate data are generated to test membership of specific dynamic system classes, referred to as null hypotheses. The three types of surrogates described by Theiler et al. (59), referred to as algorithms 0, 1, and 2, address the three null hypotheses: (0) i.i.d. (independently and identically distributed) noise; (1) linearly filtered noise; and (2) monotonic nonlinear transformation of linearly filtered noise. Each of these null hypotheses should be rejected for data generated by a nonlinear (nearly) periodic system, which is what we believe the respiratory system to be. Hence, we expect that each of these null hypotheses would be rejected for our respiratory data.
In a recent paper, Theiler (58) [and also Theiler and Rapp (61)] addresses this problem and proposes that a logical choice of surrogate for strongly periodic data, such as respiratory traces or electrocardiogram signals, should also be periodic. To achieve this, he decomposes the signal into cycles and shuffles the individual cycles. Our data, however, do not have a convenient point at which to break the cycles. To break the cycles at peak inspiration (or expiration) would introduce a degree of subjectivity, since the point at which to separate the cycles is not clearly defined, and it would introduce points of nondifferentiability (that is, sharp points, cusps, or corners) that are not present in the original data.
Theiler’s (58) alternative null hypothesis for strongly periodic signals is not a characterization of great interest to our present discussion. Theiler proposes that surrogates generated by shuffling the cycles address the hypothesis that there is no dynamical correlation between cycles. The null hypothesis that the data are generated by a noisedriven periodic orbit is of far greater significance to our analysis. Cycletocycle variation in breath structure is a problem distinct from our present inquiry but also of great interest (55). We have considered the application of cycleshuffled surrogates to similar data elsewhere (52).
NoiseDriven Nonlinear System Surrogates
Hypothesis testing with surrogate data is, essentially, a modeling process. To test whether the data are consistent with a particular hypothesis, one first builds a model that is consistent with that hypothesis and has the same properties as the original data, then one generates surrogate data from the model and checks that the original data are typical under the hypothesis by comparing them with the surrogate data. For surrogates generated by algorithm 0, 1, or 2, the model used is linear. Each of these surrogate tests addresses a null hypothesis that the data are either linear or some (linear or monotonic nonlinear) transformation of a linear process.
To address the null hypothesis that the data come from a noisedriven nonlinear system, which, one supposes, will have a focus or stable periodic orbit, we built a nonlinear model and generated noisedriven simulations from that model. These simulations are our surrogate data. The nonlinear model that we built from the data is a radialbasis model based on the methods of Judd and Mees (21, 22), and the test statistic we used was correlation dimension. Radialbasis models are used because they are known to be effective in modeling a variety of nonlinear dynamic systems, and the authors have at their disposal a sophisticated software implementation of this modeling method. To be precise, the null hypothesis we tested is that the data are consistent with a nonlinear system that can be described by a radialbasis model and that the data of such a system can be modeled adequately by using the algorithms we used. Rejection of the null hypothesis could imply that the data cannot be described by a radialbasis model or that the modeling algorithm failed to build an accurate model.
Building a nonlinear model of data is a decidedly nontrivial process. All that really concerns us here is that we were somehow able to generate a radialbasis model of the data. Figure 4 gives an example of some data generated by the methods we used. A more precise definition of the model is given in appendix .
The conclusions that can be drawn from testing with these nonlinear models are several. First, the surrogate data can indicate that our data are not consistent with a nonlinear system of the type generated by our modeling procedure. Second, it will determine whether correlation dimension can distinguish between data and surrogates. Finally, and most importantly, this is a test of the modeling procedure itself. If the null hypothesis cannot be rejected on the basis of our analysis, then this will indicate that the model we have built is an accurate model of the data, with respect to correlation dimension. Even if correlation dimension cannot distinguish between data and surrogate, other measures, for example, largest Lyapunov exponent, may.
There is one important caveat. Our methods do not test for the presence of a general nonlinear periodic orbit. They only test for the presence of a nonlinear periodic orbit that can be accurately modeled as the sum of radialbasis functions of the form described in appendix. This is not particularly restrictive, since experience has shown that such functions can model a wide range of phenomena.
NONLINEAR MODELS
This appendix provides a brief description of the form of the nonlinear model utilized in our analysis. The models used are a slight generalization of the methods described by Judd and Mees (21).
From a scalar time series y
_{t}, we build a model of the form
To build these models, we apply the methods described in Refs. 21 and22 and, more recently, by Small and Judd (51). Of the modifications discussed in Ref 51, we use the broader class of basis functions (as shown in Eq. 4 ), the directed basis selection, and the improved description length calculation.
The selection of the parameters a _{0},b _{i}, ς_{j}, ν_{j}, λ_{j}, and μ_{j} and the projectionsP _{j} is a very difficult nonlinear optimization problem. For each value of n, the parameters ς_{j}, ν_{j}, μ_{j}, and P _{j} are selected randomly (we guess a large number of possibilities and select the best). The a _{0}, b _{i}, and λ_{j} are then selected by using least sum of squares optimization. The optimal value of n is chosen by minimizing the description length (35) of the resulting model. The functions f and g are fitted separately but iteratively (first fit f to the data, then fit g to the prediction errors, refit f and then refit g, repeat several times).
RARM
In this appendix, we briefly describe RARM, their important features in the context of this study, and an algorithm for constructing them. We begin this appendix with an overview of the important ideas, and in the later sections we provide a more mathematical treatment.
The distinguishing feature of a linear autoregressive model of a time series is the assumption that at any time the next observation of the system is always a specific weighted average of selected past observations of the system, plus some quantity of noise. For example, the tidal volume of the next breath of an infant is always a specific weighted average of the tidal volume of the last seven breaths plus some unknown noise component. Determining the most appropriate weighted average from an observed time series is a classic problem of applied mathematics, statistics, and signal processing, for which there are wellestablished algorithms.
Whether one should have a model that depends on the last 7 observations, or the last 3, or the last 20 is a problem that has received a lot of attention recently, and the attention has culminated in the MDL principle. The essential problem is that to obtain a useful phenomenological model one should avoid overfitting observed data. A practical procedure is to build the best weightedsum model by using one previous observation, then two previous observations, then three, and so on, and stop building models when looking back further provides no significant improvement in predictions that the models make; that is, one should trade off improvements in modeling errors against the number of previous observations used. Akaike (3) was the first to propose a concrete algorithm to do this, and Schwarz (44) improved on this theme with deeper theoretical considerations. In the last decade, there have been significant advancements in understanding through the application of information theory by Rissanen (35) and Wallace and Freeman (66): these new approaches are often termed “minimum description length” or “minimum message length” models. The new algorithms have shifted focus from the number of previous observations used to (in essence) the total number of significant figures needed to specify the weights in the weighted sum.
RARMs are obtained when one strictly applies MDL principles to linear autoregressive models. Typically, in a weightedsum model, some previous observations have smaller weights than others and are, hence, less significant and may be eliminated. An MDL algorithm may result in a model comprising just a weighted sum of, say, the previous breath, the breath previous to that, and the breath seven breaths previous. An early suggestion of the idea of a RARM was by Haggan and Oyetunji (16) in 1984. The first suggestion of data compression as a modelselection criterion was made by Wallace and Boulton in 1968 (65). The ideas have been rediscovered since, but have never taken off because MDL algorithms were not available. Recently, Judd and Mees (21) have described a general MDL algorithm from which RARMs follow naturally.
The key feature of RARMs that makes them important to the present study is that they provide an objective means of identifying periodicities in time series. It is clear how this is possible: if a RARM selects that the current breath is best predicted by a weighted average of the previous breath, the breath previous to that, and the breath seven breaths previous, then there is evidence of a sevenbreath period. Extensive studies by some of the authors (54) have demonstrated that whenever periodicities are apparent when classic techniques are used, such as spectral analysis or autocorrelation, then the same periods are selected in a RARM. Furthermore, when known periodicities are present (in artificial data), RARM can detect periodicities when the classic techniques do not. On the other hand, extensive surrogatedata analysis has shown that RARMs do not falsely identify periodicities that are not present. In summary, RARMs perform as the wellestablished MDL theory predicts and provide a useful, reliable tool for detecting periodicities in breathing (54).
To describe the RARM algorithm, we will first discuss linear modeling. Afterward, we will describe the description length criteria of Rissanen (35) and our implementation of a modelselection algorithm.
Autoregressive Models
The traditional autoregressive model of order n [an AR(n) model] attempts to model a time series
However, a time series exhibiting periodic behavior with period τ would have strong dependence of x _{t} onx _{t − τ}. Hence, by building an AR(n) model and determining which parameters are most significant, it may be possible to estimate the period of some periodic behavior or, more significantly, several different periods within the same series.
Deciding which parameters are most significant requires sophisticated methods. To do so just on the basis of the size of the coefficientsa _{1}, a _{2},a _{3}, … , a _{n} will rarely be useful. We discuss the selection problem in Description Length below.
What we are attempting to do is to fit the best model to the data. A traditional AR(n) model has n parameters, but it may be the case that only some of these are necessary. Essentially then, we are looking to find the best model of the form
Using the concept of description length, we have a method of deciding which parameters offer a substantial improvement to the model. Rissanen’s description length (35) is just one way to compare the size of a model to its accuracy, other methods include the Schwarz (44) and Akaike (3) information criteria.
Description Length
Roughly speaking, the description length of a particular model of a time series is proportional to the number of bytes of information required to reconstruct the original time series.1 That is, the compression of the data gained by describing the model parameters (a
_{0}, a
_{1}, a
_{2}, . . . a
_{k}, l
_{1}, l
_{2}, … , l
_{k}, k) and the modeling prediction error (
Obviously, if the time series does not suit the class of models being considered, then the most economical way to do this would be to simply transmit the data. If, however, there is a model that fits the data well, then it is better to describe the model to the receiver in addition to the (minor) deviations of the time series from that predicted by the model (see Fig. 9). Thus description length offers a way to tell which model is most effective. Our encoding of description length is identical to that outlined by Judd and Mees (21) and follows the ideas described by Rissanen (35). For a model of the form shown in Eq. 6, the description length is bounded by (see Ref. 21)
Model Selection
With an appropriate measure of description length, we are able to measure the appropriateness of a given model. To determine which model is best, we apply the modelselection algorithm of Judd and Mees (21) to the trivial case, i.e., the case in which only linear models are required.
Let V
_{i} = [x
_{n+1−i}, … , x
_{N−i}] for i = 1, 2, … , n and define V
_{0} = [1, 1, … , 1]. To select the best RARM, one must select the set B ⊆ {0, 1, 2, … , M} such that the model
Our modeling method implicitly requires a parameter M, the maximum number of breaths in the past to examine. To overcome this, we examine the models produced for a range of different maximum model sizes M and ensure that we include values of M that exceed the period of any periodicities we expect to observe. The RARM procedure will produce a (possibly changing) indication of period as a function of M. We then look for the stage at which the previous breaths used to predict the next do not change by increasing the maximum model size. From this we deduce the period of any periodic behavior.
For example, by applying this analysis to one of our tidal volume data sets we observe the following. For M = 2, 3, 4, the MDL model uses information in the last two breaths to predict this breath. ForM = 5, 6, 7, 8, the model is extended to include the breath 5 breaths ago. For M = 9, 10, 11, the model also includes the breath 9 breaths ago. Finally, for M = 12, 13, 14, … , 60, the model predicts the magnitude of a breath using the previous breath and the one before that, the breath 5 breaths ago, the breath 9 breaths ago, and the breath 12 breaths previous. We conclude that for M ≥ 12 the model is