Rによるデータ分析 | Pythonによる変数の類似度の分析

Rによる変数の類似度の分析

変数の類似度の分析 です。

もともとの方法は、量的変数だけのグループや、質的変数だけのグループのための手法ですが、 いずれに対しても、このページのコードは量的、質的の両方が使えるようにしてあります。

量的変数の類似度の分析方法をベースにした方法

基本的に量的変数を扱う方法ですが、質的変数は ダミー変換 して質的・量的が混合していたり、質的変数だけでも使えるようにしてあります。

この方法の場合、質的変数同士の分析結果は、 Rによる個々のカテゴリの類似度の分析 と似て来ます。

総当たりの散布図

量的変数の類似度を見る方法として散布図は基本なので、 総当たりの散布図を見ることは大事です。

ただし、下記は変数が9個の場合ですが、たった9個でもこの後の分析が進めにくいグラフができあがります。 他の方法を使って、グラフにする変数は絞った方が良いです。

library(dummies) # ライブラリを読み込み
library(psych)
# ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- dummy.data.frame(Data)
# 質的変数はダミー変換
pairs.panels(Data1)
# 総当たりの散布図を作る

層別の折れ線グラフ

1つの量的変数と、他のすべての変数の組み合わせのグラフ

因果関係の分析をする時などは、総当たりの散布図は必要なく、1つの変数と他のすべての変数の組み合わせについて、散布図を見たいことがあります。 その時の作り方です。

下記のコードでは、1つの変数というのは、「Y」という名前にしておく必要があります。 他の変数の名前は、特に指定はありません。

まずは、注目したい変数が量的変数の場合です。

library(dummies) # ライブラリを読み込み
library(ggplot2)
# ライブラリを読み込み
library(tidyr)
# ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data2 <- dummy.data.frame(Data)
# 質的変数はダミー変換
Data_long <- tidyr::gather(Data2, key="X", value = Xs, -Y)
# 縦型に変換(Yの列以外を積み上げる)
ggplot(Data_long, aes(x=Xs,y=Y)) + geom_point() + facet_wrap(~X,scales="free")
# 散布図を大量に描く
層別の散布図

1つの質的変数と、他のすべての変数の組み合わせのグラフ

注目したい変数が質的変数の場合です。

library(dummies) # ライブラリを読み込み
library(ggplot2)
# ライブラリを読み込み
library(tidyr)
# ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- Data
# 加工用のデータを作る
Y <- Data1$Y
# Yの列を別に作っておく
Data1$Y <- NULL
# Yの列を消す
Data2 <- dummy.data.frame(Data1)
# 質的変数はダミー変換
Data3 <- cbind(Data2, Y)
# 加工したデータとYの列を合わせる
Data_long <- tidyr::gather(Data3, key="X", value = Xs, -Y)
# 縦型に変換(Yの列以外を積み上げる)
ggplot(Data_long, aes(x=Y,y=Xs)) + geom_point() + facet_wrap(~X,scales="free")
# 1次元散布図
層別のグラフ

ggplot(Data_long, aes(x=Y,y=Xs)) + geom_jitter(size=1, position=position_jitter(0.1)) + facet_wrap(~X,scales="free")# ジター散布図
層別のグラフ

ggplot(Data_long, aes(x=Y,y=Xs)) + geom_boxplot() + facet_wrap(~X,scales="free")# 箱ひげ図
層別のグラフ

相関係数+ネットワークグラフ

相関係数を計算して、相関の大きいものだけをネットワークグラフにする方法です。

library(dummies) # ライブラリを読み込み
library(igraph)
#ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- dummy.data.frame(Data)
# 質的変数はダミー変換
DataM <- as.matrix(Data1)
# テーブルデータを行列にする
GM1 <- cor(DataM)
# 相関係数行列を計算
diag(GM1) <- 0
# 対角成分は0にする
GM2 <- abs(GM1)
# 絶対値にする
GM2[GM2<0.9] <- 0
# 相関係数の絶対値が0.9未満の場合は0にする(非表示にするため)
GM3 <- GM2*10
# 一番大きな値が10になるように修正(パスの太さを指定するため)
GM4 <- graph.adjacency(GM3,weighted=T, mode = "undirected")
# グラフ用のデータを作成
plot(GM4, edge.width=E(GM4)$weight)
# グラフを作成

