How Analytic Choices Can Affect the Extraction of Electromagnetic Form Factors from Elastic Electron Scattering Cross Section Data

TL;DR

This paper demonstrates how different analytic choices in regression, such as constraints and data binning, significantly influence the extraction of proton electromagnetic form factors from electron scattering data.

Contribution

It highlights the impact of analytic choices on regression outcomes in form factor extraction and emphasizes the need for physical or mathematical grounding.

Findings

Adding a theory-motivated constraint removes oscillations in the magnetic form factor.

Data binning can dramatically alter the regression results.

Analytic choices can lead to substantially different physical parameter estimates.

Abstract

Scientists often try to incorporate prior knowledge into their regression algorithms, such as a particular analytic behavior or a known value at a kinematic endpoint. Unfortunately, there is often no unique way to make use of this prior knowledge, and thus, different analytic choices can lead to very different regression results from the same set of data. To illustrate this point in the context of the proton electromagnetic form factors, we use the Mainz elastic data with its 1422 cross section points and 31 normalization parameters. Starting with a complex unbound non-linear regression, we will show how the addition of a single theory-motivated constraint removes an oscillation from the magnetic form factor and shifts the extracted proton charge radius. We then repeat both regressions using the same algorithm, but with a rebinned version of the Mainz dataset. These examples illustrate how analytic choices, such as the function that is being used or even the binning of the data, can dramatically affect the results of a complex regression. These results also demonstrate why it is critical when using regression algorithms to have either a physical model in mind or a firm mathematical basis

Tables

Table 1. Table 1: The 34 different combinations of the 31 normalization parameters, N j , found in Ref.  Bernauer  et al. ( 2014 ) which link the data together along with the number of data points and the Q 2 superscript 𝑄 2 Q^{2} range of each dataset.

Energy	Spec.	$normA$	$normB$	Points	$Q^2$ Range [GeV2]
180 MeV	B	N1	N3	106	0.0038 – 0.0129
	B	N1	N4	41	0.0101 – 0.0190
	A	N3	-	102	0.0112 – 0.0658
	B	N1	N5	19	0.0190 – 0.0295
	C	N2	N4	38	0.0421 – 0.0740
	C	N2	N5	17	0.0740 – 0.0834
315 MeV	B	N6	N9	104	0.0111 – 0.0489
	A	N7	N9	38	0.0430 – 0.1391
	A	N9	-	40	0.0479 – 0.1441
	C	N8	N9	62	0.1128 – 0.2131
450 MeV	B	N10	N13	77	0.0152 – 0.0572
	B	N10	N15	52	0.0572 – 0.1175
	A	N13	-	42	0.0586 – 0.2663
	B	N10	N14	17	0.0589 – 0.0851
	A	N11	N13	36	0.0670 – 0.2744
	C	N12	N15	50	0.2127 – 0.3767
	A	N14	-	2	0.2744 – 0.2744
585 MeV	B	N16	N18	41	0.0255 – 0.0433
	B	N16	N19	47	0.0433 – 0.1110
	A	N18	-	27	0.0590 – 0.0964
	B	N16	N20	21	0.0920 – 0.1845
	A	N19	-	37	0.0964 – 0.4222
	C	N17	N20	20	0.3340 – 0.5665
720 MeV	B	N21	N25	47	0.0711 – 0.1564
	A	N25	-	46	0.1835 – 0.6761
	C	N24	N26	28	0.6536 – 0.7603
	B	N23	N26	27	0.2011 – 0.2520
	C	N22	N26	37	0.4729 – 0.7474
	B	N21	N26	36	0.1294 – 0.2435
855 MeV	B	N27	N31	35	0.3263 – 0.4378
	C	N28	N31	31	0.7300 – 0.9772
	A	N29	N30	32	0.3069 – 0.5011
	A	N29	-	13	0.5274 – 0.7656
	B	N27	N29	54	0.0868 – 0.3263

