2012-07-14

ケースだけ検査してある

Outcome-dependent sampling カイ二乗分布指数分布 Fisher's combined probability test

こちらから
状況
- ケース・コントロールサンプルがあって、発病SNPに興味がある
- ケースには、その重症度を表す(かのような）量的検査がある
- 具体的には、関節リウマチとその自己抗体検査であるACPA、とか
- 特徴的なのは
  - ケースにはACPAが(ほぼ)もれなく検査されていること
  - コントロールは一般集団だとすると、ACPAはほぼもれなく検査されていないこと
- このようなときに、発病に関与しつつ、ACPAを増やすのにも関与するような、「病気を起こして押し進める因子」を探すにはどうすればよいか、と言う話
データ
- 人数行ｘ３カラム(ケース・コントロールの別、SNPの３ジェノタイプ、検査値)
発病とジェノタイプの検定をしよう
検査値とジェノタイプの検定をしよう
- 検査値のあるケースだけでやる方法もある
- 検査値のないコントロールに適当に値を埋める方法もある(論文が提唱している方法)
  - 発病と検査値増大に効いているSNPを探す(その向きがそろっているものを探すと言うモデルに基づいて検定する)、という場合には、コントロール全員にケースの検査値の最小値を仮定して検定してもタイプ１エラー的にはよいことを示している
- 検査値とジェノタイプに関しては、ジェノタイプに順序を入れたF統計量の算出(自由度が(1,送検対数-2)で返ってくる→Fを自由度１のカイ二乗分布で評価してｐ値化する
発病に関する検定と検査値との検定の結果を併せる
- ２種類のｐ値がある
- このｐ値を統計量化したい
- Fisher's method to combine p-values(Wiki)
- 「合算に便利な『距離』」的な扱いをするのにはカイ二乗分布がよい
- ｐ値として出ているものが帰無仮説のもとで一様分布である、ということを利用するのであれば、ｐ値に対応するカイ二乗値として自由度は何でもよいのだが…(実際、p値を任意の自由度で統計量に変え、その和を取って、自由度の２倍を自由度としてp値に戻すと、その結果、得られるp値は一様分布になる。一様分布になるが、 $-2\log(p)$ を使って(自由度２にして)やる場合と何が違うかと言えば、併せた統計量の大小の順序に入れ替わりが生じること。自由度２の場合には、二つのp値の積の大小を統計量の順序にしており、それは、非常に納得しやすい順序である
- $-2 \log(p)$ を使って自由度２のカイ二乗統計量にするのは、「便利」
- これができるのは、自由度２のカイ二乗分布が指数分布だから
- ２つのp値を自由度２のカイ二乗統計量に対応付けて、その２つのカイ二乗統計量を単純に足せば、それは自由度４のカイ二乗分布で評価すればいいでしょう…という話なのは、２つのp値が独立な検定から出てきたときの場合
非独立なp値を併せたときには、カイ二乗分布ではなくガンマ分布、ただしガンマ分布のパラメタを推定しよう
- 独立な場合
  - Brown's method
    - $\sum_{i=1}^k -2\log(p_i) \sim \chi^2(2k)$
- 非独立な場合
  - Kost's method(こちら)
    - $-2\log(p_1)-2\log(p_2) \sim c \chi^2_d = \Gamma(2c,\frac{d}{2})$ , $c$ はスケールパラメタ、 $d$ は自由度
  - ガンマ分布のスケールパラメタと自由度とがわからないのが問題だが、統計量がたくさん取れれば、そのサンプルの平均と分散とが合うようにパラメタ推定してやりましょう、と言う話
  - その推定もいやだったら、 $\sum -2\log(p)$ の順序だけを気にするのでも「順序をつけるだけのスクリーニング」の用途には使えたり
データを作ろう

# outcome-dependent sampling
n.case <- 1000
n.control <- 1000

af <- 0.3
gf <- c(af^2,2*af*(1-af),(1-af)^2)
# 発病リスク(rsk = 1 は帰無仮説)
rsk <- 1
risk <- c(rsk, rsk^0.5,1)
gf.case <- gf*risk
gf.case <- gf.case/sum(gf.case)

# disease status
k <- c(rep(1,n.case),rep(0,n.control))
# genotype
g <- c(sample(2:0,n.case,prob=gf.case,replace=TRUE),sample(2:0,n.control,prob=gf,replace=TRUE))
# ACPA
# q.r = 0 は帰無仮説
q.r <- 0
acpa <- c(rnorm(n.case) + q.r*g[1:n.case],rep(NA,n.control))

df.all <- data.frame(k=k,g=g,acpa=acpa)
head(df.all)

df.case <- data.frame(k=k[1:n.case],g=g[1:n.case],acpa=acpa[1:n.case])
df.control <- data.frame(k=k[(n.case+1):(n.case+n.control)],g[(n.case+1):(n.case+n.control)],acpa[(n.case+1):(n.case+n.control)])

library(ggplot2)
gp <- ggplot(df.case, aes(x=factor(g),y=acpa)) + geom_boxplot()
gp

F統計量をカイ二乗分布(自由度１)で評価してよいのかどうかは、「帰無仮説を前提にやってみれ」ばよい
- まず、F統計量を論文の式でも出せるし、ただ、線形回帰してもできることを確かめる

lm.out <- lm(data = df.fill, acpa ~ g)
summary(lm.out)
snp.F <- function(g,q){
	x <- c(0,0.5,1)
	x. <- mean(g)/2
	#x. <- mean(g)
	Y. <- mean(q)
	n <- length(g)
	ret0 <- (sum((g/2-x.)*(q-Y.)))^2
	#ret0 <- (sum((g-x.)*(q-Y.)))^2
	ret1 <- (n-2) * ret0
	ret2.1 <- (sum((g/2-x.)^2)) * (sum((q-Y.)^2))
	#ret2.1 <- (sum((g-x.)^2)) * (sum((q-Y.)^2))
	ret2.2 <- ret0
	ret2 <- ret2.1 - ret2.2
	
	return(list(F=ret1/ret2,df.1 = 1, df.2 = n-2))
}

snp.F(df.case$g,df.case$acpa)

- ぐるぐる回して、Fの分布を取る

n.iter<- 10000
# outcome-dependent sampling
n.case <- 1000
n.control <- 1000

af <- 0.3
gf <- c(af^2,2*af*(1-af),(1-af)^2)
rsk <- 1
risk <- c(rsk, rsk^0.5,1)
gf.case <- gf*risk
gf.case <- gf.case/sum(gf.case)

# disease status
k <- c(rep(1,n.case),rep(0,n.control))
# genotype
g <- c(sample(2:0,n.case,prob=gf.case,replace=TRUE),sample(2:0,n.control,prob=gf,replace=TRUE))
# ACPA
q.r <- 0
acpa <- c(rnorm(n.case) + q.r*g[1:n.case],rep(NA,n.control))

df.all <- data.frame(k=k,g=g,acpa=acpa)
df.fill <- df.all
df.fill$acpa[which(is.na(df.fill$acpa))] <- min(df.case$acpa)
snp.T <- function(t){
	a <- 0:2
	m1 <- apply(t,1,sum)
	m2 <- apply(t,2,sum)
	n <- sum(t)
	ret <- sum((t[1,]-m2*m1[1]/n)^2/(m2*m1[1]/n))+ sum((t[2,]-m2*m1[2]/n)^2/(m2*m1[2]/n))
	ret
}


fs <- rep(0,n.iter)
for(ii in 1:n.iter){
	

# genotypeをシャッフル
fs[ii] <- snp.F(sample(df.fill$g),df.fill$acpa)$F

#snp.F(df.case$g,df.case$acpa)
}

plot(sort(pchisq(fs,1)))

２つの独立なpを併せて自由度４のカイ二乗値にする、その確かめ

p1<-runif(1000)
p2<-runif(1000)

k1 <- 2
k2 <- 2
g <- -k1*log(p1)-k2*log(p2)

p3 <- pchisq(g,4)

plot(sort(p3))

カイ二乗分布とガンマ分布との関係 $c\chi^2_d = \Gamma(2c,\frac{d}{2})$ の確かめ

n.iter <- 10
n.pt <- 1000
R.chisq <- c()
R.gamma <- c()
cols <- c()
for(i in 1:n.iter){
	k <- runif(1)*10
	#d <- sample(1:10,1)
	d <- runif(1)*10
	r.chisq <- rchisq(n.pt,d)*k
	r.gamma <- rgamma(n.pt,scale=2*k,shape=d/2)
	R.chisq <- c(R.chisq,sort(r.chisq))
	R.gamma <- c(R.gamma,sort(r.gamma))
	cols <- c(cols, rep(i,n.pt))
}
xlim <- ylim <- range(c(R.chisq,R.gamma))
plot(R.chisq,R.gamma,col=cols,cex=0.1,xlim=xlim,ylim=ylim)

n.pt <- 10000
p.chisq <- p.gamma <- rep(0,n.pt)
for(i in 1:n.pt){
	d <- runif(1)*10
	x <- runif(1)*50
	k <- runif(1)*10
	p.chisq[i] <- pchisq(x/k,d)
	p.gamma[i] <- pgamma(x,scale=2*k,shape=d/2)
}
plot(p.chisq,p.gamma)