PLSR Partial Least Square Regression

  • 昨日の続き
  • PLSモデルはパッケージparcorでも
    • lasso, adaptive lasso, PLS, and Ridge Regression, model selection for lasso, adaptive lasso and Ridge regression based on cross-validation.
  • 重回帰→主成分回帰→PLSRという展開
  • Y=XB+\epsilon
  • 重回帰
    • Y=X((X^TX)^{-1}X^TY)+\epsilon=XB+\epsilon
    • B=(X^TX)^{-1}X^TY
  • 主成分回帰
    • Y=T((T^TT)^{-1}T^TY)+\epsilon=TC+\epsilon
    • B=V(T^TT)^{-1}T^TY
      • ここでVはloadings行列と呼ばれるとともに、X=UDV^{-1}なる特異値分解を構成する行列でもある
      • またTはスコア行列と呼ばれるとともに、特異値分解の構成行列を用いてT=UDと書き表せる
  • PLSRでも同じ形を作る。そしてその構成行列を算出するのに少し手間がかかる
  • PCRの行列がXのみの情報を使っているのに対して、PLSRの行列は、X^TY特異値分解から計算を始め、その後もX,Yの両方に関する残差を説明する因子を加えていくことで最終的に残差をなくす繰り返し処理をしていることに反映されているように「X,Y対称性を持った回帰」を目指したものになっている
    • Y=T_s((T_s^TT_s)^{-1}T_s^TY)+\epsilon
      • PCRでのT(スコア行列)に相当するものとしてPSLRにおけるスコア行列(PCRのそれと区別するためにT_sと書くことにする)を求めたい
    • B=R(T_s^TT_s)^{-1}T_s^TY, Rを射影行列(projection行列)と言う
      • PCRではB=V(T^TT)^{-1}T^TYでありVをloadings行列と呼んだ
      • PCRではT=XVであったが、PLSRではT_s=XRである
    • ここでRがほしい
      • PCRではVX特異値分解から得られたが、PLSRでは繰り返し計算を要する
    • R=W(P^TW)^{-1},Rは2つの行列P((Xの)loadings 行列)とW(weight loadings 行列)とから計算される
      • これらのW,Pが繰り返し計算によって得られるのでRもそうなり、結局T_sもそうなる
    • Rのplsパッケージのplsr()関数の出力では、これらのB,T_s,P,W,Rはそれぞれ、"coefficients","scores","loadings","loading.weights","projection"と呼ばれ、リストの第1,2,3,4,7番要素として出力される
    • それを確かめてみる
  • 式を再掲する
    • Y=XB+\epsilon
    • 重回帰
      • Y=X((X^TX)^{-1}X^TY)+\epsilon=XB+\epsilon
    • 主成分回帰
      • Y=T((T^TT)^{-1}T^TY)+\epsilon=TC+\epsilon
    • PLSR
      • Y=T_s((T_s^TT)^{-1}T_s^TY + \epsilon
    • PCRでは、回帰がT_aYのみで計算できることを以下で確かめた
for(i in 1:maxa){
	tmpT<-Ts[[i]]
	print(pcr.out$fitted.values[,,i]-tmpT%*%(solve(t(tmpT)%*%tmpT)%*%t(tmpT)%*%Y))
}
    • PLSRではT_aT_{s,a}で置き換えよう
pls.out<-plsr(Y~X)
for(i in 1:maxa){
	tmpT<-pls.out$scores[,1:i]
	print(pls.out$fitted.values[,,i]-tmpT%*%(solve(t(tmpT)%*%tmpT)%*%t(tmpT)%*%Y))
}
    • 同様にPCRでは、回帰係数行列B_aB_a=V_a(T_a^TT)^{-1}T_a^T Yであることを用いて次のソースのように計算した
for(i in 1:maxa){
	tmpT<-Ts[[i]]
	print(pcr.out$coefficients[,,i]-pcr.out$loadings[,1:i]%*%(solve(t(tmpT)%*%tmpT)%*%t(tmpT)%*%Y))
}
    • PSLRではB_a=R_a(T_{s,a}^TT_{s,a})^{-1}T_{s,a}^TYとして計算しよう
for(i in 1:maxa){
	tmpT<-pls.out$scores[,1:i]
	print(pls.out$coefficients[,,i]-pls.out$projection[,1:i]%*%(solve(t(tmpT)%*%tmpT)%*%t(tmpT)%*%Y))
}
    • PCRにおけるX_{st}V=TとPLSRにおけるX_{st}R=T_sとも確かめておく
pcr.out$scores-X.st%*%pcr.out$loadings
pls.out$scores-X.st%*%pls.out$projection
    • PLSR特有の関係を次に確かめておく
pls.out$projection-pls.out$loading.weights%*%solve(t(pls.out$loadings)%*%pls.out$loading.weights)