In recent years, time series studies with Poisson regression using generalised additive models (GAMs)^1–4 have been the reference method for analysing the short term health effects of air pollution. Lately, however, these models have been shown to suffer from an important limitation, namely, that there is underestimation of the standard error (SE) of the estimated effect of any given pollutant in those cases where concurvity is present in the data.^5–8

Briefly, concurvity is the non-parametric analogue of multicollinearity, in which a function of one of the independent variables can be approximated by a linear combination of functions of the remaining variables, with these functions being estimated in the same way as the corresponding functions in the original model.⁶ It has been seen that, in the presence of concurvity in GAMs, confidence intervals (CI) are too narrow, p values are understated, and type I error rate is greater than that established a priori^6,8 (for example, for a concurvity of 0.6—a value based on real data—type I error has shown a rise of almost threefold⁶). One way of eliminating the problem of underestimation of SE in the presence of concurvity is to have recourse to bootstrap resampling techniques. Accordingly, in this paper we propose the use of a conditional bootstrap method in GAMs to calculate valid CI for the estimated effect of any given pollutant. For the purpose of performing such calculations, we have developed a routine, termed gam.boot.

METHODS

Conditional bootstrap

In this type of bootstrap,^9–13 B bootstrap replicates are generated. In each of these, the values of the independent variables are the same as those of the observed data, with only the values of the response variable being varied from replicate to replicate. The value assumed by the outcome in each observation is conditional (hence the technique’s name) upon the values of the set of independent variables in said observation. To this end, a Poisson model is first applied to the data for the sample, and the coefficient and predictions of the model are then obtained for each observation (also known as pilot estimates, as these are taken as reference).

To illustrate this procedure, it is as well to start with a simple model (see expression 1) having only two covariates X = (X₁,X2), one of which, X₁, is linearly related to the response via parameter β₁ (the parameter of interest), and the other, X₂, is related to the response via an unknown smooth function, f. Based on n observations (X_1,Y₁),…, (X_n,Y_n), with X_i = (X_i1,X_i2), the following model is then fitted

and thereafter the estimated coefficient β̂₁ and predictions μ̂₁,…μ̂_n are obtained for each of the observations.

In a second step, the B conditional bootstrap replicates from b = 1 to B (for example, B = 1000) are generated, so that the values of the dependent variable, Y_i, in each observation follow a random Poisson distribution with a mean of μ̂_i, that is to say, Y_i^*b∼Poisson (μ̂_i). In each of these B replicates, an estimate of the coefficient β̂₁^*b is obtained, using a GAM with the same number of degrees of freedom as the original model.

Finally, the 100percent×(1−α) limits for the confidence interval of β₁ are given by (2β̂₁−β̂₁^1−α/2, 2β̂₁−β̂₁^α/2) where β̂^p represents the percentile p of the bootstrapped estimates β̂₁^*1,…,β̂₁^*B.

Bootstrap validation: a simulation study

Scenario

A three year (1096 day) series was constructed in which the number of events per day followed a Poisson distribution, in line with model (2) below,

where Y_t denotes the number of events per day, BS_t is the black smoke, and temp_t is the temperature recorded on day t. The coefficient β₁ denotes the log relative rate of Y associated with increase in black smoke (BS). The BS is genuine and was drawn from data recorded for the city of Vigo (north west Spain) over the period 1996–1998 (see fig 1A). By default, the true β₁ had a value of 0.001, although other values for this parameter (from 0.001 to 0.008) were also considered to ascertain the effect of coefficient magnitude on the 95% CI coverage and bias in the estimated coefficient, β̂₁.

In model (2), trend is a function that introduces trend and seasonality into the simulated data to simulate an unobserved confounding variable in the data. In ecological time series studies, smooth functions (smoothing splines or LOWESS) of the time variable ^1–4 are used to control for the confounding effect of these unobserved confounding variables. The trend function is a sinusoidal function (see fig 1B) that had already been used in other simulation studies.^14–16 In (2), as in such studies, β₂ assumes the value of 0.1. The function f_temp represents the functional form of the effect of temperature on mortality in Vigo (see fig 1C).

Concurvity generation

In our simulation study, the effect of interest was BS. The concurvity measure proposed by Ramsay et al⁶ was based on the correlation between the fitted values obtained from the GAM with pollution, BS_t, as the response, and time and temperature as the predictors. Specifically, to assess the degree of concurvity in our data, one should: (a) fit the GAM

where the partial functions g₁ and g₂ are estimated using smoothing splines with 7 and 4 degrees of freedom respectively, and ε_t is a zero mean error variable; and (b) then compute the squared correlation (R²) between BS_t and, using the fitted values from (3), namely:

In our real data, concurvity was 0.29. To vary the concurvity and thereby assess its influence on CI coverage and bias in the parameter estimate, β̂₁, in model (2), the original BS was replaced by a new variable

