MultiCenter Interrupted Time Series Analysis: Incorporating Within and Between-Center Heterogeneity

Introduction

Interrupted time series (ITS) design involves repeated measurements of an outcome at several time points before and after implementation of interventions or programs and can be designed to investigate the impact of interventions and programs in healthcare settings. Segmented regression (SR) is the most common statistical method used in the analysis of ITS data.^1–4 A recent scoping review, aimed at identifying methods for analyzing ITS data, showed that SR models were the most utilized methods, having been used in approximately 45% of the included studies to analyze their ITS data.⁵ These methods have been applied in various health research areas (eg clinical research, public health, and health services) to assess the impact of interventions on patient important outcomes or clinical practice.^6–12

SR is a special case of multiple linear regression with an indicator variable representing the intervention periods, a continuous variable representing the time at which observations are taken, and an interaction variable.¹³ For ITS analysis employing SR, a regression line is fitted to each segment of the time series, for instance, the pre- and post-intervention periods.¹⁴ This approach assesses the magnitude of the effect of the intervention on the outcome of interest. Using this approach, both the immediate effect of the intervention as well as the effect of the intervention over time can be evaluated. SR analysis controls for potential confounders (such as age and gender), for example, by including covariates in the model to account for the confounders.¹⁵ The model can also be used to investigate seasonal or cyclical patterns that occur over time, as well as used to assess the presence of autocorrelation in the ITS data by including autoregressive parameters in the SR model.^3,16

Despite the advantages, the SR model has a major limitation when it is applied to aggregated data.^2,17,18 Aggregated data refer to summary statistics (eg mean, percentage, median) calculated across study participants. In health research involving ITS, aggregated data are very common since interventions are often implemented to evaluate the effect at healthcare facility level and data at a given time point are often summarized across participants within the healthcare facility. For multisite ITS studies, further aggregation occurs, where the data are summarized across several sites (eg hospital units) to evaluate the overall impact of interventions. Additionally, administrative routine data collected from different subjects across different regions are increasingly being used in ITS studies to perform post-hoc evaluation of nationwide policies and programs.⁵ Such data are usually only available in aggregated forms.

Aggregated data are a major cause of heteroskedasticity, which is a violation of the assumption of constant error variances, required during inference involving regression models.¹⁹ When data are aggregated, variability is introduced due to factors such as the differences in participant characteristics, sample size, and settings. Hence, the variance across time, and consequently the error variance, is no longer constant. To account for the variability introduced by aggregating data, Gebski et al²⁰ proposed the pooled time series analysis method, where the intercepts and slopes are pooled from individual site analysis, using meta-analytic approaches, to calculate the overall effect of the intervention. Their method may not be optimal for all types of aggregated data, since the variability across participants within the same site was not accounted for. Further, the use of summary estimates (intercepts and slopes) instead of using all the available information from the data may lead to loss of power, and hence suboptimal results.^21,22

In this study, we propose the extended weighted SR method (wSR), with the aim of incorporating both participant-level and site-level variability. The precision, accuracy, and statistical significance of our method are compared to the traditional SR method as well as the pooled time series method proposed by Gebski et al using empirical data.

Materials and Methods

Empirical Data

We use data from the Mobility of Vulnerable Elders in Ontario (MOVE-ON) project for our empirical evaluation.^23,24 The original MOVE-ON study involved implementation and evaluation of the impact of an evidence-based educational intervention, which aimed to promote early mobilization and prevent functional decline among older patients admitted to acute-care academic hospitals across Ontario. An interrupted time series design was used to evaluate the impact of the intervention. The study was conducted across 14 hospitals consisting of 30 units that provide care to inpatients aged 65 and older. Data were collected over a period of 38 weeks, with a 10-week pre-implementation period followed by an 8-week implementation period, where the intervention was rolled out, and a 20-week post-implementation period was also considered. The primary outcome of the study was the mobilization status of patients, who were assessed on twice-weekly visual audits, which occurred three times per day (in the morning, at lunch, and in the afternoon). Patients were considered mobilized if they were observed out of bed (mobility score > 2).^23,24 Further details of the study can be found elsewhere.^23,24 In this paper, we will consider the primary outcome and present the proposed method using these data as an illustration. We also compare our method with previous methods empirically.

Weighted ITS Analysis

Consider the primary outcome, patient mobility, in the MOVE-ON project. Within each site, and at each time point, patient mobility was summarized as the percentage of patients who were observed (during audits) out of bed at least once in a day. That is, if we let represent the mobility status for patient at time point , then the percentage of patients who were observed out of bed () at a time point is calculated as:

(1)

where is an indicator function dichotomizing the primary outcome into two categories, whether a patient is out of bed or not out of bed; and N represents the total number of patients observed at each time point. To examine the impact of the intervention on patient mobility, a SR of the form described in equation (2) below was fitted and comparisons in level and trend of mobility were made among the pre-, during, and post-implementation periods.²