Table 2. Table 2: The values of the polynomial terms for the unbound and bound regressions of the 1422 cross section points following the notation of Eq.  4 and  5 . If one wishes to interpret the charge form factor slope term ( i = 1 ) 𝑖 1 (i=1) in terms of charge radius using Eq.  2 , one finds the unbound fit gives a charge radius of 0.882 fm while the bound fit gives a charge radius of 0.854 fm.

	unbound		bound
i	$a_i^E$	$a_i^M$	$a_i^E$	$a_i^M$
1	-3.331	-2.523	-3.124	-2.800
2	13.05	-0.7081	8.821	5.188
3	-63.68	40.16	-25.74	-5.742
4	249.4	-176.7	60.06	2.806
5	-658.6	380.3	-89.41	-0.000
6	1099	-392.6	72.48	0.01034
7	-987.6	11.53	-24.23	-0.2766
8	57.38	442.4	0.0000	0.0000
9	853.4	-492.1	-0.0061	-0.0009
10	-810.5	230.3	0.0081	0.0013
11	250.4	-40.92	-0.0000	-0.0000

Table 3. Table 3: The values of the polynomial terms for the unbound and bound regressions of the 658 cross section points of rebinned data  Lee  et al. ( 2015 ) following the notation of Eq.  4 and  5 . If one wishes to interpret the charge form factor slope term ( i = 1 ) 𝑖 1 (i=1) in terms of charge radius using Eq.  2 , one finds the unbound fit gives a charge radius of 0.863 fm while the bound fit gives a radius of 0.845 fm.

	unbound		bound
i	$a_i^E$	$a_i^M$	$a_i^E$	$a_i^M$
1	-3.191	-2.465	-3.061	-2.760
2	10.83	-0.7271	8.413	4.979
3	-44.59	35.32	-24.46	-5.196
4	157.2	-136.4	58.23	2.193
5	-404.2	228.8	-89.36	-0.000
6	712.6	-98.11	74.77	0.5035
7	-733.1	-234.1	-25.73	-0.5330
8	133.4	349.9	0.0000	0.0000
9	632.8	-122.4	-0.0000	-0.0000
10	-695.5	-56.49	0.0000	0.0000
11	232.9	35.80	-0.0000	-0.0000

Table 4. Table 4: Four datasets of (x i ,y i ,dy i ) values adapted from  Anscombe ( 1973 ) .

x1,2,3	y1	y2	y3	x4	y4	dy1,2,3,4
10.0	8.04	9.14	7.46	8.0	6.58	1.235
8.0	6.95	8.14	6.77	8.0	5.76	1.235
13.0	7.58	8.74	12.74	8.0	7.71	1.235
9.0	8.81	8.77	7.11	8.0	8.84	1.235
11.0	8.33	9.26	7.81	8.0	8.47	1.235
14.0	9.96	8.10	8.84	8.0	7.04	1.235
6.0	7.24	6.13	6.08	8.0	5.25	1.235
4.0	4.26	3.10	5.39	19.0	12.50	1.235
12.0	10.84	9.13	8.15	8.0	5.56	1.235
7.0	4.82	7.26	6.42	8.0	7.91	1.235
5.0	5.68	4.74	5.73	8.0	6.89	1.235

Equations

$$\sigma =$$

$$\left[\frac{G_E^2\left(Q^2\right)+\tau G_M^2\left(Q^2\right)}{1+\tau }+2\tau G_M^2\left(Q^2\right){\tan }^2\left(\frac{\theta }{2}\right)\right]\text{,}$$

$$r_p\equiv {\left(-6{\frac{d{G}_E(Q^2)}{d{Q}^2}|}_{Q^2=0}\right)}^{1/2}.$$

$${\sigma }_{exp}={\sigma }_p\cdot norm{A}_p\cdot norm{B}_p,$$

$$G_{E,\mathrm{polynomial}}(a_i^E,Q^2)$$

$$G_{M,\mathrm{polynomial}}(a_i^M,Q^2)$$

