Technical note: Effects of uncertainties and number of data points on line fitting – a case study on new particle formation

Fitting a line to two measured variables is considered one of the simplest statistical procedures researchers can carry out. However, this simplicity is deceptive as the line-fitting procedure is actually quite a complex problem. Atmospheric measurement data never come without some measurement error. Too often, these errors are neglected when researchers make inferences from their data.

To demonstrate the problem, we simulated datasets with different numbers of data points and different amounts of error, mimicking the dependence of the atmospheric new particle formation rate (J_1.7) on the sulfuric acid concentration (H₂SO₄). Both variables have substantial measurement error and, thus, are good test variables for our study. We show that ordinary least squares (OLS) regression results in strongly biased slope values compared with six error-in-variables (EIV) regression methods (Deming regression, principal component analysis, orthogonal regression, Bayesian EIV and two different bivariate regression methods) that are known to take errors in the variables into account.

1 Introduction

Atmospheric measurements always come with some measurement error. Too often, these errors are neglected when researchers make inferences based on their data. Describing the relationship between two variables typically involves making deductions in a more general context than that in which the variables were directly studied. If the relationship is not defined correctly, the inference is also not valid. In some cases, the bias in the analytical method is even given a physical meaning.

When analysing the dependencies of two or more measured variables, regression models are usually applied. Regression models can be linear or non-linear, depending on the relationship between the datasets that are analysed. Standard regression models assume that the independent variables of the model have been measured without error and that the model only accounts for errors in the dependent variables or responses. In cases where the measurements of the predictors contain error, estimating using standard methods (usually ordinary least squares, OLS) does not tend to provide the true parameter values, even when a very high number of data points is used. In linear models, the coefficients are underestimated (e.g. Carroll et al., 2006); however, in non-linear models, the bias is likely to be more complicated (e.g. Schennach, 2004). If predictor variables in regression analyses contain any measurement error, methods that account for errors should be applied – particularly when errors are large. Thus, test variables in this study were chosen such that they included significant uncertainties in both the independent and dependent variables.

The sulfuric acid concentration (H₂SO₄) is known to strongly affect the formation rates (J) of aerosol particles (Kirkby et al., 2016; Kuang et al., 2008; Kulmala et al., 2006; Kürten et al., 2016; Metzger et al., 2010; Riccobono et al., 2014; Riipinen et al., 2007; Sihto et al., 2006; Spracklen et al., 2006). The relationship between J (cm⁻³ s⁻¹) and H₂SO₄ (molec cm⁻³) is typically assumed to be in the following form: log${}_{\mathrm{10}}\left(J\right)=\mathit{\beta }\times {\log }_{\mathrm{10}}$(H₂SO₄)+α (Seinfeld and Pandis, 2016). In addition, parameterisations based on the results from these fits have been implemented in global models (e.g. in Dunne et al., 2016, Metzger et al., 2010 and Spracklen et al., 2006) to estimate the effects of new particle formation on global aerosol amounts and characteristics. Theoretically, in homogeneous nucleation, the slope of this relationship is related to the number of sulfuric acid molecules in the nucleating critical cluster, based on the first nucleation theorem (Vehkamäki, 2006).

Some published results have shown discrepancies in the expected J vs. H₂SO₄ dependence. Analysing data from Hyytiälä in 2003, Kuang et al. (2008) used an unconstrained least squares method, which was not specified in the paper, and obtained a β value of 1.99 for the slope, whereas Sihto et al. (2006) reported a β value of 1.16 using OLS from the same field campaign. The studies had some differences in pre-treatment of the data and used different time windows, but a significant proportion of this inconsistency is very likely due to the use of different fitting methods. The problem regarding the relationship of H₂SO₄ and J was previously acknowledged in Paasonen et al. (2010), who noted that the bivariate fitting method, as presented in York et al. (2004), should be applied; however, this method could not be used due to the lack of proper error estimates for each quantity. They were not aware of methods that did not require knowledge of the errors in advance, and instead made use of estimated variances. Here, we present the appropriate tools required for the above-mentioned approach.

Multiple attempts have been made to present methods that account for errors in the predictor variables for regression-type analyses, going back to Deming (1943). However, the traditional least squares fitting method is still the de facto line-fitting method due to its simplicity and common availability in frequently used software. In atmospheric sciences, Cantrell (2008) drew attention to the method introduced by York (1966) and York et al. (2004) and listed multiple other methodological papers utilising similar methodology. Pitkänen et al. (2016) raised awareness of the fact that errors are not accounted for in the predictor variables in the remote-sensing community, and this study partly follows their approach and introduces multiple methods to account for the errors in predictors. Cheng and Riu (2006) studied methods involving heteroscedastic errors, whereas Wu and Yu (2018) approached the problem with measurement errors via weighted regression and applied some techniques that are also used in our study.

Measurement errors in each variable must be taken into account using approaches known as errors-in-variables (EIV) regression. EIV methods simply mean that errors in both variables are accounted for. In this study, we compared OLS regression results to six different regression methods (Deming regression, principal component analysis regression, orthogonal regression, Bayesian EIV regression and two different bivariate regression methods) that are known to be able to take errors in variables into account and provide (at least asymptotically) unbiased estimates. In this study, we focus exclusively on linear EIV methods, but it is important to acknowledge that non-linear methods also exist, e.g. ORDPACK introduced in Boggs et al. (1987) and implemented in Python SciPy and R (Boggs et al., 1989; Spiess, 2015). ORDPACK is a somewhat improved version of classical orthogonal regression, in that arbitrary covariance structures are acceptable, and it is specifically set up so that a user can specify measurement error variance and covariance point by point; some of the methods in this study carry out the same process in linear analysis.

