(R) グラフィカルモデリングする
学科の実験?実習?で多変量解析の一環としてRでグラフィカルモデリングをしたのでめも
考え方とかコードとか載せる
まずは概念?の説明
偏相関係数
参考になるやつ 参考にしたやつ - [26-4. 偏相関係数 | 統計学の時間 | 統計WEB] (https://bellcurve.jp/statistics/course/9593.html) - 偏相関係数の意味と式の導出 | 高校数学の美しい物語
自分なりの理解
ようは、2つの統計量に相関係数を計算して、相関関係を測れるが、それは第3の要素が絡んで相関関係が見えてるかもしれないから、その「第3の要素」の影響を消してその2つの相関関係を見ようといことで出てくるのが「偏相関係数」.(調べてみたら1年の時に基礎統計の教科書になってた「入門統計」に載ってた.やっぱり基礎って大事だね!!)
偏相関行列
参考にしたのはこれ.正直全部読んでない.あとで理解する(予定)
4変数以上になってくると、からについての影響を取り除いて、からについての影響を取り除いて...ってやらないと偏相関係数が出ない.
そして、と自体が多変量で、
$$ 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) $$
となっていたら、このについては、偏相関係数を集めた偏相関行列が必要になってくる...ということだと思う.
グラフィカルモデリング
この記事で理解してくれ
ようは、偏相関係数(偏相関行列の要素)が一番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
とかで推定された相関行列が見れたりする.
とりま今はここでさよなら
(コードの説明は今後気が向いたときにする予定)