Henrique Higor

Cotta • Aranda

Anselmo Valdério

Reisen

Pascal Bondon

Paulo Roberto

Prezotti Filho

H Cotta

V Reisen

P Prezotti

Identification of redundant air quality monitoring stations using robust principal component analysis

Keywords: Air quality, Monitoring networks, time series, Robust principal component analysis, Outliers

Air quality monitoring stations are essentials for monitoring air pollutants and, therefore, are essential to protect the public health and the environment from the adverse effects of air pollution. Two or more stations may monitor the same pollutant behavior. In this scenario, the equipment must be reallocated to provide a better use of public resources and to enlarge the monitored area. The identification of redundant stations can be carried out by the application of principal component analysis (PCA) as a grouping technique. The principal component analysis

Introduction

The concern about air pollution problems has increased considerably in the last 50 years. Especially in developing countries, the air quality has been degraded as a result of industrialization, population growth, high rates of urbanization, and inadequate or nonexistent policies to control air pollution. The problems caused by air pollution produce local, regional, and global impacts. In this context, the particulate matter (PM), especially the PM 10 , which has an aerodynamic diameter less than 10 µm, is one of the most important pollutants with natural and anthropogenic sources.

Its adverse impacts on humans health may lead to an increment of mortality rates, respiratory and cardiovascular problems for short and longterm exposure at high concentrations [START_REF] Beelen | Effects of long-term exposure to air pollution on natural-cause mortality: an analysis of 22 european cohorts within the multicentre escape project[END_REF][START_REF] Cesaroni | Long term exposure to ambient air pollution and incidence of acute coronary events: prospective cohort study and meta-analysis in 11 european cohorts from the escape project[END_REF][START_REF] Hoek | Longterm air pollution exposure and cardio-respiratory mortality: a review[END_REF][START_REF] Rückerl | Health effects of particulate air pollution: a review of epidemiological evidence[END_REF].

The primary purpose of air quality management is to protect public health and the environment from the adverse effects of air pollution. Adequate control of air quality involves several activities such as risk management, setting standards for emissions and air quality, implementation of control measures, and risk communication (WHO 2005). The monitoring of air quality is essential for any air pollution control policy. The realization of efficient management of air quality is important for identifying and quantifying the pollutants found in a region and their sources. This is accomplished by using stations to monitor different pollutants according to the needs of the regions where the stations are installed.

In Brazil, although the limits for pollutant concentrations are established by the federal legislation CONAMA 003/90 (Conselho Nacional do Meio Ambiente 1990), this decree does not contemplate guidelines on how to construct or how to manage monitoring networks and, thus, entrusting this task to each one of the 27 federative units. In this scenario, an actual overview of Brazil's air quality monitoring networks is given in a recent publication of the Instituto de Energia e Meio Ambiente coauthored by the Brazilian Ministry of the Environment. This publication highlights although essential the air quality monitoring in Brazil is far from being a reality.

Due to the dimensions of the country, the non-prioritization of air quality policies and the amount of financial resources destined to the monitoring activities, only 12 out of 27 unity members have an operational air quality monitoring network (IEMA 2014).

The installation and continuous operation of an air quality monitoring station are cost-intensive as it requires finding a suitable place for the installation and personnel for its maintenance. Only one monitoring station should operate in an area characterized by a specific pattern of air pollution. Pires et al. (2008a) indicate the number of stations that constitute a monitoring network must be optimized to reduce costs and expenses. If there are stations with similar patterns of pollution for a specific pollutant, the monitoring equipment could be properly relocated to another area of interest.

In this context, the principal component analysis has been successfully used in air pollution for managing a network of monitoring stations in several studies, for instance, Zamprogno (2013) studied PCA with time series models in many different applications related to air pollution data, [START_REF] Zhao | A performance assessment and adjustment program for air quality monitoring networks in shanghai[END_REF] applied PCA to verify redundant air quality monitoring networks in Shanghai (China). [START_REF] Dominick | Spatial assessment of air quality patterns in malaysia using multivariate analysis[END_REF] used PCA and Cluster Analysis (CA) to check the pattern of behavior of the pollutants carbon monoxide (CO), ozone (O 3 ), particulate matter of diameter < 10µm (PM 10 ), sulfur dioxide (SO 2 ), nitric oxide (NO) and nitrogen dioxide (NO 2 ) in five different stations in Malaysia. [START_REF] Pires | Identification of redundant air quality measurements through the use of principal component analysis[END_REF]Pires et al. ( , 2008a,b) ,b) applied PCA to identify monitoring sites with similar concentrations of pollutants for PM 10 , SO 2 , CO, NO 2 and O 3 in the metropolitan area of Porto (Portugal). [START_REF] Lu | Performance assessment of air quality monitoring networks using principal component analysis and cluster analysis[END_REF] employed PCA to study the air quality monitoring network of Hong Kong for the pollutants of SO 2 , NO 2 and Respirable Suspended Particulate (RSP). The authors found that the monitoring stations located in nearby areas are characterized by the same specific air pollution characteristics and suggested that redundant equipment should be transferred to other monitoring stations allowing for further enlargement of the monitored area. Other studies include [START_REF] Hall Lau | Interpretation of air quality in relation to monitoring station's surroundings[END_REF] and [START_REF] Gramsch | Examination of pollution trends in santiago de chile with cluster analysis of PM 10 and ozone data[END_REF].