2 Materials and methods

2.1 Data illustrating the phenomenon

Measurement data contain different types of errors. Usually, the errors are divided to two main class: systematic error and random error.

Systematic errors, commonly referred as bias, in experimental observations usually come from the measurement instruments. They may occur because something is wrong with the instrument or its data handling system, or due to operator error. In line fitting, bias cannot be taken into account and needs to be minimised by way of careful and regular instrument calibrations and zeros or data preprocessing. The random error, in comparison, may have different components; the two components discussed here are the natural error and the measurement error. In addition, one should note the existence of equation error (discussed in Carroll and Ruppert, 1996), which refers to using an inappropriate form of a fitting equation. Measurement error is more generally understood; it is where measured values do not fully represent the true values of the variable being measured. This also contains sampling error (e.g. in the case of H₂SO₄ measurement, the sampled air in the measurement instrument is not a representative sample of the outside air due to losses of H₂SO₄ occurring in the sampling lines, among other factors). Natural error is the variability caused by natural or physical phenomenon (e.g. a specific amount of H₂SO₄ does not always cause the same number of new particles to be formed).

In the analysis of the measurement data, some amount of these errors are known or can be estimated, but some of the error will usually remain unknown; this should be kept in mind when interpreting fits. Even though the measurement error is taken into account, the regression fit may be biased due to unknown natural error. In this study, we assume that the errors of the different variables are uncorrelated, but in some cases this has to be accounted for, as noted, for example, in Trefall and Nordö (1959) and Mandel (1984). The correlation between the errors of two variables, measured with separate instruments, independent of each other, such as formation rate and H₂SO₄, may come from factors such as environmental variables that affect both of the variables at the same time. Factors affecting the formation of sulfuric acid have been studied in various papers, e.g. in Weber et al. (1997) and Mikkonen et al. (2011). New particle formation rates, in turn, have been studied in works such as Boy et al. (2008) and Hamed et al. (2011) and similarities between the affecting factors can be seen. In addition, factors like room temperature in the measurement space and atmospheric pressure may affect the performance of instrumentation, thereby causing additional error.

