Evaluation of linear regression techniques for atmospheric applications: the importance of appropriate weighting

Linear regression techniques are widely used in atmospheric science, but they are often improperly applied due to lack of consideration or inappropriate handling of measurement uncertainty. In this work, numerical experiments are performed to evaluate the performance of five linear regression techniques, significantly extending previous works by Chu and Saylor. The five techniques are ordinary least squares (OLS), Deming regression (DR), orthogonal distance regression (ODR), weighted ODR (WODR), and York regression (YR). We first introduce a new data generation scheme that employs the Mersenne twister (MT) pseudorandom number generator. The numerical simulations are also improved by (a) refining the parameterization of nonlinear measurement uncertainties, (b) inclusion of a linear measurement uncertainty, and (c) inclusion of WODR for comparison. Results show that DR, WODR and YR produce an accurate slope, but the intercept by WODR and YR is overestimated and the degree of bias is more pronounced with a low R² XY dataset. The importance of a properly weighting parameter λ in DR is investigated by sensitivity tests, and it is found that an improper λ in DR can lead to a bias in both the slope and intercept estimation. Because the λ calculation depends on the actual form of the measurement error, it is essential to determine the exact form of measurement error in the XY data during the measurement stage. If a priori error in one of the variables is unknown, or the measurement error described cannot be trusted, DR, WODR and YR can provide the least biases in slope and intercept among all tested regression techniques. For these reasons, DR, WODR and YR are recommended for atmospheric studies when both X and Y data have measurement errors. An Igor Pro-based program (Scatter Plot) was developed to facilitate the implementation of error-in-variables regressions.

1 Introduction

Linear regression is heavily used in atmospheric science to derive the slope and intercept of XY datasets. Examples of linear regression applications include primary OC (organic carbon) and EC (elemental carbon) ratio estimation (Turpin and Huntzicker, 1995; Lin et al., 2009), MAE (mass absorption efficiency) estimation from light absorption and EC mass (Moosmüller et al., 1998), source apportionment of polycyclic aromatic hydrocarbons using CO and NO_x as combustion tracers (Lim et al., 1999), gas-phase reaction rate determination (Brauers and Finlayson-Pitts, 1997), inter-instrument comparison (Bauer et al., 2009; Cross et al., 2010; von Bobrutzki et al., 2010; Zieger et al., 2011; Wu et al., 2012; Huang et al., 2014; Zhou et al., 2016), inter-species analysis (Yu et al., 2005; Kuang et al., 2015), analytical protocol comparison (Chow et al., 2001, 2004; Cheng et al., 2011; Wu et al., 2016), light extinction budget reconstruction (Malm et al., 1994; Watson, 2002; Li et al., 2017), comparison between modeling and measurement (Petäjä et al., 2009), emission factor study (Janhäll et al., 2010), retrieval of shortwave cloud forcing (Cess et al., 1995), calculation of pollutant growth rate (Richter et al., 2005), estimation of ground-level PM_2.5 from MODIS data (Wang and Christopher, 2003), distinguishing OC origin from biomass burning using K⁺ as a tracer (Duan et al., 2004) and emission type identification by the EC ∕ CO ratio (Chen et al., 2001).

Ordinary least squares (OLS) regression is the most widely used method due to its simplicity. In OLS, it is assumed that independent variables are error-free. This is the case for certain applications, such as determining a calibration curve of an instrument in analytical chemistry. For example, a known amount of analyte (e.g., through weighing) can be used to calibrate the instrument output response (e.g., voltage). However, in many other applications, such as inter-instrument comparison, X and Y (from two instruments) may have comparable degrees of uncertainty. This deviation from the underlying assumption in OLS would produce biased slope and intercept when OLS is applied to the dataset.

To overcome the drawback of OLS, a number of error-in-variables regression models (also known as bivariate fittings; Cantrell, 2008) or total least-squares methods (Markovsky and Van Huffel, 2007) arise. Deming (1943) proposed an approach by minimizing sum of squares of X and Y residuals. A closed-form solution of Deming regression (DR) was provided by York (1966). Method comparison work of various regression techniques by Cornbleet and Gochman (1979) found significant error in OLS slope estimation when the relative standard deviation (RSD) of measurement error in “X” exceeded 20 %, while DR was found to reach a more accurate slope estimation. In an early application of the EC tracer method, Turpin and Huntzicker (1995) realized the limitation of OLS since OC and EC have comparable measurement uncertainty and thus recommended the use of DR for (OC ∕ EC)_pri (primary OC to EC ratio) estimation. Ayers (2001) conducted a simple numerical experiment and concluded that reduced major axis regression (RMA) is more suitable for air quality data regression analysis. Linnet (1999) pointed out that when applying DR for inter-method (or inter-instrument) comparison, special attention should be paid to the sample size. If the range ratio (max ∕ min) is relatively small (e.g., less than 2), more samples are needed to obtain statistically significant results.

In principle, a best-fit regression line should have greater dependence on the more precise data points rather than the less reliable ones. Chu (2005) performed a comparison study of OLS and DR specifically focusing on the EC tracer method application and found that the slope estimated by DR is closer to the correct value than OLS but may still overestimate the ideal value. Saylor et al. (2006) extended the comparison work of Chu (2005) by including a regression technique developed by York et al. (2004). They found that the slope overestimation by DR in the study of Chu (2005) was due to improper configuration of the weighting parameter, λ. This λ value is the key to handling the uneven errors between data points for the best-fit line calculation. This example demonstrates the importance of appropriate weighting in the calculation of best-fit line for error-in-variables regression model, which is overlooked in many studies.

In this study, we extend the work by Saylor et al. (2006) to achieve four objectives. The first is to propose a new data generation scheme by applying the Mersenne twister (MT) pseudorandom number generator for evaluation of linear regression techniques. In the study of Chu (2005), data generation is achieved by a varietal sine function, which has limitations in sample size, sample distribution, and nonadjustable correlation (R²) between X and Y. In comparison, the MT data generation provides more flexibility, permitting adjustable sample size, XY correlation and distribution. The second is to develop a nonlinear measurement error parameterization scheme for use in the regression method. The third is to incorporate linear measurement errors in the regression methods. In the work by Chu (2005) and Saylor et al. (2006), the relative measurement uncertainty (γ_Unc) is nonlinear with concentration, but a constant γ_Unc is often applied on atmospheric instruments due to its simplicity. The fourth is to include weighted orthogonal distance regression (WODR) for comparison. Abbreviations and symbols used in this study are summarized in Table B1 for quick reference.

2 Description of regression techniques compared in this study

Ordinary least squares (OLS) method

OLS only considers the errors in dependent variables (Y). OLS regression is achieved by minimizing the sum of squares (S) in the Y residuals (i.e., distance of AB in Fig. S1 in the Supplement):

$$\begin{array}{}\text{(1)}& S=\sum _{i=\mathrm{1}}^N(y_i-Y_i{)}^{\mathrm{2}},\end{array}$$

where Y_i are observed Y data points, while y_i are regressed Y data points of the regression line. N represents the number of data points that is used for regression.

Orthogonal distance regression (ODR)

ODR minimizes the sum of the squared orthogonal distances from all data points to the regressed line and considers equal error variances (i.e., distance of AC in Fig. S1):

$$\begin{array}{}\text{(2)}& S=\sum _{i=\mathrm{1}}^N\left[{\left(x_i-X_i\right)}^{\mathrm{2}}+{\left(y_i-Y_i\right)}^{\mathrm{2}}\right].\end{array}$$

Weighted orthogonal distance regression (WODR)

