大規模データに対する統計手法
大規模データと一口に言っても、実は2種類のデータがあります。
1.対象者(サンプル数)が大規模
2.変数(項目)が大規模
1.は健診、コホート等のデータなどが相当しますし、2.は遺伝子等のデータが相当します(場合によってはQOL等も)。
サンプル数をn、変数をmとすると前者はn>>m、後者はm>>nなので全く違う性質を持っています。
それぞれのデータに対しての統計的な注意点をまとめてみます。
- サンプル数が大規模なデータ
サンプル数が数千~数億といったデータでは、α水準が5%等で検定を行うこと自体が無意味なことが多いです。
なぜかというと、サンプル数が多いと必然的に平均値の標準誤差が小さくなるので有意になりやすくなるから。
例えば2群間で身長の差を検定するとすると、各群10人の場合は1cmの差は検出できませんが、各群10,000人居れば差が1cmしかなくても有意になります。
もう少し厳密にいえば身長の分散が10cmだとすると、差が1cmのとき各群10人のデータでt検定して有意になる確率(検出力)はわずか4%しかないのに対して、各群10,000人の場合は99%となります。
各群100,000人居る場合は、1mmの差も60%の確率で検出します。
検定は事前にサンプルサイズ設計を計画して始めて意味がでてくるのだと思います。
サンプルサイズ設計とは「分散がこれくらいで差がこれくらいだと、80%の確率で検出する(有意になる)ためには何例以上必要」というものです。
大体の場合では数百例居れば検出力80%を満たしますので、数万人データで検定をすることはナンセンスに思えます。
大規模サンプルデータの場合は、差の絶対値と標準偏差(標準誤差ではない)自体が重要なのではないでしょうか。
検定をするよりも、疫学情報として分布(平均と分散等)をチェックする方が重要だと思います。
あと相関係数の検定も同様のことが言えます。
相関係数が有意になるかどうかも、サンプル数に依存します。
例えば、相関係数が0.1しかなくても、サンプル数が400以上あれば有意になります(相関係数の有意水準の表が教科書やネット上にあります)。
相関係数が0.1というのは次の図のような感じ。
この図は400例で相関係数は0.139、P値は0.0054なのですが、これで有意に相関がありますって言われても、ねぇ?
相関係数の場合も有意云々よりも絶対値が重要で、厳密な基準は無いものの0.7くらいあれば相関が高いと言えます。
相関係数0.7というのは単回帰したときのR2乗が0.49ですので、ばらつきの半分を1つの変数で説明できていることになります。
- 変数が大規模なデータ
遺伝子データなどは変数が数万もあるのに、サンプルが数十から数百しかないということが多く見受けられます。
この場合に統計的に気をつけなくてはならないのが第一種の過誤(疑陽性)です。
普通の検定のようにαエラーが5%水準で数万の変数を検定すると、全ての変数で差が無かったとしても5%は有意になります。
つまり10,000遺伝子を普通に検定すると、全く意味の無いデータであっても500個は有意になるのです。
これを解消するために、FDRという概念が生まれました。
FDRは「有意になった変数のうち、間違っている変数の割合をある水準未満に抑える」というものです。
例えば10,000遺伝子をFDRが5%水準で検定して500個有意になった場合は、そのうち25個が期待的に間違っていることになります。
単純に考えると、この結果はα水準では25/10,000*100=0.25%に相当します(5%の20分の1)ので、FDRで調整すると非常に厳しい水準であることがわかります。
このような超多変量データを使って予測を行うときには、他にもover fitting(過適合)という問題があり、これを解消するためにはcross validation(クロスバリデーション)を行わなくてはなりません。
クロスバリデーション(もしくはブートストラップ)を行うと、統計的には過適合が防げるというシミュレーション結果が出ています(Efron先生の論文)。
しかしクロスバリデーションも行い方によっては、またバイアスが入ることも確認されています(Simon先生の論文)。
遺伝子データで有意な結果が出ると新発見かのように取りざたされることが多いですが、統計的に正しい方法を使っているのかを吟味しなくてはいけませんよね。
今回の内容に関するRプログラムはこちら↓
n <- 10 a <- data.frame(group = rep(c(0, 1), c(n, n)), Height = c((rnorm(n)*10 + 170) ,(rnorm(n)*10 + 171))) t.test(Height ~ group, data=a) n <- 10000 a <- data.frame(group = rep(c(0, 1), c(n, n)), Height = c((rnorm(n)*10 + 170) ,(rnorm(n)*10 + 171))) t.test(Height ~ group, data=a) power.t.test(n=10, delta=1, sd=10) power.t.test(n=10000, delta=1, sd=10) power.t.test(n=100000, delta=0.1, sd=10) library(MASS) set.seed(2) b <- mvrnorm(400, mu=c(0, 0), Sigma=matrix(c(1, 0.1, 0.1, 1), 2, 2)) cor(b) b <- data.frame(V1 = b[, 1], V2 = b[, 2]) plot(b) summary(lm(V1 ~ V2, data=b)) cor.test(b$V1, b$V2)