読書日記

« 読了: Szepannek (2018) RのclustMixTypeパッケージ | メイン | 読了：「ゆりあ先生の赤い糸」「おとなになっても」「セクシー田中さん」「邦キチ！映子さん」「よちよち文藝部」「百姓貴族」「不浄を拭う人々」 »

2019年12月30日 (月)

Burnhan, A.J., MacGregor, J.F., Viveros, R. (2001) Interpretation of regression coefficients under a latent variable regression model. Journal of Chemometrics, 15, 265-284.
　仕事の都合で、PLS回帰の係数の解釈について説明する必要が生じ、あれこれ探し回って見つけた論文。
　PLS回帰係数の解釈についての論文は珍しく、やったぜついに見つけたぜと喜んだんだけど、いざ読んでみると、疑問がいっぱい浮かんできて...

　いわく。
　主成分回帰(PCR)、PLS回帰といった潜在変数重回帰(LVMR)モデルでは、
　$y = t Q + f$
　$x = t P + e$
と考える。目的変数$y$は$m$次元, 説明変数$x$は$k$次元、潜在変数$t$は$a \leq k$次元とする。
　PCR, PLSでは
　$\hat{t} = x W$
とする。[$\hat{y} = \hat{t} \hat{Q}$に代入して$\hat{y} = x W \hat{Q}$だから] $\hat{B} = W \hat{Q}$と書けば
　$\hat{y} = x \hat{B}$
と、標準的な回帰モデルと同じ形になる。つまり$B$は"回帰係数"である。
　しかし、もともとのLVMRモデルに$B$が明示的に出ているわけではないので、これをどう定義し、どう解釈するかがはっきりしない。
　予測子のランクを減らしている回帰モデルにおいて回帰係数をどう解釈するかという問題は、我々の知る限り、これまで検討されてこなかった。本論文ではLVMRモデルの回帰係数の解釈について考える。

　まずはフルランクの回帰モデルについて考える。3つのモデルについて考えます。
　モデル1、標準的な回帰モデル。
　$y = xB_1 + f$
　モデル2、予測子を確率変数とした標準的回帰モデル。
　$E(y|x) = xB_2, \ \ V(y|x) = \Sigma_f$
　モデル3、予測子が誤差を含む回帰モデル。
　$E(y|u) = E(x|u) B_3, \ \ V(y|u) = \Sigma_f, \ \ x = u + e$

　以下では次のように定義します。
　$E(y|x) = x B^\circ$
　$E(y|u) = E(x|u) B^*$
モデル2の下では$B_2 = B^\circ = B^*$である。モデル3の下では$B_3 = B^*$で、$u, e, f$が互いに独立に平均0, 共分散$\Sigma_u, \Sigma_e, \Sigma_f$のMVNに従うとすれば、
　$\Sigma_u (\Sigma_u + \Sigma_e)^{-1} B_3 = B^\circ$
である。
　[恥ずかしながらなぜこうなるのか、考えたんだけどわからない。モデル3は
　$y = u B_3 + f$
だという理解で合っているだろうか？ すると
　$V(y) = \Sigma_u B_3 + \Sigma_f$
だよね？ いっぽう$B^\circ$についていえば
　$V(y) = V(y|x) + V(x)B^\circ = V(y|x) + (\Sigma_e + \Sigma_u) B^\circ$
だと思うんですよね。だから
　$(\Sigma_e + \Sigma_u)^{-1} (\Sigma_u B_3 + \Sigma_f - V(y|x)) = B^\circ$
じゃないかなと。この左辺が$\Sigma_u (\Sigma_u + \Sigma_e)^{-1} B_3$に化ける理由がわからない。嗚呼...]

　さあ、いよいよLVMRモデル
　$y = t Q + f$
　$x = t P + e$
について考えよう。以下では$e, f, t$がそれぞれ平均0のMVNに従い互いに独立だとする。

　まず$B^\circ$について。$x$と$y$の同時分布は
　$\displaystyle \Sigma_{xy} = \left[ \begin{array}{cc} \Sigma_e + P^T \Sigma_t P & P^T \Sigma_t Q \\ Q^T \Sigma_t P & \Sigma_f + Q^T \Sigma_t Q \end{array} \right]$
　これより
　$E(y|x) = x (\Sigma_e + P^T \Sigma_t P)^{-1} P^T \Sigma_t Q$
つまり$B^\circ = (\Sigma_e + P^T \Sigma_t P)^{-1} P^T \Sigma_t Q$となるわけである。

