[R] software comparison

[Ben Bolker]

Philippe Grosjean <phgrosjean <at> sciviews.org> writes:

> 
> Hello,
> 
> For a paper published in 2007, and submitted in April 2005, this is 
> still surprising. If I my calculation is correct, in 2004, they would 
> have used R 2.2.x, or something,... not 1.9.1?
> 
> Anyway, does someone know if there is a chance 2.5.0 provides some 
> improvements in some of the "difficult cases" for R 1.9.1 (for instance, 
> improvement of the algorithms for calculating autocorrelation, ANOVA, 
> linear regression or non linear regression)?
> 

  I did a bit of testing (with 2.6.0 unstable).
Two of the tests that R 1.9.1 did worst on were
ANOVAs of NIST data sets SiRstv and AgWt 
(LREs listed in the paper for these tests are
 5.7 and 6.7 respectively, compared
with SAS 9.1's 12.7 and 8.8, Splus 6.2's 13.4
and 9.7).

  The code below produces LREs of
13.3 for SiRstv and 9.18 for AgWt -- close to Splus
6.2, better than SAS.

  If anyone has a whole zoo of R versions lying around
they could go back and try this code on them ...

  Ben Bolker


scanstr = function(con,target,verbose=FALSE) {
  s = readLines(con,n=1)
  while (length(s)>0 && length(grep(target,s))==0) {
    if (verbose) cat(s,"\n")
    s = readLines(con,n=1)
  }
  pushBack(s,con)
}
LRE = function(target,val) {
  -log10(abs(val-target)/abs(target))
}
u1 = url("http://www.itl.nist.gov/div898/strd/anova/SiRstv.dat",open="r")
scanstr(u1,"F Statistic")  ## skip ahead to ANOVA table
r = readLines(u1,3)
## grab last value in ANOVA table
cert.F = as.numeric(strsplit(r[3]," ")[[1]][7])  
scanstr(u1,"Resistance")                         ## skip ahead to data
x = read.table(u1,skip=1)
names(x)=c("Instrument","Resistance")
x$Instrument=factor(x$Instrument)

a1 = anova(lm(Resistance~Instrument,data=x))     ## R anova(lm)
a2 = summary(aov(Resistance~Instrument,data=x))     ## R aov
R.F = a1$"F value"[1]
R.F2 = a2[[1]]$"F value"[1]
LRE(cert.F,R.F)
LRE(cert.F,R.F2)
close(u1)

u2 = url("http://www.itl.nist.gov/div898/strd/anova/AtmWtAg.dat",open="r")
scanstr(u2,"F Statistic")  ## skip ahead to ANOVA table
r = readLines(u2,4)
## grab last value in ANOVA table
cert.F = as.numeric(strsplit(r[4]," ")[[1]][7])  
scanstr(u2,"AgWt",verbose=TRUE)  ## skip ahead to data
x = read.table(u2,skip=1,header=FALSE)
names(x)=c("Instrument","AgWt")
x$Instrument=factor(x$Instrument)

a1 = anova(lm(AgWt~Instrument,data=x))     ## R anova(lm)
a2 = summary(aov(AgWt~Instrument,data=x))     ## R aov
R.F = a1$"F value"[1]
R.F2 = a2[[1]]$"F value"[1]
LRE(cert.F,R.F)
LRE(cert.F,R.F2)
close(u2)