Email updates

Keep up to date with the latest news and content from BMC Medical Research Methodology and BioMed Central.

Open Access Research article

Using a generalized additive model with autoregressive terms to study the effects of daily temperature on mortality

Lei Yang1, Guoyou Qin1, Naiqing Zhao1*, Chunfang Wang2 and Guixiang Song2

Author Affiliations

1 Department of Biostatistics, School of Public Health, Fudan University, Shanghai, China

2 Vital Statistics Division, Shanghai Municipal Center for Disease Control & Prevention, Shanghai, China

For all author emails, please log on.

BMC Medical Research Methodology 2012, 12:165  doi:10.1186/1471-2288-12-165

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2288/12/165


Received:30 January 2012
Accepted:9 October 2012
Published:30 October 2012

© 2012 Yang et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided original work is properly cited.

Abstract

Background

Generalized Additive Model (GAM) provides a flexible and effective technique for modelling nonlinear time-series in studies of the health effects of environmental factors. However, GAM assumes that errors are mutually independent, while time series can be correlated in adjacent time points. Here, a GAM with Autoregressive terms (GAMAR) is introduced to fill this gap.

Methods

Parameters in GAMAR are estimated by maximum partial likelihood using modified Newton’s method, and the difference between GAM and GAMAR is demonstrated using two simulation studies and a real data example. GAMM is also compared to GAMAR in simulation study 1.

Results

In the simulation studies, the bias of the mean estimates from GAM and GAMAR are similar but GAMAR has better coverage and smaller relative error. While the results from GAMM are similar to GAMAR, the estimation procedure of GAMM is much slower than GAMAR. In the case study, the Pearson residuals from the GAM are correlated, while those from GAMAR are quite close to white noise. In addition, the estimates of the temperature effects are different between GAM and GAMAR.

Conclusions

GAMAR incorporates both explanatory variables and AR terms so it can quantify the nonlinear impact of environmental factors on health outcome as well as the serial correlation between the observations. It can be a useful tool in environmental epidemiological studies.

Background

People have long been interested in potential impacts of temperature on human health [1-3]. Recently, many studies have been done to analyse the way and the extent temperature influences health outcomes [4-6]. Such studies have practical value for the following reasons: knowing the association between environmental factors and health outcomes will help to identify at-risk populations, benefit health department in resource allocation, and provide support for stake holders in prevention [7].

Environmental epidemiology, the investigation of the health risks related to environmental exposures, becomes the main research approach [8]. There are various study designs in environmental epidemiology for estimating health effects of temperature. One of them is model evaluation of time series data to quantify the association between temperature (or other weather/environmental factors etc.) and daily mortality (or hospitalization etc.). These are a type of ecologic study because they analyse daily population-averaged health outcomes and exposure levels [9]. This approach is suitable to study transient acute effects which are due to time-varying exposures [10-14], and the choice of model can have a large influence on the interpretation of environmental effects.

Generalized Linear Model (GLM) [15] and Generalized Additive Model (GAM) [16] are the main models used in environmental epidemiology. When the response variable represents counts (e.g. number of deaths), the model often takes the form of a Poisson regression/additive model with a log link function. The outcome Yt is assumed to follow a Poisson distribution with mean μt which is linked either to a linear combination (in this case, the model is a GLM) or smoothed functions (in this case, the model is a GAM) of environmental explanatory variables via a log function.

Many studies have identified the nonlinear relations (U, J, or V shaped) of daily mortality on temperature [17-24], which are mainly represented by piecewise linear terms [18-20] or natural cubic splines [21-24]. An assumption is often made that temperature only affects the mortality of the same day or on a single day of lag l [18,20-22], and denoted by the term tempt-l in the model. When l=0, it represents the effect on the same day [18,20-23], when l>0, it represents a lagged effect [18,20,21,23]. An alternative to single lag models is a distributed lag model [25]. It includes multiple lags of temperature over a period of time based on the assumption that the effect of time is distributed over several days into the future. In all the above models, a spline function of time is often used to explain the long time trend, and confounding effects like humidity and air pollution are also controlled by splines.

However, a thorough time series analysis should consider the order of data points and correlation of adjacent points in time [26]. In environmental epidemiological studies, the response variable may also be correlated and it is necessary to embody autocorrelation of the response variable when modelling. But all the above models are standard GLM/GAM that describes how the response variable is stochastically related to explanatory variables without considering how the response can be dependent also on its past values.