The data used in this study consist of simulated new particle formation rates at 1.7 nm (J_1.7) and sulfuric acid (H₂SO₄) concentrations mimicking observations of pure sulfuric acid in nucleation experiments from the CLOUD chamber at CERN (Kürten et al. 2016; https://home.cern/about/experiments/cloud, last access: 16 August 2019), including the corresponding expected values, their variances and covariance structures. The Proton Synchrotron at CERN provides an artificial source of “cosmic rays” that simulates the natural ionisation conditions between the ground level and the stratosphere. The core is a large (volume 26 m³) electro-polished stainless-steel chamber with temperature control (temperature stability better than 0.1 K) at any tropospheric temperature, precise delivery of selected gases (SO₂, O₃, NH₃ and various organic compounds) and ultra-pure humidified synthetic air, as well as very low gas-phase contaminant levels. The existing data on new particle formation include what are believed to be the most important formation routes that involve sulfuric acid, ammonia and water vapour (Kirkby et al., 2011); sulfuric acid and amine (Almeida et al., 2013); and ion-induced organic nucleation (Kirkby et al., 2016). The actual nucleation of new particles occurs at a slightly smaller size. After formation, particles grow by condensation to reach the detection limit (1.7 nm) of the instrument; thus, J_1.7 refers to the formation rate of particles as the instrument detects them, accounting for the known particle losses due to coagulation and deposition on the chamber walls. The relationships between precursor gas-phase concentrations and particle formation rates were chosen because they are both known to have considerable measurement errors and their relationship has been well-studied using regression-based analyses (Kirkby et al., 2016; Kürten et al., 2016; Riccobono et al., 2014; Tröstl et al., 2016). Additionally, many of the published papers on this topic do not describe how they accounted for the uncertainties in the analysis, which casts doubt on the fact that errors were treated properly. However, it should be kept in mind that the data could be any set of numbers assumed to have a linear relationship, but, in order to raise awareness in the aerosol research community, in this study we relate our analysis to the important problem of understanding new particle formation.

2.2 Regression methods

We carried out fits for the linear dependency of the logarithms of the two study variables, such that the equation for the fit was given by

$$\begin{array}{}\text{(1)}& y={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}x+\mathit{\epsilon },\end{array}$$

where y represents log₁₀(J_1.7), x is log₁₀(H₂SO₄), β values are the coefficients estimated from the data and ε is the error term. In order to demonstrate the importance of taking the measurement errors into account in the regression analysis, we tested seven different line-fitting methods: ordinary least squares (OLS), not taking the uncertainty in x variable into account; orthogonal regression (ODR; Boggs et al., 1987); Deming regression (DR; Deming, 1943); principal component analysis (PCA; Hotelling, 1957) regression; Bayesian EIV regression (Kaipio and Somersalo, 2005); and two different bivariate least squares methods by York et al. (2004) and Francq and Govaerts (2014, BLS), respectively, which are known to be able to account for errors in variables and provide (at least asymptotically) unbiased estimates. The differences between the methods stem from the criterion they minimise when calculating the coefficients and how they account for measurement errors. The minimising criteria for all methods are given in Appendix A1, but in the following we give the principles of the methods.

OLS minimises the squared distance of the observation and the fit line either in the x or y direction, but not both at the same time, whereas ODR minimises the sum of the squared weighted orthogonal distances between each point and the line. DR was originally an improved version of orthogonal regression, accounting for the ratio of the error variances, λ_xy, of the variables, (in classical non-weighted ODR λ_xy=1), and it is the maximum likelihood estimate (MLE) for the model (given in Eq. 1) when λ_xy is known. The PCA approach is the same as in ODR, but the estimation procedure is somewhat different as can be seen in Table S1 in the Supplement. The bivariate algorithm by York et al. (2004) provides a simple set of equations for iterating the MLE of the slope and intercept with weighted variables, which makes it similar to ODR in this case. However, using ODR allows for a regression to be performed on a user-defined model, whereas the York (2004) solution only works on linear models. This, for instance, enables the use of linear-scale uncertainties in ODR in this study, whereas the York (2004) approach could only use log-scale uncertainties. In Bayes EIV, statistical models for the uncertainties in the observed quantities are used and probability distributions for the line slope and intercept are computed according to Bayes' theorem. In this study, we computed the Bayesian maximum a posteriori (MAP) estimates for the slope and intercept that are the most probable values given the likelihood and prior models (see Appendix A1 for more details on models used in Bayes EIV). BLS takes errors and heteroscedasticity into account, i.e. unequal variances, in both variables; thus, it is a more advanced method than DR (under normality and equal variances, BLS is exactly equivalent to DR). PCA only accounts for the observed variance in data, whereas ODR, Bayes EIV and York bivariate regression require known estimates for measurement errors, although for Bayes EIV the error can be approximated with a distribution. DR and BLS can be applied with both errors given by the user and measurement variance-based errors. In this study, we applied measurement variance-based errors for these methods. The analyses for OLS and PCA were calculated with the “lm” and “prcomp” R functions (R Core Team, 2018), respectively, DR was calculated with the “deming” package (Therneau, 2018) and BLS was calculated with the “BivRegBLS” package (Francq and Berger, 2017) in R. The ODR-based estimates were obtained using the “scipy.odr” Python package (Jones et al., 2001), while the “PyStan” Python package (Stan Development Team, 2018) was used for calculating the Bayesian regression estimates. Finally, the York bivariate estimates were produced with a custom Python implementation of the algorithm presented by York et al. (2004).

3 Data

3.1 Simulated data

In measured data, the variables that are observed are not x and y, but (x+e_x) and (y+e_y), where e_x and e_y are the uncertainty in the measurements, and the true x and y cannot be exactly known. Thus, we used simulated data, where the true, i.e. noise-free, x and y were known to illustrate how the different line-fitting methods perform in different situations.

We simulated a dataset mimicking the new particle formation rates (J_1.7) and sulfuric acid concentrations (H₂SO₄) reported from CLOUD-chamber measurements at CERN. Both variables are known to have substantial measurement error and, thus, they are good test variables for our study. Additionally, the relationship of the logarithms of these variables is quite often described with linear OLS regression and, thus, the inference may be flawed.

We generated 1000 random noise-free H₂SO₄ concentration values assuming a log-normal distribution with a median of 2.0×10⁶ (molec cm⁻³) and a standard deviation of 2.4×10⁶ (molec cm⁻³). The corresponding noise-free J_1.7 was calculated using model log₁₀($J_{\mathrm{1.7}})=\mathit{\beta }\times {\log }_{\mathrm{10}}$(H₂SO₄)+α with the noise-free slope β=3.3 and $\mathit{\alpha }=-\mathrm{23}$, which are both realistic values presented by Kürten et al. (2016, Table 2 in their paper, for the no added ammonia cases).

Simulated observations of the noise-free H₂SO₄ concentrations were obtained by adding random errors $e_x=e_{\mathrm{rel},x}x+{\mathit{\sigma }}_{\mathrm{abs},x}$ that have a random absolute component $e_{\mathrm{abs},x}\sim$ normal($\mathrm{0},{\mathit{\sigma }}_{\mathrm{abs},x}$) and a random component relative to the observation x itself e_rel,xx, where $e_{\mathrm{rel},x}\sim$ normal($\mathrm{0},{\mathit{\sigma }}_{\mathrm{rel},x}$). Similar definitions apply for the noise-free J_1.7, e_y, σ_abs,y and σ_rel,y. The standard deviations of the measurement error components were chosen as ${\mathit{\sigma }}_{\mathrm{abs},x}=\mathrm{4}\times {\mathrm{10}}^{\mathrm{5}}$, ${\mathit{\sigma }}_{\mathrm{rel},x}=\mathrm{0.3}$, σ_abs,y = $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{3}}$ and ${\mathit{\sigma }}_{\mathrm{rel},y}=\mathrm{0.5}$, which are subjective estimates based on measurement data. The resulting total errors were occasionally about as large as the data values themselves; however, they are not unusually large error values with respect to corresponding real datasets, where overall uncertainties may reach 150 % for H₂SO₄ concentrations and 250 % for nucleation rates (e.g. Dunne et al., 2016).