(2)

where represents the time at which the outcome is measured, and Int1 and are indicator variables representing the intervention periods. Interactions between time and intervention periods are also included in the model. The parameter estimates the baseline percentage of patients who were out of bed; represents the slope (trend) of mobilization prior to the intervention; and and represent the changes in intercepts and slopes between the pre- and during intervention periods; while and represent the changes in intercepts and slopes between the pre and post intervention. The model can be presented in matrix format as:

(3)

where Y denotes a vector representing outcomes measured at different times; is a matrix of the independent variables such as time and intervention; denotes a vector representing the regression coefficients (parameters); and is the error vector.

Recall that at each site, the data at each time point are aggregated across patients, which also involves summary statistics based on different sample sizes. Further, being a multisite study, data are also aggregated across sites. Thus, there are two levels of heterogeneity introduced; within site (ie across patients within hospitals and units) and across sites (across hospitals). Due to this unaccounted heterogeneity, the SR approach will be associated with increased bias and decreased precision when used on aggregated data. It has also been previously established that the violation of the homoscedasticity assumption (due to the difference in levels of variability) leads to inflated type I error and decreased power.^2,25 To overcome this limitation, we propose the extended wSR approach, where the weighted least squares method is used to estimate the parameters of the model in equation (3). The estimators are given by:

(4)

where represents a diagonal matrix, which consists of the vector as its diagonal elements. The weights (w) are calculated as the inverse of the sum of the within-site and between-site variances. That is, if represents the variance within a site at time point , represents the variance between sites, and represents the variance associated with the random error, then the weights are given as:

(5)

The variance of the weighted regression estimates is calculated as:

(6)

The standard errors and confidence intervals can therefore be constructed accordingly, using established formulas found in the literature.²⁶

Empirical Evaluation

Our empirical evaluation consists of comparative analysis involving our wSR approach, the traditional SR approach, as well as the pooled analysis (PA) proposed by Gebski et al.²⁰ To examine the performance of the proposed method to changes in the number of sites and various levels of within- and between-site heterogeneity, we performed analyses for different scenarios, where we varied the number of sites included, the amount of within-site variability, and the magnitude of between-site heterogeneity (low, moderate, and high) as measured by the statistic used in the meta-analysis literature to quantify the amount of between-study heterogeneity.^26,27 The heterogeneity is defined as low when moderate when , and high when For all scenarios considered, the results (estimates of slopes and intercepts) and the associated 95% confidence intervals (95% CIs) obtained from the three methods were compared to assess the precision and accuracy of the estimates respectively. For time series data, the Durbin–Watson statistic is the most common measure used to evaluate the presence of autocorrelation in the data.²⁵ However, the data used in this study are from a previous project where the authors reported there was no significant autocorrelation in the data. Hence, we did not perform any further assessment of autocorrelation. All analysis was performed using R version 3.4.3 statistical software.

Results

During the 38-week study, a total of 125,025 observations were made from a total of 14,540 patients across the 14 sites, of which 3943 patients were audited pre-intervention, 3216 during intervention, and 7381 patients post intervention. For each site, the total number of patients per time point ranged from 30 to 146.

For analysis involving overall patient mobility across the 14 hospitals, the estimates of the regression coefficients obtained using each of the three methods are presented in Table 1. The wSR estimates have the narrowest confidence interval for all estimates, while the traditional SR has the widest confidence interval for most estimates. Nonetheless, the parameter estimates are relatively comparable for the SR and wSR methods in most cases.

Table 1 Estimates for Intervention Effects and 95% CIs Obtained Using Segmented Regression (SR), Pooled Analysis (PA), and Weighted Segmented Regression (wSR)

All of the models produced a negative coefficient for the week () variable, indicating a decline in patient mobilization prior to the educational intervention. The positive coefficient of the intervention () variable implies that there was an immediate increase in the percentage of patient mobilization during the implementation of the intervention (Figure 1; Table 1). Similarly, there was a substantial increase in the weekly percentage of patients who were mobile during the intervention period, = 0.64 (95% CI: −0.55, 1.83), from the wSR estimate. We again observed that although all three models showed an increase in the trend of patient mobilization during the intervention period, this increase was not statistically significant (Figure 2, Table 1).

Figure 1 Forest plot for change in level of mobilization during the implementation of the intervention compared to the pre-intervention period.

Figure 2 Forest plot for change in trend of mobilization during the implementation of the intervention compared to the pre-intervention period.