Unlike ODR, which considers even error in X and Y, weightings based on measurement errors in both X and Y are considered in WODR when minimizing the sum of squared orthogonal distance from the data points to the regression line (Carroll and Ruppert, 1996) as shown by AD in Fig. S1:

$$\begin{array}{}\text{(3)}& S={\sum }_{i=\mathrm{1}}^N\left[{\left(x_i-X_i\right)}^{\mathrm{2}}+{\left(y_i-Y_i\right)}^{\mathrm{2}}/\mathit{\eta }\right],\end{array}$$

where η is the error variance ratio, which determines the angle θ shown in Fig. S1. Implementation of ODR and WODR in Igor Pro (WaveMetrics, Inc. Lake Oswego, OR, USA) was done by the computer routine ODRPACK95 (Boggs et al., 1989; Zwolak et al., 2007).

Deming regression (DR)

Deming (1943) proposed the following function to minimize both the X and Y residuals as shown by AD in Fig. S1,

$$\begin{array}{}\text{(4)}& S={\sum }_{i=\mathrm{1}}^N\left[\mathit{\omega }\left(X_i\right)(x_i-X_i{)}^{\mathrm{2}}+\mathit{\omega }\left(Y_i\right)(y_i-Y_i{)}^{\mathrm{2}}\right],\end{array}$$

where X_i and Y_i are observed data points and x_i and y_i are regressed data points. Individual data points are weighted based on errors in X_i and Y_i,

$$\begin{array}{}\text{(5)}& \mathit{\omega }\left(X_i\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{X_i}^{\mathrm{2}}}},\mathit{\omega }\left(Y_i\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{Y_i}^{\mathrm{2}}}},\end{array}$$

where ${\mathit{\sigma }}_{X_i}$ and ${\mathit{\sigma }}_{Y_i}$ are the standard deviation of the error in measurement of X_i and Y_i, respectively. The closed-form solutions for slope and intercept of DR are shown in Appendix A.

York regression (YR)

The York method (York et al., 2004) introduces the correlation coefficient of errors in X and Y into the minimization function.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }S={\sum }_{i=\mathrm{1}}^N\left[\mathit{\omega }\left(X_i\right)(x_i-X_i{)}^{\mathrm{2}}-\mathrm{2}{r}_i\sqrt{\mathit{\omega }\left(X_i\right)\mathit{\omega }\left(Y_i\right)}\right.\\ \text{(6)}& {\displaystyle }& {\displaystyle }\left.(x_i-X_i)(y_i-Y_i)+\mathit{\omega }\left(Y_i\right)(y_i-Y_i{)}^{\mathrm{2}}\right]{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-r_i^{\mathrm{2}}}},\end{array}$$

where r_i is the correlation coefficient between measurement errors in X_i and Y_i. The slope and intercept of YR are calculated iteratively through the formulas in Appendix A.

Summary of the five regression techniques is given in Table S1 in the Supplement. It is worth noting that OLS and DR have closed-form expressions for calculating slope and intercept. In contrast, ODR, WODR and YR need to be solved iteratively. This need to be taken into consideration when choosing regression algorithm for handling huge numbers of data.

A computer program (Scatter Plot; Wu, 2017a) with a graphical user interface (GUI) in Igor Pro (WaveMetrics, Inc. Lake Oswego, OR, USA) was developed to facilitate the implementation of error-in-variables regression (including DR, WODR and YR). Two other Igor Pro-based computer programs, Histbox (Wu, 2017b) and Aethalometer data processor (Wu, 2017c), are used for data analysis and visualization in this study.

3 Data description

Two types of data are used for regression comparison. The first type is synthetic data generated by computer programs, which can be used in the EC tracer method (Turpin and Huntzicker, 1995) to demonstrate the regression application. The true “slope” and “intercept” are assigned during data generation, allowing quantitative comparison of the bias of each regression scheme. The second type of data comes from ambient measurement of light absorption, OC and EC in Guangzhou for demonstration in a real-world application.

3.1 Synthetic XY data generation

In this study, numerical simulations are conducted in Igor Pro (WaveMetrics, Inc. Lake Oswego, OR, USA) through custom codes. Two types of generation schemes are employed: one is based on the MT pseudorandom number generator (Matsumoto and Nishimura, 1998) and the other is based on the sine function described by Chu (2005).

The general form of linear regression on XY data can be written as

$$\begin{array}{}\text{(7)}& Y=kX+b,\end{array}$$

where k is the regressed slope and b is the intercept. The underlying meaning is that, Y can be decomposed into two parts. One part is correlated with X, and the ratio is defined by k. The other part of Y is constant and independent of X and regarded as b.

To make the discussion easier to follow, we intentionally avoid discussion using the abstract general form and instead opt to use a real-world application case in atmospheric science. Linear regression had been heavily applied on OC and EC data, here we use OC and EC data as an example to demonstrate the regression application in atmospheric science. In the EC tracer method, OC (mixture) is Y and EC (tracer) is X. OC can be decomposed into three components based on their formation pathway:

$$\begin{array}{}\text{(8)}& \mathrm{OC}={\mathrm{POC}}_{\mathrm{comb}}+{\mathrm{POC}}_{\text{non-comb}}+\mathrm{SOC},\end{array}$$

where POC_comb is primary OC from combustion. POC_non-comb is primary OC emitted from non-combustion activities. SOC is secondary OC formed during atmospheric aging. Since POC_comb is co-emitted with EC and well correlated with each other, their relationship can be parameterized as

$$\begin{array}{}\text{(9)}& {\mathrm{POC}}_{\mathrm{comb}}=(\mathrm{OC}/\mathrm{EC}{)}_{\mathrm{pri}}\times \mathrm{EC}.\end{array}$$

By carefully selecting an OC and EC subset when SOC is very low (considered as approximately zero), the combination of Eqs. (8) and (9) becomes

$$\begin{array}{}\text{(10)}& \mathrm{POC}=(\mathrm{OC}/\mathrm{EC}{)}_{\mathrm{pri}}\times \mathrm{EC}+{\mathrm{POC}}_{\text{non-comb}}.\end{array}$$

The regressed slope of POC (Y) against EC (X) represents (OC ∕ EC)_pri (k in Eq. 7). The regressed intercept become POC_non-comb (b in Eq. 7). With known (OC ∕ EC)_pri and POC_non-comb, SOC can be estimated by

$$\begin{array}{}\text{(11)}& \mathrm{SOC}=\mathrm{OC}-\left(\right(\mathrm{OC}/\mathrm{EC}{)}_{\mathrm{pri}}\times \mathrm{EC}+{\mathrm{POC}}_{\text{non-comb}}).\end{array}$$

The data generation starts from EC (X values). Once EC is generated, POC_comb (the part of Y that is correlated with X) can be obtained by multiplying EC by a preset constant, (OC ∕ EC)_pri (slope k). Then the other preset constant POC_non-comb is added to POC_comb and the sum becomes POC (Y values). To simulate the real-world situation, measurement errors are added on X and Y values. Details of synthesized measurement error are discussed in the next section. Implementation of data generation by two types of mathematical schemes is explained in Sect. 3.1.2 and 3.1.3, respectively.

3.1.1 Parameterization of synthesized measurement uncertainty

Weighting of variables is a crucial input for errors-in-variables linear regression methods such as DR, YR and WODR. In practice, the weights are usually defined as the inverse of the measurement error variance (Eq. 5). When measurement errors are considered, measured concentrations (Conc._measured) are simulated by adding measurement uncertainties (ε_Conc.) to the true concentrations (Conc._true):

