Multiple regression models for predicting total daily pollen concentration in Cartagena

Abstract

The use of meteorological autocorrelation variables and pollen concentrations from previous days, coupled with classification of meteorological data according to multivariate analysis techniques, is shown to improve the predictive power of multiple regression models for daily pollen forecasts. This paper presents an investigation of the meteorological and autocorrelation variables which influence pollen counts in Cartagena, from 1995 to 1999, as a basis for the development of predictive models. The analysis of total pollen concentrations, and especially Chenopodiaceae‐Amaranthaceae, was determined. Initially, forecasting models for total pollen counts were developed, using data from 1995 to 1998, and autocorrelation and meteorological variables. Secondly, predictive models were developed for different meteorological situations, which improved the results by decreasing the number of predictive parameters. Finally, data from 1999 were used to validate the predictive models.

Knowledge of airborne pollen concentration has proved to be of considerable importance for humans in fields such as agronomy, ecology, biology, allergy and occupation hygiene (Galán et al. 1995, Fornaciari et al. 1998). An accurate forecast of these concentrations would benefit pollinosis patients (Bringfelt et al. 1982, Antépara et al. 1995, Norris‐Hill 1995, Dahl & Strandhede 1996), as well as a prediction tool for different crops (Fornaciari et al. 1998).

Daily pollen counts vary as a function of multiple parameters, e.g., flowering patterns or meteorological conditions (Hyde & Adams 1960, Bringfelt et al. 1982, Andersen 1991, Alba et al. 2000, Emberlin et al. 2002). Therefore modelling these systems may be more complex than for other pollution data, where characteristics of emission source can be clearly established. Two approaches to model the movement of particulate matter in the atmosphere can be used: source and receptor orientated models (Di‐Giovanni et al. 1989, Norris‐Hill 1995).

Source orientated models use mathematical algorithms to analyse the dispersal of pollen far away from the source, or to define the delivery process. These models simplify the dispersal process, and require some knowledge of the pollen emission for each single source, or use a well‐defined area accounting for all the regions where it is possible to get pollen to the forecasting place. Emission and dispersion models have been developed to calculate mesoscale distribution of pollen types, where both vegetation maps and pollen sampling data are available (McCartney 1994, Kawashima & Takahashi 1995).

Receptor orientated models are designed to predict concentrations without knowledge of source conditions, which is often the case in aerobiological systems (Norris‐Hill 1995).

Receptor models for pollen concentrations have been used to develop forecasting models using previous time series (Moseholm et al. 1987, Emberlin et al. 1993), and to predict the severity or the start of the pollen season (Driessen et al. 1990, Andersen 1991, Ong et al. 1997, Fornaciari et al. 1998).

Recently considerable attention has been given to statistical models which use mathematical techniques to forecast the dispersion and distribution of airborne particles (Kolehmainen et al. 2001).

The aim of the research reported in this paper is to develop forecasting receptor orientated models using multiple lineal regression techniques for pollen concentration in Cartagena.

In a previous paper (Angosto et al. 2002); we have classified the different wind patterns in our city by means of a two‐step cluster analysis, as previously reported by other authors (Kalkstein et al. 1987). This classification has been used for developing forecasting models, together with other meteorological parameters, such as temperature, relative humidity, pressure, rain, and daily solar radiation. In order to evaluate long‐term trends, a method of time‐series analysis was also used.

MATERIAL AND METHODS

Total pollen concentrations were measured from 1995 to 1999 using a Hirst‐type volumetric sampler (Hirst sampler, Lanzoni VPSS‐2000), with an adjusted flow of 10 l/min. Particles with a size range between 2 and 200 μm in diameter were captured by means of a Melinex tape (cod. 200.700), impregnated with a silicon mixture, and with a driving speed of 2 mm per hour. The sample was processed as previously described (Moreno‐Grau et al. 1998), following the recommended standards of the Spanish Aerobiological Network (REA) (Domínguez et al. 1991).

Data collected in Cartagena (Fig. 1) were recorded as average daily values, expressed as total pollen grains per cubic meter of air.

Fig. 1. Map of Spain with location of Cartagena.

In order to avoid the annual variations in pollen concentrations which could mask the real patterns (Emberlin et al. 1993), the daily pollen content was standardized following Moseholm et al. (1987), using the following equation:

where:

n_p =standardized pollen count

p_r =daily pollen count

p_y =total annual pollen concentration

