where wc = vc + τˆ2. Let vc be the “true” statistical variance, and let vc* be the “estimated” statistical variance under a miss-specified model. We assumed vc* < vc, and we considered 3 cases of underestimation of vc:
- additive bias: vc* = vc − b, for arbitrary vc and b > 0;
- multiplicative bias, equal variances: vc* = k × vc, vc ≡ v; v for arbitrary v and 0 < k < 1;
- multiplicative bias, unequal variances: vc* = k × vc, for arbitrary vc and 0 < k < 1.
In the first and second case, underestimation of the variance does not affect the pooled estimate of α. The mathematical proof is given in the appendix.
In the third case, we explored the impact of multiplicative bias on the estimation of α by a simulation study. We considered 33 = 27 scenarios consisting of 3 components:
- underestimation of 50%, 30%, and 10% which corresponds to k = 0.5, 0.7, 0.9;
- number of sites: C = 15, 20, 90;
- amount of heterogeneity: τ = 0.05, 0.5, 1 corresponding, respectively, to small, medium, and large between-city standard deviations.
For each scenario, we generated 250 replicates of βˆ = (βˆ1,…, βˆC) from model(1). More specifically, we chose values for α, τ2, and vc and we sequentially generated βˆc from N(α, τ2), and then βˆc from N(βc, vc). For C = 90, we set vc equal to the estimates obtained from the NMMAPS reanalysis.28 For C = 15 and C = 20 we took a random sample from the 90 vc NMMAPS estimates. We also set α equal to 0.21 (the pooled NMMAPS estimate for total mortality at lag 128). In summary, each scenario (biased vc*, sample size C, and amount of heterogeneity τ), led to 250 simulated values of βˆ.
For each replicate βˆ, and for each scenario, we calculated the empirical Bayes estimates of τ2 and α, using both vc and vc*, by fitting model (1) and using equation (2). This led to 2 sets of 250 estimates of α, here denoted by â(vc) and â(vc*). Figure 1 shows boxplots of the 250 standardized differences between the 2 α estimates defined as (â(vc) − â(vc*))/SD, where SD is the standard deviation of the 250 estimates of â(vc). Thus, a standardized difference larger than 2 suggests high sensitivity of the pooled estimate to the use of the biased vc* instead of “true” vc.
For all scenarios, the 95% confidence intervals of the standardized differences were always within 2 standard deviations of 0, suggesting that underestimation of vc did not affect â substantially. In 8 scenarios, the distributions of the standardized differences showed more variability, with their maximum absolute differences larger than 2 standard deviations. Seven of those scenarios corresponded to extreme underestimation of 50%. The other scenario was characterized by small sample size, small heterogeneity, and a bias of 30%.
We also investigated the impact of variance underestimation on the NMMAPS national average relative rate estimate. We recalculated the national average relative rate of mortality for a 10 unit increase in PM10 by varying the underestimation parameter k from 0.1 (90% variance underestimation) to 1 (no variance underestimation). Figure 2 shows the NMMAPS national average and heterogeneity estimates as function of the underestimation parameter k. Note that for values of k increasing from 0.1 to 0.6 (variance underestimation from 90 to 40%), τˆ decreased; also, â moved from the un-weighted average of the βˆc(equal to 0.28%; 95% posterior interval equal to .05 – .51) toward the weighted average of the βˆc defined in equation (2) (equal to 0.21%; 95% posterior interval equal to .10 – .31). The latter occurred because when k was very small (large variance underestimation), τˆ2 was large relative to vc* = kvc, and roughly the same weight was assigned to all the city-specific estimates. For values of k increasing from 0.6 to 1 (variance underestimation from 40% to 0%), â was constant because τˆ2 = 0.
In summary, little or no effect was observed when the variance underestimation was less than 40% (k ≥ 0.6) (national average estimate equal to 0.21%). When underestimation was larger than 40% (k < 0.6), â gradually increased with the degree of underestimation toward the un-weighted pooled estimate which is equal to 0.28%. However, for all values of k, the t-ratio âSE(â) remained larger than 2 indicating “statistical significance” of the national average relative rate estimate.
The results of this paper indicate that, in multisite time series studies of air pollution and health, underestimation of the statistical variances vc did not affect the estimate of the pooled effect α substantially. Some sensitivity was observed when the number of sites was small (less than 20), the between-city variability was close to zero, and the underestimation was larger than 40%. In addition, our results suggest that variance underestimation of as much as 40% in the NMMAPS studies would not change the conclusions.
The mathematical results reported in this paper are generalizable to any two-stage normal-normal hierarchical model, and therefore can be applied to any meta-analysis of study results or clinical trials that uses such modeling approach. The simulation results reported in this paper can also be applied to any meta-analysis when the size of the overall parameter α with respect to the heterogeneity parameter τ and the city specific standard deviations, vc, is similar to the ones considered here. However, even in a meta-analysis with very different numerical values for α, τ, and vc, our simulation studies can be easily modified to reflect these other situations.
Our underestimation parameter k, here specified a priori, may depend upon site-specific characteristics. In each city, data characteristics may lead to different autocorrelation structures of the residuals, and therefore to different degrees of adjustment for confounding factors. Therefore, variance underestimation attributable to model mis-specification might be more or less severe in the different cities, leading to different values of k. Our simulation study considers both minor (10%) and major (50%) underestimations of standard errors and therefore should encompass the more severe cases of model mis-specification.
The robustness of air pollution effect estimates to variance underestimation applies only to pooled-effect α. Empirical Bayes estimates of the site-specific relative rates, here denoted by β˜c, are affected by the underestimation of the vc. To see this, consider their definition:
where Sc = vc−1/(τˆ−2 + vc−1) is the shrinkage factor. Note that Sc is not invariant to underestimation of the vc. In this case, underestimation of vc leads to overestimation of Sc and therefore leads to β˜c which rely too heavily on βˆc. Thus, variance underestimation leads to an overestimation of the heterogeneity of the air pollution effects, and therefore to under-shrinkage of the city-specific empirical Bayes estimates toward their overall mean â.
Unfortunately, because the true statistical variances are unknown, the distinction between additive and multiplicative bias is not straightforward. One possible approach is to specify 2 alternative but comparable site-specific regression models and to compare the estimates of vc under the 2 models. For example, Dominici et al27 estimated βˆc and vc under a GLM with natural cubic splines and under a GAM with smoothing splines; they then compared the vc under these 2 modeling approaches. In time series studies of air pollution and health, a more robust estimate of the statistical variance under a GAM can be obtained by using gam.exact.27 In other situations, robust variance estimates can be obtained by using asymptotic theory29 or bootstrap30 or sandwich estimates.31 Robust and model-based variance estimates can then be compared for a small number of sites to gain insight on the form of bias.
The results in this paper are presented using an empirical Bayesian approach to estimation, in which the Bayesian estimate of α is obtained by plugging a point estimate of τ2 into Equation (2). However our results apply also if a fully Bayesian version of the above model is fit with noninformative priors on τ2 and α.
Finally, the results of this paper apply under the assumption that the normal approximation to the likelihood function at the first stage of the hierarchical model is appropriate. Asymptotically, this approximation has an accuracy proportional to the number of days with available data in each city.32 In time series studies of air pollution and health the asymptotic normal approximation is generally accurate; however, additional work is needed to extend such results to distributions other than the normal and to examine the sensitivity of inferences if the normal approximation is not accurate.
Details on the estimation of τ2 and α. The empirical Bayes estimate of α was obtained by first computing the restricted maximum likelihood estimate of τ2, which was then plugged into equation (2). The restricted maximum likelihood estimate of τ2 was obtained by maximizing the following likelihood function
where â(τ2) = [Σcβˆcwc(τ2)−1]/[Σcwc(τ2)−1] and wc(τ2) = vc+τ2. This maximization can be done iteratively using the E-M algorithm.18
Mathematical arguments concerning the robustness of the pooled estimate to variance underestimation under cases 1 and 2. Under scenario 1, we start by showing that additive bias does not affect estimation of α for arbitrary vc. Let τ2* be the maximum likelihood estimate (mle) of τ2 obtained by maximizing the likelihood (3) with vc* = vc − b instead of vc. Note that estimation of α as defined in equation (2) depends on vc only through vc + τ2. Therefore, we simply need to show that the estimate for τ2* conditional on vc* is b units more than the estimate for τ2 conditional on vc. In this case, vc + τˆ2 = vc* + τˆ2* and therefore â(τˆ2) is unaffected by underestimation of vc. We define τ2* = τ2 + b and we maximize the likelihood (3) (with vc replaced by vc*) with respect to τ2*. By the invariance property of the maximum likelihood estimates, if τˆ2* is the mle then τˆ2 is also the mle. Therefore, vc + τˆ2 = (vc − b) + (τˆ2 + b) = vc* + τˆ2*, τˆ2* = τˆ2 + b.
Under scenario 2, we assume multiplicative bias with equal variances vc = v. Here vc* = k × vc can be rewritten as vc* = vc − (1 − k) × vc. Because vc = v, then (1 − k) × vc = (1 − k) × v = b and therefore vc* = vc − b and so the results above for additive bias apply here also.
The authors thank Jonathan M. Samet for comments.
1.Katsouyanni K, Touloumi G, Spix C, et al. Short term effects of ambient sulphur dioxide and particulate matter on mortality in 12 european cities: results from time series data from the APHEA project. Brit Med J
2.Samet JM, Zeger SL, Dominici F, Dockery D, Schwartz J. The National Morbidity, Mortality, and Air Pollution Study Part I: Methods and Methodological Issues
. Cambridge, MA: Health Effects Institute, 2000.
3.Thomas A, Spiegelhalter DJ, Gilks WR. Bugs: A program to perform Bayesian inference using Gibbs sampling. In: Bernardo JM, Berger JO, Dawid AP, Smith AFM, eds. Bayesian Statistics 4. Proceedings of the Fourth Valencia International Meeting. Oxford: Clarendon Press; 1992;837–842.
4.Gilks WR, Richardson S, Spiegelhalter DJ, eds. Markov Chain Monte Carlo in Practice
. London: Chapman & Hall, 1996.
5.Burnett R, Krewski D. Air pollution effects of hospital admission rates: A random effects modelling approach. Can J Stat
6.Roemer W, Hoek G, Brunekreef B, Kalandidi A, Pekkanen J. Daily variations in air pollution and respiratory health in a multicenter study: the peace project. Eur Respir J
7.Dominici F, Samet JM, Zeger SL. Combining evidence on air pollution and daily mortality from the twenty largest US cities: A hierarchical modeling strategy (with discussion). J R Stat Soc [Ser A]
8.Zanobetti A, Schwartz J, Dockery D. Airborne particles are a risk factor for hospital admissions for heart and lung disease. Environ Health Perspect
9.Schwartz J. Assessing confounding, effect modification, and thresholds in the associations between ambient particles and daily deaths. Environ Health Perspect
10.Dominici F. Air pollution and health: what can we learn from an hierarchical approach? [Invited commentary]. Am J Epidemiol
11.Lindley DV, Smith AFM. Bayes estimates for the linear model (with discussion). J R Stat Soc [Ser B]
12.Samet JM, Zeger SL, Berhane K. The Association of Mortality and Particulate Air Pollution
. Cambridge, MA: Health Effects Institute, 1995.
13.Kelsall J, Samet JM, Zeger SL. Air pollution, and mortality in Philadelphia, 1974–1988. Am J Epidemiol
14.Hastie TJ, Tibshirani RJ. Generalized Additive Models
. New York: Chapman & Hall, 1990.
15.McCullagh P, Nelder JA. Generalized Linear Models (Second Edition)
. New York: Chapman & Hall, 1989.
16.Hedges LV, Olkin I. Statistical Methods for Meta-analysis
. Orlando, FL: Academic Press, 1985.
17.DerSimonian R, Laird N. Meta-analysis in clinical trials. Control Clin Trials
18.Carlin BP, Louis TA. Bayes and Empirical Bayes Methods for Data Analysis
. New York: Chapman & Hall, 1996.
19.Morris CN, Normand SL. Hierarchical models for combining information and for meta-analysis. In: Bernardo JM, Berger JO, Dawid AP, Smith AFM, eds. Bayesian Statistics 4
. Oxford: Oxford University Press; 2003;321–344.
20.Gelman A, Carlin J, Stern H, Rubin D. Bayesian Data Analysis
. London: Chapman & Hall, 1995.
21.Buja A, Hastie T, Tibshirani R. Linear smoothers and additive models. Ann Stat
22.Green PJ, Silverman BW. Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach
. London: Chapman & Hall, 1994.
23.Ramsay T, Burnett R, Krewski D. The effect of concurvity in generalized additive models linking mortality to ambient particulate matter. Epidemiology
24.Chambers J, Hastie T. Statistical Models in S
. London: Chapman and Hall, 1992.
25.Lumley T, Sheppard L. Time series analyses of air pollution and health: Straining at gnats and swallowing camels? Epidemiology
26.Samet J, Dominici F, McDermott A, Zeger S. New problems for an old design: Time series analyses of air pollution and health. Epidemiology
27.Dominici F, McDermott A, Hastie T. Improved semi-parametric time series models of air pollution and mortality. Johns Hopkins University, Dept. of Biostatistics, Technical Report
28.Dominici F, McDermott A, Zeger SL, Samet JM. Airborne particulate matter and mortality: Time-scale effects in four us cities. Am J Epidemiol
29.van der Vaart AW. Asymptotic Statistics.
New York: Cambridge University Press, 1998.
30.Efron B, Tibshirani RJ. An Introduction to the Bootstrap
. New York: Chapman & Hall, 1993.
31.Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika
© 2004 Lippincott Williams & Wilkins, Inc.
32.Daniels M, Kass R. A note on first-stage approximation in two-stage hierarchical models. Sankhya Series B