怪しい研究成果の発表について（その2）

わざわざ言う必要もないのだけど，件の研究発表は，以下のようなくだらない解析結果だったのですよ。

> set.seed(1)
> d = data.frame(y = rnorm(20), x = letters[1:20])
> d
             y x
1  -0.62645381 a
2   0.18364332 b
    :
20  0.59390132 t

つまり，データの発生源を独立変数（説明変数）にしてしまった。

分析者は，1個の独立変数を使ったつもりだけど，実際は19個のダミー変数を使ってしまった。20 個のデータポイントを予測するのに，19個の変数を使うと，100%の予測が出来てしまう。それ以外の独立変数を使おうが使うまいが，100%に変わりはない。

> ans = lm(y ~ x, d)
> summary(ans) # 解は求まるが，何の意味もない。独立変数と従属変数は 1:1 の関係値にあるのだから，独立変数がある値を取るときは，対応する従属変数を予測値として返すと言うだけの話。

Call:
lm(formula = y ~ x, data = d)

Residuals:
ALL 20 residuals are 0: no residual degrees of freedom!

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.626454         NA      NA       NA
xb           0.810097         NA      NA       NA
xc          -0.209175         NA      NA       NA
    :
xt           1.220355         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:      1,    Adjusted R-squared:    NaN
F-statistic:   NaN on 19 and 0 DF,  p-value: NA

> cbind(d$y, predict(ans))
          [,1]        [,2]
1  -0.62645381 -0.62645381
2   0.18364332  0.18364332
3  -0.83562861 -0.83562861
    :
20  0.59390132  0.59390132

> all.equal(d$y,unname(predict(ans))) # 予測値は実測値と完全に一致する!!!
[1] TRUE

＝＝＝＝＝＝＝＝

R など，Factor（カテゴリー） 変数を自動的に処理してくれる統計処理システムだからこそ起こる悲劇。

もし，data.frame の x が自動的に Factor にされないとすると，d$x は

> d$x = as.integer(d$x)

として扱われたのと同じになる。

これだと，

> ans2 = lm(y ~ x, d)
> summary(ans2)

Call:
lm(formula = y ~ x, data = d)

Residuals:
    Min      1Q  Median      3Q     Max
-2.4808 -0.5168  0.1875  0.4833  1.5450

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03609    0.43158  -0.084    0.934
x            0.02158    0.03603   0.599    0.557

Residual standard error: 0.9291 on 18 degrees of freedom
Multiple R-squared:  0.01955,    Adjusted R-squared:  -0.03492
F-statistic: 0.3589 on 1 and 18 DF,  p-value: 0.5566

となるので，逆に疑わしい結果が得られる。

つまり，本来ダミー変数 a，b，c，... であるべきものが 1, 2, 3,  ... として使われるのだ。

この違いを知らないと，地獄に落ちる。

おお，神よ!!（^_^）