2014-01-16
SNP分割表〜オミックス統計学入門2014

こちらの続き
SNPのケースコントロール分割表のための配布資料
これで90分１回分か2回分と予想
以下はRmdフォーマット。knitrするとhtml化、epub化できます(もちろんフリー)→こちら)
その手間を惜しむ人向けにKindleサイトを利用中(こちらにほどなく現れる予定。1US＄(処理日の為替レートで円価決定))
SNP分割表の検定とその結果の臨床活用?オミックス統計学2014_2? 統計遺伝学
作者: ryamada
発売日: 2014/01/15
メディア: Kindle版
この商品を含むブログを見る
# SNPの関連検定

## 3行2列の分割表

一塩基多型(SNP; Single Nucleotide Polymorphism、もしくはSNV; Single Nucleotide Variation)は、DNA配列上のある１塩基が染色体によって異なる箇所のことを言う。

ヒトは常染色体を１対もつので、DNA配列のある位置には1対(２本)の染色体の塩基がある。
この個々の染色体が持つ塩基をSNPのアレルと言う。
塩基はA,T,G,Cの４塩基のいずれかである。
二種類の塩基をM,mで表せば、１対の染色体が作るアレルの対はMM,Mm,mmのいずれかとなる。
AとGとをアレルとするSNPではAA,AG,GGの３タイプがあり、この３タイプをMM,Mm,mmと書くことにする、ということである。
アレルの対をディプロタイプと呼んだり、個人が持つ遺伝子多型の型という意味でジェノタイプと呼んだりする。

同じアレルのペアとなっているディプロタイプはホモ接合型、異なるアレルのペアとなっているディプロタイプはヘテロ接合型と言う。

MMとmmはホモ接合型でありMmはヘテロ接合型である。

| |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  $n_{MM,A}$ | $n_{MM,B}$ |$n_{MM}$ |
|**Mm**| $n_{Mm,A}$ | $n_{Mm,B}$ | $n_{Mm}$ |
|**mm**| $n_{mm,A}$ | $n_{mm,B}$ | $n_{mm}$ |
|**列の和** | $n_{A}$  | $n_{B}$ | $n$ |

# 原因と因果
## ジェノタイプとフェノタイプ
遺伝因子の影響を考えるとき、遺伝子型をジェノタイプ、個体に現れる特徴をフェノタイプと言う。
上の3行2列の分割表では、ジェノタイプが行、フェノタイプが列となっている。
遺伝疫学研究では、フェノタイプ(A,B)を「病気である(A)」「病気でない(B)」とすることが多い。

ジェノタイプが行、フェノタイプが列である。
2行2列の分割表のときに、因果関係や条件とその結果を問題にするときには、原因・条件を行とし、結果を列にすることにする、と書いた。
ここでもそのルールを守っている。

ジェノタイプは個体が受精卵という単一の細胞であるときにすでに決まっているタイプであり、すべてのフェノタイプ(個体の特徴)は受精卵が発生して死亡するまでに現れるものであるから、原因と結果の関係がもしあるとすれば、ジェノタイプが原因でフェノタイプが結果である。

従って、ジェノタイプ・フェノタイプの分割表ではジェノタイプが原因・条件を表す行に、フェノタイプが結果を表す列に来る。

『関連因子解析とその臨床応用の基本』にて2行2列の表を扱ったときには、原因・条件として疾患(X,Y)を行に、検査結果(A,B)を列に取った。
この表では、ジェノタイプを原因・条件として行に、疾患フェノタイプを列に取りたいので、列の疾患フェノタイプを(X,Y)ではなく(A,B)で表していることに注意する。

## ゲノムとそれ以外の"オーム"
DNA配列を全体としてとらえてゲノムと呼ぶ。同様にDNA・染色体が修飾された状態の総体をエピゲノムと言い、転写物の総体をトランスクリプトーム、タンパク質の発現状態全体はプロテオーム、代謝系全体を指してメタボローム、表現型全体をすべて併せてフェノームと言う。

このうちDNA塩基配列以外は、原則として受精卵が発生して個体が死亡するまでの間に変化する。DNA配列も体細胞変異と呼ばれる変化がある。体細胞変異のうち、ある特定のものは癌化と関係している。
このように個体の生涯を通じて変わらないゲノムと、変わるそれ以外とを分けて、前者の変化しないゲノムを狭い意味でのゲノムと言い、その型を(狭い意味での)ジェノタイプと言うことにする。
変化するDNA配列・ゲノムは、癌の場合には「癌ゲノム」と呼ぶが、癌細胞以外にも体細胞変異はあるから、「体細胞変異ゲノム」とでも呼ぶことにする。同様の論理で「体細胞変異ジェノタイプ」と呼ぶことができる。