$${\chi }^2=\sum _{p=1}^{p_{max}}{\left(\frac{{\sigma }_{Model}(E_p,{\theta }_p)-{\sigma }_p\cdot norm{A}_p\cdot norm{B}_p}{\Delta {\sigma }_p\cdot norm{A}_p\cdot norm{B}_p}\right)}^2\text{,}$$

$$G_{E,\mathrm{dipole}}(Q^2)$$

$$G_{M,\mathrm{dipole}}(Q^2)$$

$${\chi }^2$$

$${\chi }^2/\mathrm{df}$$

$$\mathrm{AIC}$$

$$\mathrm{BIC}$$

$$r_p=(-6a_1^E{)}^{1/2}\text{.}$$

$${\chi }^2\equiv \sum _{i=1}^N{\left(\frac{y_i-y\left(x_i|a_1,a_2,...,a_M\right)}{{\sigma }_i}\right)}^2$$

$${\chi }^2\equiv \sum _{i=1}^N{\rho }_i\left(z\right)\text{and}z={\left(\frac{y_i-y\left(x_i|a_1,a_2,...,a_M\right)}{{\sigma }_i}\right)}^2$$

$$\rho \left(z\right)=2\left(\sqrt{1+z}-1\right)$$

Full text

form factors; proton radius; confirmation bias; regression; robust methods

I Introduction

Silberzahn et al. Silberzahn et al. (2018) points out that there is often little appreciation for how different analytic strategies can affect a reported result. In this work, we illustrate how analytic choices can impact the extraction of the electromagnetic form factors and the associated charge radii from electron scattering data. These extractions are frequently done with complex non-linear regression algorithms and tend to make use of prior information about the limiting behavior of the electromagnetic form factors to help constrain the value of experimental normalization parameters. Also, while many tend to look at all regressions as being the same, in fact there are different types of regressions such as descriptive, predictive, and explanatory.

A descriptive model is used to capture the features of a dataset in a compact manner without reliance on an underlying theory. A predictive model is any statistical model which tries to generalize beyond the data that is being fitted. Finally, explanatory modeling takes a theory based model and tests that model’s hypothesis by applying it to data. Further details about these differences can be found in Ref. Shmueli (2010). Though the type of regression model being developed is not always clearly stated, it is yet another choice that affects how scientists design their regression algorithms.

II Proton Elastic Scattering

There has been renewed interest in proton elastic scattering data due to muonic hydrogen Lamb shift results which determined the charge radius of the proton to be 0.84078(39) fm Pohl et al. (2010); Antognini et al. (2013), a result in stark contrast to the CODATA-2014 recommended value of 0.8751(61) fm Mohr et al. (2016). This systematic difference was known as the proton radius puzzle Pohl et al. (2013); Carlson (2015); Miller (2018). We will show that determining the proton’s charge radius is highly dependent on the analytic choices made when selecting a model to describe the world data.

In the plane-wave Born approximation, the cross section for elastic electron scattering on a proton is given by:

[TABLE]

where $\sigma_{\mathrm{Mott}}$ is the Mott cross section, $G_{E}$ and $G_{M}$ are the electric and magnetic Sachs form factors respectively, $\tau=\frac{Q^{2}}{4m_{p}^{2}}$, $Q^{2}=4E_{\mathrm{Beam}}E^{\prime}\sin^{2}\left(\frac{\theta}{2}\right)$, $E_{\mathrm{Beam}}$ is the energy of the electron beam, $E^{\prime}$ is the energy of the outgoing electron, $\theta$ is the scattering angle of the outgoing electron, and $m_{p}$ is the mass of the proton.

The proton charge radius, $r_{p}$, is extracted from the cross sections by determining the slope of the electric form factor, $G_{E}$, in the limit of four-moment transfer, $Q^{2}$, approaching zero Miller (2019):

[TABLE]

