2013年1月31日木曜日

共分散構造分析[R編]ってないの?[第3章9]その2

「共分散構造分析[Amos編]-構造方程式モデリング-」69頁、図3.27のファントム変数を含む因果モデルを計算するために、次のテキストファイル company327.txt を用意しました。
ファントム変数は母親の統制的態度因子の符号が反転しただけのものなので、ファントム変数へのパスを-1に固定するだけでなく、ファントム変数そのものの誤差分散を0に固定しています。
#観測方程式
母親の統制的態度 -> v1, NA, 1
母親の統制的態度 -> v2, p2, NA
母親の統制的態度 -> v3, p3, NA
母親の同一化傾向 -> v4, NA, 1
母親の同一化傾向 -> v5, p5, NA
母親の同一化傾向 -> v6, p6, NA
母親の絶対視 -> v7, NA, 1
母親の絶対視 -> v8, p8, NA
母親の絶対視 -> v9, p9, NA
#構成概念間の関係
母親の統制的態度 <-> 母親の同一化傾向, cor1, NA
母親の統制的態度 -> F1, NA, -1 #ファントム変数へのパスを-1に制約 
F1 -> 母親の絶対視, path1, NA  #ファントム変数からのパスは自由に推計

母親の同一化傾向 -> 母親の絶対視, path2, NA
#因子の分散
母親の統制的態度 <-> 母親の統制的態度, var1, NA
母親の同一化傾向 <-> 母親の同一化傾向, var2, NA
母親の絶対視 <-> 母親の絶対視, var3, NA
F1 <-> F1, NA, 0 #ファントム変数の誤差分散は0

そして、次のスクリプトを実行しました。
と、最近のsemパッケージのアップデートで、適合度指標の出力にはおまじないが一行、またはオプションの設定が必要になったようです。一行目がそれです。sem のヘルプのサンプルコードやsummry.semのヘルプを参照してください。
今回はshift-jisが使われているWindowsではなくて、utf-8が標準エンコードのubuntuで作業するので、pathDiagram()もいきなり.dotファイル -> .pdfファイルに出力です。
opt <- options( fit.indices = c("GFI","AGFI","CFI","RMSEA"))
c3mother327.mdl <- specifyModel("model/c3mother327.txt")
summary(c3mother327.sem <- sem(c3mother327.mdl, data=c3mother))
stdCoef(c3mother327.sem)
pathDiagram(c3mother327.sem, edge.labels="values", standardize=TRUE,
            node.font=c("IPAGothic",12), edge.font=c("IPAMincho",10),
            rank.direction="LR", ignore.double=FALSE,
            same.rank=c("母親の統制的態度,母親の同一化傾向"), #並べるノード
            min.rank=c("v1,v2,v3,v4,v5,v6"),  #最も左に配置するノード
            file="graphviz/c3mother327")       #Windowsではshift-jifで処理される


さて、path1とpath2の平均の差の統計量は、推計されたパラメータの分散共分散行列を利用して、
z = (path1 - path2 ) / sqrt ( path1の誤差分散 - 2 * path1とpath2の誤差共分散 + path2の誤差分散)
で計算するらしい。
semパッケージでは推計されたパラメータの分散共分散行列は、
vcov(c3mother.sem) で得られる行列だと思うので、この行列を使って統計量を計算してみましょう。
vc <- vcov(c3mother327.sem)
coeffs <- c3mother327.sem$coef
(coeffs[8] - coeffs[9]) / sqrt(vc[8,8] + vc[9,9] - 2 * vc[8,9])
     path1
-1.276916