　次に$B^*$について。もともと$E(y|u) = E(x|u) B^*$と定義してたわけだけど、ここでは$u = tP$だから、$E(y|t) = E(x|t) B^*$である。ってことは$t Q = tPB^*$である。これが任意の$t$について成り立つんだから$Q = PB^*$である。これを解くと[...本文に説明が載ってないし考えてもよくわかんないんだけど...]
　$B^* = GQ + (GP-I) Z$
ただし$G$は一般化逆行列で$PGP =P$であり、$Z$は任意の[←えっ？] $k \times m$実数行列である。このように、LVMRモデルでは$B^*$は一意に決まらない。
　[一生懸命説明してくれてんのにすまんけど、いまいちぴんとこない...]

　わかりにくいんで具体例で示そう[←おう、そうしてくれ]。
　$a=1, m=1, k=2$とし$\Sigma_e$は対角とする。
　$P$はいまや列ベクトル、その要素を$p_1, p_2$とする。$Q$はいまやスカラー$q$。$B$はいまや列ベクトル、その要素を$b_1, b_2$とする。潜在変数は$t \sim N(0, \sigma^2_t)$。$e_1, e_2, f$はiid正規で、平均は0, SDはそれぞれ$\sigma_1, \sigma_2, \sigma_f$。もはや$t, p, q$はスケーリングしているだけなので、$\sigma_t = 1$とする。
　[はいはい。要するに
　$y = t q + f, \ \ f \sim N(0, \sigma_f^2)$
　$x_1 = t p_1 + e_1, \ \ e_1 \sim N(0, \sigma_1^2)$
　$x_2 = t p_2 + e_2, \ \ e_2 \sim N(0, \sigma_2^2)$
　$t \sim N(0, 1)$
であるとき、$x$と$E(y|x)$の比$b^\circ$はどうなるの、潜在変数$u$の下で条件付けた$E(y|u)$と$E(x|u)$の比$b^*$はどうなるの、ってことね]
　このとき、
　$\displaystyle b^\circ_j = \frac{qp_j}{\sigma^2_j(p^2_1/\sigma^2_1 + p^2_2/\sigma^2_2 + 1)}$
であり、$b^*$は
　$b^*_1 p_1 + b^*_2 p_2 = q$
を満たす任意のベクトルである。[←ああ、なるほど...]

　では本題。PLS回帰の$\hat{b}$の分布について考えよう。さきほどと同様、$a=1, m=1, k=2$とする。
　1因子PLSでは$w = X^T y$なので、
　$\hat{b}^{PLS} = X^T y (y^T XX^T y)(y^T XX^T XX^T y)^{-1}$
となる[ついていけないけど、信じます...]。
　上の式は、ベクトル$X^T y (y^T XX^T y)$をスカラー$(y^T XX^T XX^T y)^{-1}$で割っている。比の期待値を期待値の比で近似するとして、$E(\hat{b}^{PLS})$は [...$p,q,n,\sigma_t,\sigma_f,\sigma_1,\sigma_2$を使った近似式が示されている。複雑な式ではないけれど、5行くらいの長い式なのでメモ省略]
　この式をよくみると、$\hat{b}^{PLS}$の要素の相対的サイズは、負荷$p_j$の相対的サイズだけでほぼ決まる。つまり
　$E(\hat{b}^{PLS}_1 / \hat{b}^{PLS}_2) \approx p_1/p_2$
　また、$\sigma_t$が大きいときには
　$E(\hat{b}^{PLS}) \approx p^T (q / pp^T)$
となり、$b^*$の性質$b^*_1 p_1 + b^*_2 p_2 = q$を満たす。つまり、$b^*$の特定の事例を近似的に不偏推定している。
　$n$が大きいときには、
　$\displaystyle E(\hat{b}^{PLS}) \approx p^T \frac{q (pp^T) \sigma^2_t}{(pp^T)^2 \sigma^2_t + p \Sigma_e p^T}$
となる($\Sigma_e$ってのは$\sigma_1, \sigma_2$を持つ対角行列)。分母に$p\Sigma_e p^T = p^2_1 \sigma^2_1 + p^2_2 \sigma^2_2$がある。つまり、誤差が大きいほど、$E(\hat{b}^{PLS}) $は$b^*$の線を滑るようにして小さくなっていく。