cor

グラフィカルラッソ+ネットワークグラフ

グラフィカルラッソ をしてから、ネットワークグラフにする方法です。

library(dummies) # ライブラリを読み込み
library(glasso)
# ライブラリを読み込み
library(igraph)
# ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- dummy.data.frame(Data)
# 質的変数はダミー変換
DataM <- as.matrix(Data1)
# テーブルデータを行列にする
COR <- cor(DataM)
# 相関係数行列を計算
RHO <- 0.2
# スパース性を調整するパラメタ
GM1 <- glasso(COR,RHO)$wi
# GMを実施(数値が大きいほどスパース。数値は結果を見ながら要調整)
# ここからは、グラフの描画用のデータの作成
diag(GM1) <- 0
# 対角成分は0にする
GM2 <- abs(GM1)
# 絶対値にする
GM2max <- max(GM2)
# 要素の最大値を求める
GM3 <- GM2/GM2max*10
# 一番大きな値が10になるように修正(パスの太さを指定するため)
GM3[GM3<5] <- 0
# 5未満の場合は0にする(非表示にするため)
rownames(GM3) <- rownames(COR)
# 行の名前を追記
colnames(GM3) <- colnames(COR)
# 列の名前を追記
GM4 <- graph.adjacency(GM3,weighted=T, mode = "undirected")
# グラフ用のデータを作成
plot(GM4, edge.width=E(GM4)$weight)
# グラフを作成

glasso

LiNGAM+ネットワークグラフ

LiNGAM をしてから、ネットワークグラフにする方法です。

ここでは、変数の影響の大きさを知りたいので、各変数は、 標準化 して、単位に依存して係数の大きさが決まらないようにしています。 そして、モデルの係数の絶対値を有向グラフにします。 そのため、スタートにしたデータを パス解析 するためのモデル式ではなく、標準化したデータに対してのパス解析のモデル式が求まります。

library(pcalg) #ライブラリの読み込み
library(igraph) #ライブラリの読み込み
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
for (i in 1:ncol(Data)) {
Data[,i] <- (Data[,i] - mean(Data[,i]))/sd(Data[,i]) #データを標準化
}
Str <- lingam(Data)$Bpruned # LiNGAM
rownames(Str) <- colnames(Data) # 行名を追記
colnames(Str) <- colnames(Data) # 列名を追記
GM2 <- abs(Str) # 絶対値にする
GM3 <- (GM2*5) # 値を適度な大きさにする
GM4 <- graph.adjacency(GM3,weighted=T, mode = "directed") # グラフ用のデータを作成
plot(GM4, edge.width=E(GM4)$weight) # グラフを作成

LiNGAM

主成分分析+多次元尺度構成法+散布図

主成分分析 をしてから、多次元の主成分を、 多次元尺度構成法 で可視化する方法です。

library(dummies) # ライブラリを読み込み
library(ggplot2)
# ライブラリを読み込み
library(MASS)
# ライブラリを読み込み
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
Data1 <- dummy.data.frame(Data)
# 質的変数はダミー変換
pc <- prcomp(Data1, scale=TRUE)
# 主成分分析
summary(pc)
# 「Cumulative Proportion」が累積寄与率。

pc2 <- sweep(pc$rotation, MARGIN=2, pc$sdev, FUN="*") # 因子負荷量を計算
MaxN = 5
# 使用する固有値の数を指定
Data11 <- pc2[,1:MaxN]
# 主成分の列を指定
Data11_dist <- dist(Data11)
# サンプル間の距離を計算
sn <- sammon(Data11_dist)
# 多次元尺度構成法
Data2 <- sn$points
# 得られた2次元データの抽出
Data2 <- transform(Data2 ,name1 = rownames(Data2))
# 行名を追加
ggplot(Data2, aes(x=Data2[,1], y=Data2[,2],label=name1)) + geom_text()
# 言葉の散布図
PCA

質的変数の類似度の分析方法をベースにした方法

基本的に質的変数を扱う方法ですが、量的変数は 1次元クラスタリング の方法で、質的変数に変換するコードが入っているので、 質的・量的が混合していたり、量的変数だけでも使えるようにしてあります。