$$\begin{array}{}\text{(12)}& \mathrm{Conc}{.}_{\mathrm{measured}}=\mathrm{Conc}{.}_{\mathrm{true}}+{\mathit{\epsilon }}_{\mathrm{Conc}.},\end{array}$$

where ε_Conc. is the random error following an even distribution with an average of 0, the range of which is constrained by

$$\begin{array}{}\text{(13)}& -{\mathit{\gamma }}_{\mathrm{Unc}}\times \mathrm{Conc}{.}_{\mathrm{true}}\le {\mathit{\epsilon }}_{\mathrm{Conc}.}\le +{\mathit{\gamma }}_{\mathrm{Unc}}\times \mathrm{Conc}{.}_{\mathrm{true}}.\end{array}$$

The γ_Unc is a dimensionless factor that describes the fractional measurement uncertainty relative to the true concentration (Conc._true). γ_Unc could be a function of Conc._true (Thompson, 1988) or a constant. The term ${\mathit{\gamma }}_{\mathrm{Unc}}\times \mathrm{Conc}{.}_{\mathrm{true}}$ defines the boundary of random measurement errors.

Two types of measurement error are considered in this study. The first type is γ_{Unc−nonlinear}. In the data generation scheme of Chu (2005) for the measurement uncertainties (ε_POC and ε_EC), γ_{Unc−nonlinear} is nonlinearly related to Conc._true:

$$\begin{array}{}\text{(14)}& {\mathit{\gamma }}_{\mathrm{Unc}-\mathrm{nonlinear}}={\displaystyle \frac{\mathrm{1}}{{\sqrt{\mathrm{Conc}.}}_{\mathrm{true}}}},\end{array}$$

and thus Eq. (13) for POC and EC becomes

$$\begin{array}{}\text{(15)}& {\displaystyle }-& {\displaystyle \frac{\mathrm{1}}{{\sqrt{\mathrm{POC}}}_{\mathrm{true}}}}\times {\mathrm{POC}}_{\mathrm{true}}\le {\mathit{\epsilon }}_{\mathrm{POC}}\le +{\displaystyle \frac{\mathrm{1}}{{\sqrt{\mathrm{POC}}}_{\mathrm{true}}}}\times {\mathrm{POC}}_{\mathrm{true}},\text{(16)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{{\sqrt{\mathrm{EC}}}_{\mathrm{true}}}}\times {\mathrm{EC}}_{\mathrm{true}}\le {\mathit{\epsilon }}_{\mathrm{EC}}\le +{\displaystyle \frac{\mathrm{1}}{{\sqrt{\mathrm{EC}}}_{\mathrm{true}}}}\times {\mathrm{EC}}_{\mathrm{true}}.\end{array}$$

In Eq. (14), the γ_Unc decreases as concentration increases, since low concentrations are usually more challenging to measure. As a result, the γ_{Unc−nonlinear} defined in Eq. (14) is more realistic than the constant approach, but there are two limitations. First, the physical meaning of the uncertainty unit is lost. If the unit of OC is µg m⁻³, then the unit of ε_OC becomes $\sqrt{{\mathrm{\mu }\mathrm{g}\mathrm{m}}^{-\mathrm{3}}}$. Second, the concentration is not normalized by a consistent relative value, making it sensitive to the X and Y units used. For example, if POC_true= 0.9 µg m⁻³, then ε_POC= ± 0.95 µg m⁻³ and γ_Unc= 105 %, but by changing the concentration unit to POC_true= 900 ng m⁻³, ε_OC= ± 30 ng m⁻³ and γ_Unc= 3 %. To overcome these deficiencies, we propose to modify Eq. (14) to

$$\begin{array}{}\text{(17)}& {\mathit{\gamma }}_{\mathrm{Unc}}=\sqrt{{\displaystyle \frac{\mathrm{LOD}}{\mathrm{Conc}{.}_{\mathrm{true}}}}}\times \mathit{\alpha },\end{array}$$

where LOD (limit of detection) is introduced to generate a dimensionless γ_Unc. α is a dimensionless adjustable factor to control the position of γ_Unc curve on the concentration axis, which is indicated by the value of γ_Unc at LOD level. As shown in Fig. 1a, at different values of α (α= 1, 0.5 and 0.3), the corresponding γ_Unc at the same LOD level would be 100, 50 and 30 %, respectively. By changing α, the location of the γ_Unc curve on x axis direction can be set, using the γ_Unc at LOD as the reference point. Then Eq. (7) for POC and EC becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }-\sqrt{{\displaystyle \frac{{\mathrm{LOD}}_{\mathrm{POC}}}{{\mathrm{POC}}_{\mathrm{true}}}}}\times {\mathit{\alpha }}_{\mathrm{POC}}\times {\mathrm{POC}}_{\mathrm{true}}\le {\mathit{\epsilon }}_{\mathrm{POC}}\le \\ \text{(18)}& {\displaystyle }& {\displaystyle }+\sqrt{{\displaystyle \frac{{\mathrm{LOD}}_{\mathrm{POC}}}{{\mathrm{POC}}_{\mathrm{true}}}}}\times {\mathit{\alpha }}_{\mathrm{POC}}\times {\mathrm{POC}}_{\mathrm{true}},{\displaystyle }& {\displaystyle }-\sqrt{{\displaystyle \frac{{\mathrm{LOD}}_{\mathrm{EC}}}{{\mathrm{EC}}_{\mathrm{true}}}}}\times {\mathit{\alpha }}_{\mathrm{EC}}\times {\mathrm{EC}}_{\mathrm{true}}\le {\mathit{\epsilon }}_{\mathrm{EC}}\le \\ \text{(19)}& {\displaystyle }& {\displaystyle }+\sqrt{{\displaystyle \frac{{\mathrm{LOD}}_{\mathrm{EC}}}{{\mathrm{EC}}_{\mathrm{true}}}}}\times {\mathit{\alpha }}_{\mathrm{EC}}\times {\mathrm{EC}}_{\mathrm{true}}.\end{array}$$

With the modified γ_{Unc−nonlinear} parameterization, concentrations of POC and EC are normalized by a corresponding LOD, which maintains unit consistency between POC_true and ε_POC and EC_true and ε_EC and eliminates dependency on the concentration unit.

Figure 1 (a) Example γ_{Unc−nonlinear} curves by different α values (Eq. 17). The x axis is concentration (normalized by LOD) in log scale and the y axis is γ_Unc. Black, blue and green line represents α equal to 1, 0.5 and 0.3, respectively, corresponding to the γ_{Unc−nonlinear} at LOD level equals to 100, 50 and 30 %, respectively. The red line represents γ_Unc−linear of 10 %. (b) Example of measurement uncertainty generation of γ_{Unc−nonlinear} and γ_Unc−linear. The blue circles represent γ_{Unc−nonlinear} following Eq. (17) (LOD_EC=1, a_EC=1). The red circles represent γ_Unc−linear (30 %).

Uniform distribution has been used in previous studies (Cox et al., 2003; Chu, 2005; Saylor et al., 2006) and is adopted in this study to parameterize measurement error. For a uniform distribution in the interval [a,b], the variance is $\frac{\mathrm{1}}{\mathrm{12}}{\left(a-b\right)}^{\mathrm{2}}$. Since ε_POC and ε_EC follow a uniform distribution in the interval as given by Eqs. (18) and (19), the weights in DR and YR (inverse of variance) become

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }\mathit{\omega }\left(X_i\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{X_i}^{\mathrm{2}}}}={\displaystyle \frac{\mathrm{3}}{{\mathrm{EC}}_{\mathrm{true}}\times {\mathrm{LOD}}_{\mathrm{EC}}\times {\mathit{\alpha }}_{\mathrm{EC}}^{\mathrm{2}}}},\text{(21)}& {\displaystyle }& {\displaystyle }\mathit{\omega }\left(Y_i\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{Y_i}^{\mathrm{2}}}}={\displaystyle \frac{\mathrm{3}}{{\mathrm{POC}}_{\mathrm{true}}\times {\mathrm{LOD}}_{\mathrm{POC}}\times {\mathit{\alpha }}_{\mathrm{POC}}^{\mathrm{2}}}}.\end{array}$$

