共分散選択法によるグラフィカルモデリング

共分散選択法により，最適なモデルを作成するためのRスクリプトです（要ggmパッケージ）。

使用法

cov_select(
    model = x, # 分析モデル
    cov = y,   # 共分散行列
    n.obs = 100, # データ数
    AIC = Inf    # 比較対象のAIC（通常は省略）
    cov_orig=NULL # 元の共分散行列（これも指定の必要なし）
)

使用例

x1<-c(8,14,6,14,6,13,6,5,18,10)
x2<-c(8,16,2,8,10,14,6,12,12,10)
x3<-c(8,2,13,7,9,12,10,12,7,8)
x4<-c(7,1,10,8,12,1,9,8,6,14)
x5<-c(3,1,8,9,2,5,9,6,5,1)

COV<-cov(cbind(x1,x2,x3,x4,x5))

# 初期モデルとして全変数間に関係をもたせたモデルを作成
g_model<-UG(~x1*x2*x3*x4*x5)

# 初期モデルの描画
drawGraph(g_model, adjust= FALSE)

# 修正モデルの描画
cov_select(g_model,COV,length(x1))
初期モデルと修正モデル

図：初期モデル（左）と修正モデル（右）

関数定義

cov_select.R

cov_select<-function(model=x, cov=y, n.obs=0, AIC=0, cov_orig=NULL){
  require(ggm)
  if(is.null(cov_orig)){cov_orig<-cov}
#  print(sprintf("AIC= %.4f", AIC ),quote=F)
  #' 偏相関行列を作成
  pmat<- (-1)*solve(cov) / sqrt(diag(solve(cov)) %*% t(diag(solve(cov))))
  diag(pmat)<- 1

  #' 偏相関係数を絶対値に変換
  amat<-abs(pmat)

  #' モデルの係数が0の箇所を無限大に設定
  amat[which(model==0)]<-Inf

  #' 偏相関の絶対値が最小の要素を0にした修正モデルを作成
  model_post<-rep(1,nrow(cov))%*%t(rep(1,nrow(cov)))
  model_post[which.min(amat)]<-0
  model_post<-model_post * t(model_post) * model

  #' モデルのフィットとAICの算出
  f<-fitConGraph(model_post,cov_orig,n.obs)
  AIC_post<-f$dev-2*f$df

  #' モデルの適合度が最大になるまで反復
  if (AIC_post<AIC){
    Recall(model_post,f$Shat,n.obs,AIC=AIC_post,cov_orig=cov_orig)
  }

  #' 最終的に得られたモデルを描画 & 偏相関行列を表示
  else{
    diag(pmat)<-1
    pmat[which(model==0)]<-0
    model.0<-model*0
    f<-fitConGraph(model,cov_orig,n.obs)
    f.0<-fitConGraph(model.0,cov_orig,n.obs)

    # GFIの算出
    S<-cov2cor(cov_orig)
    H<-cov2cor(f$Shat)
    p<-nrow(cov_orig)
    num<-sum(diag((solve(H)%*%(S-H))%*%(solve(H)%*%(S-H))))
    den<-sum(diag((solve(H)%*%S)%*%(solve(H)%*%S)))
    GFI<-1-(num/den)
    AGFI<-1-(p*(p+1)*(1-GFI))/(2*f$df)
    RMSEA<-sqrt(max(((f$dev-f$df)/(f$df*(nrow(cov_orig)-1))),0))
    CFI<-(f.0$dev*f.0$df-f$dev*f$df)/(f.0$dev*f.0$df)
    res<-f
    res<-c(f,GFI=GFI, AGFI=AGFI,RMSEA=RMSEA,CFI=CFI)
    return(list(fit=res,model=model, covmat=pmat))
  }
}