How meaningful, or useful, are parameter errors produced when perfroming an unweighted LinearModelFit or NonlinearModelFit?

Question

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?

Answer 1

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]

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

lm2["ParameterTable"]

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.