Dlaczego istnieje wartość R ^ 2 (i co ją determinuje), gdy lm nie ma wariancji w przewidywanej wartości?

Rozważ następujący kod R:

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)

Spojrzenie na http://svn.r-project.org/R/trunk/src/appl/dqrls.f ) nie pomogło mi zrozumieć, co się dzieje, ponieważ nie znam Fortrana. W innym pytaniu odpowiedziano, że błędy tolerancji maszyny zmiennoprzecinkowej są winne za współczynniki dla X, które są bliskie, ale niezupełnie 0.

$R^2$ jest większe, gdy wartość dla coef(example(n))["X"]jest bliższa 0. Ale ...

Dlaczego w ogóle istnieje wartość ? $R^2$
Co (konkretnie) to determinuje?
Dlaczego pozornie uporządkowany postęp NaNwyników?
Dlaczego naruszenia tego postępu?
Co z tego jest „oczekiwanym” zachowaniem?

r regression russellpierce
źródło

Uwaga: R ^ 2 dla 7 powinno wynosić 0,4542, aby zobaczyć coś bardziej konstruktywnego, patrz moja odpowiedź. :-)

Cóż, szczerze mówiąc, użytkownik powinien wiedzieć coś o metodach statystycznych przed użyciem narzędzi (w przeciwieństwie do, powiedzmy, użytkowników Excela (ok, przepraszam za tani strzał)). Ponieważ jest raczej oczywiste, że R ^ 2 zbliża się do 1, gdy błąd zbliża się do zera, wiemy lepiej niż mylić wartość NaN z granicą funkcji. Gdyby pojawił się problem z rozbieżnością R ^ 2 jako ynoise -> 0 (powiedzmy, zastąp powyższą instrukcję Y Y <- rep(1,n)+runif(n)*ynoise), byłoby to interesujące :-)

Carl Witthoft

@eznme: Myślę, że wyniki są specyficzne dla maszyny lub co najmniej 32 lub 64-bitowe; Mam 32-bitową maszynę, która daje 0,1963 dla 7, ale moja 64-bitowa maszyna daje NaN. Co ciekawe, na komputerze 64-bitowym wartości R ^ 2, które nie są NaN, są bardzo zbliżone do 0,5. Kiedy o tym myślę, ma to sens, ale na początku mnie to zaskoczyło.

Aaron opuścił przepełnienie stosu

Badasz błąd zaokrąglania podwójnej precyzji. Spójrz na współczynniki; np apply(as.matrix(2:17), 1, function(n){example(n)$coefficients[-1]}). (Moje wyniki na Xeonie Win 7 x64 wynoszą od -8e-17 do + 3e-16; około połowa to prawdziwe zera.) BTW, źródło Fortran nie pomaga: to tylko opakowanie dla dqrdc; to kod, na który chcesz spojrzeć.

whuber

(Ciąg dalszy) Ale jako użytkownik wybór CV jest lepszą witryną, z tego prostego powodu, że za staranną analizę statystyczną odpowiedzialny jest użytkownik, a nie programista. Jeśli użytkownik zobaczy błędny stosunku do wielkości kanału RSS, powinien wykonać własne przetwarzanie końcowe przed dalszym zgłaszaniem. Jeśli chodzi o programowanie, chciałbym wiedzieć, jak uniknąć tych problemów numerycznych w jak największym stopniu, ale myślę, że nie da się ich uniknąć, i dlatego ważne jest, aby mieć pilnego użytkownika i edukować innych.

R^{2}

$R^2$

Iterator

Odpowiedzi:

Jak mówi Ben Bolker, odpowiedź na to pytanie znajduje się w kodzie dla summary.lm().

Oto nagłówek:

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

A więc spójrzmy x <- 1:1000; y <- rep(1,1000); z <- lm(y ~ x)na ten nieco zmodyfikowany wyciąg:

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

Zauważ, że ans r.squared to ... $0.4998923$

Aby odpowiedzieć na pytanie pytaniem: co z tego czerpiemy? :)

Wierzę, że odpowiedź leży w tym, jak R obsługuje liczby zmiennoprzecinkowe. Myślę, że mssi rsssą to sumy bardzo małych (kwadratowych) błędów zaokrąglania, stąd przyczyna wynosi około 0,5. Co do progresji, podejrzewam, że ma to związek z liczbą wartości, które trzeba czekać na +/- przybliżenia odwołać się do 0 (zarówno i , jak to prawdopodobne, że źródłem tych wartości). Nie wiem jednak, dlaczego wartości różnią się od progresji. $R^2$ mssrss0/0NaN2^(1:k)

Aktualizacja 1: Oto fajny wątek z pomocy R, który omawia niektóre powody, dla których ostrzeżenia o niedopełnieniu nie są adresowane w R.

Ponadto w tym SO Q&A znajduje się wiele interesujących postów i przydatnych linków dotyczących niedomiaru, arytmetyki o wyższej precyzji itp.

Iterator
źródło

