naive "collinear" weighted linear regression
David Winsemius <dwinsemius <at> comcast.net> writes:
Which means those x, y, and "error" figures did not come from an experiment, but rather from theory???
The fact is I am trying to compare the results of: (1) lm under R and (2) the Java applet at http://omnis.if.ufrj.br/~carlos/applets/reta/reta.html (3) the Fit method of the ROOT system used by CERN, (4) the Numerical Recipes functions for weighted linear regression The three last ones all provide, for the "fake" data set I furnished in my first post in this thread, the same results; particularly they give erros or uncertainties in the estimated coefficients of intercept and slope which, as seems intuitive, are not zero at all, but of the order 0.1 or 0.2, whereas the method lm under R issues a "Std. Error", which is zero. Independently of terminology, which sure is of utmost importance, the data I provided should give rise to a best fit straight line with intercept zero and slope 2, but also with non-vanishing errors associated with them. How do I get this from lm???? I only want, for instance, calculation of the so-called covariance matrix for the estimated coefficients, as given, for instance, in Equation (2.3.2) of the second edition of Draper and Smith, "Applied regression analysis"; this is a standard statistical result, right? So why does R not directly provide it as a summary from an lm object???
Of course the best fit coefficients should be 0 for the intercept and 2 for the slope. Furthermore, it seems completely plausible (or not?) that, since the y_i have associated non-vanishing ``errors'' (dispersions), there should be corresponding non- vanishing ``errors'' associated to the best fit coefficients, right? When I try:
fit_mod <- lm(y~x,weights=1/error^2)
I get Warning message: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : extra arguments weigths are just disregarded.
(Actually the weights are for adjusting for sampling, and I do not see any sampling in your "design".)
Keeping on, despite the warning message, which I did not quite understand, when I type:
summary(fit_mod)
I get
Call:
lm(formula = y ~ x, weigths = 1/error^2)
Residuals:
1 2 3 4
-5.067e-17 8.445e-17 -1.689e-17 -1.689e-17
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.000e+00 8.776e-17 0.000e+00 1
x 2.000e+00 3.205e-17 6.241e+16 <2e-16 ***
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Residual standard error: 7.166e-17 on 2 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 3.895e+33 on 1 and 2 DF, p-value: < 2.2e-16
Naively, should not the column Std. Error be different from zero??
What I have in mind, and sure is not what Std. Error means, is that
if I carried out a large simulation, assuming each response y_i a
Gaussian random variable with mean y_i and standard deviation
2*error=0.6, and then making an ordinary least squares fitting of
the slope and intercept, I would end up with a mean for these
simulated coefficients which should be 2 and 0, respectively,
Well, not precisely 2 and 0, but rather something very close ... i.e, within "experimental error". Please note that numbers in the range of 10e-17 are effectively zero from a numerical analysis perspective.
http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f
> .Machine$double.eps ^ 0.5
[1] 1.490116e-08
I know this all too well and it is obviously a trivial supernewbie issue, which I have already overcome a long time ago...
and, that's the point, a non-vanishing standard deviation for these fitted coefficients, right?? This somehow is what I expected should be an estimate or, at least, a good indicator, of the degree of uncertainty which I should assign to the fitted coefficients; it seems to me these deviations, thus calculated as a result of the simulation, will certainly not be zero (or 3e-17, for that matter). So this Std. Error does not provide what I, naively, think should be given as a measure of the uncertainties or errors in the fitted coefficients...
You are trying to impose an error structure on a data situation that you constructed artificially to be perfect.
What am I not getting right??
That if you input "perfection" into R's linear regression program, you get appropriate warnings?
Thanks and sorry for the naive and non-expert question!
You are a Professor of physics, right? You do experiments, right? You replicate them. S0 perhaps I'm the one who should be puzzled.
Unfortunately you eschewed answering objectively any of my questions; I insist they do make sense. Don't mention the data are perfect; this does not help to make any progress in understanding the choice of convenient summary info the lm method provides, as compared to what, in my humble opinion and in this specific particular case, it should provide: the covariance matrix of the estimated coefficients...
-- ####################################### Prof. Mauricio Ortiz Calvao Federal University of Rio de Janeiro Institute of Physics, P O Box 68528 CEP 21941-972 Rio de Janeiro, RJ Brazil