Meteorological parameters selected for this study were daily minimum, average, and maximum values of temperature, relative humidity, air pressure, and wind speed, accumulated daily rainfall, and accumulative sunshine. Wind speed and wind direction were combined to create two categorical variables; wind range and wind course.

Wind course was obtained by grouping wind direction values into 12 classes, all with 30° of amplitude: North (N), North‐Northeast (NNE), East‐Northeast (ENE), East (E), East‐Southeast (ESE), South‐Southeast (SSE), South (S), South‐Southwest (SSW), West‐Southwest (WSW), West (W), West‐Northwest (WNW), and North‐Northwest (NNW).

Average wind speed values were grouped into four wind ranges: Range 1=0–4.9 km/h, Range 2=5–10.9 km/h, Range 3=11–15.9 km/h and Range 4≥16 km/h.

All statistical analyses were carried out using SPSS 10.0 software. First, Pearson's correlation coefficients and bilateral significance tests were used to explore the relationship between meteorological variables and standardized pollen concentrations. The values of these variables recorded during the previous days were also used in this correlation, which could influence in daily pollen counts.

Multiple regression models were developed using the Backward option in SPSS, which removes at each step any variable with a p‐value less than 0.10, 90% significance level (Vinacua 1997).

RESULTS AND DISCUSSION

A decrease in total pollen counts was observed from 1995 to 1999, except for 1998, which showed a slight increase relative to 1997 (Table I). The same behavior is seen in Chenopodiaceae‐Amaranthaceae (Table I). The highest total pollen counts were recorded during the month of May for 1995, 1996, and 1997 (Table I), but in 1998 and 1999, the peak day occurred earlier (Table I), probably due to extensive flowering of Cupressaceae.

Table I. Values for total pollen and Chenopodiaceae‐Amaranthaceae counted from 1995–1999, and the main meteorological parameters.

Years	1995	1996	1997	1998	1999	1995–1999
Total (pollen grains/m^3 )	16885	11326	11746	13391	11370	64518
N° grains/m^3 of peak day	479	423	659	767	595	‐
Date of peak day	12 May	28 May	6 May	9 Mar	23 Feb	‐
Chenop. (pollen grains/m^3 )	2860	2044	1234	1527	1182	8847
N° grains/m^3 of peak day	171	95	103	99	103	‐
Date of peak day	17 Sep	3 Sep	6 Sep	19 Sep	9 Sep	‐
Mean Pressure (mBar.)	1019	1017	1019	1021	1020	1019
Mean Temperature (°C)	19.91	19.56	19.93	19.16	19.79	19.68
Accumulated Rainfall (mm)	93.4	225.1	229.3	112.8	194.7	171.06
Humidity (%)	79.03	78.21	82.19	80.11	80.12	79.93
Accumulated Sunshine (W7m^2 )	399.61	398.27	389.19	393.37	373.28	390.74
Mean wind speed (km/h)	13.86	13.34	13.25	12.84	14.02	13.47

The highest peak for Chenopodiaceae‐Amaranthaceae was always recorded during the month of September (Table I), coincident with the second flowering of this pollen type in our region, which is the most representative of that period (Moreno‐Grau et al. 1998).

The decrease observed for pollen counts from 1995 to 1999 can be attributed to the high drought suffered in this area, with a maximum accumulated rainfall of 229.3 mm registered during 1997 (Table I), whilst average value for rainfall is 354.2 mm for the period 1975 to 1994. Other factors affecting the decrease could be changing crop species or substantial increase in the main source area.

The other values for meteorological parameters (Table I) are typical of an arid Mediterranean and subtropical climate (Moreno‐Grau et al. 1998).

Daily average total pollen concentrations from 1995 to 1999 are depicted with a dashed line (Fig. 2). Two distinct periods can be observed. The first, from January to August, is the main pollen period. It has saw‐toothed structure and includes taxa with pre‐spring and spring flowering. The second, from September to October, is coincident with the second flowering of Chenopodiaceae‐Amaranthaceae (represented by a continuous line).

Fig. 2. Daily average pollen count from 1995 to 1999. Total pollen (dashed line) ‐ Chenopodiaceae‐Amaranthaceae (continuous line).

On the other hand, temperatures show an increasing tendency from the beginning of the year until summer time, and a decreasing trend from summer to December.

These trends led us to develop two multiple regression models for pollen forecasting, one for total pollen concentrations, from January to August, and a second one for the highest peak of Chenopodiaceae‐Amaranthaceae (September and October).

As a first step, and following recommendations of different authors (Emberlin et al. 1993, Díaz de la Guardia et al. 1998, Stark et al. 1997, Galán et al. 1995, 2000), the relationships between different independent variables and the dependent variable for each model were explored.