The application of PCA is not exclusive to the management of air quality monitoring networks. Recently, [START_REF] Villas-Boas | Assessment of the water quality monitoring network of the piabanha river experimental watersheds in rio de janeiro, brazil, using autoassociative neural networks[END_REF] used PCA and nonlinear PCA to assess the redundancy of the parameters and monitoring locations of the Piabanha water quality network in Brazil. [START_REF] Maronna | Temporal and spatial assessment of river surface water quality using multivariate statistical techniques: a study in can tho city, a mekong delta area, vietnam[END_REF] applied PCA and other multivariate statistical tools to assess the river surface water quality and also redundant monitoring stations in Can Tho City (Vietnam).

At this point, PCA is one of the main multivariate statistical techniques. The goal of PCA is to explain the covariance structure of the data through auxiliary variables called components. These components are constructed from linear combinations of the original variables and are uncorrelated. Briefly, PCA calculates the eigenvalues and eigenvectors of the covariance or correlation matrix. The main application of PCA is to reduce the dimensionality of a correlated data matrix of n dimension to a m dimension, where m < n. The reduction is performed so that the new set of variables captures most of the variability contained in the original data. A review of the fundamentals of PCA using R (R Core Team 2018) can be found in [START_REF] Sergeant | A practitioners guide for exploring water quality patterns using principal components analysis and procrustes[END_REF].

Besides the use for dimensionality reduction, the PCA technique can be used for the clustering of the variables of a data matrix. [START_REF] Cadima | Loading and correlations in the interpretation of principle compenents[END_REF] discusses the clustering of variables considering the eigenvectors of the PCA. The grouping of variables consists of choosing variables that have similar values for its eigenvectors in absolute value and are highly correlated to the principal component.

In the air pollution context, outliers may arise from different scenarios such as human-made disasters and natural catastrophes, measurement errors due to the fail-ure of equipment or a sudden change in the atmosphere conditions, human errors, among others. Another critical situation arises when the observed pollutant is under the concentration limits established by legislation standards, but it may be considered as an atypical observation during the statistical analysis.

Furthermore, PCA is sensitive to outliers since the estimation of the mean vector, the covariance matrix, and the correlation matrix are directly influenced by outliers.

As a consequence, the estimation of the eigenvalues and eigenvectors of the covariance or correlation matrix will be influenced by outliers present in the data, see, e.g., [START_REF] Filzmoser | Robust principal component and factor analysis in the geostatistical treatment of environmental data[END_REF]. It is worthwhile to mention that even a single outlier may affect classical statistics methods. [START_REF] Croux | Principal component analysis based on robust estimators of the covariance or correlation matrix: Influence functons and efficiencies[END_REF] indicate that conclusions obtained from principal component analysis calculated from a dataset with outliers may be misleading.

Under these circumstances, the common choice made by a wide range of scientists and practitioners to mitigate this problem is to delete the observations suspected to be outliers. As pointed out by Maronna et al. (2006, Chapter 1), the removal of an outlier observation may lead to many issues since the deletion is based on a subjective decision. A viable option to attenuate these problems is to use robust statistical methods since these methods still work well even when the presence of outliers is uncertain. Among the methods for robust estimation of the covariance or correlation matrix with time-independent datasets, there is the estimator proposed by [START_REF] Ma | Highly robust estimation of dispersion matrices[END_REF]. This estimator uses the so-called Q n (.) estimator proposed by [START_REF] Rousseeuw | Alternatives to the median absolute deviation[END_REF], which is independent of the location parameter of the dataset.

In this paper, the central idea is to robustify the estimation of the covariance matrix before calculating its eigenvalues and eigenvectors in PCA. The methodology proposed consists of the application of robust principal component analysis and selecting the stations which presented higher correlations to the selected PCs. Then, a decision rule is to be applied to decide to keep the redundant station in the same place or to move it to a new area. This proposed methodology is also adequate when outliers are presented in the dataset. The PM 10 data of the metropolitan area of the Greater Vitória Region (GVR), Brazil, is analyzed as an illustrative example.