The parameter λ in Deming regression is then determined:

$$\begin{array}{}\text{(22)}& \mathit{\lambda }={\displaystyle \frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}}={\displaystyle \frac{{\mathrm{POC}}_{\mathrm{true}}\times {\mathrm{LOD}}_{\mathrm{POC}}\times {\mathit{\alpha }}_{\mathrm{POC}}^{\mathrm{2}}}{{\mathrm{EC}}_{\mathrm{true}}\times {\mathrm{LOD}}_{\mathrm{EC}}\times {\mathit{\alpha }}_{\mathrm{EC}}^{\mathrm{2}}}}.\end{array}$$

Besides the γ_{Unc−nonlinear} discussed above, a second type measurement uncertainty parameterized by a constant proportional factor, γ_Unc−linear, is very common in atmospheric applications:

$$\begin{array}{}\text{(23)}& {\displaystyle }& {\displaystyle }-{\mathit{\gamma }}_{\mathrm{POCunc}}\times {\mathrm{POC}}_{\mathrm{true}}\le {\mathit{\epsilon }}_{\mathrm{POC}}\le +{\mathit{\gamma }}_{\mathrm{POCunc}}\times {\mathrm{POC}}_{\mathrm{true}},\text{(24)}& {\displaystyle }& {\displaystyle }-{\mathit{\gamma }}_{\mathrm{ECunc}}\times {\mathrm{EC}}_{\mathrm{true}}\le {\mathit{\epsilon }}_{\mathrm{EC}}\le +{\mathit{\gamma }}_{\mathrm{ECunc}}\times {\mathrm{EC}}_{\mathrm{true}}.\end{array}$$

where γ_POCunc and γ_ECunc are the relative measurement uncertainties, e.g., for relative measurement uncertainty of 10 %, γ_Unc= 0.1. As a result, the measurement error is linearly proportional to the concentration. An example comparison of γ_{Unc−nonlinear} and γ_Unc−linear is shown in Fig. 1b. For γ_Unc−linear, the weights become

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }\mathit{\omega }\left(X_i\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{X_i}^{\mathrm{2}}}}={\displaystyle \frac{\mathrm{3}}{{\left({\mathit{\gamma }}_{\mathrm{ECunc}}\times {\mathrm{EC}}_{\mathrm{true}}\right)}^{\mathrm{2}}}},\text{(26)}& {\displaystyle }& {\displaystyle }\mathit{\omega }\left(Y_i\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{Y_i}^{\mathrm{2}}}}={\displaystyle \frac{\mathrm{3}}{{\left({\mathit{\gamma }}_{\mathrm{POCunc}}\times {\mathrm{POC}}_{\mathrm{true}}\right)}^{\mathrm{2}}}},\end{array}$$

and λ for Deming regression can be determined:

$$\begin{array}{}\text{(27)}& \mathit{\lambda }={\displaystyle \frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}}={\displaystyle \frac{{\left({\mathit{\gamma }}_{\mathrm{POCunc}}\times {\mathrm{POC}}_{\mathrm{true}}\right)}^{\mathrm{2}}}{{\left({\mathit{\gamma }}_{\mathrm{ECunc}}\times {\mathrm{EC}}_{\mathrm{true}}\right)}^{\mathrm{2}}}}.\end{array}$$

3.1.2 XY data generation by the Mersenne twister generator following a specific distribution

The Mersenne twister (MT) is a pseudorandom number generator (PRNG) developed by Matsumoto and Nishimura (1998). MT has been widely adopted by mainstream numerical analysis software (e.g., MATLAB, SPSS, SAS and Igor Pro) as well as popular programing languages (e.g., R, Python, IDL, C$++$ and PHP). Data generation using MT provides a few advantages: (1) frequency distribution can be easily assigned during the data generation process, allowing straightforward simulation of the frequency distribution characteristics (e.g., Gaussian or lognormal) observed in ambient measurements; (2) the inputs for data generation are simply the mean and standard deviation of the data series and can be changed easily by the user; (3) the correlation (R²) between X and Y can be manipulated easily during the data generation to satisfy various purposes; and (4) unlike the sine function described by Chu (2005), which has a sample size limitation of 120, the sample size in MT data generation is highly flexible.

In this section, we will use POC as Y and EC as X as an example to explain the data generation. Procedure of applying MT to simulate ambient POC and EC data can be found in our previous study (Wu and Yu, 2016). Details of the data generation steps are shown in Fig. 2 and described below. The first step is generation of EC_true by MT. In our previous study, it was found that ambient POC and EC data follow a lognormal distribution in various locations of the Pearl River Delta (PRD) region. Therefore, lognormal distributions are adopted during EC_true generation. A range of average concentration and relative standard deviation (RSD) from ambient samples is considered in formulating the lognormal distribution. The second step is to generate POC_comb. As shown in Fig. 2, POC_comb is generated by multiplying EC_true with (OC ∕ EC)_pri. Instead of having a Gaussian distribution, (OC ∕ EC)_pri in this study is a single value, which favors direct comparison between the true value of (OC ∕ EC)_pri and (OC ∕ EC)_pri estimated from the regression slope. The third step is generation of POC_true by adding POC_non-comb onto POC_comb. Instead of having a distribution, POC_non-comb in this study is a single value, which favors direct comparison between the true value of POC_non-comb and POC_non-comb estimated from the regression intercept. The fourth step is to compute ε_POC and ε_EC. As discussed in Sect. 3.1.1, two types of measurement errors are considered for ε_POC and ε_EC calculation: γ_{Unc−nonlinear} and γ_Unc−linear. In the last step, POC_measured and EC_measured are calculated following Eq. (12), i.e., applying measurement errors on POC_true and EC_true. Then POC_measured and EC_measured can be used as Y and X, respectively, to test the performance of various regression techniques. An Igor Pro-based program with a GUI was developed to facilitate the MT data generation for OC and EC. A brief introduction is given in the Supplement.

Figure 2 Flowchart of data generation steps using MT.

3.1.3 XY data generation by the sine function of Chu (2005)

Besides MT, inclusion of the sine function data generation scheme in this study mainly serves two purposes. First, the sine function scheme was adopted in two previous studies (Chu, 2005; Saylor et al., 2006), the inclusion of this scheme can help to verify whether the codes in Igor for various regression approaches yield the same results from the two previous studies. Second, the crosscheck between results from sine function and MT provides circumstantial evidence that the MT scheme works as expected.

In this section, XY data generation by sine functions is demonstrated using POC as Y and EC as X. There are four steps in POC and EC data generation as shown by the flowchart in Fig. S2. Details are explained as follows. (1) The first step is to generate POC and EC (Chu, 2005):

