# Estimating the uncertainty in the slope of a straight line with no associated uncertainty

I am just starting out with mathematica and cannot for the life of me find how to find the uncertainty in a slope. I am attempting to find the exponent in a logarithmic relationship between two variables, with no associated uncertainty.

M = Log /@ {70, 3.4, 21, 21, 600, 352, 17.5, 21.5, 148, 2620, 467000,
1820000, 22, 75, 88, 157, 125, 39.7, 2719000, 3672000, 407000,
57800}
R = Log /@ {50, 800, 168, 186, 33, 58, 140, 120, 55, 33, 27.4, 16,
135, 91, 65, 126, 54, 120, 4.9, 4.5, 16, 26}

I had little difficulty using this data to find the slope of best fit

data = Transpose@{M, R};
Fit[data, {1, x}, x];

returned

5.65474 - 0.242763 x

My hypothesis for this slope was -0.25, and so I would be tempted to just call this a success, however to be truly certain I need to know what the uncertainty is in this slope.

I have tried using a nonlinear model, however the best I managed to get was

Number of coordinates (1) is not equal to the number of variables (2).

when I ran the code

nlm = NonlinearModelFit[data, ax + b, {a, b}, {x,y}]
nlm[{"BestFit", "ParameterTable"}]

How do I find the uncertainty in the slope that is formed by my data? (without having to use tools outside mathematica)

EDIT: for integrity's sake, the source of my data is Table 2 in

[cal68] William A. Calder. “Respiratory and Heart Rates of Birds at Rest”. In The Condor 70.4 (1968), pp. 358–365.doi:http://dx.doi.org/10.2307/1365930

Not of importance to the question however not adding this source would be plagiarism.

• Your original code nlm = NonlinearModelFit[data, ax + b, {a, b}, {x,y}] should be corrected as nlm = NonlinearModelFit[data, a*x + b, {a, b}, x]. Then you get the estimates table. (I did edited your question by then reverted the edit -- please re-edit.) Dec 12, 2020 at 14:29
• I get a NonlinearModelFit::fitc error when I run your code. Just to be clear: Did you get the same errors, too? Dec 12, 2020 at 14:36
• @MichaelE2 Same here -- NonLinearModelFit::fitc. OP mentioned that message (but not the label) in the post. Dec 12, 2020 at 14:40
• @MichaelE2 yes, I also got a NonlinearModelFit::fitc Dec 12, 2020 at 14:43
• @AntonAntonov thank you for that suggestion, I added that edit to my code and everything works fine now. If you were to past that suggestion as an answer, I would gladly accept it:) Dec 12, 2020 at 14:44

I disagree with @AntonAntonov 's assessment of

My hypothesis for this slope was -0.25, and so I would be tempted to just call this a success, however to be truly certain I need to know what the uncertainty is in this slope.

One can directly construct the 95% confidence intervals for the slope:

M = Log /@ {70, 3.4, 21, 21, 600, 352, 17.5, 21.5, 148, 2620, 467000,
1820000, 22, 75, 88, 157, 125, 39.7, 2719000, 3672000, 407000,
57800};
R = Log /@ {50, 800, 168, 186, 33, 58, 140, 120, 55, 33, 27.4, 16,
135, 91, 65, 126, 54, 120, 4.9, 4.5, 16, 26};
data = Transpose@{M, R};
nlm = NonlinearModelFit[data, a x + b, {a, b}, x]
nlm["ParameterConfidenceIntervalTable"]

One sees that -0.25 is included in the 95% confidence interval. That suggests that -0.25 is still a plausible value for the unknown true slope (assuming that the model appropriately describes the data generation process). However, the width of the 95% confidence interval is about 0.1. Whether that confidence interval is tight enough to make the desired decision is a subject matter issue rather than a statistical issue.

Using @AntonAntonov 's quantile regression approach

An alternative is to use quantile regression which fits a not necessarily the same relationship as linear regression. (The choice depends on your objective and NOT on which results you like better.)