The paper is structured as follows: Section 2 describes the data and the statistical model introducing the proposed estimation method and how to identify monitoring stations that present similar behavior; Section 3 presents the data analysis and its discussion comparing robust PCA to the standard one. Finally, Section 4 presents the closing remarks.

Data and methods

Sampling stations in the Greater Vitoria Region

The Greater Vitória Region is located on the southeast coast of Brazil (latitude 20ž 19S, longitude 40ž20W) with a population of approximately 1.900.000 inhabitants.

The climate is tropical humid, with average temperatures ranging from 24žC to 30žC. The region has many ports being an important cargo transport hub in Brazil. Also, there are many industries presented in the region, such as steel plants, iron ore pellet mill, stone quarrying, cement and food industry, and asphalt plant.

The automatic air pollution monitoring network (AAQMN) of GVR is consisted by eight monitoring stations distributed in the cities of this region as follows: two stations in Serra (Laranjeiras and Carapina), three stations in Vitória (Jardim Camburi, Enseada do Suá and Vitória Centro), two stations in Vila Velha (Vila Velha Centro and Ibes) and one station in Cariacica (at the regional food distribution center, CEASA). The PM 10 , in µg/m 3 , is monitored in all stations. Figure 1

Principal component analysis

Most of the practitioners employ the standard PCA, which is based on the sample covariance matrix and is summarized in the sequel. Let X 1 , . . . , X n be a sample of size n of an independent and identically distributed multivariate distribution with dimension p, mean vector µ, and covariance matrix Σ. The method of moment estimator (MME) of Σ is

Σn = 1 n n i=1 (X i -μ)(X i -μ) , (1)
where μ = 1 (2)

Then, we have

Var( Ŷi ) = â i Σn âi = λi , i = 1, 2, . . . , p, (3)
Cov( Ŷi , Ŷk ) = â i Σn âk = 0, i = k, i, k = 1, 2, . . . , p, . (4)
If some λi are equal, the choice of the corresponding eigenvectors âi is not unique.

Associated with (2), it can be shown that

p i=1 Var(X i ) = λ1 + λ2 + • • • + λp = p i=1
Var( Ŷi ).

(5)

Equation 5 states that the whole variability of X is retained by the principal components Ŷ . Therefore, if the main goal of the use of PCA is to reduce the number of variables, the scientist may choose to retain only part of the total original variability.

Robust PCA

Outliers affect the estimation of the location (mean) and the scale (variance) of random variables. To address this problem, [START_REF] Rousseeuw | Alternatives to the median absolute deviation[END_REF] proposed a robust estimator, Q n , for the dispersion of a dataset. Let X 1 , . . . , X n be n i.i.d.

copies of a random variable X, the estimator Q n is the k-th order

Q n (x) = d {|X i -X j |; i < j} {k} , (6)
where i, j = 1, . . . , n, and d is a value for consistency of the estimator. The k-th order statistic is the integer value k = ( n 2 + 2)/4 + 1.

It is known that for any univariate second-order random variables X and Y it is possible to compute the covariance between them as follows

Cov(X, Y ) = αβ 4 Var(X/α + Y /β) -Var(X/α -Y /β) , (7)
for any α, β ∈ R, see, [START_REF] Huber | Robust Statistics. Wiley Series in Probability and Statistics-Applied Probability and Statistics Section Series[END_REF]. In order to robustify (7), [START_REF] Ma | Highly robust estimation of dispersion matrices[END_REF] proposed to use the estimator Q n instead of the sample variance obtaining

σQ n (X, Y) = αβ 4 Q 2 n X α + Y β -Q 2 n X α - Y β , (8)
where α = Q n (X) and

β = Q n (Y ).
The correlation between the univariate second-order random variables X and Y can be estimated by

ρQ n (X, Y) = Q 2 n X α + Y β -Q 2 n X α -Y β Q 2 n X α + Y β + Q 2 n X α -Y β , (9)
where X, Y , α and β are defined in (8).

Let X be a random vector of p ≥ 2 variables. The robust sample covariance and correlation matrix of the random vector X, namely, ΣQ n and P Q n , respectively, are obtained by estimating every covariance or correlation pairs between X i and X j , i, j = 1, . . . , p. In this work, the robustified principal component analysis is achieved by replacing the standard covariance (or correlation matrix) with ΣQ n and P Q n .

It is worthwhile to mention that the robust estimation procedure discussed above will provide similar results to the ones estimated using the standard sample estimator when there are no outliers presented in the dataset. Therefore, its usage is recommended.

PCA clustering and station selection