These choices regarding generating simulated data reflect what real dataset can often be like: the bulk of the data approximates a log-normal distribution with one of the tails possibly being thinned or cut close to a limit of detection of an instrument or close to a limit of the data filtering criterion. In our simulated data, each negative observation and each negative noise-free value was replaced with a new random simulated value, which only slightly offsets the final distribution from a perfectly symmetric log-normal shape.

Simulating the observations tends to generate infrequent extreme outlier observations from the infinite tails of the normal distribution. We discarded outliers with an absolute error larger than 3 times the combined standard uncertainty of the observation in order to remove the effect of outliers from the regression analysis. This represents the quality control procedure in data analysis and also improves the stability of our results between different simulations.

3.2 Case study on measured data

In order to show that the results gained with simulated data are also applicable in real measurement data, we applied our methods to data measured in the CLOUD chamber and published by Dunne et al. (2016). Fig. 1 in Dunne et al. (2016) shows nucleation rates (J) at a 1.7 nm mobility diameter as a function of the sulfuric acid concentration. We used their measurements with no added ammonia at two different temperatures, 278 and 292 K, as shown in their Fig. 1D and E and given in their Supplementary material.

Figure 1 Regression lines fitted to the simulated data comparing all of the respective methods. The whiskers refer to the measurement error used for simulation.

4 Results

4.1 Fits for simulated data

Differences between the regression methods are illustrated in four different ways: firstly, by showing line fits on a scatterplot of simulated data; secondly, by illustrating how the slopes change when the uncertainty in the measured variables increase; thirdly, by showing the sensitivity of the fits on number of observations; and finally, by showing how the fits are affected by adding outliers in the data. Regression fits using each of the respective methods are shown in Fig. 1.

As we know that the noise-free slope β_true is 3.30, we can easily see how the methods perform. The worst performing method was OLS, with a β_ols value of 1.55, which is roughly half of the β_true. The best performing methods that displayed equal accuracy, i.e. within 2 % range, were ODR (β_ODR=3.27), Bayes EIV (β_BEIV=3.24) and BLS (β_BLS=3.22), whereas York (β_York=3.15) was within a range of 5 %; Deming (β_DR=2.95) and PCA (β_PCA=2.92), in comparison, slightly underestimated the slope.

The sensitivity of the methods was first tested by varying the uncertainty in the H₂SO₄ observations. We simulated six datasets with 1000 observations and with varying absolute and relative uncertainties (as listed in Table 1), and then performed fits with each method on all of these datasets. The performance of the methods is shown in Fig. 2, with the results corresponding to Fig. 1 marked in black. The results show that when the uncertainty is small, the bias in the OLS fit is smaller, but when more uncertainty is added to data, the bias increases significantly. A decrease in performance can also be seen with ODR, which overestimates the slope, and PCA, DR and Bayes EIV, which all underestimate the slope. The bivariate methods, BLS and York, seem to be quite robust with increasing uncertainty, as the slopes do not change significantly.

Table 1 The uncertainties used in the simulation for the sensitivity test for increasing uncertainty.

Figure 2 Sensitivity test for increasing uncertainty in simulated data. Black markers show the initial dataset described in Sect. 3. The dashed line indicates the noise-free slope.

The sensitivity of the methods to the decreasing n (number of observations) was tested by picking 100 random samples from the 1000-sample simulation dataset with n of 3, 5, 10, 20, 30, 50, 70, 100, 300 and 500 and carrying out fits for all samples using all methods. The average slopes and their standard errors are shown in Fig. 3. It is clear that when the n≤10, the variation in the estimated slopes can be considerably high. When n≥30 the average slopes stabilised close to their characteristic levels (within 5 %), except for Bayes EIV and York bivariate, which needed more than 100 observations. The most sensitive methods for a small n were Bayes EIV, ODR and PCA; thus, these methods should not be applied for data with a small n and a similar type of uncertainty to that presented here. However, the reader is reminded that number of points needed for a good fit depends on the uncertainties in the data.

Figure 3 Effect of sample size on the uncertainty of different fits. Lines show the median and shading illustrates the 1 standard deviation range of slope estimates for 40 repeated random samples. The dashed line indicates the noise-free slope.

The sensitivity of the predictor variable H₂SO₄ to outliers was tested using two different scenarios. In the first scenario, outliers were randomly allowed to be at either the high or low end of the distribution. In the second scenario, outliers were only allowed to be large numbers, which is often the case in H₂SO₄ and aerosol concentration measurements as numbers are removed from the data when they are smaller than the detection limit of the measurement instrument. Five cases with an n of 1000 were simulated with an increasing number of outliers (0, 5, 10, 20 and 100) and 10 repetitions of H₂SO₄ values with a different set of outliers. Outliers were defined such that $x_{\mathrm{obs}}-x_{\mathrm{true}}\mathit{>}\mathrm{3}\times$ the combined standard uncertainty. The methods most sensitive to outliers in both scenarios were OLS and Bayes EIV. A high number of outliers caused underestimations in PCA and DR, especially when using the outliers with high values (second scenario mentioned above), and a slight overestimation in BLS in the random outlier case (first scenario mentioned above). York bivariate and ODR were not affected in either case, and BLS only showed small variation between the 10 replicates in the estimated slope. We did not explore how large a number of outliers would needed to be to seriously disrupt the fits for the various methods. We felt that it is likely not realistic to have situations that have more than 10 % outliers.

