Rによるデータ分析

Rによる行と列の項目の、項目同士の類似度の分析

データの表の値が、行と列の個々の項目の関係を表している場合の分析方法です。

関係を表す値が量的データで、大きいほど関係が強いことを表す場合

表の値は、頻度や関係の強さを表していて、大きいほど関係が強い場合です。 0が一番小さな値で、「関係なし」の意味になっているとします。 下記のデータの場合は、BとTが一番強い関係を表しているとします。
コードの例では、一番左の列はサンプル名として使われます。

2部グラフ

2部グラフ です。

library(igraph) # パッケージを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
DM <- read.csv("Data.csv", header=T, row.names=1)
# データを読み込み。拡張子は、列名、行名の読み込み
DM.mat = as.matrix(DM)/max(DM)*5
# データフレームを行列として読み込む。辺の太さを調整
DM.mat[DM.mat<1] <- 0
# 細い線は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, layout = layout_as_bipartite)
# グラフを作る

多次元同時付置図(コレスポンデンス分析+多次元尺度構成法+同時付置図)

データが0以上の整数の時に使える方法です。 例えば、頻度の時は使えます。

コレスポンデンス分析 の出力を 多次元尺度構成法 で処理して 多次元同時付置図 を作る方法です。
コレスポンデンス分析の結果が多次元になっている場合は、普通の同時付置図では、分析で間違いが起きる可能性があるので、 この方法を使うと便利です。

library(MASS) # パッケージを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T, row.names=1)
# データを読み込み
pc <- corresp(Data,nf=min(ncol(Data),nrow(Data)))
# コレスポンデンス分析
pc1 <- pc$cscore
# スコアを読み取り
pc1 <- transform(pc1 ,name1 = rownames(pc1), name2 = "A")
# 行名を追加
pc2 <- pc$rscore
# スコアを読み取り
pc2 <- transform(pc2 ,name1 = rownames(pc2), name2 = "B")
# 列名を追加
Data1 <- rbind(pc1,pc2)
# データを結合
round(pc$cor^2/sum(pc$cor^2),2)
# 寄与率を求める。

# 上記の結果から、3番目までの固有値は寄与率が高いことがわかったので、この後の解析にいれることにします。
# ちなみに、寄与率の低いところを入れると、それがノイズになるらしく、きれいに分離されなくなります。

MaxN = 3# 使用する固有値の数を指定
Data11 <- Data1[,1:MaxN]
# 項目のある列を指定
Data11_dist <- dist(Data11)
# サンプル間の距離を計算
sn <- sammon(Data11_dist)
# 多次元尺度構成法
output <- sn$points
# 得られた2次元データの抽出
Data2 <- cbind(output, Data1)  # 元データと多次元尺度構成法の結果を合わせる。
library(ggplot2)
# パッケージの読み込み#
ggplot(Data2, aes(x=Data2[,1], y=Data2[,2],label=name1)) + geom_text(aes(colour=name2))
# Nameを使った言葉の散布図

3部グラフ

3部グラフ を作る方法です。

この例では、Cドライブの「Rtest」というフォルダに、 「Data.csv」という名前で分割表になっているデータが入っている事を想定しています。
QM

library(MASS) # パッケージを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T, row.names=1)
# データを読み込み
pc <- corresp(Data,nf=min(ncol(Data),nrow(Data)))
# コレスポンデンス分析
pc1 <- pc$cscore
# スコアを読み取り
pc1 <- transform(pc1 ,name1 = rownames(pc1), name2 = "A")
# 行名を追加
pc2 <- pc$rscore
# スコアを読み取り
pc2 <- transform(pc2 ,name1 = rownames(pc2), name2 = "B")
# 列名を追加
Data1 <- rbind(pc1,pc2)
# データを結合
round(pc$cor^2/sum(pc$cor^2),2)
# 寄与率を求める。
QM 

「Data1」というデータは、下のようになっています。これをこの後のグラフ作成に使います。
QM