このように考えると、時間変化しない、「狭義のゲノム」と「時間変化のあるその他」に2分される。同様に「(狭義の)ジェノタイプ」と「フェノタイプ」と分けることができる。

「ジェノタイプ」と「フェノタイプ」とを比較するときには、時間的前後関係がわかるから、ジェノタイプが原因・条件であり、フェノタイプが結果となる。

「フェノタイプ」と「フェノタイプ」とを比較する場合には、時間的前後関係が必ずしも明らかではないから、関連を調べることはできても、因果関係と判断することは難しくなる。

## 3行2列の分割表
3行2列の分割表に話を戻す。



2行2列の分割表の色々な計算値を表にすると以下の通りであった。
3行2列の分割表でも基本的には同じことができるはずであるが、表が異なるので、そのままと言うわけにはいかない。

そのまま適用できるのか、修正・拡張が必要なのかを、個別に検討する。

## 行と列とが独立であるかどうかを評価する場合
### 期待値表

行と列とが独立であると仮定したときの期待値の計算式は

$ne_{i,j} = \frac{n_i n_j}{n}$

と2行2列の表の場合と変わらず、表にすれば

|期待値表 |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  $ne_{MM,A} = \frac{n_{MM}n_A}{n}$ | $ne_{MM,Y}=\frac{n_{MM}n_B}{n}$ |$n_{MM}$ |
|**Mm**|   $ne_{Mm,A} = \frac{n_{Mm}n_A}{n}$ | $ne_{Mm,Y}=\frac{n_{Mm}n_B}{n}$ | $n_{Mm}$ |
|**mm**|   $ne_{mm,A} = \frac{n_{mm}n_A}{n}$ | $ne_{mm,Y}=\frac{n_{mm}n_B}{n}$ | $n_{mm}$ |
|**列の和** | $n_{A}$  | $n_{B}$ | $n$ |

となる。
### $\chi^2$値

この観測表を期待値表に引き比べて、$\chi^2$値を計算することができる。
その式は2行2列の場合と同じで

$\chi^2 = \sum_{i=MM,Mm,mm;j=A,B}\frac{(n_{i,j}-ne_{i,j})^2}{ne_{i,j}}$

である。

### p値、自由度は2

$\chi^2$値があれば、p値を出すこともできる。
$\chi^2$値からp値を得るには、自由度というものを指定する必要がある。

分割表がa行b列であるときには、自由度が$(a-1)(b-1)$であるので、ここで得られた$\chi^2$値は自由度を$(3-1)(2-1)=2$としてp値に換算することができる。

### 2行2列の式が使えるかどうか

『関連因子解析とその臨床応用の基本』にて2行2列の表で計算された諸指標は以下の通り。
今3行2列の表にしたので(X,Y)→(MM,Mm,mm)となる。

各指標の式にXとYとが入り乱れて現れている場合には、X,YにMM,Mm,mmのどれを対応付けたらよいのかをすぐには決め難い。
決め難いということは、そのまま使うことができないことを意味する。
逆に、X,Yの片方しか現れない式を持つ指標はそのまま使える。


|そのまま使える|名称 |式 | 意味 | 関係 |
|:---|:---:|:---:|:---:|:----:|
|〇|$\chi^2$|$\sum_{i=MM,Mm,mm;j=A,B}\frac{(n_{i,j}-ne_{i,j})^2}{ne_{i,j}}$ | 帰無仮説からの距離 |p値に換算できる |
|〇|p値| 省略 | 帰無仮説を信じない程度 | $\chi^2$から計算できる |
|×|オッズ比 | $\frac{n_{X,A}n_{Y,B}}{n_{X,B}n_{Y,A}}$  | 行と列との関連の強さ |$\frac{LL(X|A)}{LL(X|B)}$  |
|〇|感度| $\frac{n_{X,A}}{n_{X}}$ | Xという条件でのAの条件付き確率|A観察のときのXの尤度と同じ|
|〇|特異度|$\frac{n_{Y,B}}{n_{Y}}$| Yという条件でのBの条件付き確率|B観察のときのYの尤度|
|〇|A観察のときのXの尤度|$\frac{n_{X,A}}{n_{X}}$|Xという条件でのAの条件付き確率|感度と同じ|
|〇|B観察のときのXの尤度|$\frac{n_{X,B}}{n_{X}}$|Xという条件でのBの条件付き確率||
|〇|A観察のときのYの尤度|$\frac{n_{Y,A}}{n_{Y}}$|Yという条件でのAの条件付き確率||
|〇|B観察のときのYの尤度|$\frac{n_{Y,B}}{n_{Y}}$|Yという条件でのBの条件付き確率|特異度と同じ|
|×|$LL(X/Y|A)$|$\frac{n_{X,A}}{n_{Y,A}}\times \frac{n_{Y}}{n_X}$|A観察のときの尤度比(X/Y)||
|×|$LL(X/Y|B)$|$\frac{n_{X,B}}{n_{Y,B}}\times \frac{n_{Y}}{n_X}$|B観察のときの尤度比(X/Y)||
|〇|Xの事前確率|$\frac{n_X}{n}$|||
|〇|A観察のときのXの事後確率|$\frac{n_{X,A}}{n_A}$||PPVと同じ|
|〇|A観察のときのYの事後確率|$\frac{n_{Y,A}}{n_A}$||1-PPV|
|〇|B観察のときのXの事後確率|$\frac{n_{X,B}}{n_B}$||1-NPVと同じ|
|〇|B観察のときのYの事後確率|$\frac{n_{Y,B}}{n_B}$||NPV|
|〇|PPV|$\frac{n_{X,A}}{n_A}$||A観察のときのXの事後確率と同じ|
|〇|NPV|$\frac{n_{Y,B}}{n_B}$||B観察のときのYの事後確率と同じ|