$$\begin{array}{}\text{(28)}& {\displaystyle }& {\displaystyle }{\mathrm{POC}}_{\mathrm{comb}}=\mathrm{14}+\mathrm{12}\left(\sin \left({\displaystyle \frac{x}{\mathit{\tau }}}\right)+\sin \left(x-\mathit{\varphi }\right)\right),\text{(29)}& {\displaystyle }& {\displaystyle }{\mathrm{EC}}_{\mathrm{true}}=\mathrm{3.5}+\mathrm{3}\left(\sin \left({\displaystyle \frac{x}{\mathit{\tau }}}\right)+\sin \left(x-\mathit{\varphi }\right)\right),\end{array}$$

where x is the elapsed hour ($x=\mathrm{1},\mathrm{2},\mathrm{3}\mathrm{\dots }n$; n ≤ 120), τ is used to adjust the width of each peak, and ϕ is used to adjust the phase of the sine wave. The constants 14 and 3.5 are used to lift the sine wave to the positive range of the y axis. An example of data generation by the sine functions of Chu (2005) is shown in Fig. 3. Dividing Eq. (28) by Eq. (29) yields a value of 4. In this way the exact relation between POC and EC is defined clearly as (OC ∕ EC)_pri= 4. (2) With POC_comb and EC_true generated, the second step is to add POC_non-comb to POC_comb to compute POC_true. As for POC_non-comb, a single value is assigned and added to all POC following Eq. (10). Then the goodness of the regression intercept can be evaluated by comparing the regressed intercept with preset POC_non-comb. (3) The third step is to compute ε_POC and ε_EC, considering both γ_{Unc−nonlinear} and γ_Unc−linear. (4) The last step is to apply measurement errors on POC_true and EC_true following Eq. (12). Then POC_measured and EC_measured can be used as Y and X, respectively, to evaluate the performance of various regression techniques.

Figure 3 POC_comb and EC_true data generated by the sine functions of Chu (2005). (a) Time series of the 120 data points for POC_comb and EC_true. (b) Scatter plot of POC_comb vs. EC_true.

3.2 Ambient measurement of σ_abs and EC

Sampling was conducted from Feb 2012 to Jan 2013 at the suburban Nancun (NC) site (23^∘0^′11.82^′′ N, 113^∘21^′18.04^′′ E), which is situated on top of the highest peak (141 m a.s.l.) in the Panyu district of Guangzhou. This site is located at the geographic center of Pearl River Delta region (PRD), making it a good location for representing the average atmospheric mixing characteristics of city clusters in the PRD region. Light absorption measurements were performed by a 7λ Aethalometer (AE-31, Magee Scientific Company, Berkeley, CA, USA). EC mass concentrations were measured by a real time ECOC analyzer (model RT-4, Sunset Laboratory Inc., Tigard, Oregon, USA). Both instruments utilized inlets with a 2.5 µm particle diameter cutoff. The algorithm of Weingartner et al. (2003) was adopted to correct the sampling artifacts (aerosol loading, filter matrix and scattering effect; Collaud Coen et al., 2010) in Aethalometer measurement. A customized computer program with GUI, Aethalometer data processor (Wu et al., 2018), was developed to perform the data correction and detailed descriptions can be found in https://sites.google.com/site/wuchengust. More details of the measurements can be found in Wu et al. (2018).

4 Comparison study using synthetic data

In the following comparisons, six regression approaches are compared using two data generation schemes (Chu sine function and MT) separately, as illustrated in Fig. 4. Each data generation scheme considers both γ_{Unc−nonlinear} and γ_Unc−linear in measurement error parameterization. In total, 18 cases are tested with different combination of data generation schemes, measurement error parameterization schemes, true slope and intercept settings. In each case, six regression approaches are tested, i.e., OLS, DR (λ=1), DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$, ODR, WODR and YR. In commercial software (e.g., OriginPro^®, SigmaPlot^®, GraphPad Prism^®), λ in DR is set to 1 by default if not specified. As indicated by Saylor et al. (2006), the bias observed in the study of Chu (2005) is likely due to λ=1 in DR. The purpose of including DR (λ=1) in this study is to examine the potential bias using the default input in many software products. The six regression approaches are considered to examine the sensitivity of regression results to various parameters used in data generation. For each case, 5000 runs are performed to obtain statistically significant results, as recommended by Saylor et al. (2006). The mean slope and intercept from 5000 runs is compared with the true value assigned during data generation. If the difference is < 5 %, the result is considered unbiased.

Figure 4 Overview of the comparison study design.

4.1 Comparison results using the dataset of Chu (2005)

In this section, the scheme of Chu (2005) is adopted for data generation to obtain a benchmark of six regression approaches. With different setup of slope, intercept and γ_Unc, six cases (Cases 1–6) are studied and the results are discussed below.

Table 1 Summary of six regression approaches comparison with 5000 runs for 18 cases.

4.1.1 Results with γ_{Unc−nonlinear}

A comparison of the regression techniques results with γ_{Unc−nonlinear} (following Eqs. 18 and 19) is summarized in Table 1. LOD_POC , LOD_EC, α_POC and α_EC are all set to 1 to reproduce the data studied by Chu (2005) and Saylor et al. (2006). Two sets of true slope and intercept are considered (Case 1: slope = 4, intercept = 0; Case 2: slope = 4, intercept = 3) to examine if any results are sensitive to the nonzero intercept. The R² (POC, EC) from 5000 runs for both Case 1 and 2 are 0.67 ± 0.03.

As shown in Fig. 5, for the zero-intercept case (Case 1), OLS significantly underestimates the slope (2.95 ± 0.14), while it overestimates the intercept (5.84 ± 0.78). This result indicates that OLS is not suitable for errors-in-variables linear regression, consistent with similar analysis results from Chu (2005) and Saylor et al. (2006). With DR, if the λ is properly calculated by weights $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$, unbiased slope (4.01 ± 0.25) and intercept (−0.04 ± 1.28) are obtained; however, results from DR with λ= 1 show obvious bias in the slope (4.27 ± 0.27) and intercept (−1.45 ± 1.36). ODR also produces biased slope (4.27 ± 0.27) and intercept (−1.45 ± 1.36), which are identical to results of DR when λ= 1. With WODR, unbiased slope (3.98 ± 0.22) is observed, but the intercept is overestimated (1.12 ± 1.02). Results of YR are identical to WODR. For Case 2 (slope = 4, intercept = 3), slopes from all six regression approaches are consistent with Case 1 (Table 1). The Case 2 intercepts are equal to the Case 1 intercepts plus 3, implying that all the regression methods are not sensitive to a nonzero intercept.

Figure 5 Regression results on synthetic data, Case 1 (Slope = 4, Intercept = 0, LOD_POC= 1, LOD_EC= 1, a_POC= 1, a_EC= 1, R² (POC, EC) = 0.67 ± 0.03). The scatter plots demonstrate regression examples from a single run. The box plots show the distribution of regressed slopes and intercepts from 5000 runs of six regression approaches. The dashed lines in orange and pink represent true slope and intercept, respectively.

For Case 3, LOD_POC= 0.5, LOD_EC= 0.5, α_POC= 0.5, α_EC= 0.5 are adopted (Table 1), leading to an offset to the left of γ_{Unc−nonlinear} (blue curve) compared to Case 1 and 2 (black curve) in Fig. 1. As a result, for the same concentration of EC and OC in Case 3, the γ_{Unc−nonlinear} is smaller than in Cases 1 and 2 as indicated by a higher R² (0.95 ± 0.01 for Case 3, Table 1). With a smaller measurement uncertainty, the degree of bias in Case 3 is smaller than in Case 1. For example, OLS slope is less biased in Case 3 (3.83 ± 0.08) compared to Case 1 (2.94 ± 0.14). Similarly, the slope (4.03 ± 0.09) and intercept (−0.18 ± 0.44) of DR (λ= 1) exhibit a much smaller bias with a smaller measurement uncertainty, implying that the degree of bias by improperly weighting in DR, WODR and YR is associated with the degree of measurement uncertainty. A higher measurement uncertainty results in larger bias in slope and intercept.

