Re: [Rd] NIST StRD linear regression
> "RobCar" == Carnell, Rob C <[EMAIL PROTECTED]> > on Sun, 30 Jul 2006 19:42:29 -0400 writes: RobCar> NIST maintains a repository of Statistical Reference RobCar> Datasets at http://www.itl.nist.gov/div898/strd/. I RobCar> have been working through the datasets to compare RobCar> R's results to their references with the hope that RobCar> if all works well, this could become a validation RobCar> package. RobCar> All the linear regression datasets give results with RobCar> some degree of accuracy except one. The NIST model RobCar> includes 11 parameters, but R will not compute the RobCar> estimates for all 11 parameters because it finds the RobCar> data matrix to be singular. RobCar> The code I used is below. Any help in getting R to RobCar> estimate all 11 regression parameters would be RobCar> greatly appreciated. RobCar> I am posting this to the R-devel list since I think RobCar> that the discussion might involve the limitations of RobCar> platform precision. RobCar> I am using R 2.3.1 for Windows XP. RobCar> rm(list=ls()) RobCar> require(gsubfn) RobCar> defaultPath <- "my path" RobCar> data.base <- "http://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA"; Here is a slight improvement {note the function file.path(); and model <- ..; also poly(V2, 10) !} which shows you how to tell lm() to "believe" in 10 digit precision of input data. --- reg.data <- paste(data.base, "/Filip.dat", sep="") filePath <- file.path(defaultPath, "NISTtest.dat") download.file(reg.data, filePath, quiet=TRUE) A <- read.table(filePath, skip=60, strip.white=TRUE) ## If you really need high-order polynomial regression in S and R, ## *DO* as you are told in all good books, and use orthogonal polynomials: (lm.ok <- lm(V1 ~ poly(V2,10), data = A)) ## and there is no problem summary(lm.ok) ## But if you insist on doing nonsense 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)" ## MM: "better": (model <- paste("V1 ~ V2", paste("+ I(V2^", 2:10, ")", sep='', collapse=''))) (form <- formula(model)) mod.mat <- model.matrix(form, data = A) dim(mod.mat) ## 82 11 (m.qr <- qr(mod.mat ))$rank # -> 10 (only, instead of 11) (m.qr <- qr(mod.mat, tol = 1e-10))$rank # -> 11 (lm.def <- lm(form, data = A)) ## last coef. is NA (lm.plus <- lm(form, data = A, tol = 1e-10))## no NA coefficients --- RobCar> reg.data <- paste(data.base, "/Filip.dat", sep="") RobCar> model <- RobCar> "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 RobCar> (V2^10)" RobCar> filePath <- paste(defaultPath, "//NISTtest.dat", sep="") RobCar> download.file(reg.data, filePath, quiet=TRUE) RobCar> A <- read.table(filePath, skip=60, strip.white=TRUE) RobCar> lm.data <- lm(formula(model), A) RobCar> lm.data RobCar> Rob Carnell A propos NIST StRD: If you go further to NONlinear regression, and use nls(), you will see that high quality statistics packages such as R do *NOT* always conform to NIST -- at least not to what NIST did about 5 years ago when I last looked. There are many nonlinear least squares problems where the correct result is *NO CONVERGENCE* (because of over-parametrization, ill-posednes, ...), owever many (cr.p) pieces of software do "converge"---falsely. I think you find more on this topic in the monograph of Bates and Watts (1988), but in any case, just install and use the CRAN R package 'NISTnls' by Doug Bates which contains the data sets with documentation and example calls. Martin Maechler, ETH Zurich __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] NIST StRD linear regression
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, [EMAIL PROTECTED] 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, UKFax: +44 1865 272595 __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] NIST StRD linear regression
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. 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) 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) A <- read.table(filePath, skip=60, strip.white=TRUE) lm.data <- lm(formula(model), A) lm.data Rob Carnell __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel