Non-Parametric Generalized Additive Models as a Tool for Evaluating Policy Interventions

The interrupted time series analysis is a quasi-experimental design used to evaluate the effectiveness of an intervention. Segmented linear regression models have been the most used models to carry out this analysis. However, they assume a linear trend that may not be appropriate in many situations. In this paper, we show how generalized additive models (GAMs), a nonparametric regression-based method, can be useful to accommodate nonlinear trends. An analysis with simulated data is carried out to assess the performance of both models. Data were simulated from linear and non-linear (quadratic and cubic) functions. The results of this analysis show how GAMs improve on segmented linear regression models when the trend is non-linear, but they also show a good performance when the trend is linear. A real-life application where the impact of the 2012 Spanish cost-sharing reforms on pharmaceutical prescription is also analyzed. Seasonality and an indicator variable for the stockpiling effect are included as explanatory variables. The segmented linear regression model shows good fit of the data. However, the GAM concludes that the hypothesis of linear trend is rejected. The estimated level shift is similar for both models but the cumulative absolute effect on the number of prescriptions is lower in GAM.


Introduction
Although well-conducted randomized control trial experiments (RCTs) provide the most reliable evidence on the effectiveness of interventions, these are not always feasible for policy intervention analysis. RCTs involves randomly allocating participant units into two groups: the treament group which includes the participants who receive the intervention, and the control group. Selection bias and confounding variables are minimized due to randomization. Thus the differences between groups can be attributed to the intervention. However, when it comes to measuring the effect of policy interventions, there may be obstacled to the use of RCTS, such as economic obstacles (impact evaluation are costly) or political constraints (to give services to some groups and not to others can generate conflicts).
As alternative to RCTs, the interrupted time series analysis (ITSA) offers a quasiexperimental research design to measure the effect of an intervention when randomization is not possible [1]. In an ITSA, the observations on the outcome before and after the intervention are used to test immediate and gradual effects of the intervention. ITSA has been used in various fields, such as financial economics [2], health policies [3] and regulatory actions [4], to name but a few.
Segmented regression analysis is the recommended approach for analysing data from an ITSA [5,6]. It requires data which are are evenly spaced and have enought information before and after the intervention. Segmented regression analysis of interrupted times series data allows us to estimate immediate and gradual effects of the intervention on the outcome. Segmented regression analysis also allows us to assess whether factors other than the intervention could explain the change, controlling for factors such as seasonality or autocorrelation.
Segmented linear regression models have been the most widely used in practice. These models allow estimating the changes both in level and trend that follow an intervention. However, the assumption of linearity often may hold only over short intervals. Trends before and/or after the intervention may follow non-linear patterns, such as curvilinear trends. Some non-linear effects can be accommodated in linear models by using polynomial trends [7] or transformations of the dependent variable such as the logarithmic transformation [8]. Other non-linear trend structures may require other alternative models such as Box-Jenkins models [9]. However, the greater complexity in the specification and interpretation of this type of model has led to their less use [10].
Generalized Additive Models (GAMs) have been proposed as an alternative to characterize general non-linear regressions, without requiring the analyst to prespecify the form of the non-linear relationship [11]. Recently, Sullivan et al. [12] showed that GAMs can be useful for characterizing trends in longitudinal data collected in Single Subject (SS) designs. SS design is most often used in applied fields of psychology, education and human behaviour. SS design is a research design in which a single individual, or very small samples, is analyzed during a baseline period followed by an intervention that can change the outcome. This period can be followed by a return to baseline due to the removal of the intervention. This design can lead to a nonlinear relationship between time and outcomes that may not be easily detected using linear models.
In this paper, we assess whether the use of GAMs can be extended to estimate the impact of an intervention in any ITSA. Simulated data are used to evaluate the performance of the segmented linear regression models and GAMs in estimating the level change and the cumulative effect of an intervention. Data were simulated assuming linear and nonlinear trends and the mean squared error (MSE) and the mean percentage error (MPE) are used to compare both methodologies. An illustrative example with real data is shown where the impact of the 2012 Spanish cost-sharing reform on pharmaceutical prescription on the volume of prescriptions dispensed from pharmacies is analyzed [13].
The rest of the paper is organized as follows. Section 2 describes the segmented linear regression models and GAMs applied to ITSA. Section 3 describes the simulation exercise where the process followed to simulate the data and the comparison of the results for both models are shown. Section 4 describes the application to real data. Finally, Section 5 deals with the discussion and the conclusion of the paper.