　これをOLS回帰と比較してみると [...面食らったんだけど、上と同じデータ生成メカニズムを仮定して、ということだと思う]
　$E(\hat{b}^{OLS}) = b^\circ$
　$E(\hat{b}^{OLS}_1 / \hat{b}^{PLS}_2) = (p_1 \sigma^2_2) / (p_2 \sigma^2_1)$
であることが示せる。

　ここからはシミュレーション。
　その1。$p_1 = p_2 = 2, q=1, \sigma_2 = 2, \sigma_f = 0.15, \sigma_t =1$に固定し, $\sigma_1$のみを4水準に動かす。つまり$x_1$のSN比だけを動かすわけである。$b^*$は$2b^*_1 + 2b^*_2 = 1$のままだが$b^\circ$は動く。[$n$が書いてない... 見落としたのだろうか...]
　$\hat{b}^{PLS}$と $\hat{b}^{OLS}$は$\sigma_1$が2から離れるにつれて離れる。$OLS$では、$\sigma_1$が小さいときに$\hat{b}^{OLS}_1$が大きくなるわけだが、PLSではそうならない。OLSは常に$b^\circ$を不偏推定しているがPLSはそうでない。$b^*$という観点から見るとどちらにもバイアスがある[それは$\sigma_t$がそんなに大きくないからだよね...]。
　$y$の予測性能を比べると、$\sigma_1$が$\sigma_2$から離れるにつれてPLSの性能が落ちる。誤差分散が違うのに回帰係数を調整しないからである。

　その2。$p_1 = 1, q=1, \sigma_1 = 0.5, \sigma_f = 0.15, \sigma_t = 1$に固定し、$p_2$と$\sigma_2$をSN比が変わらないようにして引き上げる[つまり$p_2$を4倍したら$\sigma_2$は2倍する]。$b^\circ$は動かないが$b^*$の線は動く。
　$\hat{b}^{PLS}$と $\hat{b}^{OLS}$は$p_2$が大きくなるにつれて離れる。こんどはOLSが$\hat{b}_1 = \hat{b}_2$の線上から動かず[SN比が同じだからね]、PLSがアンバランスになっていく[そりゃそうだ、$p_2$が大きいんだから]。予測性能はほとんど変わらない。[ああ、そりゃそうだな... 結局$X$のSN比は変わってないわけだから...]

　考察。
　結局、LVMRモデルの$\hat{B}$が推定しているのはなんだろうか。その答えは研究の目的によって異なる。
　Myers(1986, 回帰分析の教科書らしい)いわく、回帰分析の目的は4つある。順にみていこう。
　目的1、予測。分析者はデータ生成システムに関心を持たず、ただ$X$のなかの誤差分散のサイズに関心を持つ。この場合、回帰係数の定義には$B^\circ$がふさわしい。$B^\circ$は$X$の誤差分散を含んでいるから。
　目的2、モデルの特定、つまりシステムの説明。この場合は$B^*$がふさわしい。ある変数の誤差分散が大きいからといってその変数をディスカウントするのはおかしいから。
　目的3、変数選択。この場合はなぜ変数を選択したいかが問題になる。予測のためなら$B^\circ$が、以後の研究のためなら$B^*$がふさわしい。
　目的4、パラメータ推定。この場合、そもそもLVMRモデルには$B$なんてないということに注意。$B^\circ$の立場からみればそれは$e, f$の分布に影響されるし、$B^*$の立場からみればそれは一意に決まらない。

　Brownという人いわく、「多くの場合、$Y$の変動の一部はなんらかの予測子と連関している可能性がある。極端な例として$X_2 = cX_1$としよう。この場合、$b_1 + b_2c$だけが問題であって、$b_1$それ自体には意味がない」。
　Cramerいわく、「$x_1$と$x_2$に高い相関があるときにはどちらかだけ使え、というアドバイスは、予測の状況では正しいが、理解という目的のためには正しくない」。
　Hultquistという人いわく、「回帰係数$b_i$は単独では無意味だ。それは$Y$の予測において、collectivelyに役割を果たす」。
　こういう引用において「無意味」と呼ばれている事態は、データがreduced-rankな構造を持っていることを直接にモデリングしなかったことの結果である。LVMRモデルは回帰係数を$B^*$として定義できる。それは$x,y,t$の線形関係という観点から有意味である。[←でも先生、それって一意には決まんないんでしょ....]
　云々。

