Abstract
CurranEverett, Douglas, Yiming Zhang, M. Douglas Jones, Jr., and Richard H. Jones. An improved statistical methodology to estimate and analyze impedances and transfer functions. J. Appl. Physiol. 83: 2146–2157, 1997.—Estimating the mathematical relationship between pulsatile time series (e.g., pressure and flow) is an effective technique for studying dynamic systems. The frequencydomain relationship between time series, often calculated as an impedance (pressure/flow), is known more generally as a frequencyresponse or transfer function (output/input). Current statistical methods for transfer function analysis 1) assume erroneously that repeated observations on a subject are independent,2) have limited statistical value and power, or 3) are restricted to use in single subjects rather than in an entire sample. This paper develops a regression model for transfer function analysis that corrects each of these deficiencies. Spectral densities of the input and output time series and the crossspectral density between them are first estimated from discrete Fourier transforms and then used to obtain regression estimates of the transfer function. Statistical comparisons of the transfer function estimates use a test statistic that is distributed as χ^{2}. Confidence intervals for amplitude and phase can also be calculated. By correctly modeling repeated observations on each subject, this improved statistical approach to transfer function estimation and analysis permits the simultaneous analysis of data from all subjects in a sample, improves the power of the transfer function model, and has broad relevance to the study of dynamic physiological systems.
 discrete Fourier transform
 frequencydomain regression
 frequencyresponse function
 mixedeffects model
 spectral analysis