The Interrupted Time Series Design
In the ITSA, we have an observed outcome variable Y t measured at each equallyspaced time point t. Y t is exposed to the intervention in periods from T 0 + 1 to T and unexposed in periods from 1 to T 0 . Suppose that Y t (1) and Y t (0) represent the outcome with and without the intervention, respectively, and Y t is given as follows: The intervention effect at time t is defined as the subtraction Y t (1) − Y t (0) for t = T 0 + 1, . . . , T, and the cumulative effect is defined as the sum of the intervention effects ∑ T t=T 0 +1 (Y t (1) − Y t (0)). However the counterfactual Y t (0) is never observed for the postintervention period, so the intervention effect is not observed in the data. To estimate the intervention effect, it is necessary to make an assumption regarding the outcomes that would have occurred in the absence of intervention. In a segmented regression analysis ( [5]), separate levels and trends are estimated in each segment before and after the intervetion.
Regression forecast for the post-intervention time period assuming the parameters of the pre-intervention period becomes an accurate option when counterfactual is not observed.

Segmented Linear Regression Models
A basic statistical method for ITSA is the segmented linear regression model. In a segmented linear regression, or a break-point model, each segment of the time series before and the intervention can have different levels and trends. The segmented linear regression allows the outcome of interest to evolve differently before and after the intervention. This approach controls for secular trends. A change in the level of the outcome after the intervention may constitute an abrupt intervention, and a change in trend shows a variation in the evolution of the series.
The specification of the linear regression is: where T t , or trend variable, is the value of the time variable at moment t (takes values from 0 to T − 1), I t is an indicator variable of the intervention that takes value 0 for the periods from 1 to T 0 and value 1 for the periods from T 0 + 1 to T, and the variable TI t is a sequential numbering of the time periods of the intervention, that takes value 0 for the pre-intervention period, and takes values from 0 to T − T 0 − 1 for the periods from T 0 + 1 to T. Following the definition in (1), the expressions Y t (0) and Y t (1) are: The linear regression analysis can accommodate additional structures that allow a more accurate estimate of the change in level and/or trend of the outcome due to the intervention, such as explanatory variables not affected by the intervention (or control variables), seasonality or serial correlation of the data. In this last case, when the errors are assumed to follow a first order autoregressive process AR(1), the linear regression model can be estimated using the Prais-Winsten method [14] which uses the generalized leastsquares method. When the order of correlation is assumed to be higher, the coefficients can be estimated using the OLS estimator but Newey-West standard errors [15] are used to handle this autocorrelation.

Generalized Additive Models
Generalized additive models (GAMs) are extensions of general linear models in which the outcome depends linearly on smooth functions of some predictor variables. GAMs link the outcome variable with the independent variables using smoothing splines, which are piecewise polynomials joined together at locations in the data known as knots. There are different methods proposed by the literature for smoothing with respect to a predictor variable (cubic regression splines, p-splines, adaptive smoothing, etc.) but the choice of the smoother has not been analyzed in this paper and we used the default thin plate regression splines. They are the default smoothing for the package mgcv in R [16] because they are an optimal smoother given basis dimension/rank ( [17]) and they are more flexible than the cubic smoothing splines [18]. Thin plate regression splines do not have knots but are more computationally costly to set up [16]. GAMs can be defined by the next equation: where f j is the smoothing spline for the independent variable x t . The more knots in the GAM, the more piecewise polynomials that will be estimated and the better the fit of the model to the data. The optimal set of smoothing functions is obtained by minimizing the penalized sum of squares criterion (PSS) in: where λ ≥ 0 is a parameter that controls the trade-off between the model's fidelity to the data and the function smoothness of f . A value of λ = 0 results in the relatively minimum smoothing, whereas large values result in an extremely smoothed (i.e., linear in the limit) function. The optimal smoothing parameter is chosen by cross-validation. GAMs can be useful for characterizing trends in longitudinal data when it is thought that change over time is non-linear but the exact nature of that nonlinearity is not known. In such a case, the independent variable x t would be a variable for time T t . This work proposes to apply GAMs for evaluating policy interventions. The expression in (2) would be replaced by the expression: where s 1 and s 2 are smoothing functions for each corresponding predictor. This model applies a smoothing function to both the secular trend (T t ) and to the sequential numbering of the time periods of the intervention TI t . It implies that there may be a nonlinear data trajectory without the implementation of the intervention, as well as a (potentially different) nonlinear data trajectory after the intervention. GAMs can also adjust for serial correlation of the data assuming a generalized additive mixed model. Smooths are specified as part of the fixed effects model formula, but the wiggly components of the smooth are treated as random effects. This approach allows correlated errors to be dealt with via random effects [19,20]. Control variables not affected by the intervention or seasonality can also be included in the equation assuming a linear or non-linear relationship with the outcome variable.

