5
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Please explain what an "uncertainty" quantifies. Is it a standard error? A confidence interval? A variance? And would these be assessing measurement error or regression error? (The phrasing in this post seems to equate the two, but whether that is appropriate is doubtful.) $\endgroup$ – whuber Sep 1 '14 at 19:11
  • $\begingroup$ Weighted regression only gives relative weights. Are your known uncertainties actually known standard deviations of a noise term (which would only require a fairly simple adjustment of the output), or are they something else? $\endgroup$ – Glen_b Sep 1 '14 at 22:34
  • $\begingroup$ Could you characterize (in as near to explain like I'm five level as feasible) what the uncertainties represent, and how they vary from observation to observation? Do they include, for example, any degree of instruments consistently slightly low or high? Are they strict upper bounds on measurement error? Some kind of weak bound on measurement error? And actual standard deviation (and if so, how is that obtained)? $\endgroup$ – Glen_b Sep 1 '14 at 22:56
  • $\begingroup$ Regarding the errors in the independent variables, I'm not aware of any R packages that solve that problem. The best solution I'm aware of is ODRpack which is available from netlib, but that's in Fortran. $\endgroup$ – shane Aug 7 '15 at 19:41
3
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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.

$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.