Transcript p値 - 京都大学
P値から考えるゲノム疫学解析 日本人類遺伝学会 2014/11/20 京都大学 医学研究科 統計遺伝学分野 山田 亮 • SNPを用いたGWASなどのスタディに関 わっている・かかわろうとしている大学院 修士の学生さん、くらいを主なターゲット とします • ご参加のみなさんには、 – SNPって何？ 検定とp値？ そもそもたくさ ん検定をするって何？ という方や – あそこの分割表の検定は、厳密には、●● しないといけないのでは • という方が混ざっている・・・と思います が、せっかくですので、何かしら得るもの がありますように 今日のスライドの大部分のデータ処理はフリーソフトRで実 施しています。そのコードを含めた資料は、分野facebookか ら入手可能です。聴講しながら、叩いていただいても結構で す。 2x2分割表で確認する、p値、オッズ比 とその信頼区間 2x2分割表で確認する、p値、オッズ比 とその信頼区間 あ り コントロール 20 28 ケース 48 な し 27 47 25 53 52 100 検定する・推定する 検定する・推定する • カイ二乗統計量 0.6825 • 自由度 1 • p-値 0.4087 検定 • p値 • 棄却するべきかどうかの情報 • p値が小さいほど、無関係という仮説（帰無仮 説）を棄却する • その心は、「帰無仮説が正しいと仮定したとき に、このテーブルと同じかそれよりもっと偏っ たテーブルを観察する確率はpくらい小さいか ら」 あり なし 推定 • オッズ比 0.6614 • 95%信頼区間 0.3 ～1.458 コ 20 27 47 ケ 28 25 53 48 52 100 推定 • オッズ比を計算する。 • オッズ比は、相対危険度の推定値 • オッズ比～相対危険度が１であるときが、「無関 係」に相当する • 推定値は – 点推定値 – 信頼区間 – 「だいたいこのくらいの範囲」の中に、「無関係の相対 危険度～１」が入っているかどうかで、およその帰無 仮説を棄却するかどうかのめども立つ OR 95% CI p値と CIが1をまたぐか の関係 緩く正しく 微妙に違う 1.0 p値 p値はなんのため？ • • • • • • 小さいp値も大きいp値も同じくらい出やすい 0.5が出やすいわけではない 1が出やすいわけではない 0.5以下のp値が出る確率が0.5 0.01以下のp値が出る確率が0.01 これがわかりやすいから、p値を出して考える • 一様分布 p値はなんのため？ • • • • • • 小さいp値も大きいp値も同じくらい出やすい 0.5が出やすいわけではない 1が出やすいわけではない 0.5以下のp値が出る確率が0.5 0.01以下のp値が出る確率が0.01 これがわかりやすいから、p値を出して考える • 一様分布 一様分布のヒストグラム 0 p値 1 順番に並べてプロットすると 直線になる １番 小さい順 １万番 小さい順に並べて 対角線を描くのは QQプロット • 対数を取っても OK • p値だけではなく • カイ二乗値でもＯ Ｋ Lanktree M B et al. Stroke. Copyright © American Heart Association, Inc. All rights 2010;41:825-832 QQプロット 理論に合っているかを確認する 観 測 値 0 1 期待値・理論値 p値の基礎、終了 SNV解析の基本、2x3分割表 M M コントロール 10 5 ケース 15 M m 10 20 30 m m 10 30 40 30 55 85 コントロールに比べて、ケースが： 0.5、2、3倍 M M コントロール 10 5 ケース 15 M m 10 20 30 m m 10 30 40 30 55 85 さっそく「関連検定」しよう さっそく「関連検定」しよう • 面倒くささの元は さっそく「関連検定」しよう • 面倒くささの元は –いくつも検定法があること 2x3分割表の検定法 ２つ • トレンド検定とロジスティック回帰検定 トレンド検定する • p 値 0.007476 M M コントロール 10 5 ケース 15 M m 10 20 30 m m 10 30 40 30 55 85 オッズ比は？ • p 値 0.007476, (10x20)/(10x5) = 4 M M コントロール 10 5 ケース 15 M m 10 20 30 m m 10 30 40 30 55 85 オッズ比その２ • p 値 0.007476, (10x30)/(10x5) = 6 M M コントロール 10 5 ケース 15 M m 10 20 30 m m 10 30 40 30 55 85 トレンド検定する • p 値 0.007476, • オッズ比 –4 (Mm vs. MM) –6 (mm vs. MM) トレンド検定のオッズ比… • トレンド検定をするというのは、線形回帰をし ているのと同じです。 線形回帰では、ジェノタ イプごとのケースの割合が直線に乗るように 推定します • この『トレンド検定をするという気持ち』に照ら すと、このオッズ比はちょっと違います • どう違う？ トレンド検定という線形回帰 フェノタイプ の現れる確 率 1.0 0.0 ジェノタイプ 線形回帰は「傾き」が大事 フェノタイプ の現れる確 率 1-p2 1.0 1-p1 0.0 p1 p2 ジェノタイプ (p2/1-p2) (p1/1-p1) 線形回帰は「傾き」が大事 フェノタイプ の現れる確 率 1-p2 1.0 (p2/1-p2) (p1/1-p1) 1-p1 0.0 p1 p2 オッズ比 2.1 4 (Mm vs. MM) 5.1 6 (mm vs. MM) ジェノタイプ 線形回帰は「傾き」が大事 フェノタイプ の現れる確 率 1-p2 1.0 (p2/1-p2) (p1/1-p1) 1-p1 0.0 p1 p2 オッズ比 2.1 4 (Mm vs. MM) 5.1 6 (mm vs. MM) ジェノタイプ ロジスティック回帰する • p 値 0.007476 (トレンド検定) • p 値 0.009525 (ロジスティック回帰検定) • 大差ない ロジスティック回帰する フェノタイプ の現れる確 率 1.0 0.0 ジェノタイプ ロジスティック回帰する フェノタイプ の現れる確 率 1-p2 1.0 1-p1 0.0 p1 p2 オッズ比 2.1 → 2.3 (Mm vs. MM) 5.1ジェノタイプ → 5.2 (mm vs. MM) ロジスティック回帰する フェノタイプ の現れる確 率 1-p2 1.0 1-p1 0.0 p1 p2 オッズ比 2.1 → 2.3 5.1 → 5.2 = 2.3 x 2.3 ジェノタイプ トレンド検定とロジスティック回帰 トレンド検定とロジスティック回帰 直線と曲線の 違いはあるけれど 大差ない ロジスティック log10(p) 大差がない トレンド log10(p) 違いは • トレンド検定は – – – – – – 相加モデルに相当する 計算が簡単 正確検定もできる(低アレル頻度・少サンプル数の場合の対処法がある) 個人別のデータが不要 「線形回帰」に基づくオッズ比が結果に付いていないことが多い 共変量を組み込めない(線形回帰にすればよいが、それならロジスティック回 帰をすればよい) • ロジスティック回帰検定は – 相乗モデルに相当する – 計算が面倒ではある(けれど計算機がやってくれるので問題はないが、検定 個数が大量になるとそれなりに影響してくる) – オッズ比が結果についてくることが多い(ただし、係数は対数で返ってくること が通例) – 年齢・性別などの共変量を組み込みやすい 遺伝形式 相加 相乗 優性 劣性 • 面倒くささの元は –いくつも検定法があること 優性 優性 劣性 優性・劣性形式の検定 するのか、しないか • することのメリットとデメリット – 優性・劣性形式に照らした結果がわかる – 1つの2x3表に複数の検定をすると、複数のp値が 得られる。複数のp値が出たら、そのp値はには 「補正」をしないといけない • しないことのメリットとデメリット – 1つのp値しか出ていなければ、複数のp値の補 正について悩む必要はない – 優性・劣性形式に照らしての判断になっていない 優性・劣性形式の検定をしなかったら • 本当は「優性、または、劣性」な影響があると すると、偽陰性が増える。 • では、「優性形式が真」のときに、「相加モデ ルだけ」で検定すると、どれくらい偽陰性にな るのかを見てみることにする。 優性モデルlog10(p) アレル頻度 0.3 RR(優性) RR=1.2 1000人 vs. 1000人 相加モデルlog10(p) 優性モデルlog10(p) アレル頻度 0.3 30% RR(優性) RR=1.2 1000人 vs. 1000人 相加モデルlog10(p) 優性モデルlog10(p) アレル頻度 0.3 RR(優性) RR=1.2 30% 24% 相加モデルlog10(p) 1000人 vs. 1000人 優性モデルlog10(p) アレル頻度 0.3 偽陰性 RR(優性) RR=1.2 30% 24% 相加モデルlog10(p) 1000人 vs. 1000人 優性モデルlog10(p) 優性モデルで は拾えずに相 加モデルで 「たまたま拾 う」こともある アレル頻度 0.3 偽陰性 RR(優性) RR=1.2 1000人 vs. 1000人 相加モデルlog10(p) 優性モデルlog10(p) 相加 優性 両方を併せると 相加モデルlog10(p) パワーが上がる パワーが上がるのはよいことだ 相加・優性どちらも『あり』にしたら 偽陽性が1.8倍になった 真の優性座位 無関係の座位 パワーが上がると 偽陽性が増える いいことがあると 悪いこともある パワーを上げつつ、偽陽性を増やさない 0.01より小さい『新たなp値基準』を作る 1% 稼ぐ 失う 『新たなp値基準』を探す • この図があれば、できる • けれど、この図はない(すぐには手に入らない) 1% 素朴なマルチプルテスティング対策 素朴なマルチプルテスティング対策 黒、赤、緑、青、の比率がわかれば、 『新たなp値基準』はわかる 素朴なマルチプルテスティング対策 黒、赤、緑、青、の比率がわかれば、 『新たなp値基準』はわかる 黒、赤、緑、青、の比率がわかれば、 素朴なマルチプルテスティング対策 『新たなp値基準』はわかる 黒、赤、緑、青、の比率は 素朴なマルチプルテスティング対策 わかる 一様分布なら １辺の長さが１の正方形 ?の長さを求めなさい 0.99 ? ? １辺の長さが１の正方形 ?の長さを求めなさい 0.99 Sidak法 ? ? 細長い白枠長方形2個の面積が 0.01になるとき ? の長さはいくつ か？ 0.99より大 きくする ボンフェリニ法 ? ? 問題は、偏りがあること 優性・相加のp値には相関がある 分布がわからないので、どうするか • 分布がわからないままに、補正する • わからない分布を調べてから、それに基づい て補正する 分布がわからないままに、補正する 分布が違うのに、それでよいのか？ ボンフェロニ・Sidakを使うと 偽陽性が少なくなる パワーが弱くなる ～ストイックであれば大丈夫～ ～保守的であれば大丈夫～ 分布がわからなければ 分布を調べればよい 正確確率法 ランダマイゼーション・ パーミュテーション法 本当に知りたいこと • 本当は「関連がない」ときに • 相互に相関のある複数の検定を実施したとき に • 最も小さいp値は、どれくらい小さければ0.01 並みに珍しいか 本当に知りたいこと • 本当は「関連がない」ときに • 相互に相関のある複数の検定を実施したとき に • 最も小さいp値は、どれくらい小さければ0.01 並みに珍しいか • これは、ちょっと面倒なので、少し変えます 本当に知りたいこと • 本当は「関連がない」ときに • 相互に相関のある複数の検定を実施したとき に • 『今、観測された分割表』の周辺度数を満足 する場合のすべてを考慮して • 最も小さいp値は、どれくらい小さければ0.01 並みに珍しいか 『今、観測された分割表』の周辺度数 を満足する場合のすべてを考慮する • ２つの方法 – 本当に「すべての場合」を考慮する • 正確確率法 – 乱数を使って「一部の場合」を考慮して代用する • モンテカルロ・ランダマイゼーション法、パーミュテー ション法 正確確率法とランダマイゼーション法 の違い • 正確確率法 – 『正確』 – すべての場合を扱えるのは、自由度2くらいまで。 それは、2x3表が１個ある場合。 • ランダマイゼーション法 – 『推定値』 – 試行ごとに少し違う – 1000の場合をやれば、最小p値は0.001、10000 回やれば、最小p値は0.0001 GWAS基準の有意p値はとても小さい けれど、 それはどうするの？ • 正確確率法 – 『正確』 – すべての場合を扱えるのは、自由度2くらいまで。 それは、2x3表が１個ある場合。 • ランダマイゼーション法 – 『推定値』 – 試行ごとに少し違う – 1000の場合をやれば、最小p値は0.001、10000 回やれば、最小p値は0.0001 GWAS基準は「デフォルト推奨値」 正確確率法って？ 正確確率法って？ • サンプル数が少ないときにカイ二乗検定の代 わりに使う方法 • サンプル数が少ない、というより、分割表のセ ルの値が小さいとき… • 分割表のセルの値が小さいとき、というより、 セルの期待値が小さいとき… どうしてか？ • セルの期待値が小さめのときには、カイ二乗 検定のp値は『不正確』だから • 正確確率検定の方が『保守的』だから • 『保守的』なことは、『よいこと』だから たくさんのp値を正確法で得ると… 正確確率検定 カイ二乗検定 たくさんのp値を正確法で得ると… 正確検定は、「保守 的」なので、その結果 をたくさん集めると、 一様分布からは随分 ずれる 連鎖不平衡とp値 連鎖不平衡とマルチプルテスティング • 1つのSNP • 複数の遺伝的 モデル • 複数の検定 • 1つの遺伝子 • 複数のSNP • 個々のSNPに １つの検定 • 複数の検定 相互に独立ではない 複数の検定 • １つのSNP • 複数のSNP • 相加・優性・劣性 の３検定は、相互 に独立ではない • 相互に連鎖不平 衡にある • 個々のSNPの検 定は独立ではな い 相互に独立ではない 複数の検定 • １つのSNP • 複数のSNP • 相加・優性・劣性 の３検定は、相互 に独立ではない • 相互に連鎖不平 衡にある • 個々のSNPの検 定は独立ではな い マルチプルテスティング補正が必要 連鎖不平衡領域 SNPごとに相加検定 20個のp値 ケース・コントロールスタディを実施 有意水準0.01で関連ありとするには、どれくらい 小さいp値が適当？ クイズ、1-6のどれ？ 1. 0.01/20 = 0.0005 2. 0.000502 = 1-(1-0.01)^(1/2) 3. 4. 5. 6. 0.00045 0.00054 0.00064 0.00087 クイズ、1-6のどれ？ 1. 0.01/20 = 0.0005 ボンフェロニ 2. 0.000502 = 1-(1-0.01)^(1/2) Sidak 3. 4. 5. 6. 0.00045 0.00054 0.00064 0.00087 クイズ、1-6のどれ？ 1. 0.01/20 = 0.0005 ボンフェロニ 2. 0.000502 = 1-(1-0.01)^(1/2) Sidak 3. 4. 5. 6. 0.00045 0.00054 0.00064 0.00087 クイズ、1-6のどれ？ 1. 0.01/20 = 0.0005 ボンフェロニ 2. 0.0005023906 Sidak 3. 4. 5. 6. 0.00045 0.00054 0.00064 0.00087 ボンフェロニや Sidakより小さ いわけがない ３つの数字、３つのLD図 1. 0.01/20 = 0.0005 ボンフェロニ 2. 0.0005023906 Sidak 3. 4. 5. 6. 0.00045 0.00054 0.00064 0.00087 ３つの数字、３つのLD図 1. 0.01/20 = 0.0005 ボンフェロニ 2. 0.0005023906 Sidak 3. 4. 中5. 6. 0.00045 0.00054 0.00064 0.00087 強 弱 強 中 弱 GWAS基準は「デフォルト推奨値」 SNPの数を十倍の 一千万個に増やしたら？ • 百万個の場合と同じ領域をカバーしながら、 個数だけ増やしたのなら、マルチプルテスティ ングで目指すべきp値はそれほど下げなくて もよい • 百万個の場合にカバーしなかった場所をカ バーすることにしたのなら、マルチプルテス ティングで目指すべきp値は増やした分だけ 下げる(十分の１にする)ことになります 連鎖不平衡にあるマーカーで 代用する • • • • LDマッピングの原理そのもの SNP 1 : A / a の２アレル SNP 2 : B / b の２アレル ハプロタイプは４種類 – AB – Ab – aB – aB 2SNPの4ハプロタイプは 2x2 分割表 A a B b 0.78 0.02 0.02 0.18 0.2 0.8 0.8 0.2 1 2x2 分割表なら カイ二乗検定しよう カイ二乗値 = 0.81 = r2 A a B b 0.78 0.02 0.02 0.18 0.2 0.8 0.8 0.2 1 LD関係にある 2 SNPのカイ二乗値の 相関の良さとLDインデックス r2 カイ二乗値 の相関係数 2 LDのr LD関係にあるSNPで 代用したときのパワー (r2 = 0.81) 本体 代用 代用SNPの場合 LD関係にあるSNPによるパワー 真のリスクSNPの場合 代用SNPの場合 p値の高低、どちらが小さい？ 代用SNPのp値 の方が小さい 真のリスクSNPの場合 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 2 r =1 代用SNPの場合 p値の高低、どちらが「本物」？ 代用SNPのp値 の方が小さい 真のリスクSNPの場合 p値の高低、どちらが「本物」？ 1 代用のp値 が本物以下 の割合 0 0.1 2 LDのr 1 アレルで、ハプロタイプで 検定する M M コントロール 10 5 ケース 15 M m 10 20 30 m m 10 30 40 30 55 85 アレル本数で考える • 2x3表の相加モデル(トレンド)検定 • 2x2表を作って普通にカイ二乗検定(独立性検 定) MM Mm mm 10 10 10 5 20 30 15 30 40 M 30 30 60 m 30 80 110 30 55 85 60 110 170 MM Mm mm 10 10 10 5 20 30 15 30 40 M 30 30 60 m 30 80 110 30 55 85 60 110 170 20 29 49 25 26 51 24 18 42 22 18 50 6 3 9 6 3 9 50 50 100 23 26 49 50 50 100 クイズ ３つの分割表 相加と2x2とが同じ表が 2個ある。どれ？ 24 18 42 6 3 9 50 50 100 20 29 49 25 26 51 24 18 42 22 18 50 6 3 9 6 3 9 50 50 100 23 26 49 50 50 100 クイズ ３つの分割表 相加と2x2とが同じ表が 2個ある。どれ？ 24 18 42 6 3 9 50 50 100 20 29 49 25 26 51 24 18 42 22 18 50 6 3 9 6 3 9 50 50 100 23 26 49 50 50 100 クイズ ３つの分割表 相加と2x2とが同じ表が 2個ある。どれ？ 24 18 42 6 3 9 50 50 100 20 29 49 25 26 51 24 18 42 22 18 50 6 3 9 6 3 9 50 50 100 23 26 49 50 50 100 ハーディ・ワイン バーグ平衡かどう か 24 18 42 6 3 9 50 50 100 SNPのアレル単位でかんがえるのも ハプロタイプで考えるのも 基本は同じ • SNPのアレルの場合は2x3表の相加モデル(ト レンド)検定がある • ハプロタイプの方は、ディプロタイプがわから ないことが多く、やりようがないかもしれない ハーディ・ワインバーグ平衡のp値 HWE検定p値が小さいとき • 『サンプルは集団構造化がある母集団を代表 している』 • 『母集団を代表していない』 • 『実験がうまく行っていない』 HWE検定p値が小さいとき • 『サンプルは集団構造化がある母集団を代表 している』 – GWASならば補正方法がある • 『母集団を代表していない』 – GWASならば個々のHWE検定を問題にする必要 はない • 『実験がうまく行っていない』 – GWASにおいて、個々のHWE検定p値を利用する べきは、これ HWE検定で実験の失敗を疑う • HWE検定p値が小さいとき – 『サンプルは集団構造化がある母集団を代表し ている』 – 『母集団を代表していない』 – 『実験がうまく行っていない』 「ずれ」を見るなら、QQプロット p値が一様じゃない • p値は、一様分布に従っているから、その値を 0.01と聞けば、「あー、0.01的に珍しいことな んだ」とわかるわけですから、p値の本領は一 様分布になっていることです。 • しかしながら、実際にGWASを実施して、数十 万個のp値を算出して、その分布を見てやる と、一様分布になっていない。 ２つのアプローチ • 「本当は一様分布」なはず。「一様分布」にな るように修正してしまおう、という作戦。 • 一様分布になるわけがない。個々の検定結 果のp値を見て、対立仮説が真なのか、帰無 仮説が真なのかを選別する情報が得られれ ばよい、という作戦 一様分布に修正する作戦 ジェノミック・コントロール 単純な１要因 • 集団が完全には均一でないときに、帰無仮説の 検定結果(カイ二乗値、p値)が理論的分布から 外れる • その外れ方は、「うまく混ざっていない」という単 純な要因で説明できると仮定すると • 観測されたカイ二乗値の中央値が理論的な中 央値になるように、割り算補正すると解決するこ とが知られている • じゃあ、そうしてしまおう、というのがジェノミック コントロール法 中央値が揃うように補正 たくさんの本物がある場合～FDR～ 「合否」の基準を 一律にせず 何番目に小さいかで 手加減する 色々方法はあるが、基礎的なFDRは 直線に照らして「合否判定」 0.05 まとめ • p値は、判断するための値 – 0 から 1 、一様分布 – 一様分布であることを使って判断したい • パワーと偽陽性とは、お互い様 • たくさんのp値があったら、その特性に応じて 補正する – マルチプルテスティング補正 • 相互に関連しあう検定があったら補正は少し甘くする – FDRを使うことも 京大 統計遺伝学