### 3行2列の表に特徴的な指標

#### 比は２つを比べるもの

##### 尤度比

2行2列の式をそのまま使えない指標について考えてみる。
それは2行2列と3行2列との違いを知る手掛かりになる。

尤度比$LL(X/Y|A) = \frac{n_{X,A}}{n_{Y,A}}\times \frac{n_{Y}}{n_X}$はX,Yの両方を式に含んでいるので、MM,Mm,mmにそのままは使えない。

なぜならA(病気である)という観察がなされたときに、それを引き起こす条件がX,Yの2種類であるならば、その２条件の条件付き確率(尤度)の比をとることはできるが、MM,Mm,mmという３条件を考えているときには、単純に２ジェノタイプを取って比を取ることにわかりやすい意味がないからである。

尤度比が２つの条件の比であるのに、3条件あるという事情がふさわしくないから、と言い換えることもできる。

##### オッズ比

オッズ比$\frac{n_{X,A}n_{Y,B}}{n_{X,B}n_{Y,A}}$も「比」である。
従って、３条件あるときに利用するには、工夫が必要になる。

ジェノタイプが疾患フェノタイプに影響を与えるかどうか、ということに興味があるときには、「疾患を最も起こしにくいジェノタイプ」を基準にして、「あるジェノタイプはどれくらい疾患を起こしやすいか」を知りたい。

今、mmジェノタイプが「疾患を最も起こしにくいジェノタイプ」であるとすれば、MM,Mm,mmの３ジェノタイプのオッズ比は

|オッズ比|式 ||
|:---:|:---:|:---:|
|MM vs. mm | $\frac{n_{MM,A}n_{mm,B}}{n_{MM,B}n_{mm,A}}$  | |
|Mm vs. mm | $\frac{n_{Mm,A}n_{mm,B}}{n_{Mm,B}n_{mm,A}}$  | |
|MM vs. mm | $\frac{n_{mm,A}n_{mm,B}}{n_{mm,B}n_{mm,A}}$  | これは1|

これは(もっともリスクの低いジェノタイプを基準にしたときの)「レラティブ・リスク(相対危険度)」と呼ばれ、分割表のオッズ比がレラティブ・リスクのよい推定値になっていることが多い。

> リスク因子を持っているために、リスク因子が無かった場合に比べて、何倍病気になりやすいか

という値である。したがって、2なら2倍なりやすい、0.5なら1/2倍なりやすい。

## 遺伝的モデルを考慮する〜ジェノタイプに順序がある〜

次の２つの表の$\chi^2$値とオッズ比を計算してみる。
オッズ比の計算に当たっては、mmジェノタイプを基準にすることにする。

| |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  15 | 10 |25 |
|**Mm**| 20 | 10 | 30 |
|**mm**| 10 | 10 | 20 |
|**列の和** | 45  | 30 | 75 |

| |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  20 | 10 |30 |
|**Mm**| 15 | 10 | 25 |
|**mm**| 10 | 10 | 20 |
|**列の和** | 45  | 30 | 75 |

この２つの表はMMの行とmmの行の値が入れ替わっているけれど、それ以外はすべて同じである。

$\chi^2$の計算式$\sum_{i=MM,Mm,mm;j=A,B}\frac{(n_{i,j}-ne_{i,j})^2}{ne_{i,j}}$は、３つのジェノタイプについてそれぞれ別々に計算してそれを足し合わせているだけなので、２つの表の$\chi^2$値は等しい。

オッズ比はどうだろうか。

