「束」同士をペアワイズで処理する - ryamadaの遺伝学・遺伝統計学メモ

行ベクトルのペアの内積を計算する

N1<-4
N2<-3
M<-matrix(runif(N1*N2),N1,N2)
# 各行を単位ベクトルにする
# まず長さを計算する
L<-sqrt(apply(M^2,1,sum))
# 各行を各行ベクトルの長さで割って単位ベクトルにする
M2<-M/L
M2
# 単位ベクトルであることの確認をする
apply(M2^2,1,sum)
# 以下の式でペアワイズの内積が計算できる
# 単位行列なのでcos(角度)になっている
M2%*%t(M2)

行ベクトルのペアに距離を定義して階層的クラスタリングする

N1<-4
N2<-20
M<-matrix(runif(N1*N2),N1,N2)
# 4つの行ベクトルが作るペアに「距離」を定めて木を作る
dM<-dist(M)
dM
library(ape)
plot(nj(dM))
# クラスター解析
plot(hclust(dM))

固有値分解して「分離の良い軸」を抽出して、その軸に合わせてプロットする
- こちらの"書籍内Rソースコード"のファイル"R7-5.R"

#偏った集団構成(100人規模の亜集団４つと10人規模の亜集団を20個)で
#100項目のデータを作成
Nm<-100 #項目数
# 亜集団別の人数発生(100人くらいの4亜集団と20人くらいの10亜集団)
Ns<-c(rpois(4,100),rpois(20,10)) 
Npop<-length(Ns) #亜集団数
M<-NULL #全ジェノタイプデータを納める行列
#亜集団別に平均を振ってシミュレーション
for(j in 1:Npop){
 tmpM<-matrix(rep(0,Nm*Ns[j]),ncol=Nm)
 for(i in 1:Nm){ # 項目ごとのループ
  af<-rnorm(1) # 項目の頻度のおよその値
  tmpM[,i]<-rnorm(Ns[j],af) # 亜集団別の頻度
 }
 #全データ行列に格納
 M<-rbind(M,tmpM)
}
# データを標準化
wholemean<-mean(M)
M<-M-wholemean # 全平均が0になるように
mu<-apply(M,2,mean) # 列平均
M<-t(t(M)-mu) # 列平均が0になるように
# 固有値分解
svdout<-svd(M)
M2<-svdout$u%*%diag(svdout$d) # 分解後データ行列
par(mfcol=c(1,2))
# 固有値分解前後をimage()プロット
image(1:sum(Ns),1:Nm,M,xlab="サンプル(大集団→小集団)",ylab="項目")
image(1:sum(Ns),1:Nm,M2,xlab="サンプル(大集団→小集団)",
ylab="PCA後eigen項目")
df1<-as.data.frame(M);df2<-as.data.frame(M2) # データフレーム化
L<-1:5;par(mfcol=c(1,1))
plot(df1[,L]) # ５軸がなす軸ペアでサンプルをプロット。分離しない
plot(df2[,L]) # 固有値分解後に分離力のあるトップ５軸でのプロットは分離する
vM1<-apply(M,2,var) # 項目ごとの分散(分解前)
vM2<-apply(M2,2,var) # 項目ごとの分散(分解後)
ylim<-c(min(vM1,vM2),max(vM1,vM2))
# 固有値分解前の各項目の分散はどれも同じ程度だが
# 固有値分解には、分散の大きいものと小さいものとのコントラストが大きくなっている
plot(vM1,ylim=ylim,type="b")
par(new=T)
plot(vM2,ylim=ylim,type="b",col="red")