In addition, autocorrelation causes trouble in estimation of GLM/GAM, since GLM/GAM essentially requires each observation to be independently distributed. Violation of this assumption can lead to problematic estimates even in very simple cases. For example, if the error terms in a linear regression model are in reality positively autocorrelated, failure to account for this may underestimate the standard errors of the estimated regression coefficients [27].

Another statistical issue often encountered is overdispersion, and many sources of autocorrelation are related to sources of overdispersion [28]. One usual approach to adjust for overdispersion is to specify a dispersion parameter ϕ for estimation [15]. However, if autocorrelation exists, this approach only inflates the variance of the estimate by ϕ but leaves the estimate unchanged, which is an inadequate solution to our issue.

Generalized Estimating Equations (GEE) [29,30] and Generalized Linear Mixed Model (GLMM)/Generalized Additive Mixed Model (GAMM) [31,32], are two extensions from GLM/GAM for grouped or clustered data. GEE and GLMM/GAMM can be implemented in popular software like SAS and R [31-34]. For all these models, a within-group variance-covariance structure can be used to account for the corresponding within-group autocorrelation [31,34]. Time series data, treated as single cluster data, can also be modelled by GAMM. Performance of this degenerative GAMM will be studied in simulation.

In this article, we introduced GAM with Autoregressive terms (GAMAR), which is derived from Generalized Autoregressive Moving Average (GARMA) models [35]. ARMA is a family of models for analysing time series. The notation ARMA(p,q) refers to a model with p autoregressive terms and q moving-average terms [26]. GARMA belongs to the class observation driven models [36] for extending Gaussian ARMA to non-Gaussian settings. GAMAR differs from GARMA in that the linear components in GARMA are generalized to natural splines and MA terms are omitted. The generalization is motivated by modelling nonlinear relationships, while the simplification is justified because AR can be used to approximate MA or ARMA. Additionally, larger AR order in GAMAR will not compromise the estimation of the effects of explanatory variables.

In all, GAMAR has two advantages over GAM: 1) it is a model for generalized time series analysis rather than a probabilistic model like GAM; 2) the AR part of GAMAR can explain the autocorrelation structure of observations. So the Pearson residuals of GAMAR can be closer to white noise than those from GAM, yielding more reliable estimation.

Methods

GARMA

In GARMA, the conditional distribution of each observation yt, for t=1,…,n, given the previous information set Ht={X1,…,Xt,y1,…,yt-1} containing past observations yt and covariate vectors Xt = (Xt1,…,Xtm), is assumed to follow the same exponential family distribution [35]. As with the standard GLM, the conditional mean μt is related to the variables by a twice-differentiable one-to-one monotonic function g, which is called the link function. However, unlike the standard GLM, the formula here allows autoregressive moving average terms to be included additively [35]:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M1">View MathML</a>

where <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M2">View MathML</a> are autoregressive terms and <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M3">View MathML</a> are moving average terms.

Count data are often assumed to follow Poisson distribution. And the Poisson GARMA submodel is:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M4">View MathML</a>

where yt* = max(yt), τ is a positive threshold parameter. Any 0 or negative values of y are replaced by τ, because ln() is not defined for 0 or negative values.

GAMAR

GAMAR is derived from GARMA with linear terms replaced by smoothers and MA terms omitted:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M5">View MathML</a>

(1)

where <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M6">View MathML</a> are smoothers of covariates, <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M7">View MathML</a> are autoregressive terms.

Compared to GAM, (1) allows autoregressive terms to be included additively in the link predictor.

For count data y, we use the Poisson submodel:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M8">View MathML</a>

(2)

here yt* = max(yt), τ is a positive threshold parameter.

In our study, we use natural cubic spline (ns) as smoother. ns is a piecewise-cubic real function on an interval a,b, where a,b is separated by a sequence of k ordered knots: α = ξ0< ξ1< ⋯ < ξk-1< ξk= b. ns is continuous at interior knots and linear beyond the boundary [16]. Degrees of freedom (df) for ns equals to the number of subintervals separated by these knots, thus df satisfies: df = k. ns is often represented by linear combination of its B-spline basis [16]. When df is specified, the standard practice is to place df-1 interior knots at evenly spaced intervals in the data to generate the B-spline basis. By using ns as smoother in (2), the Poisson GAMAR becomes:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M9">View MathML</a>

Algorithm

Maximum partial likelihood estimator (MPLE)

