Почему lm возвращает значения, если в прогнозируемом значении нет дисперсии? - PullRequest
Новая
       26

Почему lm возвращает значения, если в прогнозируемом значении нет дисперсии?

6 голосов
/ 12 февраля 2012

Рассмотрим следующий код R (который, я думаю, в конечном итоге вызывает некоторый Фортран): 

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

Почему значения возвращаются в виде сводки? Не должна ли эта модель не соответствовать, так как нет Y-дисперсии? Более важно, почему модель R ^ 2 ~ = .5?

Редактировать

Я отследил код от lm до lm.fit и вижу этот вызов: 

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

Вот тут-то и происходит фактическая подгонка. Просмотр http://svn.r -project.org / R / trunk / src / appl / dqrls.f ) не помог мне понять, что происходит, потому что я не знаю фортран.

Ответы [ 3 ]

5 голосов
/ 12 февраля 2012

Статистически говоря, что мы должны ожидать (я хотел бы сказать «ожидать», но это очень специфический термин ;-))? Коэффициенты должны быть (0,1), а не "не соответствовать". Ковариация (X, Y) предполагается пропорциональной дисперсии X, а не наоборот. Поскольку X имеет ненулевую дисперсию, проблем нет. Поскольку ковариация равна 0, предполагаемый коэффициент для X должен быть равен 0. Таким образом, в пределах допуска машины это ответ, который вы получаете.

Здесь нет статистической аномалии. Там может быть статистическое недоразумение. Есть также проблема машинной толерантности, но коэффициент порядка 1E-19 довольно незначителен, учитывая масштаб предиктора и значения ответа.

Обновление 1: краткий обзор простой линейной регрессии можно найти на этой странице Википедии . Ключевым моментом, который стоит отметить, является то, что Var(x) находится в знаменателе, Cov(x,y) в числителе. В этом случае числитель равен 0, знаменатель не равен нулю, поэтому нет оснований ожидать NaN или NA. Однако можно спросить, почему не получается результирующий коэффициент для x a 0, и это связано с проблемами численной точности QR-разложения.

2 голосов
/ 12 февраля 2012

Я согласен, что проблема может быть с плавающей точкой.но я не думаю, что это единственное число.

Если вы проверяете, используя solve(t(x1)%*%x1)%*%(t(x1)%*%Y) вместо QR, (t(x1)%*%x1) не является единственным числом

, используйте x1 = cbind(rep(1,1000,X), поскольку lm(Y~X) включает перехват

2 голосов
/ 12 февраля 2012

Я считаю, что это просто потому, что декомпозиция QR реализована с использованием арифметики с плавающей запятой.

Параметр singular.ok фактически относится к проектной матрице (т.е. только X). Попробуйте 

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

против. 

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