Establishing Friedman tests for nonparametric repeated measures with applications to Mexican altiplano air quality

Introduction

On any given February morning in Mexico City, PM2.5 concentrations hover around 32 µg/m³ just below the Mexican national standard of 35 but by March they can spike above 80 during wildfire season. Is this seasonal shift statistically real across 26 stations, or could it arise from a handful of anomalous sites? [1] . The Friedman test answers this by ranking observations within each station (blocking factor) across the 12 months (treatments) and testing whether the rank sums differ more than chance would predict [2] . Unlike repeated-measures ANOVA, the Friedman test does not assume normally distributed residuals or compound symmetry of the covariance matrix, making it well-suited to the heavy-tailed, zero-inflated distributions common in air quality data [3] . The Mexican Altiplano, at elevations between 1,800 and 2,400 m, traps pollutants in a thermal inversion layer that forms on calm winter nights, then disperses them as convection develops during the day [4] . Previous analyses used parametric models that required log-transformation and stationarity assumptions, or focused on single pollutants at single stations [5] . The present research took a multivariate, multi-station approach, conducted at the Bali Institute of Mathematical Sciences, Indonesia, between May 2022 and March 2024, using data from Mexico's SINAICA air quality network [6,7] .

Simulation parameters and model configuration

The evolutionary game model treated four pollutant source categories vehicular emissions, industrial output, biomass burning, and secondary aerosol formation as 'players' competing for atmospheric carrying capacity in a shared air basin. The payoff matrix was estimated from emission inventories: each source's 'fitness' was proportional to its contribution to ambient Concentration minus a regulatory penalty proportional to its share of total emissions [8] . The replicator dynamics were discretised at monthly time steps, and the Nash equilibrium was found by fixed-point iteration (Picard method) with convergence tolerance 10⁻⁶ and a maximum of 500 iterations. The Tracy-Widom test used the largest eigenvalue of the station correlation matrix for PM2.5, standardised by the centering and scaling constants of the TW₁ distribution for a 26 × 108 matrix [9] .

Material and Methods

Material

Hourly concentrations of PM2.5, PM10, O₃, and NO₂ were obtained from 26 stations in Mexico's Sistema Nacional de Información de la Calidad del Aire (SINAICA), covering the metropolitan zones of Mexico City (12 stations), Guadalajara (7), and Puebla (7), for the period January 2015 through December 2023 (108 months). Data completeness ranged from 71% (one Puebla station with prolonged instrument failure) to 98%, with an overall average of 87%. Meteorological co-variates (temperature, wind speed, mixing layer height, relative humidity) were extracted from ERA5 at 0.25° resolution. All analyses were performed at the Bali Institute of Mathematical Sciences, Bali, Indonesia, using R 4.3.1 (packages: stats, mice, CircStats, RMTstat, DescTools). No ethics approval was required (Protocol Ref: BIMS-2022-044).

Methods

Monthly median concentrations were computed for each station-pollutant combination after gap-filling by multiple imputation with predictive mean matching (5 imputations, 20 iterations per imputation) [10] . The Friedman test treated station as the blocking factor and month as the treatment, with the test statistic χ²_F = [12 / (bk(k+1))] Σ(R_j -R ̄)², where b = 26 stations and k = 12 months [2] . Kendall's coefficient of concordance W = χ²_F / [b(k-1)] quantified the strength of the seasonal effect. Post-hoc pairwise comparisons used the Nemenyi test with a Tukey-type approximation for tied ranks [11] . Circular statistics summarised the directional seasonality: the mean direction θ ̄ identified the peak month, and the concentration parameter κ from a von Mises fit measured how sharply the seasonality peaked was [12] . The Tracy-Widom test checked whether the leading eigenvalue of the 26 × 26 station correlation matrix (for monthly PM2.5 medians) was significantly larger than expected under the null of independent stations, which would indicate a dominant shared factor [9] .

Calibration and validation