For the jointly distributed time series {Xt,yt}, t = 1,…n, the parameters of GAMAR can be estimated by maximum partial likelihood. The partial likelihood based on yt for {Xt,yt}, t = 1,…n can be expressed as the product of a sequence of conditional likelihoods f(yt| X(t),y(t-1)), t = 1,…n, where X(t)= X1,…Xt, y(t-1) = y1,…yt-1, that is:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M10">View MathML</a>

This function is referred to as a partial likelihood instead of likelihood because Xt are stochastic [35]. For a Poisson GAMAR, the partial likelihood is:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M11">View MathML</a>

and μt can be expanded as <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M12">View MathML</a>. Since a natural cubic spline is a linear combination of its B-spline basis, it is linear with respect to its parameters, so natural cubic spline functions are like linear terms from the perspective of computation. Therefore,

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M13">View MathML</a>

(3)

where θ = (β1m,c1,⋯,cp)T is the model parameter vector.

Modified Newton’s method

Maximum partial likelihood estimators are solved by modified Newton’s method. This procedure is described in detail in Appendix.

Evaluation of GAMAR by simulation

We conducted two simulation studies, each with 1000 samples, to compare the performance of GAM and GAMAR for estimation. In simulation study 1, the performance of GAMM is also studied. The R code for simulations and simulation output is included in Additional file 1.

Additional file 1. R code for this study. This file contains core R codes for this study, including the function of GAMAR, data generation in simulation studies, data fitting in simulation studies, real case analysis (including the procedure to choose the parameters) and generation of tables and figures. This file also contains a brief description of every R program.

Format: ZIP Size: 28KB Download fileOpen Data

Simulation study 1

The first model:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M14">View MathML</a>

(4)

Here xt represents a daily averaged temperature series from year 2000–2008 obtained from Shanghai's Meteorological Bureau, and the terms si5(xt), i = 1, 2, 3, 4, 5 form the B-spline basis for the natural cubic spline. The parameters are:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M15">View MathML</a>

We used data ranging over time points 1828–3288 (year 2005–2008) to eliminate the impact of the starting value, where the time points stand for the days since Jan 1st 2000.

For two models on Sample 1, the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) [26] of the Pearson residuals [15] are plotted against different lag periods to examine the presence of residual autocorrelation. Estimates of temperature effects from the two models and the true effect were plotted against temperature to show which model provided a better fit.

Besides showing the analysis for a single sample, the averaged results from GAM, GAMM and GAMAR were given, and the statistics calculated were mean estimates, bias, relative error and coverage, short for “coverage rate of confidence interval (CI) on true value”. The confidence level was chosen to be 95%, so the coverage of a correct model should be around 95%.

In addition, GAMAR with various AR orders were studied to explore the effects of AR order on estimation. Finally, to verify the consistency of partial maximum likelihood estimation, 1000 samples in the time range of 2558–3288 (2 years), 1828–3288 (4 years), 1097–3288 (6 years) were used for calculation.

Simulation study 2

A second simulation study was performed to study whether the proposed method could approximate a nonlinear curve. We simulated 1000 samples from a model where the covariate is a cosine function of temperature and nonlinear with respect to the parameters.

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M16">View MathML</a>

(5)

Here, the AR terms parameters are identical to (4): c1 = 0.5, c2 = 0.25, c3 = 0.12.

As in simulation study 1, sample statistics were calculated over the time points 1828–3288 (year 2005–2008) for each of the 1000 realizations by GAM and GAMAR. The first sample was analysed just the same way as in simulation study 1, but here the mean and standard deviation of the parameter estimates from the two models for the 1000 samples are presented. To compensate for the absence of true parameters, we compared the two models visually, plotting the mean predicted values from the two models and the true effect. Moreover, the Pearson correlation coefficients of the fitted linear predictors from the two models and true effect were also calculated.

Application to temperature-mortality research

The real example came from three sources: daily mortality of Shanghai in 2001–2004 from Shanghai Municipal Center for Disease Control & Prevention; daily average temperature and humidity in 2001–2004 from Shanghai's Meteorological Bureau; and air pollution data (NO2, SO2, PM10) in 2001–2004 from (http://www.envir.gov.cn/airnews/ webcite), which was announced by the Shanghai Environmental Monitoring Centre.

Autocorrelation of the Pearson residuals from the two models are compared via figure. We also compared their parameter estimates, and their estimated temperature effects via figures. The R code for real example modelling and its outputs is included in Additional file 1.

Results

Simulation results

Simulation study 1