PCA technique can also be used for clustering of the variables. A method for clustering variables using PCA is discussed in [START_REF] Cadima | Loading and correlations in the interpretation of principle compenents[END_REF]. The grouping of variables consists of choosing variables that have similar values for its eigenvectors in module and are highly correlated to the principal component. The correlation between a retained PC group and the related full PC (containing all the variability of Ŷi ) is given by

rk = λ1/2 j (â k j Σ-1 n,k âk j ) 1/2 , ( 10
)
where λj is eigenvalue of j-th component, âk j is the clustered vector of âj containing k variables and Σ-1 n,k is the sub-matrix of Σn , which lines and columns correspond to the k grouped variables.

The main idea behind the method is to address monitoring stations which present similar behaviors for the PM 10 pollutant (the technique is easily expanded to any other pollutant or meteorological parameter). Thus, a decision rule can be applied In this context, the following methodology for addressing monitoring stations which present similar behavior for a given pollutant is proposed:

1. Perform a descriptive statistical analysis of the data to verify the occurrence of possible outliers and to check for different scale of the measured variables;

2. Compute the robust PCA using the covariance or the correlation matrix;

Data analysis and discussion

In this study, the robust PCA was applied as a classification tool to group monitoring sites with redundant measurements of PM 10 concentrations from January 1st of 2005

to December 31 of 2009 (n = 1826). All the plots and analysis were performed using the computing environment R. ΣQ n and P Q n are available in the package tsqn (Cotta et al. 2017). The dataset and the R codes are available upon request. The boxplot of the data and the series of PM 10 are shown in Figures 2 and3, respectively. From the boxplot and the plots of the series, one can observe higher levels of PM 10 pollutant compared to WHO's guidelines, where the established limit is 50 µg/m 3 for 24-hour concentrations. Although the high levels of PM 10 are essential information that should be considered in the context of the air pollution and its impact on human health, these observations can be identified, from a statistical point of view, as being outliers. Therefore, the high levels of PM 10 presented in the series justify the use and comparison of the robust PCA.

Tables 2 and3 show the correlations and the robust correlations (as in Section 2.2) between the monitoring stations in the study. From both tables, we observe strong correlations between the variables, e.g., 0.78 for Ibes and Enseada do Suá stations.

The grouping of stations with redundant measurements for the PM 10 pollutant was carried out following the methodology proposed in Section 2.2.2. That is, stations having the same contribution in a given component will have similar values for their eigenvectors, and they will also be correlated to the component.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Laranjeiras Carapina J. Camburi Ens. do Suá Vitória Centro Ibes Vila V. Centro Cariacica
In the PCA tool, the estimates of the eigenvalues and their corresponding eigenvectors using P and P Q n are given in P . However, base on boxplot (Figure 2) and descriptive statistics (Table 1), the grouping based on P Q n is suggested here.

To visually confirm the grouping results for both estimators, the daily averages of PM 10 for the groups are shown in Figure 4. It is seen that the grouping using P Q n is superior since, for the first component, the grouped stations have similar concentrations.

To end this analysis and continuing with the procedure of the methodology discussed in Section 2.2.2, the stations of Vitória Centro and Enseada do Suá may be selected to be moved to a new area to enlarge the total monitored area. It is highlighted that although Cariacica has no important contribution to the robust cluster, it is the only station located in Cariacica municipality and, therefore, must be kept.

Conclusions

This article proposed and applied a grouping methodology to identify monitoring stations that present similar behavior for a given pollutant. As a case of study, the AAQMN of GVR (Brazil), which monitors the PM 10 pollutant was considered in order to enable better management of the local monitoring network.

The methodology proposed consists of the application of robust principal component analysis and selecting the stations which presented higher contributions to the chosen PCs. Then, a decision rule is to be applied to decide to keep the redundant station in the same place or to move it to a new area.

In the case study, it was found the occurrence of possible outliers observations during the descriptive analysis of the PM 10 data, which justified the comparison between the robust and standard PCA. It was found that Ibes, Enseada do Suá, and

Fig. 1

 1 Fig. 1 Geographical location of the stations.

  i . As stated byJollife (2002), the big drawback of PCA tool based on covariance matrices is the sensitivity of the PCs to the units of measurement of the variables. Therefore, if large differences in the variances of variables are found, the variables with large variances will tend to dominate the first PCs. To avoid this problem, the use of PCA based on the correlation matrix is suggested. To this end, the sample correlation matrix P can be obtained as P = D Σn D, where D = diag(1/ √ σ11 , . . . , 1/ σpp ), where σii , for i = 1, . . . , p, is the sample covariance. It is straightforward to see that even one outlier will affect the sample mean, and, thus, the whole covariance (or correlation matrix). Now, consider the random vector X = [X 1 , X 2 , . . . , X p ] with sample covariance matrix Σn and its associated sample eigenvalues λ1 ≥ λ2 ≥ . . . ≥ λp ≥ 0 with corresponding normed eigenvectors â = [â 1 , â2 , . . . , âp ]. Let Ŷi = â i X.

  to decide to keep the redundant station in the same place or to move it to a new area. As a possible decision rule, Pires et al. (2009) suggested three criteria: (i) sites should be monitoring the highest possible pollutant concentrations; (ii) the number of pollutants being monitored at each site should be maximized; and (iii) the distribution should maximize distances between locations.

