From ripley at stats.ox.ac.uk  Mon Jul 31 08:52:01 2006
From: ripley at stats.ox.ac.uk (Prof Brian Ripley)
Date: Mon, 31 Jul 2006 07:52:01 +0100 (BST)
Subject: [Rd] NIST StRD linear regression
In-Reply-To: <0F923E9EBBF1554ABA641ACA04BEBEC0051BD9FE@WS-BCO-MSE7.milky-way.battelle.org>
References: <0F923E9EBBF1554ABA641ACA04BEBEC0051BD9FE@WS-BCO-MSE7.milky-way.battelle.org>
Message-ID: <Pine.LNX.4.64.0607310735001.20216@gannet.stats.ox.ac.uk>

That is not an appropriate way to fit a degree-10 polynomial (in any 
language, if fitting a degree-10 polynomial is in fact an appropriate 
statistical analysis, which seems unlikely).

On Sun, 30 Jul 2006, Carnell, Rob C wrote:

> NIST maintains a repository of Statistical Reference Datasets at
> http://www.itl.nist.gov/div898/strd/.  I have been working through the
> datasets to compare R's results to their references with the hope that
> if all works well, this could become a validation package.

What does it validate?  The R user's understanding of numerical methods?

> All the linear regression datasets give results with some degree of
> accuracy except one.  The NIST model includes 11 parameters, but R will
> not compute the estimates for all 11 parameters because it finds the
> data matrix to be singular.
> 
> The code I used is below.  Any help in getting R to estimate all 11
> regression parameters would be greatly appreciated.
> 
> I am posting this to the R-devel list since I think that the discussion
> might involve the limitations of platform precision.
> 
> I am using R 2.3.1 for Windows XP.
> 
> rm(list=ls())
> require(gsubfn)

That is not needed.

> defaultPath <- "my path"
> 
> data.base <- "http://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA"
> 
> reg.data <- paste(data.base, "/Filip.dat", sep="")
> 
> model <-
> "V1~V2+I(V2^2)+I(V2^3)+I(V2^4)+I(V2^5)+I(V2^6)+I(V2^7)+I(V2^8)+I(V2^9)+I
> (V2^10)"
> 
> filePath <- paste(defaultPath, "//NISTtest.dat", sep="")
> download.file(reg.data, filePath, quiet=TRUE)

filePath <- 
url("http://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/Filip.dat")

will suffice.

> A <- read.table(filePath, skip=60, strip.white=TRUE)

> lm.data <- lm(formula(model), A)
> 
> lm.data

lm(V1 ~ poly(V2, 10), A)

works.

> kappa(model.matrix(V1 ~ poly(V2, 10, raw=TRUE), A), exact=TRUE)
[1] 1.767963e+15

shows the design matrix is indeed numerically singular by the naive 
method.

-- 
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



[Rd] NIST StRD linear regression
Prof Brian Ripley