The ACF and PACF plot of the GAM Pearson residuals from the first sample indicate obvious autocorrelation (Figure 1). The Pearson estimate of the dispersion parameter [37] is <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M17">View MathML</a>, indicating data overdispersion. Since the ACF tails off and PACF cuts off after lag 3, it is natural to use GAMAR(3) instead. The model is given below:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M18">View MathML</a>

(6)

thumbnailFigure 1. ACF and PACF of GAM and GAMAR(3) for Sample 1 from simulation study 1. These statistics are calculated on Pearson residuals <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M19">View MathML</a> from both models.

The ACF and PACF plot of the GAMAR(3) Pearson residuals suggest no autocorrelation, and the dispersion parameter estimate is now <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M20">View MathML</a>, indicating no overdispersion. Apparently, GAMAR(3) controls autocorrelation and overdispersion simultaneously, and Figure 2 indicates that the predicted spline function <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M21">View MathML</a> from GAMAR(3) is much closer to the real model than that from GAM.

thumbnailFigure 2. The temperature effects in link scale for Sample 1 from simulation study 1. Black: the true effect, Red: <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M22">View MathML</a> from GAM, Blue: <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M23">View MathML</a> from GAMAR(3).

Concerning the general performance of GAM and GAMAR, Table 1 shows that the biases of the mean parameter estimates from GAM are almost the same as GAMAR. However, the coverages of the 95% confidence intervals from GAM are far less than 95%, while those from GAMAR are very close to 95%. Also, relative errors of the parameter estimates from GAM are larger than that of GAMAR.

Table 1. Results from GAM, GAMM, GAMAR(3) in simulation study 1

Coverages for GAM corrected for overdispersion [15] are also shown in Table 1 as coverage2. This approach is the same as a standard GAM except that its standard error has been expanded by the dispersion parameter ϕ. As a result, the bias and relative errors remain unchanged, while the confidence interval is broadened, coverage improved. However, the coverage is still far less than 95%.

The results of GAMM with AR(3) correlation structure and GAMAR(3) are very similar, but the former model is much more time consuming. This disparity will be described in discussion. Hence results from only 50 samples fitted by GAMM are summarized in Table 1.

To examine how the AR order in GAMAR influences the results, we also fitted the GAMAR(1), GAMAR(2), GAMAR(4), and GAMAR(5). As Table 1 in Additional file 2 shows, the coverage grows and relative error drops as the AR order increases from 1 to 3. For models with AR order larger than 3, there is little difference in coverage, relative error and bias among them.

Additional file 2. Additional Tables. This file contains 2 tables related to simulation study 1.

Format: DOCX Size: 94KB Download fileOpen Data

The results for different time ranges are listed in Table 2 in Additional file 2, from which we can see that the larger the number of time-points, the lower the bias and relative errors, and the higher the coverage. This illustrates the asymptotic unbiased (consistency) property of MPLE.

Table 2. Results from GAM and GAMAR(3) in simulation study 2

Simulation study 2

In this simulation study, The ACF and PACF plot of the GAM Pearson residuals from the first sample also reveal obvious autocorrelation (Figure 3). Just as in simulation study 1, we can see that ACF tails off and PACF cuts off after lag 3 for GAM, suggesting AR(3) would be suitable. In contrast, the ACF and PACF of GAMAR(3) on the same data are both very close to 0. And Figure 4 indicates that the predicted spline function <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M34">View MathML</a> from GAMAR(3) is much closer to the real model than that from GAM.

thumbnailFigure 3. ACF and PACF of GAM and GAMAR(3) for Sample 1 from simulation study 2. These statistics are calculated on Pearson residuals <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M35">View MathML</a> from both models.

thumbnailFigure 4. The temperature effects in link scale for Sample 1 from simulation study 2. Black: the true effect, Red: <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M36">View MathML</a> from GAM, Blue: <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M37">View MathML</a> from GAMAR(3).

In Figure 5, we can see the mean estimated temperature effects from GAMAR(3) is closer to the true effect than that from GAM. Meanwhile, the Pearson correlation coefficients between estimated temperature effect and true effect are also different (GAM: 0.9341; GAMAR(3): 0.9980), which means GAMAR(3) provides a better fit. Also, the standard deviations of estimates from GAM are larger than those from GAMAR (Table 2).

thumbnailFigure 5. The averaged temperature effects in link scale from Simulation study 2. Black: the true effect, Red: <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M38">View MathML</a> from GAM, Blue: <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M39">View MathML</a> from GAMAR(3).

Application to temperature-mortality research

In the real case analysis, GAM is first used. The long time trend of original observations is unobvious as shown in the left side of Figure 6. So a rather small df: 2, is used to control this secular trend. The daily mortality after adjustment is in the right side of Figure 6.

