Simple lm/regression question

I am trying to use lm for a simple linear fit with weights. The results 
I get from IDL (which I am more familiar with) seem correct and 
intuitive, but the "lm" function in R gives outputs that seem strange to 
me.

Unweighted case:

x<-1:4
y<-(1:4)^2
summary(lm(y~x))
Call:
lm(formula = y ~ x)

Residuals:
1  2  3  4
1 -1 -1  1

Coefficients:
           Estimate Std. Error t value Pr(>|t|)
(Intercept)  -5.0000     1.7321  -2.887   0.1020
x             5.0000     0.6325   7.906   0.0156 *

So far, so good - IDL does much the same:

IDL> vec=linfit(x,y,sigma=sig)
IDL> print,vec,sig
    -5.00000      5.00000
     1.73205     0.632456

Now, if the dependent variable has known (measurement) uncertainties (10, 
say) then it is appropriate to use weights defined as the inverse of the 
variances, right?
If you think that in y = x'b + e the error has

E(e^2) = sigma^2 * w

then you should use

weights = 1/w

in R because that is how the weights enter the objective function

sum 1/w (y - x'b)^2
summary(lm(y~x,weights=c(.01,.01,.01,.01)))
Call:
lm(formula = y ~ x, weights = c(0.01, 0.01, 0.01, 0.01))

Residuals:
  1    2    3    4
0.1 -0.1 -0.1  0.1

Coefficients:
           Estimate Std. Error t value Pr(>|t|)
(Intercept)  -5.0000     1.7321  -2.887   0.1020
x             5.0000     0.6325   7.906   0.0156 *

But here the *residuals* have changed and are no longer the actual 
response minus fitted values. (They seem to be scaled by the sqrt of the 
weights).
The summary() shows under "Residuals" the contributions to the objective 
function, i.e. sqrt(1/w) (y - x'b) in the notation above.

However, by using the residuals extractor function you can get the 
unweighted residuals:

residuals(lm(y~x,weights=c(.01,.01,.01,.01)))
The uncertainties on the parameter estimates, however, have 
*not* changed, which seems very odd to me.
lm() interprets the weights as precision weights, not as case weights.

Thus, the scaling in the variances is done by the number of (non-zero) 
weights, not by the sum of weights.
The behaviour of IDL is rather different and intuitive to me:

IDL> vec=linfit(x,y,sigma=sig,measure_errors=[1,1,1,1])
IDL> print,vec,sig
    -5.00000      5.00000
     1.22474     0.447214

IDL> vec=linfit(x,y,sigma=sig,measure_errors=[10,10,10,10])
IDL> print,vec,sig
    -5.00000      5.00000
     12.2474      4.47214
This appears to use sandwich standard errors. If you load the sandwich and 
lmtest package you can do

coeftest(m, vcov = sandwich)

where 'm' is the fitted lm object.

Hope that helps,
Z
Note that the uncertainties are 10* larger when the errors are 10* larger.

My question is really how to replicate these results in R. I suspect there is 
some terminology or convention that I'm confused by. But in the lm help page 
documentation, even the straightforward definition of "residual" seems 
incompatible with what the code actually returns:

	residuals: the residuals, that is response minus fitted values.

But, in the weighted case, this is not actually true...

James
-- 
James D Annan jdannan at jamstec.go.jp Tel: +81-45-778-5618 (Fax 5707)
Senior Scientist, Research Institute for Global Change, JAMSTEC
Yokohama Institute for Earth Sciences, 3173-25 Showamachi,
Kanazawa-ku, Yokohama City, Kanagawa, 236-0001 Japan
http://www.jamstec.go.jp/frcgc/research/d5/jdannan/

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Simple lm/regression question

Thread (4 messages)