2011-08-31

ちょっとメモ

# ケースとコントロールの人数
Nca<-Nco<-1000
# ケースコントロール別のレアアレル保有者数
nca<-20
nco<-0

# 計算するための、ケース・コントロールのアレル頻度(パラメタとして振ってある)
dp<-0.001
pca<-seq(from=dp,to=0.05,length=100)
pco<-pca

# 確率計算、アレル頻度比計算
pcaco<-expand.grid(pca,pco)
L<-rep(0,length(pcaco[,1]))
R<-L
for(i in 1:length(pcaco[,1])){
	tmpca<-pcaco[i,1]
	tmpco<-pcaco[i,2]
	R[i]<-tmpca/tmpco
	L[i]<-lchoose(Nca,nca)*lchoose(Nco,nco)*nca*log(tmpca)+nco*log(tmpco)+(Nca-nca)*log(1-tmpca)+(Nco-nco)*log(1-tmpco)
}

Ans<-cbind(pcaco,L,R)
SortedAns<-Ans[order(Ans[,4]),]

# GWASレベルでの、Effective number of independent tests
a<-0.01
pt<-10^(-8)
Nm<-log(1-a)/log(1-pt)

#Lm<-1-(1-exp(L))^Nm
Lm<-1-exp(Nm*log(1-exp(L)))

# 同じ計算をプロットするために行列型オブジェクトで再計算
M<-matrix(0,length(pca),length(pco))
for(i in 1:length(pca)){
	for(j in 1:length(pco)){
		M[i,j]<-lchoose(Nca,nca)+lchoose(Nco,nco)+nca*log(pca[i])+nco*log(pco[j])+(Nca-nca)*log(1-pca[i])+(Nco-nco)*log(1-pco[j])
	}
}

#image(exp(M))
# ケース・コントロールのアレル頻度の比の逆数
as<-seq(from=1,to=3,by=0.25)

# マルチプルテスティング補正なし
contour(exp(M))
for(i in 1:length(as)){
	abline(0,1/as[i],col=i)

}

# マルチプルテスティング補正した
Mm<-1-exp(Nm*log(1-exp(M)))

#image(Mm)
contour(Mm)

for(i in 1:length(as)){
	abline(0,1/as[i],col=i)

}

# これはフィッシャーの正確確率検定
x<-c(rep(1,nca),rep(2,Nca-nca),rep(1,nco),rep(2,Nco-nco))
y<-c(rep(1,Nca),rep(2,Nco))
fisher.test(x,y)

N<-100
ns<-0:20
ps<-matrix(0,length(ns),length(ns))
for(i in 1:length(ns)){
	for(j in 1:length(ns)){
		
	n<-ns[i]
if(ns[i]*ns[j]!=0){
x<-c(rep(1,n),rep(2,N-n),rep(1,ns[j]),rep(2,N-ns[j]))
y<-c(rep(1,N),rep(2,N))
ps[i,j]<-fisher.test(x,y)$p
}
}
}
matplot(log(ps),type="l")

Ns<-2^(5:14)
n<-20
Ps<-rep(0,length(Ns))
for(i in 1:length(Ns)){
x<-c(rep(1,n),rep(2,Ns[i]-n),rep(1,0),rep(2,Ns[i]))
y<-c(rep(1,Ns[i]),rep(2,Ns[i]))
Ps[i]<-fisher.test(x,y)$p
}

plot(Ns,Ps)


Nperm<-1000

ps<-rep(0,Nperm)

for(i in 1:Nperm){
	
	ps[i]<-fisher.test(x,y[sample(1:length(y))])$p
}

plot(sort(ps))


x<-c(rep(1,20),rep(2,980),rep(1,0),rep(2,1000))
y<-c(rep(1,1000),rep(2,1000))

fisher.test(x,y)

Nca<-Nco<-1000
nca<-10
nco<-0

dp<-0.001
pca<-seq(from=dp,to=0.1,length=100)
pco<-pca

pcaco<-expand.grid(pca,pco)
L<-rep(0,length(pcaco[,1]))
R<-L
for(i in 1:length(pcaco[,1])){
	tmpca<-pcaco[i,1]
	tmpco<-pcaco[i,2]
	R[i]<-tmpca/tmpco
	L[i]<-nca*log(tmpca)+nco*log(tmpco)+(Nca-nca)*log(1-tmpca)+(Nco-nco)*log(1-tmpco)
}

Ans<-cbind(pcaco,L,R)
SortedAns<-Ans[order(Ans[,4]),]




M<-matrix(0,length(pca),length(pco))
for(i in 1:length(pca)){
	for(j in 1:length(pco)){
		M[i,j]<-nca*log(pca[i])+nco*log(pco[j])+(Nca-nca)*log(1-pca[i])+(Nco-nco)*log(1-pco[j])
	}
}

image(M)

Effective No. tests

N<-log(0.99)/log(1-10^(-8))

Nca<-Nco<-1000
nca<-20
nco<-1


dp<-0.001
pca<-seq(from=dp,to=0.05,length=100)
pco<-pca

pcaco<-expand.grid(pca,pco)
L<-rep(0,length(pcaco[,1]))
R<-L
for(i in 1:length(pcaco[,1])){
	tmpca<-pcaco[i,1]
	tmpco<-pcaco[i,2]
	R[i]<-tmpca/tmpco
	L[i]<-lchoose(Nca,nca)*lchoose(Nco,nco)*nca*log(tmpca)+nco*log(tmpco)+(Nca-nca)*log(1-tmpca)+(Nco-nco)*log(1-tmpco)
}

Ans<-cbind(pcaco,L,R)
SortedAns<-Ans[order(Ans[,4]),]

a<-0.01
pt<-10^(-8)
Nm<-log(1-a)/log(1-pt)

#Lm<-1-(1-exp(L))^Nm
Lm<-1-exp(Nm*log(1-exp(L)))


M<-matrix(0,length(pca),length(pco))
for(i in 1:length(pca)){
	for(j in 1:length(pco)){
		M[i,j]<-lchoose(Nca,nca)+lchoose(Nco,nco)+nca*log(pca[i])+nco*log(pco[j])+(Nca-nca)*log(1-pca[i])+(Nco-nco)*log(1-pco[j])
	}
}

#image(exp(M))
# ケース・コントロールのアレル頻度の比の逆数
as<-seq(from=1,to=3,by=0.25)

contour(exp(M))
for(i in 1:length(as)){
	abline(0,1/as[i],col=i)

}

Mm<-1-exp(Nm*log(1-exp(M)))

#image(Mm)
contour(Mm)

for(i in 1:length(as)){
	abline(0,1/as[i],col=i)

}
x<-c(rep(1,nca),rep(2,Nca-nca),rep(1,nco),rep(2,Nco-nco))
y<-c(rep(1,Nca),rep(2,Nco))
fisher.test(x,y)