■

source("MST.R")
par(ask=F)

Np<-100
HN<-500
df<-2000
N<-280							#バリアントの数 i
ka<-1000
coa<-c(0,1000,2000,3000,4000,5000)
lcoa<-length(coa)
syu<-matrix(200,Np,lcoa)
syu2<-matrix(0,Np,lcoa)
colnames(syu)<-coa
nrvcv<-matrix(200,Np,lcoa)


afb<-0
for (i in 1:round(df/20)){
	afd<-(i-0.5)/df
	afb<-c(afb,rep(afd,round(N*(pexp(i/df,ka)-pexp((i-1)/df,ka)))))
}
raf<-afb[-1]
#caf<-c(0.02,0.02,0.03,0.04,0.05,0.06,0.08,0.10,0.20,0.30)
#af<-c(raf,caf)
af<-raf
hist(af,breaks=seq(0,max(af)+1/df,1/df),xlab="model allele frequency",main="allele frequency model")
N<-length(af)

for (inp in 1:Np){

ir<-al<-matrix(0,HN,N)		#各ローカスにおけるリスク、持っているアレル
P<-sr<-matrix(0,HN,lcoa)			#フェノタイプ、リスクの総和
iaf<-rep(0,N)		#真のアレル頻度
risk<-matrix(0,N,lcoa)

kijun<-sum(af)
plmi<-sign(rnorm(N, mean=0, sd=1))
for (icoa in 1:lcoa){
kr<-coa[icoa]
risk[,icoa]<-rnorm(N,mean=exp(-kr*af),sd=exp(-kr*af)*0.1)*plmi
kjn<-kijun/sum(abs(risk[,icoa])*af)
risk[,icoa]<-risk[,icoa]*kjn
plot(af,risk[,icoa],xlab="model allele frequency",main="allele frequency - risk")
}

#個人のアレルの有無振り分け iさんについて
for(i in 1:HN){
	a1<-a2<-tmp1<-tmp2<-rep(0,N)
	a1<-runif(N)					#ＲＶの有無を決定
	a2<-runif(N)
	tmp1[which(a1<af)]<-1
	tmp2[which(a2<af)]<-1
	al[i,]<-tmp1+tmp2				#アレルをいくつもつか
	for (icoa in 1:lcoa){
		sr[i,icoa]<-sum(risk[,icoa]*al[i,])		#リスク総和
	}
}


iaf<-apply(al,2,mean)/2			#アレル頻度確かめ
#plot(af,iaf,xlab="model allele frequency",ylab="real allele frequency")						#グラフを描画

for(icoa in 1:lcoa) plot(density(sr[,icoa]),xlab="sum(risk)",main="sum(risk) = Phenotype",xlim=c(-4,4),ylim=c(0,8))
P<-sr
syu2[inp,]<-apply(sr,2,var)
#sdfrac<-0.5			#遺伝率
#P<-sr*sdfrac+rnorm(N,sd=sqrt(var(sr))*(1-sdfrac))
#plot(sr,P)

#----------------ＭＳＴ・量的形質----------------


asdf<-apply(al,2,sum)
qwer<-which(asdf==0)
ala<-al[,-qwer]
jaf<-apply(ala,2,mean)/2

hist(jaf,breaks=seq(0,max(jaf)+1/HN,1/HN/2),xlab="real allele frequency",main="real allele frequency")

Nrv<-length(jaf[which(jaf<0.05)])
Ncv<-length(jaf[which(jaf>=0.05)])
for (icoa in 1:lcoa){
nrvcv[inp,icoa]<-Nrv
X<-cbind(P[,icoa],ala)
MST<-StatsMST(X,perm=TRUE,Nperm=100,tobeshuffled=c(1),dist.method="euclid")
plot(sort(MST$PermL),ylim=range(c(MST$L,MST$PermL)))
abline(h=MST$L,col=2)
title(paste(as.character(coa[icoa]),as.character(inp),as.character(Nrv),as.character(Ncv)))
ll<-MST$L
cat1<-seq(from=0,to=1,by=0.01)
cat2<-quantile(MST$PermL,cat1)
cat3<-cat2-ll
syu[inp,icoa]<-length(cat3[cat3<0])
}
}

boxplot(syu,xlab="kr (risk=rsign*exp(-kr*allelefrequency))")
syu

save.image("10051.RData")

なるべく、リスク総和の分布が似た形になると、検出力がリスクの与え方によって変わるのか、分布の仕方によって変わるのかの区別をつけられる、と考えた。

レアであるほどリスクが高くなり、その高くなる割合を変える、という前提は変えられないので、リスク総和が小さくなりがちなリスクの与え方になっているとき、リスクの値自体を増やすことにする。

sum(risk * allele frequency)の値が一緒になるようにすると、リスク総和の分布が近くなる。