Circular a posteriori projection - ryamadaの遺伝学・遺伝統計学メモ

いくつかの基準カテゴリがあって、それらに遠近関係があるという。多次元だと視覚化しにくいので、２次元のしかも円周上に配置する。その上で、たくさんあるその他もろもろの点を、その基準カテゴリに帰属する「割合」を基に、その２次元円内にプロットする方法の話
昨日の記事でも触れた血球細胞の１細胞RNAseqの論文(こちら)から、基準カテゴリを定め、全細胞をプロットしよう、という話

# 適当にデータ次元を定める
d <- 2
# 基準点の数
Ng <- 7

# 基準点の座標を遠近取り混ぜて決める
X <- matrix(rnorm(d*Ng,0,3),ncol=d)
X[2,] <- X[1,]+rnorm(d,0,0.4)
X[4,] <- X[3,] + rnorm(d,0,0.6)
X[6,] <- X[5,] + rnorm(d,0,0.1)
# 各基準点の周囲に標本点を発生させる
n.per.gr <- sample(10:1000,Ng)
Y <- matrix(0,0,d)
for(i in 1:Ng){
	tmp <- matrix(0,n.per.gr[i],d)
	for(j in 1:d){
		tmp[,j] <- rnorm(n.per.gr[i],X[i,j])
	}
	Y <- rbind(Y,tmp)
}
# その散らばり具合を見てみる
plot(as.data.frame(Y))
# 各点の基準点への帰属配分を距離を基準に決めておく。あくまでもこの配分を適当に作りたいだけ
U <- matrix(0,length(Y[,1]),Ng)
for(i in 1:length(Y[,1])){
	tmp <- rep(0,Ng)
	for(j in 1:Ng){
		tmp[j] <- 1/exp(0.5*sum((Y[i,]-X[j,])^2))
	}
	U[i,] <- tmp/sum(tmp)
}
# 帰属具合のばらつきを表示
plot(as.data.frame(U))
# この帰属具合のばらつきに応じて、基準点の距離行列を以下の要領で定める
# この部分からは論文の手法
V <- t(U) %*% U
V. <- V
for(i in 1:Ng){
	for(j in 1:Ng){
		V.[i,j] <- V[i,j]/sum(V[i,])/sum(V[,j])
	}
}
D <- exp(-1000*V.)
D
# 基準点間距離行列から、巡回セールスマン問題を実行し
# 周回性の道を作る
# create a TSP
library(TSP)
tsp <- TSP(D)
tsp

tour <- solve_TSP(tsp)

# 巡回路の道順を登録する。巡回路が「円周」となる
ord <- c()
for(i in 1:Ng){
	ord <- c(ord,tour[[i]])
}
ord <- c(ord,ord[1])
# 巡回路の順序ごとにノード間距離を取り出す
# それが円周上の距離に相当するようにする(角度が距離に比例することになる)
D.tour <- rep(0,Ng)
for(i in 1:Ng){
	D.tour[i] <- D[ord[i],ord[i+1]]
}

a <- cumsum(D.tour)/sum(D.tour)*2*pi

# これにより、基準点の円周上座標が定まる
x.gr <- cos(a)
y.gr <- sin(a)

# そのほかのすべての点は、基準点の割合と基準点の座標とを使って以下の要領で計算する
xy.all <- matrix(0,sum(n.per.gr),2)
for(i in 1:Ng){
	xy.all[,1] <- xy.all[,1] + U[,i] *x.gr[i]
	xy.all[,2] <- xy.all[,2] + U[,i] *y.gr[i]
}
t <- seq(from=0,to=2*pi,length=1000)
x <- cos(t)
y <- sin(t)

# プロットする
plot(x,y,type="l")
points(x.gr,y.gr,col=1:Ng,pch=20,cex=5)
points(xy.all,pch=20,cex=0.1)