[R] Straight

Prof Brian D Ripley ripley at stats.ox.ac.uk
Wed Nov 1 19:09:33 CET 2000


On Wed, 1 Nov 2000, Zsombor Cseres-Gergely wrote:

> On Wed, Nov 01, 2000 at 08:23:36AM +0000, Prof Brian Ripley wrote:
> 
> > What do you want there: estimation under restrictions or tests of linear
> > combinations?  The latter is a simple R programming exercise. I don't know
> 
> Sorry for bad wording, I meant the latter. And this is why I did not ask this
> before: easy, _provided_ we have the covariance matrix of the estimators. Now
> I have the answer also for this, so it is really straight. (I still do not
> know, how summarize() gets it, since in a lm.model, I have not found the
> underlying data for it [or overlooked]. The reason I did not do a X'X was that
> I believe that this is expensive. I have not tried, but my previous
> experiences (again: Stata) tells me this. Or am I wrong, and this is cheap in R?

Modern linear models theory does not use X'X ^-1 but R'R ^-1 for the QR
decomposition, and R is there in the object. Further, all you have to do is
invert R and crosspoduct the inverse. That's cheap and stable.  Here's the
exact code from summary.lm:

    chol2inv(Qr$qr[p1, p1, drop = FALSE])

which is X'X^-1.  chol2inv is a Fortran routine to invert an matrix from
its Cholesky decomposition, and R (or moe precisely the top p1 rows of R)
is the (possibly pivoted) Cholesky decomposition of X'X.  On the other
hand, S just inverts R by solve().


-- 
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 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list