library(psych) # ライブラリを読み込み
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
Data1 <- dummy.data.frame(Data)# 質的変数はダミー変換
pairs.panels(Data1) # 総当たりの散布図を作る
変数の類似度の分析 です。
もともとの方法は、量的変数だけのグループや、質的変数だけのグループのための手法ですが、 いずれに対しても、このページのコードは量的、質的の両方が使えるようにしてあります。
基本的に量的変数を扱う方法ですが、質的変数は ダミー変換 して質的・量的が混合していたり、質的変数だけでも使えるようにしてあります。
この方法の場合、質的変数同士の分析結果は、 Rによる個々のカテゴリの類似度の分析 と似て来ます。
量的変数の類似度を見る方法として散布図は基本なので、 総当たりの散布図を見ることは大事です。
ただし、下記は変数が９個の場合ですが、たった９個でもこの後の分析が進めにくいグラフができあがります。 他の方法を使って、グラフにする変数は絞った方が良いです。
library(dummies) # ライブラリを読み込み
library(psych) # ライブラリを読み込み
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
Data1 <- dummy.data.frame(Data)# 質的変数はダミー変換
pairs.panels(Data1) # 総当たりの散布図を作る
因果関係の分析をする時などは、総当たりの散布図は必要なく、１つの変数と他のすべての変数の組み合わせについて、散布図を見たいことがあります。 その時の作り方です。
下記のコードでは、１つの変数というのは、「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")# 散布図を大量に描く
注目したい変数が質的変数の場合です。
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")# １次元散布図
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) # グラフを作成
グラフィカルラッソ をしてから、ネットワークグラフにする方法です。
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) # グラフを作成
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) # グラフを作成
主成分分析 をしてから、多次元の主成分を、 多次元尺度構成法 で可視化する方法です。
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() # 言葉の散布図
基本的に質的変数を扱う方法ですが、量的変数は １次元クラスタリング の方法で、質的変数に変換するコードが入っているので、 質的・量的が混合していたり、量的変数だけでも使えるようにしてあります。
上記の量的変数の類似度の分析をベースにした方法では、線形の関係しか見ていませんが、 質的変数の類似度の分析をベースにした方法を使うと、量的変数の非線形の関係も見ることができるようになります。 数値の情報が粗くなる弱点はありますが、探索的データ分析では、あまり気にしなくて良いと思います。
連関係数 を計算して、相関の大きいものだけをネットワークグラフで抽出する方法です。
library(igraph) #ライブラリを読み込み
library(vcd) #ライブラリを読み込み
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T) # データを読み込み
Data1 <- Data # １次元クラスタリングの出力先を作る
n <- ncol(Data1) # データの列数を数える
for (i in 1:n) { # ループの始まり。データの列数を数えて同じ回数繰り返す
if (class(Data1[,i]) == "numeric") { # 条件分岐の始まり
Data1[,i] <- droplevels(cut(Data1[,i], breaks = 5,include.lowest = TRUE))# ５分割する場合。量的データは、質的データに変換する。
} # 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) # グラフを作成
対数線形分析 です。
library(dplyr) #ライブラリを読み込み
library(MASS) # 作ライブラリの読み込み
setwd("C:/Rtest") # 作業用ディレクトリを変更
Data <- read.csv("Data.csv", header=T, stringsAsFactors=TRUE) # # データを読み込み
Data1 <- Data # １次元クラスタリングの出力先を作る
nc <- ncol(Data1) # データの列数を数える
for (i in 1:nc) { # ループの始まり。データの列数を数えて同じ回数繰り返す
if (class(Data1[,i]) == "numeric") { # 条件分岐の始まり
Data1[,i] <- droplevels(cut(Data1[,i], breaks = 5,include.lowest = TRUE))# ５分割する場合。量的データは、質的データに変換する。
} # if文の処理の終わり
} # ループの終わり
Data2 <- count(group_by(Data1,Data1[,1:nc],.drop=FALSE)) # 分割表を作る
gm <- step(glm(n~.^2, data=Data2,family=poisson)) # 対数線形分析
summary(gm) # 結果の出力
この結果では、A、B、Cという変数のグループと、D、E、Fという変数のグループがあることがわかりました。
上のコードの場合は、２変数の交互作用項がすべて入るモデルをスタートにして、変数の絞り込みが始まります。
この場合は、２つの組み合わせを見ている連関係数を使う方法と、あまり変わらないことをしています。
例えば、３変数の交互作用項も入れたい場合は、下記にします。
この場合は、計算時間が相当かかるだけでなく、n数不足が原因と思われるエラーが発生しやすいです。
変数が数個なら、問題なくできます。
筆者が試した時は、変数が６個の時はダメで、３個だとできました。
gm <- step(glm(n~.^3, data=Data2,family=poisson)) # 対数線形分析
summary(gm) # 結果の出力
祐子 武田 氏の記事
質的変数の総当たりの組み合わせから、連関係数を求める方法は、このページを参考にさせていただきました。
https://qiita.com/ytakeda/items/058e83ebdd721f87ceb4
Tweet
