中心化後特異値分解と固有値分解 - ryamadaの遺伝学・遺伝統計学メモ

参考こちら
$X=USV^T$ が特異値分解。
変形して
。
- これを解いて、S,V->Uが得られる。
今、Xを中心化すると
$XX^T$ 　はXの分散共分散行列に比例した値になるので、中心かした特異値分解と固有値分解は、同じようなもの。

#構造化集団をシミュレート
Nm<-1000 #マーカー数
Npop<-4 #亜集団数
Ns<-c(100,150,200,250) #集団別人数
M<-NULL #全ジェノタイプデータを納める行列
#亜集団別にアレル頻度を振ってシミュレーション
for(j in 1:Npop){
 tmpM<-matrix(rep(0,Nm*Ns[j]),nrow=Nm)
 for(i in 1:Nm){
  af<-runif(1)*0.8+0.1
  f<-rnorm(1,sd=0.01)
  if(abs(f)>1) f=0
  df<-c(af^2,2*af*(1-af),(1-af)^2)
  df[1]<-df[1]+f/2*df[2]
  df[3]<-df[3]+f/2*df[2]
  df[2]<-1-df[1]-df[3]
  tmpM[i,]<-sample(c(0,1,2),Ns[j],replace=TRUE,prob=df)
 }
 #全データ行列に格納
 M<-cbind(M,tmpM)
}

##########################
##PCA (Eigenstratの方式)##
##########################


#マーカー別平均
mu<-apply(M,1,mean)
#マーカー平均で標準化
M2<-M-mu
#分散で標準化
M2<-M2/sqrt(mu/2*(1-mu/2))
#個人間の分散・共分散行列
X<-1/Nm*t(M2)%*%M2
#固有値分解
eiout<-eigen(X)
#固有値をプロット
plot(eiout$values)
#意味のある固有値の数は、亜集団の数-1
#PLoS "Population Structure and Eigenanalysis by Nick Patterson"にも記述がある通り
#plotしてみる
eivect<-as.data.frame(eiout$vectors)
eilist<-1:(Npop)
plot(eivect[,eilist])

##########################
##Mをいきなりsvdしてみる##
##########################
svdout<-svd(M)
#固有値をプロット
plot(svdout$d)
#意味のある固有値の数は、亜集団の数

plot(eiout$values,svdout$d)
#svdのトップ固有値がeigenにはない

#標準化しないときに出る、トップの固有値は何？
svdvect<-as.data.frame(svdout$v)

plot(svdvect[,eilist])


#############################
##標準化したM2をsvdしてみる##
#############################
svdout2<-svd(M2)
#固有値をプロット
plot(svdout2$d)
#意味のある固有値の数は、亜集団の数-1

plot(eiout$values,svdout2$d)
#eigenの固有値と同じ

svdvect2<-as.data.frame(svdout2$v)

plot(svdvect2[,eilist])

################################
#３つのeigenevectorの関係を見る#
################################

comparison<-as.data.frame(cbind(eivect[,eilist],svdvect[,eilist],svdvect2[,eilist]))
plot(comparison)