2
$\begingroup$

This may be a question which teeters on the brink of belonging more to the realm of statistics and Cross-Validated SE, but I'm also specifically interested in Mathematica fitting routines.

Usually, if I want to fit a model to some data using either NonlinearModelFit or LinearModelFit I will have some errors associated with my $y$-data that I use to weight the fits. These errors might be simply the standard error acquired from repeated measurements, or I might know something about the physical processes and can assign weights.

For example Weights->1/YDataErrors^2 and I always set my variance estimator as VarianceEstimatorFunction -> (1 &). I can then get my parameter errors from the covariance matrix, or simply with MyFit["ParameterErrors"].

However in some cases one might not have any errors for the data you want to fit, meaning one cannot provide weights in the way that I described above. My question is then, how reliable -- or more importantly -- how physically/statistically meaningful are the parameter errors for an unweighted fit in Mathematica?

$\endgroup$
2
  • $\begingroup$ You continue to challenge me to think in ways that I would never had considered. (At least I think that's a good thing.) Most folks (outside of physics it seems) don't use weighted data. But when they do, they also specify an error structure appropriate to the weighting. The holy triplet is the model ("the equation"), the data, and the error structure. Once you specify all three of those, then the approach to fitting and whether or not to weight becomes (relatively) clear. $\endgroup$ – JimB Aug 16 '20 at 16:42
  • $\begingroup$ When there are measurement errors found from repeated measurements, one also always needs to consider the "lack-of-fit" error to the model. (As John Tukey states: "All models are wrong but some are useful.") Currently (and it would be great if this would change) Mathematica only allows for a single type of error in the fitting process (unless you write your own functions). Mathematica needs to add "mixed models" where there are multiple error sources. (These can be very difficult and numerically unstable to fit which is likely why folks at Wolfram, Inc. haven't added this yet.) $\endgroup$ – JimB Aug 16 '20 at 16:50
4
$\begingroup$

For example, if one has two sources of error, say a measurement error and a lack-of-fit error, then using the weights based on the measurement errors can result in gross underestimates of the standard errors. Consider the following model:

$$y=a+b x +\gamma + \epsilon$$

where $y$ is the measured response, $x$ is the predictor, $a$ and $b$ are constants to be estimated, $\gamma$ is the repeated measurement error with $\gamma \sim N(0,\sigma_{ME})$, and $\epsilon$ is the lack-of-fit error with $\epsilon \sim N(0,\sigma)$ and all errors are assumed to be independent.

First set some specific parameters:

(* Measurement error standard deviation *)
σME = 10;

(* Lack-of-fit error standard deviation *)
σ = 20;

(* Regression coefficients *)
a = 1;
b = 1;

Generate and plot some data:

n = 100;
x = Range[n];
SeedRandom[12345];
measurementError = RandomVariate[NormalDistribution[0, σME], n];
lackOfFitError = RandomVariate[NormalDistribution[0, σ], n];
y = a + b x + measurementError + lackOfFitError;
data = Transpose[{x, y}];
data2 = {#[[1]], Around[#[[2]], σME]} & /@ data;
ListPlot[data2]

Data with associated measurement error bars

Now consider two different linear model fits where lm1 is what you suggest and lm2 is what I suggest:

lm1 = LinearModelFit[data, z, z, Weights -> 1/ConstantArray[σME^2, n],
   VarianceEstimatorFunction -> (1 &)];
lm2 = LinearModelFit[data, z, z];
lm1["ParameterTable"]

Parameter table for lm1

lm2["ParameterTable"]

Parameter table for lm2

The estimates of the parameters are identical but the standard errors for lm1 are less than half the size as those for lm2. Which one is correct?

The "true" covariance matrix of the least squares estimators of a and b for this model is

$$\left(\sigma ^2+\sigma_{ME}^2\right) \left(X^T.X\right)^{-1}$$

where $X$ is the design matrix. In Mathematica code the standard error for b is

X = Transpose[{ConstantArray[1, n], Range[n]}]
Sqrt[(σME^2 + σ^2) Inverse[Transpose[X].X][[2, 2]]] // N
(* 0.0774635 *)

That matches pretty well with lm2.

This is a slightly contrived example in that I have all of the measurement standard errors identical because Mathematica's regression functions only allow a single error term. And by having the measurement standard errors identical, that results in an equivalent model with a single error.

However, even when the measurement standard deviations vary considerably, the issue about weighting improperly such that it doesn't match the error structure of the model remains.

Mathematica's regression routines are not yet adequate for models with more than one source of error. I wish they were.

$\endgroup$
2
  • $\begingroup$ Okay this is really interesting, I was always working under the assumption that errors from weightless fits were more-or-less meaningless. I need to think about this a bit more. Thanks for the detailed answers as always. $\endgroup$ – Q.P. Aug 16 '20 at 20:16
  • 1
    $\begingroup$ So the error variance being estimated by lm2 is $\sigma^2+\sigma^2_{ME}$ and not the individual variances. If the measurement standard errors are really known, then the LogLikelihood function will help with that. But what I mean as "known" is not the square root of the measurement as (at best) that just gives one an estimate of the standard deviation. By "known" I mean "really known" and I find that it's rare that anything is "really known." ("Believed to be true", yes. That's why there are so many conspiracy theories - even in science. But "really known", no.) $\endgroup$ – JimB Aug 16 '20 at 20:24

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.