主成分分析 のRによる実施例です。
主成分分析 の使い道はいろいろですが、以下は順を追って説明しています。
Rの実施例は下記になります。 下記は、コピーペーストでそのまま使えます。
サンプルのデータは下記を使っています。「Name」の列はなくても動きます。
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
DataName <- Data$Name # Nameの列を別名で保管する
Data$Name <- NULL # データからNameの列を消して、Xの列だけにする
pc <- prcomp(Data, scale=TRUE) # 主成分分析
summary(pc) # 「Cumulative Proportion」が累積寄与率。
上記のコードの続きです。 「他の多変量解析で、データの前処理に使う」をするために必要な主成分得点を得る方法です。
行番号でグラフを作ります。
pc1 <- pc$x # 主成分得点を得る
例えば、累積寄与率を見て２つ目の主成分までを使うことにするのなら、「pc1」の左の２列が次元削減後のデータということになります。
上記のコードの続きです。 「総合的にサンプルを見る」や「サンプルの仲間分け」の方法です。 これを進めるには、ggplot2のインストールが事前に必要になります。
行番号でグラフを作ります。
pc1 <- transform(pc1 ,name1 = DataName,name2 = "A")# サンプル名を追加
library(ggplot2)# ライブラリを追加
ggplot(pc1, aes(x=PC1, y=PC2,label=rownames(pc1))) + geom_text()# 第１主成分と第２主成分で言葉の散布図
「Name」の列があるデータを使う場合、「Name」でグラフを作ることもできます。
ggplot(pc1, aes(x=PC1, y=PC2,label=name1)) + geom_text()# 第１主成分と第２主成分で言葉の散布図
pc2 <- sweep(pc$rotation, MARGIN=2, pc$sdev, FUN="*") # 因子負荷量を計算
pc2 <- transform(pc2,name1=rownames(pc2),name2="B")# 言葉の散布図
ggplot(pc2, aes(x=PC1, y=PC2,label=name1)) + geom_text()# 言葉の散布図
３つの変数で似ているものがないことがわかりました。
「変数の仲間分け」１の続きです。
変数が多いと、上から２つの主成分で作った散布図上では重なって見える点が、別の主成分で見ると離れているということがおきます。
どの主成分と関係しているかではなく、変数の類似度を調べたい場合は、
多次元尺度構成法
を使って、
多数の主成分を２つに凝縮すると良いです。
MaxN = 5# 使用する固有値の数を指定
library(MASS)# ライブラリを読み込み
Data11 <- pc2[,1:MaxN]# 主成分の列を指定
Data11_dist <- dist(Data11)# サンプル間の距離を計算
sn <- sammon(Data11_dist) # 多次元尺度構成法
output <- sn$points# 得られた2次元データの抽出
Data2 <- cbind(output, pc2) # 元データと多次元尺度構成法の結果を合わせる。
ggplot(Data2, aes(x=Data2[,1], y=Data2[,2],label=name1)) + geom_text() # 言葉の散布図
左が、第１、第２の主成分で作った散布図で、右が、５番目までの主成分を２個の変数に凝縮して作った散布図です。
このデータは、「X1とX2」、「X3とX4」、「X5とX6」、「X7とX8」の組が高い相関になるように作ったものなので、変数の仲間分けとしては、
右側の散布図が欲しかったものになります。
元の変数と主成分の関係は、
２部グラフ
にするとわかるようになります。
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) # グラフを作る
因子負荷量は、プラスとマイナスがあり、絶対値が大きいほど、元の変数と相関が高いことを表します。
例えば、PC1という因子負荷量については、元の変数には、プラス側で相関が高い場合、マイナス側で相関が高い場合、相関が低い場合の３つに分かれます。
このことがわかるように、PC1という変数から、「PC1+」と「PC1-」という変数を作るようにして、どちらと相関が高いのかがわかるようにしてみました。
biplotを使うと、サンプルと変数を同時に描いた図が作れます。
biplot(pc) # グラフを作る
このグラフを作るためのデータは、下記で作れます。
Data1 <- rbind(pc1,pc2)# データを結合
広義の数量化Ⅲ類 と、このサイトで呼んでいる方法のひとつになります。
質的変数の前処理の方法が２種類あるので、下記には２種類を示しています。
この例では、Cドライブの「Rtest」というフォルダに、
「Data.csv」という名前でデータが入っている事を想定しています。
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
crs <- table(Data$X,Data$Y)# 分割表を作る
pc <- prcomp(crs, scale=TRUE) # 主成分分析
pc1 <- pc$x # 主成分得点を得る
pc1 <- transform(pc1 ,name = rownames(pc1))# 行名を追加
library(ggplot2) # パッケージの読み込み
ggplot(pc1, aes(x=PC1, y=PC2,label=name)) + geom_text()# 第１主成分と第２主成分で言葉の散布図
# ここまでは、X側のカテゴリの仲間分けでした。この後がY側のカテゴリの仲間分けです。
pc2 <- sweep(pc$rotation, MARGIN=2, pc$sdev, FUN="*") # 因子負荷量を計算
pc2 <- transform(pc2,name=rownames(pc2))# 名前を書き替え
ggplot(pc2, aes(x=PC1, y=PC2,label=name)) + geom_text()# 言葉の散布図
# ２つの結果を結合します。
pc3 <- rbind(pc1,pc2) # ２つの結果を結合
ggplot(pc3, aes(x=PC1, y=PC2,label=name)) + geom_text()# 言葉の散布図
分割表とグラフを見比べると、数字が大きいところの関係は、近くに配置されていることがわかります。
サンプルのデータは下記を使っています。「Name」の列はなくても動きます。
library(dummies) # ライブラリを読み込み
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
DataName <- Data$Name # Nameの列を別名で保管する
Data$Name <- NULL # データからNameの列を消して、Xの列だけにする
Data_dmy <- dummy.data.frame(Data)# ダミー変換
pc <- prcomp(Data_dmy, scale=TRUE) # 主成分分析
summary(pc) # 「Cumulative Proportion」が累積寄与率。
pc1 <- pc$x # 主成分得点を得る
pc1 <- transform(pc1 ,name = DataName)# サンプル名を追加
pc1$Index <-row.names(Data) # Indexという名前の列を作り、中身は行番号にする
library(ggplot2) # パッケージの読み込み
ggplot(pc1, aes(x=PC1, y=PC2,label=name)) + geom_text()# 第１主成分と第２主成分で言葉の散布図
Z4とZ7が重なっています。 数量化Ⅲ類は、主成分分析がベースなため、 サンプルの類似度の分析 をした時に、３次元以上にサンプルが割り振られていることがあり、散布図でうまく見れません。 ２次元で見たい場合は、 多次元尺度構成法 の方が良いです。
pc2 <- sweep(pc$rotation, MARGIN=2, pc$sdev, FUN="*") # 因子負荷量を計算
pc2 <- transform(pc2,nameCol=rownames(pc2))# 言葉の散布図
ggplot(pc2, aes(x=PC1, y=PC2,label=nameCol)) + geom_text()# 言葉の散布図
Tweet
