• 以下のRmdに書いていないこと
    • トレンド検定は、いわゆる「トレンドカイ二乗検定」と「Cockran-Armitage検定」とが使えて、それは\frac{n-1}{n}(nは総サンプル数)で相互の統計量を換算し合う関係の、自由度1のカイ二乗分布に基づく検定であること
  • Rmdファイルの処理の仕方はこちらを参照

# 遺伝子多型のためのp値の話

## 2x2分割表で確認する、p値、オッズ比とその信頼区間
こんな2x2分割表があったとする。
```{r}
n <- 100
g <- sample(1:2,n,replace=TRUE)
g[1:(n/2)] <- sort(g[1:(n/2)],decreasing=TRUE)
p <- sample(0:1,n,replace=TRUE)
tab <- table(data.frame(p,g))
tab
```

### 検定する・推定する
カイ二乗検定をして、縦軸と横軸が無関係かどうかの検定をする。
```{r}
chisq.test(tab)
```
p値が帰ってきて、棄却するべきかどうかの情報が得られる。

p値が小さいほど、無関係という仮説(帰無仮説)を棄却することになる。
「基準を定めて」p < 0.05 なら棄却する、ということにすれば、棄却するかどうかの決断ができる。
「その基準を有意水準として、帰無仮説を棄却する」と説明し、「関連があった」と主張する。
その逆なら「関連がなかった」と主張する。
その心は、「帰無仮説が正しいと仮定したときに、このテーブルと同じかそれよりもっと偏ったテーブルを観察する確率はpくらい小さいから(pくらい大きいから)」。

オッズ比を計算する。
オッズ比は、相対危険度の推定値として使えるから。
オッズ比〜相対危険度が1であるときが、「無関係」に相当する。
推定値なので、「この値らしい」という点推定値があり、また、推定には幅があるから、だいたいこのくらいの範囲という意味で信頼区間がある(細かいことはここでは省略する)。
「だいたいこのくらいの範囲」の中に、「無関係の相対危険度〜1」が入っているかどうかで、およその帰無仮説を棄却するかどうかのめども立つ。

```{r}
library(vcd)
or <- oddsratio(tab,log=FALSE)
or
confint(or)
```

このようにデータが得られたら、たいていの場合、「検定」と「推定」をする。

「検定」するのは、「帰無仮説」を捨てるため。
「推定」するのは、相対危険度など、「因子の強さ」の程度に見当をつけるため。

たいていの場合は、「検定」をして、「関連あり」を選んだ後、「その関連の強さを数値として理解する」というプロセスをとる。

### p値とオッズ比の信頼区間との大まかな関係
たくさんの2x2表を適当に作り、カイ二乗検定p値とオッズ比の95%信頼区間を出して、その関係をプロットしてみる。
横軸にp値、縦軸に95%信頼区間をとると、オッズ比の推定値とその信頼区間は1より大きい側になることもあれば、小さい側になることもあり、また、p値の大小の順番がオッズ比の信頼区間の上下限の値の大小の順番と完璧に一致するわけではないので、多少のでこぼこはあるが、総じて、p値の大小とオッズ比の信頼区間の大小には滑らかな関係があり、p値が0.05のあたりで、オッズ比の95%信頼区間に1が入るかどうかが分離できることが見て取れる。


```{r,warning=FALSE}
library(MCMCpack)
n <- 100
n.iter <- 1000
fracs <- rdirichlet(n.iter,rep(50,4))
ps <- rep(0,n.iter)
ors <- matrix(0,n.iter,2)
for(i in 1:n.iter){
  tmp.tab <- matrix(rmultinom(1,n,fracs[i,]),2,2)
  ps[i] <- chisq.test(tmp.tab)$p.value
  ors[i,] <- confint(oddsratio(tmp.tab,log=FALSE))
}

plot(ps,ors[,1],pch=20,cex=0.5,xlim=c(0,0.2),xlab="p-value",ylab="95%CI",main="p値とオッズ比の95%信頼区間の関係")
points(ps,ors[,2],pch=20,cex=0.5)
for(i in 1:n.iter){
  segments(ps[i],ors[i,1],ps[i],ors[i,2])
}
abline(h=1,col=3)
abline(v=0.05,col=3)
```


## p値と一様分布、QQプロット
起きやすさを0から1までの数値で表したもの。
色々な検定方法などでは、異なる名前の統計量が計算されるけれど、最終的には、p値に変換することで、解釈がしやすくなる。

p値が小さければ小さいほど、観測データと同じかそれよりも起こりにくいデータは生じにくいことを意味する。
したがって、p値が小さいとき、その観測データを引き起こす前提として設定した仮説(帰無仮説)は信用しにくいことになる。

p値が小さければ小さいほど、「観測データと同じかそれよりも起こりにくいデータが得られる確率の総和」は小さくなる。
しかしながら、「観測データそのものが起きる確率」自体が小さいわけではない。
実際、p値は0から1のどの値も同程度に起きやすく起きにくい。
それを『帰無仮説が成り立つとき、p値は一様分布に従う』と言う。
p値が仮説を信じるかどうかにおいて、わかりやすいのは、この『一様分布に従う』という性質のお蔭。