Figure 6 Slope and intercept biases by different regression schemes in two test scenarios (A and B) in which the assumed error for one of the regression variables deviates from the actual measurement error. In Test A data generation, γ_{Unc_X} is fixed at 30 % and γ_{Unc_Y} is varied between 1 and 50 %. In Test B, γ_{Unc_X} varies between 1 and 50 % and γ_{Unc_Y} is fixed at 30 %. The assumed measurement error for regression is 10 % for both X and Y. (a) Slope biases as a function of γ_{Unc_Y} in Test A. (b) Intercept biases as a function of γ_{Unc_Y} in Test A. (c) Slope biases as a function of γ_{Unc_X} in Test B. (d) Intercept biases as a function of γ_{Unc_X} in Test B.

An uneven LOD_POC and LOD_EC is tested in Case 4 with LOD_POC= 1, LOD_EC= 0.5, α_POC= 0.5, α_EC= 0.5, which yield an R²(POC, EC) of 0.78 ± 0.02. The results are similar to Case 1. For DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$ unbiased slope and intercept are obtained. For WODR and YR, unbiased slopes are reported with a small bias in the intercepts. Large bias values are observed in both the slopes and intercepts in Case 4 using OLS, DR (λ=1) and ODR.

4.1.2 Results with γ_Unc−linear

Cases 5 and 6 represent the results from using γ_Unc−linear and are shown in Table 1. γ_Unc is set to 30 % to achieve an R² (POC, EC) of 0.7, a value close to the R² in studies of Chu (2005) and Saylor et al. (2006). In Case 5 (slope = 4, intercept = 0), unbiased slopes and intercepts are determined by DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$, WODR and YR. OLS underestimates the slope (3.32 ± 0.20) and overestimates intercept (3.77 ± 0.90), while DR (λ=1) and ODR overestimate the slopes (4.75 ± 0.30) and underestimate the intercepts (−4.14 ± 1.36). In Case 6 (slope = 4, intercept = 3), results similar to Case 5 are obtained. It is worth noting that although the mean intercept (3.05 ± 1.22) of DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$ is closest to the true value (intercept = 3), the deviations are much larger than for WODR (2.72 ± 0.74).

4.2 Comparison results using data generated by MT

In this section, MT is adopted for data generation to obtain a benchmark of six regression approaches. Both γ_{Unc−nonlinear} and γ_Unc−linear are considered. With different configuration of slope, intercept and γ_Unc, 12 cases (Cases 7–18) are studied and the results are discussed below.

4.2.1 γ_{Unc−nonlinear} results

Cases 7 and 8 use data generated by MT and γ_{Unc−nonlinear} with results shown in Table 1. In Case 7 (slope = 4, intercept = 0, LOD_POC= 1, LOD_EC= 1, α_POC= 1, α_EC= 1), unbiased slope (4.00 ± 0.03) and intercept (0.00 ± 0.17) is estimated by DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$. WODR and YR yield unbiased slopes (3.96 ± 0.03) but overestimate the intercepts (1.21 ± 0.13). DR (λ=1) and ODR report slightly biased slopes (4.17 ± 0.04) with biased intercepts (−0.94 ± 0.18). OLS underestimates the slope (3.22 ± 0.03) and overestimates the intercept (4.30 ± 0.14). In Case 8 (slope = 4, intercept = 3, LOD_POC= 1, LOD_EC= 1, α_POC= 1, α_EC= 1), DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$ provides unbiased slope (4.00 ± 0.03) and intercept (3.00 ± 0.18) estimations. WODR and YR report unbiased slopes (3.97 ± 0.03) and overestimate intercepts (4.11 ± 0.13). OLS, DR (λ=1) and ODR report biased slopes and intercepts.

To test the overestimation/underestimation dependency on the true slope, Case 9 (slope = 0.5, intercept = 0, LOD_POC= 1, LOD_EC= 1, α_POC= 1, α_EC= 1) and Case 10 (slope = 0.5, intercept = 3, LOD_POC= 1, LOD_EC= 1, α_POC= 1, α_EC= 1) are conducted and the results are shown in Table 1. Unlike the overestimation observed in Cases 1–8, DR (λ=1) and ODR underestimate the slopes (0.46 ± 0.01) in Case 9. In Case 10, DR (λ=1), DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$ and ODR report unbiased slopes and intercepts. Cases 11 and 12 test the bias when the true slope is 1, as shown in Table 1. In Case 11 (intercept = 0), all regression approaches except OLS can provide unbiased results. In Case 12, all regression approaches report unbiased slopes except OLS, but DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$ is the only regression approach that reports unbiased intercept.

These results imply that if the true slope is less than 1, the improper weighting (λ=1) in Deming regression and ODR without weighting tends to underestimate slope. If the true slope is 1, these two estimators can provide unbiased results. If the true slope is larger than 1, the improper weighting (λ=1) in Deming regression and ODR without weighting tends to overestimate slope.

4.2.2 γ_Unc−linear results

Cases 13 and 14 (Table 1) represent the results from using γ_Unc−linear (30 %) and data generated from MT. For Case 13 (slope = 4, intercept = 0), DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$, WODR and YR provide the best estimation of slopes and intercepts. DR (λ=1) and ODR overestimate slopes (4.53 ± 0.05) and underestimate intercepts (−2.94 ± 0.24). For Case 14 (slope = 4, intercept = 3), DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$, WODR and YR provide an unbiased estimation of slopes. But DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$ is the only regression approach reporting unbiased intercept (3.08 ± 0.23). Cases 15 and 16 are tested to investigate whether the results are different if the true slope is smaller than 1. As shown in Table 1, the results are similar to Cases 13 and 14, i.e., that DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$ can provide unbiased slope and intercept while WODR and YR can provide unbiased slopes but biased intercepts. Cases 17 and 18 are tested to see if the results are the same for a special case when the true slope is 1. As shown in Table 1, the results are similar to Cases 13 and 14, implying that these results are not sensitive to the special case when the true slope is 1.

Figure 7 Regression results using ambient σ_abs520 and EC data from a suburban site in Guangzhou, China.

4.3 The importance of appropriate λ input for Deming regression

As discussed above, inappropriate λ assignment in the Deming regression (e.g., λ= 1 by default for much commercial software) leads to biased slope and intercept. Besides λ= 1, inappropriate λ input due to improper handling of measurement uncertainty can also result in bias for Deming regression. An example is shown in Fig. S3. Data are generated by MT with following parameters: slope = 4, intercept = 0, and γ_Unc−linear (30 %). Figure S2a and b demonstrate that when an appropriate λ is provided (following γ_Unc−linear, $\mathit{\lambda }=\frac{{\mathrm{POC}}^{\mathrm{2}}}{{\mathrm{EC}}^{\mathrm{2}}})$, unbiased slopes and intercepts are obtained. If an improper λ is used due to a mismatched measurement uncertainty assumption $\left({\mathit{\gamma }}_{\mathrm{Unc}-\mathrm{nonlinear}},\mathit{\lambda }=\frac{\mathrm{POC}}{\mathrm{EC}}\right)$, the slopes are overestimated (Fig. S3c, 4.37 ± 0.05) and intercepts are underestimated (Fig. S3, −2.01 ± 0.24). This result emphasizes the importance of determining the correct form of measurement uncertainty in ambient samples, since λ is a crucial parameter in Deming regression.