thumbnailFigure 6. Mortality before and after adjustment for secular trend. Left: the original mortality, Right: the mortality adjusted for secular trend.

And the complete GAM is given below:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M40">View MathML</a>

(7)

Lagged days of two temperature terms and df of all remaining natural spline functions are first determined in a sequence to minimize AIC. Then this set of parameters is used as starting value to find the final parameters which minimize AIC locally. The final parameters are: lag1 = 4, lag2 = 10, dftempt − lag1 = 5, dftempt − lag2 = 4, dfpres = 2, dfhumi = 2, dfno2 = 3, dfso2 = 2, dfpm10 = 3

The two lagged temperature terms separately represent short term temperature effect and long term temperature effect. The term <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M41">View MathML</a> stands for week effect, where weekt represents corresponding day in a week for date t, and Ii (weekt), i = 1,2,3,4,5,6 are indicator functions for the number of a day in a week. For each function, if weekt = i, Ii (weekt) = 1; if weekti, Ii (weekt) = 0.

Figure 7 shows that PACF all exceed 95% CI bounds for the autocorrelations (blue dashed line) [38] lags less than 5, and are contained within the bounds for lags larger than 4. So GAMAR(4) is then chosen to fit the data as described below:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M42">View MathML</a>

(8)

thumbnailFigure 7. ACF and PACF of GAM and GAMAR(4) for real case. These statistics are calculated on Pearson residuals <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M43">View MathML</a> from both models.

All ACF and PACF of Pearson residuals from (8) are now below 0.1, thus Pearson residuals from GAMAR(4) are close to white noise. The estimated coefficients of two lagged temperature terms and AR terms are given in Table 3. In Figure 8 and Figure 9, we can see that the estimated effects of tempt-4, tempt-10 from two models are different.

Table 3. temperature effects and AR estimates from GAM and GAMAR(4)

thumbnailFigure 8. Effects of tempt − 4. Black: <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M44">View MathML</a> from GAM, Red: <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M45">View MathML</a> from GAMAR(4).

thumbnailFigure 9. Effects of tempt − 10. Black: <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M46">View MathML</a> from GAM, Red: <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M47">View MathML</a> from GAMAR(4).

Discussion

In environmental epidemiological studies, the meteorological/environmental influences on human health indicators are often explained in a modelling framework. And the predominately used models are GLM/GAM. In our research, the drawbacks of GLM/GAM are studied, and their melioration: GAMAR, is given. Simulation studies reveal that GAMAR is more suitable when observations are autocorrelated.

While the true AR order is known in simulation studies, AR order needs to be determined in real case. Three enlightening perspectives into this issue are given below:

1. In time series analysis, ACF and PACF are good indicators of the order of process. CIs of an uncorrelated series are often treated as criterion for selecting orders [38]. For example, when PACF exceed the limit of CI at certain lag, autoregressive term at that lag needs to be modelled. This approach is used to preliminarily determine the AR order in the real case study.

2. Significance of AR terms can also be considered in modelling. If p-value of AR(n0+1) from GAMAR(n0+1) is larger than 0.05, while p-value of AR(n0) from GAMAR(n0) is less than 0.05, then GAMAR(n0+1) is unnecessary and GAMAR(n0) is a favorable model. To illustrate it, we also use GAMAR(5) to the real data, and find p-value of AR(5) to be 0.1460>0.05.

3. The goal of introducing AR terms is to control autocorrelation. So if Pearson residuals of GAMAR(n0) show little sign of autocorrelation, or the PACF is within the CI for all lags, then the AR terms are adequate. In real data case, this requirement is justified.

In application, the first and second advices are practical methods to determine AR order preliminarily, and the third can be used to finally justify the choice. There have been extensive discussions about choosing AR orders [39].

We also want the covariate parameter estimation to be robust with respect to different AR orders. Simulation study provides some information for this issue: while the real model is GAMAR(3) in simulation study 1. GAMAR(1), GAMAR(2), GAMAR(4), GAMAR(5) are also used for estimation. From Table 1 in Additional file we can see that GAMAR(1), GAMAR(2) fit data much better than GAM; while they differ little from GAMAR(3). Since the true AR(3) parameter is rather small (0.12), probably when the effect of neglected AR term is small, covariate parameter estimation won’t differ much. When the AR order is larger than the real, the estimation results are almost the same as GAMAR(3). In all, covariate parameter estimation is robust with respect to disturbance of AR order.

