乱点

• 格子が交互に充てんされる模様の２次元版が市松模様
• k次元一様乱点を発生させ、そのうち、0を含む偶数要素の値を負にすると、k次元市松模様乱点が発生できる
• ４次元以上にすると、ただの一様分布だけれど、３次元までなら、「市松模様」を「見る」ことができる
• さて、k=2のときは、明らかに、第１軸と第２軸とに「関係」が見えるけれども、k>2のときには、どの２軸を選んでも、その２軸の間に「関係」が見えない
• しかしながら、明らかに、２軸を考慮すると「関係〜ルール」がある
• これは、SNPのハプロタイプで言えば、３ＳＮＰの作る８ハプロタイプの頻度が、(ABC,ABc,AbC,Abc,aBC,aBc,abC,abc)=(0.25,0,0,0.25,0,0.25,0.25,0)となっているときに、A/a,B/b,C/cには連鎖不平衡があるのに、(A/a,B/b)で見るとそれが連鎖平衡に見える((B/b,C/c)で見ても、(C/c,A/a)で見ても同じこと)のと同じ話
• この３ＳＮＰ〜３軸の「関係〜ルール」の存在をp値にすることを考える
• 最小全域木を使って「見える」ようにしてみる(後半のソース)

# 多次元市松模様乱点発生

IchimatsuK<-function(n,k){
	r<-matrix(runif(n*k),n,k)
	loop<-TRUE
	ok<-0
	tmps<-NULL
	while(loop){
		s<-matrix(sample(c(0,-1),n*k*3,replace=TRUE),ncol=k)
		tmps<-rbind(tmps,s[which(apply(s,1,sum)%%2==0),])
		if(length(tmps[,1])>=n){
			loop<-FALSE
		}
	}
	tmps<-tmps[1:n,]
	r*((tmps+0.5)*2)
}
n<-10000
k<-4
IchimatsuOut<-IchimatsuK(n,k)
library(rgl)
plot3d(IchimatsuOut[,1:3])

plot(as.data.frame(IchimatsuOut))
n<-10000
k<-3
IchimatsuOut<-IchimatsuK(n,k)
plot3d(IchimatsuOut[,1:3])

plot(as.data.frame(IchimatsuOut))

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))
	}
}
# Mの列sをシャッフルして返す
ShuffleForPerm<-function(M,s=NULL){
	if(is.null(s)){
		s<-2:length(M[1,])
	}
	sM<-matrix(1:length(M),nrow(M),ncol(M))
	sM<-apply(sM[,s],2,sample)
	sM2<-M
	sM2[,s]<-matrix(M[c(sM)],nrow(M),length(s))
	sM2
}
n<-100
k<-4
IchimatsuOut<-IchimatsuK(n,k)
plot3d(IchimatsuOut[,1:3])


sMST<-StatsMST(IchimatsuOut,Nperm=100)
plot(sort(sMST$PermL),ylim=range(c(sMST$L,sMST$PermL)))
abline(h=sMST$L,col=2)
print(length(which(sMST$PermL<sMST$L))/sMST$Nperm)