[R] Change in 'solve' for r-patched

[Prof Brian Ripley]

GLS is usually solved by taking a matrix square root and converting to 
least squares (and BTW you might want to use lm.fit not solve for least 
squares to allow for aliased columns).

So your idea is right (and is the one given in my 1981 Spatial Statistics 
book).  These days I would usually use an eigendecomposition instead of 
Cholesky as it will enable you to cope better with nearly 
non-positive-definite Omega.  See lm.gls in package MASS for an outline 
implementation.

On Fri, 31 Oct 2003, Ole F. Christensen wrote:

> Dear R help
> 
> Thanks to Professor Bates for his information about how to calculate 
> least square estimates [not using solve] in ``the right way''.
> This is very useful indead, I am clearly one of of the package 
> maintainers who is not using using solve in a proper way at the moment.
> However, the calculations in my code look more like GLS than LS.
> 
> ## GLS could in principlpe be implemented like this :
> betahat <- solve(t(X) %*% solve(Omega)%*% X) %*% t(X)%*%solve(Omega)%*% y
> ## where Omega is a strictly p.d. symmetric matrix
> 
> Does someone have a recommendation on how to do this in ``the right way'' ?
> 
> My first attempt (trying to imitate the LS solution recommended by Prof. Bates) is :
> 
> temp <- backsolve(chol(Omega),cbind(X,y))
> betahat <- qr.coef(qr(temp[,1:ncol(X)]), temp[,ncol(X)+1])
> 
> 
> 
> Thank you in advance for any help
> 
> 
> Cheers Ole
> 
> 

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595