一様分布とは、次のヒストグラムに示すように、等確率の分布のこと。

ヒストグラムを描いてみる。
```{r}
n <- 10000
ps <- runif(n)
hist(ps,xlab="p value")
```

分布が一様であることを示すには、小さい順に並べてプロットすることでも示せる。
```{r}
plot(sort(ps))
```

このように、「観察されたp値は一様分布である」ということを、「グラフにすると対角線上に並ぶ」ということを利用したのがQQプロット。
p値・一様分布の場合には、ただ、小さい順に並べてプロットすればよいけれど、それ以外(たとえばカイ二乗値)の場合には、小さい順に並べてプロットすると直線ではなくてカーブを描く。

「たくさんの観測されたカイ二乗値が、カイ二乗分布になるはず」と思っているときは、横軸に「カイ二乗分布からの値が同じ数だけランダムに得られたとして、それを小さい順に並べたときの理論的な値」を、縦軸に「観察された値」を取ってプロットしたものが、QQプロット。
対角線上に並んでいれば、「理論通り」と言える。

QQプロットは、「帰無仮説が成り立つときに、たくさんのp値が得られたとしたときに、得られるべきp値の理論値を横軸に、得られたp値の観測値を縦軸にしたプロット」。
観測が理論通りなら、QQプロットは対角線上、に一直線になる。
下の図はその様子。
```{r,warning=FALSE}
qqplot(qunif(ppoints(n)),ps,xlab="Theoretical",ylab="Observed",main="一様分布に従うp値のQQ plot")
```

## SNV解析の基本、2x3分割表
### 2x3分割表とは
SNVのほとんどは、2アレル(M,m)。
常染色体上のSNVは3ジェノタイプ(MM,Mm,mm)。
病気 vs. 非病気のとき、フェノタイプは0/12値。

```{r}
table2columns <- function(tab){
  m <- matrix(0,0,2)
  for(i in 1:length(tab[,1])){
    for(j in 1:length(tab[1,])){
      m <- rbind(m,cbind(rep(i,tab[i,j]),rep(j,tab[i,j])))
    }
  }
  m
}
tab <- matrix(c(10,10,10,20,10,40),2,3)
tab
gp <- table2columns(tab)
p <- gp[,1] - 1
g <- gp[,2]
#gp
#n <- 1000
#a <- 0.3
#g.f <- c(a^2,2*a*(1-a),(1-a)^2)
#g <- sample(0:2, n, replace = TRUE,prob=g.f)
#p <- rep(0,n)
#pr <- c(0.35,0.4,0.45)
#for(i in 1:3){
#  p[which(g==i-1)] <- sample(0:1,length(which(g==i-1)),replace=TRUE,prob=c(pr[i],1-pr[i]))
#}
#g[1:(n/8)] <- sort(g[1:(n/8)], decreasing = TRUE)
#p[1:(n/8)] <- sort(p[1:(n/8)], decreasing = TRUE)
#tab <- table(data.frame(p,g))
#tab
```

### いくつもある関連検定
#### トレンド検定とロジスティック回帰検定

##### トレンド検定
3つのジェノタイプはあるけれど、アレルがリスクを生み出しており、アレル本数が0,1,2本と増えるにしたがってリスクが増えると考えるのは、よい考え方だ。

その考え方に沿った検定がトレント検定とロジスティック回帰検定です。

まず、テーブルを作ってみる。
3ジェのタイプのコントロールの数は等しく、ケースの数が増えている。
```{r}
tab <- matrix(c(10,5,10,20,10,30),2,3)
#tab <- matrix(c(40,80,50,50,50,50),2,3)*3
tab
```

トレンド検定して、オッズ比も計算する。
```{r}
gp <- table2columns(tab)
p <- gp[,1] - 1
g <- gp[,2]
prop.trend.test(tab[1,],apply(tab,2,sum))$p.value
# オッズ比
or2vs1 <- oddsratio(tab[,1:2],log=FALSE)
or2vs1
#confint(or2vs1)

or3vs1 <- oddsratio(tab[,c(1,3)],log=FALSE)
or3vs1
#confint(or3vs1)
```
オッズ比は、1.5、2と計算できましたが、『トレンド検定をするという気持ち』に照らすと、このオッズ比はちょっと違います。

トレンド検定をするというのは、線形回帰をしているのと同じです。
線形回帰では、ジェノタイプごとのケースの割合が直線に乗るように推定します。
```{r,warning=FALSE}
# 線形回帰と同じこと
lm.out <- lm(p~g)
lm.out$coefficients
a <- lm.out$coefficients
plot(matrix(c(1,0,2,0,3,0,1,1,2,1,3,1),byrow=TRUE,ncol=2),pch=20,cex=t(tab)*0.5,xlim = c(-2,6),ylim=c(-1,2),xlab="genotype",ylab="phenotype",main="トレンド検定と線形回帰")
abline(a[1],a[2],col=2)
p.prop <- a[1] + (1:3)*a[2]
points(cbind(1:3,p.prop),col=2,pch=20,cex=3)
abline(h=c(0,1))
abline(v=c(1,2,3))
```

