3
$\begingroup$

I have a set of measurements y taken at 17 different values of x, with 50 repeated measurements at each value of x. They follow a simple linear relationship y = mx + c, and I am fitting the parameters m (slope) and c (intercept) using fitting functions in MatLab (though I believe this is more of a statistics than programming issue).

I have attempted the fit in two ways:

1) Sending every single (17*50 = 850) measurement to the fitting function.

2) First calculating the mean measurement at each value of x (17 mean values), and the variance of the measurements at each value of x. I then send the mean values to the fitting function, with 1/variance of each measurement used as the weight of each point (as stated to be optimal in MatLab documentation and reading around linear regression that I've done).

Both methods give almost identical values for m and c, but the calculated error on the parameters is very different. For example, on the fitted slope of magnitude 416 I get a standard error of 0.2 with method 1 (sending each individual measurement to the fitting function) but an error of 1.8 with method 2 (using the mean measurements with weights).

Naively, I expected both methods to give similar errors. Could anyone give an insight about whether this difference is expected, and which method is more robust/which errors are more realistic if it is?

For reference I have tried using (happy to provide code snippets/more information if the above isn't enough):

MatLab "fit" function (via cftool), which gives 95% confidence intervals that I am dividing by 2 to give the above standard errors.

MatLab LinearModel.fit function, which gives very similar results to the above (although it gives a standard error directly).

EDIT: A code snippet in response to whuber's reply below.

% create rough linear distribution with some random normally distributed errors
nPoints = 11;
nRep = 50;
x = NaN(nRep,nPoints);
y = NaN(nRep,nPoints);
for i=1:nRep
    x(i,:)= 1:nPoints;
    y(i,:) = (1:nPoints) + random('norm',0,3,1,nPoints);
end

% calculate mean and variance at each value of x
yMean = mean(y,1);
yVar = var(y,0,1);

fprintf('---------------------------------------------------------------\n');
fprintf('fit using 1/variance as weights\n');
fprintf('---------------------------------------------------------------\n');
LinearModel.fit(1:nPoints, yMean, 'Weights', 1./yVar)

fprintf('\n---------------------------------------------------------------\n');
fprintf('fit using sqrt(nRep-1)/variance as weights\n');
fprintf('---------------------------------------------------------------\n');
LinearModel.fit(1:nPoints, yMean, 'Weights', sqrt(nRep-1)./yVar)

fprintf('\n---------------------------------------------------------------\n');
fprintf('fit using all points rather than mean\n');
fprintf('---------------------------------------------------------------\n');
LinearModel.fit(x(:), y(:))

This code gives the following fitter parameters:

mean values with 1/variance as weights:

intercept: estimate = 0.029077, error = 0.27823

slope: estimate = 1.0291, error = 0.040575

mean values with sqrt(nRep-1)/variance as weights:

intercept: estimate = 0.029077, error = 0.27823

slope: estimate = 1.0291, error = 0.040574

all measurements:

intercept: estimate = -0.076889, error = 0.27133

slope: estimate = 1.0427, error = 0.040006

It does not recreate my problem with large differences in the error using all measurements or using the mean with weights. However, in the two examples using the mean value the different absolute values of the weights has no effect on the standard error of the parameters.

$\endgroup$
  • $\begingroup$ Using the variances as relative weights was the right idea. However, they are not the variances of the means: those need to be divided by $50-1$ if you want meaningful standard errors to be estimated. The resulting standard errors in the fit will therefore be approximately $\sqrt{50-1}=7$ times too great. (I presume that $7\times 0.2\approx 1.8$ due to imprecision in reporting these values.) Could you confirm that this indeed is the case? $\endgroup$ – whuber Mar 22 '16 at 14:08
  • $\begingroup$ Thanks for the reply! That makes complete sense but in fact I already tried that in the past and it makes no difference to the standard error on the parameters returned by the MatLab functions. I just did it again to make sure and there is no difference - as long as the relative weights are maintained the absolute value seems to change nothing. $\endgroup$ – Jack Roberts Mar 22 '16 at 16:56
  • $\begingroup$ There are various ways in which software can accommodate and interpret weights in regression. (Stata, for instance, offers four distinct ones.) Possibly, Matlab might not be using these weights as you intended. It would help if you could provide a small reproducible example, because there are many conceivable explanations for the differences in output. $\endgroup$ – whuber Mar 22 '16 at 17:06
  • 1
    $\begingroup$ I added a MatLab code snippet to the original post above. It seems to give similar errors fitting using both the mean measurements and all measurements, unlike my real data set, so I will see if I can understand why that is the case. However, it recreates the issue with the absolute value of the weights having no effect on the returned error of the fit parameters. $\endgroup$ – Jack Roberts Mar 22 '16 at 17:47
  • $\begingroup$ Thank you. The majority of readers of this question will not have access to MatLab, so it would help also to include the key parts of its output. That would open your question up to consideration by a greater number of experts. $\endgroup$ – whuber Mar 22 '16 at 17:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.