Since the scattering data is measured at finite $Q^{2}$, an extrapolation is required to extract the charge radius. A purely mathematical fit to the scattering data would be a descriptive model and would generally only be valid in the region of the data making extrapolation risky. When extrapolating back to $Q^{2}$=0 it is desirable to use a predictive model. This requires extra care such as adding physics considerations to the model and/or mathematical requirements to keep the fit well behaved (e.g. not diverging in the region) and not unduly complex. Authors have taken many different approaches to this extraction, yielding different outcomes Horbatsch and Hessels (2016); Griffioen et al. (2016); Higinbotham et al. (2016); Lee et al. (2015); Graczyk and Juszczak (2014); Lorenz and Meißner (2014); Horbatsch et al. (2017); Alarcón et al. (2019, 2020); Zhou et al. (2019); Mihovilovič et al. (2020).

III Regressions

To illustrate how analytic choice can strongly affect the extracted radius, we first use the Mainz dataset of 1422 cross section points along with its 31 normalization parameters Bernauer et al. (2010). As noted in the work of Bernauer et al. Bernauer et al. (2014), knowledge of the absolute value of cross sections is limited by the determination of the absolute luminosity, which in turn is limited by the uncertainty of the target thickness and beam current. In order to compensate for these uncertainties, normalization parameters were introduced to the original fits of this data which were constrained only by the value of the charge and magnetic form factors in the limit of $Q^{2}=0$. While there is no debate about the value of the form factors at $Q^{2}=0$, how the data at finite $Q^{2}$ is connected to the known end-point brings a model dependence to the analysis that is not easily understood.

These parameters are taken in combinations to link sets of data together, with the final value of each cross section point defined by:

[TABLE]

where $normA_{p}$ and $normB_{p}$ are the two normalization parameters associated with that data point. A complete list of the 31 different normalization parameters, Nj, that are taken in 34 unique combinations for the 1422 points, is shown in Table 1. Further details of how these parameters connect to each of the 1422 cross section points can be found in the supplemental material of Bernauer et al. Bernauer et al. (2014).

For our example of how analytic choice can effect the outcome, we first fit the Mainz dataset with an unbound complex non-linear regression with the form factors parameterized in terms of polynomials:

[TABLE]

where $\mu_{p}$ is the magnetic moment of the proton and $n$ is the order of the polynomial. This is one of the many form factors parameterizations described in Ref. Bernauer et al. (2014) that has been used for extracting the proton radius from the Mainz dataset.

For these regressions, we perform a weighted least squares minimization with a $\chi^{2}$ function defined as follows:

[TABLE]

where for each data point $p$ there is a cross-section, $\sigma_{p}$, with energy $E_{p}$, angle $\theta_{p}$, and normalization parameters $normA$ and $normB$ as shown in Table 1. As was done in the Mainz fits, the normalization parameters are allowed to float freely.

We repeat this same regression adding one requirement: that the terms of the polynomial have successively alternating signs. This makes the polynomial more closely resemble a “completely monotone” function and is referred to as a bound regression. A true completely monotone function, $f$, would possess derivatives, $f^{n}$, of all orders such that $(-1)^{n}f^{n}(x)\geq 0,x>0$ Merkle (2012). Of course for a finite order polynomial being fit to experimental data, we are simply approximating a completely monotone function over the range of the data by alternating the signs of the terms. This seemingly simple constraint imposes an analytic behavior to the form factors that is constant with nuclear physics calculations such as chiral effective field theory Alarcón and Weiss (2018). In statistics terms, adding the condition that the polynomial approximate a physically motivated function would be classified as creating a robust regression model Press et al. (2007). Robust regression models are designed such that they are not unduly affected by outliers, whereas least squares estimates are highly sensitive to outlying points as illustrated in Appendix A.

As a further check of how sensitive these two functions are to the handling of the data in the fit, we use the rebinned version of the Mainz data that is provided in the supplemental material of Ref. Lee et al. (2015). These authors carefully rebinned and re-weighted the full Mainz dataset and provided a new set of 658 cross section points, though with the same 31 normalization parameters as the original set. By simply replacing the original Mainz dataset with this set, we can repeat our unbound and bound regressions.