estimation of the mathematical relationship between pulsatile time series is a timehonored approach to the study of dynamic systems. Using this strategy, physiologists have explored respiratory mechanics (11, 15, 23, 30, 39), pulmonary ventilation (29,33), cardiovascular function (3, 7, 24, 27, 32), and cardiorespiratory regulation (4, 25, 34, 35). The relationship between time series is calculated most often in the frequency domain (5, 24): this entails first transforming the time series, originally functions of time, into their equivalent Fourier coefficients, written as functions of frequency. When two time series represent the input and output of a dynamic system, the mathematical relationship between them is often designated the transfer function (h), defined in terms of frequency ( f ) as h( f ) = output( f )/input( f ), where the output and input functions are the Fourier transforms of the original time series. When input is flow and output is pressure, the transfer function is called impedance. A transfer function is a complex expression, with real and imaginary parts, that can be described by amplitudes (output amplitude divided by input amplitude) and phase angles (timing of output with respect to input).
In experimental situations, transfer function analysis is complicated by betweensubject variability as well as by random measurement error. Current approaches to transfer function analysis 1) assume erroneously that repeated observations on a subject are independent,2) have limited statistical value and power, or 3) because they are unable to account for betweensubject variability, are restricted to use in single subjects rather than in an entire sample. This paper develops a regression model for transfer function analysis that corrects each of these deficiencies and has broad relevance to the analysis of dynamic physiological systems. Before we derive and illustrate this improved methodology, we review previous approaches to applied transfer function analysis.
APPLIED TRANSFER FUNCTION ANALYSIS
Description.
Many studies (3, 7, 11, 15, 23, 26, 28, 34, 42) simply describe the effect of an experimental intervention on a transfer function: they report, for example, changes in the frequencies at which maxima in impedance amplitude occur. But to draw conclusions about an underlying population, one must use more than rudimentary description: one must employ the inferential techniques of confidence intervals and significance tests (38).
Usual indexes of betweensubject variability.
In most physiological studies (see Ref. 24), betweensubject variability is estimated, first by deriving the amplitude and phase angle of a transfer function for all subjects and then by calculating SDs. This standard approach of handling betweensubject variability fails to account for the fact that repeated observations on a subject (e.g., observations made during baseline and then during an experimental intervention) are correlated (21). Because of this correlation, the true error variability is underestimated, and the reported values for the SDs of amplitude and phase underestimate the true variabilities (21). The shows how correlation decreases variability.
Separate analyses of amplitude and phase.
Some researchers (4, 35, 40) analyze the amplitude and phase of a transfer function as if they were unrelated variables. But the real and imaginary parts of a transfer function, from which amplitude and phase are derived, are determined simultaneously. Therefore, the components of a transfer function, either its real and imaginary parts or its amplitude and phase, must also be analyzed simultaneously. There is a quantitative benefit: the simultaneous analysis of jointly derived data improves statistical power (10).
Derivation and analysis of analog parameters.
Because transfer function estimation simply reduces a relationship between pulsatile time series to its mathematical form, the correspondence of a transfer function to physical characteristics of a system is unclear. To circumvent this limitation, some investigators (15, 39, 42) first obtain an estimate of the transfer function and then, from that estimate, derive analog parameters that represent general characteristics (e.g., total compliance) or specific components (e.g., lung tissue elastance) of the system. It is the estimates of the analog parameters, rather than the estimate of the transfer function, that are subjected to statistical analysis. Although this approach gives physiological meaning to a transfer function, it does have a statistical drawback: the estimates of the analog parameters are correlated (38). This means that the analysis of only one (preselected) parameter is valid: the statistical outcomes of the remaining analyses will be related to the outcome of the initial analysis.
Timedomain analysis.
The timedomain technique of bivariate autoregression has been used also to estimate a transfer function (22). Bivariate autoregressions (20) can be fit to the original time series by using the YuleWalker equations, which define the relationships between the covariance functions, the crosscovariance function, and the autoregression matrices (18, 41). The resulting bivariate autoregressive estimates are related to the autocorrelation and crosscorrelation functions, often used in physiology to estimate a transfer function (24, 25, 32). From the estimated autoregressive matrices, estimates of the spectral densities of an input and output and the crossspectral density between them can be calculated; it is from these autoregressive spectra that a transfer function can be computed (see preliminary data analysis). The advantage of autoregression is that it can produce smooth spectral estimates of a transfer function. The disadvantage is that the autoregressive estimates are correlated across frequency (20,37): therefore, the statistical properties required to compare transfer functions at specific frequencies are unknown.
Regression techniques.
Since the 1940s, regression analysis in the frequency domain has been used in transfer function estimation (see Ref. 6). In this approach, the real and imaginary parts of the Fourier transforms of two time series are used to obtain a regression estimate of the transfer function (6, 16, 31, 37). The classic regression model (see Eq.3 ) assumes, however, that estimates of the transfer function during a given condition are virtually identical for all subjects. This is unrealistic: physiological transfer function estimates vary considerably among subjects (7, 8, 28, 32, 34, 42). Indeed, it was explicitly because of betweensubject differences that O’Rourke and Taylor (Ref. 28, p. 368–369) presented only individual impedance spectra: “Although statistical analysis of a large number of experimental results is desirable, it is obvious that features of the impedance curves can be lost if data are pooled from [many subjects] …”
Because it does not allow for betweensubject variability, the classic regression model can estimate a transfer function in single subjects only (37); this confines statistical and biological inferences to single individuals rather than to an underlying population. The methodology we develop here extends this classic regression approach in such a way that all subjects in a sample can be analyzed simultaneously. This statistical enhancement has important theoretical and applied advantages. In the next section, we review the experimental problem that motivated the development of this methodology.
ORIGINAL EXPERIMENTAL PROBLEM
The development of this methodology was driven by our interest in the effects of systemic circulatory perturbations on dynamic properties of the cerebral circulation (see Ref. 8). In this context, we considered systemic arterial blood pressure to be the input and cerebral cortical blood flow to be the output of a linear vascular system. In 10 rabbits, arterial blood pressure, measured by highfidelity transducer catheter, and cortical blood flow, measured by laserDoppler flowmeter, were sampled for 1 min during two conditions: baseline and combined hypoxiahypercapnia. Our experimental objective was to estimate the transfer function between arterial blood pressure and cortical blood flow during each condition by using data from the entire sample. To do this, we selected a 20s segment from each time series, giving 800 pairs of observations (sampling interval 0.025 s) for each subject during each condition (Fig. 1). The individual transfer functions derived from each of these bivariate time series could be estimated reliably at the frequencies of heart rate and its first harmonic.1 For the frequency of heart rate, Fig. 2 depicts the betweensubject variability inherent to transfer function estimation. The methodology that follows was devised to incorporate betweensubject variability into the regression model for the transfer function between arterial blood pressure and cortical blood flow.
PRELIMINARY DATA ANALYSIS
Before we derived the regression model for the transfer function, however, we verified an important statistical assumption: that the error term of the model, i.e., the spectral density of the residual series, was nearly constant around the frequencies of heart rate and its first harmonic. We did this using bivariate autoregression (20), identifying the order of the model by Akaike’s information criterion (1, 2). The estimated autoregressive spectra of pressure and flow (ŝ
_{xx} andŝ
_{yy}, respectively) showed striking peaks at the frequencies of heart rate and its harmonics (Fig.3). The estimated spectral density of the residual series was then derived in a twostep procedure analogous to that used in regression analysis. First, the transfer function estimateĥ was calculated from the bivariate autoregressive spectral estimates as
Next, we review time series regression in the frequency domain. As we develop the full regression model for the transfer function, we illustrate the process for the transfer function (evaluated at heart rate) between arterial blood pressure and cortical blood flow using data from Ref. 8.
TIME SERIES REGRESSION
Hannan (16) pioneered time series regression in the frequency domain (see also Refs. 6, 31, and 37). Ifx
_{j} and y
_{j}represent the input and output time series, respectively, each with n observations at time intervals of δ, then the discrete Fourier transforms (z
_{ν}) of these series are
The example.
Because we study dynamic vascular properties at discrete frequencies, the first task is to identify the frequency of heart rate; we do this by using the spectral density of arterial pressure. Then, for each of 15 0.05Hz frequency bands2centered at heart rate, we calculate the real and imaginary components (Eq. 2 ) of the discrete Fourier transforms of both the input and output time series (Table 1).
Transfer function estimation.
Using the real and imaginary parts (Eq. 2
) of the Fourier transforms, the transfer function h can be approximated in the frequency domain as
where * denotes the complex conjugate. These estimates are smoothed periodogram estimates with uniform weights, a smoothing span of 2w + 1, and 4w + 2 degrees of freedom (19).3 Next, the real and imaginary parts of the crossspectral density, the cospectrumĉ
_{yx} and the quadrature spectrumq̂
_{yx}, are estimated as
The example.
From the real and imaginary components of the discrete Fourier transforms within the 15interval frequency band (see Table 1), we next obtained a single smoothed estimate of the transfer function (Table2; see Transfer function estimationabove). The phase estimate is meaningful only if there is high coherence between the input and output: that is, if changes in the output correspond closely to changes in the input. Based on the estimated autoregressive spectrum of squared coherence (see Fig. 3), we expect to obtain a reliable estimate of the transfer function at the frequency of heart rate. Although smoothing increases the precision of the transfer function estimate, it can also decrease squared coherence. Therefore, after we get the single transfer function estimate, we verify that squared coherence across the 15interval frequency band remains high
RANDOM SUBJECT EFFECT
In the preceding section, we reviewed the classic regression model (Eq. 3 ) that assumes, for some conditions, that estimates of a transfer function are virtually identical for all subjects. In this section, we first add a random parameter γ to the transfer function model, thereby transforming the classic fixedeffects model (Eq.3 ) into a mixedeffects model that can handle inherent betweensubject variability in the transfer function estimates. Next, we use this enhanced statistical model to estimate the transfer function during b experimental conditions for an entire sample of N subjects. Finally, we use a likelihood ratio test to evaluate the improvement gained by adding the random subject effect to the transfer function model.
Development of the mixedeffects model.
When a transfer function is estimated in a group of subjects, the estimates will differ because of random betweensubject variation. This variation can be included in the transfer function model, thereby improving its power, by using a random subject effect (see Ref. 21). With this revision, the classic regression model (Eq. 3
) becomes, for subject k
When the transfer function is estimated over a band of mfrequencies, the Fourier transforms
For a sample of N subjects, the estimate ĥ of the transfer functions is
Analysis of the random effect.
Statistical evaluation of the random subject effect is based on a likelihood ratio test, an essential component of which is the maximum likelihood estimate of ς^{2}. [A maximum likelihood estimate is the value of a population parameter that most likely would have produced the sample observations (see
).] The maximum likelihood estimate of ς^{2} is obtained from a multivariate complex Gaussian distribution (12)
By substituting Eq. 11
into Eq. 10, the likelihood can now be written as a function of the single variablec
^{2}, the ratio of the error variances
The example.
Before we compare the b estimates of the transfer function, we test whether the random effect for betweensubject variability,
TRANSFER FUNCTION ANALYSIS
In this section, we first review statistical contrasts and contrast coefficients; contrast coefficients are part of the test statistic used to compare transfer function estimates. Next, we detail the test statistic itself. Last, we calculate 100(1 − α)% confidence intervals for the amplitude and phase of the transfer function.
Contrasts.
Comparisons of the transfer function estimates require contrast coefficients, which are used to compute statistical contrasts: a contrast represents the size of a comparison among treatment means (see Ref. 38). The sample contrast D̂estimates the population contrast D and is any linear combination
The example.
The transfer function estimates, obtained by using the mixedeffects model, are
Test statistic.
A statistical comparison between transfer function estimates requires the test statistic
The example.
The enhanced ability of the mixedeffects model to detect transfer function differences becomes clear with the example presented in Table 3. Using each model, we first estimate the two transfer functions and then compute the same contrast:ĥ
_{1} (hypoxiahypercapnia) vs.ĥ
_{0} (baseline). To compareĥ
_{1} toĥ
_{0}, C, the row vector of contrast coefficients in Eq. 14, must be
The explanation for the improved power of the mixedeffects model can be appreciated most easily by comparing the fixedeffect and mixedeffects models for subject k. When b transfer functions are estimated for subject k, the fixedeffect model is
Confidence intervals for amplitude and phase.
Confidence intervals for the amplitude and phase of the transfer function—or, by analogy, a transfer function contrast—can be calculated by using procedures derived by Groves and Hannan (14). The approximate 100(1 − α)% confidence interval for the amplitudeA is
The example.
Last, we calculate 95% confidence intervals for the amplitude and phase of the transfer function during baseline by using the inequalities in Eqs. 15 and 16 (Table4). In addition to the estimates of amplitude and phase, Â and θ̂, respectively, and the estimated SD of the transfer functionς̂_{h}, we use the following information: α = 0.05, m = 15, N = 10, z _{α/2} = 1.96, and t _{α} = 1.65 (298 degrees of freedom).
DISCUSSION
Previous approaches to transfer function estimation and analysis1) assume erroneously that repeated observations on a subject are independent, 2) have limited statistical value and power, or 3) are restricted to use in single subjects rather than in an entire sample. The methodology we develop here is founded on a regression approach to time series analysis and corrects each of these deficiencies. By virtue of the random effect for betweensubject variability in the transfer function estimates, our approach 1) can correctly model repeated observations on a subject, 2) has general relevance and increased statistical power, and 3) enables the simultaneous analysis of data from all subjects in a sample. Before discussing advantages and applications of this methodology, we review the use of a random subject effect in another statistical model developed for time series analysis.
Another random subject effect.
Diggle and AlWasel (9) incorporate a random effect for betweensubject variability in their model for the spectral density of an output time series measured during multiple experimental conditions. In our approach, it is the model for the transfer function between the input and output time series that has the random betweensubject effect. When the transfer function is of primary experimental interest, our methodology is preferable.
Advantages of this new transfer function model.
Our addition of the random effect γ_{k} to the classic transfer function model has important theoretical and applied advantages. Three theoretical advantages pertain to the model itself. First, the new transfer function model (Eq. 6 ) accounts explicitly for biological differences that exist between subjects. Second, the model handles correctly the correlation among repeated observations on each subject. Third, the model is better quantitatively (see Table 3).
Three applied advantages relate to the use of this new transfer function model for the study of dynamic systems. First, the model provides for the simultaneous analysis of multiple subjects; the ability to properly analyze an entire sample is essential to draw conclusions about an underlying population. Second, the model provides for more accurate estimates of the mean and SD of the transfer function; the impact on the estimate of the SD,ς̂_{h}, is especially pronounced (see Table 4). Third, the model provides for an improved ability to detect differences between transfer functions (see Table 3).
Application to physiology.
This improved statistical approach to transfer function estimation and analysis has broad relevance to the study of dynamic physiological systems. It can be used to estimate the mathematical relationship between virtually any pulsatile time series. In our animal model (Ref.8), we estimate and analyze a transfer function at discrete frequencies: spontaneous heart rate and its first harmonic. Many physiological investigations, however, use random pacing to generate nearly continuous transfer function spectra. The methodology we present here accommodates this common situation also: transfer function analyses at multiple frequencies will be uncorrelated as long as nonoverlapping frequency bands are used.
Summary.
The transfer function methodology developed here furthers regression techniques to analyze the mathematical relationship between pulsatile time series. By correctly modeling repeated observations on each subject, this improved methodology can account for inherent betweensubject variability in the transfer function estimates; this makes possible the simultaneous analysis of all subjects in a sample and improves the power of the transfer function model. This methodology can be used to estimate and analyze the dynamic properties of a variety of physiological systems.
Acknowledgments
We thank Dr. Mary Anne Farrell Epstein and the anonymous reviewers for their constructive comments and suggestions.
Footnotes

