|Home | About | Journals | Submit | Contact Us | Français|
Statistical methods for detecting changes in longitudinal time series of psychophysiological data are limited. ANOVA and mixed models are not designed to detect the existence, timing, or duration of unknown changes in such data. Change Point (CP) analysis was developed to detect distinct changes in time series data. Preliminary reports using CP analysis for fMRI data are promising. Here, we illustrate the application of CP analysis for detecting discrete changes in ambulatory, peripheral physiological data leading up to naturally occurring panic attacks (PAs). The CP method was successful in detecting cardio-respiratory changes that preceded the onset of reported PAs. Furthermore, the changes were unique to the pre-PA period, and were not detected in matched non-PA control periods. The efficacy of our CP method was further validated by detecting patterns of change that were consistent with prominent respiratory theories of panic positing a relation between aberrant respiration and panic etiology.
Longitudinal assessment of physiological data has long been a focus in basic and applied research. Such assessments can offer insight into both chronic and discrete events. Examples for the latter include the examination of recurrent syncopal episodes (Inamdar, Mehta, Juang, & Cohen, 2006), ventricular arrhythmias (Krahn, Klein, Skanes, & Yee, 2004), obstructive sleep apnea (Collop, 2008), and/or seizures (Chang, Ives, Schomer, & Drislane, 2002). Examples in the area of psychopathology can include the onset of a manic or psychotic episode, or the often dramatic and abrupt occurrence of a panic attack (PA).
The apparent “out-of-the-blue” nature of PAs is intriguing, both from a psychological and biological perspective. Consequently, attempts have been made to capture their physiological profile in laboratory and real-life settings. However studies to date remain sparse and results are mixed. While some studies report physiological arousal immediately prior to and during the reported PAs (Goetz et al. 1993; Taylor et al., 1986; Freedman, Ianni, Ettedgui, & Puthezhath, 1985), others were unable to detect such changes (Gaffney, Fenton, Lane, & Lake, 1988; Hibbert & Pilsbury, 1988). The examination of physiological precursors to naturally occurring PAs poses a unique challenge due to their sudden occurrence. Thus, extended recording times and sophisticated ambulatory monitoring that do not greatly interfere with patients’ activities are required. At the same time, as advancements in ambulatory multi-channel technology continue to improve our abilities to record longitudinal data, new challenges of processing and analyzing these data surface. Questions such as the following arise: “When do discrete events occur and how long do they last?”, “Are there detectable precursors of psychophysiological events?”, and “How do we know that potential events are not simply random noise?” Adding to these challenges is the fact that the obtained biosignals (e.g., heart rate [HR] and respiration rate [RR]) are impacted by confounding variables, such as movement and speech, whose influence must be controlled.
The difficulty in addressing such questions arises, in part, from the apparent lack of appropriate statistical tools to provide meaningful answers to these questions. Most inferential statistics are designed to test specific hypotheses concerning the nature of the data and are not designed to answer questions related to when, and if, changes occur in a time series of data.
In the following we describe commonly used statistical methods for the analysis of physiological data. We then illustrate the application of a statistical approach for analyzing longitudinal data that until recently has rarely been used to examine psychophysiological measures.
Psychophysiologists have typically analyzed longitudinal data from a group (or groups) of individuals using repeated measures ANOVA (Keselman, 1998; Vasey & Thayer, 1987; for a specific example, see Goetz et al., 1993). In the context of longitudinal data, this approach, however, suffers from several disadvantages: First, repeated measures ANOVA requires complete data for each participant. Secondly, traditional repeated measures ANOVA does not allow for time varying covariates. Controlling for movement at each point in the time series could help establish if a change (in cardiac activity, for example) is due to psychological arousal over and above the influence of physical activity. Thirdly, ANOVA is insensitive to discrete, time limited events that occur at some unknown point during a long time series. Finally, repeated measures ANOVA is subject to a number of very restrictive assumptions (e.g., sphericity) which are rarely met (Tabachnick & Fidell, 2007). Although various methods can be used to adjust for violations of these assumptions, such as Greenhouse-Geisser or multivariate approaches to repeated measures ANOVA, the power of these analyses is often reduced by these methods. In sum, ANOVA is not optimally designed to detect a “signal” that might exist in the midst of background noise; consequently the question of relevant change often remains unanswered.
A more general mixed model approach, such as hierarchical linear modeling (HLM), can effectively deal with some of the limitations of repeated measured ANOVA (Keselman, 1998; Tabachnick & Fidell, 2007). In particular, HLM can easily handle data that are missing at random (Hall et al., 2001; Hamer & Simpson, 2009) and can incorporate time varying covariates (Singer & Willett, 2003). However, mixed models, are designed to detect specific, hypothesized, overall trends in the time series beginning at specified time points. They often cannot detect changes of unknown timing and duration, especially if the change is time limited. Different analytical approaches are necessary to answer questions related to discrete, time delimited changes occurring in time series data (Lindquist, Waugh, & Wager, 2007; Robinson, Wager, & Lindquist, in press)
Change Point (CP) analysis (e.g., Basseville & Nikiforov, 1993; Lindquist, Waugh, & Wager, 2007) is an approach adapted from statistical control theory that was developed to detect significant changes in time series data. The CP method also allows the examination of various kinds of hypotheses through the imposition of different constraints on the CP parameters. There now exist numerous variations of CP analysis emanating from various fields, such as statistical control theory, econometrics (e.g., Fomby & Lin, 2006), and psychophysiology (Lindquist et al., 2007; Robinson et al., in press). These various analytic strategies are largely comparable and are designed to detect change points in time series data while at the same time addressing questions unique to their respective fields. Here, we present an application of the CP method to ambulatory, psychophysiological data. Our method is based on statistical control theory, and is designed to be a straightforward, simple application that is aimed at addressing a diverse array of psychophysiological questions. Our application accommodates multiple time series and corrects for the potentially confounding effects of extraneous third variables.
Two applications of CP analysis that could be particularly relevant to psychophysiology are: 1) the detection of the timing and magnitude of a single change point in a time series (single CP analysis), and 2) the detection of the timing, magnitude, and duration of multiple changes in a time series (multiple CP analysis) (see Lindquist et al., 2007). We discuss these applications below:
Single CP analysis is designed to detect a single sustained change in a physiological measure. For example, it can detect a physiological event that might be a trigger of a psychological phenomenon (e.g., an increase in PCO2 that might lead to hyperventilation and precipitate a PA; or a drop in blood pressure that might result in fainting in patients with blood phobia). Single CP analysis examines each time point in a time series of data to determine if the average level of the outcome variable up to a specific time point (β0, the mean of the distribution from time 1 to time t0) is significantly different from the average level of the variable after that time point (β1, the mean of the distribution from time t0 to the final time, T). The significance of the difference between β0 and β1 is evaluated using the two-sample pooled t-test (Montgomery, 2001, pp. 117–119). The time point at which the CP occurred is calculated using maximum likelihood estimation (MLE) (Basseville & Nikiforov, 1993). The MLE calculates β0, β1, and the time t0 (the MLE estimate of the time of the change point) to maximize the probability that the actual distribution of the data came from distributions with the calculated values for these parameters (that is, distributions with a mean of β0 that persists through time point t0, followed by a distribution with mean β1 which persists to the end of the time series). The formula for the maximum likelihood estimator of the change time can be found in the appendix of this manuscript.
Multiple CP analysis can be used to examine processes that might turn on and off (e.g., evoked potentials) or that might vary substantially over time (e.g., respiratory instability before PAs). While single CP analysis does not allow for the possibility that the data may return to their initial states, or for the possibility that the time series might change more than once, multiple CP analysis allows for these events. In multiple CP analysis, the time series is scanned from the beginning to detect the first CP (the first point at which the distribution changes from one level to another). Then, the remaining time series is scanned for the next change point. This process continues until no more CPs are detected
While Lindquist and colleagues (Lindquist et al., 2007; Robinson, Wager, & Lindquist, in press) have effectively demonstrated how CP analysis can be applied to complex functional magnetic resonance imaging data, the aim of this report was to describe an alternative approach to CP analysis applied to ambulatory, peripheral psychophysiological data, while still drawing from their important work. We analyzed a selection of peripheral parameters that are thought to play a central role in the onset of PAs. Additionally, we present the CP method within the context of accounting for time varying “confounding” variables, such as movement, and providing for comparisons between “experimental” time series (pre-PA) and those from a “control” time period (not preceding a PA). By doing so, we hope to expand the range of application of CP analysis, while at the same time addressing additional methodological issues that are relevant to its application.
To illustrate the application of CP analysis to ambulatory data, we chose data from a set of 24-hr cardio-respiratory recordings, during which real-life PAs were observed. The data seem particularly germane as they address a long-lasting debate about the relation of autonomic changes, in particular respiratory changes, to the onset of PAs (Roth, Wilhelm, & Petitt, 2005). Prominent respiratory theories of PD (Klein, 1993, Ley, 1985; for a review see Roth, Wilhelm, & Petitt, 2005) have posited respiratory instability as a trigger of PAs. More specifically, the hyperventilation theory of panic (Ley, 1985) postulates that a sudden or sustained reduction in PCO2 precedes the onset of a PA. On the other hand, the false suffocation alarm theory (Klein, 1993) assumes that an increase in PCO2 might be responsible for the development of a PA. Experimental studies have shown that changes in PCO2 in either direction can provoke panic-like symptoms (Papp, Klein, & Gorman, 1993; Wilhelm, Gerlach, & Roth, 2001; Hegel & Ferguson, 1997) that can accumulate into what are described as PAs. Even though experimental manipulation of certain aspects of physiology in the laboratory can create PAs, it remains unknown whether these particular physiological changes reflect actual precursors to naturally occurring PAs. Similarly, although the period preceding the onset of a PA may be marked by substantial fluctuations in physiology, it is unknown whether these fluctuations reflect underlying changes in physiology or whether they are merely just random fluctuations in physiology. Thus, the application of CP analysis could offer a much needed test of the predictions of such theories and vitally enhance our understanding on potential biological risk factors for PAs. However, since CP analysis does not address causality, we can only demonstrate whether or not the data are consistent with hypothesized theories.
The patients included in this study were selected from a larger sample of patients (N=43) that had undergone a 24-hr ambulatory monitoring prior to entering a biobehavioral treatment study for PD (for more details see Meuret, Wilhelm, Ritz, & Roth, 2008). Their inclusion for analysis was based on the presence of reported PAs during the 24-hr ambulatory monitoring. All patients suffered from panic disorder with (N=6) or without agoraphobia (N=5) (PD/A) (Structural Clinical Interview for DSM-IV, First et al., 1994) and were recruited through local advertisements. Patients were informed that the ambulatory recordings would serve as extended baselines to assess day and nighttime autonomic and respiratory functions. The sample included 4 men and 7 women, averaging 41.6 years of age (SD = 12.0, range 23–61). The majority of the participants were Caucasians (81.8%), 9.1% were African American and 9.1% were Asian. Patients had suffered an average of three years from PD/A (range: 0.5–16). The study was approved by the Stanford University Institutional Review Board and informed consent was obtained prior to the study from all participants.
Physiological parameters representing respiratory (PCO2, RR, tidal volume [VT]) and cardiac (HR) activity were selected, in addition to control parameters (motion (ACC), vocalization (SND)). Sixty minutes of continuous minute-sby-minute data preceding the onset of 13 naturally occurring, spontaneous PAs were analyzed. The attacks were captured during 24-h ambulatory monitoring in the patient’s natural environment. Of the 11 patients reporting a PA, two patients experienced two PAs.1 Patients indicated the onset of the PA by either pressing a marker button attached to the equipment, or by recording its inception in an electronic log that was part of the monitoring equipment (see Wilhelm, Alpers, Meuret, & Roth (2008) for detailed information on ambulatory protocols and measurement techniques). To verify the existence of significant event-related changes beyond naturally occurring random fluctuations, we contrasted the signal patterns of the data that preceded the PA with data from matched non-panic periods. Thus, each pre-PA epoch was matched with an alternate, panic-free period (non-PA) during the 24-hr recording period. These non-PA periods were randomly selected out of samples of periods which were subject to the following criteria: They had to match the duration of the corresponding PA (e.g., if the PA lasted 20 minutes, the corresponding non-PA period also lasted 20 minutes), and they had to be a minimum of 3 hours apart from the time of the PA. As in the case of the PAs, we extracted the 60 minutes of data preceding the matching non-PA period.
The data was collected using portable ambulatory devices, the Vitaport II System (Meditec, Karlsruhe, Germany) with attached respiratory inductive plethysmography bands (Ritz et al., 2002). End-tidal PCO2 was assessed using a portable, battery-operated capnometry device (Capnocount Mini, Weinmann GmbH, Hamburg, Germany) that samples exhaled gas through a nasal cannula.2 The capnometer was connected to the ambulatory recorder. Signals were then analyzed off-line and averaged for 1-min periods using an integrated set of biosignal analysis programs (for more details see Conrad, Isaac, & Roth, 2008) written by the fourth author [AC] in MATLAB 7.0 (Mathworks, Inc., Natick, MA).
In the following we outline the steps for using CP analysis on the above described longitudinal data set. The algorithm to perform the CP analyses was implemented in MATLAB by the second author [EZ]:
Physiological data is often impacted by events unrelated to the variables of interest. For example, previous research indicates that physical exertion and movement can substantially contribute to non-emotional changes in autonomic variables (e.g., HR activity) (Grossman, Wilhelm, & Spoerle, 2004). Thus, failure to control for physical activity when analyzing physiological data may result in misleading conclusions caused by changes in activity level rather than psychophysiological changes (e.g., Taylor et al., 1982; 1986). Similarly, speech can be a source of signal bias in certain respiratory measures (e.g., frequency and depth). Since it is possible that CPs might be found that are simply artifacts of movement or speech, controlling for both movement and vocalization data is essential.
We employed HLM to remove variations in the outcomes that might have been due to changes in these covariates (activity and speech). Although HLM is not well suited for detecting unknown change points, it is a viable approach to control for time varying covariates. HLM can easily handle 1) correlated data (e.g., longitudinal data), 2) missing observations, 3) time varying covariates, and 4) it allows for estimating the variance over individuals.
Level 1 of the HLM analysis consisted of the 60 data points from each of the 13 different PAs during the one hour preceding the PAs or during the matched non-PA period. The level 1 regression equation included movement and vocalization as predictors of the cardio-respiratory outcome over the 60 data points within each PA. The level 1 residual then became the value of the outcome after the effects of movement and vocalization were removed (i.e., the residual is the difference, at each data point, between the actual value of the outcome and the value that would have been expected given the level of movement and vocalization at that time point).
Level 2 consisted of the 13 different PAs (or non-PAs for the matched time series). The level 2 equation allowed these relations (and the intercept) to vary between individuals.
The residuals from the HLM analyses provided the adjusted values of the outcomes, corrected for movement and vocalization, at each data point for each PA (or non-PA). We combined the data from each individual PA to form an “aggregated” time series, in order to detect CPs for the “average” PA. Hence, we averaged the value of each outcome at each data point across the 13 PAs, yielding an “average” time series for each outcome measure. We then performed the CP analyses on these aggregated time series. This approach is similar to that of other investigators that examined the physiological profile of PAs (e.g., Goetz et al., 1993), who used analyses of variance to test changes in the average HR of patients immediately prior to PAs. It is also similar to the group control chart approach used in statistical quality control (Montgomery, 2008), in which limits are calculated using the mean of multiple streams of data, and to “pooling” data in econometrics. This approach is different from the multilevel approach presented in Lindquist et al. (2007), in which time series from the individual PAs are considered nested within individuals. Please refer to the Discussion for a comparison of the present approach and the one by Lindquist and colleagues.
To minimize the possibility of erroneously detecting random, transient changes as CPs, we added two constraints to the CP analysis. First, in our single CP analysis we only examined changes that persisted through to the onset of the PAs. We reasoned that physiological changes that returned to earlier levels prior to the onset of the PAs were unlikely to be a cause of, or directly connected to, the PAs. Second, to avoid detecting transient changes that may only represent random variations in the physiological measures, we added the constraint that any CP must last a minimum of 5 minutes. We choose this time restriction to insure the detection of substantial change instead of randomly occurring changes. Finally, to further confirm that the detected CPs were not merely typical random variations, we performed identical CP analyses for the matched 60 minute non-PA time series. The absence of CPs during the non-PA epochs strengthens our conviction that detected CPs are indeed related to the evolving PAs and not simply due to random variation. These criteria were designed to conservatively insure that we would not mistakenly misinterpret random events as precursors of the PA.
Finally, since the minimum required CP length was set at 5 minutes, we only examined minutes 5 to 55 for change points. However, examining so many time points for change will inflate Type 1 error, requiring an adjustment to the significance tests. The Bonferroni correction is a conservative method to compensate for multiple significance tests. Since we tested 50 time points for four different physiological indices, a Bonferroni correction of 4*50 to the standard p level of .05 was used. Thus, CPs were only considered significant if they reached a significance level of p < .05/200 (i.e., p < .00025).
We first used HLM to test whether there were baseline differences between the pre-PA period and the non-PA period to verify that the matched non-PA period was initially equivalent (over the first 10 minutes) to the pre-PA period.3 We found no differences between pre-PA and non-PA on initial levels of any of the measured variables: RR (Meanpre-PA [Mpre-PA])= 18.18, Mnon-PA = 17.94, pdiff = .78, σ2 (level 1 variance in HLM) = 6.15), VT (Mpre-PA = 537.67, Mnon-PA = 501.50, pdiff = .37, σ2 =10682.2), PCO2 (Mpre-PA = 28.16, Mnon-PA = 28.28, pdiff = .92, σ2 = 4.50), or HR (Mpre-PA = 86.02, Mnon-PA = 87.29, pdiff = .53, σ2 = 12.83).
We then performed the HLM analyses which used the covariates to control for the effects of movement and vocalization on the outcome variables. As expected, greater movement was related to higher RR (b = 28.49, t (12) = 3.92, p < .005), higher VT, (b = 1718.69, t (12) = 8.59, p < .001), and faster HR (b = 125.07, t (12) = 6.73, p < .001), within subjects over time. Furthermore, it was also related to lower PCO2 (b = −9.00, t (12) = −3.22, p < .01). Similarly, vocalization was related to higher VT (b = 2.28, t (12) = 2.30, p < .05) and lower PCO2 (b = −.004, t (12) = −2.46, p < .05).4 The residuals for each outcome at each time point for each PA was calculated, and used in the CP analyses.
The single change point analysis was used to detect significant and sustained changes in the aggregated time series leading up to the onset of the PAs. Figure 1 (labeled Pre-PA) illustrates the findings from these analyses. The analyses resulted in the following sequence of significant changes preceding the PAs: 1) a significant increase in VT about 21 minutes before onset of the PAs, p <.0001, followed by 2) a significant decrease in RR about 19 minutes prior to the PAs, p <.0001, and 3) a significant increase in PCO2 about 16 minutes before the onset, p <.0001. HR did not significantly change during this time period (a CP at 42 minutes before the PA had p = .001, not significant after the Bonferroni correction). To verify that these changes were not merely random or characteristic of changes that regularly occur in these variables, we performed the single CP analyses on the matched non-PA time series. These analyses showed only one significant CP (for PCO2 at time point 8, p = .0001) for the non-PA time series (Figure 1, labeled Non-PA).
Although we controlled for the effects of movement and speech on the above physiological measures, thereby removing any changes due to these control variables, we performed CP analyses on movement and speech to further understand changes that may have occurred in these parameters during the pre-PA period. No significant single CPs were detected for movement or speech either during the pre-PA or non-PA time period.
As described above, the purpose of the multiple CP analysis was to detect substantial changes prior to, but not necessarily sustaining through, the onset of the PAs. This analysis consistently found activation across all physiological measures, showing significant CPs during the 60 minutes before the PA (see Figure 2, labeled Pre-PA) on each of our physiological measures. To determine if this high degree of instability reflected typical variations in these physiological measures, we again compared the results of the multiple CP analyses of the pre-PA data to the non-PA time series. The results for non-PA revealed a very different pattern: in the non-PA time series: only one of the four physiological measures (PCO2) showed a CP, and it showed only one CP (Figure 2, labeled Non-PA).
The changes detected in the multiple CP analyses for the pre-PA series showed an interesting pattern of activation and deactivation (see Figure 2). An initial spike in HR (p <.0001) about 47 minutes before the onset of the PAs was followed soon after by a substantial drop in HR (p = .0001). Then, HR increases again about 30 minutes before the beginning of the PAs (p = .0001). This spike in HR (which occurred before changes in activity) seems to be followed by a cascade of respiratory changes. About 27 minutes prior the onset of the PAs, PCO2 dropped (p < .0001). A few minutes later, VT dropped (p = .0001), followed by a compensatory increase in RR (p = .0001). These changes culminated in a final increase in PCO2 (p = .0001) about 15 minutes before the onset of the PA. CPs were also detected for physical activity. There was a decrease in activity about 27 minutes before the onset of the PAs, followed by an increase in activity about 16 minutes before the onset of the PAs. There were no CPs detected for physical activity during the non-PA period and none detected in speech for either the pre-PA or non-PA time series.
The purpose of this study was to investigate the application of CP analysis for detecting discrete physiological events using cardio-respiratory data that was recorded prior to the onset of self-reported PAs. Two sets of analyses, single and multiple CP analyses, were performed to demonstrate the steps required to satisfy the multiple criteria for the detection of CPs. Our first set of analyses was designed to detect changes that persisted through the beginning of the PAs (and lasted at least 5 minutes), reasoning that changes which return to baseline before the onset of an event are unlikely to be directly connected to (or responsible for) the onset of that event. It needs to be said that persistence through the onset of the event cannot be regarded as evidence that the change caused the event (e.g., the occurrence of the PAs). But, it seems doubtful that a change that returns to baseline long before an event is likely to be the proximal cause of that event, although it could be part of a chain of events leading up to the event.
As demonstrated by the single CP analyses, the onset of the PAs was preceded by sustained changes in respiratory physiology, compared to the time-matched non-PA control periods, where only one isolated decrease in pCO2 was detected. The observed respiratory pattern in the pre-PA time period is consistent with respiratory theories of panic disorder discussed in the literature. For instance, Klein (1993) proposed that spontaneous PAs are triggered by an oversensitive “monitor” in the brain stem that sends a “false suffocation alarm”. The alarm is triggered by certain conditions, such as raising levels of PCO2. The observed pattern of physiological change detected by the single CP analysis follow a sensible sequence of respiratory changes, namely a decrease in VT followed by a compensatory increase in RR, which in turn was succeeded by an increase in PCO2. This increase in PCO2 indicates that the previous decrease in VT apparently overrode the effect of increased RR. This finding supports recent reports on the role of depth of breathing rather than rate of breathing as a contributing factor to changes in PCO2 in patients with anxiety disorders (Ritz, Wilhelm, Meuret, Gerlach, & Roth, 2009; Ayala, Meuret, & Ritz, in revision).
Our second set of analyses searched for all significant points of change (or systematic instabilities) in the time series, with the only criterion being that the change was maintained for at least 5 minutes. This criterion was imposed to avoid detecting transient and/or random changes in these highly variable data streams. We found that all four of our physiological measures displayed discrete changes preceding the onset of the PAs. With the exception of a change in PCO2, such changes were not found in the matched control periods. While the findings of the multiple CP analysis reinforced the findings of the single CP analysis that suggested the role of ventilation depression (reduction in VT and elevation in PCO2) before panic onset, the additional change points identified by multiple CP analysis added a dynamic perspective to these findings. The PCO2 increase was preceded by a sustained drop in PCO2, which may have, in turn, been a result of the preceding instabilities in RR (although these instabilities were not found to be significant).
The results from both CP analyses (single and multiple) seem to indicate that, before PA onset, various physiological systems exhibit deviations from a set point. These deviations may be a part of a cascade of systemic dysregulation escalating into a PA. Thus, our CP analysis appears to be an effective method to detect events of unknown time or duration in a time series of highly variable psychophysiological data. These findings are especially remarkable given the very conservative Bonferroni correction applied to the statistics. In identifying up-and down-regulations of physiological processes, CP analysis can help demonstrate dynamic temporal patterns of dysregulation that could be relevant to models of psychopathology, such as respiratory theories of panic disorder in the present example, or that could be applied to more long-term time series data detecting precursors of depressive episodes.
While these findings do not prove that the observed patterns of instabilities were the cause of the PAs, they may reflect the presence of physiological or psychological processes that can lead to a PA. For example, the changes detected by the single CP analysis, which persisted through the beginning of the PAs, could have been sufficient to trigger interoceptive awareness (either conscious or unconscious). Such awareness may have then triggered further physiological arousal. Similar feedback loops between symptoms and cognitions have been discussed extensively in the PD literature (e.g., Clark, 1986). Although the PAs reported in this study are categorized as spontaneous (in that there was little or no outside provocation), the pattern of physiological instability preceding the onset does not appear random or spontaneous but suggests an accumulation of marked activity beginning roughly 45 minutes prior to the onset.
As noted earlier, we chose to aggregate the 13 PAs together and analyze the resultant “mean” time series, rather than adopting the multilevel approach developed by Lindquist et al. (2007). The present approach is similar to the group control chart approach used in statistical quality control (Montgomery, 2008), and the “pooling” of data in the econometrics literature. In our view, our approach is accessible, and the change points estimated by it are the same as those estimated by the Lindquist method when the number of observations per individual is equal and when the errors of the individual subjects are similar. On the other hand, the Lindquist et al. approach makes fewer assumptions concerning the distributions of data from the individuals, could allow for predicting individual differences in change points, and is more powerful and less conservative than using the Bonferroni correction. Further, Robinson, Wager, & Lindquist (in press) illustrate a procedure that does not assume a specific functional form for the “activated state” and allows for the possibility that some subjects may show no activation. We offer our “averaging” approach as a straightforward method of calculating CPs when individual differences in the CPs are not the focus.
Despite the ability of CP analysis to detect changes of unknown time and duration, there are a number of limitations in the present analysis that deserve consideration. First, we analyzed the aggregated time series across multiple PAs. As in all psychological processes, the individual PAs that comprised this data showed a wide variety of patterns. Thus, the timing of these CPs may not precisely reflect the timing of the changes that occur in each individual. Second, our CP analysis did not statistically compare CPs in pre-PA to the non-PA time series, so we cannot conclude that these time series were significantly different from each other. Third, we cannot conclude that these changes were directly related to the onset of the PAs. They may rather have been the result of the same process that produced the PAs, instead of being a causal link in the process. Finally, the complexity of the multiple CP analysis makes it difficult to calculate the most appropriate Bonferroni correction. However, adopting a p level of .05/200 (.00025) in our study appears to have face validity since we were able to detect significant CPs in the pre-PA epoch but found few in the non-PA epoch.
We illustrated two applications of CP analysis for examining peripheral psychophysiological data. The analysis provided a unique test of current theories proposing that respiratory changes precede the occurrence of PAs. The findings demonstrate that a cumulative pattern of cardio-respiratory instability does occur prior to impending PAs. This observation adds significantly to our understanding of this clinical phenomenon and can inform treatment approaches. Indeed, alterations of sustained levels of hypocapnia have proven to be successful in reducing panic symptoms (Meuret et al., 2008) and are linked causally to improvement (Meuret et al., 2009; Meuret et al., revisions). Future studies using CP analysis are needed to further validate the detected patterns of cardiorespiratory instability and to examine their specificity. Other applications of CP analysis include the examination of the physiological profile of subtypes of PAs (e.g., situationally-bound, nocturnal) and the investigation of changes occurring prior to fainting in blood phobia patients (Ayala et al., revisions). CP analysis of time-series data, whether recorded in the patients’ natural environment by ambulatory monitoring, or obtained in a laboratory setting, is a promising tool for detecting discrete changes in autonomic patterns even in the context of complex and highly variable random “noise”.
While the time series in our analyses involved aggregated data from multiple PAs, the CP technique can also be used to analyze individual time series data. In such, the following issues await further exploration and refinement: Individual data are inherently more variable than aggregated data. Thus, the power to detect individual CPs is expected to be lower than in aggregated data. It is feasible that techniques such as exponential smoothing could be used to reduce the random variation inherent in the data (see Lindquist et al., 2007). Such methods would require the proper selection of the smoothing technique and smoothing parameters, which can be challenging without prior knowledge of the nature of the discrete event to be detected.
In addition, the handling of missing data in individual time series requires further exploration. It is not unusual for a longitudinal series of psychophysiological data to have “skips”, where data is missing because of equipment failure, artifacts, missed assessments, etc. Such missing data, if missing at random, does not pose a particular problem for aggregated data. For example, in our study if there was missing data from some subjects scattered throughout the 60 minutes before the onset of the PAs, we still had the aggregated data from the other subjects at each time point to provide data at all 60 time points. However, when analyzing data from a single PA, missing data causes “skips” in the time series. It is unclear whether traditional choices such as imputing the last valid observation into the time series, interpolating the missing data, or “closing up” the skips in the time series are viable options for treating the missing data CP analyses. Closing the skips can lead to discontinuities in the series due to the amount of time separating the adjacent “closed-up” time points, which can be mistaken for CPs. Interpolating the data may not be a better choice, since this can lead to underestimation of the variance of the distribution, thereby increasing Type I error.
It is therefore clear that further research is needed to understand and develop effective applications of CP analysis to individual data series with missing data. But despite its challenges, the application of CP analysis to individual time series could be particularly important since there are presently few statistical tools available to analyze such data and since many questions relating to temporal dynamics of biopsychological organismic regulation still need to be addressed.
In conclusion, CP analysis is a very promising method to detect changes in longitudinal time series data, with multiple areas of application. The method reaches far beyond traditional methods of statistical analysis that are not well suited to detect discrete changes in time series. Further testing and refinement is needed to fully develop its promising capabilities.
The exact procedure for the change point analysis was as follows:
For detecting a single change point at time = t0, the data time series was modeled as:
where yi is the value of the outcome (measurement) at time i, ui is the measurement noise (error) at time i, βj is the underlying data, j = 0 if i < t0, and j = 1 if i ≥ t0. We want to know if there is a significant change point in the data such that β0 ≠ β1, and, if there is such a change, we want to estimate the change time t0. We assume that u = [u1, …, uT] is normally distributed with mean zero and covariance matrix ∑. Then the maximum likelihood estimator of the change time point is
where y = [y1, …, yT]’, ‘ denotes transpose, and (k) is a vector where the first to the (k-1)th elements are 0 and the kth to the Tth elements are 1 where
The covariance matrix ∑ can be estimated using the non PA data.
If we further assume that ui is identically and independently distributed, i.e., each ui is normally distributed with zero mean and the same variance (Basseville & Nikiforov, 1993, pp. 66–69), then the above formulae can be simplified to:
which further simplifies to:
Thus, 0 and 1 are the sample means of the distribution, before and after the change point respectively. One uses (3) to test each time point, and the time point that maximizes the argument is the CP. That time point is the point that minimizes the sum of the two summations in (3). The first summation is the sum of squared deviations around β0 (the sum of squares prior to t0) and the second summation is the sum of squared deviations around β1. Thus, the time point that produces the smallest sum of squared deviations around the sample means before and after that point is the change point.
In our case, the errors (ui) are not independent. Therefore, we used (2) to calculate our CPs. After the change point was identified, we performed a permutation analysis to examine the significance of the CP since the errors were autocorrelated. In the permutation analysis, we randomly shuffle the data and calculate the maximum likelihood of the shuffled data series. We repeat it multiple times (we chose 10000 times in our analysis) to obtain many likelihoods, and use all the likelihoods to form an empirical distribution. The p-value of the CP is then determined by where the likelihood associated with the CP lies in this empirical distribution.
If we assume the errors are independent, we can instead use a 2 sample pooled t-test to examine the significance of the CP using the following formula from Montgomery (2001, pp. 117–119):
where is the pooled sample standard deviation.
This Z statistic follows a Student’s t distribution with T-2 degrees of freedom, and hence the p-value and confidence interval can be calculated according to the Student’s t distribution.
The detection of multiple change points is based on the algorithm for detecting a single change point. We assume that each change lasts a minimum duration of a minutes (a = 5 minutes in our study). Initially, a single change point test is carried out for the data segment from time 0 to time a. If there is no significant change point detected, we increased the length of the data segment by 1, i.e., from time 0 to a + 1, and then performed the single change point test again for this segment. This step was repeated until a significant change point was detected, or the data segment reached the end of the time series. In the latter case, the time series had no change points. In the former case, a new time series was formed by taking the data segment from the latest change point to the end of the original time series, and the new time series was tested using the same approach as above. The algorithm terminates either when no more change point can be found or when the remaining time series from the latest change point is less than a minutes.
Our algorithm to detect single and multiple CPs using MATLAB is available upon request.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1A total of 15 PAs were reported. We excluded the two nocturnal PAs from the analysis due to missing control variables (vocalization and motion).
2Mouth breathing can be a concern during end-tidal PCO2 ambulatory capnometry recordings, in that it can result in lower levels of PCO2. Precautionary steps were taken to minimize fluctuation in PCO2 due to mouth breathing. Patients were instructed to remind themselves to breathe through their nose as much as possible. In addition, verbal reminders in the form of taped instructions for standard tests throughout the days were provided. Furthermore, recordings were rescheduled if nose breathing was limited (e.g., stuffy or runny nose). Occasional mouth breathing in this study was likely to affect both the pre-PA and non-PA data in an equally random fashion. Furthermore, no significant differences for speech were observed for the pre-PA and non-PA recordings. Most importantly, the first step in the data analysis was the visual inspection of the analog PCO2 curves. Only breaths with full plateaus were considered in the data analysis. The capnometry device used in this study adheres to international accuracy standard (Biedler et al., 2003).
3In this HLM analysis, the first 10 assessments of the pre-PA and non-PA periods comprised the level 1 data that was nested within subjects at level 2. A dummy variable which coded pre-PA/non-PA was included as a level 1 predictor of outcome, and motion and vocalization were also included as level 1 control variables (to avoid detecting differences between pre-PA and non-PA due merely to differences in these control variables). The significance of the dummy variable indicated whether there were significant differences between the pre-PA and non-PA data.
4Vocalization is usually accompanied by slight hyperventilation (Ritz et al., 2002). Lower PCO2 during physical activity could result from panic patients hyperventilating in response to activity-induced physical symptoms. Thus, the presented CP analyses employ the residuals from the HLM analyses in order to control for the effects of speech and movement on the physiological measures.