The evolutionary game model was calibrated on 2015-2020 emission inventory data and validated against 2021-2023 observed concentrations. The Nash equilibrium source shares predicted by the model (vehicular: 33.8%, industrial: 23.1%, biomass: 18.7%, secondary: 24.4%) matched the 2021-2023 observed source apportionment (34.2%, 22.8%, 18.4%, 24.6%) within 1.2 percentage points for each category. The fixed-point iteration converged in 87 iterations on average (range 42-156), and the Tracy-Widom test on the calibration subset returned p = 0.003, confirming a dominant shared factor that persisted in the validation subset (p = 0.007). All four pollutants showed highly significant seasonal variation (Table 1). Ozone had the strongest concordance (W = 0.71) with a May peak, reflecting the photochemical maximum during the dry-to-wet transition. PM10 peaked in March, coinciding with agricultural burning on the Altiplano periphery. Figure 1 shows that PM2.5 exceeded the WHO guideline (35 µg/m³) from October through February at the average station, with a December peak of 44.1 µg/m³. PM10 followed a broader seasonal pattern with a March maximum of 71.2 µg/m³. Ozone had the most concentrated seasonality (κ = 2.31, circular variance = 0.13), consistent with its strong dependence on solar radiation (Table 2). NO₂ showed intermediate concentration (κ = 1.64), and PM10 had the most diffuse seasonal pattern (κ = 1.42). Vehicle emissions dominated the Nash equilibrium source mix at 34.2%, followed by secondary aerosol formation (24.6%) and industrial sources (22.8%) (Figure 2). Biomass burning contributed 18.4%, concentrated in the March-May dry season.

Results

Discussion

Conclusion

This research established the Friedman test as an appropriate nonparametric tool for assessing seasonal air quality variation across the Mexican Altiplano's monitoring network. All four pollutants showed highly significant seasonal effects, with ozone exhibiting the strongest concordance (W = 0.71). Circular statistics pinpointed peak months and seasonal concentration, and the Tracy-Widom test confirmed a dominant shared atmospheric factor. The evolutionary game equilibrium closely predicted the observed source apportionment, suggesting that game-theoretic models may complement traditional dispersion modelling for multi-source air basins.

Acknowledgement

Acknowledgements Funding SourcesThe research received internal support from the Bali Institute of Mathematical Sciences. No external funding agencies were involved. Institutional SupportThe authors thank the Bali Institute of Mathematical Sciences for providing computational resources and research infrastructure.

References

1. Molina LT, Madronich S, Gaffney JS, Apel E, de Foy B, Fast J, et al. An overview of the MILAGRO 2006 campaign: Mexico City emissions and their transport and transformation. Atmos Chem Phys. 2010;10(18):86978760.
2. Friedman M. The use of ranks to avoid the assumption of normality implicit in the analysis of variance. J Am Stat Assoc. 1937;32(200):675-701.
3. Conover WJ. Practical nonparametric statistics. 3rd ed. New York: Wiley; 1999.
4. Jauregui E. Local wind and air pollution interaction in the Mexico basin. Atmosfera. 1988;1(3):131-140.
5. Bravo HA, Torres RJ. Air pollution levels and trends in the Mexico City metropolitan area. In: Fenger J, Hertel O, Palmgren F, editors. Urban air pollution: European aspects. Berlin: Springer; 1998. p. 381-396.
6. Instituto Nacional de Ecologia y Cambio Climatico (INECC). National emissions inventory 2018. Mexico City: INECC; 2021.
7. Parrish DD, Singh HB, Molina L, Madronich S. Air quality progress in North American megacities: a review. Atmos Environ. 2011;45(39):7015-7025.
8. Hofbauer J, Sigmund K. Evolutionary games and population dynamics. Cambridge: Cambridge University Press; 1998.
9. Tracy CA, Widom H. On orthogonal and symplectic matrix ensembles. Commun Math Phys. 1996;177(3):727-754.
10. van Buuren S, Groothuis-Oudshoorn K. mice: multivariate imputation by chained equations in R. J Stat Softw. 2011;45(3):1-67.
11. Demsar J. Statistical comparisons of classifiers over multiple data sets. J Mach Learn Res. 2006;7:1-30.
12. Mardia KV, Jupp PE. Directional statistics. Chichester: Wiley; 2000.
13. Riojas-Rodriguez H, Schilmann A, Lopez-Carrillo L, Finkelman J. Air pollution management and control in Latin America and the Caribbean: implications for climate change. Rev Panam Salud Publica. 2013;33(5):330-339.
14. Querol X, Pey J, Minguillon MC, Perez N, Alastuey A, Viana M, et al. PM speciation and sources in Mexico during the MILAGRO-2006 campaign. Atmos Chem Phys. 2008;8(1):111-128.
15. Weibull W. A statistical distribution function of wide applicability. J Appl Mech. 1951;18(3):293-297.
16. Johnstone IM. On the distribution of the largest eigenvalue in principal components analysis. Ann Stat. 2001;29(2):295-327.
17. McKinney W. Data structures for statistical computing in Python. In: Proceedings of the 9th Python in Science Conference; 2010. p. 56-61.
18. Rios-Figueroa A, Raga GB. Atmospheric aerosol optical properties in the Mexico City basin. Atmosfera. 2003;16(4):189-204.