[R] naive "collinear" weighted linear regression

[Mauricio O Calvao]

David Winsemius <dwinsemius <at> comcast.net> writes:


> 
> It's really not that difficult to get the variance covariance matrix.  
> What is not so clear is why you think differential weighting of a set  
> that has a perfect fit should give meaningfully different results than  
> a fit that has no weights.

Again, David, what I have in mind is: since there are errors or uncertainties in
the response variables (despite the perfect collinearity of the data), which I
assume are Gaussian, if I make a large enough number of simulations of four
response values, there will undoubtedly be a dispersion in the best fit
intercept and slope obtained from a usual unweighted least squares procedure,
right? Then, if I calculate the arithmetic mean of these simulated intercept and
slope, I would certainly check that they would be 0 and 2, respectively.
However, and THAT IS THE POINT, there will also be a standard deviation
associated with each one of these two coefficients, right??, and that is what I
would assign as the measure of uncertainty in the estimation of the
coefficients. This is not, as Dalgaard has called attention to, what the simple
command summary(lm(y~x,weights=1/err^2)) provides in its Std. Error. However, as
Dalgaard also recalled, the command
summary(glm(y~x,family=gaussian,weights=1/err^2),dispersion=1) does provide Std.
Errors in the coefficients which look plausible (at least to me) and, at any
rate, which do coincide with results from other packages (Numerical Recipes,
ROOT and possibly GSL...)

> 
> ?lm
> ?vcov
> 
>  >  y <- c(2,4,6,8) # response vect
>  > fit_mod <- lm(y~x,weights=1/error^2)
> Error in eval(expr, envir, enclos) : object 'error' not found
>  > error <- c(0.3,0.3,0.3,0.3)
>  > fit_mod <- lm(y~x,weights=1/error^2)
>  > vcov(fit_mod)
>                (Intercept)             x
> (Intercept)  2.396165e-30 -7.987217e-31
> x           -7.987217e-31  3.194887e-31
> 
> Numerically those are effectively zero.
> 
>  > fit_mod <- lm(y~x)
>  > vcov(fit_mod)
>              (Intercept) x
> (Intercept)           0 0
> x                     0 0
>