residual - ¿Por qué lm devuelve valores cuando no hay varianza en el valor predicho?

regresion lineal

Question

Considere el siguiente código R (que, creo, finalmente llama Fortran):

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

¿Por qué los valores son devueltos por un 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

Seguí el código de lm a lm.fit y puedo 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")

Ahí es donde parece que ocurre el ajuste real. Ver http://svn.r-project.org/R/trunk/src/appl/dqrls.f ) no me ayudó a entender lo que está pasando, porque no sé por qué.

Answer 1

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 del predictor y los valores de respuesta.

Actualización 1: se puede encontrar una revisión rápida de la regresión lineal simple en esta página de Wikipedia . 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 razón 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.

Answer 2

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

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

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

vs.

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

Answer 3

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

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

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