lm weights and the standard error

Question

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.

Answer 1

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.

Answer 2

It seems that the u values represent known standard errors of the y values. However, as pointed out by Glen_b, with lm(y ~ x, data = d.f, weights = 1/u^2), the variances (u^2) are treated as if they are known only up to a proportionality constant. If you want to fit a model with known standard errors, you can approach this from a meta-analytic perspective (where we are in the same situation where the standard errors of the effect size estimates are known). Hence:

library(metafor)
rma(y ~ x, u^2, method="FE", data=d.f, digits=5)

will yield an estimate of the slope (not intercept!) equal to $1.9$ with standard error equal to $0.04472$, which is what the OP was expecting.

Answer 3

The reported standard error for fit2 is not correct, indeed. The reason for this is that the residual standard error for the inverse-variance weighted regression should be 1. Hence, to obtain the correct value of the uncertainty in the slope for fit2 you have to divide the reported standard error by the reported residual standard error.

       fit2 <- lm(y ~ x, data = d.f, weights = 1/u^2)
       SE <- summary(fit2)$coef[2,2]/summary(fit2)$sigma

SE = 0.04472136 - the correct SE of the slope.