¿Por qué hay un valor R ^ 2 (y qué lo determina) cuando lm no tiene variación en el valor predicho?

10

Considere el siguiente código 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)

Mirar http://svn.r-project.org/R/trunk/src/appl/dqrls.f ) no me ayudó a entender lo que está sucediendo, porque no conozco Fortran. En otra pregunta , se respondió que los errores de tolerancia de máquina de coma flotante eran los culpables de los coeficientes para X que están cerca, pero no del todo 0.

es mayor cuando el valor deestá más cerca de 0. Pero ...R2coef(example(n))["X"]

  1. ¿Por qué hay un valor ? R2
  2. ¿Qué (específicamente) lo está determinando?
  3. ¿Por qué la progresión aparentemente ordenada de los NaNresultados?
  4. ¿Por qué las violaciones de esa progresión?
  5. ¿Qué de esto es el comportamiento 'esperado'?
russellpierce
fuente
Nota: R ^ 2 de 7 debería ser 0.4542 para ver algo más constructivo, vea mi respuesta. :-)
1
Bueno, para ser justos, se supone que el usuario realmente debe saber algo acerca de los métodos estadísticos antes de usar herramientas (a diferencia de, por ejemplo, los usuarios de Excel (ok, perdón por el tiro barato)). Como es bastante obvio que R ^ 2 se acerca a 1 cuando el error se acerca a cero, sabemos mejor que confundir un valor de NaN con el límite de una función. Ahora, si hubiera un problema con R ^ 2 divergiendo como ynoise -> 0 (digamos, reemplace la declaración Y anterior con Y <- rep(1,n)+runif(n)*ynoise), eso sería interesante :-)
Carl Witthoft
@eznme: Creo que los resultados son específicos de la máquina, o al menos de 32 o 64 bits específicos; Tengo una máquina de 32 bits que da 0.1963 por 7, pero mi máquina de 64 bits da NaN. Curiosamente, en la máquina de 64 bits, los R ^ 2 que no son NaN están muy cerca de 0.5. Tiene sentido cuando lo pienso, pero al principio me sorprendió.
Aaron dejó Stack Overflow
1
Estás estudiando un error de redondeo de doble precisión. Echa un vistazo a los coeficientes; por ejemplo, apply(as.matrix(2:17), 1, function(n){example(n)$coefficients[-1]}). (Mis resultados, en un Win 7 x64 Xeon, varían de -8e-17 a + 3e-16; aproximadamente la mitad son ceros verdaderos). Por cierto, la fuente Fortran no es de ayuda: es solo un contenedor para dqrdc; ese es el código que quieres mirar.
whuber
1
R2

Respuestas:

6

Como dice Ben Bolker, la respuesta a esta pregunta se puede encontrar en el código para summary.lm().

Aquí está el encabezado:

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

Entonces, x <- 1:1000; y <- rep(1,1000); z <- lm(y ~ x)eche un vistazo a este extracto ligeramente modificado:

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

0.4998923

Para responder una pregunta con una pregunta: ¿qué extraemos de esto? :)

mssrssR2mssrss0/0NaN2^(1:k)


Actualización 1: Aquí hay un buen hilo de R-help que aborda algunas de las razones por las cuales las advertencias de subflujo no se abordan en R.

Además, este SO Q&A tiene una serie de publicaciones interesantes y enlaces útiles sobre flujo inferior, aritmética de mayor precisión, etc.

Iterador
fuente
8

Tengo curiosidad acerca de su motivación para hacer la pregunta. No puedo pensar en una razón práctica por la que este comportamiento debería importar; La curiosidad intelectual es una razón alternativa (y la OMI mucho más sensata). Creo que no necesita comprender FORTRAN para responder esta pregunta, pero creo que necesita saber sobre la descomposición de QR y su uso en la regresión lineal. Si lo trata dqrlscomo un cuadro negro que calcula una descomposición QR y devuelve información diversa al respecto, entonces puede seguir los pasos ... o simplemente ir directamente summary.lmy rastrear para ver cómo se calcula el R ^ 2. En 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)

Luego debe volver a lm.fitver que los valores ajustados se calculan como r1 <- y - z$residuals(es decir, como la respuesta menos los residuos). Ahora puede averiguar qué determina el valor de los residuos y si el valor menos su media es exactamente cero o no, y a partir de ahí averiguar los resultados computacionales ...

Ben Bolker
fuente
La curiosidad intelectual es la mayor parte de la razón de mi pregunta. Un colega informó sobre el comportamiento y quería curiosear y ver si podía resolverlo. Después de rastrear el problema más allá de mi conjunto de habilidades, decidí hacer la pregunta. Como cuestión práctica, a veces los análisis se realizan por lotes, u ocurren otros errores, y este comportamiento me parece decididamente "extraño".
russellpierce
1
mms y rss son ambos resultados de z, que es el nombre del objeto lm dentro de summary.lm. Por lo tanto, una respuesta probablemente requiera una explicación de la descomposición QR, su uso en regresión lineal y, específicamente, algunos detalles de la descomposición QR como se instancia en el código subyacente R para explicar por qué la descomposición QR termina con aproximaciones de 0 en lugar de 0 .
russellpierce
mssrssR2R2
R2
0

R2R2=1SSerrSStot

Bernd Elkemann
fuente
1
¿Puedes dar una situación práctica donde este comportamiento importaría?
Ben Bolker
3
@Brandon: ¡Iterator puso la carita allí y todavía te dieron una paliza!
Carl Witthoft
2
@eznme Si bien un error es bueno, es bastante difícil detectar todo tipo de lugares donde surgen problemas de coma flotante, especialmente en el mundo de la aritmética IEEE-754. La lección aquí es que incluso los cálculos de pan y mantequilla con R deben manejarse con delicadeza.
Iterator
2
Estas consideraciones son especialmente importantes porque en sus escritos, John Chambers (uno de los creadores de S y, por lo tanto, un "abuelo" de R) enfatiza fuertemente el uso de R para computación confiable. Por ejemplo, vea Chambers, Software for Data Analysis: Programming with R (Springer Verlag 2008): "los cálculos y el software para el análisis de datos deben ser confiables: deben hacer lo que dicen y ver que lo hagan". [En la p. 3.]
whuber
2
El problema es que, para bien o para mal, R-core es resistente (como lo ven) a engañar el código con muchas, muchas comprobaciones que interceptan todos los casos de esquina y posibles errores de usuarios extraños; tienen miedo (creo) de que (a) tomarán grandes cantidades de tiempo, (b) hará que la base del código sea mucho más grande y difícil de leer (porque hay literalmente miles de estos casos especiales), y (c) ralentizará la ejecución al forzar dichos controles todo el tiempo incluso en situaciones donde los cálculos se repiten muchas, muchas veces.
Ben Bolker