Rによるデータ分析

Rによる主成分回帰分析

主成分回帰分析 の使い方の紹介になります。

主成分回帰分析の基本

CドライブのRtestというフォルダに、Data.csvというファイルがある状況を想定しています。 目的変数は「Y」という名前になっている必要があります。 Y以外のすべての変数を主成分分析するようになっていますので、説明変数の名前はなんでも良いです。

抽出する主成分は第３位までにしてあります。 主成分を何位まで使うのかは、累積寄与率を考えて決める必要があります。

setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
DataY <- Data$Y # Yの列を別名で保管する
Data$Y <- NULL # データからYの列を消して、Xの列だけにする
pc <- prcomp(Data, scale=TRUE) # Xを主成分分析
summary(pc) # 「Cumulative Proportion」が累積寄与率。これを使って、何個目の主成分までを使うのかを決められる。

この例の場合は、累積寄与率がPC3までで99%になっているので、PC3までの主成分で回帰分析をすることにします。
pc1 <- predict(pc, Data)[,1:3] # 主成分を第３位まで作る
DataPCR <- transform(pc1, Y = DataY)　 # Yと主成分で新しいデータセットを作る
pcr <- lm(Y ~ . ,DataPCR) # 主成分回帰分析
summary(pcr) #　結果を出力

この例の場合は、t-Valueを見ると、PC1だけがとても大きいので、PC1だけでほぼ決まっているモデルと考えられます。

主成分分析による説明変数の分析（パターン１）

上記でYと、主成分のPC1に相関があることがわかりましたが、これではYと説明変数の関係がわかりません。 そこで、PC1と説明変数の関係を主成分分析で調べることにします。 この方法は、 Rによる主成分分析 のページにあります。

MaxN = 3# 使用する固有値の数を指定
library(MASS)# ライブラリを読み込み
pc2 <- sweep(pc$rotation, MARGIN=2, pc$sdev, FUN="*") # 因子負荷量を計算
pc2 # 因子負荷量を表示

この結果を見ると、 PC1は、X01、X02、X03と相関が高いことがわかります。

主成分と元の変数の相関の強さの全体像を見るためのグラフを作ってみます。
Data11 <- pc2[,1:MaxN]# 主成分の列を指定
Data11_dist <- dist(Data11)# サンプル間の距離を計算
library(igraph) # パッケージを読み込み
library(sigmoid) # パッケージを読み込み
Data1p = Data11# 分析対象のデータを抽出
colnames(Data1p) = paste(colnames(Data1p),"+",sep="")# 列名を変更
DM.matp = apply(Data1p,c(1,2),relu)# 0以下の値は0に変換
Data1m = -Data11# 分析対象のデータを抽出。符号を反転
colnames(Data1m) = paste(colnames(Data1m),"-",sep="")# 列名を変更
DM.matm = apply(Data1m,c(1,2),relu)# 0以下の値は0に変換
DM.mat =cbind(DM.matp,DM.matm)# プラス側とマイナス側を合体
DM.mat <- DM.mat / max(DM.mat) * 3 # 値が0から3になるように変換（線の太さになる）
DM.mat[DM.mat < 1] <- 0 #　1未満の場合は0にする（非表示にする）
DM.g<-graph_from_incidence_matrix(DM.mat,weighted=T) # グラフ用のデータを作る
V(DM.g)$color <- c("steel blue", "orange")[V(DM.g)$type+1] # 色を変える
V(DM.g)$shape <- c("square", "circle")[V(DM.g)$type+1] # マークの形を変える
plot(DM.g, edge.width=E(DM.g)$weight) # グラフを作る

６つあった変数は、３つの主成分に分類できていることがわかります。

主成分分析による説明変数の分析（パターン２）

上記のパターン１は、主成分回帰分析による分析で、簡単な場合です。 たいていの場合は、パターン１が当てはまります。

下記のパターン２は、主成分回帰分析で、考察が難しくなる場合です。 難しいのですが、主成分を使わない普通の回帰分析では、もっと考察が難しくなるデータのパターンです。

上記のコードで、回帰分析が以下のようになったとします。 Ｙと相関の高い変数はPC2であることがわかります。

因子負荷量を見ると、どの変数もPC1と相関が高く、PC2と相関の高い変数はないです。 こうなって来ると、パターン１のような考察ができないです。

Rによる隠れ変数の分析 にある方法を使って、散布図を作ると、全体的な相関関係がよくわかります。 Yと相関が高いのは、PC2だけで、その他の主成分も元データの変数もYとは相関が低いことがわかります。

よくよくPC2の因子負荷量を見ると、X01、X02は符号がプラスで、X03、X04は符号がマイナスなことがわかります。

主成分回帰分析からは、変数はX01、X02と、X03、X04の２つのグループに分かれ、 このグループの違いにある何かが、Ｙと非常に高い相関をしていることがわかります。 ちなみに、このサイトでは、この「何か」を、 隠れ変数 と呼んでいます。

主成分分析では、隠れ変数がどのように説明変数に入っているのかがわからないので、 因子分析 をしてみます。 主成分回帰分析では、Yを元データから分けましたが、ここでは一緒にするのがポイントです。
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
factanal(Data,2)$loadings # 潜在変数を２個で仮定して、因子分析

すると、Ｙと共通の潜在変数が含まれているのは、X03とX04の２つであることがわかります。

ちなみに、psychというライブラリでも同じ結論を出せます。
library(psych) # ライブラリの読み込み
library(GPArotation) # ライブラリの読み込み
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
fa(Data, nfactors = 2, fm = "ml", rotate = "varimax")$loadings　# 潜在変数を２と仮定して、因子負荷量を出力

こうすると、X01、X02のグループとは違う何かが、X03、X04のグループには少し含まれていて、それがＹと関係していることまでわかるようになります。

ちなみに、Ｙを一緒にするやり方を、主成分分析や、独立成分分析でしても、少なくともこのデータについては、すっきりとした結果になりませんでした。

なお、目的変数と説明変数の関係にある変数を、同列に考えて、因子分析をしてしまうやり方は、 因果関係の構造がおかしいので、因子分析の本来の使い方すれば誤用になります。 本来は、 SEM・共分散構造分析 をするものになります。

ただ、 SEM・共分散構造分析 は手間がかかりますので、 因果関係の分析をするための簡単な方法として、ここでは目的変数と説明変数を同列にして扱う方法を使っています。