R で主成分回帰分析を行う方法
概要
主成分回帰分析pCRとは、主成分分析と重回帰分析を合わせたもので、主成分分析を通じて得られた主成分を新たな独立変数とした回帰分析をいう。実際、統計学の観点から主成分分析それ自体はあまり必要なく、通常は回帰分析に用いられるときにこそ意味を持つ。
実習
(多重共線性を見つける方法の続き)
主成分を作り出す過程は行列分解を含む複雑な計算を伴うが、Rではprcomp()関数を通じて簡単に求めることができる。
- (1) Rotation: 実際に元の独立変数に掛けて主成分を作る係数を表す。例えばPC1を作るためには、populationに$-0.098$、whiteに$0.056$、… 、povertyに$-0.187$を掛けてすべて足すことで作られる。またデータによっては、どの主成分がどの変数についての情報を多く持っているかを一目で確認できる指標となることもある。
- (2) Importance of components: 文字通り各主成分の重要度を表す。Proportion of Varianceが高いほどデータの多くの部分を説明し、それだけ重要な主成分であるという意味になる。これを視覚的に表したものが右側のScree Plotである。Scree plotはただplot()関数に主成分分析の結果を入れるだけでも簡単に描くことができる。この図を通じて主成分をどこまで使うのかを決定する。
- (3) Elbow Point: Scree plotで急激に傾きが変わる地点であり、これ以降は主成分をさらに使っても得られる利益が少なくなる部分をいう。Elbow point以降の主成分は分析に使わなくても特に問題ない。例題のような場合には、PC1、PC2、PC3までのみを使うのが望ましいだろう。
ところが結局これをどのように回帰分析に組み込むかが問題である。独立変数をRotationと掛けてPCを作り、回帰分析をした後に元の独立変数についての式を立てるのは非常に面倒なことなので、pls パッケージのpcr()関数を使ってみよう。
summary()関数を通じて出力された結果では、いくつの主成分を使ったときにcollgeGradをどれだけ説明するかを教えてくれる。
上の主成分分析では主成分の個数は$3$個で十分だった。jack.test()関数のncompオプションに3を入れて実行してみると、主成分についての回帰係数ではなく元の独立変数についての回帰係数を教えてくれるだけでなく、t検定まで一度に行ってくれる。
また、右側の残差プロットを見ると、特別な問題はないように見える。
win.graph(4,4); plot(scale(PCRout$residuals[,,3]),main="Residual")
残差プロットを出力する際に注意すべきことは、pcr()関数が返す結果の$residualsを参照することで得られるのは正しいが、こうするだけでは主成分をいくつ使ったかに応じた残差がすべて含まれた3次元配列として出てくるということである。この例題では[,,3]でスライシングすることで、主成分が3個まで使われたPCRの残差を参照した。もし主成分をk個だけ使用した場合は[,,k]のようにスライシングしなければならない。
コード
以下は例題コードである。
install.packages('pls')
install.packages('car')
library(pls)
library(car)
X<-MplsDemo[,-1]
Y<-X[,7]; X<-X[,-7]
PCAout<-prcomp(X[,-7],center=T,scale.=T); PCAout; summary(PCAout)
win.graph(3,3); plot(PCAout, type='l',main="screeplot")
PCRout<-pcr(collegeGrad~.,data=MplsDemo[,-1],validation="LOO",jackknife=T)
summary(PCRout)
jack.test(PCRout,ncomp=3)
win.graph(4,4); plot(scale(PCRout$residuals[,,3]),main="Residual")