where ε̂_t = BS_t−▒BS_t, k₁ and k₂ were constants, such that 0⩽k_1,k₂⩽1 and 0⩽k₁+k₂⩽1. Note that BSs_t is a combination of the original BS_t of the estimate ▒BS_t and the errors ε̂_t = BS_t−▒BS_t obtained from the model (3). Thus, by varying the constants k₁ and k₂, the degree of concurvity between BSs and the remaining independent variables may easily be modified, from 0 (k₁ = k₂ = 0) through 1 (k₁ = 0, k₂ = 1).

Data analysis

In our study, a number of scenarios were defined using different values for the coefficient β₁ (from 0.001 to 0.009) and different values for k₁ and k₂, so that the degree of concurvity varied from 0.05 to 0.65. In each of the scenarios, 1000 samples {t,temp_t, Y_t} were generated with 1096 points in time (days) each (t = 1,…,1096) based on the model (2). In each sample, β̂₁ was obtained on the basis of the estimated model (4)

Lastly, the corresponding bootstrap and asymptotic 95% CI were calculated for β₁. The CI coverages were calculated as the proportion of samples in which these included the true β₁.

The estimates in (4) were obtained by using smoothing splines with 7 degrees of freedom per year for the estimated trend effect, f̂_trend (following Dominici et al),⁵ and 4 degrees of freedom for estimated temperature effect, f̂_temp.

RESULTS

Figure 1A shows the BS sequence for the city of Vigo, which was used as the independent variable to generate the replicates; fig 1B shows the unobserved confounding variable having trend and seasonality; fig 1C shows the functional form of the relation between temperature and mortality, which was obtained from previous estimates in a relation observed for Vigo.

In table 1, the CI coverages are compared for different degrees of concurvity using asymptotic GAM and GAM bootstrap. It will be seen that, with the asymptotic approach, increases in concurvity were accompanied by a progressive decline in coverage (see fig 2A), which decreased to as low as 80% when concurvity was 0.6. With the bootstrap method, however, coverage remained above 94% throughout, irrespective of the degree of concurvity. The coverage was also evaluated for different values of the true coefficient. With regard to the coefficient’s influence on coverage (see fig 2B), coverage proved seemingly independent of the magnitude of the coefficient, for the asymptotic and bootstrap methods alike.

DISCUSSION

The results of our study show that application of conditional bootstrap techniques could enable the CI for the effect of pollutants in time series studies to be calculated with optimal coverage, regardless of the degree of concurvity, number of covariates in the model or magnitude of the coefficient to be estimated.

The publication of two scientific papers a few months ago, the first by Dominici et al⁵ and the second by Ramsay et al,⁶ showed that standard errors are considerably underestimated in the presence of concurvity.^7,8 To eliminate the problem of concurvity, flexible regression methods can be used, such as parametric cubic splines⁶ or penalised splines, which do not require the use of backfitting to estimate model components. However, parametric cubic splines might not provide sufficient flexibility to capture the trend or seasonality of the series, as the knots have to be established a priori. In this respect, penalised splines^17,18 might afford greater flexibility, although, as far as we are aware, no assessment has yet been made of the performance of such smoothers in the context of time series studies.

To eradicate the problem of concurvity, consideration has been given to the application of alternative designs, such as that of case-crossover, which control trend and seasonality by means of design and thereby eliminate the effect of concurvity on the estimation of standard errors.¹⁶ Dominici et al¹⁹ have introduced a closed form estimate of the asymptotically exact covariance matrix of the linear component of a GAM, and offer a development of the gam.exact package, which is an extended version of gam.²⁰ However, it has also been observed that, in given conditions (with many smoothing functions included in the model), the routine gam.exact can give rise to standard errors very much greater than those that would result from application of parametric techniques.

In this paper, we propose an alternative method based on conditional bootstrap resampling techniques, which provides optimal coverages and affords the following additional advantages: (1) unlike the gam.exact method, which requires the smooth functions to be of the smoothing-spline type, gam.boot is general, in that it is applicable to models with any type of smoother, such as smoothing splines, LOWESS, natural splines, and any combination of same; (2) although our bootstrap method was initially developed to construct CI for coefficients, this approach is easily generalised to non-linear dose-risk relations. Thus, non-linear dose-response relations can be evaluated for pollutants or climatological variables; and lastly, (3) our method could also be extended to the calculation of CI for possible interactions between pollutants and/or climatological variables.

Furthermore, a recent study seems to show that, in cases where the results for a number of cities are combined using hierarchical models, underestimation of standard errors in each of the cities seems to have little effect on the overall results under most conditions.²¹ In this regard, our method could be generalised to a hierarchical bootstrap method that would start on the basis of standard errors correctly calculated at the first level (city), and then proceed to estimate the overall effects at a multicity level.

While the main limitation of our method might possibly lie in the computational cost entailed, we do not regard this as too high. Furthermore, there is no need to apply the bootstrap method at each stage of building the model: suffice to apply it once the basal model has been constructed and the final effect of each pollutant is to be ascertained.

Indeed, the conditional bootstrap method proposed in this paper could enable the community of air pollution researchers to make a satisfactory calculation of CI in dose-response relations, whether linear or non-linear. To this end, the S-plus gam.boot software, a user friendly application, will be made freely available to readers on request (rocauvigo.es).