While regressions that are linear in terms can be solved exactly, this is not the case with non-linear regressions where algorithms can converge in a local or non-physical minimum; thus choosing reasonable initialization parameters is an important step when developing non-linear regression algorithms. To have reasonable initialization parameters for our complex non-linear regressions, we first perform a regression with dipole functions for $G_{E}$ and $G_{M}$ and use the resultant normalization parameters as initialization parameters in the more complex regressions.

In Table 2 and Table 3 we show the results of fitting with both the unbound and bound regressions for the 1422 Mainz cross section points and the rebinned 658 Mainz cross section points respectively. The regression results and residuals are shown graphically in Fig. 1 and Fig. 2 for the 1422 Mainz cross section points and the rebinned 658 Mainz cross section points respectively. For clarity we divide $\sigma_{exp}$ by $\sigma_{dipole}$, where $\sigma_{dipole}$ is simply Eq. 1 with standard dipole form factors:

[TABLE]

The fits of the unbound and bound regressions clearly differ significantly for both the original and the rebinned Mainz cross section points, but the residuals of each fit are quite similar. The reason the locations of the cross section data points differ between the unbound and bound regressions in Fig. 1 and Fig. 2 is due to the choice of regression model shifting the 31 normalization parameters to maintain agreement with our prior knowledge of the values of the electromagnetic form factors in the limit of $Q^{2}$ = 0. The magnitude of individual normalizaton shifts is a few tenths of a percent, which is much smaller than absolute normalizatons can be determined, but can have a clear effect on the results as the point-to-point uncertainties are also just few tenths of a percent.

IV Model Selection

For a fixed number of fit parameters, the unbound regressions presented in this work will always have a total $\chi^{2}$ equal to or lower than a bound regression as shown in Fig. 3 and Fig. 4 which plot the total $\chi^{2}$ and charge radius given by the unbound and bound regressions of the 1422 Mainz cross section points and the 658 rebinned cross sections respectively. Since adding parameters will always either decrease or keep total $\chi^{2}$ the same, $\chi^{2}$ by itself is not a valid model selection criterion. More appropriate model selection techniques include using an F-Test for nested models that are linear in terms Higinbotham et al. (2016); Širca (2016) or model selection methods like the Akaike Information Criterion (AIC) Akaike (1974) or the Bayesian Information Criterion (BIC) Schwarz (1978) which can be used with non-nested non-linear models (see Ref. Higinbotham et al. (2018) for more details).

Since the regressions herein are non-linear, we use the AIC and BIC to determine the most appropriate number of parameters for these regressions. These statistical criteria, along with the frequently quoted $\chi^{2}$ per degree of freedom $\chi^{2}/{\mathrm{df}}$ Andrae et al. (2010), are defined as follows:

[TABLE]

where $N$ is the number of data points, $data_{i}$ and $\sigma_{i}$ are the measured values and their estimated uncertainties respectively, $model_{i}$ is the model value, and $N_{\mathrm{var}}$ is the number of model parameters. Starting from lowest order fits, at first the value of these criteria will decrease as a parameter is added indicating an underfitting of the data while eventually the criteria will start to increase as parameters are added indicating an overfitting of the data; thus, with these criteria the model with the lowest AIC or BIC value should be selected as most appropriate.

For the 1422 original Mainz data points, we find with both AIC and BIC the most appropriate of the bound fits is the 7th order fit with a $\chi^{2}/{\mathrm{df}}$ of 1.21, while for the unbound fits the 10th order fit with a $\chi^{2}/{\mathrm{df}}$ of 1.14 is most appropriate as previously found Bernauer et al. (2010, 2014). For the 658 rebinned Mainz data points, we find with both AIC and BIC the most appropriate of the bound fits is the 7th order fit with a $\chi^{2}/{\mathrm{df}}$ of 0.865, while for the unbound fits the 9th order fit with a $\chi^{2}/{\mathrm{df}}$ of 0.830 is most appropriate. Ideally a criteria would have been defined prior to the analysis of the data, since we are doing a re-analysis, we are presenting two of the most common techniques so one can judge how even the criteria can effect the outcome of the analysis. Further details about model selection techniques can be found in Wit et al. (2012).

