Hatena::ブログ(Diary)

アイアナ:データ分析や人工知能(AI)などの技術雑記

2009 | 11 |
2010 | 02 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2011 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2012 | 01 | 03 | 04 | 05 | 10 | 12 |
2013 | 01 | 02 | 04 |
2014 | 03 | 12 |
2016 | 03 |
2017 | 10 |

2010-09-25

メモリの確認

RはPCのメモリ上にデータ(オブジェクト)を一時的に記録しています。

現在どの程度メモリを使用しているかをチェックするには、

memory.size(T)

をタイプし、利用可能な最大メモリを調べるためには、

memory.limit()

をタイプします。

つい最近私は64bitマシンを買い替えてRも64bit版をインストールしまして、memory.limit()をタイプ実行すると、

3959

と返ってきました。

以前のPCでは1535だったので、64bitになったんだなぁと実感。


いや、まてよ以前のPCにも4Gのメモリを積んでたよ?

win32だと3.2Gまでしかメモリを使えないとしても、今までも半分くらいまでしか使ってなかったのかな?


今回のPCはメモリが8Gまでいけるみたいなので、時間があるときに増設するつもりです。

64bitのRでどこまで大規模データを扱えるのか。。。人柱になってきます!

2010-09-24

大規模データに対する統計手法

大規模データと一口に言っても、実は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というのは次の図のような感じ。

f:id:isseing333:20100924170411j:image


この図は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)

2010-09-21

多重比較

臨床データの解析、特に経時データの解析をする際に気を付けなければいけないのが検定の多重性です。

検定のαエラーは「ほんとうは差が無い(帰無仮説が真)データなのに間違って差があるといってしまう確率」です。

5%の有意水準(α=0.05)の下で帰無仮説が真のデータを100個検定すると、5個は有意になります。

例1.例えば5時点の経時データで、全ての時点で差がないとします。

その5時点全てを対比較すると、5C2の10通りあります。

この10通りを検定すると期待的に0.5個は有意になり、αエラーに換算すると1-(0.95)^10=0.4になります(つまり40%の確率で第一種の過誤が起こる)。

例2.さらに、変数(検査値)が10個あれば、10検査値×経時10通り=100(回)検定することになりますので、期待的に5個は有意になります。