In any event, one needs to determine a measure of precision associated with an estimate of the slope from quantile regression. I'm not aware of a built-in estimate of precision (like a standard error) for the quantile regression package but that could very well be just ignorance on my part. However, one can perform a bootstrap to obtain approximate 95% confidence intervals for the underlying slope.

Below is a brute-force (i.e., not efficient and certainly not elegant) approach:

(* Get predicted and residuals from quantile regression *)
result = QRMonUnit[N@data]⟹QRMonQuantileRegressionFit[{1, x},
0.5, Method -> NMinimize];
predicted = (List @@ result[[1]])[[1]] /@ data[[All, 1]];
residuals = data[[All, 2]] - predicted;

(* Perform bootstrap samples *)
n = 1000;
slopes = ConstantArray[0, n];
SeedRandom[123];
Do[boot = Transpose[{data[[All, 1]], predicted + RandomChoice[residuals, Length[data]]}];
out = QRMonUnit[N@boot]⟹QRMonQuantileRegressionFit[{1, x},
0.5, Method -> NMinimize];
(* Grab just the slope from the output *)
slopes[[i]] = (List @@ out[[1]] /. Function -> List)[[1, 1]] /. Plus[a_, Times[b_, Slot[1]]] -> b,
{i, n}]

(* Histogram of results *)
Histogram[slopes, "FreedmanDiaconis", "PDF"]

(* Approximate 95% confidence interval for the underlying slope *)
Quantile[slopes, {0.025, 0.975}]
(* {-0.33217059970198626,-0.214088495992208} *)

One can see that -0.25 falls within the 95% confidence interval. Again, that is not evidence that the slope is -0.25 but rather that -0.25 is still a plausible value to consider also for the quantile regression slope.

You only have 22 data points and there is certainly a fair amount of variability around the fitted lines. So a large confidence interval should not be surprising.

• Thanks for posting this. The thing is that "mean" and "median" disagree so much there might be some problem with the fits and/or data. See the update of my answer. If we remove the second data point, {1.22378, Log[800]}, we get agreeing "mean" and "median", and a slope estimate ~(-0.22). That value is also "too far away" from -0.25. Dec 12, 2020 at 23:27
• The "estimates" do look different between the two methods. But the associated precisions show that sampling error is larger than the difference. Also, comparing estimates of two different methods that have different objectives is a best providing only indirect evidence of differences and/or similarities.
– JimB
Dec 13, 2020 at 5:19
• I like the brute force bootstrap approach! I consider making an add-on to QRMon to do that. Dec 13, 2020 at 15:15
• "(The choice depends on your objective and NOT on which results you like better.)" -- Agreed. Should be in bold not parenthesized. Dec 13, 2020 at 15:25
• LinearModelFit[data, x, x] would work too, producing same results. Feb 15, 2023 at 8:00

### Code correction

NonlinearModelFit[data, ax + b, {a, b}, {x, y}]

should be corrected as

NonlinearModelFit[data, a*x + b, {a, b}, x]

Then you get the estimates table.

### Slope estimate

In the computations below I am using QRMon, but WFR’s QuantileRegression can be also used.

My hypothesis for this slope was -0.25, and so I would be tempted to just call this a success, however to be truly certain I need to know what the uncertainty is in this slope.

You should not believe it -- using Quantile Regression (QR) the slope is $$0.28$$:

qrObj =
QRMonUnit[N@data]⟹
QRMonLeastSquaresFit[{1, x}]⟹
QRMonQuantileRegressionFit[{1, x}, 0.5, Method -> NMinimize]⟹
QRMonPlot[AspectRatio -> Automatic, PlotTheme -> "Detailed"];