One should also keep in mind whether one is trying to do a descriptive fit of the data or, by adding physical constraints, building a predictive or explanatory model of the data Shmueli (2010). One must also keep in mind that none of these model selection techniques will prevent the use of completely inappropriate functions nor do they ensure that the best type of function has been selected. As noted in Ref. Anscombe (1973); Andrae et al. (2010), it is essential to plot the fit functions and residuals to ensure a reasonable regression as $\chi^{2}$ minimization alone is insufficient as illustrated in Appendix B.

As the changes we have presented in these four fits are larger than the statistical parameter uncertainties, we have limited ourselves to a discussion of the shifts of the mean values of the points. For non-linear regressions such as these, statistical bootstrapping which makes use of sampling with replacement can be used to find the statistical parameter uncertainties Press et al. (2007).

V Results

In Fig. 5 and Fig. 6, we show the individual electric and magnetic form factors obtained from the unbound and bound regressions for the 1422 Mainz cross section points. In Fig. 5 we see that $G_{M}$ for both the unbound and bound fits remains well behaved at high $Q^{2}$. However, for $G_{E}$ both the unbound and bound fits begin to diverge at high $Q^{2}$. This is due to the dominance of the magnetic form factor in the cross section at high $Q^{2}$ which leads the electric form factor to become unconstrained in this region. Due to using a high-order polynomial regression, along with the unconstrained nature of $G_{E}$ in the high $Q^{2}$ region, the divergence of $G_{E}$ is to be expected.

In Fig. 6 the ratios of the electromagnetic form factors to the standard dipole are shown. This ratio reveals that the unbound $G_{E}$ ratio has a small oscillation and the unbound $G_{M}$ ratio has a large oscillation. Whereas, the bound $G_{E}$ ratio has mostly removed the oscillation and the bound $G_{M}$ ratio has almost no oscillation. The unbound ratios are descriptive models of the scattering data without any physics considerations, but the bound ratios are more akin to predictive models as the terms alternate sign as one would expect from chiral effective field theory. By adding this one physical constraint the oscillations are removed as the model becomes more predictive.

Though it is beyond the range of the data used in the regression, the results of regressions like these are frequently used to extract the charge radius of the proton by using Eq. 2 to relate the fit function to the charge radius of the proton. Ideally, these extractions would use a predictive model to fit the data as finding the charge radius requires extrapolating beyond the range of the experimental data. For the case of a polynomial regression, this is simply:

[TABLE]

Using Eq. 2 one finds a charge radius of 0.882 fm from the unbound regression and 0.854 fm from the bound regression of the original Mainz cross section points (Table 2), and a charge radius of 0.863 fm from the unbound regression and 0.845 fm from the bound regression of the rebinned Mainz data (Table 3). These results show that even just rebinning the data can shift the result of a high-order polynomial regression significantly.

In the end, the radii extracted from the more descriptive unbound regressions are closer to the CODATA-2014 value, while the radii extracted from the more physically justifiable and predictive bound regressions are closer to the muonic results as well as the most recent atomic result Bezginov et al. (2019). With freedom to make analytic choices that so strongly affect the results, there is the potential for unconscious confirmation bias, and for researchers to select and report the regressions that confirm their expectations Koehler (1993); Ioannidis (2005).

VI Conclusions

We have shown that small changes in analytic functions and binning choices applied to a complex non-linear regression can result in significantly different results. In particular, using the Mainz dataset of elastic cross section points to extract a proton charge radius, we have shown results consistent with the CODATA-2014 value when using high-order unbounded polynomial fits and values close to the muonic results when using bounded polynomial regressions. Thus, by simply trying different functions, limits, and bounds, one can easily extrapolate to different results which can lead to confirmation bias and/or inappropriate rejection of certain results. Enrico Fermi noted that these types of problems should be addressed using either a firm mathematical basis or a physical model Dyson (2004). Ideally, regression models should be carefully developed prior to obtaining experimental data, as was done by the PRad collaboration Yan et al. (2018); Xiong et al. (2019); otherwise, one must be exceedingly careful to avoid confirmation bias though the rigorous use of model selection techniques Koehler (1993); Ioannidis (2005).

