1
$\begingroup$

Basically I have some datapoints on the form $(x_i,y_i)$

data = {{1, -12.3989}, {2, -13.3147}, {3, -14.528}, {4, -15.2962}, {5,
-16.3706}, {6, -16.9926}, {7, -17.7229}, {8, -18.603}, {9, -18.8811},
{10, -19.6465}, {11, -20.5828}, {12, -20.9566}, {13, -20.8691}, {14,
-21.695}, {15, -22.5606}, {16, -22.1918}, {17, -22.8421}, {18,
-23.3825}, {19, -23.6358}, {20, -23.8957}}

I also have the variances of $y_i$ in a vector

var = {0.30746, 3.76794, 11.1175, 15.7907, 21.2136, 23.0218, 25.8991,
26.9632, 28.9219, 29.369, 34.2043, 36.2889, 33.9908, 36.3291,
39.5869, 31.1223, 32.9673, 34.1042, 38.1544, 39.4407}

I expect the relationship between $x_i$ and $y_i$ to be linear. To find the line of best fit I write

nlm = NonlinearModelFit[data, a x + b, {a, b}, {x}, 
Weights -> 1/error^2]

I can get the standard error in the slope by writing

nlm[{"BestFit", "ParameterTable"}]

and obtain that the standard error in the slope is is

0.0261997

However, this error seems a bit small to me. Is there a better way to estimate the error in the slope that is easily implemented in mathematica? For instance is there a simple way to find the maximum and minimum possible lines and say that the error is the difference divided by two?

$\endgroup$
  • 3
    $\begingroup$ What is the value of error in nlm? Any information about the distribution of $y$? $\endgroup$ – MikeY May 9 at 16:50
  • $\begingroup$ The error I obtain is in the question 0.0261997. The distribution of y is unknown $\endgroup$ – MOOSE May 9 at 17:13
  • $\begingroup$ After plotting +/- 2 standard deviations (based on the values of var), the variability associated with data seems way, way too low to have independent errors. Either the values in var are way to high or there is likely a large amount of serial correlation among the errors. $\endgroup$ – JimB May 9 at 17:16
  • $\begingroup$ There is indeed a trend in the variance. But at the moment I'm just concerned with finding a decent measure for the uncertainty in the slope of the line $\endgroup$ – MOOSE May 9 at 17:52
  • $\begingroup$ Who said anything about a trend in the variance? You likely have serial correlation in the errors about the line which violates the assumption of independent errors. $\endgroup$ – JimB May 9 at 18:46
1
$\begingroup$

This is straight from the manual:

Using VarianceEstimatorFunction->(1&) and Weights->{1/Δy1^2,1/Δy2^2,…}, Δyi is treated as the known uncertainty of measurement yi, and parameter standard errors are effectively computed only from the weights.

So you need to do

nlm = NonlinearModelFit[data, a x + b, {a, b}, x, 
  VarianceEstimatorFunction -> (1 &), Weights -> 1/var]

and you get a standard error on the slope of 0.125726:

nlm[{"BestFit", "ParameterTable"}]

enter image description here

In this linear case you can also use LinearModelFit:

lm = LinearModelFit[data, {1, x}, x, 
  VarianceEstimatorFunction -> (1 &), Weights -> 1/var]
(* same results as nlm above *)

To see if the results make sense:

P1 = Plot[Evaluate@lm["MeanPredictionBands"], {x, Min[data[[All, 1]]], Max[data[[All, 1]]]}];
P2 = ListPlot[MapThread[{#1[[1]], Around[#1[[2]], Sqrt[#2]]} &, {data, var}]];
Show[P1, P2]

enter image description here

$\endgroup$
1
$\begingroup$

If all of the regressions assumptions are met (or at least not greatly violated) the answer that @Roman gives is exactly what you want.

However, your data and/or the estimated variances for each point grossly violates some of the regression assumptions. As seen from the plot of the data and fit from @Roman's answer the error bars (for each data point) far exceeds the variability one sees in the plot of the data points.

While there are many causes of this behavior, two of the causes are (1) lack of independence among the errors about the line (maybe some sort of serial correlation among the errors about the line) and (2) the variances in var are wrong.

Here's one way to check on the assumption of independent errors. Suppose that the estimated fit is correct and the variances in var are correct. Then, if the errors for each data point were independent, we would see the following kind of data:

new = data;
SeedRandom[12345]; 
new[[All, 2]] = (nlm[#] & /@ data[[All, 1]]) + 
  (RandomVariate[NormalDistribution[0, #], 1][[1]] & /@ (var^0.5));
ListPlot[{data, new}, PlotLegends -> {"Original data", "Simulated data"}]

Original and simulated data

We see far more variability in the simulated data than in the original data. This implies that there is a least one major regression assumption violation.

While the estimates of the slope and intercept won't change much however this issue is reconciled, the estimates of the associated precision of the estimates can be wildly different (which is the basis of your question).

If you're doing something important with the results of the regression analysis (or if your inferences might be challenged in court), I strongly recommend reconciling the deviations from the standard regression assumptions. Possibly the way the data was collected induced serial correlation. (That's just a guess.) If so, then there are ways to deal with that and get more appropriate estimates of precision (for the predictions and the estimated parameters). Just not with NonlinearModelFit or LinearModelFit which assume independence of observations.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.