We also applied an alternative method for simulating the data to different testing methods. The main difference compared with our method was that the distribution of noise-free H₂SO₄ followed a uniform distribution in log-space. With this assumption, it could be seen that OLS works almost as well as the EIV methods introduced here if the range of the data is wide (H₂SO₄ concentration in the range of 10⁶–10⁹). However, when scaled to concentrations usually measured in the atmosphere (10⁴–10⁷), the high uncertainties caused similar behaviour to the data seen in our previous simulations. Details of these results can be seen in Supplement S1.

Figure 4 Effect of outliers in the data. (a) The random outliers case and (b) the only high positive values case. Lines show the median and the shading shows ±1 standard deviation of the slope estimates in the 10 repeated studies. The dashed line indicates the noise-free slope.

Figure 5 Regression lines fitted to data from Dunne et al. (2016), in a similar fashion to Fig. 1. In (a) the observations are from a temperature of 292 K, and in (b) they are from 278 K.

4.2 Results of the case study

Figure 5 shows the fits of the data from Dunne et al. (2016). As expected, the fit using OLS is underestimated at both temperatures: β_ols(278 K) was 2.4 and β_ols(292 K) was 3.0. The regression equations for all methods are shown in Fig. 5. Dunne et al. (2016) did not use a linear fit in their study and instead applied a non-linear Levenberg–Marquardt algorithm (Moré, 1978) on function $J_{\mathrm{1.7}}=k\times [{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}{]}^{\mathit{\beta }},$ where k is a temperature-dependent rate coefficient with a non-linear function including three estimable parameters (see Sect. 8 of their Supplement for details). Thus, the results are not directly comparable as, for simplicity, we fit the data measured at different temperatures separately. However, their β value for the fit (β=3.95) is quite close to our results using EIV methods, especially as slopes from Bayes EIV at 292 K and BLS and PCA at both temperatures were within a range of 5 %. We also carried out some tests on data measured at lower temperatures (results not shown here). However, the slopes did not vary drastically from those at β_ols(278 K) and β_ols(292 K) when the other conditions were similar, even though the lower number of observations at lower temperatures increased uncertainty in the data. Nevertheless, the intercepts β₀(T) varied between temperatures.

5 Conclusions

Ordinary least squares regression can be used to answer some simple questions regarding data, such as “How is y related to x?”. However, if we are interested in the strength of the relationship and the predictor variable X contains some error, then error-in-variables methods should be applied. There is no single correct method to make the fit, as the methods behave slightly differently with different types of error. The choice of method should be based on the properties of data and the specific research question. There are usually two types of error in the data: natural and measurement error, where natural error refers to stochastic variation in the environment. Even if the natural error in the data is not known, taking the measurement error into account improves the fit significantly. Weighting the data based on some factor, typically the inverse of the uncertainty, reduces the effect of outliers and makes the regression more depend on the data that are more certain (see e.g. Wu and Yu, 2018), but it does not solve the problem completely.

As a test study, we simulated a dataset mimicking the dependence of the atmospheric new particle formation rate on the sulfuric acid concentration. We introduced three major sources of uncertainty when establishing inference from scatterplot data: the increasing measurement error, the number of data points and the number of outliers. In Fig. 1, we showed that for simulations where errors are taken from real measurements of J_1.7 and H₂SO₄ four of the methods gave slopes within 5 % of the known noise-free value: BLS, York bivariate, Bayes EIV and ODR. Estimates from BLS and York bivariate even remained stable when the uncertainty in the simulated H₂SO₄ concentration was increased drastically, as seen in Fig. 2. The main message to take away from Fig. 3, in comparison, is that if the data contain some error, all fit methods are highly uncertain when small numbers of observations are used. BLS was the most accurate with the smallest sample sizes (10 or less), ODR stabilised with 20 observations, and York bivariate and Bayes EIV needed 100 or more data points to become accurate. After that, these methods approached the noise-free value asymptotically, whereas the OLS slope converged towards an incorrect value. With an increasing number of outliers (Fig. 4), ODR and York bivariate were the most stable methods, even when 10 % of observations were classified as outliers in both test cases. BLS remained stable in the scenario with only high outlier values. Bayes EIV was the most sensitive to outliers after OLS.

From this, we can recommend that if the uncertainty in the predictor variable is known, York bivariate, or another method able to use known variances, should be applied. If the errors are not known, and they are estimated from data, BLS and ODR were found to be the most robust in cases with increasing uncertainty (relative error, rE > 30 % in Fig. 2) and with a high number of outliers. In our test data, BLS and ODR remained stable up to rE > 80 % (Fig. 2), whereas DR and PCA began to become more uncertain when rE > 30 % and Bayes EIV when rE > 50 %. If the number of observations is less than 10 and the uncertainties are high, we would recommend considering if a regression fit is appropriate at all. However, with the chosen uncertainties in our simulation tests, BLS was found to be the most robust with small numbers of data points. Bayes EIV displayed significant advantages if the number of observations was high enough and there were not too many outliers, as it did not require an explicit definition of the errors and could treat them as unknown parameters given their probability distributions.