Pearson's correlation coefficients, comparing standardized daily average total pollen counts for the first 30 weeks and second flowering of Chenopodiaceae‐Amaranthaceae with different meteorological parameters and variables for the previous days, provides information on which variables correlate with pollen production (Table II). A negative correlation coefficient was found for mean atmospheric pressure (MAP), with a 95% significance level, implying that an increase atmospheric in MAP leads to a decrease in airborne pollen counts. Glassheim et al. (1995) have reported a positive and significant correlation between MAP and pollen concentrations. In our study, the settlement of air masses would reduce the flotation ability for the biggest pollen grains, thus having a stronger effect during the second flowering of Chenopodiaceae‐Amaranthaceae than for total pollen.

Table II. Pearson's correlation coefficients between daily average pollen count meteorological parameters during 1995–1999.

Period (1995–1999)	Total pollen (30 weeks) (n=985)	Chenopod/Amaranth pollen (2^nd peak) (n=295)
Mean Pressure (MP)	−0.078*	−0.177**
Mean temperature (MT)	0.079*	−0.051
Average relative humidity (RH)	−0.217**	−0.087
Mean wind speed (WS)	0.153**	0.152**
Accumulative sunshine (AS)	0.293**	−0.012*
Accumulated daily rainfall (RF)	−0.048	0.155**
Wind range (RA)	0.182**	0.135*
Minimum temperature (MIT)	0.074*	−0.260
Maximum temperature (MT)	0.086*	−0.101
Minimum humidity (RHMIN)	−0.022**	−0.091
Minimum pressure (MIP)	−0.079*	−0.142*
Maximum pressure (MXP)	−0.075*	−0.173**
Season (ST)	0.026	−0.281**
ΔT	0.017	−0.098
ΔT1	−0.042	−0.071
ΔT2	0.004	−0.053
ΔT3	0.024	−0.061
ΔT4	0.053	−0.028
ΔT5	0.024	−0.083
ΔT6	−0.029	−0.046
ΔH1	0.084**	0.049
ΔH2	0.077**	0.078
ΔH3	0.053	0.089
ΔH4	0.067*	0.0130*
ΔH5	0.038	0.091
ΔH6	0.027	0.110
P1ST	0.720**	0.769**
P2ST	0.580**	0.664**
P3ST	0.490**	0.608**

**p≤0.01; *p≤0.05; ΔT=variation between maximum and minimum temperature; ΔT1=thermal variation for the previous day; ΔT2=thermal variation two days before; ΔT3=thermal variation three days before; ΔT4=thermal variation four days before; ΔT5=thermal variation five days before; ΔT6=thermal variation six days before; ΔH1=humidity variation for the previous day; ΔH2=humidity variation two days before; ΔH3=humidity variation three days before; ΔH4=humidity variation four days before; ΔH5=humidity variation five days before; ΔH6=humidity variation six days before; P1ST=total standardized pollen concentration for the previous day; P2ST=total standardized pollen concentration two days before; P3ST=total standardized pollen concentration three days before

For total pollen, relative humidity (RH) and its minimum daily value (RHMIN) show significant negative correlation. This behavior has already been reported by other authors (Herrero & Fraile 1997), and is explained by variation in shape and size of pollen grains according to their hydration level (Blackmore & Barnes 1986). As humidity increases, pollen grains increase their ability for wet deposition.

Wind speed (WS) showed a positive and 99% significant correlation with both total and Chenopodiaceae‐Amaranthaceae pollen, which indicates an increase in pollen concentrations with the increase of wind speed. This fact could be conditioned to the distance between the pollen sampler and plant masses: with a shorter distance, wind could act as a dilution factor.

Accumulative sunshine (AS) is positively correlated with total pollen counts, because of its influence with formation and delivery of pollen grains (Akers et al. 1979). Similar results have been reported by Bringfelt et al. (1982) and Glassheim et al. (1995). However, the second peak of Chenopodiaceae‐Amaranthaceae displayed a low negative Pearson's correlation coefficient with this meteorological parameter, as shown by a gradual decay of solar radiation during the month of September and an increase of this pollen type. Similar results have been reported by Muñoz et al. (2000) for mean temperature and relative humidity.

For both pollen peaks, there is a 99% significant correlation with standardized pollen counts in the previous three days. This fact is highly influenced by pollen delivery pattern.