> 上の表では、$OR(MM/mm)=\frac{15\times 10}{10\times 10} =1.5$,$OR(Mm/mm)=\frac{20\times 10}{10\times 10}=2$、$OR(mm/mm)=1$。


> 下の表では、$OR(MM/mm)=\frac{20\times 10}{10\times 10} =2$,$OR(Mm/mm)=\frac{15\times 10}{10\times 10}=1.5$、$OR(mm/mm)=1$。

オッズ比だけを抜き出すと、1.5, 2, 1 と、2, 1.5, 1 とである。それぞれは次のように言い換えられる。

> ジェノタイプmmを基準にして、２つのアレルの片方をMに変えたときのレラティブ・リスクが2倍となり、２つのアレルの両方をMに変えるとレラティブ・リスクは少し小さくなり1.5倍となる。


> ジェノタイプmmを基準にして、２つのアレルの片方をMに変えたときのレラティブ・リスクが1.5倍となり、２つのアレルの両方をMに変えるとレラティブ・リスクはさらに大きくなり2倍となる。

どちらの文も真実かもしれないけれど、どちらかと言えば、下の文「2,1.5,1」というオッズ比の並びの方が納得が行くように感じられる。
実際、３ジェノタイプMM,Mm,mmのレラティブ・リスクに順序があることを前提にして解析することの方が普通である。


### ジェノタイプの順序を考慮したモデル

順序があると仮定して解析をすることを、順序があるモデルをあてはめて解析するともいう。
ジェノタイプの順序のモデルにはいくつか代表的なものがあり、モデルにはそれに対応する解析手法がある。


|モデル名|レラティブ・リスクの設定($r>1$)|解析手法名|
|:---:|:---:|:---:|
|優性モデル|RR(MM)=RR(Mm)=$r$,RR(mm)=$1$|2x2表化して$\chi^2$検定|
|劣性モデル|RR(MM)=$r$,RR(Mm)=RR(mm)=$1$|2x2表化して$\chi^2$検定|
|相加モデル|RR(MM)=$2r+1$,RR(Mm)=$r+1$,RR(mm)=$1$|トレンド検定・コクラン・アーミテージ検定|
|相乗モデル|RR(MM)=$r^2$,RR(Mm)=$r$,RR(mm)=$1$|ロジスティック回帰で検定|

このうち、優性モデルと劣性モデルは、いわゆるメンデル型遺伝病に見られるように、責任アレルを１本でも持っていれば発病し(優性モデル)、責任アレルを１本しかもっていない場合には保因者ではあるが疾患フェノタイプは示さず、2本持っている場合には発病する(劣性モデル)というものである。

このようにはっきりと区別できる場合とは異なり、なんらかの影響力を持つアレルがあり、1コピーあるよりも２コピーある方がその影響力が大きいと考えれば、優性モデルと劣性モデルの中間のモデルが必要になる。大多数の遺伝的影響がこのモデルに合致すると思われる。
相加モデルも相乗モデルもどちらもジェノタイプMmに中間的なレラティブ・リスクを設定している。実際、リスクアレルがフェノタイプの発現にどのように影響するかは、個々の分子の生物学的機能の大きさ、発現量、複雑で精密な分子ネットワークによる発現制御、時と場所による発現パターンなどの影響など、複数の要素の複合として現れるので、相加モデルにせよ相乗モデルにせよ、リスクの大きさがモデル式通りになる保証はなく、そうなっていないと考える方が自然である。

### ジェノタイプの順序を考慮するのはなぜか

上述のように、ジェノタイプの順序を考慮したモデル(相加モデル、相乗モデル)はあるが、そのモデルが必ずしも正しいわけではない。
では、なぜそのようなモデルを設定して解析するのだろうか。

端的に言えば、モデルに近いデータを優先して取りあげ、モデルから遠いデータに重きを置かないためである。

例を示す。
次のような３つのデータが得られたとする。

|データ１ |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  15 | 10 |25 |
|**Mm**| 20 | 10 | 30 |
|**mm**| 10 | 10 | 20 |
|**列の和** | 45  | 30 | 75 |

|データ２ |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  20 | 10 |30 |
|**Mm**| 15 | 10 | 25 |
|**mm**| 10 | 10 | 20 |
|**列の和** | 45  | 30 | 75 |

|データ３ |A | B | 行の和 |
|:---:|:---:|:---:|:----:|
|**MM**|  16 | 9 |25 |
|**Mm**| 19 | 11 | 30 |
|**mm**| 10 | 10 | 20 |
|**列の和** | 45  | 30 | 75 |

それぞれのオッズ比はmmを基準として、

|データ番号 |OR(MM/mm)|OR(Mm/mm)  |
|:--:|:--:|:--:|
|データ１|1.5|2|
|データ２|2|1.5|
|データ３|1.78|1.73|