3.

  Select a desirable number of PCs to be retained, e.g., 80% or more of the total variability; 4. Arbitrarily choose a cutoff point for the absolute values of the eigenvectors; 5. Create a group of variables whose coefficient of eigenvectors are equal or greater than the cutoff point in the component; 6. Using (10) compute the correlation between the selected variables in the PC and the full component. If the chosen variables and the component are not correlated, verify the cutoff point and redo the steps 4-6; 7. Apply the decision criteria of Pires et al. (2009) to decide to keep or to move to a new area the monitoring equipment of the pollutant considered in the study.

Fig. 2 Fig. 3

 23 Fig. 2 Boxplot of PM 10 's concentrations of the AAQMN of the GVR.

Table 1

 1 shows the descriptive statistics (i.e., the averages, standard deviations, and quantile values, among others) of the variables considered. The concentrations of PM 10 pollutants exceeded hourly and annually, the guidelines suggested by the World Health Organization (WHO 2005). It is observed a high range for all stations.

Table 1

 1 Descriptive statistics of PM

		Vila Velha Centro Cariacica	28.94 44.16	11.33 13.12	5.92 8.92	21.51 36.14	27.21 43.33	33.92 50.79	94.75 106.30
		Ibes	28.13	9.20	7.00	22.01	27.29	32.91	88.12
		Vitória Centro	26.01	7.37	5.62	21.46	25.25	29.78	70.42
		Jardim Camburi Enseada do Suá	28.97 28.08	8.01 8.29	8.67 7.50	23.64 22.71	28.33 27.00	33.46 32.46	78.08 74.58
	10 data	Laranjeiras Carapina	32.26 24.13	11.29 7.67	6.08 5.75	24.50 19.33	31.27 23.00	38.07 27.71	86.46 88.25
			Mean	Std. Dev	Min.	25th perc.	50th perc.	75th perc.	Max.

Table 3

 3 Robust correlation matrix ( P

		Vila Velha Centro Cariacica							1.00	0.51 1.00
		Vitória Centro Ibes					1.00	0.66 1.00	0.60 0.56	0.73 0.69
	Qn ) between the stations	Carapina Jardim Camburi Enseada do Suá		1.00	0.57 1.00	0.58 0.59 1.00	0.65 0.61 0.71	0.61 0.62 0.78	0.55 0.48 0.54	0.70 0.55 0.60
		Laranjeiras	1.00	0.40	0.59	0.59	0.45	0.66	0.44	0.46
			Laranjeiras	Carapina	Jardim Camburi	Enseada do Suá	Vitória Centro	Ibes	Vila Velha Centro	Cariacica

  Table4where, for each component, the grouped stations are highlight by boldface. For both estimators, four components could explain approximately 85 % of the total variability of the dataset, leading to a dimension reduction of the data. It is observed that PCA computed by using P Q n preserved a higher percentage of variability in the components.For both PCAs, the cutoff point was selected to be 0.37 in absolute value, which led to the highest correlation values. In the standard PCA, this cutoff led to a correlation between the selected PCs groups and the original PCs of 0.96, 0.88, 0.66, and 0.96, for the four PCs, respectively. In the case of robust PCA, correlations of 0.96, 0.89, 0.66, and 0.95 were found. The values are close in both standard and robust PCA.Thus, for the method of moments estimator for the first component, it is possible to visualize the existence of a group of stations formed by Ibes, Vila Velha Centro, and Cariacica. In the second component, the group is formed by Laranjeiras and Carapina. For the third component, Vila Velha Centro forms a group. Finally, the fourth component is the group formed by Jardim Camburi, and Enseada do Suá.For the grouping through robust PCA, in the first component, Ibes, Enseada do Suá, and Vitória Centro can be grouped. For the second component, Laranjeiras and Cariacica form a group. In the third component, Vila Velha Centro is the only station in the group. For the fourth component, the group is formed by Enseada do Suá and Jardim Camburi. Therefore, the proposed method allocated groups differently from