The next step was to develop multiple lineal regression models. The first regression models were developed using those meteorological parameters with a significant correlation with pollen counts from 1995 to 1998.

The two best regression models for total pollen and the best for each cluster, (first 30 weeks of the year), are from 1995 to 1998 (Table III). Variables were included in the model following a step forward method. Models 1 and 2 for total pollen yield a determination coefficient of 0.54, when the following variables are included: pollen count of the previous day, minimum relative humidity, pollen count of three days before, wind speed and direction, and the variations of maximum and minimum temperature for four and six days before. Although the value of unexplained variance is still high, determination coefficients are higher than those previously obtained where for single variables and for models and data reported by other authors (Bringfelt et al. 1982, Galán et al. 1995).

Table III. Results of multiple regression analysis for the two best models for total standardized pollen (first 30 weeks of each year) (1995–1998) and for each cluster.

Models	r^2	Equations
1	0.539	PST=1.110+0.635*P1ST−0.008*HRMIN−0.007*ΔH1+0.108*P3ST+0.011*WS−1.480*RUM11−0.016*ΔT6
2	0.541	PST=0.995+0.634*P1ST−0.009*HRMIN−0.007*ΔH1+0.111*P3ST+0.011*WS−1.470*RUM11−0.019*ΔT6+0.015*ΔT4
C1	0.685	PST=−0.611+0.772*P1ST+0.079*ΔT5+0.013*ΔH1
C2	0.520	PST=0.615+0.640*P1ST−0.010*ΔT6+0.136*P3ST
C3	0.575	PST=0.856+0.651*P1ST−0.104*ΔT1+0.015*ΔH2
C4	0.435	PST=0.916+0.649*P1ST−0.015*ΔH1
C5	0.592	PST=−0.395+0.527*P1ST−0.009*ΔH1+0.013*ΔH+0.179*P3ST+0.003WS+0.027*ΔT4

RUM11=Categorical variable for wind direction, corresponding to WNW direction

Wind direction appears to be an important parameter for the amount and type of pollen arriving at the pollen sampler (Herrero & Fraile 1997, Moreno‐Grau 1998). For that reason, and in order to achieve a better result for our models, a two‐step cluster multivariate analysis has been used for a daily classification of wind fluxes, previously reported in Angosto et al. (2002). The average and 95% confidence interval data for standardized total pollen concentration includes the period 1995 to 1999 (Fig. 3).

Fig. 3. Average values and 95% confidence interval data for total standardized pollen concentrations for each cluster, from 1995 to 1999.

Statistically significant differences were observed between different clusters, which imply that each wind cluster picks up different pollen source areas. The highest pollen counts were obtained for cluster 1, noted for a predominance of NNW and N wind direction, and including high pressures with intermediate wind speed (Angosto et al. 2002). This fact could be explained due to the situation of the main vegetation masses in respect to the pollen sampler.

The masses of air coming from the sea are represented by clusters 3 and 4. As reported in a previous paper (Angosto et al. 2002), cluster 3 groups low pressure situations with a maximum wind speed in the city and a predominant SSW wind direction. On the other hand, cluster 4 includes high pressure situations, with a maximum of temperature and N and NNE predominant wind directions. Both clusters showed a relatively low pollen count.

The lowest pollen count was achieved for cluster 5, representing the minimum wind speeds without a predominant wind direction.

For these reasons, we consider that cluster variable stands for the influence of wind fluxes, and therefore can be applied to the prediction of pollen concentrations. It is easier to handle this variable for statistical studies and modelling procedures than wind direction variable alone.

The next step was to develop forecasting models for each wind situation (for each cluster). Standardized total pollen concentration was obtained for these models by dividing the original data according to the cluster group corresponding to each day. The best model was determined for each cluster situation (Table III). Except for cluster 4 (C4), all models show stronger determination coefficient than those passed on the full data set. The highest determination coefficient was for cluster 1 (C1), 0.685. Except for cluster 5 (C5), all the cluster models also used a much lower number of predictive variables.

The results of the two best multiple regression models for the second peak of Chenopodiaceae‐Amaranthaceae yielded a similar determination coefficient (Table IV), but the first one uses only two variables, standardized pollen count of two and three previous days, and the second model the same two variables plus the maximum pressure value for this day.

Table IV. Results of multiple regression analysis for the two best models obtained for second peak of Chenopodiaceae/Amaranthaceae (1995–1998) and for each cluster.