この３つの表に、３ジェノタイプの順序を考慮しない検定($\chi^2$検定(自由度2))と、３ジェノタイプの順序を考慮した検定(トレンド検定)を実行してみます

```{r}
d1 <- matrix(c(15,10,20,10,10,10),byrow=TRUE,ncol=2)
d2 <- d1[c(2,1,3),] # 行を入れ替えているだけ
d3 <- matrix(c(16,9,19,11,10,10),byrow=TRUE,ncol=2)
d1
d2
d3
# 順序を考慮しないchi-squared test
chisq.test(d1,correct=FALSE)
chisq.test(d2,correct=FALSE)
chisq.test(d3,correct=FALSE)
# トレンド検定。引数は第1列と、列の和
prop.trend.test(d1[,1],apply(d1,1,sum))
prop.trend.test(d2[,1],apply(d2,1,sum))
prop.trend.test(d3[,1],apply(d3,1,sum))
```

ジェノタイプの順序を考慮しない検定では、データ１とデータ２との$\chi^2$が等しく、データ３のそれはより小さくなっている(X-squared = 1.389,1.389,1.139)。
これは、ジェノタイプ順序を考慮しないのであれば、データ１とデータ２は同程度にジェノタイプがフェノタイプに影響していると考えるのがよく、データ３はそれらに比べて影響があると考えにくい、と読み取れる。
一方、順序を考慮したトレンド検定の方では、統計量が0.3731,1.37,0.8396とあるように、もっとも関連がありそうなのはデータ２、次がデータ３、最後にデータ１が来ることがわかる。

このようにジェノタイプに順序を入れたモデルを適用することで、モデルに合致しないデータには小さい関連指標値が出るようになる。
このような指標値を用いることで、モデルに合致したデータが優先的にスクリーニングすることができる。

## 活用する
以上が、ディプロタイプ型ジェノタイプと2値型フェノタイプに観察される3行2列分割表の基本事項である。
以下は、その知識を活用するための練習である。

### データを作ってみる
#### Hardy-Weinberg 平衡とランダム・メイティング

今、ある集団にアレルMとアレルmとが割合$P(M)=0.7$,$P(m)=0.3$で存在しているとする。
これはアレルMの染色体が全体の7割、アレルmの染色体が３割である、という意味である。
この集団がランダム・メイティングしており、このアレルが出生・成長・生殖に影響がないとすると、この集団の個体ジェノタイプMM,Mm,mmの頻度は$P(MM)=P(M)^2$,$P(Mm)=2P(M)P(m)$,$P(mm)=P(m)^2$になると考えられる。
アレル頻度とジェノタイプ頻度との関係がこのような関係にあることをHardy-Weinberg平衡(HWE)にある、と言う。

リスクアレルmの頻度が$P(m)=0.2$であり、ジェノタイプ頻度がHWEになっているような集団から、$N_s=100000$人をランダムにサンプリングし、そのジェノタイプを長さ$N_s$のベクトルにし、それを集計せよ。ただし、ジェノタイプMM,Mm,mmにはリスクアレルmの本数0,1,2を対応づけること。

```{r}
P.m <- 0.2
N.s <- 100000
P.M <- 1-P.m
P.MM <- P.M^2; P.Mm <- 2*P.M*P.m; P.mm <- P.m^2
P.g <- c(P.MM,P.Mm,P.mm)
G <- sample(0:2,N.s,replace=TRUE,prob=P.g)
table(G)
```

#### ジェノタイプ別に有病率を設定する

リスクアレルを持たない場合、この病気の有病率を$Q(MM)=0.1$とし、レラティブ・リスクには相乗モデルを仮定し、ジェノタイプMmのレラティブ・リスクを$r=1.5$とし、$N_s$人に確率的にフェノタイプXを割り当てよ。
```{r}
Q.MM <- 0.1
r <- 1.5
Q.Mm <- Q.MM * r
Q.mm <- Q.MM * r^2
Q.G <- c(Q.MM,Q.Mm,Q.mm)
Q.G
Q.s <- Q.G[G+1] # これで各個人のRRが得られる
plot(G,Q.s) # 念のため確認する
plot(jitter(G),jitter(Q.s)) # 点が重なって不安なので少しぶれを入れると安心するかもしれない
R <- runif(N.s) # 各個人用の乱数
X <- R < Q.s # 病気がTRUE,非病気がFALSE
X <- as.numeric(X) # TRUEを1、FALSEを0に換える
table(data.frame(G=G,X=X))
```

#### データを作る処理を関数にする