We also carried out a case study on data measured in the CLOUD chamber and published by Dunne et al. (2016). In these analyses, we saw that our above-mentioned recommended methods also performed best for these data. Our tests indicated that the slope β₁ for the fit is not highly sensitive to changes in temperature in the chamber but the intercept β₀ in linear fit is. This dependency was also seen, and taken into account, in Dunne et al. (2016).

Appendix A: Minimising criteria for the regression methods applied in the paper

In this appendix, we introduce the minimising criteria (C_method) for all methods applied in the main text. We also give the equations for the regression coefficients (${\widehat{\mathit{\alpha }}}_{\mathrm{method}}$ and ${\widehat{\mathit{\beta }}}_{\mathrm{method}}$) for the methods.

A1 Ordinary least squares (OLS)

OLS minimises the sum of squares vertical distances (residuals) between each point and the fitted line. OLS regression minimises the following criterion:

$$\begin{array}{}\text{(A1)}& C_{\mathrm{OLS}}={\sum }_{i=\mathrm{1}}^N{\left(y_i-{\widehat{\mathit{\alpha }}}_{\mathrm{OLS}}-{\widehat{\mathit{\beta }}}_{\mathrm{OLS}}{x}_i\right)}^{\mathrm{2}},\end{array}$$

where ${\widehat{\mathit{\alpha }}}_{\mathrm{OLS}}$ and ${\widehat{\mathit{\beta }}}_{\mathrm{OLS}}$ refer to estimators calculated from the data. These estimations are given by

$$\begin{array}{}\text{(A2)}& {\widehat{\mathit{\beta }}}_{\mathrm{OLS}}={\displaystyle \frac{S_x}{S_y}},{\widehat{\mathit{\alpha }}}_{\mathrm{OLS}}=\stackrel{\mathrm{‾}}{x}-{\widehat{\mathit{\beta }}}_{\mathrm{OLS}}\stackrel{\mathrm{‾}}{y},\end{array}$$

where observed variances for $x{S}_x={\sum }_{i=\mathrm{1}}^N{\left(x_i-\stackrel{\mathrm{‾}}{x}\right)}^{\mathrm{2}}$ and for $y{S}_y={\sum }_{i=\mathrm{1}}^N{\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}$, and observed covariance for x and y $S_{xy}={\sum }_{i=\mathrm{1}}^N\left(x_i-\stackrel{\mathrm{‾}}{x}\right)\left(y_i-\stackrel{\mathrm{‾}}{y}\right)$.

A2 Orthogonal regression (ODR)

ODR (https://docs.scipy.org/doc/external/odrpack_guide.pdf, last access: 16 August 2019, https://docs.scipy.org/doc/scipy/reference/odr.html, last access: 27 July 2018) minimises the sum of the square of the orthogonal distances between each point and the line. The criteria is given by

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& C_{\mathrm{ODR}}={\sum }_{i=\mathrm{1}}^N\left({\left(x_i-{\displaystyle \frac{y_i+x_i/{\widehat{\mathit{\beta }}}_{\mathrm{ODR}}-{\widehat{\mathit{\alpha }}}_{\mathrm{ODR}}}{{\widehat{\mathit{\beta }}}_{\mathrm{ODR}}+\mathrm{1}/{\widehat{\mathit{\beta }}}_{\mathrm{ODR}}}}\right)}^{\mathrm{2}}\right.\\ & \left.+{\left(y_i-{\widehat{\mathit{\alpha }}}_{\mathrm{ODR}}-{\displaystyle \frac{{\widehat{\mathit{\beta }}}_{\mathrm{ODR}}{y}_i+x_i-{\widehat{\mathit{\alpha }}}_{\mathrm{ODR}}{\widehat{\mathit{\beta }}}_{\mathrm{ODR}}}{{\widehat{\mathit{\beta }}}_{\mathrm{ODR}}+\mathrm{1}/{\widehat{\mathit{\beta }}}_{\mathrm{ODR}}}}\right)}^{\mathrm{2}}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A4)}& {\widehat{\mathit{\beta }}}_{\mathrm{ODR}}={\displaystyle \frac{S_y-S_x+\sqrt{{\left(S_y-S_x\right)}^{\mathrm{2}}+\mathrm{4}{S}_{xy}^{\mathrm{2}}}}{\mathrm{2}{S}_{xy}}}\end{array}$$

and

$$\begin{array}{}\text{(A5)}& {\widehat{\mathit{\alpha }}}_{\mathrm{ODR}}=\stackrel{\mathrm{‾}}{y}-{\widehat{\mathit{\beta }}}_{\mathrm{ODR}}\stackrel{\mathrm{‾}}{x}.\end{array}$$

ODR accounts for the fact that errors exist in both axes but does not account for the exact values of the variances of variables. Thus, only the ratio between the two error variances (λ_xy) is needed to improve the methodology. With notation of Francq and Govaerts (2014), this ratio is given by the following:

$$\begin{array}{}\text{(A6)}& {\mathit{\lambda }}_{xy}={\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{{\mathit{\sigma }}_x^{\mathrm{2}}}},\end{array}$$

where the numerator of the ratio is the error variance in the data in the y axis and the denominator is the error variance in the data in the x axis.

A3 Deming regression (DR)

DR is the ML (maximum likelihood) solution of Eq. 1 when λ_xy is known. In practice, λ_xy is unknown and it is estimated from the variances of x and y calculated from the data.