Thus, while one can argue that the bounded non-linear regression is the more physical function, it would be more appropriate to approach the analysis such that the analytic choices do not so strongly affect the results. To do this, one can either fit only lower $Q^{2}$ data where fewer free parameters are required Hand et al. (1963); Murphy et al. (1974); Borkowski et al. (1975); Simon et al. (1980); Higinbotham et al. (2016); Griffioen et al. (2016); Yan et al. (2018); Hagelstein and Pascalutsa (2019); Higinbotham et al. (2018); Hayward and Griffioen (2018) and the results are not sensitive to the magnetic form factor, as shown explicitly in Ref. Higinbotham et al. (2018); or, as Fermi preferred, use a physical model, such as that of Bernard et al. Bernard et al. (1998) or Alarcón and Weiss Alarcón and Weiss (2018) to constrain the fits Hohler et al. (1976); Lorenz et al. (2015); Horbatsch et al. (2017); Alarcón et al. (2019). There are also the physically motivated functions such as rational functions Higinbotham et al. (2016); Kelly (2004), continued fractions Griffioen et al. (2016); Sick (2003); Lorenz et al. (2012), or the z-expansion fits Hill and Paz (2010); Lee et al. (2015) though these still require model selection techniques to determine the appropriate number of regression parameters. We hope to have illustrated that by using extremely complex non-linear regressions and deep searches, one can find nearly any radius in a wide range of radii from a single dataset Higinbotham (2019). To quote Nobel laureate Ronald H. Coase, “if you torture the data long enough, it will confess.”

Acknowledgements.

Many thanks to Franziska Hagelstein and Vladimir Pascalutsa for their questions about normalization parameters that prompted this work. Thanks to Dave Meekins for his many useful comments and discussions. Thanks to Nigel Tufnel for suggesting going to eleven. And thanks to Marcy Stutzman and Simon Širca for their editorial comments. This work was supported by the U.S. Department of Energy contract DE-AC05-06OR23177 under which Jefferson Science Associates operates the Thomas Jefferson National Accelerator Facility.

Appendix A Robust Regressions

Ordinary least squares regression (OLSR) is one of the most commonly employed techniques used to fit a given model and its parameters to a dataset. However, OLSR is commonly misunderstood and misapplied by researchers. For OLSR, fit parameters are determined by minimizing the square of the differences between real-world data and model predictions. This is known as a $\chi^{2}$ minimization. Eq. 14 shows this minimization for $N$ data points with $M$ fit parameters.

[TABLE]

Here $y_{i}$ are the measured data values, $y\left(x_{i}|a_{1},a_{2},...,a_{M}\right)$ are the values given by the model with fit parameters $a_{1}$ to $a_{M}$ when evaluated at the $x_{i}$ of the measured data, and $\sigma_{i}$ are the uncertainties on each measured data point.

While OLSR via $\chi^{2}$ minimization is often a useful initial method for checking the ‘goodness’ of a fit, it can fail if the dataset being fit does not meet certain conditions. OLSR is based on the core assumptions that the errors are random variables that are normally distributed, the errors are uncorrelated to each other, and the errors are homoscedastic, which is to say they have the same variance. Unfortunately, these assumptions often do not hold true in the case of real-world data. When OLSR’s assumptions are not met, such as when the dataset has significant outliers, OLSR is not sufficient for fitting the data and can yield misleading results. Even a singular outlier can skew the results of an OLSR pulling the fit away from the data’s true behavior Press et al. (2007); Draper (1988). To avoid these pitfalls, robust methods such as robust least squares regression (RLSR) should be used instead of OLSR techniques.