上記の量的変数の類似度の分析をベースにした方法では、線形の関係しか見ていませんが、 質的変数の類似度の分析をベースにした方法を使うと、量的変数の非線形の関係も見ることができるようになります。 数値の情報が粗くなる弱点はありますが、探索的データ分析では、あまり気にしなくて良いと思います。

連関係数+ネットワークグラフ

連関係数 を計算して、相関の大きいものだけをネットワークグラフで抽出する方法です。

library(igraph) #ライブラリを読み込み
library(vcd)
#ライブラリを読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T)
# データを読み込み
Data1 <- Data
# 1次元クラスタリングの出力先を作る
n <- ncol(Data1)
# データの列数を数える
for (i in 1:n) {
# ループの始まり。データの列数を数えて同じ回数繰り返す
if (class(Data1[,i]) == "numeric") {
# 条件分岐の始まり
Data1[,i] <- droplevels(cut(Data1[,i], breaks = 5,include.lowest = TRUE))
# 5分割する場合。量的データは、質的データに変換する。
}
# if文の処理の終わり
}
# ループの終わり
GM2 <- matrix(0,nrow=n,ncol=n)
# 出力先を作る
for (i in 1:(n-1)) {
# ループの始まり
for (j in (i+1):n) {
# ループの始まり
cross<-xtabs(~Data1[,i]+Data1[,j],data=Data1)
# 分割表の作成
res<-assocstats(cross)
# 連関分析
cramer_v<-res$cramer
# クラメールの連関係数の抽出
GM2[i,j] <- cramer_v
# クラメールの連関係数の書き込み
GM2[j,i] <- cramer_v
# クラメールの連関係数の書き込み
}
# ループの終わり
}
# ループの終わり
rownames(GM2)<-colnames(Data1)
# 行名をつける
colnames(GM2)<-colnames(Data1)
# 列名をつける
GM2[GM2<0.8] <- 0
# 相関係数の絶対値が0.8未満の場合は0にする(非表示にするため)
GM3 <- GM2*10
# 一番大きな値が10になるように修正(パスの太さを指定するため)
GM4 <- graph.adjacency(GM3,weighted=T, mode = "undirected")
# グラフ用のデータを作成
plot(GM4, edge.width=E(GM4)$weight)
# グラフを作成

cramer

対数線形分析

対数線形分析 です。

library(dplyr) #ライブラリを読み込み
library(MASS)
# 作ライブラリの読み込み
setwd("C:/Rtest")
# 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T, stringsAsFactors=TRUE) #
# データを読み込み
Data1 <- Data
# 1次元クラスタリングの出力先を作る
nc <- ncol(Data1)
# データの列数を数える
for (i in 1:nc) {
# ループの始まり。データの列数を数えて同じ回数繰り返す
if (class(Data1[,i]) == "numeric") {
# 条件分岐の始まり
Data1[,i] <- droplevels(cut(Data1[,i], breaks = 5,include.lowest = TRUE))
# 5分割する場合。量的データは、質的データに変換する。
}
# if文の処理の終わり
}
# ループの終わり
Data2 <- count(group_by(Data1,Data1[,1:nc],.drop=FALSE))
# 分割表を作る
gm <- step(glm(n~.^2, data=Data2,family=poisson))
# 対数線形分析
summary(gm)
# 結果の出力
cramer
この結果では、A、B、Cという変数のグループと、D、E、Fという変数のグループがあることがわかりました。

上のコードの場合は、2変数の交互作用項がすべて入るモデルをスタートにして、変数の絞り込みが始まります。 この場合は、2つの組み合わせを見ている連関係数を使う方法と、あまり変わらないことをしています。
例えば、3変数の交互作用項も入れたい場合は、下記にします。 この場合は、計算時間が相当かかるだけでなく、n数不足が原因と思われるエラーが発生しやすいです。 変数が数個なら、問題なくできます。 筆者が試した時は、変数が6個の時はダメで、3個だとできました。
gm <- step(glm(n~.^3, data=Data2,family=poisson))
# 対数線形分析
summary(gm)
# 結果の出力



参考文献

祐子 武田 氏の記事
質的変数の総当たりの組み合わせから、連関係数を求める方法は、このページを参考にさせていただきました。
https://qiita.com/ytakeda/items/058e83ebdd721f87ceb4


Tweet データサイエンス教室