GAMM can also give good estimates as shown in Table 1. However, while the speed of GAMM in fitting short time series is acceptable, say 2.05 seconds for sample with 100 observations and 10.17 seconds for sample with 200 observations (both executed in server), the time grows nonlinearly with respect to the length of time series and thus GAMM would become computationally formidable when time series data is long. For example, one sample of 1461 observations in simulation study 1 would take about 1 hour and 7 minutes for GAMM, while only 0.02 seconds for GAMAR (both executed in server): the former is 2×105 times the latter. Thus GAMM is much more computing intensive than GAMAR when sample size is large. Therefore, we only calculated 50 samples for GAMM in simulation study 1.

While we use natural cubic spline as smoother in the model, penalized spline and smoothing spline can also be implemented. With natural splines, one constructs a spline basis with knots at fixed locations throughout the range of the data. Smoothing splines and penalized splines have circumvented the problem of choosing the knot locations by constructing a very large spline basis and then penalizing the spline coefficients to reduce the effective number of df [40]. Despite the flexibility the penalized way provided, [40] identified both fully parametric and nonparametric methods can perform well in similar studies. Thus natural spline smoother can still address many practical problems. In addition, our model can be straightforwardly extended to the nonparametric way by including the penalty.

For natural spline with B-spline basis, selecting the df is of essential importance for application. A general approach is to use a data-driven method and to select the number of df which optimizes a particular criterion [40], like AIC. We used AIC to determine df of all splines except that of time in the real case. Another strategy is to use a df based on background knowledge or previous studies. For example, natural spline of time is chosen to represent the long term trend among different years, and itself shouldn’t contain any yearly fluctuation. Whether df meets the requirement is judged by comparing the trend of observations before and after adjustment, as well as the shape of the spline function visually in Figure 6. Section 2.1 of [40] gives a full treatment of this issue.

Besides including lagged temperature effects additively in the model, [41,42] have developed a family of distributed lag non-linear models (DLNM), which can simultaneously represent non-linear exposure-response dependencies and delayed effects. Still, AR terms can be aptly added to DLNM just as what we’ve done to GAM.

Finally, the model can also be applied in other areas. Researchers can use GAMAR for their specific research purposes. The most immediate extensions are other environmental epidemiological studies, specifically studies quantifying relationships between air pollution and mortality. Air pollution functions in a different way from temperature to impact human health, and there are also some differences between their modelling: 1) the air pollution-mortality relation is often simplified as linear, for such simplicity facilitates parameter estimation and interpretation; 2) distributed lag model is often used, since the effects of air pollution can last for many days; 3) cumulative effect: the total impact of pollution of a certain day over a period of following days, represented by the sum of parameters at all lags, is of great interest [9]. Such differences pose new questions for further study, like assessing the impact of AR terms on estimated cumulative effect rather than a single parameter.

Conclusions

This article proposes GAMAR for fitting time series data with explanatory variables and autoregressive terms. Two simulation studies with functions that approximate the response from the real example showed that GAMAR performed better than GAM. In the real example, there was residual autocorrelation with GAM, but little sign of autocorrelation with GAMAR. Also, different estimates of the temperature effects were obtained with GAMAR and GAM.

Appendix

Modified Newton’s method

Given partial likelihood, maximum partial likelihood estimators are solved by a modified Newton’s method. For Newton’s method, the iteration goes:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M48">View MathML</a>

(9)

until convergence.For a modified Newton’s method, the iteration goes:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M49">View MathML</a>

until it convergence.Here <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M50">View MathML</a>, which is the information matrix, and Γ*− 1(θm) is a modified version of Γ− 1(θm).In (9):

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M51">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M52">View MathML</a>

Since <a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M53">View MathML</a>, then:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M54">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M55">View MathML</a>

So

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M56">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M57">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M58">View MathML</a>

And the modified inverse matrix Γ*− 1(θm) is defined as: if Γ(θ) is reversible, then Γ*− 1(θ) = Γ− 1(θ);If Γ(θ) is irreversible. And its eigenvalue are λ1, λ2,…,λmp, then we can find orthogonal matrix P, which satisfies:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M59">View MathML</a>

Let λi* = max(λi, δ), δ = 0.01, i = 1, 2, ⋯, mp, then:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/12/165/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/12/165/mathml/M60">View MathML</a>

Such procedure ensures Γ*− 1(θm) to be positive definite.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

LY conducted the analysis and drafted the manuscript. GQ participated in discussions and provided advice. NZ inspired and tutored LY, and reviewed the paper. CW and GS collected the mortality data for real case study. All authors read and approved the final manuscript.