Models	r^2	Equations
1	0.603	PST=0.315+0.651*P2ST+0.183*P3ST
2	0.603	PST=32.116+0.639*P2ST+0.178*P3ST−0.031*MXP
C1	0.401	PST=0.946+0.459*P2ST
C2	0.682	PST=0.144+0.658*P1ST+0.354*RF+0.228*P2ST
C3	0.839	PST=92.971+0.828*P1ST−0.961*ST+0.332*R−0.090*MXP
C4	0.554	PST=0.753+0.719*P1ST
C5	0.535	PST=0.322+0.534*P1ST+0.235*P3ST

When the data for this second peak of Chenopodiaceae‐Amaranthaceae were subdivided according to the cluster (see Table IV), stronger r² values were obtained for C2 and C3, 0.839 and 0.682, and for clusters 1, 4, and 5. The predictive equations explain a similar variance to that explained for the whole sample.

Our next step was validation of the best models with data from 1999. The validation was performed for the best model of total pollen count and second peak of Chenopodiaceae‐Amaranthaceae (model 1) (Tables III & IV). Models representing total pollen count by clusters (Table III) were also validated. Models by clusters for the second peak of Chenopodiaceae‐Amaranthaceae were not validated because sample size was too small for cluster subsets.

The dispersion plot of predictions and observed data of the best model for total standardized pollen (Fig. 4), with an adjusted determination coefficient of 0.541 is concordant between forecast and observed pollen counts, as shown by the Pearson's correlation coefficient of 0.97 (99% significant−R² =0.94).

Fig. 4. Validation of the best model for total standardized pollen concentration using data from 1999.

For the calculated and critical values (α=0.05) of variance (F) and mean values (t) (Table V), the critical value is always higher than calculated value. Thus, we cannot reject the null hypothesis of equal means and variances.

Table V. Calculated and critical values (α=0.05) for variance (F) and mean values (t) for both models.

Model	r	r^2	F calculated	F critical	t calculated	t critical
Total pollen	0.97**	0.94	0.53	1.00	0.36	1.65
Chenopodiaceae/Amaranthaceae	0.83**	0.68	1.62	1.90	0.99	1.65

** p≤0.01

Using data from 1999 a dispersion plot of the validation of the best model corresponding to the second peak of Chenopodiaceae‐Amaranthaceae was constituted (Fig. 5). Although the number of samples is not very high, because of a short flowering period, they seem to fit the model, with a correlation coefficient between predicted and observed data of 0.83 (99% statistically significant). As for the previous model for total pollen concentrations, we cannot reject the null hypothesis for equal means and variances (see Table V).

Fig. 5. Validation of the best model for Chenopodiaceae/Amaranthaceae pollen using data from 1999.

The results of the best model validation for total pollen count for each cluster shows that the highest and 99% statistically significant coefficients were achieved for clusters 1, 2, 3, and 5, meanwhile cluster 4 obtained a lower Pearson's correlation coefficient and only 95% significance (Table VI). This fact could be explained by the lower number of days which were gathered by that cluster.

Table VI. Determination and correlation coefficients between predicted and calculated data for 1999, for each cluster situation. Calculated and critic values (α=0.05) for variance (F) and mean values (t).

Model	r	r^2	F calculated	F critic	T calculated	t critic
C1	0.81**	0.66	0.59	2.03	0.45	1.65
C2	0.69**	0.47	0.58	2.05	0.05	1.65
C3	0.87**	0.75	1.80	4.45	0.47	1.71
C4	0.51*	0.26	2.84	2.65	0.24	1.70
C5	0.67**	0.45	0.47	1.60	1.35	1.65

**p≤0.01 *p≤0.05

Statistic parameters F and t for test of variance and mean differences showed that null hypothesis could not be rejected for clusters 1, 2, 3, and 5, meanwhile for cluster 4 null hypothesis for variance differences has to be rejected.

CONCLUSIONS

The use of meteorological autocorrelation variables and pollen counts in previous days clearly improved the results of pollen predictive models, both for total and Chenopodiaceae‐Amaranthaceae counts.

Prior classification of meteorological data (a two‐step cluster analysis) and the development of a separate model for each meteorological situation, also improved the forecasting power of the models by decreasing the number of predictive variables.

Models 1 and 2 developed for the forecasting of total pollen concentration could be used in other cities of the Mediterranean area, with similar meteorological characteristics. Meanwhile, models obtained for each cluster could only be used in those cities of the Mediterranean area with wind patterns similar to those identified for the city of Cartagena.

ACKNOWLEDGEMENTS

Authors acknowledge the Environmental Service Department of the City Town of Cartagena for their support with this study.