Table 4

 4 PCA results for PM 10 of AAQMN of the GVR Stations

		4	0.0683	0.2701	0.7037	-0.5144	-0.1276	-0.3395	0.1475	-0.0962	0.4612	5.77	86.14
	PCA -P Qn	2 3	0.6998 0.0533	-0.4144 -0.1961	0.2356 -0.2115	0.1519 -0.0045	-0.2867 -0.0211	0.1902 -0.0881	-0.1838 0.8942	-0.3283 -0.3259	0.7568 0.5334	9.46 6.67	73.71 80.38
		1	-0.3123	-0.3488	-0.3446	-0.3722	-0.3745	-0.3863	-0.3203	-0.3625	5.146	64.25	64.25
		4	0.1460	0.1750	0.7019	-0.6118	-0.1629	-0.2271	0.0141	0.0507	0.4973	6.22	84.97
	PCA -P	2 3	0.7193 -0.1756	-0.4004 0.2628	0.1700 0.0502	0.2163 0.0406	-0.2265 -0.1026	0.1787 0.2359	-0.2942 -0.8391	-0.2766 0.3542	0.7744 0.6282	9.68 7.85	70.90 78.75
		1	-0.3002	-0.3554	-0.3472	-0.3632	-0.3864	-0.3869	-0.3055	-0.3721	4.8971	61.22	61.22
			Laranjeiras	Carapina	Jardim Camburi	Enseada do Suá	Vitória Centro	Ibes	Vila Velha Centro	Cariacica	Eigenvalue	Proportion	Cumulative

Acknowledgements

The authors would like to thank the support from CNPq, ERASMUS, CAPES and FAPES. Part of this paper was revised when Professor Valdério Reisen was visiting CentraleSupélec in July 2018. This author is indebted to CentraleSupélec for its financial support.

PM 10 (µg/m 3 ) -average

Vitória Centro Ibes Cariacica Vila Velha Centro

0 10 20 30 40 50 60

Enseada do Suá Vitória Centro Ibes

Laranjeiras Carapina Vitória Centro presented similar behavior and thus can be grouped. Also, Jardim

Camburi and Enseada do Suá form another group. Therefore, two stations, Ibes and Enseada do Suá, are the candidates to be moved to a new site to enlarge the monitored area.


François Blanquart
email: francois.blanquart@college-de-france.fr

Clémence Abad

Joevin Ambroise

Mathieu Bernard

Gina Cosentino

Jean-Marc Giannoli

Florence Débarre
email: florence.debarre@normalesup.org

Characterization of vaccine-breakthrough infections of SARS-CoV-2 Delta and Alpha variants and within

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

individuals, Ct values at the first positive test were higher in fully vaccinated individuals (> 2 weeks after final dose) than non fully vaccinated individuals (1.7 [1, 2.3], p < 1e-6). In symptomatic individuals however, Ct values at the time of symptoms were not significantly different in vaccinated compared to unvaccinated individuals (p = 0.26). This was true both for infections by Delta and non-Delta (essentially Alpha in France at the time) variants. These results suggest that some infected vaccinated individuals, especially if symptomatic, may transmit the virus as much as unvaccinated individuals.

The SARS-CoV-2 variant of concern Delta, first detected in India, spread across the world in 2021, and in particular in Europe in late spring -early summer 2021, where it displaced the previously dominant Alpha variant. Delta was shown to be fitter than Alpha, [START_REF] Campbell | Increased transmissibility and global spread of SARS-CoV-2 variants of concern as at June 2021[END_REF][START_REF] Mlcochova | SARS-CoV-2 B.1.617.2 Delta variant emergence and vaccine breakthrough[END_REF][START_REF] Alizon | Rapid spread of the SARS-CoV-2 Delta variant in some French regions[END_REF] and may be associated with higher virulence [START_REF] Fisman | Progressive Increase in Virulence of Novel SARS-CoV-2 Variants in Ontario, Canada [Internet]. Infectious Diseases (except HIV/AIDS)[END_REF][START_REF] Ong | Clinical and Virological Features of SARS-CoV-2 Variants of Concern: A Retrospective Cohort Study Comparing B.1.1.7 (Alpha), B.1.315 (Beta)[END_REF][START_REF] Sheikh | SARS-CoV-2 Delta VOC in Scotland: Demographics, risk of hospital admission, and vaccine effectiveness[END_REF] and lesser vaccine effectiveness [START_REF] Sheikh | SARS-CoV-2 Delta VOC in Scotland: Demographics, risk of hospital admission, and vaccine effectiveness[END_REF][START_REF] Bernal | Effectiveness of Covid-19 Vaccines against the B.1.617.2 (Delta) Variant[END_REF] against symptomatic disease, especially after just one dose.