人数、アレル頻度、レラティブ・リスクを与えて、個人別のフェノタイプを発生させ、さらにそれを3行2列の表にした。この処理を関数にする。
```{r}
# 入力変数
N.s <- 100000
P.m <- 0.2
Q.MM <- 0.1
r <- 1.5

# 出力
# ジェノタイプG,フェノタイプX,そのテーブルTa
# まずはその枠だけ作ります
my.make.SNPdata <- function(N.s,P.m,Q.MM,r){
  
  return(list(G=G,X=X,Ta=Ta)) # GをGとして、XをXとして、TaをTaとして返す
}
```

関数の入力と出力ができたら、あとは入力の変数と出力の変数をつなぐ処理を関数の{ }の内側に書き込む。
```{r}
# 入力変数
N.s <- 100000
P.m <- 0.2
Q.MM <- 0.1
r <- 1.5

# 出力
# ジェノタイプG,フェノタイプP,そのテーブルTa
# まずはその枠だけ作る
my.make.SNPdata <- function(N.s,P.m,Q.MM,r){
  P.M <- 1-P.m
  P.MM <- P.M^2; P.Mm <- 2*P.M*P.m; P.mm <- P.m^2
  P.g <- c(P.MM,P.Mm,P.mm)
  G <- sample(0:2,N.s,replace=TRUE,prob=P.g)
  Q.Mm <- Q.MM * r
  Q.mm <- Q.MM * r^2
  Q.G <- c(Q.MM,Q.Mm,Q.mm)
  Q.s <- Q.G[G+1] # これで各個人のRRが得られる
  R <- runif(N.s) # 各個人用の乱数
  X <- R < Q.s # 病気がTRUE,非病気がFALSE
  X <- as.numeric(X) # TRUEを1、FALSEを0に換える
  Ta <- table(data.frame(G=G,X=X))
  return(list(G=G,X=X,Ta=Ta))
}
# 使ってみます
test.out <- my.make.SNPdata(N.s,P.m,Q.MM,r)
test.out$Ta
```

### 関連を見出す

データをシミュレーション作成できるようになったので、これを用いて解析することにする。
大きく分けると次の２点をになる。
> 関連があることを、「関連がないという帰無仮説を棄却」するに足るp値が得られたことで示し、

> 関連の強さをジェノタイプ別のレラテイィブ・リスクの推定値として示す。その推定値は分割表から計算するオッズ比であったり、ロジスティック回帰の係数から算出したものであったりする。推定値は点推定値と信頼区間とで示す。

では、上で作成した集団においてケース・コントロール関連解析を実施したものとして、解析をしてみる。
ケース・コントロール関連解析では、フェノタイプごとにサンプリングする。
今、ケースとコントロールをそれぞれ$N_{ca}=500,N_{co}=500$人サンプリングして相加モデルで検定してみる。

```{r}
N.ca <- 500
N.co <- 500
# フェノタイプtest.out$Xの値が1,0の番号を取りだす
cases <- which(test.out$X==1)
conts <- which(test.out$X==0)
# N.ca,N.co人をサンプリングする
cases.sample <- sample(cases,N.ca)
conts.sample <- sample(conts,N.co)
samples <- c(cases.sample,conts.sample)
# サンプルのジェノタイプとフェノタイプ
G.sample <- test.out$G[samples]
X.sample <- test.out$X[samples]
# サンプルの分割表
Ta.sample <- table(data.frame(G.sample,X.sample))
Ta.sample
tr.out <- prop.trend.test(Ta.sample[,1],apply(Ta.sample,1,sum))
tr.out
```

大きな$\chi^2$値と小さなp値が得られる。

> ジェノタイプとフェノタイプには関連がないと考えにくい。
> その判断は相加モデルによって検定した。
> その意味は、３ジェノタイプのリスクが$2r+1,r+1,1$となる場合を最優先としてジェノタイプとフェノタイプの関連を評価・検定した、という意味である。

次にMMを基準として、オッズ比を計算する。
オッズ比は２つのジェノタイプが作る2行2列の表を用いて計算する。
そのためにパッケージをインストールしてそれを読み込んだ上で、オッズ比とその95%信頼区間を計算する関数を使う。
その関数は、1行目が因子あり、1列目が病気ありを前提にしているので、関数への2行2列データの受け渡しに際して、それに合うようにしている。
```{r}
# "epiR"パッケージのインストールが済んでいなければ、以下のコメントアウトした行を実行します
#install.packages("epiR")
library(epiR)
OR.Mm <- epi.2by2(Ta.sample[c(2:1),c(2:1)])
OR.mm <- epi.2by2(Ta.sample[c(3,1),c(2:1)])
OR.Mm
OR.mm
```

複数の指標が示されるが、ここでは、Odds ratioに着目する。
点推移定値と信頼区間とが x (y,z) という形で示されている。