The DR minimises the criterion C_DR, which is the sum of the square of (weighted) oblique distances between each point to the line:

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}& C_{\mathrm{DR}}={\sum }_{i=\mathrm{1}}^N\left({\mathit{\lambda }}_{xy}{\left(x_i-{\displaystyle \frac{y_i+{\mathit{\lambda }}_{xy}{x}_i/{\widehat{\mathit{\beta }}}_{\mathrm{DR}}-{\widehat{\mathit{\alpha }}}_{\mathrm{DR}}}{{\widehat{\mathit{\beta }}}_{\mathrm{DR}}+{\mathit{\lambda }}_{xy}/{\widehat{\mathit{\beta }}}_{\mathrm{DR}}}}\right)}^{\mathrm{2}}\right.\\ & \left.+{\left(y_i-{\widehat{\mathit{\alpha }}}_{\mathrm{DR}}-{\displaystyle \frac{{\widehat{\mathit{\beta }}}_{\mathrm{DR}}{y}_i+{\mathit{\lambda }}_{xy}{x}_i-{\widehat{\mathit{\alpha }}}_{\mathrm{DR}}{\widehat{\mathit{\beta }}}_{\mathrm{DR}}}{{\widehat{\mathit{\beta }}}_{\mathrm{DR}}+{\mathit{\lambda }}_{xy}/{\widehat{\mathit{\beta }}}_{\mathrm{DR}}}}\right)}^{\mathrm{2}}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A8)}& {\widehat{\mathit{\beta }}}_{\mathrm{DR}}={\displaystyle \frac{S_y-{\mathit{\lambda }}_{xy}{S}_x+\sqrt{{\left(S_y-{\mathit{\lambda }}_{xy}{S}_x\right)}^{\mathrm{2}}+\mathrm{4}{\mathit{\lambda }}_{xy}{S}_{xy}^{\mathrm{2}}}}{\mathrm{2}{S}_{xy}}}\end{array}$$

and

$$\begin{array}{}\text{(A9)}& {\widehat{\mathit{\alpha }}}_{\mathrm{DR}}=\stackrel{\mathrm{‾}}{y}-{\widehat{\mathit{\beta }}}_{\mathrm{DR}}\stackrel{\mathrm{‾}}{x}.\end{array}$$

A4 Bivariate least squares regression (BLS)

BLS is a generic name but here we refer to the formulation described in Francq and Govaerts (2014) and references therein. BLS takes errors and heteroscedasticity in both axes into account and is usually written in matrix notation. BLS minimises the criterion C_BLS, which is the sum of weighted residuals W_BLS given by the following:

$$\begin{array}{}\text{(A10)}& C_{\mathrm{BLS}}={\displaystyle \frac{\mathrm{1}}{W_{\mathrm{BLS}}}}{\sum }_{i=\mathrm{1}}^N{\left(y_i-{\widehat{\mathit{\alpha }}}_{\mathrm{BLS}}-{\widehat{\mathit{\beta }}}_{\mathrm{BLS}}{x}_i\right)}^{\mathrm{2}}\end{array}$$

with

$$\begin{array}{}\text{(A11)}& W_{\mathrm{BLS}}={\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{n_y}}+{\widehat{\mathit{\beta }}}_{\mathrm{BLS}}^{\mathrm{2}}{\displaystyle \frac{{\mathit{\sigma }}_x^{\mathrm{2}}}{n_x}}.\end{array}$$

Estimators for the parameters are computed by iterations using the following formulas:

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{W_{\mathrm{BLS}}}}\left(\begin{array}{cc}N& {\sum }_{i=\mathrm{1}}^N{x}_i\\ {\sum }_{i=\mathrm{1}}^N{x}_i& {\sum }_{i=\mathrm{1}}^N{x}_i^{\mathrm{2}}\end{array}\right)\left(\begin{array}{c}{\widehat{\mathit{\alpha }}}_{\mathrm{BLS}}\\ {\widehat{\mathit{\beta }}}_{\mathrm{BLS}}\end{array}\right)={\displaystyle \frac{\mathrm{1}}{W_{\mathrm{BLS}}}}\\ & \left({\sum }_{i=\mathrm{1}}^N\left(x_i{y}_i+{\widehat{\mathit{\beta }}}_{\mathrm{BLS}}{\displaystyle \frac{{\mathit{\sigma }}_x^{\mathrm{2}}}{n_x}}{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{\left(y_i-{\widehat{\mathit{\alpha }}}_{\mathrm{BLS}}-{\widehat{\mathit{\beta }}}_{\mathrm{BLS}}{x}_i\right)}^{\mathrm{2}}}{W_{\mathrm{BLS}}}}\right)\right),\end{array}\end{array}$$

where known uncertainties ${\mathit{\sigma }}_x^{\mathrm{2}}$ and ${\mathit{\sigma }}_y^{\mathrm{2}}$ are replaced with estimated variances S_x and S_y in this study.

A second bivariate regression method that was used in this study is an implementation of the regression method described by York et al. (2004, Section III). The minimisation criterion is described in York (1968):

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{\mathrm{C}}_{\mathrm{york}}& ={\sum }_{i=\mathrm{0}}^N{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-r_i^{\mathrm{2}}}}\left\{w\left(x_i\right){\left(x_{i,\mathrm{adj}}-x_i\right)}^{\mathrm{2}}\right.\\ & -\mathrm{2}r\sqrt{w\left(x_i\right)w\left(y_i\right)}\left(x_{i,\mathrm{adj}}-x_i\right)\left(y_{i,\mathrm{adj}}-y_i\right)\\ & \left.+w\left(y_i\right){\left(y_{i,\mathrm{adj}}-y_i\right)}^{\mathrm{2}}\right\},\end{array}\end{array}$$

