偏相関係数による有向相関分析 の実施例です。
下のようなデータを例にします。
5つの変数がありますが、X3に誤差を足すことで、その他の変数を作っています。
library(ppcor)
setwd("C:/Rtest")
Data <- read.csv("Data.csv", header=T)
cor2<-cor(Data)^2
pcor2 <- (pcor(Data)$estimate)^2
n <- ncol(cor2)
for (j in 1:n) {
for (i in 1:n) {
if(i != j){
if(cor2[i,j] < pcor2[i,j]){
if(i == 2 & j == 1){
directed_cor2 <- cor2
use_cor <- 1
}
if(cor2[i,j] < 0.01){
directed_cor2[,j] <- 0
}
} else {
if(i == 2 & j == 1){
directed_cor2 <- pcor2
use_cor <- 0
}
if(pcor2[i,j] < 0.01){
directed_cor2[i,] <- 0
}
}
}
}
}
round(cor2,3) # 相関係数の2乗
round(pcor2,3) # 偏相関係数の2乗
round(directed_cor2,3) # 有向相関分析のデータ
library(igraph)
GM1 <- directed_cor2
diag(GM1) <- 0
GM3 <- GM1*10 # 一番大きな値が10になるように修正
GM3[GM3<1] <- 0 # 1未満の場合は0にする(非表示にするため)
if(use_cor == 1){
GM4 <- graph.adjacency(GM3,weighted=T, mode = "directed") # グラフ用のデータを作成
} else {
GM4 <- graph.adjacency(GM3,weighted=T, mode = "undirected") # グラフ用のデータを作成
}
plot(GM4, edge.width=E(GM4)$weight) # グラフを作成
library(ppcor)
library(tidyr)
library(dplyr)
library(igraph)
setwd("C:/Rtest")
Data <- read.csv("Data.csv", header=T)
for (k in 1:100) { # ランダムサンプリングを100回実行
Data1 <- Data[,sample(ncol(Data), 3)]# ランダムに3列を選択する
cor2<-cor(Data1)^2
pcor2 <- (pcor(Data1)$estimate)^2
n <- ncol(cor2)
for (j in 1:n) {
for (i in 1:n) {
if(i != j){
if(cor2[i,j] < pcor2[i,j]){
if(i == 2 & j == 1){
directed_cor2 <- cor2
}
if(cor2[i,j] < 0.1){
directed_cor2[,j] <- 0
}
} else {
if(i == 2 & j == 1){
directed_cor2 <- pcor2
}
if(pcor2[i,j] < 0.1){
#directed_cor2[i,] <- 0
}
}
}
}
}
if(min(cor2, pcor2) > 0.01){
directed_cor2[,1:3] <- 0
}
diag(directed_cor2) <- 0
directed_cor22 <- data.frame(directed_cor2)
directed_cor22$Name <- row.names(directed_cor22)
directed_cor23 <- tidyr::gather(directed_cor22, key="X", value = Xs, -Name)
if(k == 1){
directed_cor24 <- directed_cor23
} else {
directed_cor24 <- rbind(directed_cor24,directed_cor23)
}
directed_cor24 <- directed_cor24[!(directed_cor24$Xs < 0.001),]
}
directed_cor25<-directed_cor24
directed_cor25[,3] <- NULL
directed_cor25 <- distinct(directed_cor25)
GM4 <- graph_from_edgelist(as.matrix(directed_cor25), directed = T)
plot(GM4)