2011-08-13

PLSR Partial Least Square Regression

回帰オミックス chemometrics 計量化学主成分回帰特異値分解 PLSR

昨日の続き
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$
- - ここで $V$ はloadings行列と呼ばれるとともに、 $X=UDV^{-1}$ なる特異値分解を構成する行列でもある
  - また $T$ はスコア行列と呼ばれるとともに、特異値分解の構成行列を用いて $T=UD$ と書き表せる
PLSRでも同じ形を作る。そしてその構成行列を算出するのに少し手間がかかる
PCRの行列がXのみの情報を使っているのに対して、PLSRの行列は、の特異値分解から計算を始め、その後もの両方に関する残差を説明する因子を加えていくことで最終的に残差をなくす繰り返し処理をしていることに反映されているように「X,Y対称性を持った回帰」を目指したものになっている
- - PCRでの $T$ (スコア行列)に相当するものとしてPSLRにおけるスコア行列(PCRのそれと区別するために $T_s$ と書くことにする)を求めたい
- , を射影行列(projection行列)と言う
  - PCRでは $B=V(T^TT)^{-1}T^TY$ であり $V$ をloadings行列と呼んだ
  - PCRでは $T=XV$ であったが、PLSRでは $T_s=XR$ である
- ここでがほしい
  - PCRでは $V$ は $X$ の特異値分解から得られたが、PLSRでは繰り返し計算を要する
- ,は２つの行列((Xの)loadings 行列)と(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_a$ と $Y$ のみで計算できることを以下で確かめた

for(i in 1:maxa){
	tmpT<-Ts[[i]]
	print(pcr.out$fitted.values[,,i]-tmpT%*%(solve(t(tmpT)%*%tmpT)%*%t(tmpT)%*%Y))
}

- PLSRでは $T_a$ を $T_{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_a$ を $B_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)