相関偏相関分析 、 有向相関偏相関分析 、 ハイブリッド有向相関分析 の実施例です。
有向相関偏相関分析 の事前分析が、 相関偏相関分析 になっています。 さらに、 ハイブリッド有向相関分析 の事前分析が、 有向相関偏相関分析 になっています。
下のようなデータを例にします。
library(ppcor)
library(tidyr)
library(dplyr)
library(igraph)
setwd("C:/Rtest")
Data <- read.csv("Data.csv", header=T)
cor2<-cor(Data)^2
pcor2 <- (pcor(Data)$estimate)^2
directed_cor2 <- cor2
n <- ncol(cor2)
for (j in 1:n) {
for (i in 1:n) {
if(i < j){
if(cor2[i,j] < pcor2[i,j]){
directed_cor2[i,j] <- cor2[i,j]
directed_cor2[j,i] <- cor2[i,j]
} else {
directed_cor2[i,j] <- pcor2[i,j]
directed_cor2[j,i] <- pcor2[i,j]
}
}
}
}
round(cor2,3) # 相関係数の２乗
round(pcor2,3) # 偏相関係数の２乗
round(directed_cor2,3) #
library(igraph)
#GM1 <- cor2
#GM1 <- pcor2
GM1 <- directed_cor2
diag(GM1) <- 0
GM3 <- GM1*10 # 一番大きな値が10になるように修正
GM3[GM3<1] <- 0 # 1未満の場合は0にする（非表示にするため）
GM4 <- graph.adjacency(GM3,weighted=T, mode = "undirected")
plot(GM4, edge.width=E(GM4)$weight)
以下も進めると、有向相関偏相関分析になります。
for (k in 1:n) {
Data1 <- Data
for (j2 in n:1) {
if(directed_cor2[k,j2] < 0.1){
Data1[,j2] <- NULL
}
}
cor2b <- cor(Data1)^2
pcor2b <- (pcor(Data1)$estimate)^2
nb <- ncol(cor2b)
for(j2 in 1:nb){
if(colnames(Data1[j2]) == colnames(Data[k])){
na <- j2
}
}
if(k == 1){
directed_cor24 <- NULL
}
if(ncol(cor2b) > 1){
directed_cor2b <- cor2b
for (j in 1:nb) {
if(j != na){
directed_cor2b[,j] <- 0
}
}
for (i in 1:nb) {
if(cor2b[i,na] > pcor2b[i,na]){
directed_cor2b[i,na] <- pcor2b[i,na]
}
}
if(min(cor2b) > 0.1){
directed_cor2b[,1:nb] <- 0
}
diag(directed_cor2b) <- 0
directed_cor22 <- data.frame(directed_cor2b)
directed_cor22$Name <- row.names(directed_cor22)
directed_cor23 <- tidyr::gather(directed_cor22, key="X", value = Xs, -Name)
directed_cor24 <- rbind(directed_cor24,directed_cor23)
directed_cor24 <- directed_cor24[!(directed_cor24$Xs < 0.001),]
}
}
directed_cor31 <- directed_cor2
directed_cor31[,1:n] <- 0
if(nrow(directed_cor24) > 1){
for (k in 1:nrow(directed_cor24)) {
for (i in 1:n) {
for (j in 1:n) {
if(directed_cor24[k,1] == rownames(directed_cor31)[i] & directed_cor24[k,2] == colnames(directed_cor31)[j]){
directed_cor31[i,j] <- directed_cor24[k,3]
}
}
}
}
}
directed_cor31[directed_cor31<0.1] <- 0
GM1 <- directed_cor31
GM1 <- ceiling(GM1)
diag(GM1) <- 0
GM5 <- graph.adjacency(GM1,weighted=F, mode = "directed")
plot(GM5)
Data2 <- Data
directed_cor201 <- directed_cor2
directed_cor201[,1:n] <-0
for (k in 1:n) {
Data2[,k] <- (Data[,k] - min(Data[,k]))/(max(Data[,k]) - min(Data[,k]))
}
for (i in 1:n) {
for (j in 1:n) {
if(i !=j){
if(directed_cor2[i,j] > 0.1){
if(directed_cor31[i,j] < 0.1 & directed_cor31[j,i] < 0.1){
directed_cor201[i,j] <- sd(Data2[,j]) / sd(Data2[,i])
if(directed_cor201[i,j] > 1){
directed_cor201[i,j] <- 0
}
if(max(sd(Data2[,j]) , sd(Data2[,i])) < 0.2){ # 有向かを判断できる前提が成り立っていなければ、判断しないようにする
directed_cor201[i,j] <- 0
}
}
} else {
directed_cor201[i,j] <- 0
}
}
}
}
directed_cor201
GM1 <- ceiling(directed_cor201+directed_cor31)
diag(GM1) <- 0
GM6 <- graph.adjacency(GM1,weighted=F, mode = "directed")
plot(GM6)