αエラーは1-0.95^100=0.99、つまり99%の確率で...(ry

このような現象をαエラーの増大などと言いますが、それを解消するのが多重性の調整です。


多重性の調整方法は大きく分けて2種類あります。

  • FWERの調整(Family Wise Error Rate)
  • FDRの調整(False Discovery Rate)

FWERはさっきの例でいうと、例1の5時点の対比較のなかでαエラーを名義水準(5%)に保つということ。

方法としては以下のようなものがあります。

    • Bonferroni法
      • P値を検定した回数でかけて、有意水準(5%)と比較する、簡単に手計算でできるけど他の方法より有意になりにくい
    • Tukey法
      • 対比較のときに利用する方法
    • Dunnett法
      • 対照群と他の群の比較に利用、経時データなどで第1時点とその後を比較したいときなど
    • Williams法
      • 対照群と他の群の比較だが、単調変化しているという仮定が入る。第1時点と比較して、第○時点以降は全て差がある場合なと

FDRはBenjaminiとHockbergが1995年に提唱した新しい概念で、主に遺伝子の分野で利用されている方法です。

定義は「有意になった仮説のうち、間違って検出された仮説の割合の期待値が有意水準以下になるようにする」というものです。

例えば10,000回FDR=0.05水準で検定して100個有意になったとしたら、期待的に5個は間違っているのです(これをFWER=0.05で調整しようとすると、有意になる変数はほとんど無いでしょう)。


Rで解析する場合はFWERではmultcompというパッケージがあります。

以下解析例です。

breaksという変数をtensionで説明したいので、まず箱ひげ図を確認。

f:id:isseing333:20100921143422j:image

水準Lを対照として、Dunnett法で調整したP値は次のようになる。

M - L == 0   -10.00       3.96  -2.525 0.027536 *  
H - L == 0   -14.72       3.96  -3.718 0.000978 ***

この場合は、調整しても両方とも有意に下がっているようです。

対比較の場合はTukeyにすれば良くて、結果は次のようになります。

M - L == 0  -10.000      3.960  -2.525  0.03847 * 
H - L == 0  -14.722      3.960  -3.718  0.00144 **
H - M == 0   -4.722      3.960  -1.192  0.46311   

HとMはそんなに統計的には差があるといえないみたい。


Rプログラムはこちら↓

library("multcomp")

#breaks変数をtensionで説明したい
plot(warpbreaks[,c(3,1)])

#モデルの指定
amod <- aov(breaks ~ tension, data = warpbreaks)

#多重性の調整をしたP値の算出
amod.du <- glht(amod, linfct = mcp(tension = "Dunnett"))
summary(amod.du)

amod.tu <- glht(amod, linfct = mcp(tension = "Tukey"))
summary(amod.tu)

2010-09-17

混合効果モデルと固定効果モデルの比較

前回の続きです。

変量効果モデルで推定した場合と、固定効果モデルで交互作用として推定した場合の違いを絵でみるとこのようになります。

成長曲線のデータで、8〜14才を追跡し身長が伸びた分を記録しているそうです。

以下のグラフは9人分の図。

  • 変量効果

f:id:isseing333:20100917182804j:image


  • 固定効果

f:id:isseing333:20100917182803j:image


まぁ、あんまり変わらないね。。。

変量効果の方が欠測値を扱いやすいとかっていう利点があった(固定効果でも推定はできるけど)。

あまりに欠測が多かったりカテゴリが多いとMCMCで推定する必要が出てきます。

結果はあまり変わらないけど、変量効果モデルを使うと「お、統計を知ってるな」って思われるかも。

プログラムはこちら↓


#------変量効果と固定効果の比較
#---変量効果
library(nlme)
formula(Orthodont)

#交互作用モデルではFactorの値が使われるので文字にしておく
Orthodont$Subject <- as.character(Orthodont$Subject)

fm1      <- lme(distance ~ age, data = Orthodont)	#subject毎にageと切片の推定値を算出
IntMix   <- fm1$coefficients$fixed[1] + fm1$coefficients$random$Subject[, 1]		#混合効果モデルでの切片
SlopeMix <- fm1$coefficients$fixed[2] + fm1$coefficients$random$Subject[, 2]		#混合効果モデルでの傾き

MixedPlot <- function(Group){
	plot(Orthodont[Orthodont[ ,3] == Group, c(2, 1)], main=Group)
	abline(a=IntMix[names(IntMix) == Group], b=SlopeMix[names(SlopeMix) == Group])
}

par(mfrow=c(3, 3), pty="m")
MixedPlot("M01")
MixedPlot("M02")
MixedPlot("M03")
MixedPlot("M04")
MixedPlot("M05")
MixedPlot("M06")
MixedPlot("M07")
MixedPlot("M08")
MixedPlot("M09")


#---固定効果
fm2      <- lm(distance ~ age + age * Subject, data = Orthodont)	#固定効果で交互作用を指定したモデル
IntFix   <- fm2$coefficients[1] + fm2$coefficients[3:28]
SlopeFix <- fm2$coefficients[2] + fm2$coefficients[29:54]

FixedPlot <- function(Group){
	plot(Orthodont[Orthodont[ ,3] == Group, c(2, 1)], main=Group)
	#---交互作用モデルでは推定パラメータ名に変数名Subjectが付いているので名前をあわせる
	SGroup <- paste("Subject", Group, sep="")
	SAGroup <- paste("age:Subject", Group, sep="")
	abline(a=IntFix[names(IntFix) == SGroup], b=SlopeFix[names(SlopeFix) == SAGroup])
}

par(mfrow=c(3, 3), pty="m")

#M01を基底にして交互作用を推定しているので、M01だけはそのまま描く
plot(Orthodont[Orthodont[ ,3] == "M01", c(2, 1)], main="M01")
SGroup <- paste("Subject", "M01", sep="")
SAGroup <- paste("age:Subject", "M01", sep="")
abline(a=IntFix[names(IntFix) == SGroup], b=SlopeFix[names(SlopeFix) == SAGroup])

FixedPlot("M02")
FixedPlot("M03")
FixedPlot("M04")
FixedPlot("M05")
FixedPlot("M06")
FixedPlot("M07")
FixedPlot("M08")
FixedPlot("M09")