このように、ジェノタイプごとのフェノタイプ率が推定できるので、基準とするジェノタイプを決めれば、それに対するオッズ比も計算できる。
トレンド検定・線形回帰では、このようにして計算したオッズ比を想定していたことになる。
トレンド検定のp値を用いて「関連あり」と信じるのであれば、トレンド検定が支持しているオッズ比は、1.5,2ではなく、線形回帰から計算されるオッズ比になる。

```{r}
# オッズとオッズ比
odds <- p.prop/(1-p.prop)
(or2vs1.2 <- odds[2]/odds[1])
(or3vs1.2 <- odds[3]/odds[1])
```

##### ロジスティック回帰とその検定

ロジスティック回帰もトレンド検定と同様によく用いられる。
やってみる。
p値と
```{r}
glm.out <- glm(p~g,family=binomial)
summary(glm.out)$coefficients[2,4]
```
2ジェノタイプのオッズ比が得られる。
```{r,warning=FALSE}
exp(summary(glm(p~g,family=binomial))$coefficients[2,1])^(1:2)
#summary(glm(p~g,family=binomial))$coefficients
```

ロジスティック回帰によって推定されるオッズ比の特徴は、ホモ接合体のそれがヘテロ接合体のそれの二乗になっていることである。

トレンド検定では直線を引いたが、ロジスティック回帰では、シグモイドカーブを引く。

```{r,warning=FALSE}
plot(matrix(c(1,0,2,0,3,0,1,1,2,1,3,1),byrow=TRUE,ncol=2),pch=20,cex=t(tab)*0.5,xlim = c(-2,6),ylim=c(-1,2),xlab="genotype",ylab="phenotype",main="ロジスティック回帰検定")

#abline(lm.out$coefficients[1],lm.out$coefficients[2],col=2)
#points(cbind(1:3,p.prop),col=2,pch=20)
#plot(matrix(c(1,0,2,0,3,0,1,1,2,1,3,1),byrow=TRUE,ncol=2),pch=20,cex=t(tab)*0.1,xlim =  c(-2,6),ylim=c(-0.5,1.5))
x <- seq(from=-2,to=6,length=100)
tmp <- exp(glm.out$coefficients[1]+glm.out$coefficients[2]*x)
y <- tmp/(1+tmp)
points(x,y,col=3,type="l")
#abline(lm.out$coefficients[1],lm.out$coefficients[2],col=2)
abline(h=c(0,1))
abline(v=c(1,2,3))
```

トレンド検定とロジスティック回帰とを併せて描いてみる。

直線と曲線はそれなりに違うけれど、3ジェノタイプに対応するx軸上の位置はほとんど同じであることが分かる。
```{r,warning=FALSE}
plot(matrix(c(1,0,2,0,3,0,1,1,2,1,3,1),byrow=TRUE,ncol=2),pch=20,cex=t(tab)*0.5,xlim = c(-2,6),ylim=c(-1,2),xlab="genotype",ylab="phenotype",main="トレンド検定とロジスティック回帰検定")

abline(lm.out$coefficients[1],lm.out$coefficients[2],col=2)
points(x,y,col=3,type="l")

abline(h=c(0,1))
abline(v=c(1,2,3))
```

式に表して整理してみる。

トレンド検定では、
$$
\text{フェノタイプ1の割合} = a \times \text{ジェノタイプ} + b
$$
という式、
ロジスティック回帰では
$$
\frac{\text{フェノタイプ1の割合}}{\text{フェノタイプ0の割合}} = b\times e^{a \times \text{ジェタノタイプ}}
$$
となるように線を引いている。

##### トレンド検定とロジスティック回帰検定の結果に大差はない
```{r,warning=FALSE}
n.iter <- 100
p.prop <- p.log <- rep(0,n.iter)
for(i in 1:n.iter){
  a <- runif(1,0.2,0.8)
  g.f <- c(a^2,2*a*(1-a),(1-a)^2)
  g <- sample(0:2, n, replace = TRUE,prob=g.f)
  p <- sample(0:1, n, replace = TRUE)
  tab <- table(data.frame(p,g))
  p.prop[i] <- prop.trend.test(tab[1,],apply(tab,2,sum))$p.value
  p.log[i] <- summary(glm(p~g,family=binomial))$coefficients[2,4]
}
plot(log(p.prop),log(p.log),xlab="トレンド",ylab="ロジスティック",main="トレンドとロジスティックのlog(p)")
abline(0,1,col=2)
```
##### 何が違うか
トレンド検定は
> 相加モデルに相当する
> 計算が簡単
> 正確検定もできる
> 個人別のデータが不要
> 「線形回帰」に基づくオッズ比が結果に付いていないことが多い
> 共変量を組み込めない(線形回帰にすればよいが、それならロジスティック回帰をすればよい)

ロジスティック回帰検定は
> 相乗モデルに相当する
> 計算が面倒ではある(けれど計算機がやってくれるので問題はない)
> オッズ比が結果についてくることが多い(ただし、係数は対数で返ってくることが通例)
> 年齢・性別などの共変量を組み込みやすい