本に書かれていた出力と随分違う。んんん。
もちろん、sqrt(diag(vc)[8:9] はsemパッケージのpath1とpath2の標準誤差の出力と一致している。何か大きな勘違いをしているのでしょうか。AMOSの非標準化解はどうなっているんでしょう。いかん!おいらがインチキだってことがバレバレです。semパッケージの標準誤差の計算方法は今後の宿題にしておこう。

図3.28。等値制約を置いたモデルを試してみます。パスに同じ名前をつければ等値制約がかかるはずなので、c3mother327のモデルのpath2をpath1 に書き換えただけのファイルをsem()に計算してもらいました。
非標準化解では、「母親の絶対視 <--- F1」の値だけが表示されています。
標準化解を出力すると、「母親の絶対視 <--- F1」と「母親の絶対視 <--- 母親の同一化傾向」は一致していません。他の推計値も「Amos編」とは微妙に違う。

2013年1月27日日曜日

共分散構造分析[R編]ってないの?[第3章9]その1

「共分散構造分析[Amos編]-構造方程式モデリング-」68頁、図3.26の因果モデルを計算してみます。
次のテキストファイル company326.txt を用意しました。

#観測方程式
母親の統制的態度 -> v1, NA, 1
母親の統制的態度 -> v2, p2, NA
母親の統制的態度 -> v3, p3, NA
母親の同一化傾向 -> v4, NA, 1
母親の同一化傾向 -> v5, p5, NA
母親の同一化傾向 -> v6, p6, NA
母親の絶対視 -> v7, NA, 1
母親の絶対視 -> v8, p8, NA
母親の絶対視 -> v9, p9, NA
#構成概念間の関係
母親の統制的態度 <-> 母親の同一化傾向, cor1, NA
母親の統制的態度 -> 母親の絶対視, path1, NA
母親の同一化傾向 -> 母親の絶対視, path2, NA
#因子の分散
母親の統制的態度 <-> 母親の統制的態度, v1, NA
母親の同一化傾向 <-> 母親の同一化傾向, v2, NA
母親の絶対視 <-> 母親の絶対視, v3, NA

そして、次のスクリプトを実行しました。
c3mother1.mdl <- specifyModel("model/c3mother1.txt")
summary(c3mother1.sem <- sem(c3mother1.mdl, data=c3mother))
stdCoef(c3mother1.sem)
pathDiagram(c3mother1.sem, edge.labels="values", standardize=TRUE,
            node.font=c("IPAGothic",12), edge.font=c("IPAMincho",10),
            rank.direction="LR", ignore.double=FALSE, 
            same.rank=c("母親の統制的態度,母親の同一化傾向"),
            min.rank=c("v1,v2,v3,v4,v5,v6") )

company326.dotファイルは、見栄えのために赤い箇所を書き換えました。
  • headport=nw tailport=sw 矢印が因子に刺さる位置を指定して矢印を曲線にする
  • 終わりの部分はタイトルを表示する

digraph "c3mother1.sem" {
  rankdir=LR;
  size="8,8";
  node [fontname="IPAGothic" fontsize=12 shape=box];
  edge [fontname="IPAMincho" fontsize=10];
  center=1;
  {rank=same "母親の統制的態度" "母親の同一化傾向"}
  "母親の同一化傾向" [shape=ellipse]
  "母親の絶対視" [shape=ellipse]
  "母親の統制的態度" [shape=ellipse]
  "v1" -> "母親の統制的態度" [label="0.8" dir=back];
  "v2" -> "母親の統制的態度" [label="0.68" dir=back];
  "v3" -> "母親の統制的態度" [label="0.6" dir=back];
  "v4" -> "母親の同一化傾向" [label="0.73" dir=back];
  "v5" -> "母親の同一化傾向" [label="0.87" dir=back];
  "v6" -> "母親の同一化傾向" [label="0.7" dir=back];

 
"母親の絶対視" -> "v7" [label="0.88"];
  "母親の絶対視" -> "v8" [label="0.7"];
  "母親の絶対視" -> "v9" [label="0.92"];
  "母親の統制的態度" -> "母親の同一化傾向" [label="0.31" dir=both headport=nw tailport=sw];
  "母親の統制的態度" -> "母親の絶対視" [label="-0.29"];
  "母親の同一化傾向" -> "母親の絶対視" [label="0.37"];
  "母親の統制的態度" -> "母親の統制的態度" [label="1" dir=both];
  "母親の同一化傾向" -> "母親の同一化傾向" [label="1" dir=both];
  "母親の絶対視" -> "母親の絶対視" [label="0.85" dir=both];
  "v1" -> "v1" [label="0.36" dir=both];
  "v2" -> "v2" [label="0.54" dir=both];
  "v3" -> "v3" [label="0.64" dir=both];
  "v4" -> "v4" [label="0.47" dir=both];
  "v5" -> "v5" [label="0.24" dir=both];
  "v6" -> "v6" [label="0.51" dir=both];
  "v7" -> "v7" [label="0.23" dir=both];
  "v8" -> "v8" [label="0.51" dir=both];
  "v9" -> "v9" [label="0.15" dir=both];

  fontname="IPAGothic" ;
  fontsize=15;
  label="構成概念間の因果モデル";
  labelloc="t";

}

以下を実行して、パス図を作成します。
 dot -Tpdf < company325.dot > company326.pdf




「おまじない」なしで、UbuntuにRをインストール

 基本的にWindowsでGUIな私ではありますが、今回購入した廉価なPCにはUbuntuを使ってみることにしました。RもUbuntuで使ってみます。

まず、cranでr-base-core-xxx.debをダウンロードします。すると、プログラム(ソフトウェアセンター)で処理するか、保存するかを聞いてくるので、ソフトウェアセンターにお任せすることにしました。同じように、r-recommended-xxx.debもインストールします。apt-getなどというおまじないは使いません。

 次に、semやggplot2をインストールしようとしました。cranにubuntu用のリンクがないので、Rからインストールします。Rを起動して、

install.packages(c("sem","ggplot2"), dep=TRUE)

なんてことをすると、エラーが出てしまいました。プログラムのコンパイルのためにgfortranやgcc、それに関連するライブラリが必要なようです。
ソフトウェアセンターでsynapticパッケージマネージャーをインストールして起動します。synapticからgfortranやgccを検索して、関連しそうなパッケージをインストールしました。
 これでもエラーがでました。どうも、数値計算関係のライブラリが不足しているようです。blasやlapackのライブラリをsynapticで検索すると、libblasとかliblapackとかいうライブラリがあったので、それもインストールしてみます。
 またまたRに戻って、

install.packages(c("sem","ggplot2"), dep=TRUE)

とやると、今度はエラーもでず、うまく行ったようです。

 乱暴だったけど、「おまじない」 なし、「知識」すらなしで、エラーメッセージを頼りになんとかsemパッケージが動くようになりました。良い子はマネしちゃいけませんね。ちゃんと勉強しましょう。

 RStudioも、ダウンロード-ソフトウェアセンターにお任せでインストールしてしまいました。
-----------------------------------------------------------------------

 時々地図も描くので、rgdalをインストールしようとすると、またもやエラーがでます。そこでsynapticでgdalを検索すると、libgdalというものがあったので、こいつもインストールすると、rgdalもインストールできたようです。

2013年1月19日土曜日

共分散構造分析[R編]ってないの?[第3章8]その2

「共分散構造分析[Amos編]-構造方程式モデリング-」の67頁、図3.25階層因子分析です。
次のテキストファイル company325.txt を用意しました。

#下位因子
V11 <- 自信, NA ,1
V10 <- 自信, s10, NA
V9 <- 自信, s9, NA
V8 <- 自信, s8, NA
V7 <- 自信, s7, NA
V6 <- 自信, s6, NA
V5 <- 自信, s5, NA
V4 <- 自信, s4, NA
V3 <- 自信, s3, NA
V2 <- 自信, s2, NA
V1 <- 自信, s1, NA
#一般因子
V11 <- 仲間関係, NA, 1
V10 <- 仲間関係, r1, NA
V9 <- 仲間関係, r2, NA
V8 <- 身体外見,  NA, 1
V7 <- 身体外見, l1, NA
V6 <- 身体外見, l2, NA
V5 <- 身体外見, l3, NA
V4 <- 身体能力, NA, 1
V3 <- 身体能力, a1, NA
V2 <- 身体能力, a2, NA
V1 <- 身体能力, a3, NA
#因子の分散
仲間関係 <-> 仲間関係, d1, NA
身体外見 <-> 身体外見, d2, NA
身体能力 <-> 身体能力, d3, NA
自信 <-> 自信, d4, NA
#因子の間の相関
仲間関係 <-> 身体外見, cor1
身体能力 <-> 仲間関係, cor2
身体外見 <-> 身体能力, cor3

そして、次のスクリプトを実行しました。
ccompany325.mdl <- specifyModel("model/company325.txt")
summary(company325.sem <- sem(company325.mdl, S = cor.company, N = 100, standardized = TRUE))
stdCoef(company325.sem)
pathDiagram(company325.sem, edge.labels = "values",standardize = TRUE,
            node.font = c("IPAGothic",12), edge.font = c("IPAMincho",10),
            rank.direction = "TB", ignore.double = FALSE,
            same.rank = c("仲間関係,身体外見,身体能力", "自信"))

company325.dotファイルは出力のママだと、下位因子が一般因子が同じ並びになるので、赤い箇所の矢印の向きを反転し、dir=back を書き加えてみました。ついでにラベルを書き加えてタイトルまでおごってしまいます。
digraph "company325.sem" {
  rankdir=TB;
  size="8,8";
  node [fontname="IPAGothic" fontsize=12 shape=box];
  edge [fontname="IPAMincho" fontsize=10];
  center=1;
  {rank=same "仲間関係" "身体外見" "身体能力"}
  "仲間関係" [shape=ellipse]
  "身体外見" [shape=ellipse]
  "身体能力" [shape=ellipse]
  "自信" [shape=ellipse]

  "V11" -> "自信" [label="0.63" dir=back];
  "V10" -> "自信" [label="0.29" dir=back];
  "V9"  -> "自信" [label="0.21" dir=back];
  "V8"  -> "自信" [label="0.51" dir=back];
  "V7"  -> "自信" [label="0.6" dir=back];
  "V6"  -> "自信" [label="-0.03" dir=back];
  "V5"  -> "自信" [label="0.02" dir=back];
  "V4"  -> "自信" [label="0.38" dir=back];
  "V3"  -> "自信" [label="0.64" dir=back];
  "V2"  -> "自信" [label="0.5" dir=back];
  "V1"  -> "自信" [label="0.2" dir=back];


  "仲間関係" -> "V11" [label="0.55"];
  "仲間関係" -> "V10" [label="0.82"];
  "仲間関係" -> "V9" [label="0.61"];
  "身体外見" -> "V8" [label="0.72"];
  "身体外見" -> "V7" [label="0.54"];
  "身体外見" -> "V6" [label="0.72"];
  "身体外見" -> "V5" [label="0.89"];
  "身体能力" -> "V4" [label="0.78"];
  "身体能力" -> "V3" [label="0.45"];
  "身体能力" -> "V2" [label="0.61"];
  "身体能力" -> "V1" [label="0.66"];
  "仲間関係" -> "仲間関係" [label="1" dir=both];
  "身体外見" -> "身体外見" [label="1" dir=both];
  "身体能力" -> "身体能力" [label="1" dir=both];
  "自信" -> "自信" [label="1" dir=both];
  "仲間関係" -> "身体外見" [label="0.69" dir=both];
  "身体能力" -> "仲間関係" [label="0.4" dir=both];
  "身体外見" -> "身体能力" [label="0.32" dir=both];
  "V11" -> "V11" [label="0.3" dir=both];
  "V10" -> "V10" [label="0.24" dir=both];
  "V9" -> "V9" [label="0.58" dir=both];
  "V8" -> "V8" [label="0.22" dir=both];
  "V7" -> "V7" [label="0.35" dir=both];
  "V6" -> "V6" [label="0.48" dir=both];
  "V5" -> "V5" [label="0.21" dir=both];
  "V4" -> "V4" [label="0.24" dir=both];
  "V3" -> "V3" [label="0.38" dir=both];
  "V2" -> "V2" [label="0.38" dir=both];
  "V1" -> "V1" [label="0.52" dir=both];

 fontname = "IPAGothic" ;
 fontsize = 15;
 label = "階層因子分析による結果";
 labelloc = "t"

}

 dot -Tpdf < company325.dot > company325.pdf


で、こんなふうな絵が描けました。

2013年1月16日水曜日

共分散構造分析[R編]ってないの?[第3章8]その1

「共分散構造分析[Amos編]-構造方程式モデリング-」の66頁、2次因子分析です。
モデルがやや複雑になってきたので、モデルファイルを別に作成します。
観測変数の誤差分散は specifyModel() が自動的に補完してくれるので、省略します。
用意したファイル company324.txt の中身は次のとおりです。
specifyModel() でファイルを読み込むときには、ファイルに空白行があっても無視されるし、
#でコメントアウトすることもできるみたい。
#2次因子
仲間関係 <- 仲間評価, NA, 1
身体外見 <- 仲間評価, p1, NA
身体能力 <- 仲間評価, p2, NA
#一次因子
V11 <- 仲間関係, NA, 1
V10 <- 仲間関係, r1, NA
V9<- 仲間関係, r2, NA
V8 <- 身体外見, NA, 1
V7 <- 身体外見, l1, NA
V6 <- 身体外見, l2, NA
V5 <- 身体外見, l3, NA
V4 <- 身体能力, NA, 1
V3 <- 身体能力, a1, NA
V2 <- 身体能力, a2, NA
V1 <- 身体能力, a3, NA

仲間関係 <-> 仲間関係, d3, NA
身体外見 <-> 身体外見, d2, NA
身体能力 <-> 身体能力, d1, NA
仲間評価 <-> 仲間評価, d0, 1

次のスクリプトで計算してもらいます。
company324.mdl <- specifyModel("model/company324.txt")
summary(company324.sem <- sem(company324.mdl, S=cor.company, N=100))
stdCoef(company324.sem)

graphvizで絵を描いてみましょう。誤差分散も出力します。英語以外のフォントを使うときには、フォントの指定をします。
pathDiagram(company324.sem, edge.labels="values",standardize=TRUE,
            node.font=c("IPAGothic",12), edge.font=c("IPAMincho",10),
            rank.direction="TB", ignore.double=FALSE,
            same.rank=c("仲間関係,身体外見,身体能力"))


こんなの出たら、company324.dot に utf-8 のエンコードで保存します。
digraph "company324.sem" {
  rankdir=TB;
  size="8,8";
  node [fontname="IPAGothic" fontsize=12 shape=box];
  edge [fontname="IPAMincho" fontsize=10];
  center=1;
  {rank=same "仲間関係" "身体外見" "身体能力"}
  "仲間関係" [shape=ellipse]
  "身体外見" [shape=ellipse]
  "身体能力" [shape=ellipse]
  "仲間評価" [shape=ellipse]
  "仲間評価" -> "仲間関係" [label="0.97"];
  "仲間評価" -> "身体外見" [label="0.83"];
  "仲間評価" -> "身体能力" [label="0.6"];
  "仲間関係" -> "V11" [label="0.8"];
  "仲間関係" -> "V10" [label="0.81"];
  "仲間関係" -> "V9" [label="0.63"];
  "身体外見" -> "V8" [label="0.91"];
  "身体外見" -> "V7" [label="0.76"];
  "身体外見" -> "V6" [label="0.58"];
  "身体外見" -> "V5" [label="0.73"];
  "身体能力" -> "V4" [label="0.84"];
  "身体能力" -> "V3" [label="0.74"];
  "身体能力" -> "V2" [label="0.8"];
  "身体能力" -> "V1" [label="0.64"];
  "仲間関係" -> "仲間関係" [label="0.07" dir=both];
  "身体外見" -> "身体外見" [label="0.31" dir=both];
  "身体能力" -> "身体能力" [label="0.64" dir=both];
  "仲間評価" -> "仲間評価" [label="1" dir=both];
  "V11" -> "V11" [label="0.37" dir=both];
  "V10" -> "V10" [label="0.34" dir=both];
  "V9" -> "V9" [label="0.6" dir=both];
  "V8" -> "V8" [label="0.18" dir=both];
  "V7" -> "V7" [label="0.43" dir=both];
  "V6" -> "V6" [label="0.67" dir=both];
  "V5" -> "V5" [label="0.47" dir=both];
  "V4" -> "V4" [label="0.3" dir=both];
  "V3" -> "V3" [label="0.45" dir=both];
  "V2" -> "V2" [label="0.36" dir=both];
  "V1" -> "V1" [label="0.59" dir=both];
}

そして、コマンドプロンプトから以下を実行。
dot -Tpdf company324.pdf

サイズの小さいラスター画像は貧相ですが、pdfはうつくしーですね。




2013年1月14日月曜日

共分散構造分析[R編]ってないの?[第3章7]その2

しばらく放置していましたが、約3ヶ月ぶりに再開です。

「共分散構造分析[Amos編]-構造方程式モデリング-」の65頁、図3.22の確認的因子分析です。
因子からのパスを1、誤差分散を0に制約する方法という前者の方法です。
test2.mdl <- specifyModel()
英語学力 -> 英語, NA, 1
学力 -> 国語, p5, NA
学力 -> 地理, p4, NA
学力 -> 数学, p3, NA
学力 -> 化学, p2, NA
学力 -> 物理, p1, NA
学力 <-> 学力, NA, 1
英語学力 <-> 英語学力, NA, 1
英語学力 <-> 学力, c1, NA
国語 <-> 国語, e1, NA
地理 <-> 地理, e2, NA
英語 <-> 英語, NA, 0
数学 <-> 数学, e4, NA
化学 <-> 化学, e5, NA
物理 <-> 物理, e6, NA

summary( test2.sem <- sem( test2.mdl, S = cor( test ), N = nrow( test )))

チラと適合度指標を表示すると
Model Chisquare = 26.604 Df = 10 Pr(>Chisq) = 0.003
GFI = 0.912
AGFI = 0.816
RMSEA = 0.130 90% CI: (0.071, 0.190)
CFI = 0.681
SEMではAGFIを0.9以上とか、RMSEAを0.05以下という基準が受け入れられているようですが、この因子分析は、少し当てはまりが良くないようにみえます。
以前あるセミナーで、因子分析も適合度を確認するので、学生さんが論文を書けなくて大変なんだということを聞きました。実際、アンケートで集めたデータで因子分析してみると、まずいなって思うことがあります。ていうか、ほとんどいつも悲惨。コストかけて、質の悪いデータを集めたと告白すべきでしょうか。

さて、信頼度係数ρ=0.81ということにして、sqrt(ρ)*Sx=9.981をパス、(1-ρ)*Sx^2=23.37を誤差分散として用いて計算してみます。図3.23のパス係数を出すためには、sem()のSオプションには共分散を用いる必要があるようです。(または、data=testとする。)
反対に図3.22の場合は、Sオプションに共分散行列(cov(test))や生データを与えると、 標準化解は図3.22と一致しません。なぜだかはわかりません。教えてください。
test3.mdl <- specifyModel()
英語学力 -> 英語, NA, 9.981
学力 -> 国語, p5, NA
学力 -> 地理, p4, NA
学力 -> 数学, p3, NA
学力 -> 化学, p2, NA
学力 -> 物理, p1, NA
学力 <-> 学力, NA, 1
英語学力 <-> 英語学力, NA, 1
英語学力 <-> 学力, c1, NA
国語 <-> 国語, e1, NA
地理 <-> 地理, e2, NA
英語 <-> 英語, NA, 23.37
数学 <-> 数学, e4, NA
化学 <-> 化学, e5, NA
物理 <-> 物理, e6, NA
 
summary( test2.sem <- sem( test2.mdl, S = cov( test ), N = nrow( test )))
stdCoef( test3.sem )