# Estimate error on slope of linear regression given data with associated uncertainty

Given a set of data, is it possible to create a linear regression which has a slope error that takes into account the uncertainty of the data?

This is for a high school class, and so the normal approach to find the uncertainty of the slope of the linear regression is to find the line that passes through the first data point minus its uncertainty and the last data point plus its uncertainty, and vice versa. Then, the slope of the line with the greater slope is subtracted from the other slope. However, this is not very accurate. Is there another way?

Both the x-coordinate and y-coordinate has an associated error. However, the error in the x-coordinate can be safely ignored without loss of marks. I would prefer a solution that takes into account both errors, but one that takes into account only the error in the y-coordinate is acceptable.

-
Welcome to Mathematica.SE! Can you please confirm that you plan to do this analysis in Mathematica? Otherwise I think it would be better to migrate your question to Cross Validated. –  Verbeia Oct 15 '12 at 0:55
It is commendable that you plan on using Mathematica for a high school project/class. –  drN Oct 15 '12 at 0:56
@J.M. Almost, there is an additional uncertainty on the x variable, so $(x_k, y_k, \text{xerr}_k, \text{yerr}_k)$, becoming $(x_k \pm \text{xerr}_k, y_k \pm \text{yerr}_k)$. –  George S Oct 15 '12 at 2:13
I'm not saying it's required; you're the only one who is supposed to determine whether they should be accounted for or ignored. My point was that weighted orthogonal regression (which accounts for errors in both coordinates) is a tougher problem to solve (and thus requires more elaborate methods) than weighted linear regression, which is easily handled by the built-in function LinearModelFit[] (via its Weights option). –  Ｊ. Ｍ. Oct 15 '12 at 2:52

Here's a method for doing weighted orthogonal regression of a straight line, based on the formulae in Krystek/Anton and York:

ortlinfit[data_?MatrixQ, errs_?MatrixQ] :=
Module[{n = Length[data], c, ct, dk, dm, k, m, p, s, st, ul, vl, w, wt, xm, ym},
(* yes, I know I could have used FindFit[] for this... *)
{ct, st, k} = Flatten[MapAt[Normalize[{1, #}] &,
NArgMin[Norm[Function[{x, y}, y - \[FormalM] x - \[FormalK]] @@@ data],
{\[FormalM], \[FormalK]}], 1]];
(* find orthogonal regression coefficients *)
{c, s, p} = FindArgMin[{
Total[(data.{-\[FormalS], \[FormalC]} -
\[FormalP])^2/((errs^2).{\[FormalS]^2, \[FormalC]^2})],
\[FormalC]^2 + \[FormalS]^2 == 1},
{{\[FormalC], ct}, {\[FormalS], st}, {\[FormalP], k/ct}}];
(* slope and intercept *)
{m, k} = {s, p}/c;
wt = 1/errs^2; w = (Times @@@ wt)/(wt.{1, m^2});
{xm, ym} = w.data/Total[w];
ul = data[[All, 1]] - xm; vl = data[[All, 2]] - ym;
(* uncertainties in slope and intercept *)
dm = w.(m ul - vl)^2/w.ul^2/(n - 2);
dk = dm (w.data[[All, 1]]^2/Total[w]);
{Function[\[FormalX], Evaluate[{m, k}.{\[FormalX], 1}]], Sqrt[{dm, dk}]}] /;
Dimensions[data] === Dimensions[errs]


ortlinfit[] expects data to contain the $(x_j,y_j)$ pairs, and errs to contain the corresponding uncertainties $(\rm{dx}_j,\rm{dy}_j)$. The routine returns the best-fit line as a pure function, as well as the uncertainties in the slope and intercept ($\sigma_m$ and $\sigma_k$).

As an example, here's some data used in York's paper:

data = {{0, 5.9}, {0.9, 5.4}, {1.8, 4.4}, {2.6, 4.6}, {3.3, 3.5}, {4.4, 3.7},
{5.2, 2.8}, {6.1, 2.8}, {6.5, 2.4}, {7.4, 1.5}};

errs = {{1000., 1.}, {1000., 1.8}, {500., 4.}, {800., 8.}, {200., 20.},
{80., 20.}, {60., 70.}, {20., 70.}, {1.8, 100.}, {1, 500.}} // Sqrt[1/#] &;

{lin, {sm, sk}} = ortlinfit[data, errs]
{Function[x, 5.47991 - 0.480533 x], {0.0710065, 0.361871}}


Now, let's look at the data, the associated error ellipses (constructed from the uncertainties), the best-fit line, and the "bounding lines" $(m-\sigma_m)x+(k-\sigma_k)$ and $(m+\sigma_m)x+(k+\sigma_k)$:

Show[
Frame -> True],
Plot[{lin[x], lin[x] - sm x - sk, lin[x] + sm x + sk}, {x, -1, 9},
PlotStyle -> {Directive[Thick, Red],
Directive[Dashed, Gray], Directive[Dashed, Gray]}]
]


-
At the bottom of your example you call the function with an undefined wts (presumably errs); also I get a different output : {Function[\[FormalX]$, 5.76267 - 0.536842 \[FormalX]$], {0.0468427, 0.198548}}. My approach gives a not too bad solution : {a -> -0.539577, b -> 5.76119} for the parameters and {0.0609037, 0.206709} for their errors. –  b.gatessucks Oct 16 '12 at 8:17
Whoops, thanks for catching that typo! (That was from a previous test, which I forgot to fix here.) Please try the code now. –  Ｊ. Ｍ. Oct 16 '12 at 8:48
A big +1. Now you got me wondering, how this problem could be solved for arbitrary curves (so for non-linear regression) –  Ajasja Oct 16 '12 at 9:48
@Ajasja:: the unweighted case of nonlinear orthogonal regression is already a bit difficult; more so if we have to account for weights/uncertainties. I'd say that should be for a different question... :) (However, the algorithm for unweighted orthogonal regression of a straight line is much simpler than the one I gave here; all that's needed is a clever application of the SVD.) –  Ｊ. Ｍ. Oct 16 '12 at 9:55