Acknowledgements

This research is supported by a project 30972551 from National Natural Science Foundation of China (NSFC).

References

  1. Ellis FP, Nelson F: Mortality in the elderly in a heat wave in New York City, August 1975.

    Environ Res 1978, 15:504-512. PubMed Abstract | Publisher Full Text OpenURL

  2. McKee CM: Deaths in winter: can Britain learn from Europe?

    Eur J Epidemiol 1989, 5(2):178-182. PubMed Abstract | Publisher Full Text OpenURL

  3. Ballester F, Corella D, Perez-Hoyos S, Saez M, Hervas A: Mortality as a function of temperature. A study in Valencia, Spain, 1991–1993.

    Int J Epidemiol 1997, 26(3):551-561. PubMed Abstract | Publisher Full Text OpenURL

  4. Ha J, Shin Y, Kim H: Distributed lag effects in the relationship between temperature and mortality in three major cities in South Korea.

    Sci Total Environ 2011, 409(18):3274-3280. PubMed Abstract | Publisher Full Text OpenURL

  5. Yu W, Guo Y, Ye X, Wang X, Huang C, Pan X, Tong S: The effect of various temperature indicators on different mortality categories in a subtropical city of Brisbane, Australia.

    Sci Total Environ 2011, 409(18):3431-3437. PubMed Abstract | Publisher Full Text OpenURL

  6. Basu R, Malig B: High ambient temperature and mortality in California: Exploring the roles of age, disease, and mortality displacement.

    Environ Res 2011, 111(8):1286-1292. PubMed Abstract | Publisher Full Text OpenURL

  7. Patz JA, Engelberg D, Last J: The effects of changing weather on public health.

    Annu Rev Publ Health 2000, 21:271-307. Publisher Full Text OpenURL

  8. Baker DANM: Environmental Epidemiology: Study Methods and Application. Oxford: Oxford University Press; 2008. OpenURL

  9. Peng RD, Dominici F: Statistical methods for environmental epidemiology in R. New York: Springer; 2008. OpenURL

  10. Ostro BD, Roth LA, Green RS, Basu R: Estimating the mortality effect of the July 2006 California heat wave.

    Environ Res 2009, 109(5):614-619. PubMed Abstract | Publisher Full Text OpenURL

  11. Hertel S, Le Tertre A, Jockel KH, Hoffmann B: Quantification of the heat wave effect on cause-specific mortality in Essen, Germany.

    Eur J Epidemiol 2009, 24(8):407-414. PubMed Abstract | Publisher Full Text OpenURL

  12. Goldberg MS, Gasparrini A, Armstrong B, Valois MF: The short-term influence of temperature on daily mortality in the temperate climate of Montreal, Canada.

    Environ Res 2011, 111(6):853-860. PubMed Abstract | Publisher Full Text OpenURL

  13. Hajat S, Kovats RS, Atkinson RW, Haines A: Impact of hot temperatures on death in London: a time series approach.

    J Epidemiol Community Health 2002, 56(5):367-372. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  14. Ballester J, Robine JM, Herrmann FR, Rodo X: Long-term projections and acclimatization scenarios of temperature-related mortality in Europe.

    Nat Commun 2011, 2:358. PubMed Abstract | Publisher Full Text OpenURL

  15. McCullagh P, Nelder JA: Generalized Linear Models. London: Chapman and Hall/CRC; 1989. OpenURL

  16. Hastie TJ, Tibshirani RJ: Generalized Additive Models. London: Chapman and Hall/CRC; 1990. OpenURL

  17. Braga AL, Zanobetti A, Schwartz J: The effect of weather on respiratory and cardiovascular deaths in 12 U.S. cities.

    Environ Health Perspect 2002, 110(9):859-863. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  18. Carson C, Hajat S, Armstrong B, Wilkinson P: Declining vulnerability to temperature-related mortality in London over the 20th century.

    Am J Epidemiol 2006, 164(1):77-84. PubMed Abstract | Publisher Full Text OpenURL

  19. Muggeo VM: Modeling temperature effects on mortality: multiple segmented relationships with common break points.

    Biostatistics 2008, 9(4):613-620. PubMed Abstract | Publisher Full Text OpenURL

  20. Kim H, Ha JS, Park J: High temperature, heat index, and mortality in 6 major cities in South Korea.

    Arch Environ Occup Health 2006, 61(6):265-270. PubMed Abstract | Publisher Full Text OpenURL

  21. Curriero FC, Heiner KS, Samet JM, Zeger SL, Strug L, Patz JA: Temperature and mortality in 11 cities of the eastern United States.

    Am J Epidemiol 2002, 155(1):80-87. PubMed Abstract | Publisher Full Text OpenURL

  22. Ren C, Williams GM, Morawska L, Mengersen K, Tong S: Ozone modifies associations between temperature and cardiovascular mortality: analysis of the NMMAPS data.

    Occup Environ Med 2008, 65(4):255-260. PubMed Abstract | Publisher Full Text OpenURL

  23. Hoffmann B, Hertel S, Boes T, Weiland D, Jockel KH: Increased cause-specific mortality associated with 2003 heat wave in Essen, Germany.

    J Toxicol Environ Health A 2008, 71(11–12):759-765. PubMed Abstract | Publisher Full Text OpenURL

  24. Anderson BG, Bell ML: Weather-related mortality: how heat, cold, and heat waves affect mortality in the United States.

    Epidemiology 2009, 20(2):205-213. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  25. Schwartz J: The distributed lag between air pollution and daily deaths.

    Epidemiology 2000, 11(3):320-326. PubMed Abstract | Publisher Full Text OpenURL

  26. Hamilton JD: time series analysis. Princeton, NJ: Princeton University Press; 1994. OpenURL

  27. Kutner M, Nachtsheim C, Neter J, Li W: Applied Linear Statistical Models Fifth Edition. New York: McGraw-Hill/Irwin; 2005. OpenURL

  28. Barron DN: The analysis of count data: overdispersion and autocorrelation.

    Sociol Methodol 1992, 22:179-220. OpenURL

  29. Zeger SL, Liang KY: Longitudinal data analysis for discrete and continuous outcomes.

    Biometrics 1986, 42(1):121-130. PubMed Abstract | Publisher Full Text OpenURL

  30. Liang KY, Zeger SL: Longitudinal data analysis using generalized linear models.

    Biometrika 1986, 73(1):13-22. Publisher Full Text OpenURL

  31. Pinheiro JC, Bates DM: Mixed-Effects Models in S and S-plus. New York: Springer; 2000. OpenURL

  32. Wood SN: Generalized Additive Model: an introduction with R. New York: Chapman and Hall/CRC; 2006. OpenURL

  33. Lin X: SAS Macro GAMM1 to fit generalized additive mixed models using smoothing splines. Harvard University;

    http://scholarjsagotskydev.iq.harvard.edu/amaity/software/gamm1 webcite

    OpenURL

  34. The GLIMMIX Procedure. SAS Institute inc; 2006.

    support.sas.com/rnd/app/papers/glimmix.pdf webcite

    OpenURL

  35. Benjamin MA, Rigby RA, Stasinopoulos DM: Generalized autoregressive moving average models.

    J Am Stat Assoc 2003, 98(461):214-223. Publisher Full Text OpenURL

  36. Cox DR, Gudmundsson G, Lindgren G, Bondesson L, Harsaae E, Laake P, Juselius K, Lauritzen SL: Statistical analysis of time series: some recent developments [with discussion and reply].

    Scand J Stat 1981, 8(2):93-115. OpenURL

  37. Ruoyan M: Estimation of Dispersion Parameters in GLMs with and without Random Effects. Stockholm University; 2004.

    http://www2.math.su.se/matstat/reports/serieb/2004/rep5/report.pdf webcite

    OpenURL

  38. Box G, Jenkins GM, Reinsel G: Time Series Analysis: Forecasting & Control. 3rd edition. Prentice Hall; 1994. PubMed Abstract | Publisher Full Text OpenURL

  39. Zhang S, Qi L: shi jian xu lie fen xi jian ming jiao cheng. Beijing: Beijing jiaotong university press; 2003. OpenURL

  40. Peng RD, Dominici F, Louis TA: Model choice in time series studies of air pollution and mortality.

    J R Stat Soc A Stat Soc 2006, 169(2):179-203. Publisher Full Text OpenURL

  41. Gasparrini A, Armstrong B, Kenward MG: Distributed lag non-linear models.

    Stat Med 2010, 29(21):2224-2234. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  42. Gasparrini A: Distributed lag linear and non-linear models in R: The package dlnm.

    J Stat Softw 2011, 43(8):1-20. PubMed Abstract | PubMed Central Full Text OpenURL

Pre-publication history

The pre-publication history for this paper can be accessed here:

http://www.biomedcentral.com/1471-2288/12/165/prepub