up vote 8 down vote favorite
1
share [g+] share [fb]

Consider the following R code:

example <- function(n) {
    X <- 1:n
    Y <- rep(1,n)
    return(lm(Y~X))
}
#(2.13.0, i386-pc-mingw32)
summary(example(7))    #R^2 = .1963
summary(example(62))   #R^2 = .4529
summary(example(4540)) #R^2 = .7832
summary(example(104))) #R^2 = 0
#I did a search for n 6:10000, the result for R^2 is NaN for
#n = 2, 4, 16, 64, 256, 1024, 2085 (not a typo), 4096, 6175 (not a typo), and 8340 (not a typo)

Looking at http://svn.r-project.org/R/trunk/src/appl/dqrls.f) did not help me understand what is going on, because I do not know Fortran. In another question it was answered that floating point machine tolerance errors were are to blame for coefficients for X that are close to, but not quite 0.

$R^2$ is greater when the value for coef(example(n))["X"] is closer to 0. But...

  1. Why is there an $R^2$ value at all?
  2. What (specifically) is determining it?
  3. Why the seemingly orderly progression of NaN results?
  4. Why the violations of that progression?
  5. What of this is 'expected' behavior?
link|improve this question
94% accept rate
Note: 7's R^2 should be 0.4542 to see something more constructive see my answer. :-) – eznme Feb 14 at 18:49
Well, to be fair, the user is supposed to actually know something about statistical methods before using tools (unlike, say, Excel users (ok, sorry about the cheap shot)). Since it's rather obvious that R^2 approaches 1 as error approaches zero, we know better than to confuse a NaN value with the limit of a function. Now, if there were a problem with R^2 diverging as ynoise-->0 (say, replace Y statement above with Y <- rep(1,n)+runif(n)*ynoise ), that would be interesting :-) – Carl Witthoft Feb 14 at 19:45
@eznme: I think the results are machine specific, or at least 32 or 64 bit specific; I have a 32-bit machine that gives 0.1963 for 7, but my 64-bit machine gives NaN. Interestingly, on the 64-bit machine, the R^2s that are not NaN are all very close to 0.5. Makes sense when I think about it, but it surprised me at first. – Aaron Feb 14 at 20:35
1  
You're studying double precision rounding error. Take a look at the coefficients; e.g., apply(as.matrix(2:17), 1, function(n){example(n)$coefficients[-1]}). (My results, on a Win 7 x64 Xeon, range from -8e-17 to +3e-16; about half are true zeros.) BTW, the Fortran source is of no help: it's just a wrapper for dqrdc; that's the code you want to look at. – whuber♦ Feb 14 at 21:30
1  
(Continued) But, as a user, the choice of CV is a better site, for the simple reason that diligent statistical analysis is the responsibility of the user, not the developer. If the user sees an erroneous $R^2$ relative to the magnitude of the RSS, then they should do their own post-processing before reporting further. Programming-wise, I'd like to know how to avoid these numerical issues as much as possible, but I think that they can't be escaped, and that's where it's important to have a diligent user and to educate others. – Iterator Feb 14 at 22:53
show 4 more comments
feedback

migrated from stackoverflow.com Feb 14 at 19:59

This question came from our site for professional and enthusiast programmers.

3 Answers

up vote 6 down vote accepted

As Ben Bolker says, the answer to this question can be found in the code for summary.lm().

Here's the header:

function (object, correlation = FALSE, symbolic.cor = FALSE, 
    ...) 
{

So, let x <- 1:1000; y <- rep(1,1000); z <- lm(y ~ x) and then take a look at this slightly modified extract:

    p <- z$rank
    rdf <- z$df.residual
    Qr <- stats:::qr.lm(z)
    n <- NROW(Qr$qr)
    r <- z$residuals
    f <- z$fitted.values
    w <- z$weights
    if (is.null(w)) {
        mss <- sum((f - mean(f))^2)
        rss <- sum(r^2)
    }
    ans <- z[c("call", "terms")]
    if (p != attr(z$terms, "intercept")) {
        df.int <- 1L
        ans$r.squared <- mss/(mss + rss)
        ans$adj.r.squared <- 1 - (1 - ans$r.squared) * ((n - 
            df.int)/rdf)
    }

Notice that ans\$r.squared is $0.4998923$...

To answer a question with a question: what do we draw from this? :)

I believe the answer lies in how R handles floating point numbers. I think that mss and rss are the sums of very small (squared) rounding errors, hence the reason $R^2$ is about 0.5. As for the progression, I suspect this has to do with the number of values that it takes for the +/- approximations to cancel out to 0 (for both mss and rss, as 0/0 is likely the source of these NaN values). I don't know why the values differ from a 2^(1:k) progression, though.