Considering the change in level from the pre-intervention phase to the post-intervention phase, there was a significant increase in mobility, where 13% more patients were mobile in the post-intervention period compared to the pre-intervention period (Figure 3, Table 1). However, this difference was not statistically significant using both the traditional SR and PA methods. For all of the methods, the weekly percentage of patients who were mobile was shown to have increased by approximately 0.1% per week post intervention, albeit not significant (Table 1).

Figure 3 Forest plot for change in level of mobilization after the implementation of the intervention compared to the pre-intervention period.

From Figures 1–3, there is substantial variability within sites and heterogeneity across sites. To perform more extensive comparative analysis with varying levels of heterogeneity, we performed subgroup analyses. The results for these subgroup analyses are presented in Table 2 as well as Figures 4–10.

Table 2 Subgroup Estimates for Intervention Effects and 95% CIs Obtained Using Segmented Regression (SR), Pooled Analysis, and Weighted Segmented Regression (wSR)

Figure 4 Forest plot for change in level of mobilization post intervention for high between-site heterogeneity.

Figure 5 Forest plot for change in level of mobilization during intervention for moderate between-site heterogeneity.

Figure 6 Forest plot for change in trend of mobilization during intervention for moderate between- and within-site heterogeneity.

Figure 7 Forest plot for change in trend of mobilization post intervention for low between-site heterogeneity.

Figure 8 Forest plot for change in level of mobilization post intervention for low within- and between-site heterogeneity.

Figure 9 Forest plot for change in level of mobilization post intervention for moderate within- and between-site heterogeneity.

Figure 10 Forest plot for change in trend of mobilization post intervention for high within-site variability.

For analysis involving four sites, when the between-site heterogeneity was very high, , the wSR method produced the most precise estimate compared to the PA method and the conventional SR method (Figure 4), where the length of the confidence intervals for our estimates are significantly smaller and hence most precise compared to estimates based on the traditional SR and the PA methods (Table 2). The same conclusions were drawn when there was moderate between-site heterogeneity (Figure 5). On the other hand, when there was substantial within-site heterogeneity and moderate between-site heterogeneity, the wSR and the PA methods gave results that were similar in precision whereas the traditional SR still had the worst performance (Figure 6). When there is low between-site heterogeneity, , but substantial within-site variability, the wSR method provided the most precise estimates while the PA and the SR methods produced similar results (Figure 7). Interestingly, in the instance when the within- and between-site heterogeneity was low, the pooled method produced the most precise estimates, while the traditional SR performed significantly worse (Figure 8).

Similarly, when the number of included sites was increased to eight, the wSR method produced estimates with the narrowest confidence intervals in most cases (Table 2). When the between-site heterogeneity was at least moderate, , the wSR method gave the most precise estimates compared to the other two methods (Figure 9). Interestingly, we observed that in some cases, where there was substantial within-site variability and moderate to high between-site heterogeneity, the PA method produced estimates associated with the widest confidence intervals (Figure 10).

Discussion

In this study, our aim was to extend the traditional SR approach, which is most commonly used in applications involving ITS studies and provide a methodological framework that allows the incorporation of within-site (between participants) and across-site (between healthcare facilities) heterogeneity. We performed extensive empirical evaluations by creating various scenarios using a multisite ITS study and compared our proposed method with that of the traditional SR and the PA method previously published.

Overall, our proposed wSR method produced estimates with the narrowest 95% confidence intervals for most of the scenarios considered, indicating that our method led to increased precision by incorporating two levels of heterogeneity in the data. The PA approach produced comparable estimates with that of the wSR in some scenarios, while the conventional SR method had the widest 95% confidence intervals compared to the other two methods. The observed wider confidence interval for the pooled method in situations with substantial within-site heterogeneity and a relatively high number of sites agrees with the reported limitation of the pooled analysis method.²⁰ However, further investigation, preferably using extensive simulations, is required to understand performance limitations for the pooled method in comparison with the weighted analysis.

For almost all of the scenarios considered, our proposed method produced the most precise estimates, with narrower confidence intervals, hence allowing meaningful differences to be detected. This finding also indicates that our method led to increased statistical power. Finally, we would like to highlight that the extensive empirical evaluations performed in this paper lay the groundwork for further extensive study involving simulations to establish more extensive performance characteristics of our method and compare performance with the other methods in terms of well-established performance measures such as bias, mean square error (MSE), type I error rate, and statistical power.

Conclusions

In conclusion, the findings in this study show that accounting for participant variability and differences across healthcare facilities in the analysis of ITS data leads to increased precision and statistical power. The results also show that simply pooling slopes and intercepts from site-level analysis has performance limitations. Taking into consideration that most ITS studies are conducted in healthcare settings and involve aggregated data, we believe this study provides findings that encourage researchers to consider differences in the participant populations and across healthcare settings. Our study also highlights the need for researchers who collect administrative routine data, which are usually presented in aggregated forms, to also report the variability associated with such data for meaningful statistical analysis to be performed.