MaxN = 5# 使用する固有値の数を指定
library(igraph)
# パッケージを読み込み
library(sigmoid)
# パッケージを読み込み
Data1p = Data1[,1:MaxN]
# 分析対象のデータを抽出
names(Data1p) = paste(names(Data1p),"+",sep="")
# 列名を変更
DM.matp = apply(Data1p,c(1,2),relu)
# 0以下の値は0に変換
Data1m = -Data1[,1:MaxN]
# 分析対象のデータを抽出。符号を反転
names(Data1m) = paste(names(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) * 10
# 値が0から10になるように変換
DM.mat[DM.mat < 4] <- 0
# 絶対値が4未満の場合は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)$color[1:9] <- "red"
# 色を変える。「9」はAの要素の数
V(DM.g)$shape <- c("square", "circle")[V(DM.g)$type+1]
# マークの形を変える
plot(DM.g, edge.width=E(DM.g)$weight)
# グラフを作る
QM

対数線形分析

対数線形分析 を使って、数値と項目の関係を数式で表す方法です。

library(tidyr) # ライブラリの読み込み
library(MASS)
# ライブラリの読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- tidyr::gather(Data, key="Name2", value = Val, -Name)
# 縦型に変換(Nameの列以外を積み上げる)
summary(step(glm(Val~.^2, data=Data1,family=poisson)))
# 対数線形分析


この例だと、「(Name + Name2)^2」がモデル式の右辺として最終的に残っているので、 行と列の項目の単純な足し合わせ以外の要因で、値が決まっていることがわかります。
データによっては、例えば「Name2」が最後に残ることがあります。 この場合は、行の項目とは関係なく、値が決まっていることになります。

関係を表す値が量的データで、値の大きさが関係の強さを表さない場合

表の値は、量的データで、項目間の関係を表す数値ですが、値の大きさが関係の強さを表していないとします。

ヒートマップ

一番左の列を列名として入力すること以外は、 Rによるデータ全体の可視化 にある方法と基本的に同じです。

ヒートマップ を使う方法は、「関係を表す値が量的データで、大きいほど関係が強いことを表す場合」でも、もちろん使えます。

library(heatmaply) # ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T, row.names=1)
# データを読み込み
heatmaply(Data, Colv = NA, Rowv = NA)
# ヒートマップを作成
グラフの例

ヒートマップ×クラスター分析

ヒートマップを作るための関数のデフォルトは、なぜかクラスター分析をして、行と列を似たものになるように並べ替える機能が付いたものになっています。 これを使います。

library(heatmaply) # ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T, row.names=1)
# データを読み込み
heatmaply(Data)
# クラスター分析をしてからヒートマップを作成
グラフの例

行列の分解1:特異値分解(svd)

Rで 特異値分解(svd) をする時の例は、下記になります。
NMF

「C:/Rtest」のフォルダに、「Data.csv」という名前のcsvファイルがある事を想定しています。 また、このページの例のように、一番上の行に列の名前、一番左の列に行の名前が書かれているデータを想定しています。

library(ggplot2) # パッケージの読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data0 <- read.csv("Data.csv", header=T, row.names=1)
# データを読み込み
Data1 = as.matrix(Data0)
# データフレームを行列として読み込み
Pt <- svd(Data1)
# パターンが2個になるようにNMF
PtB = Pt$u
# 行側のパターンを作成
PtA = Pt$v
# 列側のパターンを作成
# 列側のパターンをグラフ化(変数の仲間分けの分析)
A = data.frame(PtA)
# 行列を転置してからフレームにする
A$Name <-colnames(Data1)
# 名前の列を作る
ggplot(A, aes(x=A[,1], y=A[,2],label=Name)) + geom_text()
# 言葉の散布図
NMF
# 行側のパターンをグラフ化(サンプルの仲間分けの分析)
B = data.frame(PtB)
# 行列をフレームにする
B$Name <-rownames(Data1)
# 名前の列を作る
ggplot(B, aes(x=B[,1], y=B[,2],label=Name)) + geom_text()
# 言葉の散布図
NMF

行列の分解2:非負行列分解(NMF)

Rで 非負行列分解(NMF) をする時の例は、下記になります。
NMF

library(NMF) # パッケージの読み込み
library(ggplot2)
# パッケージの読み込み

setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data0 <- read.csv("Data.csv", header=T, row.names=1)
# データを読み込み
Data1 = as.matrix(Data0)
# データフレームを行列として読み込み
Pt <- nmf(Data1, 2)
# パターンが2個になるようにNMF
PtB = basis(Pt)
# 行側のパターンを作成
PtA = coef(Pt)
# 列側のパターンを作成
# 列側のパターンをグラフ化(変数の仲間分けの分析)
A = data.frame(t(PtA))
# 行列を転置してからフレームにする
A$Name <-row.names(A)
# 名前の列を作る
ggplot(A, aes(x=A[,1], y=A[,2],label=Name)) + geom_text()
# 言葉の散布図
NMF
# 行側のパターンをグラフ化(サンプルの仲間分けの分析)
B = data.frame(PtB)
# 行列をフレームにする
B$Name <-row.names(B)
# 名前の列を作る
ggplot(B, aes(x=B[,1], y=B[,2],label=Name)) + geom_text()
# 言葉の散布図
NMF

なお、筆者が最初に上記のコードを使った時は、何も問題がなかったのですが、しばらくしてから使ったら、NMFのインストールでエラーが出るようになりました。 エラーメッセージを見て、「Biobase」と「BiocGenerics」というパッケージをインストールしたところ、エラーはなくなりました。 これらの2つのパッケージはCRANにはないので、インストールはBioconductorというサイトから圧縮ファイルをダウンロードして、 それを展開して使うことになりました。

層別の1次元散布図

層別1次元散布図 を使って、体系的なグラフにする方法です。

library(tidyr) # ライブラリの読み込み
library(ggplot2)
# ライブラリの読み込み#
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- tidyr::gather(Data, key="Name2", value = Val, -Name)
# 縦型に変換(Nameの列以外を積み上げる)
ggplot(Data1, aes(x=Name, y=Val)) + geom_point() +facet_grid(.~Name2)
# 層別1次元散布図

NameとName2を入れ替えると、グラフの分け方が入れ替わります。
ggplot(Data1, aes(x=Name2, y=Val)) + geom_point() +facet_grid(.~Name) # 層別1次元散布図


数量化T類

数量化T類 を使って、数値と項目の関係を数式で表す方法です。

library(tidyr) # ライブラリの読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- tidyr::gather(Data, key="Name2", value = Val, -Name)
# 縦型に変換(Nameの列以外を積み上げる)
summary(step(lm(Val~., data=Data1)))
# 数量化T類


決定係数(R-squared)が0.9991なので、非常に高いモデル式ができていることがわかります。
結果の見方ですが、「TとB」の組み合わせの推論の値は、
5.010 + 6.120 -3.100 = 8.03
ということになります。 多重共線性 の回避のために、AとSがダミー変数から抜けていますが、例えば、 「AとT」の組み合わせの値は、
5.010 + 6.120 = 11.13
と求めます。

なお、結果として、有意なモデル式が作られた場合は、複雑な規則で値が決まっていることは推察できますが、 どのような規則なのかはわかりません。

関係を表す値が質的データの場合

表の値は、質的データの場合です。

いずれの方法も、行の項目名が1列目、列の項目が2列目、値が3列目になるようにデータを並べ変えてから使います。

3部グラフ

library(tidyr) # ライブラリの読み込み
library(dplyr)
# ライブラリの読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- tidyr::gather(Data, key="Name2", value = Val, -Name)
# 縦型に変換(Nameの列以外を積み上げる)
Data21 <- Data1
# 縦型に変換(Nameの列以外を積み上げる)
Data21[,2] <- NULL #
# データから2列目を消す
Data22 <- Data1[,2:3] #
# 1列目のないデータを作る
names(Data22)[1] <- "Name"
# 列名を変更
Data31 <- count(group_by(Data21,Data21[,1:2],.drop=FALSE)) #
# 重複している組み合わせを集計
Data32 <- count(group_by(Data22,Data22[,1:2],.drop=FALSE)) #
# 重複している組み合わせを集計
Data4 <- rbind(Data31,Data32)
# データを縦に結合
Data5 <- Data4/max(Data4)*5
# 辺の太さを調整
DM.g<- graph.data.frame(Data5[,c(1,2)], directed = F)
# グラフ用のデータを作る
E(DM.g)$weight <-Data5[[3]]
# 辺のデータを取る
plot(DM.g, edge.width=E(DM.g)$weight)
# グラフを作る


この結果は、解釈しにくいですが、とりあえずデータがグラフになりました。

コレスポンデンス分析(多重対応分析)+多次元尺度構成法

library(tidyr) # ライブラリの読み込み
library(dplyr)
# ライブラリの読み込み
library(dummies)
# ライブラリの読み込み
library(MASS)
# ライブラリの読み込み
library(ggplot2)
# ライブラリの読み込み#
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- tidyr::gather(Data, key="Name2", value = Val, -Name)
# 縦型に変換(Nameの列以外を積み上げる)
Data_dmy <- dummy.data.frame(Data1)
# ダミー変換
pc <- corresp(Data_dmy,nf=min(ncol(Data_dmy)))
# コレスポンデンス分析
pc1 <- pc$cscore
#
pc1 <- transform(pc1 ,name1 = rownames(pc1))
# 行名を追加
round(pc$cor^2/sum(pc$cor^2),2)
# 寄与率を求める。

上の例では、9番目以降の固有値は寄与率が低いので、この後の解析から外すことにします。

MaxN = 8# 使用する固有値の数を指定
Data11 <- pc1[,1:MaxN]
# 項目のある列を指定
Data11_dist <- dist(Data11)
# サンプル間の距離を計算
sn <- sammon(Data11_dist)
# 多次元尺度構成法
output <- sn$points
# 得られた2次元データの抽出
Data2 <- cbind(output, pc1)  # 元データと多次元尺度構成法の結果を合わせる。
ggplot(Data2, aes(x=Data2[,1], y=Data2[,2],label=name1)) + geom_text()
# Nameを使った言葉の散布図

これも、解釈しにくいですが、とりあえずデータがグラフになりました。 Redが外れたところに配置されるのは、上の方法と同じになりました。



Tweet データサイエンス教室