In the λ calculation, different representations for POC and EC, including mean, median and mode, are tested as shown in Fig. S4. The results show that when X and Y have a similar distribution (e.g., both are lognormal), any of mean, median or mode can be used for the λ calculation.

4.4 Caveats of regressions with unknown X and Y uncertainties

In atmospheric applications, there are scenarios in which a priori error in one of the variables is unknown, or the measurement error described cannot be trusted. For example, in the case of comparing model prediction and measurement data, the uncertainty of model prediction data is unknown. A second example is the case in which measurement uncertainty cannot be determined due to the lack of duplicated or collocated measurements and as a result, an arbitrarily assumed uncertainty is used. Such a case was illustrated in the study by Flanagan et al. (2006). They found that in the Speciation Trends Network (STN), the whole-system uncertainty retrieved by data from collocated samplers was different from the arbitrarily assumed 5 % uncertainty. Additionally, the discrepancy between the actual uncertainty obtained through collocated samplers and the arbitrarily assumed uncertainty varied by chemical species. To investigate the performance of different regression approaches in these cases, two tests (A and B) are conducted.

In Test A, the actual measurement error for X is fixed at 30 %, while γ_Unc for Y varies from 1 to 50 %. The assumed measurement error for regression is 10 % for both X and Y. Results of Test A are shown in Fig. 6a and b. For OLS, the slopes are underestimated (−14 to −12 %) and intercepts are overestimated (90–103 %) and the biases are independent of variations in γ_{Unc_Y}. ODR and DR (λ=1) yield similar results with overestimated slopes (0–44 %) and underestimated intercepts (−330–0 %). The degree of bias in slopes and intercepts depends on the γ_{Unc_Y}. WODR, DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$ and YR perform much better than other regression approaches in Test A, with a smaller bias in both slopes (−8–12 %) and intercepts −98–55 %).

Figure 8 The user interface of the Scatter Plot Igor program. The program and its operation manual are available from https://doi.org/10.5281/zenodo.832417.

In Test B, γ_{Unc_Y} is fixed at 30 % and γ_{Unc_X} varies between 1 and 50 %. The results of Test B are shown in Fig. 6c and d. The assumed measurement error for regression is 10 % for both X and Y. OLS underestimates the slopes (−29 to −0.2 %) and overestimates the intercepts (2–209 %). In contrast to Test A in which slope and intercept biases are independent of variations in γ_{Unc_Y}, the slope and intercept biases in Test B exhibit dependency on γ_{Unc_X}. The reason for this is that OLS only considers errors in Y and X is assumed to be error-free. ODR and DR (λ=1) yield similar results with overestimated slopes (11–18 %) and underestimated intercepts (−144 to −87 %). The degree of bias in slopes and intercepts is relatively independent on the γ_{Unc_X}. WODR, DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$ and YR performed much better than the other regression approaches in Test B, with a smaller bias in both slopes (−14–8 %) and intercepts (−59–106 %).

The results from these two tests suggest that, if one of the measurement errors described cannot be trusted or a priori error in one of the variables is unknown, WODR, DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$ and YR should be used instead of ODR and DR (λ=1) and OLS. This conclusion is consistent with results presented in Sect. 4.1 and 4.2. This analysis, albeit crude, also suggests that, in general, the magnitude of bias in slope estimation by these regression approaches is smaller than those for intercept. In other words, slope is a more reliable quantity compared to intercept when extracting quantitative information from linear regressions.

5 Regression applications to ambient data

This section demonstrates the application of the six regression approaches on a light absorption coefficient and EC dataset collected in a suburban site in Guangzhou. As mentioned in Sect. 4.4, measurement uncertainties are crucial inputs for DR, YR and WODR. The measurement precision of Aethalometer is 5 % (Hansen, 2005), while EC by the RT-ECOC analyzer is 24 % (Bauer et al., 2009). These measurement uncertainties are used in DR, YR and WODR calculation. The dataset contains 6926 data points with an R² of 0.92.

As shown in Fig. 7, the y axis is light absorption at 520 nm (σ_abs520) and the x axis is EC mass concentration. The regressed slopes represent the mass absorption efficiency (MAE) of EC at 520 nm, ranging from 13.66 to 15.94 m² g⁻¹ by the six regression approaches. OLS yields the lowest slope (13.66 as shown in Fig. 7a) among all six regression approaches, consistent with the results using synthetic data. This implies that OLS tends to underestimate regression slope when mean Y to X ratio is larger than 1. DR (λ=1) and ODR report the same slope (14.88) and intercept (5.54); this equivalency is also observed for the synthetic data. Similarly, WODR and YR yield identical slope (14.88) and intercept (5.54), in line with the synthetic data results. The regressed slope by DR (λ=1) is higher than DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$, and this relationship agrees well with the synthetic data results.

Regression comparison is also performed on hourly OC and EC data. Regression on OC ∕ EC percentile subset is a widely used empirical approach for primary OC ∕ EC ratio determination. Figure S5 shows the regression slopes as a function of OC ∕ EC percentile. OC ∕ EC percentile ranges from 0.5 to 100 %, with an interval of 0.5 %. As the percentile increases, SOC contribution in OC increases as well, resulting in decreased R² between OC and EC. The deviations between six regression approaches exhibit a dependency on R². When percentile is relatively small (e.g., < 10 %), the differences between the six regression approaches are also small due to the high R² (0.98). The deviations between the six regression approaches become more pronounced as R² decreases (e.g., < 0.9). The deviations are expected to be even larger when R² is less than 0.8. These results emphasize the importance of applying error-in-variables regression, since ambient XY data more likely has an R² less than 0.9 in most cases.

As discussed in this section, the ambient data confirm the results obtained in comparing methods with the synthetic data. The advantage of using the synthetic data for regression approaches evaluation is that the ideal slope and intercept are known values during the data generation, so the bias of each regression approach can be quantified.

6 Recommendations and conclusions

This study aims to provide a benchmark of commonly used linear regression algorithms using a new data generation scheme (MT). Six regression approaches are tested, i.e., OLS, DR (λ=1), DR $\left(\mathit{\lambda }=\frac{\mathit{\omega }\left(X_i\right)}{\mathit{\omega }\left(Y_i\right)}\right)$, ODR, WODR and YR. The results show that OLS fails to estimate the correct slope and intercept when both X and Y have measurement errors. This result is consistent with previous studies. For ambient data with R² less than 0.9, error-in-variables regression is needed to minimize the biases in slope and intercept. If measurement uncertainties in X and Y are determined during the measurement, measurement uncertainties should be used for regression. With appropriate weighting, DR, WODR and YR can provide the best results among all tested regression techniques. Sensitivity tests also reveal the importance of the weighting parameter λ in DR. An improper λ could lead to biased slope and intercept. Since the λ estimation depends on the form of the measurement errors, it is important to determine the measurement errors during the experimentation stage rather than making assumptions. If measurement errors are not available from the measurement and assumptions are made on measurement errors, DR, WODR and YR are still the best option that can provide the least bias in slope and intercept among all tested regression techniques. For these reasons, DR, WODR and YR are recommended for atmospheric studies when both X and Y data have measurement errors.

Application of error-in-variables regression is often overlooked in atmospheric studies, partly due to the lack of a specified tool for the regression implementation. To facilitate the implementation of error-in-variables regression (including DR, WODR and YR), a computer program (Scatter Plot) with a GUI in Igor Pro (WaveMetrics, Inc. Lake Oswego, OR, USA) was developed (Fig. 8). It is packed with many useful features for data analysis and plotting, including batch plotting, data masking via GUI, color coding in the z axis, data filtering and grouping by numerical values and strings. The Scatter Plot program and user manual are available from https://sites.google.com/site/wuchengust and https://doi.org/10.5281/zenodo.832417.

