3

私は R の関数をいじっていましたが、決定係数 R^2lm()の計算で奇妙なことに気付きました。おもちゃの回帰問題があるとします。

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?

どうしたの?ウィキペディアの記事の定義を決定係数(または必要に応じて説明された分散のパーセンテージ) に使用してみましょう。

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 # ???

どうやら全体のトリックは function によって当惑するようsummary.lm()です。R^2 は次のように計算されますmss/ (rss+ mss)rssは残差平方和でsum(residuals(lmX_)^2)ありmss、切片の存在に依存する平均平方和 (フィット) です (summary.lm() からの cp):

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 は 2 番目のモデル行列が切片を持つとは識別しませんでした。私の場合のように、誰かがモデルマトリックスを動かしている場合、それは重大な問題です。切片なしでモデル マトリックスを取得することが明らかな解決策であることはわかっていますが、これは、これが少し壊れやすい設計上の選択であるという事実を緩和するものではありません。より一般的な定義では、わずかな調整だけでこの問題を解決できるのではないでしょうか? これらの問題にぶつからないように使用できる他の明白な修正はありますか?

この問題は、インターセプトがなく、交互作用のある因子変数を使用している場合、非常に深刻になります。実際には、R はレベルの 1 つをインターセプトとして処理しようとする場合があります。その後、全体の状況は悲劇的です。

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 を使用しています。(Linux 32 ビット)

4

1 に答える 1