分散が集約されて共分散が小さくなる - ryamadaの遺伝学・遺伝統計学メモ

ryamada本のR7-5.Rでは、固有値分解を実施している
固有値分解では、オリジナルの軸が説明する分散が10人並みなのに対して、固有値分解を施すことにより、説明する分散が大きい方から軸をとっている
その図がこれだが

そのようにして軸を取り直すことで、分散共分散行列が
- 対角成分(分散)の和は変わらないが、大きな値を持つものと、ほとんど値を持たないものとに分けなおして、大きい順に並べていること
- 非対角成分(共分散)が小さくなっていることが見て取れる

#偏った集団構成(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=TRUE)
plot(vM2,ylim=ylim,type="b",col="red")
library(rgl)
cD<-cov(M) # 分解前の行列の分散共分散行列
cD2<-cov(M2) # 分解後の行列の分散共分散行列
# 対角成分の和(分散の和)は等しい
sum(diag(cD)) 
sum(diag(cD2))
# 非対角成分の和(共分散の和)は違う
sum(cD)-sum(diag(cD))
sum(cD2)-sum(diag(cD2))

open3d()
plot3d(row(cD),col(cD),cD,type="l")
open3d()
plot3d(row(cD2),col(cD2),cD2,type="l")

par(mfcol=c(1,2))
persp(cov(M))
persp(cov(M2))
par(mfcol=c(1,1))