Solved – lm weights and the standard error

regression

I am fairly new to R, but am experienced in other languages and also in data analysis in the physical sciences. I have a problem and will illustrate with a straight line fit using lm()

In the physical sciences, if the dependent variable (y) has explicit uncertainties (also called errors) and those uncertainties all have the same value, then if we fit to a straight line then the values of the slope and intercept are unchanged but the uncertainties in the slope and intercept will be different than if the data do not have explicit uncertainties. I can not find an R fitting routine that does the correct thing. I've written fitters that do the right thing in C, Mathematica, Maple, and LabView, but hope to avoid doing so with R.

Here's an example:

d.f = data.frame( x = c(1, 2, 3, 4), y = c(2, 3.9, 6.1, 7.9), u = c(.1, .1, .1, .1))

The variable u is the uncertainties in the values of y. Do a straight line fit ignoring the uncertainties:

fit1 <- lm(y ~ x, data = d.f)

This gives exactly what I think is the correct fit, including the "Std. Error" in the slope and intercept.

The weights are conventionally $1 / (\text{uncertainty}^2)$, so I give that to lm()

fit2 <- lm(y ~ x, data = d.f, weights = 1/u^2)

The Residual standard error is different for these 2 fits, as it should be, but the Std. Error in the fitted parameters is not.

For example, the intercept for both fits is $1.99$, which is correct, and the Std. error in the intercept for fit1 is $0.05196$, which is also correct. However, these values are exactly the same for fit2. The correct value for the uncertainty in the intercept for fit2 is actually $0.04472$.

I've searched around a fair amount, and have failed to find a least-square regression fitter in R that does this correctly. In fact, the ideal such fitter would also accept uncertainties in both the x and y coordinates, and fit using an "effective variance" technique. Or perhaps I just don't understand how to use lm() or interpret its results.

Best Answer

I think you are looking for:

fit2 = lm(y ~ x, data = d.f, weights = 1/u^2)
parameter.covariance.matrix = vcov(fit2)/summary(fit2)$sigma^2

there is probably a better way to get this information from the fit2 object. This works even if it isn't elegant.

Related Question