Jestem ciekawy twojej motywacji do zadania pytania. Nie mogę wymyślić praktycznego powodu, dla którego takie zachowanie powinno mieć znaczenie; ciekawość intelektualna jest alternatywnym (i IMO o wiele bardziej rozsądnym) powodem. Myślę, że nie musisz rozumieć FORTRAN, aby odpowiedzieć na to pytanie, ale myślę, że musisz wiedzieć o rozkładzie QR i jego zastosowaniu w regresji liniowej. Jeśli traktujesz dqrlsjak czarną skrzynkę, która oblicza rozkład QR i zwraca różne informacje na jego temat, być może będziesz w stanie prześledzić kroki ... lub po prostu przejdź prosto summary.lmi prześledzić, aby zobaczyć, jak oblicza się R ^ 2. W szczególności:

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)

Następnie musisz wrócić lm.fiti zobaczyć, że dopasowane wartości są obliczane jako r1 <- y - z$residuals(tj. Jako odpowiedź minus reszty). Teraz możesz dowiedzieć się, co determinuje wartość reszt i czy wartość minus jej średnia wynosi dokładnie zero, czy nie, a następnie dowiedzieć się, jakie są wyniki obliczeń ...

Ben Bolker
źródło

Intelektualna ciekawość to większość powodów mojego pytania. Kolega zgłosił to zachowanie, a ja chciałem się z nimi pogrzebać i sprawdzić, czy uda mi się to rozgryźć. Po prześledzeniu problemu poza moim zestawem umiejętności postanowiłem zadać pytanie. Jako praktyczny problem, czasami analizy są wykonywane partiami lub występują inne błędy, a to zachowanie wydaje mi się zdecydowanie „dziwne”.

russellpierce

mms i rss są wynikami z, który jest nazwą obiektu lm w pliku Summary.lm. Tak więc odpowiedź prawdopodobnie wymaga wyjaśnienia rozkładu QR, jego zastosowania w regresji liniowej, a konkretnie niektórych szczegółów dekompozycji QR przedstawionych w kodzie leżącym u podstaw R, aby wyjaśnić, dlaczego rozkład QR kończy się na przybliżeniu 0 zamiast samego 0 .

russellpierce

@drknexus Nie zgadzam się. Rozkład QR jest jednym z wielu algorytmów numerycznych; jeśli podstawowym problemem jest precyzja liczbowa, pojawi się to w QR, mnożeniu macierzy, rozwiązaniach nieliniowych i wielu innych miejscach. Zasadnicza sekwencja jest prosta: współczynniki są bardzo nieznacznie wyłączone (powinny wynosić (0,1)); nie jest to nierozsądne, ale produkuje mssi rss„szumu”. Jest to zasada GIGO, która zapewnia, że jest dokładna, ale niepoprawna. Wolę wstawić „wykrywacz śmieci” przed obliczeniem niż zmodyfikować algo QR, ponieważ wątpię, aby jego poprawność mogła zostać poprawiona.

R^{2}

$R^2$

R^{2}

$R^2$

Iterator

Wydaje mi się, że wykrywacz śmieci powinien znajdować się na QR lub tuż przed nim. Prosta kontrola poczytalności wariancji Y i ostrzeżenie, że Y nie ma wariancji, byłaby w porządku (mogę napisać opakowanie lm dla moich przyjaciół, które to właśnie robi). Wydaje mi się, że zanim obliczysz , jest już zbyt daleko od obliczeniowej nory królika, aby wiedzieć, czy patrzysz na śmieci, czy nie.

R^{2}

$R^2$

russellpierce

$R^2$ jest zdefiniowane jako ( http://en.wikipedia.org/wiki/R_squared ), więc jeśli suma kwadratów-suma wynosi 0, to jest niezdefiniowana. Moim zdaniem R powinien pokazać komunikat o błędzie. $R^2 = 1-\frac{\textrm{SS}_{err}}{\textrm{SS}_{tot}}$

Bernd Elkemann
źródło

Czy możesz podać praktyczną sytuację, w której takie zachowanie miałoby znaczenie?

Ben Bolker

@Brandon - Iterator umieścił tam buźkę, a ty wciąż krzyczysz!

Carl Witthoft

@eznme Chociaż błąd jest dobry, dość trudno jest uchwycić wszelkiego rodzaju miejsca, w których pojawiają się problemy z liczbą zmiennoprzecinkową, zwłaszcza w świecie arytmetyki IEEE-754. Lekcja polega na tym, że nawet obliczenia R na chleb i masło powinny być traktowane delikatnie.

Iterator

Te rozważania są szczególnie ważne, ponieważ w swoich pismach John Chambers (jeden z twórców S, a zatem „dziadek” R) mocno podkreśla użycie R do niezawodnego przetwarzania. Np. Patrz Chambers, Oprogramowanie do analizy danych: Programowanie w języku R (Springer Verlag 2008): „obliczenia i oprogramowanie do analizy danych powinny być godne zaufania: powinny robić to, co twierdzą i być do tego zobowiązani”. [W p. 3.]

whuber

Problem polega na tym, że na lepsze lub gorsze, rdzeń R jest odporny na (jak to widzą) festonowanie kodu wieloma, wieloma kontrolami, które przechwytują wszystkie przypadki narożne i możliwe dziwne błędy użytkowników - obawiają się (myślę), że to (a) zajmie dużo czasu, (b) sprawi, że baza kodu będzie znacznie większa i trudniejsza do odczytania (ponieważ dosłownie są tysiące takich specjalnych przypadków), oraz (c) spowolni wykonanie poprzez wymuszanie takich kontroli przez cały czas nawet w sytuacjach, w których obliczenia są powtarzane wiele, wiele razy.

Ben Bolker