2
$\begingroup$

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.

$\endgroup$
8
  • 2
    $\begingroup$ 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.) $\endgroup$ Dec 12, 2020 at 14:29
  • $\begingroup$ I get a NonlinearModelFit::fitc error when I run your code. Just to be clear: Did you get the same errors, too? $\endgroup$
    – Michael E2
    Dec 12, 2020 at 14:36
  • $\begingroup$ @MichaelE2 Same here -- NonLinearModelFit::fitc. OP mentioned that message (but not the label) in the post. $\endgroup$ Dec 12, 2020 at 14:40
  • $\begingroup$ @MichaelE2 yes, I also got a NonlinearModelFit::fitc $\endgroup$
    – Poseidaan
    Dec 12, 2020 at 14:43
  • $\begingroup$ @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:) $\endgroup$
    – Poseidaan
    Dec 12, 2020 at 14:44

2 Answers 2

4
$\begingroup$

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"]

Parameter confidence interval table

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"]

Histogram of bootstrap slope estimates

(* 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.

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ Dec 12, 2020 at 23:27
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Dec 13, 2020 at 5:19
  • $\begingroup$ I like the brute force bootstrap approach! I consider making an add-on to QRMon to do that. $\endgroup$ Dec 13, 2020 at 15:15
  • 1
    $\begingroup$ "(The choice depends on your objective and NOT on which results you like better.)" -- Agreed. Should be in bold not parenthesized. $\endgroup$ Dec 13, 2020 at 15:25
  • $\begingroup$ LinearModelFit[data, x, x] would work too, producing same results. $\endgroup$
    – kirma
    Feb 15, 2023 at 8:00
3
$\begingroup$

Code correction

Your original code

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"];

enter image description here

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"];

enter image description here

enter image description here

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]

enter image description here

Here are some of the conditional CDFs:

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

enter image description here

enter image description here

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]]

enter image description here

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}]

enter image description here

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).

enter image description here

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.)

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Dec 13, 2020 at 4:08
  • 1
    $\begingroup$ 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. $\endgroup$
    – JimB
    Dec 13, 2020 at 4:14
  • $\begingroup$ @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. $\endgroup$ Dec 13, 2020 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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