「統計学関連なんでもあり」の過去ログ--- 042

No.08562　R について　　【青木繁伸】　2008/12/08(Mon) 21:38

R は，フリーソフトということで，プログラムの間違いなどについても，一ユーザーが異論を唱えるということも自由に行えるというのが「うり」のはず。まあ，誰も彼もが自分の価値観で説を唱え始めると収集がつかないということはあるかも知れない。異論を唱えても，「私たちはそうは思わない」で片付けられることも多いようだ（実際そのようなあしらいを受けたこともあるが）。
先日，欠損値を含むデータセットでスピアマンの順位相関係数を求めるとき，欠損値を含まないデータだけを使って相関係数行列を求める際に，バグがあることに気づいた。まあ，先日の件もあるので，自分が注意すれば済むだけで，他人がバグの影響を受けようが知ったことではないという気がするので，CRAN に報告するのはやめにする。
ただ，日本語を解する我が同胞には，注意換気したいので，ここに記事をまとめておきたいと思う次第。

問題は，R 2.8.0 からの cor 関数を使って，スピアマンの順位相関係数を求める場合に限定される模様。しかも，データセットの内で，欠損値を含むケースを除外した（つまり，完全なデータセットについて計算する）場合についての問題のようだ。
計算例を示すような目的で，R に用意されているデータセットのうち，airquality というデータセットがある。このデータセットは，大気汚染に関する測定値の毎日のデータである。必然的に欠損値を多く含むデータセットである。
さて，一日の各種測定値に一つたりとも欠損値のないデータについて各測定値間の順位相関係数行列を求めようとして，バグの存在が発覚した。
欠損値が含まれない月日データは，na.omit(データセット）で得られる。cor 関数で，欠損値を含まないデータだけを使って相関係数行列を求めるには，use="complete.obs" を指定すればよいと，オンラインヘルプは語っている。じゃあ，二通りのやりかたで計算してみようじゃないか。。。。

> cor(na.omit(airquality[1:4]), method="spearman")
             Ozone     Solar.R        Wind       Temp
Ozone    1.0000000  0.34818647 -0.60513642  0.7729319
Solar.R  0.3481865  1.00000000 -0.06169636  0.2095369
Wind    -0.6051364 -0.06169636  1.00000000 -0.4993228
Temp     0.7729319  0.20953692 -0.49932278  1.0000000
> cor(airquality[1:4], use="complete.obs", method="spearman")
             Ozone     Solar.R        Wind       Temp
Ozone    1.0000000  0.35846998 -0.60797155  0.7807127
Solar.R  0.3584700  1.00000000 -0.07022037  0.2157419
Wind    -0.6079715 -0.07022037  1.00000000 -0.5034014
Temp     0.7807127  0.21574193 -0.50340142  1.0000000

同じ答えになるはずが違う。どちらが間違えているのか。
結論は，
cor(airquality[1:4], use="complete.obs", method="spearman")
が，間違えている。つまり，cor 関数の定義が間違えている。R 2.8.0 から，use 引数が拡張された。そのときに，対処を誤った（実は，cor 関数は，しばらく前までは，信じられない結果を生み出す実装をしていた。「こんな風にするしかないんじゃないの？」なんて，とんでもない居直りをしていたんだ(^_^;)）やっと，直ったと思ったら，また，前のようなとんでもない実装をしてしまったんだね。
取りあえず，スピアマンの順位相関係数を求めたいときで，SPSS で言うならリストワイズ欠損値除去による場合には，バグが修正されるまでは（いつの事やら），
cor(na.omit(airquality[1:4]), method="spearman")
のようにした方が宜しかろうということです。