where $w\left(x_i\right)=\mathrm{1}/{\mathit{\sigma }}_x^{\mathrm{2}}$ and $w=\left(y_i\right)\mathrm{1}/{\mathit{\sigma }}_y^{\mathrm{2}}$ are the weight coefficients for x and y, respectively, and r is the correlation coefficient between x and y. x_i,adj and y_i,adj are adjusted values of x_i, y_i, which fulfil the requirement

$$\begin{array}{}\text{(A14)}& y_{i,\mathrm{adj}}={\widehat{\mathit{\alpha }}}_{\mathrm{york}}+{\widehat{\mathit{\beta }}}_{\mathrm{york}}{x}_{i,\mathrm{adj}}.\end{array}$$

The solution for ${\widehat{\mathit{\alpha }}}_{\mathrm{york}}$ and ${\widehat{\mathit{\beta }}}_{\mathrm{york}}$ is found iteratively following the 10-step algorithm presented in York et al. (2004, Section III).

A5 The principal component analysis-based regression (PCA)

PCA can be applied for bivariate and multivariate cases.

For one independent and one dependent variable, the regression line is

$y={\widehat{\mathit{\alpha }}}_{\mathrm{PCA}}+{\widehat{\mathit{\beta }}}_{\mathrm{PCA}}x$ where the error between the observed value y_i and the estimated value a+bx_i is minimum. For n data points, we compute a and b using the method of least squares that minimises

$$\begin{array}{}\text{(A15)}& C_{\mathrm{PCA}}={\sum }_{i=\mathrm{1}}^N{\left(y_i-{\widehat{\mathit{\alpha }}}_{\mathrm{PCA}}-{\widehat{\mathit{\beta }}}_{\mathrm{PCA}}{x}_i\right)}^{\mathrm{2}}.\end{array}$$

This is a standard technique that gives the regression coefficients α and β.

$$\begin{array}{}\text{(A16)}& \left[\begin{array}{c}{\widehat{\mathit{\alpha }}}_{\mathrm{PCA}}\\ {\widehat{\mathit{\beta }}}_{\mathrm{PCA}}\end{array}\right]={\displaystyle \frac{\left[\begin{array}{cc}{S}_x& -\stackrel{\mathrm{‾}}{x}\\ -\stackrel{\mathrm{‾}}{x}& \mathrm{1}\end{array}\right]}{S_x-{\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}}}\left[\begin{array}{c}\stackrel{\mathrm{‾}}{y}\\ S_{xy}\end{array}\right]\end{array}$$

A6 Bayesian error-in-variables regression (Bayes EIV)

Bayes EIV regression estimate applies Bayesian inference using the popular Stan software tool (http://mc-stan.org/users/documentation/, last access: 27 July 2018), which allows the use of prior information of the model parameters. We assumed

$$\begin{array}{c}{\displaystyle }{\mathit{\beta }}_{\mathrm{BEIV}}\sim \mathrm{student}\mathit{\_}t(\mathrm{5},\mathrm{0.0},\mathrm{100.0})\\ {\displaystyle }{\mathit{\alpha }}_{\mathrm{BEIV}}\sim \mathrm{student}\mathit{\_}t(\mathrm{5},\mathrm{0.0},\mathrm{100.0})\\ {\displaystyle }{x}_{\mathrm{true}}\sim \log \mathrm{normal}({\mathit{\mu }}_x,{\mathit{\sigma }}_x)\\ {\displaystyle }{y}_{\mathrm{true}}={\mathrm{10.0}}^{\left({\mathit{\alpha }}_{\mathrm{BEIV}}+{\mathit{\beta }}_{\mathrm{BEIV}}\times {\log }_{\mathrm{10}}\left(x_{\mathrm{true}}\right)\right)},\end{array}$$

where μ and σ are the respective mean and standard deviation of x_true and y_true and are treated as unknowns. The observations x_obs and y_obs of x_true and y_true, respectively, were defined as follows:

$$\begin{array}{c}{\displaystyle }{x}_{\mathrm{obs}}\sim \mathrm{normal}\left(x_{\mathrm{true}},{\mathit{\sigma }}_{\mathrm{rel},x}\times x_{\mathrm{true}}+{\mathit{\sigma }}_{\mathrm{abs},x}\right)\\ {\displaystyle }{y}_{\mathrm{obs}}\sim \mathrm{normal}\left(y_{\mathrm{true}},{\mathit{\sigma }}_{\mathrm{rel},y}\times y_{\mathrm{true}}+{\mathit{\sigma }}_{\mathrm{abs},y}\right),\end{array}$$

where σ_rel and σ_abs are the relative and absolute components of standard uncertainties, respectively.

The Stan tool solved regression problems using 1000 iterations, and it provided a posteriori distributions for the model parameters β_BEIV and α_BEIV. For the definitions of given Student t, log-normal and normal probability distributions, see Stan documentation. In our regression analysis, we used the maximum a posteriori estimates for β_BEIV and α_BEIV provided by the software tool.