To avoid outliers having too much influence over a fit, we desire a method by which outliers can be identified and then re-weighted such that they do not skew the overall fit. The least squares minimization found in Eq. 14 can be generalized to Eq. 15 by introducing the function $\rho(z)$ Triggs et al. (2000). OLSR is then simply the case where $\rho(z)=z$. Many functions can be used for $\rho(z)$ to introduce robustness, but for the following examples the ‘soft loss’ (softl1) function given in Eq. 16 was selected and implemented using the Python package SciPy Jones et al. (01); Oliphant (2007); Triggs et al. (2000).

[TABLE]

[TABLE]

With soft loss, as a $z_{i}$ gets larger, the magnitude of $\rho_{i}(z)$ is increasingly reduced with respect to OLSR. A RLSR with soft loss essentially re-weights the outliers of a dataset, decreasing their influence when fitting. Note that if a dataset meets all of the above assumptions inherent to OLSR (i.e. errors are normally distributed, uncorrelated, and have the same variance) then OLSR and RLSR techniques should both produce the same fit results since the dataset, by definition, does not contain excessive outliers.

A simple example reproduced from a classic statistics paper Anscombe (1973) is shown in Fig. 7. The dataset has a single clear outlier which pulls the fit considerably away from the bulk of the data when using OLSR. However, when RLSR with soft loss is used to fit the data the influence of the outlier is greatly reduced, and the fit returns to the bulk of the data yielding superior results.

For a real-world example using RLSR we can study the full Initial State Radiation (ISR) dataset found in the supplemental material of Ref. Mihovilovič et al. (2019, 2017). Fig. 8 shows the results of two regressions of the ISR proton electric form factor data using the theory model of Alarcón and Weiss Alarcón and Weiss (2018), with the proton radius as its single free parameter. The light gray curve uses an OLSR to fit the dataset and the dark curve uses a RLSR with soft loss.

There is a clear separation between the results of the two regressions, with the RLSR fitting the higher $Q^{2}$ data better. The OLSR finds a proton radius of 0.874 fm, and the RLSR finds a significantly smaller radius of 0.844 fm. Again, the purely analytic choice of the regression type significantly influences the fit results. Further, the fact that the OLSR and the RLSR fit results differ significantly is evidence that there are outliers in the dataset that are not following a normal distribution, otherwise the OLSR and RLSR fits would have better agreement.

Appendix B Anscombe’s Quartet

With the power of modern computing, one can be tempted to blindly assume the results of the regression are correct if the $\chi^{2}/{\mathrm{df}}$ is close to unity; but this can be extremely misleading Širca and Horvat (2012). To illustrate this point we use the Anscombe quartet Anscombe (1973), though as nuclear physicists tend to use $\chi^{2}/{\mathrm{df}}$ instead of R2. We have taken the 1973 example problem and added uncertainties to the points as shown in Table 4.

When fit with a linear function, these four sets of (x,y,dy) values give the same statistical quantities to three significant figures: mean, variance, $\chi^{2}$, $\chi^{2}/{\mathrm{df}}$, etc. So if one fails to make graphical checks, one can be completely fooled into thinking the fits are all equally good; but by simply graphing (see Fig. 9) one can see that only the first set of data is distributed in an ideal way around the fit function.

Dataset two clearly has a curved residual yet has exactly the same mean, errors, and $\chi^{2}$ as the first fit. This suggests that the fitter should likely add a quadratic term to the regression as well as check that the uncertainties have been correctly reported.

Dataset three illustrates the effect of an outlier on the regression. Of course, the scientist doesn’t simply throw out an outlier. Instead one should report on the outlier’s influence on the result. For example, in dataset three it would be worth noting that if the outlier is removed, the data exactly follow a line of $y=4+0.346x$ and that that measurement should be repeated.

Dataset four looks unsatisfactory, since all the information about the slope comes from one observation. This is very different from dataset one where any one point can be removed and one will obtain nearly the same result. Thus, for dataset four it should be pointed out that a single observation plays a critical role in the result.