Appendix A: Equations of regression techniques

Ordinary least squares (OLS) calculation steps.

First calculate average of observed X_i and Y_i.

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{X}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{X}_i}{N}}\text{(A2)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{Y}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{Y}_i}{N}}\end{array}$$

Then calculate S_xx and S_yy.

$$\begin{array}{}\text{(A3)}& {\displaystyle }& {\displaystyle }{S}_{xx}={\sum }_{i=\mathrm{1}}^N{\left(X_i-\stackrel{\mathrm{‾}}{X}\right)}^{\mathrm{2}}\text{(A4)}& {\displaystyle }& {\displaystyle }{S}_{yy}={\sum }_{i=\mathrm{1}}^N{\left(Y_i-\stackrel{\mathrm{‾}}{Y}\right)}^{\mathrm{2}}\end{array}$$

OLS slope and intercept can be obtained from

$$\begin{array}{}\text{(A5)}& {\displaystyle }& {\displaystyle }k={\displaystyle \frac{S_{yy}}{S_{xx}}},\text{(A6)}& {\displaystyle }& {\displaystyle }b=\stackrel{\mathrm{‾}}{Y}-k\stackrel{\mathrm{‾}}{X}.\end{array}$$

Deming regression (DR) calculation steps (York, 1966).

Besides S_xx and S_yy as shown above, S_xy can be calculated from

$$\begin{array}{}\text{(A7)}& {\displaystyle }{S}_{xy}={\sum }_{i=\mathrm{1}}^N\left(X_i-\stackrel{\mathrm{‾}}{X}\right)\left(Y_i-\stackrel{\mathrm{‾}}{Y}\right),\end{array}$$

DR slope and intercept can be obtained from

$$\begin{array}{}\text{(A8)}& {\displaystyle }& {\displaystyle }k={\displaystyle \frac{S_{yy}-\mathit{\lambda }{S}_{xx}+\sqrt{{\left(S_{yy}-\mathit{\lambda }{S}_{xx}\right)}^{\mathrm{2}}+\mathrm{4}\mathit{\lambda }{S}_{xy}^{\mathrm{2}}}}{\mathrm{2}{S}_{xy}}},\text{(A9)}& {\displaystyle }& {\displaystyle }b=\stackrel{\mathrm{‾}}{Y}-k\stackrel{\mathrm{‾}}{X}.\end{array}$$

York regression (YR) iteration steps (York et al., 2004).

Slope by OLS can be used as the initial k in W_i calculation.

$$\begin{array}{}\text{(A10)}& {\displaystyle }& {\displaystyle }{W}_i={\displaystyle \frac{\mathit{\omega }\left(X_i\right)\mathit{\omega }\left(Y_i\right)}{\mathit{\omega }\left(X_i\right)+k^{\mathrm{2}}\mathit{\omega }\left(Y_i\right)-\mathrm{2}k{r}_i\sqrt{\mathit{\omega }\left(X_i\right)\mathit{\omega }\left(Y_i\right)}}}\text{(A11)}& {\displaystyle }& {\displaystyle }{U}_i=X_i-\stackrel{\mathrm{‾}}{X}=X_i-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{W}_i{X}_i}{{\sum }_{i=\mathrm{1}}^N{W}_i}}\text{(A12)}& {\displaystyle }& {\displaystyle }{V}_i=Y_i-\stackrel{\mathrm{‾}}{Y}=Y_i-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{W}_i{Y}_i}{{\sum }_{i=\mathrm{1}}^N{W}_i}}\end{array}$$

Then calculate β_i.

$$\begin{array}{}\text{(A13)}& {\displaystyle }{\mathit{\beta }}_i=W_i\left[{\displaystyle \frac{U_i}{\mathit{\omega }\left(Y_i\right)}}+{\displaystyle \frac{k{V}_i}{\mathit{\omega }\left(X_i\right)}}-\left[k{U}_i+V_i\right]{\displaystyle \frac{r_i}{\sqrt{\mathit{\omega }\left(X_i\right)\mathit{\omega }\left(Y_i\right)}}}\right]\end{array}$$

Slope and intercept can be obtained from

$$\begin{array}{}\text{(A14)}& {\displaystyle }& {\displaystyle }k={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{W}_i{\mathit{\beta }}_i{V}_i}{{\sum }_{i=\mathrm{1}}^N{W}_i{\mathit{\beta }}_i{U}_i}},\text{(A15)}& {\displaystyle }& {\displaystyle }b=\stackrel{\mathrm{‾}}{Y}-k\stackrel{\mathrm{‾}}{X}.\end{array}$$

Since W_i and β_i are functions of k, k must be solved iteratively by repeating Eqs. (A11) to (A15). If the difference between the k obtained from Eq. (A15) and the k used in Eq. (A11) satisfies the predefined tolerance ($\frac{k_{i+\mathrm{1}}-k_i}{k_i}<e^{-\mathrm{15}})$, the calculation is considered as converged. The calculation is straightforward and usually converged in 10 iterations. For example, the iteration count on the dataset of Chu (2005) is around 6.

Appendix B: Summary of abbreviations and symbols

Abbreviation/symbol	Definition
α	a dimensionless adjustable factor to control the position
	of γUnc curve on the concentration axis
b	intercept in linear regression
${\mathit{\beta }}_i,U_i,V_i,W_i$	intermediates in York regression calculations
γUnc	fractional measurement uncertainties
	relative to the true concentration (%)
DR	Deming regression
εEC,εPOC	absolute measurement uncertainties of EC and POC
EC	elemental carbon
ECtrue	numerically synthesized true EC concentration
	without measurement uncertainty
ECmeasured	EC with measurement error (ECtrue+εEC)
λ	ω(Xi) to ω(Yi)
	ratio in Deming regression
k	slope in linear regression
LOD	limit of detection
MT	Mersenne twister pseudorandom number generator
OC	organic carbon
OC ∕ EC	OC to EC ratio
(OC ∕ EC)pri	primary OC ∕ EC ratio
OCnon-comb	OC from non-combustion sources
ODR	orthogonal distance regression
OLS	ordinary least squares regression
POC	primary organic carbon
POCcomb	numerically synthesized true POC from combustion
	sources (well correlated with ECtrue),
	measurement uncertainty not considered
POCnon-comb	numerically synthesized true POC from non-combustion
	sources (independent of ECtrue)
	without considering measurement uncertainty
POCtrue	sum of POCcomb and POCnon-comb
	without considering measurement uncertainty
POCmeasured	POC with measurement error (POCtrue+εPOC)
${\mathit{\sigma }}_{X_i},{\mathit{\sigma }}_{Y_i}$	the standard deviation of the error in
	measurement of Xi and Yi
ri	correlation coefficient between errors in Xi and Yi in YR
S	sum of squared residuals
SOC	secondary organic carbon
τ	parameter in the sine function of Chu (2005)
	that adjusts the width of each peak
ϕ	parameter in the sine function of Chu (2005)
	that adjusts the phase of the curve
WODR	weighted orthogonal distance regression
$\stackrel{\mathrm{‾}}{X},\stackrel{\mathrm{‾}}{Y}$	average of Xi and Yi
YR	York regression
ω(Xi),ω(Yi)	inverse of ${\mathit{\sigma }}_{X_i}$ and ${\mathit{\sigma }}_{Y_i}$,
	used as weights in DR calculation.