2012-02-11 24 views
6

Considere el siguiente código R (que, creo, finalmente llama a algunos Fortran):¿Por qué lm devuelve valores cuando no hay varianza en el valor predicho?

X <- 1:1000 
Y <- rep(1,1000) 
summary(lm(Y~X)) 

¿Por qué valores devueltos por resumen? ¿No debería este modelo no encajar ya que no hay varianza en Y? Más importante aún, ¿por qué el modelo R^2 ~ = .5?

Editar

Rastreé el código de lm a lm.fit y puede ver esta llamada:

z <- .Fortran("dqrls", qr = x, n = n, p = p, y = y, ny = ny, 
    tol = as.double(tol), coefficients = mat.or.vec(p, ny), residuals = y, 
    effects = y, rank = integer(1L), pivot = 1L:p, qraux = double(p), 
    work = double(2 * p), PACKAGE = "base") 

Eso es donde el ajuste real parece suceder. En cuanto a http://svn.r-project.org/R/trunk/src/appl/dqrls.f) no me ayudó a entender lo que está pasando, porque no sé fortran.

+1

Ah, la R^2 de 0.5 es una pregunta bastante interesante. – Iterator

+0

Creo que voy a girarlo como una pregunta separada, entonces ... – russellpierce

Respuesta

5

Estadísticamente hablando, ¿qué deberíamos anticipar (me gustaría decir "esperar", pero ese es un término muy específico ;-))? Los coeficientes deben ser (0,1), en lugar de "no ajustarse". La covarianza de (X, Y) se supone proporcional a la varianza de X, y no al revés. Como X tiene una varianza distinta de cero, no hay problema. Como la covarianza es 0, el coeficiente estimado para X debe ser 0. Entonces, dentro de la tolerancia de la máquina, esta es la respuesta que está obteniendo.

No hay anomalía estadística aquí. Puede haber un malentendido estadístico. También está la cuestión de la tolerancia de la máquina, pero un coeficiente del orden de 1E-19 es bastante insignificante, dada la escala de los valores de predicción y respuesta.

Actualización 1: se puede encontrar una revisión rápida de la regresión lineal simple en this Wikipedia page. La clave a tener en cuenta es que Var(x) está en el denominador, Cov(x,y) en el numerador. En este caso, el numerador es 0, el denominador no es cero, por lo que no hay motivos para esperar un NaN o NA. Sin embargo, uno puede preguntarse por qué no es el coeficiente resultante para x a 0, y eso tiene que ver con problemas de precisión numérica de la descomposición QR.

+0

Veo su punto (s). La tolerancia de la máquina está más cerca de 1E-17 para problemas de N más pequeños, pero aún "insignificante". Supongo que esperaba que la función simplemente fallara como lo hace cuando N = 4 (pero, nuevamente, (para mí) extrañamente no falla para N = 3). – russellpierce

2

Creo que esto es simplemente porque la descomposición QR se implementa con aritmética de punto flotante.

El parámetro singular.ok en realidad se refiere a la matriz de diseño (es decir, solo a X). Trate

lm.fit(cbind(X, X), Y) 

vs

lm.fit(cbind(X, X), Y, singular.ok=F) 
2

Estoy de acuerdo en que el problema podría ser de punto flotante. pero no creo que sea singularidad

Si marca usando solve(t(x1)%*%x1)%*%(t(x1)%*%Y) en lugar de QR, (t(x1)%*%x1) no es singular

uso x1 = cbind(rep(1,1000,X) porque lm(Y~X) incluye la intersección.

Cuestiones relacionadas