Update 1: Here is a nice thread from R-help addressing some of the reasons that underflow warnings are not addressed in R.

In addition, this SO Q&A has a number of interesting posts and useful links regarding underflow, higher precision arithmetic, etc.

link|improve this answer
feedback
up vote 8 down vote

I'm curious about your motivation for asking the question. I can't think of a practical reason this behavior should matter; intellectual curiosity is an alternative (and IMO much more sensible) reason. I think you don't need to understand FORTRAN to answer this question, but I think you do need to know about QR decomposition and its use in linear regression. If you treat dqrls as a black box that computes a QR decomposition and returns various information about it, then you may be able to trace the steps ... or just go straight to summary.lm and trace through to see how the R^2 is calculated. In particular:

mss <- if (attr(z$terms, "intercept")) 
          sum((f - mean(f))^2)
       else sum(f^2)
rss <- sum(r^2)
## ... stuff ...
ans$r.squared <- mss/(mss + rss)

Then you have to go back into lm.fit and see that the fitted values are computed as r1 <- y - z$residuals (i.e. as the response minus the residuals). Now you can go figure out what determines the value of the residuals and whether the value minus its mean is exactly zero or not, and from there figure out the computational outcomes ...

link|improve this answer
Intellectual curiosity is the majority of the reason for my question. A colleague reported the behavior and I wanted to poke around and see if I could figure it out. After I traced the issue beyond my skill-set, I decided to ask the question. As a practical issue, sometimes analyses are done by batch, or other errors occur, and this behavior strikes me as decidedly 'odd'. – drknexus Feb 14 at 21:47
1  
mms and rss are both results of z, which is the name of the lm object inside of summary.lm. So, an answer probably does require an explanation of QR decomposition, its use in linear regression, and specifically some details QR decomposition as instantiated in the code underlying R in order to explain why the QR decomposition ends up with approximations of 0 rather than 0 itself. – drknexus Feb 14 at 21:50
@drknexus I disagree. QR decomp is one of many numerical algorithms; if the underlying issue is numerical precision, this'll crop up in QR, matrix multiplication, non-linear solvers, and so many other places. The essential sequence is simple: the coefficients are very slightly off (should be (0,1)); this isn't unreasonable, yet produces the mss and rss "noise". It's the GIGO principle that assures that $R^2$ is precise, but incorrect. I'd rather insert a "garbage detector" before calculating $R^2$ than to modify the QR algo, because I doubt its validity could be improved. – Iterator Feb 14 at 23:02
It seems to me that the garbage detector should be at the QR or right before it. A simple sanity check on the variance of Y and warning that Y lacks variance would be fine (I may write an lm wrapper for my friends that does just this). It seems to me that by the time you are calculating $R^2$, one is already too far down the computational rabbit hole to know if one is looking at garbage or not. – drknexus Feb 15 at 22:00
feedback
up vote 0 down vote

$R^2$ is defined as $R^2 = 1-\frac{\textrm{SS}_{err}}{\textrm{SS}_{tot}}$ ( http://en.wikipedia.org/wiki/R_squared ), so if the sum-of-squares-total is 0 then it is undefined. In my opinion R should show an error-message.

link|improve this answer
1  
Can you give a practical situation where this behavior would matter? – Ben Bolker Feb 14 at 20:12
3  
@Brandon -- Iterator put the smiley in there and you still got whooshed! – Carl Witthoft Feb 14 at 20:45
2  
@eznme While an error is good, it is quite difficult to catch all kinds of places where floating point issues arise, notably in the world of IEEE-754 arithmetic. The lesson here is that even the bread and butter calculations with R should be handled delicately. – Iterator Feb 14 at 21:08
2  
These considerations are especially important because in his writings, John Chambers (one of the originators of S and therefore a "grandfather" of R) strongly emphasizes the use of R for reliable computing. E.g., see Chambers, Software for Data Analysis: Programming with R (Springer Verlag 2008): "the computations and the software for data analysis should be trustworthy: they should do what they claim, and be seen to do so." [At p. 3.] – whuber♦ Feb 14 at 21:37
2  
The problem is that for better or worse, R-core is resistant to (as they see it) festooning the code with many, many checks that intercept all corner cases and possible weird user errors -- they are afraid (I think) that it will (a) take huge amounts of their time, (b) make the code base much larger and harder to read (because there are literally thousands of these special cases), and (c) slow down execution by forcing such checks all the time even in situations where computations are being repeated many, many times. – Ben Bolker Feb 14 at 22:12
show 9 more comments
feedback

Your Answer

 
or
required, but never shown
discard

Not the answer you're looking for? Browse other questions tagged or ask your own question.