SPSS ときど記(121〜130）

SPSSを使っていてトラぶったところや変な出力や裏技表技の便利な使い方を中心に書き留めてみる。何回話題があるかわからですが，時々書きます。（Keizo Hori）
最終更新日: (2002/11/25から)
(111)〜(120)　ときど記（メニュー）へ　 (131)〜(140)
  1. SPSS ときど記(130)　2003/03/06　対数線形モデル　loglinear と genlog のモデル指定
  2. SPSS ときど記(129)　2003/03/02　行列言語　Error # 12354：*>添え字が範囲外にあります。
  3. SPSS ときど記(128)　2003/02/04　対数線型モデル　モデル選択
  4. SPSS ときど記(127)　2003/01/30　対数線型モデル　ANOVA codingの巻
  5. SPSS ときど記(126)　2003/01/30　対数線型モデル　dummy codingの巻
  6. SPSS ときど記(125)　2003/01/17　SPSS11.5　インストールできた
  7. SPSS ときど記(124)　2003/01/16　SPSS11.5　インストールがうまくいかない
  8. SPSS ときど記(123)　2002/12/10　因子分析　適合度指標計算 χ2 ULS編
  9. SPSS ときど記(122)　2002/11/28　因子分析　適合度指標計算 χ2
  10. SPSS ときど記(121)　2002/11/28　SPSS　SPSS 11.5 日本語版のリリース

    SPSS ときど記(130)　2003/03/05

    対数線形モデル　loglinear と genlog のモデル指定

    ＳＰＳＳの対数線形モデルのloglinear と genlog にはロジットモデルを計算するようになっている。そしてloglinear はメニューから消えてシンタックスでしか使えなくなっている。メニュー「対数線型」の一般的→genlog, ロジット→genlog, モデル選択→hiloglinear となっている。

    genlogの ロジットモデルの場合、変数指定のところで by を使う。メニューでは従属変数を指定する。「ロジット対数線型分析」のSPSSの説明はヘルプ参照のこと。

    GENLOG
    brand BY soft prev temp
    /MODEL = MULTINOMIAL
    /PRINT = all
    /CRITERIA = CIN(95) ITERATE(20) CONVERGE(.001) DELTA(.5)
    /DESIGN .
    loglinear でも同様である。下の例(1)(2)のようにする。通常の対数線形モデルの場合(3)のようにする。
    反応変数がある場合、ロジットモデルで指定するほうが楽である。下の(1)と等価の対数線形モデル(3)の場合、design で指定が多くなっているのがわかる。
    /DESIGN brand*prev*soft brand*soft*temp brand*temp*prev
    brand*prev brand*soft brand*temp brand.
    これに次の部分が増える。
    soft prev temp soft*prev soft*temp prev*temp soft*prev*temp

    loglinear 同士，またはgenlog同士の結果が同じになることを確かめて下さい。

    ところで、対数線形モデル(loglinear model)、ロジットモデル(logit model)、ロジスティック回帰分析(logistic regression)のことばが紛らわしい。

    (a)対数線形モデルは反応変数と説明変数の区別をしない。ある意味ではすべて反応変数といってもいい。数学的にはシンプルなので、理解しやすいでしょう。Everrit の本でもこれが説明されているのはそのような点を評価したのでしょう。ただし、Everrit(1991)の第２版ではロジットモデルも説明されています。

    (b)ロジットモデルはも反応変数が一つであとは説明変数となる。
    (イ)反応変数はカテゴリー変数，従属変数(これはSPSSの言い方)，応答変数　これはどのような言い方をしても必須
    (ロ)説明変数はカテゴリー変数だけ

    一般には、説明変数に連続変数が含まれる場合にロジスティック回帰分析になる。SPSSのロジットモデルはすべての変数がカテゴリカル変数（質的変数）ということになる。

    ロジスティック回帰分析はメニュー分析では回帰→二項ロジスティック　を選ぶ。

    LOGISTIC REGRESSION VAR=brand
    /METHOD=ENTER soft prev temp
    /CRITERIA PIN(.05) POUT(.10) ITERATE(20) CUT(.5) .

    説明変数に連続変数を含み，反応変数が多項の場合は多項ロジスティック回帰モデルを使う。メニュー分析では回帰→二項ロジスティック　を選ぶ。

    NOMREG
    brand BY soft prev temp
    /CRITERIA CIN(95) DELTA(0) MXITER(100) MXSTEP(5) CHKSEP(20) LCONVERGE(0)
    PCONVERGE(0.000001) SINGULAR(0.00000001)
    /FULLFACTORIAL
    /INTERCEPT INCLUDE
    /PRINT CLASSTABLE FIT PARAMETER SUMMARY LRT .

    ＳＡＳの解説書（Categorical data analysis using SAS system）ではロジスティック回帰分析で両方の一般名称として使っている。

    このあたり，ロジスティック回帰分析を含め日本語で学ぶには Agresti,A. (渡邊裕之ほか訳)(2003).カテゴリカルデータ解析入門サイエンティスト社
    がいいだろう。そのほかの文献については，    対数線形モデル文献案内を見てください。

    他にもGEE なんかもチェックしたほうがいいかな。

    *----------------------------------------------------------.
    TITLE 'Ries and Smith(1963) 新しい洗剤と古い洗剤とどちらを好む'.
    DATA LIST FREE /Soft prev temp brand FREQ.
    VAR LABELS
    Soft '使用水の硬度'/
    prev'以前にＭを使用していたか'/
    Temp '使用温度'/
    brand '新しい洗剤ＸまたはＭを好むか'.
    VALUE LABELS
    Soft 1 '硬水' 2 '中硬水' 3 '軟水'/
    prev 1 '使用なし' 2 '使用あり'/
    Temp 1 '低温' 2 '高温'/
    brand 1 '洗剤 Ｘ' 2 '洗剤 Ｍ'.
    BEGIN DATA.
    1 1 1 1 68
    1 1 1 2 42
    1 2 1 1 37
    1 2 1 2 52
    1 1 2 1 42
    1 1 2 2 30
    1 2 2 1 24
    1 2 2 2 43
    2 1 1 1 66
    2 1 1 2 50
    2 2 1 1 47
    2 2 1 2 55
    2 1 2 1 33
    2 1 2 2 23
    2 2 2 1 23
    2 2 2 2 47
    3 1 1 1 63
    3 1 1 2 53
    3 2 1 1 57
    3 2 1 2 49
    3 1 2 1 29
    3 1 2 2 27
    3 2 2 1 19
    3 2 2 2 29
    END DATA.
    weight by freq.
    *(1).
    loglinear
    brand (1 2) by soft(1 3) prev temp (1 2)
    /PRINT ESTIM
    /PLOT NONE
    /DESIGN brand*prev*soft brand*soft*temp brand*temp*prev
    brand*prev brand*soft brand*temp brand.
    *(2).
    loglinear
    brand (1 2) by soft(1 3) prev temp (1 2)
    /PRINT ESTIM
    /PLOT NONE
    /contrast (soft)=deviation(2)
    /DESIGN brand*prev*soft brand*soft*temp brand*temp*prev
    brand*prev brand*soft brand*temp brand.
    *(3).
    loglinear
    brand (1 2) soft(1 3) prev temp (1 2)
    /PRINT ESTIM
    /PLOT NONE
    /DESIGN brand*prev*soft brand*soft*temp brand*temp*prev
    brand*prev brand*soft brand*temp brand
    soft prev temp soft*prev soft*temp prev*temp soft*prev*temp.

    *(4)Genlog 例も一つ。loglinear と比較して見よう.
    GENLOG
    brand BY soft prev temp
    /MODEL = MULTINOMIAL
    /PRINT = all
    /CRITERIA = CIN(95) ITERATE(20) CONVERGE(.001) DELTA(.5)
    /DESIGN .


    メニューへ　　トップへ　(131)へ

    SPSS ときど記(129)　2003/03/02

    行列言語　Error # 12354：*>添え字が範囲外にあります。


    久しぶりにマクロをいじったら、エラー続出である。そのなかで比較的わかりにくかったエラーメッセージ。

    Error # 12354
    *>添え字が範囲外にあります。

    行列の指示に計算式をそのまま使ったためにでたメッセージだった。
    たとえば、

    print data(#nds-1,:).

    というような行列記述がエラーを起こす。これを

    print data((#nds-1),:).

    とするとエラーは生じない。

    さらにデバッグ中にわかったのだが、loop の中でこのエラーを起こすと、そのままloop の外に飛び出し、end loop の次から実行する。


    メニューへ　　トップへ　(130)へ

    SPSS ときど記(128)　2003/02/04

    対数線型モデル　モデル選択


    対数線形モデルの場合，モデル選択がもっとも重要である。

    (1)χ2値を使う方法であるが，spssのhiloglinearでは偏連合（偏χ2値）(PARTIAL associations)が出力されるのでこれを使う。

    松田(1988,p122)の    表4.8のデータを使った例を示す。部分連合と訳している例もある。

    度数
    C B合計
    12
    1A13176179
    24293297
    合計7469476
    2A117197214
    222325
    合計19220239
    このようなクロス表データを入力するためには，
    変数名　B1 B2として，次のデータをデータ窓に直接入力する。

    3 176
    4 293
    17 197
    2 23

    そして，マクロ「多重クロス表データ→コード＋頻度生成マクロ（http://www.ec.kagawa-u.ac.jp/~hori/spss/spss.html#x2wedata）」をはしらせます。

    x2wedata var=B1 B2/evar=C A/nvar=2 2/cvar=B.


    そうすると次のようにデータが生成される。

    B C A FRQ
    1 1 1 3
    1 1 2 4
    1 2 1 17
    1 2 2 2
    2 1 1 176
    2 1 2 293
    2 2 1 197
    2 2 2 23
    このマクロを使うとそのまま対数線形モデルの処理ができます。
    *weight by frq. /*macro の中で処理済み*/.

    この処理が面倒なら，上の表をデータ窓で打ち込み，
    weight by frq.
    を走らせます。これは重み付けをする処理です。このデータをspssファイルに保存すると重み付けの情報は残っています。

    そのあと次の処理を行います。

    HILOGLINEAR
    a(1 2) b(1 2) c(1 2) /METHOD=BACKWARD
    /CRITERIA MAXSTEPS(10) P(.01) ITERATION(100) CONVERGE(.001) DELTA(0)
    /PRINT=FREQ RESID ASSOCIATION ESTIM
    /DESIGN .

    /CRITERIA P(.01) によって有意水準を１％にしています。既定値は５％水準

    /print association のオプションにより次の結果が出ます。


    * * * * * * * *　H I E R A R C H I C A L 　L O G 　L I N E A R　* * * * * * * *

    Tests that K-way and higher order effects are zero.

    　　　　K 　　DF 　L.R. Chisq　　Prob　Pearson Chisq　　Prob 　Iteration

    　　　　3　　　1 　　　　.043 　.8352 　　　　　.044 　.8338 　　　　　9
    　　　　2　　　4　　　211.482 　.0000　　　　199.646 　.0000 　　　　　2
    　　　　1　　　7 　　1066.428 　.0000 　　　1035.836 　.0000 　　　　　0

    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


    Tests that K-way effects are zero.

    　　　　K 　　DF 　L.R. Chisq　　Prob　Pearson Chisq　　Prob 　Iteration

    　　　　1　　　3　　　854.946 　.0000　　　　836.191 　.0000 　　　　　0
    　　　　2　　　3　　　211.439 　.0000　　　　199.602 　.0000 　　　　　0
    　　　　3　　　1 　　　　.043 　.8352 　　　　　.044 　.8338 　　　　　0



    p値＝.8338となっているので３次の交互作用（A*B*C)は有意でないことがわかる。２次の交互作用までは考慮しないといけない。

    * * * * * * * *　H I E R A R C H I C A L 　L O G 　L I N E A R　* * * * * * * *

    Tests of PARTIAL associations.

    　Effect Name　　　　　　　　　　　　　　　　　　DF　Partial Chisq　　Prob　Iter


    　A*B 　　　　　　　　　　　　　　　　　　　　　　1 　　　　　.039 　.8434 　　2
    　A*C 　　　　　　　　　　　　　　　　　　　　　　1　　　　188.081 　.0000 　　2
    　B*C 　　　　　　　　　　　　　　　　　　　　　　1 　　　　12.173 　.0005 　　2
    　A 　　　　　　　　　　　　　　　　　　　　　　　1　　　　　7.062 　.0079 　　2
    　B 　　　　　　　　　　　　　　　　　　　　　　　1　　　　767.820 　.0000 　　2
    　C 　　　　　　　　　　　　　　　　　　　　　　　1 　　　　80.064 　.0000 　　2

    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    この偏連合から1%水準で考慮すればA*B（p値＝.843)はいらない。松田(1988,P129)表4.10参照。

    /method backward 指定による出力
    高次の交互作用から有意であるかチェックし，有意でなければstepdown していく。最後に残ったものを報告する。


    * * * * * * * *　H I E R A R C H I C A L　 L O G　 L I N E A R　* * * * * * * *

    Backward Elimination (p = .010) for DESIGN 1 with generating class

    　A*B*C

    Likelihood ratio chi square =　　　.00000　　DF = 0　P =　-INF

    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


    If Deleted Simple Effect is　　　　　　　　 DF　 L.R. Chisq Change　　Prob　Iter

    A*B*C　　　　　　　　　　　　　　　　　　　 1　　　　　　　　.043　 .8352　　 9

    Step 1

    　The best model has generating class

    　　　A*B
    　　　A*C
    　　　B*C

    　Likelihood ratio chi square =　　　.04329　　DF = 1　P =　.835

    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


    If Deleted Simple Effect is　　　　　　　　 DF　 L.R. Chisq Change　　Prob　Iter

    A*B　　　　　　　　　　　　　　　　　　　　 1　　　　　　　　.039　 .8434　　 2
    A*C　　　　　　　　　　　　　　　　　　　　 1　　　　　　 188.081　 .0000　　 2
    B*C　　　　　　　　　　　　　　　　　　　　 1　　　　　　　12.173　 .0005　　 2

    Step 2

    　The best model has generating class

    　　　A*C
    　　　B*C

    　Likelihood ratio chi square =　　　.08230　　DF = 2　P =　.960

    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


    If Deleted Simple Effect is　　　　　　　　 DF　 L.R. Chisq Change　　Prob　Iter

    A*C　　　　　　　　　　　　　　　　　　　　 1　　　　　　 193.654　 .0000　　 2
    B*C　　　　　　　　　　　　　　　　　　　　 1　　　　　　　17.746　 .0000　　 2



    * * * * * * * *　H I E R A R C H I C A L　 L O G　 L I N E A R　* * * * * * * *

    Step 3

    　The best model has generating class

    　　　A*C
    　　　B*C

    　Likelihood ratio chi square =　　　.08230　　DF = 2　P =　.960

    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


    * * * * * * * *　H I E R A R C H I C A L　 L O G　 L I N E A R　* * * * * * * *

    The final model has generating class

    　　A*C
    　　B*C


    この表に出てくる
    Likelihood ratio chi square =　　　.08230　　DF = 2　P =　.960
    というのが松田の表4.10のG2， df， p である。

    という訳で残ったのは[AC][BC]と松田の結果と同じになった。このようにどれが反応変数かを考慮しないで階層をたどればいい場合はHiloglenear にまかせればよい。しかし，普通は反応変数が決まっているのでこの処理に任せることはできない。しかし，その場合でも参考資料にはなる。

    (2) AICを使う。spssではAICは出力されないが，χ2値，自由度が出力されるので，計算すればよい。

    上のbackward の出力のstep 3 にLikelihood ratio chi square　が出力されています。このデータを使いましょう。

    Step 3

    　The best model has generating class

    　　　A*C
    　　　B*C

    　Likelihood ratio chi square =　　　.08230　　DF = 2　P =　.960


    松田(1988, p45)(1.34)式を使っているということですね。
    AIC=飽和モデルのAIC−2*df+G2 という計算をしてます。 40.864-2*2+0.082＝36.946

    しかし，「飽和モデルのAIC」は共通なので，実際には坂元慶行(1985, p47)にあるように
    AIC=χ2−2*df　（χ2≒G2）

    を使って求めたほうが簡単である。しかも0以下の値をとると飽和モデルよりもいいということになります。もちろん一番小さいAICとなるモデルを選択する。ただし，２．０以内の違いは意味のある違いと解釈しない方がいい(坂元慶行, 1985)。

    [AB][BC]のAICは
    -2*2+0.082＝-3.918
    となり，飽和モデルよりよい選択であることもわかります。

    このようにAICを計算して，一番小さいモデルを選びます。

    (3)最近，共分散構造分析（SEM）においていろんな適合度指標が開発されている。それらのうちなにか使えるかもしれない。

    また，重回帰分析や因子分析においてAICはあまりよい指標とはいえないという議論も増えてきている。BICにしたほうがいいかもしれないが，このタイプの研究はいまのところみたことがない。

    (4)対数線形モデルと類似の方法として，SASにある proc GENMOD や坂元(1985)のＣＡＴＤＡＰ０２（ｓｐｓｓのマクロ）およびグラフィカルモデリングがある。


    グラフィカルモデリング関連文献
    宮川 雅巳 (1997)グラフィカルモデリング 朝倉書店
    Edwards, D.(2000).Introduction to Graphical Modelling, 2nd Edition. Springer
    catdap02関連文献
    坂元 慶行 (1985)カテゴリカルデータのモデル分析 共立出版（絶版のようです）


    メニューへ　　トップへ　(129)へ

    SPSS ときど記(127)　2003/01/30

    対数線型モデル　ANOVA codingの巻


    SPSS ときど記(126)の続きである。同じデータを使う。
    飽和モデルのANOVAコーディングの推定値は p97　表3.12(3.4　効果量の推定　pdf ファイル)にある。

    今度はloglinear を使用する。
    ところで，anovaコーディングでは基準セルがないのにSPSSではrefcat がある。おかしい。SPSSの既定値は最後のセルである。

    /CONTRAST(a)=deviation
    /CONTRAST(b)=deviation
    を入れてやるとrefcat はa,b とも3になる。
    明示的に
    /CONTRAST(a)=deviation(3)
    /CONTRAST(b)=deviation(3)
    としたほうがあとで問題が起こらないかもしれない。
    *-------------------------.
    title "松田(1988),p70,表3.2 の分析　飽和モデルANOVAコーディング".
    data list free/A B FRQ.
    begin data.
    1 1 5
    1 2 9
    1 3 7
    2 1 7
    2 2 3
    2 3 9
    3 1 8
    3 2 21
    3 3 6
    end data.
    weight by frq.
    LOGLINEAR a(1,3) b(1,3)
    /PRINT = ALL
    /CONTRAST(a)=deviation
    /CONTRAST(b)=deviation
    /CRITERIA=DELTA(0)
    /DESIGN.
    *-------------------------.
    推定値は次の通り。
    Estimates for Parameters

    A

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　1 　-.072480576　　　.18400 　　-.39391 　　-.43313　　　.28817
    　　　　2 　-.242755784　　　.19697　　-1.23247 　　-.62881　　　.14330

    B

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　3 　-.111741588　　　.18540 　　-.60272 　　-.47512　　　.25164
    　　　　4 　.1234483121　　　.18780　　　.65734 　　-.24464　　　.49153

    A BY B

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　5 　-.196344712　　　.27094 　　-.72468 　　-.72738　　　.33470
    　　　　6 　.1562520523　　　.25383　　　.61559 　　-.34125　　　.65375
    　　　　7 　.3104027321　　　.26832 　　1.15683 　　-.21551　　　.83631
    　　　　8 　-.772085028　　　.30669　　-2.51749　　-1.37319 　　-.17098
     genlog よりもさらに不親切になっている。parameter に対応するものはややこしい書き方だ。ま，deviation(x)で指定したｘがあるセルは省略されているということだ。ANOVAコーディングなら出力されるべき，refcat のデータも示されない。前のマニュアル advanced models 11.0およびJの付録 にあったカテゴリカ変数コード化方式 (categorical variables coding scheme)は11.5からはどこを見ればいいのだろう。わからん。

    Parameter 　Aliased　Term

    　　　　1　　　　　　Constant
    　　　　2　　　　　　[A = 1.00]
    　　　　3　　　　　　[A = 2.00]
    　　　　　　　　x　　[A = 3.00]
    　　　　4　　　　　　[B = 1.00]
    　　　　5　　　　　　[B = 2.00]
    　　　　　　　　x　　[B = 3.00]
    　　　　6　　　　　　[A = 1.00]*[B = 1.00]
    　　　　7　　　　　　[A = 1.00]*[B = 2.00]
    　　　　　　　　x　　[A = 1.00]*[B = 3.00]
    　　　　8　　　　　　[A = 2.00]*[B = 1.00]
    　　　　9　　　　　　[A = 2.00]*[B = 2.00]
    　　　　　　　　x　　[A = 2.00]*[B = 3.00]
    　　　　　　　　x　　[A = 3.00]*[B = 1.00]
    　　　　　　　　x　　[A = 3.00]*[B = 2.00]
    　　　　　　　　x　　[A = 3.00]*[B = 3.00]

     であることに気づかなければならない。
    おっと，本当はparameter に対応する表は出力されていた。

    Design Matrix

    　1-A　2-B

    Factor　　　　　　　　　　　　　　　　　　Parameter

    　　1　 2　　　　　　　1　　　　 2　　　　 3　　　　 4　　　　 5

    　　1　 1　　　　 -1.000　　-1.000　　-1.000　　-1.000　　 1.000
    　　1　 2　　　　 -1.000　　-1.000　　 1.000　　　.000　　-1.000
    　　1　 3　　　　 -1.000　　-1.000　　　.000　　 1.000　　　.000
    　　2　 1　　　　　1.000　　　.000　　-1.000　　-1.000　　-1.000
    　　2　 2　　　　　1.000　　　.000　　 1.000　　　.000　　 1.000
    　　2　 3　　　　　1.000　　　.000　　　.000　　 1.000　　　.000
    　　3　 1　　　　　 .000　　 1.000　　-1.000　　-1.000　　　.000
    　　3　 2　　　　　 .000　　 1.000　　 1.000　　　.000　　　.000
    　　3　 3　　　　　 .000　　 1.000　　　.000　　 1.000　　　.000

    　　1　 2　　　　　　　6　　　　 7　　　　 8

    　　1　 1　　　　　1.000　　 1.000　　 1.000
    　　1　 2　　　　　 .000　　-1.000　　　.000
    　　1　 3　　　　 -1.000　　　.000　　-1.000
    　　2　 1　　　　 -1.000　　　.000　　　.000
    　　2　 2　　　　　 .000　　　.000　　　.000
    　　2　 3　　　　　1.000　　　.000　　　.000
    　　3　 1　　　　　 .000　　-1.000　　-1.000
    　　3　 2　　　　　 .000　　 1.000　　　.000
    　　3　 3　　　　　 .000　　　.000　　 1.000


    出力されていないセルはどうしたらいいのか。

    cntrast に refcat を指定して出力するのである。
    /CONTRAST(a)=deviation(1)
    /CONTRAST(b)=deviation(1)
    LOGLINEAR a(1,3) b(1,3)
    /PRINT = ALL
    /CONTRAST(a)=deviation(1)
    /CONTRAST(b)=deviation(1)
    /CRITERIA=DELTA(0)
    /DESIGN.

    Estimates for Parameters

    A

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　1 　-.242755784　　　.19697　　-1.23247 　　-.62881　　　.14330　　　　(2)
    　　　　2 　.3152363604　　　.17208 　　1.83196 　　-.02203　　　.65250　　　　(3)

    B

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　3 　.1234483121　　　.18780　　　.65734 　　-.24464　　　.49153　　　　(2)
    　　　　4 　-.011706724　　　.18062 　　-.06481 　　-.36572　　　.34231　　　　(3)

    A BY B

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　5 　-.772085028　　　.30669　　-2.51749　　-1.37319 　　-.17098　　　　(2,2)
    　　　　6 　.4616822963　　　.25830 　　1.78737 　　-.04459　　　.96796　　　　(2,3)
    　　　　7 　.6158329761　　　.23049 　　2.67184　　　.16407 　　1.06759　　　　(3,2)
    　　　　8 　-.501774956　　　.25240　　-1.98799 　　-.99649 　　-.00706　　　　(3,3)
    同様に
    2も指定して，a by b のところを完成させる。

    /CONTRAST(a)=deviation(2) /CONTRAST(b)=deviation(2)
    LOGLINEAR a(1,3) b(1,3)
    /PRINT = ALL
    /CONTRAST(a)=deviation(2)
    /CONTRAST(b)=deviation(2)
    /CRITERIA=DELTA(0)
    /DESIGN.


    Estimates for Parameters

    A

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　1 　-.072480576　　　.18400 　　-.39391 　　-.43313　　　.28817　　　　(1)
    　　　　2 　.3152363604　　　.17208 　　1.83196 　　-.02203　　　.65250　　　　(3)

    B

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　3 　-.111741588　　　.18540 　　-.60272 　　-.47512　　　.25164　　　　(1)
    　　　　4 　-.011706724　　　.18062 　　-.06481 　　-.36572　　　.34231　　　　(3)

    A BY B

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　5 　-.196344712　　　.27094 　　-.72468 　　-.72738　　　.33470　　　　(1,1)
    　　　　6 　.0400926601　　　.25556　　　.15688 　　-.46080　　　.54099　　　　(1,3)
    　　　　7 　-.114058020　　　.24663 　　-.46246 　　-.59746　　　.36934　　　　(3,1)
    　　　　8 　-.501774956　　　.25240　　-1.98799 　　-.99649 　　-.00706　　　　(3,3)
    　　　
    これを継ぎ合わせて，松田(1988)表3.12のようにする。

    こんな面倒なこと誰がするのだろう。SPSSはANOVAコーディングの意味を理解できていないのだろうか。

    心理学のデータの場合，ANOVAコーディングのほうが解釈しやすいことが多いであろう。ダミーコーディングの場合，仮説がより鮮明なときに使うことになる。
    hiloglinear もANOVAコーディングと同じ結果となる。ただし，飽和モデルのみ推定値を出力する。

    *-----------------------.
    title "松田(1988),p70,表3.2 の分析　飽和モデルhilog".
    recode a,b(1=4).
    HILOGLINEAR a(2,4) b(2,4)
    /PRINT = ALL
    /method=backward
    /CRITERIA=DELTA(0) CONVERGE(0.001) ITERATE(100)
    /DESIGN
    .
    *-----------------------.
    最後のカテゴリー（セル）の出力が省略される．あらかじめ，最初のカテゴリーを最後に移動した。
    結果は次のようになる。
    Estimates for Parameters.

    A*B

    　Parameter　　　Coeff. 　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　1　-.7720850284　　　.30669　　-2.51749　　-1.37319 　　-.17098
    　　　　2 　.4616822963　　　.25830 　　1.78737 　　-.04459　　　.96796
    　　　　3 　.6158329761　　　.23049 　　2.67184　　　.16407 　　1.06759
    　　　　4　-.5017749564　　　.25240　　-1.98799 　　-.99649 　　-.00706

    A

    　Parameter　　　Coeff. 　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　1　-.2427557841　　　.19697　　-1.23247 　　-.62881　　　.14330
    　　　　2 　.3152363604　　　.17208 　　1.83196 　　-.02203　　　.65250

    B

    　Parameter　　　Coeff. 　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　1 　.1234483121　　　.18780　　　.65734 　　-.24464　　　.49153
    　　　　2　-.0117067240　　　.18062 　　-.06481 　　-.36572　　　.34231


    using an iterative proportional-fitting algorithm.とはあるがコーディングについては特に書いていない。このアルゴリズムを言っただけで分かるのか。algorithmのファイルを調べなくては。といっても実り少ないのですぐにはやる気がおこらない。

    hiloglinear のいいところは，偏連合を出すことができる。

    * * * * * * * *　H I E R A R C H I C A L 　L O G 　L I N E A R　* * * * * * * *

    Tests that K-way and higher order effects are zero.

    　　　　K 　　DF 　L.R. Chisq　　Prob　Pearson Chisq　　Prob 　Iteration

    　　　　2　　　4 　　　11.299 　.0234 　　　　10.466 　.0333 　　　　　2
    　　　　1　　　8 　　　20.874 　.0075 　　　　25.200 　.0014 　　　　　0

    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


    Tests that K-way effects are zero.

    　　　　K 　　DF 　L.R. Chisq　　Prob　Pearson Chisq　　Prob 　Iteration

    　　　　1　　　4　　　　9.575 　.0482 　　　　14.734 　.0053 　　　　　0
    　　　　2　　　4 　　　11.299 　.0234 　　　　10.466 　.0333 　　　　　0

    ２次効果（A*B)，１次効果(a,b)とも５％水準で有意である。

    * * * * * * * *　H I E R A R C H I C A L 　L O G 　L I N E A R　* * * * * * * *

    Tests of PARTIAL associations.

    　Effect Name　　　　　　　　　　　　　　　　　　DF　Partial Chisq　　Prob　Iter


    　A 　　　　　　　　　　　　　　　　　　　　　　　2　　　　　5.802 　.0550 　　2
    　B 　　　　　　　　　　　　　　　　　　　　　　　2　　　　　3.773 　.1516 　　2

    しかし，主効果（１次の項）A,B を単独で調べると5%水準で有意ではない。つまり，１次効果はA,Bあわせて５％水準に届いたのである。

    SASのgenmod も変数が有意かどうか調べることができる。が，今SASを使える環境になっていないのでこのデータのチェックができない。誰か確かめて欲しい。

    hiloglinear はこのほかstepdown の処理ができる。

    hiloglinearとloglinear にちてはStevens(2001)の１４章参照のこと。この本はまだgenlog の説明がない。

    PS. 2003.2.2
    hiloglinear について書いた件に間違いがあったので，一部削除した。カテゴリーの値の範囲が１〜３でもちゃんとANOVAコーディングの出力がされた。

    また，出力はドラフトビューに出力する方が結果は見やすい。
    ドラフトビューにするには
    ファイル→新規作成→ドラフト出力

    ANOVAコーディングの場合，カテゴリが２つの場合，一方は省略できる。これは交互作用の場合でも同じ。カテゴリが３以上のときは省略できない。spss のような出力は納得できない。ecta はカテゴリが２つの場合のみ省略している。どうせなら全部出力してほしいが。


    メニューへ　　トップへ　(128)へ

    SPSS ときど記(126)　2003/01/30

    対数線型モデル　dummy codingの巻


    対数線型モデルについてはNifty のspss フォーラムに随分書き込んだ。
    また，文献案内は    対数線形モデル文献案内に書いた。

    質的情報の多変量解析 (松田 紀之 著)のデータについての質問があった( [fpr 2308] SPSSでの対数線形モデルの出力の見方)ので，このデータを使って分析してみよう。 　p70の表3.2 ロケット発射実験のデータ(pdf ファイル)である。

    ロケットの発射試験（松田,1988,p70,表3.2)（度数）
    横方向ずれ（ヤード） 合計
    B1(-250,-51) B2(-50-+49) B3(+50,+199)
    航続距離（ヤード） a1(0-1199) 5 9 7 21
    a2(1200-1799) 7 3 9 19
    a3(1800-2699) 8 21 6 35
    合計 20 33 22 75

    基準セル(spss ではrefcat 参照カテゴリーか？) を(1,1)にした場合の飽和モデルのダミーコーディングの推定値は p97　表3.11(3.4　効果量の推定　pdf ファイル)にある。

    ついでなので、度数データを入れて、genlogの処理をするシンタックスを記す。

    (a)松田の推定値においては基準セル(1,1)であるが、ＳＰＳＳの既定値は(最後のセル,最後のセル）なので基準セル a=1,b=1 を最後のセルにする。
    recode a,b(1=4).

    (b) /CRITERIA = CIN(95) ITERATE(20) CONVERGE(.001) DELTA(0)
    CRITERIA において delta(0)としている。delta は度数に加える数を意味する。このdelta は0のセルがある場合に対数がとれなくなるのを防ぐものである。このデータでは0のセルはない。0にしておかないと松田の結果と比較できない。他の「CIN(95) ITERATE(20) CONVERGE(.001)」はメニューから貼り付けたら出力されたものである。既定値なのでなくてもいい。ITERATE, CONVERGEはITERATE(100) CONVERGE(.0001)と反復回数を増やし、収束基準を厳しくしたほうがいいだろう。このデータでは大した問題ではない。

    (c) /MODEL = MULTINOMIAL
    genlog の既定値はポアソン分布(/MODEL = POISSON)なので松田の処理と合わせるために多項分布（MULTINOMIAL）を指定しなければならない。

    次のシンタックスをシンタックス窓に入れて走らせてみよう。

    *------------------------------.

    title "松田(1988),p70,表3.2 の分析　飽和モデルダミーコーディング".
    DATA LIST FREE/A B freq.
    BEGIN DATA
    1 1 5
    1 2 9
    1 3 7
    2 1 7
    2 2 3
    2 3 9
    3 1 8
    3 2 21
    3 3 6
    END DATA.

    WEIGHT BY freq.
    recode A,B(1=4).
    GENLOG
    a b
    /MODEL = MULTINOMIAL
    /PRINT = FREQ ESTIM CORR COV
    /CRITERIA = CIN(95) ITERATE(20) CONVERGE(.001) DELTA(0)
    /DESIGN .
    *------------------------------. 出力において次の表があるのに注意しよう。番号と効果の対応が記されている。xがあるセルは基準にしているので推定値は0である。

    Correspondence Between Parameters and Terms of the Design

    Parameter 　Aliased　Term

    　　　　1　　　　　　Constant
    　　　　2　　　　　　[A = 2.00]
    　　　　3　　　　　　[A = 3.00]
    　　　　4 　　　x　　[A = 4.00]
    　　　　5　　　　　　[B = 2.00]
    　　　　6　　　　　　[B = 3.00]
    　　　　7 　　　x　　[B = 4.00]
    　　　　8　　　　　　[A = 2.00]*[B = 2.00]
    　　　　9　　　　　　[A = 2.00]*[B = 3.00]
    　　　10 　　　x　　[A = 2.00]*[B = 4.00]
    　　　11　　　　　　[A = 3.00]*[B = 2.00]
    　　　12　　　　　　[A = 3.00]*[B = 3.00]
    　　　13 　　　x　　[A = 3.00]*[B = 4.00]
    　　　14 　　　x　　[A = 4.00]*[B = 2.00]
    Parameter 　Aliased　Term
    　　　15 　　　x　　[A = 4.00]*[B = 3.00]
    　　　16 　　　x　　[A = 4.00]*[B = 4.00]

    Note: 'x' indicates an aliased (or a redundant) parameter.
    　　　These parameters are set to zero.

    推定値の出力は
    　　　　　　　　　　　　　　　　　　　　　　　Asymptotic 95% CI Parameter 　Estimate 　　　　SE　　Z-value　　　Lower　　　Upper

    　　　　2　　　.3365　　　.5855　　　　.57 　　　-.81 　　　1.48
    　　　　3　　　.4700　　　.5701　　　　.82 　　　-.65 　　　1.59
    　　　　4　　　.0000　　　.　　　　　　.　　　　　.　　　　　.
    　　　　5　　　.5878　　　.5578 　　　1.05 　　　-.51 　　　1.68
    　　　　6　　　.3365　　　.5855　　　　.57 　　　-.81 　　　1.48
    　　　　7　　　.0000　　　.　　　　　　.　　　　　.　　　　　.
    　　　　8　　-1.4351　　　.8873　　　-1.62　　　-3.17　　　　.30
    　　　　9 　　-.0852　　　.7725 　　　-.11　　　-1.60 　　　1.43
    　　　10　　　.0000　　　.　　　　　　.　　　　　.　　　　　.
    　　　11　　　.3773　　　.6955　　　　.54 　　　-.99 　　　1.74
    　　　12 　　-.6242　　　.7966 　　　-.78　　　-2.19　　　　.94
    　　　13　　　.0000　　　.　　　　　　.　　　　　.　　　　　.
    　　　14　　　.0000　　　.　　　　　　.　　　　　.　　　　　.
    　　　15　　　.0000　　　.　　　　　　.　　　　　.　　　　　.
    　　　16　　　.0000　　　.　　　　　　.　　　　　.　　　　　.
    となり，松田(1988) p97表3.11と一致する。

    ところで，loglinear でも出るはずと思っていたが，結果が違っている。

    /CONTRAST(a)=simple(1)
    /CONTRAST(b)=simple(1)
    loglinear の場合refcat を指定できる。上の１のところである。
    *-----------------------.
    title "松田(1988),p70,表3.2 の分析　飽和モデルダミーコーディング".
    data list free/A B FRQ.
    begin data.
    1 1 5
    1 2 9
    1 3 7
    2 1 7
    2 2 3
    2 3 9
    3 1 8
    3 2 21
    3 3 6
    end data.
    weight by frq.
    LOGLINEAR a(1,3) b(1,3)
    /PRINT = ALL
    /CONTRAST(a)=simple(1)
    /CONTRAST(b)=simple(1)
    /CRITERIA=DELTA(0) CONVERGE(0.0001) ITERATE(100)
    /DESIGN.
    *-----------------------.
    結果は

    Estimates for Parameters

    A

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　1 　-.170275208　　　.34014 　　-.50060 　　-.83695　　　.49640
    　　　　2 　.3877169366　　　.29688 　　1.30596 　　-.19417　　　.96961

    B

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　3 　.2351899002　　　.32659　　　.72015 　　-.40492　　　.87530
    　　　　4 　.1000348642　　　.31420　　　.31838 　　-.51580　　　.71587

    A BY B

    　Parameter 　　Coeff.　　Std. Err. 　　Z-Value Lower 95 CI Upper 95 CI

    　　　　5 　-1.43508453　　　.88730　　-1.61736　　-3.17419　　　.30402
    　　　　6 　-.085157808　　　.77254 　　-.11023　　-1.59935 　　1.42903
    　　　　7 　.3772942311　　　.69551　　　.54247 　　-.98590 　　1.74049
    　　　　8 　-.624154309　　　.79657 　　-.78355　　-2.18543　　　.93712

    A BY Bの４つはgenlog と一致しているが，主効果はA,Bとも違っている。なんでだろう。

    refarence manual(p631 ,ver 11.5用) によると。
    loglinear の場合
    (1)再パラメタ化法をとっている。(genlog の場合は回帰）。デザインが冗長の場合，解釈に困る。
    (2)計算結果が間違っていることがあるとしている。しかし，不飽和モデルのときの自由度がらみの問題だからこれと違っている。

    再パラメタ化法が原因だろうか？？　確かに松田のダミーコーディングの例とSPSSのdesign は違っている。

    spss の design matrix は次のようになる。松田(1988, p82-p84)のダミーコーディングと比較してみよう。付録p200-203も参照のこと
    *-------------------------------
    Design Matrix

    　1-A　2-B

    Factor　　　　　　　　　　　　　　　　　　Parameter

    　　1 　2　　　　1 　　　　2 　　　　3 　　　　4 　　　　5 　　　　6 　　　　7 　　　　8

    　　1 　1　　 -.333 　　-.333 　　-.333 　　-.333　　　.111　　　.111　　　.111　　　.111
    　　1 　2　　 -.333 　　-.333　　　.667 　　-.333 　　-.222　　　.111 　　-.222　　　.111
    　　1 　3　　 -.333 　　-.333 　　-.333　　　.667　　　.111 　　-.222　　　.111 　　-.222
    　　2 　1　　　.667 　　-.333 　　-.333 　　-.333 　　-.222 　　-.222　　　.111　　　.111
    　　2 　2　　　.667 　　-.333　　　.667 　　-.333　　　.444 　　-.222 　　-.222　　　.111
    　　2 　3　　　.667 　　-.333 　　-.333　　　.667 　　-.222　　　.444　　　.111 　　-.222
    　　3 　1　　 -.333　　　.667 　　-.333 　　-.333　　　.111　　　.111 　　-.222 　　-.222
    　　3 　2　　 -.333　　　.667　　　.667 　　-.333 　　-.222　　　.111　　　.444 　　-.222
    　　3 　3　　 -.333　　　.667 　　-.333　　　.667　　　.111 　　-.222 　　-.222　　　.444
    *----------------------. spssには各自がコーディングするオプションもあるのでそれで記述すればいいということになる。が，genlog でできるのでそこまでやる必要はない。

    メニューへ　　トップへ　(127)へ

    SPSS ときど記(125)　2003/01/17

    SPSS11.5　インストールできた


    SPSS11.5日本語版が来たがインストールがうまくいかないと前の    SPSS ときど記(124)で書いたが，できました。

    SPSS社に問い合わせのメールをうったところ，質問が返ってきた。

    1.　お使いのPCのメモリ容量と、SPSSインストール前のハードディスクの空き容量、CPUのスペックを教えてください。以下の方法で確認できます。（略）

    2. 下記インストール1)を試したときに、SPSS製品をインストールするフォルダが作成されていると思います。続けて2) 3)を試すにあたり、前のインストールで作成されたフォルダを削除したり、または、アンインストールを試すような処理を行っていますでしょうか。

    3.インストールに際し、ウィルス対策ソフトなどの常駐プログラムや、他のアプリケーションは稼動していなかったでしょうか。

    4.インストールは、ローカルPCにAdministrator権限のあるユーザで行っているでしょうか。

    １については，253.0MB，20.6GB
    ２については，アンインストールはしていません。
    ３については，atokは常駐してました。
    と答えたがこれは大嘘だった。ウィルスバスターやNetscape7Jなどが常駐していた。
    後で試したところ，Netscape7 はかなりでかいメモリーを喰っていた。
    ４についてはWindows 98なので関係ないだろう。

    さらに「卓駆」というファイラーで調べたところ，最初のAtok15 が立ち上がる当たりで他のソフトを起動すると，メチャ　CPUメモリを失うことがわかった。私のパソコンの設定において正常なら83％〜85%使用だが，93%以上のメモリを使用する。原因の第１はこれだろう。
    いくつか実験したところ，どうもメモリーの解放に時間がかかるものがあるらしく，結果は不安定である。

    ある実験の場合
    (1)起動＋卓駆 85%使用
    (2)Dell touch をはずす　84% (ちょっと怪しい他の場合には2%変動することがある）
    (3)Netscape7.0 をはずす　77%（大きすぎる。2%しか減らないときもあった。なお，常駐をはずすだけで，Netscape そのものは立ち上げていない）
    (4)Real Player 76%
    (5)スーパーニッポニカライト 75%
    (6)Book Shelf 2.0クイック　74%
    (7)ウィルスバスター2002 73%

    というわけで，Netscape 7.0, と Dell Touch(こいつは何かに使っているのか？）が最悪のようだ。

    ところで，常駐のものをはずさす，Windows 98がきちんと起動し終わるまって，Base だけをとりあえずインストールしたところ，正常にインストールできた。
    その後，Netscape 7などの常駐を一部はずして，他をインストールしたところ正常にインストールできた。とりあえず分析をする。

    ところで，11.5から紙のマニュアルをなくすということだ。それらはhelp で参照できるということだ。だが，日本語化していなかった 統計の説明をしている User Manualはhelp には見あたらないようだがどうなっているの。ケーススタディがそれだというのだろうか。

    きちんとした使用例がないと分からない場合も結構あるのだけどね。それともalgorithm の参考文献を読みまくれということだろうか。

    なんか，ますますSPSS難民が増えそうですね。SAS難民はどこにいけばいいの。

    ps. 2003年1月29日
    ウィルスバスター2002を2003に入れ替えたためか，メモリーが随分と楽になった。73%に減っている。常駐のnetscape 7.0をはずした状態だと，63%になった。もちろん条件によっていろいろ変わるが，最小値はこの数字のようだ。

    メニューへ　　トップへ　(126)へ

    SPSS ときど記(124)　2003/01/16

    SPSS11.5　インストールがうまくいかない


    SPSS11.5日本語版が来たがインストールがうまくいかない

    Windows 98　を使ってます。
    SPSSの11.5日本語版をインストールに問題が生じます。いくつかやってそれぞれ反応が違っています。

    (1)カスタムで持っているのを全部インストールしてみようとする。
    (a)BASE　(b)Regression　(c)advanced (d)exact test (e)categories
    の５つ。すべてをインストールすることにする。
    インストール中 77%のところ(Fatsurface.save ？「Fat surfactant.sav？か」とかのファイルのところ）でトラブル発生。
    ×マークだけのダイアログ。×で対応。

    すると「メモリーが足りないので展開（解凍？）できない」とのエラーメッセージのメッセージ。無視を押す。何度かこの繰り返し。
    その後インストールを継続する。

    11.5を立ち上げようとすると，「spssprnt.dll が壊れてます」のダイアログボックスがでて，立ち上がらない。以下のインストールでも同じ。

    (2)BASEだけをカスタムでインストール。前の文句はなく，スムースにいくが，立ち上げようとすると，「spssprnt.dll が壊れてます」のダイアログボックスがでて，立ち上がらない。

    (3)標準を選ぶ。インストール中。すぐに「0による割り算がスクリプト中で検出されました。エラー101」のダイアログ。４回。スムースにいくが，立ち上げようとすると，「spssprnt.dll が壊れてます」のダイアログボックスがでて，立ち上がらない。

    PS.解決編は    (125)

    メニューへ　　トップへ　(125)へ

    SPSS ときど記(123)　2002/12/10

    因子分析　適合度指標計算 χ2 ULS編

    SPSS ときど記(122)の続きである。uls（最小２乗法）がもっとも簡単なはずなのに意外と手間取った。EQS, AMOS の出力と一致しないからである。spss，Mplus では適合度を出力しない。SteigerがSTATISTICA の日本語マニュアルp3579-3580でも望ましくないと言っているように多くのものは統計学的に望ましくないことを指摘している。EQS ではそのためχ2を求めるのに別の計算法を使っている(EQS program manual1995, p217)。複雑なのでよくわからないが，値は大きくなるようなので有意になりにくいのか。

    AMOS はマニュアルには特に何もしていないようだが(Amos 4.0 User's Guid, p392)，値が小さい。Bartlett の修正をしているのだろうか？しかし，Bartlett の修正した値とも微妙に違う。

    macro は素直である。

    *---------------------------. *適合度 uls.
    matrix.
    compute S={1.0, .3, .2, .1;
    .3,1., .2, .2;
    .2, .2,1., .3;
    .1, .2, .3,1}.
    print S/title "元の相関行列 S"/format f7.3.
    * 'M: 変数の数'.
    compute m=nrow(S).
    print m/title "M: 変数の数".

    compute N=240.
    print N/title "N サンプルサイズ".
    * "１因子解のときの重み".
    *spss の結果.
    compute A={.42541;
    .50750;
    .50750;
    .42541
    }.

    print A/title "１因子解のときの重み係数 A"/format f7.3.
    * "重みから再生相関行列を作る。".
    compute nfac=ncol(A).
    print nfac/title "nfac: 因子数"/format f5.0.
    compute C=A*t(A).
    loop i=1 to m.
    compute C(i,i)=1.
    end loop.

    print C/title "C: 重みから再生相関行列を作る。"/format f7.3.
    *ULS.
    compute d0=(S-C).
    print d0/title "d0: 残差行列"/format f7.3.
    compute Fk=trace(d0*d0)/2.
    print Fk/title 'Fk: F 式(1) (ULS)'/format f7.5.
    compute chi2=Fk*(N-1).
    print chi2/title 'Chi2: Chi square (ULS)'/format f7.3.
    compute df=((m-nfac)**2-(m+nfac))/2.
    print df/format f7.0.
    compute p=1-CHICDF(chi2,df).
    print p/title "p: p-value"/format f7.4.
    end matrix.
    *---------------------------.

    出力

    Run MATRIX procedure:

    元の相関行列 S
    　1.000　　.300　　.200　　.100
    　　.300 　1.000　　.200　　.200
    　　.200　　.200 　1.000　　.300
    　　.100　　.200　　.300 　1.000

    M: 変数の数
    　4

    N サンプルサイズ
    　240

    １因子解のときの重み係数 A
    　　.425
    　　.507
    　　.507
    　　.425

    nfac: 因子数
    　　1

    C: 重みから再生相関行列を作る。
    　1.000　　.216　　.216　　.181
    　　.216 　1.000　　.258　　.216
    　　.216　　.258 　1.000　　.216
    　　.181　　.216　　.216 　1.000

    d0: 残差行列
    　　.000　　.084 　-.016 　-.081
    　　.084　　.000 　-.058 　-.016
    　-.016 　-.058　　.000　　.084
    　-.081 　-.016　　.084　　.000

    Fk: F 式(1) (ULS)
    　.02452

    Chi2: Chi square (ULS)
    　5.861

    DF
    　　　2

    p: p-value
    　.0534

    ------ END MATRIX -----
    EQS のχ2=9.302
    Amosのχ2=5.812
    Bartlettの修正しているものとして元に戻すと
    5.812/(240-1-(2*4+5)/6-2*1/3)*(240-1)=5.881727594

    メニューへ　　トップへ　(124)へ

    SPSS ときど記(122)　2002/11/28

    因子分析　適合度指標計算 χ2

    SEMの影響を受けて探索的因子分析でも適合度を求める機会がでてきた。すでに一部スクリプトは    発表している。服部環さん＠筑波大学の因子分析ソフトでは多くの指標を求めることができる。これをexcel上のデータを使って，データ処理を行い，結果をexcelに読みとって，見やすく表示するvba を開発中。β版はできた。あとは使用説明書をつくるまでになっている。

    それだけでは淋しいので，SPSSのmatrix 言語を使って適合度指標を求めるプログラムを少し作ってみた。

    まずは基本中の基本χ2を求める。χ2値を求める前の距離（乖離度）つまり適合度関数（F)を Loehlin(1998)に従って示す。

    一般的には次の式で表すことができる。

    F=1/2 trace((S-C)V)2　…(1)

    S は測定されたデータの相関行列（R）または共分散行列（S)　ここではS
    C は推測されたモデルから再生された相関行列または共分散行列　Σ，Σ(θ＾)などと表される。

    ここでは簡単のため相関行列で考える。

    因子分析の因子パタン行列をAとすると

    C=(AA' の対角に１を入れる）　……(2)　 ('は転置を示す　SPSSではT）


    Vは次のようになる。

    最小２乗法
    OLS　のときは　V=I　　…(3)
    重み付き最小２乗法
    GLS　のときは　V=S-1　…(4)
    最尤法
    ML 　のときは　V=C-1　…(5)

    (1)式により３つの方法を統一して理解できることを示している。S-Cつまり残差行列に重みをつける付け方が、OLS,GLS、mlによって違う。OLSは単位行列をかけているのでそのまま。そのtrace なので、残差行列の非対角要素の２乗和÷２つまり、下三角の非対角要素の乗和となる。GLSは標本相関行列によって重みをつけている。(5)の方法はRLS（Reweighted Least Squares)と呼ばれている。

    MLの場合実際は

    F=ln|C|-ln|S|+trace(SV-1)　…(6)

    を使って求める。(1)式は独立モデルのときに使われる。また、(1)の式と(6)式は値は多少異なるが、比例関係にある。

    χ2=F*(N-1)　　　…(7)

    N=サンプルサイズ

    としてχ2値が求められる。

    MLの場合


    *--------------------------------.
    *適合度χ2　 ml.
    matrix.
    compute S={1.0, .3, .2, .1;
    .3,1., .2, .2;
    .2, .2,1., .3;
    .1, .2, .3,1}.
    print S/title "標本相関行列 S"/format f7.3.
    compute m=nrow(S).
    print m/title "M: 変数の数".

    compute N=240.
    print N/title "N サンプルサイズ".
    * "１因子解のときの重み".
    compute A={ 0.430 ;
    0.502 ;
    0.502 ;
    0.430}.
    print A/title "１因子解のときの重み係数 A"/format f7.3.
    * "重みから再生相関行列を作る。".
    compute C=A*t(A).
    loop i=1 to m.
    compute C(i,i)=1.
    end loop.

    print C/title "C: 重みから再生相関行列を作る。"/format f7.3.
    *ML.
    compute Fk=ln(det(C))-ln(det(S))+trace(S*inv(C))-m.
    print Fk/title 'Fk 適合度関数：本来の方法(6)'/format f7.4.
    compute chi2=Fk*(N-1).
    print chi2/title 'Chi2: Chi square　(ML)'/format f7.4.
    compute d0=(S-C)*inv(C).
    compute Fk0=trace(d0*d0)/2.
    print Fk0/title 'Fk0: F 式(1)'/format f7.4.
    print (Fk0*(n-1))/title "(1)式から求めたχ2値"/format f7.4.
    end matrix.

    matrix data var=v1 to v4/format=full/N=240/contents=corr.
    begin
    1.0, .3, .2, .1
    .3,1., .2, .2
    .2, .2,1., .3
    .1, .2, .3,1
    end data.
    '*---------------------------マクロ終わり

    実行結果
    ---------------------------------
    Run MATRIX procedure:

    標本相関行列 S
    　1.000　　.300　　.200　　.100
    　　.300 　1.000　　.200　　.200
    　　.200　　.200 　1.000　　.300
    　　.100　　.200　　.300 　1.000

    M: 変数の数
    　4

    N サンプルサイズ
    　240

    １因子解のときの重み係数 A
    　　.430
    　　.502
    　　.502
    　　.430

    C: 重みから再生相関行列を作る。
    　1.000　　.216　　.216　　.185
    　　.216 　1.000　　.252　　.216
    　　.216　　.252 　1.000　　.216
    　　.185　　.216　　.216 　1.000

    Fk 適合度関数：本来の方法(6)
    　.0390

    Chi2: Chi square　(ML)
    　9.3100

    Fk0: F 式(1)
    　.0397

    (1)式から求めたχ2値
    　9.4953

    ------ END MATRIX -----
    ---------------------------------

    *----------------------------spss で求めてみる.
    matrix data var=v1 to v4/format=full/N=240/contents=corr.
    begin
    1.0, .3, .2, .1
    .3,1., .2, .2
    .2, .2,1., .3
    .1, .2, .3,1
    end data.
    FACTOR
    /matrix=in(cor=*)
    /PRINT UNIVARIATE INITIAL DET INV REPR EXTRACTION
    /CRITERIA FACTORS(1) ITERATE(25)
    /EXTRACTION ML
    /ROTATION NOROTATE .
    *---------------.

    おや，χ2値が違ってますね。
    SPSSの結果
    適合度検定
    カイ2乗自由度 有意確率
    9.199 2 .010
    SPSSは上の方法でもとめていませんね。Bartlettの修正をしている。
    修正を戻すには，
    compute chi2=chi2B/(n-1-(2*nvar+5)/6-2*nfac/3)*(n-1).
    とするとほぼ一致した。
    9.199/(240-1-(2*4+5)/6-2*1/3)*(240-1)=9.309362032


    ＡＭＯＳは行列言語と一致しました(χ2=9.310)。
    -------------------------------
    AMOS出力(カラムずれご容赦)

    カイ2乗 = 9.310
    自由度 = 2
    確率水準 = 0.010


    最尤(ML)推定値
    --------------

    係数：　　　　　　　　　　　　　　　　　　推定値　　標準誤差 検定統計量 ラベル 　
    ------ 　　　　　　　　　　　　　　　　　-------- 　------- 　------- 　-------

    　　　　　　　　　　　v1 <--- 1因子 　　　0.430 　　0.091 　　4.726 　　　　　
    　　　　　　　　　　　v2 <--- 1因子 　　　0.502 　　0.094 　　5.326 　　　　　
    　　　　　　　　　　　v3 <--- 1因子 　　　0.502 　　0.094 　　5.326 　　　　　
    　　　　　　　　　　　v4 <--- 1因子 　　　0.430 　　0.091 　　4.726 　　　　　
    --------------------


    モデルの概要
    ----------

    　　　　　　　モデル　NPAR　　　　CMIN　　DF 　　　　　P 　　CMIN/DF
    　　----------------　---- 　---------　　-- 　--------- 　---------
    　　　Default model 　　8 　　　9.310 　　2 　　　0.010 　　　4.655
    　　　　　飽和モデル　　10 　　　0.000 　　0
    　　　　　独立モデル 　　4　　　65.685 　　6 　　　0.000　　　10.947

    　　　　　　　モデル 　　　　RMR 　　　　GFI　　　　AGFI　　　　PGFI
    　　----------------　----------　----------　----------　----------
    　　　Default model 　　　0.049 　　　0.980 　　　0.902 　　　0.196
    　　　　　飽和モデル 　　　0.000 　　　1.000　　　　　　　　　　　　
    　　　　　独立モデル 　　　0.175 　　　0.866 　　　0.776 　　　0.519

    　　　　　　　　　　　　　DELTA1　　　　RHO1　　　DELTA2　　　　RHO2
    　　　　　　　モデル 　　　　NFI 　　　　RFI 　　　　IFI 　　　　TLI 　　　　CFI
    　　----------------　----------　----------　----------　----------　----------
    　　　Default model 　　　0.858 　　　0.575 　　　0.885 　　　0.633 　　　0.878
    　　　　　飽和モデル 　　　1.000 　　　　　　　　　1.000 　　　　　　　　　1.000
    　　　　　独立モデル 　　　0.000 　　　0.000 　　　0.000 　　　0.000 　　　0.000

    　　　　　　　モデル　　　PRATIO　　　　PNFI　　　　PCFI
    　　----------------　----------　----------　----------
    　　　Default model 　　　0.333 　　　0.286 　　　0.293
    　　　　　飽和モデル 　　　0.000 　　　0.000 　　　0.000
    　　　　　独立モデル 　　　1.000 　　　0.000 　　　0.000

    　　　　　　　モデル 　　　　NCP 　　　LO 90 　　　HI 90　　　　　　
    　　----------------　----------　----------　----------
    　　　Default model 　　　7.310 　　　1.289　　　20.804
    　　　　　飽和モデル 　　　0.000 　　　0.000 　　　0.000
    　　　　　独立モデル　　　59.685　　　37.265　　　89.560

    　　　　　　　モデル　　　　FMIN　　　　　F0 　　　LO 90 　　　HI 90
    　　----------------　----------　----------　----------　----------
    　　　Default model 　　　0.039 　　　0.031 　　　0.005 　　　0.087
    　　　　　飽和モデル 　　　0.000 　　　0.000 　　　0.000 　　　0.000
    　　　　　独立モデル 　　　0.275 　　　0.250 　　　0.156 　　　0.375

    　　　　　　　モデル 　　　RMSEA 　　　LO 90 　　　HI 90　　　PCLOSE
    　　----------------　----------　----------　----------　----------
    　　　Default model 　　　0.124 　　　0.052 　　　0.209 　　　0.046
    　　　　　独立モデル 　　　0.204 　　　0.161 　　　0.250 　　　0.000

    　　　　　　　モデル 　　　　AIC 　　　　BCC 　　　　BIC　　　　CAIC
    　　----------------　----------　----------　----------　----------
    　　　Default model　　　25.310　　　25.651　　　64.245　　　61.155
    　　　　　飽和モデル　　　20.000　　　20.427　　　68.669　　　64.806
    　　　　　独立モデル　　　73.685　　　73.856　　　93.152　　　91.607

    　　　　　　　モデル　　　　ECVI 　　　LO 90 　　　HI 90 　　　MECVI
    　　----------------　----------　----------　----------　----------
    　　　Default model 　　　0.106 　　　0.081 　　　0.162 　　　0.107
    　　　　　飽和モデル 　　　0.084 　　　0.084 　　　0.084 　　　0.085
    　　　　　独立モデル 　　　0.308 　　　0.214 　　　0.433 　　　0.309

    　　　　　　　　　　　　HOELTER 　　HOELTER
    　　　　　　　モデル 　　　　.05 　　　　.01
    　　----------------　----------　----------
    　　　Default model 　　　　154 　　　　237
    　　　　　独立モデル　　　　　46　　　　　62
    --------------------------------------------
    EQS でも行列言語の結果(χ2=9.310)と一致しました。
    MPLUS は9.349
    9.349/240*239＝9.310045833
    となるので，Fk*N でχ2値を求めているようだ。


    メニューへ　　トップへ　(123)へ

    SPSS ときど記(121)　2002/11/28

    SPSS　SPSS 11.5 日本語版のリリース

    SPSS11.5Jのリリースが    http://www.spss.co.jp/でアナウンスされた。2002年12月25日発売だそうだ。

    http://www.spss.co.jp/product/ALL/base/11_5.htm

    に新しい機能を紹介している。

    メニューへ　　トップへ　(122)へ
    堀　啓造ホームページへ