Играех си с функцията lm() на R и забелязах нещо странно, когато стигнах до изчисляването на коефициента на детерминация R^ 2. Да предположим, че имам проблем с регресията на играчката:
set.seed(0); N= 12;
x_ <- 1:N; y_ <- 2* x_+ rnorm(N, sd=2);
lm1_ <- lm( y_ ~ x_ ); S1 <- summary(lm1_);
lm2_ <- lm( y_ ~ -1 + model.matrix(lm1_) ); S2 <- summary(lm2_);
Така че в този момент всъщност задаваме lm1_ и lm2_ като едно и също нещо. Използва се същата зависима променлива, използва се същата матрица на модела и трябва да получим същото точно напасване. Нека проверим тези:
identical( fitted(lm2_),fitted(lm1_))
[1] TRUE #SURE!
identical( residuals(lm2_), residuals(lm1_))
[1] TRUE #YEAP!
identical( lm1_$df.residual, lm2_$df.residual)
[1] TRUE #SORTED!
identical( as.numeric(model.matrix(lm2_ )), as.numeric(model.matrix(lm1_ )))
[1] TRUE #WHY NOT?
identical( S2$r.squared, S1$r.squared)
[1] FALSE #HUH?
Какво стана? Нека да използваме определението от статията в wikipedia за коефициента на определяне (или обяснения процент на отклонение, ако желаете):
1 - ((sum( residuals(lm1_)^2))/( sum( (y_ -mean(y_))^2)))
[1] 0.8575376
1 - ((sum( residuals(lm2_)^2))/( sum( (y_ -mean(y_))^2)))
[1] 0.8575376
identical(1 - ((sum( residuals(lm1_)^2))/(sum( (y_ -mean(y_))^2))), S1$r.squared)
[1] TRUE
identical(1 - ((sum( residuals(lm2_)^2))/(sum( (y_ -mean(y_))^2))), S1$r.squared)
[1] TRUE
identical(1 - ((sum( residuals(lm2_)^2))/(sum( (y_ -mean(y_))^2))), S2$r.squared)
[1] FALSE
S2$r.squared
[1] 0.9700402 # ???
Сега очевидно целият трик се обърква от функцията summary.lm(). R^2 се изчислява като mss/ (rss+ mss), където rss е остатъчната сума на квадратите sum(residuals(lmX_)^2), а mss е средната сума на квадратите (побиране), която зависи от наличието на прихващане или не (c-p от summary.lm()):
mss <- if (attr(z$terms, "intercept"))
sum((f - mean(f))^2)
else sum(f^2)
или в нашия случай:
sum(scale(fitted(lm2_),scale=F)^2) /
(sum(scale(fitted(lm2_),scale=F)^2)+sum(residuals(lm2_)^2))
[1] 0.8575376
sum(fitted(lm2_)^2) / (sum(fitted(lm2_)^2) + sum(residuals(lm2_)^2))
[1] 0.9700402
Така че ясно R не идентифицира матрицата на втория модел като имаща пресечка. Това обаче е сериозен проблем, ако, както в моя случай, някой мести моделни матрици наоколо. Знам, че очевидното решение ще бъде да се вземе матрицата на модела без прихващането, но това не облекчава факта, че това е дизайнерски избор, който е малко предразположен към счупване. Няма ли по-общата дефиниция да ни спести този проблем само с малка корекция? Има ли някакво друго очевидно решение, което мога да използвам, ако не се натъкна на тези проблеми?
Този проблем става изключително сериозен, ако човек не прихваща и използва факторни променливи с взаимодействия. На практика R може да се опита да третира едно от нивата като прихващане. Тогава цялата ситуация е просто трагична:
library(MASS)
lm_oats <- lm( Y ~ -1+ V*N , oats)
S_oats <- summary(lm_oats)
S_oats$r.squared
[1] 0.9640413 # :*D
1- ( sum( residuals(lm_oats)^2) / sum( ( oats$Y- mean(oats$Y))^2 ) )
[1] 0.4256653 # :*(
#For the record:
sum((fitted(lm_oats))^2 )/(sum( residuals(lm_oats)^2) +sum((fitted(lm_oats))^2))
[1] 0.9640413
sum((scale(scale=F,fitted(lm_oats)))^2 ) /( sum( residuals(lm_oats)^2)
+sum((scale(scale=F,fitted(lm_oats)))^2 ))
[1] 0.4256653
Това черта на дизайна на R, която по някакъв начин може да бъде заобиколена? (Оценявам, че въпросът ми беше доста отворен, благодаря на всички, които стигнаха дотук!)
Използвам R версия 3.0.1, Matrix_1.0-12, lattice_0.20-15, MASS_7.3-26. (32-битов Linux)
