myscratchpad

Rで確認的因子分析を行う。
Indiana Universityのサイトとデータを参考にする。感謝。

library(foreign)
dat <- data.frame(read.spss("http://www.indiana.edu/~statmath/stat/all/cfa/values.sav"))

library(sem)
model <- specify.model()
economic   ->  privtown,  p1,   1
economic   ->  govtresp,  p2,   NA
economic   ->  compete,   p3,   NA
moral      ->  homosex,   p4,   NA
moral      ->  abortion,  p5,   NA
moral      ->  euthanas,  p6,   NA
privtown  <->  privtown,  e1,   NA
govtresp  <->  govtresp,  e2,   NA
compete   <->  compete,   e3,   NA
homosex   <->  homosex,   e4,   NA
abortion  <->  abortion,  e5,   NA
euthanas  <->  euthanas,  e6,   NA
economic  <->  economic,  NA,   1
moral     <->  moral,     NA,   1
economic  <->  moral,     em,   NA

semres <- sem(model, cov(dat), nrow(dat))
(sumsem <- summary(semres))
std.coef(semres) # 標準化係数
# このサイトとこのサイトにある関数を利用してAIC, BIC, R二乗を出す。感謝
source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1266106324")
sem.aic(semres) # AIC
rsquare.sem(semres) #  R二乗
sumsem$chisq - sumsem$df * log(semres$N)  # BIC


## 作図。Rgraphvizが上手くインストールできないので (libcdt-4.dllがないとかいわれる) のでGraphvizで書く
 # とりあえずインストールはできた
# 以下のコードで出力されたdigraph...のテキストをGraphvizで描画する。Graphvizを起動し、File -> Newで出たウィンドウにコピペしてF5を押す。設定ダイアログは何もせずにOKとした
path.diagram(semres, ignore.double=FALSE, edge.labels="values", digits=3, standardize=T) # file="filename"と指定するとpdfが作業ディレクトリに出力される。だけど表示がおかしい

出力されたテキスト
# 望んだ作図ではないが、今日はこの辺で勘弁してやる


元サイトと全然結果が一致しない…まいったのう

## specify.model の覚書
変数間の関係、共分散・係数の名前、初期値

変数間の関係は -> か <-> で指定。前後に変数名を入れる
 sem関数に与える相関・共分散行列にない名前は潜在変数になる
その次は共分散の名前。p1とか文字列で指定する。
その次が初期値。NAこれを推定する。
一方向 -> は回帰係数
両方向 <-> で変数名が違う場合は共分散
両方向で変数名が同じ場合はその変数の分散

係数をNAとすると固定母数として扱われる
潜在変数 (因子) の分散は固定母数1にすること
F1 <-> F1, NA, 1
観測変数の分散はNAにする
O1 <-> O1, e1, NA

同じ名前のパスが2つあると同値制約
同じ潜在変数同士　F1 <-F1 は負荷量にNA, 初期値に1を与えて固定母数にする
初心者の方へ: モデル式中にはパス係数、分散、共分散の全てを表記しなければならない。パス係数だけ指定してもモデルが推定されないので、要注意である (Rjiwikiより)

psychパッケージのstructure.sem関数でもモデルが書ける