Address for reprint requests and software inquiries: D. CurranEverett, Dept. of Pediatrics, B195, Univ. of Colorado Health Sciences Center, 4200 East 9th Ave., Denver, CO 80262 (Email:dcurranevere{at}castle.cudenver.edu).

The development and analysis of the statistical methodology was funded by National Institute of General Medical Sciences Grant GM–38519 (to R. H. Jones). The physiological studies cited were funded by the Academic Enrichment Fund of the University of Colorado School of Medicine.

The authors are pleased to assist interested readers with the software that executes this methodology.

↵1 Although our statistical methodology can be used to estimate a transfer function at any frequency of interest, we use only the frequency of heart rate to illustrate the procedure.

↵2 The number of frequency bands is determined empirically. The choice involves a tradeoff between increased precision of the transfer function estimate (i.e., greater equivalent degrees of freedom associated with the error variance) and decreased resolution of the frequency interval within which the transfer function is estimated.

↵3 These spectral density estimates are proportional to their true spectral densities: to be a proper estimate, each would have to be multiplied by a constant so that the area under its spectral density curve estimated the variance of its respective time series.

↵4 A complex Gaussian distribution has random real and imaginary parts that are statistically independent and have identical variance ς^{2}/2.
 Copyright © 1997 the American Physiological Society
