Feature pubs.acs.org/est
Environmental Statistics Revisited: Is the Mean Reliable? Chantal de Fouquet* Mines ParisTech, Centre de Géosciences-Géostatistique, Ecole Nationale Supérieure des Mines de Paris, 35, rue saint-Honoré, 77305 Fontainebleau, France Statistics are commonly used in order to assess the state of environment and its evolution,1−4 and environmental policy5 makes reference to annual mean concentration or fractiles over the year.6−9 However are elementary statistics such as the mean always reliable, for example when the data are collected irregularly in time or in space? First of all, is the definition of the “mean” always clear?
The error standard deviation is
σ(m − m*) = s / n
Expressed in %, the relative error standard deviation is 100 × σ(m − m*)/m*. Minimum, maximum, mean, and standard deviation are expressed in the same units as the data. The variance has the dimension of data square. Relative standard deviation is dimensionless. The statistical summary of the daily concentrations expressed in μg·m−3 at Arles station is as follows: • size: nday = 364, minimum: 4.9 and maximum: 61.9;
■
EXAMPLE FROM AN EXHAUSTIVE TIME-SERIES Let us consider the time-series of daily nitrogen dioxide concentrations (NO2), an air pollutant, measured at a monitoring station in Arles (France) during the year 2005 (Figure 1). Measurements are derived from the air flux through the sensor. Data represent thus the averaged concentration during one day. As December 31 measurement is missing, the considered “year” will henceforth involve 364 days. Let z denote the concentration. Data are indexed by the day number i. Statistical software gives the following summary of the “population”: • size n, minimum and maximum;
• sample mean: zday = 26.5; • variance sday2 = 109.01 and standard deviation sday = 10.4. The estimated distribution mean is mday* = 26.5 (eq 3) and the error standard deviation is σday(m − m*) = 0.55(eq 6), with a relative error standard deviation of 2.06%. Weekly concentrations are the average of seven successive daily concentrations. The statistical summary of the 52 week concentrations is as follows: • size: nweek = 52, minimum: 18.5, maximum: 44.9;
• sample mean:
z̅ =
1 n
• sample mean: zweek = 26.5;
n
∑ zi
• variance Sweek2 = 39.6 and standard deviation sweek = 6.3. Weekly concentrations are less variable than the daily ones: the minimum increases, the maximum decreases, and variance or standard-deviation also decreases. The fluctuations on the weekly time-series (Figure 2) appear effectively more reduced than on the daily ones (Figure 1). The sample mean remains unchanged because the 52 weekly concentrations correspond exactly with the 364 daily concentrations, averaged by interval of seven days. Thus zday = zweek and from eq 3 the estimated mean remains unchanged. But calculated from the weekly concentrations, the associated error standard deviation σweek(m − m*) is now 0.87. The precision is worse than that from the daily concentrations. Let us now consider the 26 2-week (i.e., 14-day) and the 13 4-week (28-day) intervals. The statistical summary of these two populations is presented in Table 1. As the 364-day “year” is exactly divided in an integer number of days, or 1-, 2-, or 4-week intervals, the four populations have exactly the same sample mean. But the precision on the estimated mean gets worse when the time interval at which the data relate (also called “support”) increases. The error standard deviation on the mean is 0.94 μg/m3 for 2-week data, and 1.14 μg/m3 for 4-week data.
(1)
i=1
• variance:
s2 =
1 n−1
(6)
n
∑ (zi − z ̅)2 (2)
i=1
The variance characterizes the dispersion of the population around its mean. Its positive root s is the standard deviation. • variation coefficient: s/z.̅ It characterizes the relative dispersion. The mean m of the underlying distribution is estimated by the mean of the data population (10), which is written: (3) m* = z ̅ The variance of the estimation error on the mean, denoted σ2(m − m*), characterizes the precision of this estimation. As the mean m is a deterministic parameter, the error variance is:
σ2(m − m*) = Var(m*)
(4)
It is calculated from the sample variance
σ2(m − m*) ≃ s 2 /n
(5) © 2011 American Chemical Society
Published: December 16, 2011 1964
dx.doi.org/10.1021/es2024143 | Environ. Sci. Technol. 2012, 46, 1964−1970
Environmental Science & Technology
Feature
Figure 1. Time-series of daily NO2 concentration expressed in μg/m3, measured at Arles (France) in 2005. In abscissa time is expressed as the day number in the year.
Figure 2. Time-series of weekly NO2 concentration. Same abscissa and ordinate scales as in Figure 1.
Table 1. Statistical Summary of Concentrations Expressed in μg·m−3 and Associated Error Standard Deviation on the Estimated Mean (Data Are Daily Values and 1-, 2-, or 4-Week Averages) time support in days
sample size
minimum
maximum
1 7 14 28
364 52 26 13
4.9 18.5 20.1 20.6
61.9 44.9 38.5 34.0
sample mean
sample standard deviation
variation coefficient %
error standard deviation
26.5
10.4 6.3 4.8 4.1
39 24 18 16
0.55 0.87 0.94 1.14
Now can the concentrations z1,...,zn be considered as a realization of random variables Z1,...,Zn with same mean and variance? The seasonal variations visible in Figure 1 and Figure 2 indicate rather that the mean and the variability vary significantly over the year. In fact in the urban areas of the northern hemisphere, the nitrogen dioxide concentrations are generally larger in winter, due to atmospheric photochemical reactions and additional fossil hydrocarbon combustion for heating. • Hypothesis 2: the random variables Z1,...,Zn are jointly independent. The link between two random variables is described by their correlation coefficient, which varies between +1 and −1. It is maximum when the variables are deterministically linked with a linear relation and vary in the same way. It is zero between independent variables and it becomes negative when the variables vary in the opposite direction. The correlation coefficient between concentrations at two successive days (i and i + 1) is 0.55, and it decreases to 0.10 for
What is the meaning of these results? If the precision on the estimated mean improves when the year is divided into smaller intervals, should the precision really get better considering hourly concentrations, or even minutes or seconds related values? Furthermore the 364 daily concentrations are known, and thus the “annual” mean is perfectly known. It is the sample mean (26.5 μg·m−3) identical for the four data sets. What represents then the nonzero error standard deviation given by eq 6? In other words, would there not be some divergence between statistics and physical reality?
■
BACK TO THE STATISTICS COURSE To solve this paradox, let us explicitly discuss on which hypotheses eq 3 to eq 6 are based. • Hypothesis 1: the n-sample is drawn from an unique distribution Z, with expectation m and variance s2. “Expectation” is another word for “probabilistic mean” of a distribution. 1965
dx.doi.org/10.1021/es2024143 | Environ. Sci. Technol. 2012, 46, 1964−1970
Environmental Science & Technology
Feature
Figure 3. Zoom on NO2 daily concentrations between January 30 and March 1. Time is expressed as the day number. Black points indicate the Sunday lower concentration.
a 4-day delay (i and i + 4). The light periodic component (coefficient of 0.17 between days i and i + 7) reflects the weekly periodicity of urban traffic. A zoom on daily concentrations (figure 3) between January 30 and March 1 shows that even during the higher winter concentrations the variations are not totally erratic. Thus the daily concentrations cannot be supposed independent along time. Another thing does not work. Let us design by “regularized” concentration its average on a fixed support. The variance of the concentration regularized on a k-day support is given by the general expression ⎛ 1 Var⎜⎜ ⎝k
k
⎞
∑ Zi⎟⎟ = i=1
⎠
1 k2
k
∑ VarZi + i=1
1 k2
concentration, which is lower than twice the variance of the 2-week concentration and so on. At least one of the previous hypotheses is not consistent with the statistical summaries. • Hypothesis 3: the “parameter” to be estimated is the expectation m of the concentrations distribution. Does the probabilistic mean represent really the pertinent quantity in order to characterize the annual average? What is the meaning of the “mean”? As the time-series of the 364 daily concentrations during the “year” is complete, there is no uncertainty on its average. The standard deviation of the associated estimation error should be zero. But things are different about the uncertainty on the expectation (the probabilistic mean) m of the random variable Z, sampled by the 364 daily values. Let the seasonal variations be neglected. In the absence of a systematic evolution of concentrations along time (and under an additional hypothesis of “ergodicity” of the random variables) the expectation m is the limit of the arithmetic average of the concentrations during an increasing number of years. Assuming a constant expectation is strictly consistent with an annual average different in 2005 and 2006 for example. The “annual mean” for year T, denoted ZT, can properly be defined as the average of all the daily concentrations during T: ZT = Zday. Thus, the difference ZT − Zday is zero and the estimation error on the “annual mean” ZT from the average of the exhaustive time-series of daily values is zero. And the same for the estimation from regularized concentrations of any support d meeting the condition d is an integer divisor of the year. Moreover the definition of the annual mean ZT as a time average is consistent with the presence of seasonal variations.
k
∑
Cov(Zi , Zj)
i,j=1 j≠i
(7)
where Cov(Zi,Zj) denotes the covariance between the variables. If the k random variables Z i have same variance s 2 (Hypothesis 1), the covariance between Zi and Zj is Cov(Zi,Zj) = s2ρij (with ρij their correlation coefficient), and the variance of the regularized concentration becomes k ⎛ k ⎞ 1 1 2 1 2 ⎟ ⎜ Var⎜ ∑ Zi⎟ = 2 ks + 2 s ∑ ρij k i,j=1 ⎝ k i=1 ⎠ k j≠i
⎛ ⎞ ⎜ ⎟ k ⎟ 1 2⎜ 1 = s ⎜ + 2 ∑ ρij⎟ k i,j=1 ⎟ ⎜k ⎜ ⎟ j≠i ⎝ ⎠
■
(8)
POSSIBLE CONFUSION FOR THE DEFINITION OF THE “MEAN”? If it seems more pertinent to characterize the environment during the year T by the time average ZT, a perfectly defined quantity, rather than by the probabilistic mean m supposed to be constant, why is then the “formula” (eq 6) of the error standard deviation on the mean so generally applied and even promoted, for example by European regulation? Because, from N regularly spaced measurements Zi over the year (e.g., one measurement a week or a month) the estimator ZT̅ * of the annual average ZT and the estimator m* of the
If the variables are independent (Hypothesis 2) or just not correlated, then all correlation coefficients are zero (for i ≠ j) and thus all correlation terms too. Then the variance of the sum is equal to the sum of the variances, and the variance of the k-day regularized concentration is s2/k. If Hypothesis One and Hypothesis Two are verified, the variance of the k-day regularized concentration should vary proportionally to 1/k. Now the sample variance is 109.0 μg2·m−6 for daily concentrations, 39.6 for 1-week, 23.0 for 2-week, and 17.0 for 4-week concentration. Thus the variance of the daily concentration is lower than 7 times the variance of the 1-week 1966
dx.doi.org/10.1021/es2024143 | Environ. Sci. Technol. 2012, 46, 1964−1970
Environmental Science & Technology
Feature
expectation m are both equal to the sample average:
1 ZT* = Z̅ and m* = Z̅ , with Z̅ = N
set in the middle of each semester (mean of campaign 1: 27.2).
N
• the mean of the campaign 3 preferential sampling is larger (30.1) and exceeds the limit value of 30 μg/m3 for vegetation protection. At the same measurement number and total duration, spreading the data throughout the year improves the estimation of the annual average (Table 2), even though the sample standard deviation increases and consequently the error standard deviation derived from eq 6. Indeed seasonal variations increase the dispersion of data and thus the predicted error standard deviation. Statistics do not consider the measurement sequence. All occurs as if the measurements are “disorderly”, without the seasonal variations taken into account. The deliberate increase of the measurement frequency during the winter large concentrations (campaign 3) induces an overestimation of the mean, called “bias”. In practice it is not always possible to get a regular sampling schema by removing data from an irregular one. Moreover removing data means losing information and during the statistical treatments it is not always possible to verify if an irregular sampling schema is deliberate, or not. Is it then possible to correct the effects of an irregular sampling without losing information?
∑ Zi i=1
(9)
However the estimation errors Z̅ T − Z̅ and m − Z̅ are different and thus the associated standard deviations are different. Moreover, the error standard deviation on the expectation m given by eq 6 supposes the absence of time correlation. A “confidence interval” around the estimated mean is commonly deduced from the error standard deviation. Of course its validity relies on the validity of the error standard deviation as well as of the normal distribution of the error.
■
NUMERICAL EXPERIMENT Let us now address another aspect of the annual average estimation. Not all survey sites are provided with a permanent monitoring station. At most air or water quality survey sites, only a few measurements are made each year. What is the consequence of the sampling timing and of a possible preferential sampling? Temporary NO2 measurements are made in complement to the fixed monitoring stations. “Passive diffusion samplers” are exposed during one week at some carefully chosen places. In the absence of measurement errors, data are the averaged weekly concentrations. Annual campaigns involve several measurement phases. To reduce the duration of installing and removing the samplers, the phases generally involve 2 or 3 successive weeks. The total measurement duration must be long enough in order to approach the annual average with a “good” precision. For the Arles station, let us simulate different campaigns with a total of 12 weeks (about one-fourth of the year). Data are the associated 12 weekly averages, and the other concentrations are forgotten. The first two campaigns are regularly distributed throughout the year with respectively two 6-week phases or four 3-week phases. The phase duration of the third campaign is reduced to 2 weeks, and two phases are added at the beginning of the year, in order to better survey large winter concentrations (Figure 4).
■
GEOSTATISTICAL ESTIMATION OF A TEMPORAL AVERAGE A simple idea consists of weighting the data in order to reduce the influence of closer measurements and increase that of less frequent ones. Weighting involves replacing eq 1 by n
* = Z̅2005
∑ SiZi i=1
(10)
where the “weights”S i are not anymore uniformly equal to 1/n, but depend on the measurement durations and dates. How are the weights determined? In geostatistics, “optimal” weights S i are derived according to the following criteria: • no bias: the expectation of the estimation error is zero; • optimality: the variance of the estimation error is minimal. The estimator thus defined is called “kriging”.11,12 Kriging enables estimating an unmeasured concentration as well as any temporal average (monthly, etc.) and provides the associated error standard deviation. Kriging also allows estimating the expectation m; the associated precision will differ from that on the annual average. Kriging needs the time covariance of concentrations, or more generally their “variogram”. The variogram is the expectation of the quadratic deviation of concentrations as a function of the time intervals τ between dates t and t + τ
1 E⎡⎣(Z(t + τ) − Z(t ))2 ⎤⎦ (11) 2 If temporal correlation is present, concentrations Z(t) at time t and Z(t + τ) at time t + τ are close at small time intervals τ and the variogram is low. At greater intervals, the amplitude of the deviation increases and the variogram is larger. A periodic component (annual, weekly, etc.) present in the concentrations is also visible on the variogram. The variogram of NO2 daily concentrations calculated for time intervals until 6 weeks shows a rather large variability between successive day measurements, a correlation range about 4 days, and a slight weekly periodicity γ(τ) =
Figure 4. Measurement phases (in dark) for the simulated campaigns. Time is indicated with the week number.
The statistics of the three surveys are presented in Table 2. Logically, • the true annual average (26.5) is better approached by shorter phases regularly distributed per quarter (mean of campaign 2: 26.7) than by only two longer phases 1967
dx.doi.org/10.1021/es2024143 | Environ. Sci. Technol. 2012, 46, 1964−1970
Environmental Science & Technology
Feature
Table 2. Statistical Summary of the Weekly Concentrations Expressed in μg·m−3 and Associated Error Standard Deviation on the Estimated Mean, for the Three Measurement Campaigns regular 2 × 6 weeks regular 4 × 3 weeks preferential
sample size
minimum
maximum
mean
sample standard deviation
error standard deviation
12 12 12
20.9 19.7 19.7
39.6 44.8 44.9
27.2 26.7 30.1
5.3 6.6 9.1
1.53 1.91 2.63
Figure 5. Sample variogram of daily NO2 concentrations. (a) Calculation until a 40-day interval shows a “range” about 4 days and a slight weekly (7 days) periodicity, and (b) until a 530-day interval (on the 2005−2006 time-series) shows the annual periodicity. Horizontal line indicates the sample variance.
are generally larger in winter (Figure 6) but the contrary occurs too. To better monitor the large nitrate concentrations some agencies double the measurement frequency in winter (Figure 6). The statistical mean assigns the same weight 1/18 to all the data, i.e., the same as a weight 2/18 to each of the six winter months and 1/18 to each of the six summer months. In the absence of time correlation kriging weights would be all equal too, whereas for a linear variogram they would be exactly proportional to the sampling interval (“segment of influence” weighting, Figure 7). A detailed study on numerous time-series13,15 showed that for nitrates or nitrites concentrations, kriging weights calculated from the time-series own variograms are mostly close to the “segment of influence”
(Figure 5a); the annual periodicity becomes visible at larger calculation intervals (Figure 5b). The variogram characterizes thus synthetically how a phenomenon is structured:12 amplitude of the different components of the variability, ranges, and periodicities.
■
APPLICATION TO WATER QUALITY REPORTING Are previous developments only sophistication or can the environmental reporting be modified? For nutrient concentrations in rivers, Bernard-Michel et al.13 and Polus et al.14 compared the usual statistical reporting with the geostatistical estimation of the annual average. In metropolitan France, seasonal concentrations in rivers are more or less marked according to substances. Nitrate concentrations 1968
dx.doi.org/10.1021/es2024143 | Environ. Sci. Technol. 2012, 46, 1964−1970
Environmental Science & Technology
Feature
the relative deviation between usual reporting and geostatistical estimation as a function of the number or annual measurements for nitrate concentrations on a set of stations on the Loire basin. The deviation was calculated rather than the difference, because the sign of the difference changes with the stations or the years. Mean relative deviation is generally of 5−10% and raises 20% for preferential sampling with seven samples a year. When sampling varies between years, the evolution of the annual average can change following the estimation method. In some cases, the actual sign of the variation is modified.15 Only a detailed analysis of the data could highlight that the evolution of the statistical mean reflects in fact the sampling changes rather than the concentration variations.
■
CONCLUSION The calculation method is important for environmental reporting. The differences between standard calculations and kriging are not negligible, and become important to characterize the temporal evolution when sampling changes with time. Kriging allows all the measurements to be retained while, in most cases, avoiding bias on the estimated mean and without concern whether oversampling affects large or low concentrations. There is therefore need for revisiting environmental statistics.
Figure 6. Nitrates and nitrites concentration over time measured at a station on the Loire (France) in 1995. From ref 15.
Figure 7. “Segment of influence” weighting. For the estimation of the annual average the weight of measurement i is equal to Si divided by the total length ∑ Si.
■ ■
ones: 1/12 to one-month spaced data and 1/24 to two-week spaced data, i.e., the same weight 1/12 to each month. But weights can be different for substances with a more erratic temporal variability. For the station of Figure 6 the usual reporting overestimates the nitrates annual average of about 13% compared to kriging: the additional precautionary winter measurements result in a bias on the estimated mean. The kriging correction of preferential oversampling is all the more important because the seasonal variation of the nitrites is here opposite. With the classical reporting, the precautionary preferential sampling for the nitrates results for the nitritesmore critical for human healthin an under-estimation of the annual average of about 10%. The error standard deviation deduced from the statistics is here 135% of the geostatistical one for the nitrates and 76% for the nitrites. The precision derived from the statistical calculations should be considered with carefulness: from the same sampling schema it can underevaluate or overevaluate the actual uncertainty, depending on the substance! Figure 8 shows
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected].
ACKNOWLEDGMENTS Thanks to INERIS for transmitting the data of the station in Arles, acquired by Airfobep, and extracted from the BDQA managed by ADEME. The geostatistical estimation of nutrient concentrations in rivers was supported by the French Environment Office. Thanks to H. Beucher and to one Reviewer for suggested improvements and to Ph. Le Caër for graphic help.
■
REFERENCES
(1) AEE. L’environnement en Europe: état et perspectives 2010 − Synthèse; Office des publications officielles de l’Union européenne, Agence Européenne pour l’Environnement: Copenhague, 2010. (2) Observation Et Statistiques; http://www.statistiques. developpement-durable.gouv.fr/theme/environnement/1097.html (3) Baumgardner, R. E.; Lavery, T. F.; Rogers, C. M.; Isil, S. S. Estimates of the Atmospheric Deposition of Sulfur and Nitrogen Species: Clean Air Status and Trends Network, 1990−2000. Environ. Sci. Technol. 2002, 36 (12), 2614−2629, DOI: 10.1021/es011146g. (4) Wang, C.; Corbett, J. J.; Firestone, J. Modelling Energy Use and Emissions from North American Shipping: Application of the Ship Traffic, Energy, and Environment Model. Environ. Sci. Technol. 2007, 41 (9), 3226−3232, DOI: 10.1021/es060752e. (5) Environmental Policy: Past, Present, and Future. A Special issue. Environ. Sci. Technol. 2011, 45 (1). (6) http://ec.europa.eu/environment/policy_en.htm. (7) http://eur-lex.europa.eu/LexUriServ/LexUriServ.do?uri= CELEX:32008L0050:EN:NOT. (8) http://ec.europa.eu/environment/water/participation/notes_ en.htm. (9) Aneja, V. P.; Schlesinger, W. H.; Erisman, J. W. Effects of Agriculture upon the Air Quality and Climate: Research, Policy, and Regulations. Environ. Sci. Technol. 2009, 43 (12), 4234−4240, DOI: 10.1021/es8024403. (10) Saporta, G. Probabilités, Analyse des Données et Statistique; Editions Technip: Paris, 1990.
Figure 8. Relative deviation between classical statistics and geostatistical estimation of the annual average as a function of the number of measurements during the year. Nitrate concentration from measurement stations on the Loire basin between 1985 and 2002. For each sample size, the cumulated number of stations is indicated. 1969
dx.doi.org/10.1021/es2024143 | Environ. Sci. Technol. 2012, 46, 1964−1970
Environmental Science & Technology
Feature
(11) Matheron, G. Traité de Géostatistique Appliquée; Editions Technip: Paris, 1962. (12) Chilès, J.-P.; Delfiner, P. Geostatistics: Modelling Spatial Uncertainty; Wiley: New York, 1999. (13) Bernard-Michel, C.; de Fouquet, C. Geostatistic indicators of waterway quality for nutrients. In Geostatistics Banff 2000, volume 2; Leuangthong, O., Deutsch, C., Eds.; Springer: New York, 2005. (14) Polus, E.; de Fouquet, C.; Flipo, N.; Poulin, M. Caractérisation spatiale et temporelle des Masses d’Eau Cours d’Eau. Rev. Sci. Eau. 2010, 23 (4), 415−429. (15) Bernard-Michel, C. Indicateurs géostatistiques de la pollution dans les cours d’eau. Thèse de Doctorat (PhD thesis); Ecole Nationale Supérieure des Mines de Paris, 2006.
1970
dx.doi.org/10.1021/es2024143 | Environ. Sci. Technol. 2012, 46, 1964−1970