### Summary

Use the Weights option to LinearModelFit, setting the weights to be inversely proportional to the variances of the error terms.

### Theory

This is a standard problem: when the errors in the individual $y$ values are expressed in a way that can be related to their variances, then use weighted least squares with the reciprocal variances as the weights. (Search our sister site Cross Validated for more about this, including references and generalizations.)

### Creating realistic data to study

To illustrate, suppose the data are given as vectors $x$ and $y$ with the "errors" expressed either as standard deviations of $y$ or as standard errors of estimate of $y$, or any other quantity that can be interpreted as a fixed, constant multiple of the standard deviations of the $y$. Specifically, the applicable model for these data is

$$y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$$

where $\beta_0$ (the intercept) and $\beta_1$ (the slope) are constants to be estimated, the $\varepsilon_i$ are independent random deviations with mean zero, and $\text{Var}(\varepsilon_i) = \sigma_i^2$ for some given quantities $\sigma_i$ assumed to be known accurately. (The case where all the $\sigma_i$ equal a common unknown constant is the usual linear regression setting.)

In Mathematica we can simulate such data with random number generation. Let's wrap this into a function whose arguments are the amount of data and the slope and intercept. I will make the sizes of the errors vary randomly, but generally they will increase as $x$ increases.

simulate[n_Integer, intercept_: 0, slope_: 0] /; n >= 1 :=
Module[{x, y, errors, sim},
x = Range[n];
errors = RandomReal[GammaDistribution[n, #/(10 n)]] & /@ x;
y = RandomReal[NormalDistribution[intercept + slope  x[[#]],  errors[[#]]]] & / Range[n];
sim["x"] = x;
sim["y"] = y;
sim["errors"] = errors;
sim
]


Here is a tiny example of its use.

SeedRandom[17];
{x, y, errors} = simulate[16, 50, -2/3][#] & /@ {"x", "y", "errors"};
ListPlot[{y, y + errors, y - errors}, Joined -> {False, True, True},
PlotStyle -> {PointSize[0.015], Thick, Thick},
AxesOrigin -> {0, Min[y - errors]}]


The simulated points are surrounded by error bands.

### Weighted least-squares estimation

To fit these data, use the Weights option of LinearModelFit. Once again, let's prepare for later analysis by encapsulating the fitting in a function. For comparison, let's fit the data both with and without the weights.

trial[n_Integer: 1, intercept_: 0, slope_: 0] :=
Module[{x, y, errors, t, fit, fit0},
{x, y, errors} = simulate[n, intercept, slope][#] & /@ {"x", "y", "errors"};
fit = LinearModelFit[y, t, t, Weights -> 1 / errors^2];
fit0 = LinearModelFit[y, t, t];
{fit[#], fit0[#]} & @ "BestFitParameters"
]


The output is a list whose elements give {intercept, slope}: the first element is for the weighted fit, the second for the unweighted.

### Monte-Carlo comparison of the weighted and ordinary least squares methods

Let's run a lot of independent trials (say, $1000$ of them for simulated datasets of $n=32$ points each) and compare their results:

SeedRandom[17];
simulation = ParallelTable[trial[32, 20, -1/2], {i, 1, 1000}];

ranges = {{18.5, 22}, {-0.65, -0.35}};
TableForm[
Table[Histogram[simulation[[All, i, j]], ImageSize -> 250,
PlotRange -> {ranges[[j]], Automatic}, Axes -> {True, False}],
{i, 1, 2}, {j, 1, 2}],
TableHeadings -> {{"Weighted", "OLS"}, {"Intercept", "Slope"}}
]


Because I specified an intercept of $20$ and slope of $-1/2$, we will want to use these values as references. Indeed, the histograms in the left column ("Intercept") display sets of estimated intercepts hovering around $20$ and the histograms in the right column ("Slope") display sets of slopes hovering around $-0.50 = -1/2$. This illustrates the theoretical fact that the estimates in either case are unbiased. However, looking more closely at the spreads of the histograms (read the numbers on the horizontal axes), we see that those in the upper row ("Weighted") have smaller widths of their counterparts in the lower row ("OLS," or "Ordinary Least Squares"). This is evidence that the weighted estimates tend, on the whole, to be better then the unweighted ones, because they tend to deviate less from the true parameter values.

When the underlying data truly conform to the hypothesized model--there is a linear relationship between the $x$'s and $y$'s, with errors in the $y$'s having known but differing standard deviations--then among all unbiased linear estimates of the slope and intercept, weighted least squares using reciprocal variances for the weights is best in the sense just illustrated.

### Obtaining the estimation error of the slope

Now, to answer the question: we would like to assess the estimation error in the slope. This can be obtained from the fit object in many ways: consult the help page for details. Here is a nicely formatted table:

fit = LinearModelFit[y, t, t, Weights -> 1 / errors^2];
fit0 = LinearModelFit[y, t, t];
TableForm[{{fit[#], fit0[#]} & @ "ParameterConfidenceIntervalTable"},


In this case, for this particular set of simulated data (as created previously), the weighted method reports a much smaller standard error for the intercept than the OLS method (because errors near $x=0$ are low according to the information in errors) but the weighted estimate of the slope has only a slightly smaller standard error than the OLS estimate of the slope.

Errors in both $x$ and $y$ can be handled, using--for instance--methods of maximum likelihood. However, this involves considerably more mathematical, statistical, and computational machinery and requires a careful assessment of the nature of those errors (such as whether the $x$ errors and $y$ errors are independent). One general result in the statistical literature is that when the $x$ errors are typically smaller than the $y$ errors, yet independent of them, it is usually safe to ignore the $x$ errors. For more about all this, good search terms include "errors-in-variables regression," "Deming regression," and even "principal components analysis (PCA)".

-

I made this implementation of York's classical (and easy to understand) method following this paper by Cameron Reed.

f[x_List, y_List, wxi_List, wyi_List] :=
Module[{n = Length@x, wi, ui, vi, wmean, d, g, a, b, set, least},

wi[i_, m_]        := wxi[[i]] wyi[[i]]/(m ^2 wyi[[i]] + wxi[[i]]);
ui[i_, m_]        := x[[i]] - wmean[x, m];
vi[i_, m_]        := y[[i]] - wmean[y, m];
wmean[q_List, m_] :=  Sum[wi[i, m] q[[i]], {i, n}]/Sum[wi[i, m], {i, n}];
d[m_]             :=  Sum[wi[i, m]^2 ui[i, m]^2/wxi[[i]], {i, n}];
g[m_]             :=- Sum[wi[i, m]   ui[i, m] vi[i, m], {i, n}]/d[m];
a[m_]             :=2 Sum[wi[i, m]^2 ui[i, m] vi[i, m]/wxi[[i]], {i, n}]/(3 d[m]);
b[m_]             := (Sum[wi[i, m]^2 vi[i, m]^2/wxi[[i]], {i, n}] -
Sum[wi[i, m] ui[i, m]^2, {i, n}])/(3 d[m]);

set = {ToRules@ Reduce[m^3 - 3 a[m] m m + 3 b[m] m - g[m] == 0 &&
c == wmean[y, m] - m wmean[x, m], {m, c},
Backsubstitution -> True]};

least = Sum[wxi[[i]] (x[[i]] - (y[[i]] - c)/m)^2 +
wyi[[i]] (y[[i]] - (m x[[i]] + c))^2, {i, Length@x}] /.
set[[Flatten@Position[m /. set, _Real]]];

set[[Flatten@Position[m /. set, _Real]]][[Position[least, Min@least][[1]]]]];


Usage

f[Range[10], 3 Range[10] + RandomReal[.2], Array[# &, 10], Array[# &, 10]]
(*
-> {{m -> 3., c -> 0.110805}}
*)

-
This solution doesn't seem to work for all data. Changing the x error list to "{1, 2, 3, 4, 5, 6, 7, 8, 9, 2}" makes f return an error and a long, weird solution. Pastebin: pastebin.com/dMzTjweu –  George S Oct 16 '12 at 3:40
@GeorgeS I forgot to filter out imaginary results. Fixed now –  belisarius Oct 16 '12 at 3:58
I still receive the same error. The data I am using: f[{2.04, 3.06, 4.08, 5.1, 6.12, 7.14, 8.16}, {1.9975, 3.216, 4.0939, 4.9878, 6.4685, 7.2003, 8.2944}, {0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08}, {0.1696, 0.1793, 0.1821, 0.2568, 0.4197, 0.1342, 0.2304}] –  George S Oct 16 '12 at 5:11
I'm probably missing something, but shouldn't it be possible to obtain the error estimates for m and c? –  Ajasja Oct 16 '12 at 9:43
@Ajasja, yes, the formulae are in York's paper... –  Ｊ. Ｍ. Oct 16 '12 at 10:00

I'd take the simplest approach : solve the usual least-squares problem to determine the coefficients of the line; they will depend on the input data $(x_i, y_i)$. If the input data has an uncertainty $(dx_i, dy_i)$ then we can propagate it to the solution for the coefficients. We can do all this symbolically and substitute numerical values at the very end.

nPoints = 10;
data = Table[{Subscript[x, i], Subscript[y, i]}, {i, nPoints}];
errors = Table[{Subscript[dx, i], Subscript[dy, i]}, {i, nPoints}];

model[a_, b_, x_] = a x + b;

(* The least-squares functional; can be different, i.e. actual distance to the line *)
objFun[a_, b_, data_] := Total[(#[[2]] - model[a, b, #[[1]]])^2 & /@ data]

(* The usual solution *)
solab = First@Solve[{D[objFun[a, b, data], a] == 0,
D[objFun[a, b, data], b] == 0}, {a, b}] ;

(* The squared deltas relative to the input data *)
deltaa = Flatten@{D[a /. solab, #]^2 & /@ data[[All, 1]] ,
D[a /. solab, #]^2 & /@ data[[All, 2]]};
deltab = Flatten@{D[b /. solab, #]^2 & /@ data[[All, 1]] ,
D[b /. solab, #]^2 & /@ data[[All, 2]]};

(* The error is the sum of delta times uncertainty, assuming independence *)
errora = Sqrt[Dot[deltaa, Flatten[errors]^2]];
errorb = Sqrt[Dot[deltab, Flatten[errors]^2]];


Now we can use our actual numeric data and substitute them everywhere in the above :

ndata = Sort@RandomReal[{-2, 2}, {nPoints, 2}];
nerrors = RandomReal[{0, 0.2}, {nPoints, 2}];

nmodel[x_] = model[a, b, x] /. (solab /. Thread[Rule[Flatten[data] , Flatten[ndata]]])
nsolab = solab /. Thread[Rule[Flatten[data] , Flatten[ndata]]]
nerrorab = {errora, errorb} //. Join[Thread[Rule[Flatten[data] , Flatten[ndata]]],

(* -0.481409 + 0.125842 x
{a -> 0.125842, b -> -0.481409}
{0.0463579, 0.0385911} *)


Finally we can plot the data together with the best fit and the model corresponding to changing the parameters by $\pm$ 1 standard deviation.

Needs["ErrorBarPlots"]

Show[Plot[{model[(a /. nsolab) - nerrorab[[1]], (b /. nsolab) - nerrorab[[2]], x],
model[(a /. nsolab) - nerrorab[[1]], (b /. nsolab) + nerrorab[[2]],x],
model[(a /. nsolab) + nerrorab[[1]], (b /. nsolab) - nerrorab[[2]], x],
model[(a /. nsolab) + nerrorab[[1]], (b /. nsolab) + nerrorab[[2]],x],
nmodel[x]}, {x, -2, 2}, PlotStyle -> {Dashed, Dashed, Dotted, Dotted, Red}, PlotRange -> All],
ErrorListPlot[Transpose[{ndata, ErrorBar @@@ nerrors}]]]


-
This seems to be conceptually contradictory: the OLS solution assumes one error structure, then you propagate different errors. Why should that work at all? –  whuber Oct 15 '12 at 15:46
@b.gatessucks I seem to be unable to run this example properly. Running it when the random data your code uses, or pasting in my data returns a lot of errors. After a long delay, a graph with all the points but no lines of best fit appears. Pastebin: pastebin.com/xG9R79UN –  George S Oct 16 '12 at 1:51
@GeorgeS You need a ; after the definition of ndata or use different cells for ndata and nerrors. Also you need to adjust the plot range. –  b.gatessucks Oct 16 '12 at 8:08

Here is a simplistic approach, but perhaps it is one step on from examining just the end points. It makes the assumption that nothing is known about the distribution of x, y errors other than they are uniform, if they are not uniformly distributed the problem becomes significantly more complex.

Data

Some experimental data:

data = {Range@(size + 1),Sin@Range[0, (2 \[Pi])/3, (2 \[Pi])/(3 size)]}//Transpose;


Errors

A vector of {xErr,yErr} pairs which represent the estimated, or measured, magnitude of mean error for x and y measurements at each value of x:

errorDeltas = ConstantArray[{.1, 0.2}, Length@data];


Here the values of mean x and y errors have been chosen to be constant, but from experimental data they might well be variable.

Uniform

The simplicity derived from the uniform assumption is that the errors are symmetric, so the aggregate effect, of rotation or translation, on the linear fit cancel out. This allows for a simple addition of the mean error vector to the data to represent a least squares solution.

model = LinearModelFit[data + errorDeltas, x, x]


Plotted

Show[ListLinePlot[data, PlotStyle -> Dashed],
Plot[model[x], {x, 1, Length@data}], PlotRange -> All]
`

-
I don't understand how one the method used; why would adding error to the data help? In addition, the solution doesn't estimate an error on the slope of the line, as is necessary. –  George S Oct 16 '12 at 2:05
@GeorgeS This is a refinement of the example method you gave in your question using end points where the errors were added to the data. If the statistical distribution of the individual errors is uniform, or symmetric, then most likely ( maximum likelihood ) error value is the mean value of the errors. Thus, the most likely linear regression is the data plus the most likely error. The difference between this linear regression and the linear regression of the data on it's own, gives the error at any point on the regression. I hope that explanation wasn't too confusing :) –  image_doctor Oct 16 '12 at 8:02