■

par(ask=TRUE)
N<-200							#バリアントの数 i
HN<-100000						#人数 j
ir<-al<-matrix(0,HN,N)			#各ローカスにおけるリスク、持っているアレル
wei<-wei2<-P<-sr<-rep(0,HN)	#重み、フェノタイプ、リスクの総和
risk<-af<-iaf<-rep(0,N)		#真のアレル頻度

Np<-50							#MST何回か
KN<-1							#何種類のKについて
knnu<-rep(0,KN)					#何種類のKについて
IDN<-10							#何種類の遺伝率か
idnnu<-rep(0,IDN)				#何種類の遺伝率か
syukei<-matrix(0,Np,KN)		#boxplot用
syukei2<-matrix(0,Np,KN)		#boxplot用

syukei2<-matrix(0,Np,IDN)		#boxplot用



thress<-seq(from=0.9,to=0.98,by=0.005)	#0.005刻みのthreshold

aa<-60.0						#アレル頻度パラメータ
af<-rexp(N,aa)					#アレル頻度を与える

for(knp in 1:KN){

mm<-60							#リスク分布のmeanのパラメータ
#knnu[knp]<-mm
risk<-exp(-mm*af)*sign(rnorm(N, mean=0, sd=1))


#risk[which(af<0.005)]<-200


#hist(af)
#hist(risk)

cccc<-as.character(knp)
ssss<-paste("manaf",cccc,".png")
png(ssss)          #デバイスドライバの用意
plot(af,risk)                      #グラフを描画
dev.off()

#個人のアレルの有無振り分け 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				#アレルをいくつもつか
	sr[i]<-sum(risk*al[i,])		#リスク総和
}


#for (i in 1:N) iaf[i]<-mean(al[,i])/2	#アレル頻度確かめ
#plot(af,iaf)
#plot(density(sr))

for(iii in 1:IDN){
sdfrac<-0.98-((iii-1)*0.1)			#遺伝率
idnnu[iii]<-sdfrac
P<-sr*sdfrac+rnorm(N,sd=sqrt(var(sr))*(1-sdfrac))    
plot(sr,P)

#----------------ＭＳＴ・量的形質----------------
par(ask=FALSE)

Ncaco<-200
for (counte in 1:Np){
	Cmst<-sample(1:HN,Ncaco,replace=FALSE)
	X<-cbind(P[c(Cmst)],al[c(Cmst),])
#	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(counte+1000*knp)
#	ll<-MST$L
#	cat1<-seq(from=0,to=1,by=0.01)
#	cat2<-quantile(MST$PermL,cat1)
#	cat3<-cat2-ll
#	syukei[counte,knp]<-length(cat3[cat3<0])
	
	MST<-StatsMST(X,perm=TRUE,Nperm=100,tobeshuffled=c(1),dist.method="manhattan")
	plot(sort(MST$PermL),ylim=range(c(MST$L,MST$PermL)))
	abline(h=MST$L,col=2)
	title(counte+1000*knp+1000000)
	ll<-MST$L
	cat1<-seq(from=0,to=1,by=0.01)
	cat2<-quantile(MST$PermL,cat1)
	cat3<-cat2-ll
	syukei2[counte,knp]<-length(cat3[cat3<0])

}
}
}
colnames(syukei2) <- idnnu	#knnu
boxplot(syukei2)

上手く回るかまだ試してないので、これは家でまわそう。多分25分くらいかかる。

#Mの列sをシャッフルして返す
ShuffleForPerm<-function(M,s=NULL){
	if(is.null(s)){				#s=NULLなら（つまりデフォでは）２列目以降をシャッフル
		s<-2:length(M[1,])
	}
	sM<-matrix(1:length(M),nrow(M),ncol(M))	#sMを列ごとにシャッフルしてそのインデックスのを取り出す感じ
	sM<-apply(as.matrix(sM[,s]),2,sample)
	sM2<-M
	sM2[,s]<-matrix(M[sM],nrow(M),length(s))
	sM2
}

library(vegan)

#MSTの値について、検定
StatsMST<-function(X,perm=TRUE,tobeshuffled=NULL,Nperm=1000,dist.method="euclid"){
	library(vegan)
	Nnode<-length(X[,1])
	distM<-dist(X,method=dist.method)

	# calculate stats
	MST<-spantree(distM)
	
	if(!perm){
		return(list(Nnode=Nnode,distM=distM,MST=MST,L=sum(MST[[2]]),PermL=NULL,perm=perm,Nperm=Nperm,dist.method=dist.method))
	}else{
		if(is.null(tobeshuffled)){
			tobeshuffled<-2:length(X[1,])
		}
		PermL<-rep(0,Nperm)

		for(i in 1:Nperm){
			tmpdistM<-dist(ShuffleForPerm(X,tobeshuffled),method=dist.method)
			tmpMST<-spantree(tmpdistM)
			PermL[i]<-sum(tmpMST[[2]])
		}
		return(list(Nnode=Nnode,distM=distM,MST=MST,L=sum(MST[[2]]),PermL=PermL,perm=perm,Nperm=Nperm,dist.method=dist.method))
	}
}

#MSTの値について、プロット
plot.MST.L<-function(mst){
	ylim<-range(c(mst$L,mst$PermL))
	plot(sort(mst$PermL),ylim=ylim)
	abline(h=mst$L)
}

読み込み用