Appendix
This shows how correlation decreases variability, details the test statistic used to evaluate the improvement gained by adding the random subject effect to the transfer function model, and reviews maximum likelihood estimation and statistical contrasts and contrast coefficients. Maximum likelihood estimation is used to obtain the transfer function estimates (Eq. 7 ). Contrast coefficients are a component of the test statistic (Eq. 14 ) used to compare transfer function estimates.
Impact of Correlation on Variability
Repeated observations on a subject are correlated (21) by virtue of that subject’s particular biological makeup. When investigators fail to consider this correlation, as they usually do, they underestimate the true variability in the amplitude and phase estimates. To illustrate the impact of correlation on variability, imagine an investigation in a sample of n subjects. In each subject, some random variable B is measured during two experimental conditions: a control period (for which B is designatedB _{1}) and a subsequent intervention period (for which B is designated B _{2}). Last, supposeB _{1} and B _{2} are distributed normally: B _{1} with SD ς_{1}, andB _{2} with SD ς_{2}.
If the random variables B
_{1} andB
_{2} are considered jointly, then the distribution of the variable pair (B
_{1}, B
_{2}) can be thought of as a joint probability density function (37), the bivariate normal distribution. For this bivariate normal distribution, ς_{21}, the SD of the conditional distribution ofB
_{2} given that B
_{1} equals a specific value, depends on the correlation ρ betweenB
_{1} and B
_{2} (17)
Maximum Likelihood Estimation
Sample observations are used to estimate population parameters: for example, the sample mean
To illustrate maximum likelihood estimation, the approach used in this methodology, imagine a sample of n independent observations,y
_{1},y
_{2}, … , y
_{n}, drawn from a population that is distributed normally with mean μ and variance 1. The probability density function f describes the normal distribution of the population
The value of μ that maximizes L also maximizes the logarithm of L. This is convenient: exponential terms are characteristic of many likelihood functions, and because products become sums, it is often simpler to maximize the logarithm of a likelihood rather than the likelihood itself. The logarithm l of the likelihood in Eq.EA2
is
Test Statistic Used to Evaluate Random Subject Effect
To test whether the random subject effect improves the transfer function model, we evaluate the corresponding change in l, the −2 ln likelihood (Eq. 12
), by using a likelihood ratio test (see Analysis of the random effect). In this situation, the null hypothesis H
_{0} is that there is no betweensubject variability in the transfer function estimates. This null hypothesis is written formally as
The first possibility is that ĉ ^{2} = 0. In this case, the values of l for the fixedeffect (random effect absent) and mixedeffects (random effect present) models will be identical: therefore, Δl = 0, which can not be distributed as χ^{2}. If ĉ ^{2} = 0, then the statistical decision is fail to rejectH _{0} and conclude that the data fail to demonstrate betweensubject variability in the transfer function estimates.
The second possibility is that ĉ
^{2} > 0. In this case, H
_{0} : c
^{2} = 0 can be restated by moving it off the edge of the parameter space by an infinitesimal distance ε, where ε > 0:
Statistical Contrasts
To illustrate contrasts and their coefficients, imagine an investigation of some biological variable during four experimental conditions: two control periods, during which the values of the sample means are
The contrast used to make inferences about the population means μ_{1} and μ_{3} (objective 1) is
The contrast used to make inferences about the average effect of the intervention on the population mean of the biological variable (objective 2) is