怪しい研究成果の発表について（その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,  ... として使われるのだ。

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

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

怪しい研究成果の発表について

怪しいんですよと解説してくれているサイトでの表示が，奥歯にものが挟まったような気がするので...

まず，データをダウンロードして実行して見ると，件のサイトには明示されていない「警告メッセージ」が現れる。

>  d.glm<-glm(Result~.,d,family=binomial)
 警告メッセージ:
1:  glm.fit: アルゴリズムは収束しませんでした  
2:  glm.fit: 数値的に 0 か 1 である確率が生じました  

つまり，「以下に示される結果は『信頼できませんよ』ということですかね

それを無視して，分析を進めると

> table(d$Result,round(predict(d.glm,newdata=d[,-3],type='response'),0))
   
      0   1
  0 251   0
  1   0 240
 警告メッセージ:
 predict.lm(object, newdata, se.fit, scale = 1, type = ifelse(type ==  で:
  prediction from a rank-deficient fit may be misleading

またしても，警告メッセージが出ていますけど。（表示される結果は，件のサイトと同じ）

以下，件のサイとで行われたとおりの分析を続けると，出てくる結果はそのサイトに示されているのと全く同じ。

違う所（というか，件のサイとの結果の表示では省略されているだけかも知れないのだけれど），「警告メッセージ」というのが出てくるのよ。

こんなのが，毎回出てくると「心穏やかにならなくなるよね」

> d.glm.true<-glm(Result~.,d[,-c(1,2)],family=binomial)
> table(d$Result,round(predict(d.glm.true,newdata=d[,-c(1,2,3)],type='response'),0))
   
      0   1
  0 238  13
  1  14 226
 警告メッセージ:
 predict.lm(object, newdata, se.fit, scale = 1, type = ifelse(type ==  で:
  prediction from a rank-deficient fit may be misleading
> library(e1071)
> d.svm<-svm(Result~.,d,kernel='linear')
> table(d$Result,predict(d.svm,newdata=d[,-3]))
   
      0   1
  0 251   0
  1   1 239
> d.svm.true<-svm(Result~.,d[,-c(1,2)],kernel='linear')
> table(d$Result,predict(d.svm.true,newdata=d[,-c(1,2,3)]))
   
      0   1
  0 238  13
  1  15 225
> nrow(d[,1:2])
[1] 491
> nrow(unique(d[,1:2]))
[1] 485
> d.glm.name<-glm(Result~.,d[,1:3],family=binomial)
 警告メッセージ:
1:  glm.fit: アルゴリズムは収束しませんでした  
2:  glm.fit: 数値的に 0 か 1 である確率が生じました  
> table(d$Result,round(predict(d.glm.name,newdata=d[,1:2],type='response'),0))
   
      0   1
  0 247   4
  1   2 238
 警告メッセージ:
 predict.lm(object, newdata, se.fit, scale = 1, type = ifelse(type ==  で:
  prediction from a rank-deficient fit may be misleading
> summary(d.glm)

Call:
glm(formula = Result ~ ., family = binomial, data = d)

Deviance Residuals:
       Min          1Q      Median          3Q         Max  
-7.228e-06  -2.409e-06  -2.110e-08   2.409e-06   8.367e-06  

Coefficients: (47 not defined because of singularities)
                                Estimate Std. Error z value Pr(>|z|)
(Intercept)                    1.184e+02  4.986e+06       0        1
Player1A.Kuznetsov            -3.999e+01  1.553e+06       0        1
Player1A.Man0rino             -2.280e+01  1.667e+06       0        1
Player1A.Ramos                -1.334e+01  2.653e+06       0        1
Player1A.Seppi                 2.574e+01  1.920e+06       0        1
Player1A.Ungur                 6.009e+01  1.841e+06       0        1
Player1Adrian Man0rino        -2.082e+01  2.717e+06       0        1
　: あまりにも多いので，途中省略

BPC.2                         -6.438e+00  9.798e+04       0        1
BPW.2                          1.426e-01  4.375e+04       0        1
NPA.2                         -1.442e-01  6.894e+04       0        1
NPW.2                          4.440e-01  7.039e+04       0        1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 6.8042e+02  on 490  degrees of freedom
Residual deviance: 3.9260e-09  on  74  degrees of freedom
AIC: 834

Number of Fisher Scoring iterations: 25

この結果見ただけで，普通の人は，「げんなりする」はず。

原論文の解析者は，この結果を，見ていないのかな（みていないんだろうな）。

ビッグダーどころか，統計学，統計データ，多変量解析について，なんの知識もない人なんだなあというのが，正直な感想（卒論レベルの評価でも，不合格ですよ）。筆者は工学分野が専門なようで，当該のような医学・薬学分野の常識が欠除していたんだろうなあと思います。

このような人を特任研究員？として任用した，京都大学の然るべき部署の責任が問われるべきでしょうね。

一度は公表して，すぐに取り下げたということであるが，責任ある京都大学の然るべき部署として，なぜ取り下げたか，どういう根拠に基づいて取り下げるを得なかったのかをしっかりと説明しないと，ダメですよ～～～～～。って。嘆かわしい。