データをシミュレーションした際に設定したレラティブ・リスクは、Mm,mmそれぞれ、$1.5,1.5^2=2.25$であった。
サンプルから計算したオッズ比はその値から遠くはないが、等しいわけではない。

また、検定モデルとしては相加モデルを用いたが、それはMmのリスクがmmのそれよりも小さいことを仮定したかったからであって、MM,Mm,mmのリスクが$1,1+r,1+2r$でなければならないと仮定したわけではないので、Mm,mmのレラティブ・リスクは、上記のような計算をしてx (y,z)と信頼区間付きで考えておけばよい。

さて、場合によってはロジスティック回帰を手法として採用するかもしれない。
以下のように実施することができる
```{r}
fit <- glm(X.sample~G.sample,family=binomial())
summary(fit) # display results
```

ここに現れるのは、相乗モデルが最も適切だと考えた上で、ジェノタイプとフェノタイプとに関連があるかどうかの評価結果である。
いくつもの指標と数値が現れるが、prop.trend.test()で算出した$\chi^2$値とp値とに相当するのはCoefficients:に引き続くG.sample のz value とPr(>|z|)の値である。

G.sampleとして与えた、ジェノタイプの値がフェノタイプに与える影響をzという値で表し、それに対応するp値をPr(|z|)と書いている。
$z^2$の値が$\chi^2$値と対応するので、prop.trend.test()の$\chi^2$値とこのz値とを、両手法のp値とともに比べてみる。
```{r}
print(c(tr.out$statistic,(summary(fit))$coefficients[2,3]^2))
print(c(tr.out$p.value,(summary(fit))$coefficients[2,4]))
```
だいたい同じであることがわかる。
関連がないという帰無仮説を棄却する上でだいたい同じである、と言うことである。

ジェノタイプMMを基準にしたオッズ比とその信頼区間はモデルによらず算出することができるので、ロジスティック・回帰モデルで解析した後にも、その値を計算し、それを利用することはもちろん可能である。
ロジスティック回帰を用いた場合には、リスク・アレル１本の寄与の程度が計算されているので、それに基づいて、Mm,mmジェノタイプのレラティブ・リスクを推定すれば以下のようになる。

Mmは回帰係数$\beta$を用いて$e^{\beta}$で計算されるから、以下のようになる。
初めに点推定値、次に信頼区間を示す。
```{r}
print(c(exp(coef(fit))[2],exp(confint(fit))[2,])) # 95% CI 信頼区間
```

mmは2本分であるから、1本分の2乗であり、次のようになる。

```{r}
print(c(exp(coef(fit))[2]^2,exp(confint(fit))[2,]^2)) # 95% CI 信頼区間
```

この値も、シミュレーションで用いた1.5,2.25と完全には一致しないが、推定値としては妥当であることが確かめられる。

#### 結果を応用する

みずからサンプルを収集し関連解析を行う場合には、上記のような手順を踏む。
では、そのような研究成果が出ているとき、それはどのように臨床応用すればよいかを考える。

3行2列の分割表とprop.trend.test()のp値とMM基準のジェノタイプ別オッズ比が報告されたとする。
```{r}
Ta.sample
tr.out$p.value
OR.Mm$rval$OR[c(1,4,5)]
OR.mm$rval$OR[c(1,4,5)]
```

ある人のジェノタイプがmmだったとき、この報告に基づいて

> あなたのジェノタイプはmmなので、あなたが病気になる確率は、MMの人に比べて

```{r}
paste(OR.mm$rval$OR[1],"倍")
```

高いです、と説明することは可能である。
しかしながら、その説明を受けても、結局、その人は自分がどれくらいの確率で病気になるのかがわからない。

これは、ケース・コントロールサンプリングをしたときに、ケースの人数とコントロールの人数の比が、
集団におけるケースとコントロールの人数の比と異なっているために生じる問題である。

今、ある人が病気かそうでないかの情報が全くないものとする。
この人のジェノタイプが判明したとき、この人が病気なのかそうでないのかを予想することを考える。

これは、疾患と臨床検査の分割表
| |A(検査が陽性である) | B(検査が陰性である) | 行の和 |
|:---:|:---:|:---:|:----:|
|**X(病気であるならば)**|  $n_{X,A}$ | $n_{X,B}$ |$n_{X}$ |
|**Y(病気でないならば)**| $n_{Y,A}$ | $n_{Y,B}$ | $n_{Y}$ |
|**列の和** | $n_{A}$  | $n_{B}$ | $n$ |
で言うところの、A検査が陽性であるをある特定のジェノタイプである、と読み換えたときの
PPVと同じである。

この計算をするためには、X,Yの人数比がわからなければ、分割表自体がPPVの計算にそぐわない。

