Warum gibt es einen R ^ 2-Wert (und was bestimmt ihn), wenn lm keine Varianz im vorhergesagten Wert aufweist?

Betrachten Sie den folgenden 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)

Ein Blick auf http://svn.r-project.org/R/trunk/src/appl/dqrls.f ) hat mir nicht geholfen zu verstehen, was los ist, da ich Fortran nicht kenne. In einer anderen Frage wurde beantwortet, dass Gleitkomma-Maschinentoleranzfehler für Koeffizienten für X verantwortlich sind, die nahe, aber nicht ganz 0 sind.

$R^2$ ist größer, wenn der Wert für coef(example(n))["X"]näher bei 0 liegt. Aber ...

Warum gibt es überhaupt einen Wert? $R^2$
Was bestimmt es (konkret)?
Warum das scheinbar geordnete Fortschreiten der NaNErgebnisse?
Warum die Verstöße gegen diesen Fortschritt?
Was ist das "erwartete" Verhalten?

r regression russellpierce
quelle

Hinweis: 7s R ^ 2 sollte 0,4542 sein, um etwas Konstruktiveres zu sehen, siehe meine Antwort. :-)

Um fair zu sein, sollte der Benutzer tatsächlich etwas über statistische Methoden wissen , bevor er Tools verwendet (im Gegensatz zu beispielsweise Excel-Benutzern (ok, Entschuldigung für den billigen Schuss)). Da es ziemlich offensichtlich ist, dass sich R ^ 2 1 nähert, wenn sich der Fehler Null nähert, wissen wir besser, als einen NaN-Wert mit der Grenze einer Funktion zu verwechseln. Wenn es ein Problem mit R ^ 2 geben würde, das als ynoise -> 0 divergiert (sagen Sie, ersetzen Sie die Y-Anweisung oben durch Y <- rep(1,n)+runif(n)*ynoise), wäre das interessant :-)

Carl Witthoft

@eznme: Ich denke, die Ergebnisse sind maschinenspezifisch oder mindestens 32- oder 64-Bit-spezifisch. Ich habe eine 32-Bit-Maschine, die 0,1963 für 7 gibt, aber meine 64-Bit-Maschine gibt NaN. Interessanterweise liegen auf der 64-Bit-Maschine die R ^ 2s, die nicht NaN sind, alle sehr nahe bei 0,5. Sinnvoll, wenn ich darüber nachdenke, aber es hat mich zuerst überrascht.

Aaron verließ Stack Overflow

Sie untersuchen Rundungsfehler mit doppelter Genauigkeit. Schauen Sie sich die Koeffizienten an. zB , apply(as.matrix(2:17), 1, function(n){example(n)$coefficients[-1]}). (Meine Ergebnisse auf einem Win 7 x64 Xeon reichen von -8e-17 bis + 3e-16; ungefähr die Hälfte sind echte Nullen.) Übrigens ist die Fortran-Quelle keine Hilfe: Es ist nur ein Wrapper für dqrdc; Das ist der Code, den Sie sich ansehen möchten.

whuber

(Fortsetzung) Als Benutzer ist die Wahl des Lebenslaufs jedoch eine bessere Website, da die sorgfältige statistische Analyse in der Verantwortung des Benutzers und nicht des Entwicklers liegt. Wenn der Benutzer ein fehlerhaftes Bezug auf die Größe des RSS sieht , sollte er seine eigene Nachbearbeitung durchführen, bevor er weiter berichtet. In Bezug auf die Programmierung würde ich gerne wissen, wie diese numerischen Probleme so weit wie möglich vermieden werden können, aber ich denke, dass sie nicht vermieden werden können, und hier ist es wichtig, einen sorgfältigen Benutzer zu haben und andere zu schulen.

R^{2}

$R^2$

Iterator

Antworten:

Wie Ben Bolker sagt, finden Sie die Antwort auf diese Frage im Code für summary.lm().

Hier ist der Header:

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

x <- 1:1000; y <- rep(1,1000); z <- lm(y ~ x)Schauen wir uns also diesen leicht modifizierten Auszug an:

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

Beachten Sie, dass Ans r.squared ist ... $0.4998923$

Um eine Frage mit einer Frage zu beantworten: Was ziehen wir daraus? :) :)

Ich glaube, die Antwort liegt darin, wie R mit Gleitkommazahlen umgeht. Ich denke das mssund rsssind die Summen sehr kleiner (quadratischer) Rundungsfehler, daher liegt der Grund für bei etwa 0,5. Ich vermute, dass dies mit der Anzahl der Werte zusammenhängt, die erforderlich sind, damit sich die +/- Näherungen auf 0 aufheben (für beide und , wie wahrscheinlich die Quelle dieser Werte ist). Ich weiß jedoch nicht, warum sich die Werte von einer Progression unterscheiden. $R^2$ mssrss0/0NaN2^(1:k)

Update 1: Hier ist ein netter Thread von R-help , der einige der Gründe behandelt, warum Unterlaufwarnungen in R nicht behandelt werden.

Darüber hinaus enthält diese SO Q & A eine Reihe interessanter Beiträge und nützlicher Links zu Unterlauf, Arithmetik mit höherer Genauigkeit usw.

Iterator
quelle