Data Generation Process
In this subsection we show how the simulated data were generated. For the sake of simplicity, we have assumed a fixed sample size of 100 for all time series, where the intervention affected the last 10 periods of the series. We have considered linear and non-linear trends before and after the intervention in the simulation process to study the performance of the segmented linear regression models and the GAMs in each of the cases. The number of simulations for each model was 500 and the possible autocorrelation in the series is considered assuming that the error term is distributed according to a first order autoregressive process AR(1) with parameter 0.3 and a standard deviation of the white noise process of 0.5.
The level change and the cumulative effect are analyzed and the performance of the models was evaluated through the comparison of estimated impacts of the intervention and the expected real impacts assumed in the simulation, using the mean squared error (MSE) and the mean percentage error (MPE). The first simulation model assumes a linear but different trend before and after the intervention. It also assumes a change in the level of 5 units. Table 1 shows the parameters of this model. Figure 1 shows the deterministic part of the time series, where it is easily observed the assumed impact of the intervention through the comparison of the series before and after the intervention. The expected level change is 5 and the expected cumulative effect which includes the change in trend is 50.9 (5 · 10 + 0.02 · (0 + 1 + 2 + . . . + 9)). Linear model Y t = 30 − 0.08 · T t + 5 · I t + 0.02 · TI t + u t 5 50.9 Quadratic model The second simulation model assumes a quadratic trend. Table 1 shows the parameters of this model and Figure 1 shows the deterministic part of the simulated data. In this case we have assumed a negative level change of −5 and an expected cumulative effect of −44.075.
Finally, the third simulation model assumes a cubic function. With this example we try to show a nonlinear model with several trend changes during the pre-intervention period. A cubic function can accommodate this behaviour. In practice a great majority of time series can be adequately fitted with polynomials with a maximum degree of 3 [21]. Table 1 and Figure 1 show the behaviour of this model, where the expected level change is assumed to be −5 and the expected cumulative effect is −59.5475.
The models have been estimated using the R statistical software, version 4.0.4. The codes are provided in the supplementary documents. Table 2 shows the results of the simulation analysis. Results include the mean and standard deviation of the level change estimated for the 500 simulations, along with the MSE and MPE obtained from the comparison with the expected level change for each simulation model. Similarly, the results for the cumulative effect are showed. As expected, the segmented linear regression model obtains the best results for the linear simulation model. The mean level change is very close to the real expected level change (5.0218 and 5, respectively). The MSE is 0.1550 and the MPE is 6.27%. The mean cumulative effect is 51.0629 versus the real expected cumulative effect of 50.9. The MPE for the cumulative effect is 3.98%. Surprisingly, the GAM also obtains good results for the linear simulation model. The mean level change is 5.0246 although the dispersion is greater than that observed for the segmented linear regression model (0.4411 and 0.3935, respectively). The MPE for the level change is slightly higher, 7.04%. The results are similar for the estimation of the cumulative effect, where the MPE is 4.49%.