今、ケース・コントロールスタディを実施したので、$n_{ca}=n_{X},n_{co}=n_Y$が一般集団のそれとは異なっている。
一般集団での$n_{ca}=n_X,n_{co}=n_Y$の割合は有病率Qによって$Q,1-Q$となるから、有病率の情報が必要であることがわかる。

今、$Q=0.12$だと考えられていて、ある人のジェノタイプがmmであると判明したとする。
Qを考慮して、mmかそうでないかで2行2列の表の値を再計算することにする。
ここではジェノタイプを列に取ることにする。ジェノタイプとの関係では、病気か否かはA,Bで表して来たことに注意する。

| |mm | mm以外 | 行の和 |
|:---:|:---:|:---:|:----:|
|**A(病気である)**| $n_1=$? | $n_2=$? |$n\times Q$ |
|**B(病気でない)**| $n_3=$? | $n_4=$? | $n \times (1-Q)$ |
|**列の和** | $n_1+n_3=$?  | $n_2+n_4=$? | $n$ |

で?の値をケース・コントロールの観察表から算出したい。

$n_1 : n_2 = n_{mm,A} : n_{MM,A}+n_{Mm,A}$ であり、$n_3 : n_4 = n_{mm,B} : n_{MM,B}+n_{Mm,B}$であることを考慮して、

| |mm | mm以外 | 行の和 |
|:---:|:---:|:---:|:----:|
|**A(病気である)**| $nQ\frac{n_{mm,A}}{n_A}$ | $nQ\frac{n_{MM,A}+n_{Mm,A}}{n_{A}}$ |$n\times Q$ |
|**B(病気でない)**| $n(1-Q)\frac{n_{mm,B}}{n_{B}}$ | $n(1-Q)\frac{n_{MM,B}+n_{Mm,B}}{n_{B}}$ | $n \times (1-Q)$ |
|**列の和** | $n_1+n_3$ | $n_2+n_4$| $n$ |

という表が得られます。

この表に基づいて、いわゆる$PPV = \frac{n_1}{n_1+n_3}$を計算すれば、ジェノタイプがmmとわかったときの病気である確率が得られます。
式で書くと

$\frac{Q\frac{n_{mm,A}}{n_A}}{Q\frac{n_{mm,A}}{n_A}+(1-Q)\frac{n_{mm,B}}{n_{B}}}$

実際にQを与えてこの値を計算してみると
```{r}
Q <- 0.125
n.mm.A <- Ta.sample[3,2]
n.A <- sum(Ta.sample[,2])
n.mm.B <- Ta.sample[3,1]
n.B <- sum(Ta.sample[,1])

(Q*n.mm.A/n.A)/(Q*n.mm.A/n.A + (1-Q)*n.mm.B/n.B)
```

となる。
ジェノタイプ情報がなかったときに$Q=0.125$と予想していたものが、かなり大きな値となった。
これが、ジェノタイプが判明したときの病気である事後確率である。

### ジェノタイプ情報を用いた疾病確率の計算に関する注意

上記の計算で明らかにしたのは、実は、「病気かどうかわからない」状態のときに、「ジェノタイプがmmと判明したとして」、その人が「有病率Q」で表される「疾病集団の一員である」(事後)確率である。

今、Qが生まれてから死ぬまでにこの病気にかかる人の割合(生涯有病率)を表し、ジェノタイプを調べ、病気になる確率を計算した対象者が新生時であれば、この計算でよい。

しかしながら、そもそも生涯有病率にしろ、ある一時点での有病率(時点有病率)も正確な数字は不明なので、上の計算は、幅のあるおよその数字である。

また、実際には「ある年齢の人」が「これから病気になる」確率を知りたいとすると、「その年齢までは病気になっていない人のうち、その年齢以降に病気になる割合」をQとしなくてはならない。
この割合は、年齢とともに下がるが、どのように下がるかがわかっている疾患はあまりない。
たとえば、20歳の人が死ぬまでに肺癌になる確率を訪ねる場合と、85歳の人が同じ質問をする場合とでは異なる考え方が必要だと言うことである。

このように、有病率Qに関する曖昧さや難しさの他に、レラティブ・リスクの信頼区間のことも考慮に入れる必要がある。
上の計算では、観測表の度数をそのまま使っているので、オッズ比で言うところの点推定値のみを使い、信頼区間は無視していることとに相当する。
有病率の曖昧さや年齢の要素の曖昧さがあるので、オッズ比に関するところだけを厳密にしても、適当とは考えない。
したがって、今の段階では、

> 「複数の理由で曖昧さのある点推定値」を「ジェノタイプ情報を得た後の発病事後確率」の目安として計算できる

ということに留める。
時間があれば、この点について回を改めて取り扱う。