　。。。いやー、自分なりに真剣に読んでみたんだけど、やはり学力が追い付かない内容であった。ケモメトリクスについて知らないから問題設定がわからないんだよな、と自分を慰めつつ...

　恥を恐れず、この論文のメッセージを要約すると、こういうことなんじゃないかと思う。
　PLS回帰を使おうとしているそこのあなた。あなたは、$X$には誤差があり、真の予測子は少数の潜在変数だと考えているわけだよね。あなたが推定したいのがなんなのかを考えるべきだ。$X$と$E(Y|X)$の関係 $B^\circ$を推定したいのか、それとも潜在変数$T$を固定したときの$E(X|T)$と$E(Y|T)$の関係 $B^*$を推定したいのか。あなたがやりたいことが予測なら、推定したいのは$B^\circ$であり、$X$の誤差分散に応じて回帰係数を割り引いてくれるOLS回帰がお勧め。あなたがやりたいことが説明なら、推定したいのは$B^*$であり、$X$の誤差分散がどうであれ回帰係数を泰然と維持してくれるPLS回帰がお勧め。

　ううううむ。わからん。
　回帰分析を説明という目的で行っているユーザが、きっとデータ生成構造は
　$y = t Q + f$
　$x = t P + e$
だ、よおしPLS回帰を使うぞ、と思うことがあるんだろうか？ 仮にそうであれば、たとえばPLS因子数は実質科学的観点から決めるべきで、(通常行われているように)CVで決めるのは筋が通らないと思う。それに、もしそういうユーザがいたら、その人は負荷行列$P,Q$にのみ関心を持ち、$B^*$になんて関心を持たないのではないだろうか？
　むしろ、回帰分析を説明という目的で行っているユーザがわざわざPLS回帰を使う、もっとメジャーな理由はこうじゃなかろうか。「私ですか？ いえ、別にreduced-rankなデータ生成構造とか$x$の誤差とか考えてません。単に$E(y|x) = xB$だと考えてます。私は$B=B^\circ = B^*$を推定したいんです。でもほら、OLSだとたしかに$\hat{B}$は$B$の一致推定量だけど、手元のデータはサイズが小さいので(or 説明変数間の相関が高いので)、$MSE(\hat{B})$が大きくなりそうなんですよね。ですからPLSを選びました」。
　こういう無節操なユーザに対して、いや流石にお目が高い、PLSの$\hat{B}$は($B^\circ$ではなく)$B^*$の特定の事例を近似的に不偏推定する推定量であり、つまり$X$の誤差分散の影響を受けませんからね...などと説明しても、ぽかーんとされるのではなかろうか。むしろ、分析者の推定対象は$B^\circ = B^*$であると考えて、

• 個々の説明変数$x_j$において$MSE(\hat{b}^{OLS}_j)$と$MSE(\hat{b}^{PLS}_j)$を比べたとき、PLSが有利になるのはどういう状況か(私の理解では、$X$がランク落ちに近づいている状況でさえ、PLSが必ず有利になるとはいえない)
• $E(\hat{b}^{PLS}_j)$は真の$b_j = b^\circ_j = b^*_j$からどっち向きにどのくらい歪むのか(私の理解では、$\hat{b}^{PLS}$はその二乗和が$\hat{b}^{OLS}$の二乗和より小さいという意味では縮小推定量だけど、個別の$\hat{b}^{PLS}_j$は必ずしも0に向かって縮小するとは言えない)

という2点を教えてあげるほうが大事なのではなかろうか...
　というか、私はそういう話を聞きたかったのだと気持ちが整理できたので、読んでよかった、ということにしよう。逆にいうと、PLS回帰モデルをデータ生成メカニズムのモデルとしてまじに捉える人がやっぱしいるんだな(その割には$\hat{B}$なんかをみちゃうんだな)、という点が勉強になった。

論文：データ解析(2018-) - 読了：Burnham et al. (2001) PLS回帰モデルの「回帰係数」をどう解釈するか