Ich bin gespannt auf Ihre Motivation, die Frage zu stellen. Ich kann mir keinen praktischen Grund vorstellen, warum dieses Verhalten von Bedeutung sein sollte. Intellektuelle Neugier ist ein alternativer (und IMO viel vernünftigerer) Grund. Ich denke, Sie müssen FORTRAN nicht verstehen, um diese Frage zu beantworten, aber ich denke, Sie müssen etwas über die QR-Zerlegung und ihre Verwendung in der linearen Regression wissen. Wenn Sie dqrlsals Black Box behandeln, die eine QR-Zerlegung berechnet und verschiedene Informationen darüber zurückgibt, können Sie möglicherweise die Schritte verfolgen ... oder gehen Sie einfach direkt zu summary.lmund verfolgen Sie, um zu sehen, wie der R ^ 2 berechnet wird. Speziell:

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)

Dann müssen Sie zurückgehen lm.fitund sehen, dass die angepassten Werte berechnet werden als r1 <- y - z$residuals(dh als Antwort abzüglich der Residuen). Jetzt können Sie herausfinden, was den Wert der Residuen bestimmt und ob der Wert minus seines Mittelwerts genau Null ist oder nicht, und von dort aus die Berechnungsergebnisse herausfinden ...

Ben Bolker
quelle

Intellektuelle Neugier ist der Hauptgrund für meine Frage. Ein Kollege berichtete über das Verhalten und ich wollte herumstöbern und sehen, ob ich es herausfinden konnte. Nachdem ich das Problem über meine Fähigkeiten hinaus verfolgt hatte, beschloss ich, die Frage zu stellen. In der Praxis werden Analysen manchmal stapelweise durchgeführt, oder es treten andere Fehler auf, und dieses Verhalten erscheint mir ausgesprochen „seltsam“.

Russellpierce

mms und rss sind beide Ergebnisse von z, dem Namen des lm-Objekts in summary.lm. Eine Antwort erfordert daher wahrscheinlich eine Erklärung der QR-Zerlegung, ihrer Verwendung bei der linearen Regression und insbesondere einiger Details der QR-Zerlegung, wie sie im Code zugrunde liegen, der R zugrunde liegt, um zu erklären, warum die QR-Zerlegung mit Annäherungen von 0 statt 0 endet .

Russellpierce

@drknexus Ich bin anderer Meinung. Die QR-Dekomposition ist einer von vielen numerischen Algorithmen. Wenn das zugrunde liegende Problem die numerische Genauigkeit ist, tritt dies in QR, Matrixmultiplikation, nichtlinearen Lösern und so vielen anderen Stellen auf. Die wesentliche Reihenfolge ist einfach: Die Koeffizienten sind sehr geringfügig abweichen (sollte (0,1) sein); Dies ist nicht unangemessen, erzeugt jedoch das mssund rss"Rauschen". Es ist das GIGO-Prinzip, das sicherstellt, dass präzise, aber falsch ist. Ich würde lieber einen "Mülldetektor" einfügen, bevor ich berechne, als das QR-Algo zu modifizieren, da ich bezweifle, dass seine Gültigkeit verbessert werden könnte.

R^{2}

$R^2$

R^{2}

$R^2$

Iterator

Es scheint mir, dass der Mülldetektor am QR oder direkt davor sein sollte. Eine einfache Überprüfung der Varianz von Y und die Warnung, dass Y keine Varianz aufweist, wäre in Ordnung (ich kann einen Film-Wrapper für meine Freunde schreiben, der genau dies tut). Es scheint mir, dass man zu dem Zeitpunkt, an dem man berechnet , bereits zu weit unten im rechnerischen Kaninchenbau ist, um zu wissen, ob man Müll betrachtet oder nicht.

R^{2}

$R^2$

Russellpierce

$R^2$ ist definiert als ( http://en.wikipedia.org/wiki/R_squared ). Wenn also die Summe der Quadrate 0 ist, ist sie undefiniert. Meiner Meinung nach sollte R eine Fehlermeldung anzeigen. $R^2 = 1-\frac{\textrm{SS}_{err}}{\textrm{SS}_{tot}}$

Bernd Elkemann
quelle

Können Sie eine praktische Situation angeben, in der dieses Verhalten von Bedeutung wäre?

Ben Bolker

@Brandon - Iterator hat den Smiley da reingelegt und du hast immer noch einen Rausch!

Carl Witthoft

@eznme Obwohl ein Fehler gut ist, ist es ziemlich schwierig, alle Arten von Stellen zu erfassen, an denen Gleitkommaprobleme auftreten, insbesondere in der Welt der IEEE-754-Arithmetik. Die Lehre hier ist, dass selbst die Brot-und-Butter-Berechnungen mit R vorsichtig gehandhabt werden sollten.

Iterator

Diese Überlegungen sind besonders wichtig, da John Chambers (einer der Urheber von S und daher ein "Großvater" von R) in seinen Schriften die Verwendung von R für zuverlässiges Rechnen nachdrücklich betont . Siehe z. B. Chambers, Software für die Datenanalyse: Programmieren mit R (Springer Verlag 2008): "Die Berechnungen und die Software für die Datenanalyse sollten vertrauenswürdig sein: Sie sollten das tun, was sie behaupten, und dies auch tun." [Auf S. 3.]

whuber

Das Problem ist, dass R-Core zum Guten oder Schlechten dagegen ist, den Code (wie sie sehen) mit vielen, vielen Überprüfungen zu versehen, die alle Eckfälle und mögliche seltsame Benutzerfehler abfangen - sie befürchten (glaube ich), dass dies der Fall ist wird (a) sehr viel Zeit in Anspruch nehmen, (b) die Codebasis viel größer und schwerer lesbar machen (weil es buchstäblich Tausende dieser Sonderfälle gibt) und (c) die Ausführung verlangsamen, indem solche Überprüfungen ständig erzwungen werden selbst in Situationen, in denen Berechnungen viele, viele Male wiederholt werden.

Ben Bolker