Delta has spread in countries with high vaccination levels, and breakthrough infections have been reported, with Ct values suggesting similar viral loads between vaccinated and unvaccinated infected individuals. [START_REF] Brown | Outbreak of SARS-CoV-2 Infections, Including COVID-19 Vaccine Breakthrough Infections, Associated with Large Public Gatherings -Barnstable County, Massachusetts[END_REF][START_REF] England | SARS-CoV-2 variants of concern and variants under investigation in England[END_REF][START_REF] Elliott | REACT-1 round 13 final report: Exponential growth, high prevalence of SARS-CoV-2 and vaccine effectiveness associated with Delta variant in England during May to July 2021[END_REF] A longitudinal study has confirmed similar Ct values between Delta-infected vaccinated and unvaccinated patients in the first week after diagnosis or symptom onset, with later faster decline in vaccinated patients. [START_REF] Chia | Virological and serological kinetics of SARS-CoV-2 Delta variant vaccine-breakthrough infections: A multi-center cohort study[END_REF] Comparisons of Ct values in infections with Delta compared to infections with previous variants require controlling for infection age when variants has different epidemiological dynamics. This is because viral load depends on infection age, and the distribution of infection ages depends on whether the number of cases is growing or shrinking. [START_REF] Alizon | Epidemiological and clinical insights from SARS-CoV-2 RT-PCR cycle amplification values [Internet[END_REF][START_REF] Hay | Estimating epidemiologic dynamics from cross-sectional viral load distributions[END_REF][START_REF] Cosentino | SARS-CoV-2 viral dynamics in infections with Alpha and Beta variants of concern in the French community[END_REF][START_REF] Althaus | A tale of two variants: Spread of SARS-CoV-2 variants Alpha in Geneva, Switzerland, and Beta in South Africa [Internet[END_REF][START_REF] Hay | Viral loads observed under competing strain dynamics[END_REF] We studied the determinants of Ct values at the time of test and, for symptomatic individuals, as a function of the time since symptoms, in data from 292284 patients tested by a large private laboratory in the community in France, from 14 June 2021 to 30 July 2021 in three regions (Bretagne, Île-de-France, Provence-Alpes-Côte d'Azur). These data include information on the result of the PCR test, the associated Ct value, the patient's self-reported vaccine status (whether fully vaccinated since at least two weeks, or not), whether the patient has been symptomatic and the time since the onset of symptoms. Positive tests were screened for the L452R mutation, which characterizes the Delta variant (9343 positive tests with mutation information). In the case of multiple tests per individual, we kept the last negative test if there were no positive tests, and the first positive test otherwise.

Consistent with the French vaccination campaign, vaccinated individuals are on average older than non-vaccinated patients in our dataset (12 years older). The proportion of vaccinated individuals in the dataset (24%) is lower than in the community (47.5% by 10 July 2021), reflecting the fact that the data are not surveillance-based.

Reasons for testing may vary between vaccinated and non-vaccinated individuals. This may especially be the case since France introduced a "sanitary passport," requiring a proof of either full vaccination or a negative test for specific events, which may artificially inflate the proportion of negative tests among non-vaccinated individuals. Conversely, vaccinated individuals may get tested only if they have good reasons to suspect an infection; these reasons may also vary if symptoms differ depending on the infecting variant.

We compared the cycle threshold (Ct; targeted at gene RdRp) values of the PCR of positive tests depending on vaccine status, the presence of symptoms, and the infecting variant (Delta: presence of the L452R mutation), for the 8437 individuals for which all pieces of information are available. The Ct is the number of PCR cycles needed to detect a target; it is negatively correlated with viral load.

We find that the presence of symptoms is associated with significantly lower Ct values (ANOVA; -2.7

Ct [-3, -2.5]; adjusted p < 1e-6 (Tukey's Honest Significant Difference)). An infection with the Delta variant is also associated with lower Ct values in these data (-6.7 Ct [-7.1, -6.3]; adjusted p < 1e-6)); note that age of infection is not controlled for here, but will be later on. Comparing vaccinated and non-vaccinated individuals, we find that vaccinated individuals have significantly higher Ct values for both Delta and non-Delta asymptomatic infections (1.7 [1, 2.3] Ct difference; adjusted p < 1e-6), but that the differences are not significant for symptomatic infections (p = 0.8). For non-Delta variants, this result may be due to too small a sample size (only N = 18 vaccinated symptomatic individuals infected with non-Delta).

Ct values also depend on the age of infection of tested individuals. It is therefore useful to control for age of infection, especially when comparing Ct values of variants with different epidemiological dynamics, [START_REF] Hay | Estimating epidemiologic dynamics from cross-sectional viral load distributions[END_REF][START_REF] Hay | Viral loads observed under competing strain dynamics[END_REF] as was the case for the Delta (increasing numbers of infections) and Alpha (decreasing numbers of infections) in the early summer in France. To this end, we add time since symptoms onset as a continuous variable in the linear model. We find that the Delta variant has Ct at day of symptoms -3.32 [-4.38, -2.25] lower than non-Delta (Alpha) (p < 1e-6, N = 3439). The slope of Ct as function of time is 0.6 [0.54, 0.66] per day for Delta and 0.92 [0.73, 1.1] for non-Delta variants (p < 1e-6). Vaccine status does not significantly alter the outcomes, whether regarding the Ct at symptom onset (p = 0.256) or the slope of Ct as function of time since symptoms (p = 0.947) and was therefore not included in the final model.

Limitations of our study stem from the way the data were collected: this is a community-based study. Reasons for seeking a PCR test are unknown and may vary among individuals and across time.

Symptom and vaccine information are self-reported. Yet our dataset is unique for France, because variant information and vaccine status data have not been linked yet in public datasets, and Ct values are not reported at the national level (only positive/negative test results are). Our results are in line with a retrospective cohort study which found lower Ct values with Delta and longer duration of infection with low Ct. [START_REF] Ong | Clinical and Virological Features of SARS-CoV-2 Variants of Concern: A Retrospective Cohort Study Comparing B.1.1.7 (Alpha), B.1.315 (Beta)[END_REF] Regarding vaccine-breakthroughs, our results confirm studies finding similar Ct values among fully vaccinated individuals and those who were not, with the majority of infections being due to Delta. [START_REF] Brown | Outbreak of SARS-CoV-2 Infections, Including COVID-19 Vaccine Breakthrough Infections, Associated with Large Public Gatherings -Barnstable County, Massachusetts[END_REF][START_REF] England | SARS-CoV-2 variants of concern and variants under investigation in England[END_REF] Another limitation of our dataset is the lack of longitudinal data (in Figure 2, each point corresponds to a single patient). A recent study in Singapore found similar Ct values among vaccinated and unvaccinated individuals infected by Delta, at the time of diagnosis or of symptom onset. [START_REF] Chia | Virological and serological kinetics of SARS-CoV-2 Delta variant vaccine-breakthrough infections: A multi-center cohort study[END_REF] After a week however, Ct values increased faster (i.e. viral load declined faster) among vaccinated individuals than unvaccinated individuals, even after excluding asymptomatic individuals (personal communication). [START_REF] Chia | Virological and serological kinetics of SARS-CoV-2 Delta variant vaccine-breakthrough infections: A multi-center cohort study[END_REF] With transmission occurring early in infection, this delayed differential decrease is expected not to have much of an impact on relative transmissibility.

Ct values are linked to viral load, and viral load has been shown to be positively associated with probability of transmission in household contacts -but there is also high inter-individual variation, and transmission from individuals with high Ct values is possible. [START_REF] Marc | Quantifying the relationship between SARS-CoV-2 viral load and infectiousness[END_REF][START_REF] Lyngse | Association between SARS-CoV-2 Transmissibility, Viral Load, and Age in Households [Internet[END_REF] In spite of these limitations, our confirmatory results indicate that epidemic control may require similar measures for symptomatic PCR-positive vaccinated individuals as non-vaccinated infected individuals, and that it is important to not to stop testing vaccinated individuals.

  Figures

Figure 1 :Figure 2 :

 12 Figure 1: Distributions of Ct values, according to vaccine status (light: unvaccinated, dark: vaccinated), infecting variant (orange: non-Delta, purple: Delta) and whether the individual was symptomatic (left: asymptomatic at least until the test, right: symptomatic). The widths of the violin plots reflect the number of tests; points and arrows represent predicted values and confidence intervals. "***" means that the corresponding comparison of unvaccinated and vaccinated individuals is statistically significant with p < 0.001.

##

  . codes: 0 *** 0.001 ** 0.01 * 0.05 . 0: aov(formula = Ct ~vacf * sympf * deltaf, data = tmp) 1-0:0 -6.5769166 -7.1265213 -6.027312 0.0000000 ## 1:1-0:0 -5.9746925 -6.6780558 -5.271329 0.0000000 :1:0-1:0:0 -9.1605078 -13.5534270 -4.76758868 0.0000000 ## 0:0:1-1:0:0 -11.0393880 -13.1792936 -8.89948235 0.0000000 ## 1:0:1-1:0:0 -9.6154935 -11.8746836 -7.35630334 0.0000000

Acknowledgements

We thank Mircea Sofonea for critical comments on a previous version of the manuscript.

FB was funded by a Momentum grant from the Centre National de la Recherche Scientifique. FD was funded by grant ANR-19-CE45-0009-01 from Agence Nationale de la Recherche.

Model for