Results of the Simulation Analysis
However, when the data are simulated from a non-linear model, the performance of the segmented linear regression model worsens. For the quadratic simulation model, the segmented linear regression model estimates a mean level change of −4.5442 when the real expected level change is −5. The MPE is 12.73%. The MPE for the cumulative effect is even greater, 46.43%. The flexibility of the GAM allows for a better fit for the non-linear simulation models. The mean estimated level change is −4.8969, very close to the real level change of −5. The good results of the GAM are maintained when estimating the cumulative effect. The mean cumulative effect is −40.5424 and the real cumulative effect is −44.075. The MPE is 11.19%.
For the more complex polynomial simulation model, the results of the GAM get worse, with an MPE of 11.33% and 16.98% for the estimation of the level change and the cumulative effect, respectively. However, these results are better than those obtained by the segmented linear regression model, where the MPEs are 15.94% and 39.54%, respectively. Figure 2 shows an illustrative example of how both models fit the same time series simulated from the polynomial model. The segmented linear regression model estimates a negative trend for the pre-intervention period. The level change is estimated at −4.3307, and the trend becomes a positive trend after the intervention, estimating a cumulative effect of −33.5915. The GAM fits a positive trend at the end of the pre-intervention period. The level change is now estimated at −4.6535. The estimated trend for the post-intervention period is lower than that observed before the intervention, and the cumulative effect is estimated at −45.1620.

Illustration with Real Data: Impact of the Cost-Sharing Reforms on Pharmaceutical Prescriptions Established in Spain
To illustrate the use of GAMs in a real-life application, we investigate the impact of the 2012 Spanish cost-sharing reforms on pharmaceutical prescription financed by the Spanish National Health Systems (SNHS) [13]. In June 2012 Spain enacted a reform of the co-payment for outpatient prescription drugs scheme, implemented early July 2012. Cost-free arrangements for all pensioners' drugs were replaced by a 10% co-payment subject to a monthly cap, depending on the income [22].
We use data relating to dispensed prescriptions for Pharmacies and financed by the SNHS. Data were collected from reports published by the Spanish Ministry of Health. We use monthly data from January 2004 to December 2015. The per-capita total prescription dispensed was calculated by dividing by the resident population of Spain.
The segmented linear regression model and GAM are applied to this dataset. Besides the level and trend before and after the intervention we have included as explanatory variables the seasonality (using indicator variables for the segmented linear regression model and a smoothing function for the GAM) and an indicator variable for the month previous to the intervention which examines the "stockpiling" effect between the announcement and the implementation of the law [23]. The codes are provided in the supplementary documents. Table 3 shows the descriptive statistics for the dependent variable for the pre-and post-intervention periods. The mean number of prescriptions decreased after the reform. The histogram plots of the per-capita prescriptions are shown in Figure 3. The Shapiro-Wilk normality tests confirm the normality hypothesis with p-values of 0.0806 and 0.1066 for the pre-and post-intervention periods, respectively.      Table 5 shows the results of the GAM. The model M1 includes three coefficients: the intercept, an indicator variable for the intervention (T t ) and the indicator variable for the month before the intervention which refers to the stockpiling effect. Besides, the model includes three terms to be smoothed, the secular trend s(T t ), the change in the trend after the intervention s(TI t ) and the seasonality s(month t ). For this last smooth term, the maximum possible dimensions of the basis used for the spline is set to 12 (k = 12), the number of months, while k is set by default as 9 for the rest of the smooth terms.
The level change is estimated in −0.2704 [IC95%: (−0.3096, −0.2312)], similar to that obtained by the segmented linear regression model, −0.3068. The stockpiling effect is statistically significant (p-value: 0.0009). The results for the smooth terms are summarized by the effective degrees of freedom (EDF), which measure the complexity of a penalized smooth term. EDF can be interpreted as an estimate of how many parameters are needed to represent the smooth [20]. If the EDF is equal to 1, a linear relationship cannot be rejected. In this analysis, the EDF is estimated at 4.167 for the secular trend showing its non-linearity. However, we cannot reject that the change in the trend after the intervention is linear. Seasonality is clearly non-linear. The cumulative effect is estimated at −9.5623 [IC95%: (−10.2504, −8.8742)], which is lower than that observed in the segmented linear regression model. Due to the linearity of the TI t covariate, an alternative model M2, where TI t is included in the linear part of the model, is also shown in Table 5. Estimates do not vary and the variable TI t is not statistically significant. The cumulative effect is also estimated in −9.5623 [IC95%: (−10.2504, −8.8742)]. The normality and unbiased error distribution were verified through four residual plots (using the command gam.check [16]). We also checked that the default maximum possible dimensions of the basis used for the trend spline (k = 9) was enough. Figure 4 shows the fitted values for both models. The fits are similar, but the GAM predicts a lower trend for the post-intervention period if the intervention had not been performed, which leads to a smaller total impact. Even in this case where the linear model fits the data well, the difference in the total impact estimated by both models is statistically significant.

