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.