(R) グラフィカルモデリングする - ナルミヤの備忘録(仮)

学科の実験?実習?で多変量解析の一環としてRでグラフィカルモデリングをしたのでめも
考え方とかコードとか載せる

まずは概念？の説明

偏相関係数

参考になるやつ参考にしたやつ - [26-4. 偏相関係数 | 統計学の時間 | 統計WEB] (https://bellcurve.jp/statistics/course/9593.html) - 偏相関係数の意味と式の導出 | 高校数学の美しい物語

自分なりの理解

ようは、2つの統計量に相関係数を計算して、相関関係を測れるが、それは第３の要素が絡んで相関関係が見えてるかもしれないから、その「第３の要素」の影響を消してその2つの相関関係を見ようといことで出てくるのが「偏相関係数」.(調べてみたら１年の時に基礎統計の教科書になってた「入門統計」に載ってた.やっぱり基礎って大事だね!!)

偏相関行列

参考にしたのはこれ.正直全部読んでない.あとで理解する(予定)

偏相関係数

4変数以上になってくると、 $x$ から $z_1,z_2,...$ についての影響を取り除いて、 $y$ から $z_1,z_2,...$ についての影響を取り除いて...ってやらないと偏相関係数が出ない.
そして、 $x$ と $y$ 自体が多変量で、

$$ x = \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \end{array} \right) ,y = \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right) $$

となっていたら、この $x,y$ については、偏相関係数を集めた偏相関行列が必要になってくる...ということだと思う.

グラフィカルモデリング

この記事で理解してくれ

グラフィカルモデリング

ようは、偏相関係数(偏相関行列の要素)が一番0に近いところを0において、もう一度相関係数と偏相関係数を出して、AICという指標を元に、その0とおく仮定が正しいのかいなかを判断し、OKなら(AICが前のAICより減少してるなら)、また次に0に一番近い偏相関係数を0とおいて....ということを繰り返しやっていき、AICが前のAICより増加した時点でこの操作をストップする.
ちなみに、偏相関行列を算出する際にするこの操作を「共分散選択アルゴリズム」という(らしい).
そして最終的に、算出された偏相関行列を用いて、無向グラフをかくといもの.
ここで、グラフを描くとき、偏相関係数が0のところは線を繋げず、非0なら線を繋げる.

実際のRのコード

必要な関数は一個のファイルにまとめて、#でコメントアウトしてるコードをコンソールでうてば実行してテキストファイルに保存してくれる.

# source("/Users/m-narumiya/OneDrive/ドキュメント/keisu_experiment/cov_selec_algo.R")
# result_list <- cov_sel_algo(df)

# 出力したいファイル先をここで指定
root_path <- "/Users/m-narumiya/OneDrive/ドキュメント/keisu_experiment/result"

#偏相関行列を計算
cor2par <- function(R) { #Rは相関行列
  X <- solve(R)
  p <- nrow(R)
  P <- matrix(0,p,p)
  for (i in (1:p)) {
    for (j in (1:p)) {
      if (i != j) P[i,j] <- -X[i,j]/sqrt(X[i,i]*X[j,j])
      if (i == j) P[i,j] <- 1
    }
  }
  return(P)
}

#消去する(i,j)成分を算出
select_ij <- function(P,amat) { #偏相関行列
  p <- nrow(P)
  minabsP <- Inf
  for (i in (2:p)) {
    for (j in (1:(i-1))) {
      if (amat[i,j] == 1 && abs(P[i,j]) < minabsP) {
        minabsP <- abs(P[i,j])
        i0 <- i
        j0 <- j
      }
    }
  }
  return(c(i0,j0)) #
}


# はじめのamatをセット
set_amat <- function(df) {
  X <- df
  n <- nrow(X)
  p <- ncol(X)
  R <- cor(X)
  amat <- matrix(1,p,p) - diag(p)
  options(digits = 3)
  dimnames(amat) <- dimnames(R)
  return(amat)
}

# amatの更新
renew_amat <- function(amat,c) {
  amat[c[1],c[2]] <- amat[c[2],c[1]] <- 0
  return(amat)
}

#AICの計算のための準備
est_cor_f <- function(df,amat) { #c = sele_ij()
  library("ggm")
  options(digits = 3)
  X <- df
  n <- nrow(X)
  p <- ncol(X)
  R <- cor(X)
  f <- fitConGraph(amat,R,n)
  return(f)
}


# AICの計算
calc_aic <- function(f){
  aic <- f$dev - 2*f$df
  return(aic)
}

# リストのテキストファイルへの出力
print_list <- function(list,out_path) {
  out <- file(out_path,"w")
  for (i in 1:length(list)) {
    writeLines(paste(list[[i]][1]), out, sep=", ")
    writeLines(paste(list[[i]][2]), out, sep="\n")
  }
  close(out) 
}

cov_sel_algo <- function(df) {
  R <- round(cor(df),3)
  P <- round(cor2par(R),3)
  amat <- set_amat(df)
  select_ij_list <- list()
  aic_list <- 0
  last_aic <- 999
  while(1) {
    c <- select_ij(P,amat) #(i,j)選択
    P[c[1],c[2]] <- P[c[2],c[1]] <- 0
    amat <- renew_amat(amat,c) #amatの更新
    f <- est_cor_f(df,amat)
    aic <- calc_aic(f) # aicを計算
    if (aic > last_aic) {
      amat <- last_amat
      aic_list <- c(aic_list,aic)
      P <- estP
      break
    }
    R <- f$Shat #round (f$Shat,3) #推定した相関行列
    P <- cor2par(R) #round(cor2par(R),3) #偏相関行列を生成
    estP <- P;
    aic_list <- c(aic_list,aic)
    select_ij_list <- append(select_ij_list,list(c))
    last_aic <- aic #aicの値を保存
    last_amat <- amat
  }
  aic_list <- c(aic_list)
  select_ij_list
  result_list <- list(round(R,3),round(P,3),select_ij_list,aic_list,amat)
  names(result_list) <- c("cor","cor2par","ijs","aics","amat")
  
  png(paste(root_path,"/plot_graph.png",sep = ""), width = 600, height = 600)  # 描画デバイスを開く
  drawGraph(amat,adjust=FALSE)
  dev.off() # 描画デバイスを閉じる
  
 #推定された相関行列をテキストファイルに保存
  write.table(round(R,3),paste(root_path,"/est_cor.txt",sep = ""),quote=F, col.names=F, append=T)

#推定された偏相関行列をテキストファイルに保存
  write.table(round(P,3),paste(root_path,"/est_cov.txt",sep = ""),quote=F, col.names=F, append=T)

#AICの変化をテキストファイルに保存
  write.table(aic_list,paste(root_path,"/aic_list.txt",sep = ""),quote=F,col.names=F, append=T)

#0にされた成分をテキストファイルに保存
  print_list(select_ij_list,paste(root_path,"/ij_list.txt",sep = ""))

  return (result_list)
}

ちなみに、amatといのは、いわゆる「隣接行列」のこと(だと思う)
ちなみに、テキストファイルを覗かなくても、result_listの中身を見ればresult_list$corとかで推定された相関行列が見れたりする.

とりま今はここでさよなら
(コードの説明は今後気が向いたときにする予定)