Map[#[x] &, qrObj⟹QRMonTakeRegressionFunctions]

(*<|"mean" -> 5.65474 - 0.242763 x, 0.5 -> 5.74435 - 0.280452 x|>*)

### Conditional CDFs

Further, I would suggest to look at QR fits for probabilities other than 0.5.

Also, it is very likely the OP analysis would benefit from examining the conditional CDFs instead of just looking at the mean (median) of the fit.

Here are QR fits additional a list of probabilities:

qrObj2 =
QRMonUnit[N@data]⟹
QRMonEchoDataSummary⟹
QRMonQuantileRegressionFit[{1, x}, {0.01, Sequence @@ Range[0.1, 0.9, 0.1], 0.99}, Method -> NMinimize]⟹
QRMonPlot[AspectRatio -> Automatic, PlotTheme -> "Detailed"];

Out of curiosity let us look the slopes summary:

aSlopes = Coefficient[#, x] & /@ Map[#[x] &, qrObj2⟹QRMonTakeRegressionFunctions]

(*<|0.01 -> -0.22157, 0.1 -> -0.228508, 0.2 -> -0.247381, 0.3 -> -0.272245, 0.4 -> -0.276046, 0.5 -> -0.280452, 0.6 -> -0.217056, 0.7 -> -0.209984, 0.8 -> -0.21576, 0.9 -> -0.193597, 0.99 -> -0.285205|>*)
ResourceFunction["RecordsSummary"][Values@aSlopes]

Here are some of the conditional CDFs:

qrObj2⟹
QRMonConditionalCDF[{5, 12}]⟹
QRMonEchoValue⟹
QRMonConditionalCDFPlot[{5, 8, 12}, ImageSize -> Medium];

### Update: Slopes jackknifed

Here we systematically drop data points, Least squares and Quantile regression, and extract the corresponding slopes:

aQRObjs =
Association@
Map[
# -> QRMonUnit[Drop[data, {#}]]⟹
QRMonLeastSquaresFit[{1, x}]⟹
QRMonQuantileRegressionFit[{1, x}, 0.5, Method -> NMinimize] &,
Range@Length@data];

aPDropSlopes =
Map[Coefficient[#, x] & /@
Map[Simplify[#[x]] &, #⟹QRMonTakeRegressionFunctions] &, aQRObjs];

Here is a grid of plots of the fits:

Multicolumn[
Map[aQRObjs[#]⟹QRMonPlot["Echo" -> False,
AspectRatio -> Automatic,
FrameLabel -> {{None, None}, {aPDropSlopes[#],
Row[{"Drop index: ", #, ", point:", data[[#]]}]}},
PlotTheme -> "Detailed",
ImageSize -> 350]⟹QRMonTakeValue &,
Keys@aQRObjs], 4, Appearance -> "Horizontal", Dividers -> All,
FrameStyle -> GrayLevel[0.8]]

Slopes summaries:

aPDropSlopes2 = Join @@ KeyValueMap[Function[{k, v}, KeyMap[{k, #} &, v]], aPDropSlopes];
Association@
Map[Function[{p},
p -> ResourceFunction["RecordsSummary"][
Values@KeySelect[aPDropSlopes2, MatchQ[#, {_, p}] &]]], {"mean", 0.5}]

We can see that removing the second data point, {1.22378,Log[800]}, produces Least squares and Quantile regression fits that agree; the slope is approximately (-0.22).

Now, I think -0.22 is too far away from the hypothesis slope -0.25, and the slope from NonLinearModelFit shown in JimB's answer. (Approximately, 1 and 1.5 standard errors away, respectively.)

• To decide if -0.25 is a plausible value for the true value of a parameter one needs a measure of precision to judge "distance" from an estimate. I see no measures of precision described here. Finding that a point that could be removed in which the least squares slope is very similar to the quantile regression estimate for medians is not a standard justification for making decisions. What is the measure of precision for the quantile regression slope? I assume one would either use the jackknife estimates or a bootstrap procedure to obtain such a measure.
– JimB
Dec 13, 2020 at 4:08
• The jackknife values for "0.5" are (essentially) -0.215 (6 values) or -0.28 (16 values) only. That's a lot of sampling variability that far exceeds the difference between -0.25 and -0.28.
– JimB
Dec 13, 2020 at 4:14
• @JimB This is a good point — I noticed the “binarization” of values, but forgot to reason with it. I think, though, that using the conditional CDFs is a better way of analyzing the data. Getting the slope is more of a "summarization courtesy" to the consumers of the analysis. Dec 13, 2020 at 15:24