Conclusions
Interrupted time series analysis (ITSA) is a useful quasi-experimental design with which to evaluate the longitudinal effects of interventions. Its design is particularly useful for evaluating policy interventions. Segmented linear regression models have been the most used models to carry out this analysis. However, it may not be appropriate when trends are not linear and they cannot be transformed to be so.
In this paper, we show how generalized additive models (GAMs) [24][25][26][27] can be useful to accommodate nonlinear trends. GAM is a non-parametric regression-based method that can estimate non-linear trends in time series and can handle the irregular structure of some time series.
Our method generalizes the widely used regression methods applied to ITSA, which explicitly models the time series observed both before and after the intervention. The projection of the pre-intervention model for the post-intervention period can then be used as a counterfactual for the post-intervention period as if the intervention had not occurred.
The analysis with simulated data showed how GAMs improve on existing methods when the trend is non-linear, but they also show good performance when the trend is linear. The analysis with simulated data also showed how the segmented linear regression model fails as the trend model gets more complex. Other intervention effects, such as a pulse intervention, or other non-linear models for the trends, are also possible but we do not expect to observe different conclusions.
A real-life application where the impact of the 2012 Spanish cost-sharing reforms on pharmaceutical prescription is analyzed allows us to observe the differences that can be achieved when applying one or another methodology even in the case where the linear model fits the data well. GAMs also allow the inclusion of other explanatory (control) variables into the analysis assuming a linear or non-linear relationship with the outcome. In our example, the seasonality was included assuming a non-linear trend. The EDF has shown how the change in trend after the intervention (TI t ) could be modeled linearly. In that case, we recommend its inclusion in the linear part of the model due to its simplicity of estimation and interpretation. The autocorrelation in the error term can also be considered with GAMs which makes this method flexible enough to be applied in most situations.
In addition to GAMs, there are other alternative statistical methods for dealing with non-linear trends. These methods can be divided into two categories: the first includes regression methods like GAM such as the autoregressive integrated moving average (ARIMA) [10], and local regression (LOESS) [28]. The second category includes computing methods such as recurrent neural networks (RNN) [29], and other artificial intelligent systems. ARIMA models are usually considered robust for a long time-series only. These models have been used more in a predictive than an explanatory approach. To use ARIMA models we must to transform a time series into stationary one. ARIMA models are backward looking and not very flexible. Besides, ARIMA models are subjective and the reliability of the chosen model depends on the skill and experience of the researcher. Like GAMs, LOESS is a non-parametric regression method that fits a smooth line through data. But unlike LOESS, GAMs use flexible smoothing functions with automatic choice of smoothing parameteres. Finally, opacity is the most important disadvantage of RNN methods. Furthermore, training of RNN models can be difficult [30]. GAMs allow for model shapes from linear to nonlinear trends, a balanced reducing of model uncertainty, and the identification of time-periods of significant events [31]. However, the propensity to overfit is the main limitation of GAMs. Another limitation is that the model will lose predictability when the smoothed variables have values outside of the range of training dataset. GAMs are also restricted to be additive, thus important interactions can be missed. However, as with regular linear regression, we can manually add the interaction effects. Funding: Financial support for this study was provided in part by Grant ECO2017-85577-P (Ministerio de Ciencia, Innovación y Universidades, Agencia Estatal de Investigación, Spain).