猫も杓子も構造化

発達障害、特別支援などについて書いています。最近は心理学関係の内容が多めです。

いきなり因子分析(その6):lavaanを使って確証的因子分析

シリーズものの記事。過去の記事はこちら。

その1(とりあえず試した編)
その2(因子数の決定編)
その3(様々な推定法編)
その4(因子軸の回転編)
その5(高次因子モデルと双因子モデル編)

前回までは基本的に探索的因子分析(Exploratory Factor Analysis)に焦点を当てていた。これは観測されたデータを説明するモデルを探索的に調べるときに使われる手法である。

一方で、観測データと因子との間の関係についての仮説的なモデルを立てて、そのモデルをデータに基づき評価する際に使われるのが、確証的因子分析(Confirmatory Factor Analysis)である。(あるいは確認的因子分析や検証的因子分析とも)

前回までのデータセットを用いてやってみる。

lavaanのパッケージ

Rで確認的因子分析をするときはlavaanというパッケージに含まれるcfa()を使えばできるようだ。さっそくやってみる。

まず最初にやるべきことは、モデルを表現する式を書くことらしい。式の記法は決まっており, lavaanパッケージの説明pdfの日本語版には次のように書いてある。

モデルのタイプ オペレータ 意味
潜在変数の定義 =∼ ~により測定される
回帰 ~への回帰
(残差) (共)分散 ∼∼  〜と相関を持つ
切片 ∼1 切片

これらを使って、潜在変数・観測変数の関係性を記述するのである。

前回までの記事では、観測変数が9つあり、それらが「視空間能力」「言語」「処理速度」と3つの因子によって説明されるとしていたのであった。 そこで、潜在変数を定義するオペレータを用いて次のようにモデルを書く。

library(lavaan)
dat <- HolzingerSwineford1939

model1 <- "
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
"

モデル式の右辺はデータセットの列名である。モデルは''もしくは""で囲む。

あとはcfa()に作ったモデルとデータを読ませるだけである。

# モデルへのあてはめ
fit1 <- cfa(model = model1,
            data = dat
            )

modelという引数には書いたモデルを、dataにはデータセットを指定する。

summary(fit1)を用いると結果の要約をみることができるのでさっそく見てみよう。

> summary(fit1)
lavaan (0.6-1) converged normally after  35 iterations

  Number of observations                           301

  Estimator                                         ML
  Model Fit Test Statistic                      85.306
  Degrees of freedom                                24
  P-value (Chi-square)                           0.000

一番最初には

  • lavaanのバージョン
  • 解が収束したか
  • 何回の繰り返しが必要だったかが
  • 分析で使用された観察数
  • 使われた推定量
  • 検定統計量、自由度、p値

などが表示されている。

次のブロックにはパラメータについての記述。得られたパラメータが期待値なのか観測値なのか、標準誤差が標準のものかブートストラップなのかなどが表示されるらしい。(詳しくはよく分からん)

Parameter Estimates:

  Information                                 Expected
  Information saturated (h1) model          Structured
  Standard Errors                             Standard

そこから先が各変数のパラメータの推定値とその標準誤差、Waldの統計量とそのp値が出力される。Waldの統計量とは推定値をその標準誤差で割ったもののようだ。

cfa()は各潜在変数の最初の因子負荷量を1に固定しているので、当然のことながら推定をおこなっていないそれらの標準誤差は存在しない。

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  visual =~                                           
    x1                1.000                           
    x2                0.554    0.100    5.554    0.000
    x3                0.729    0.109    6.685    0.000
  textual =~                                          
    x4                1.000                           
    x5                1.113    0.065   17.014    0.000
    x6                0.926    0.055   16.703    0.000
  speed =~                                            
    x7                1.000                           
    x8                1.180    0.165    7.152    0.000
    x9                1.082    0.151    7.155    0.000

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)
  visual ~~                                           
    textual           0.408    0.074    5.552    0.000
    speed             0.262    0.056    4.660    0.000
  textual ~~                                          
    speed             0.173    0.049    3.518    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .x1                0.549    0.114    4.833    0.000
   .x2                1.134    0.102   11.146    0.000
   .x3                0.844    0.091    9.317    0.000
   .x4                0.371    0.048    7.779    0.000
   .x5                0.446    0.058    7.642    0.000
   .x6                0.356    0.043    8.277    0.000
   .x7                0.799    0.081    9.823    0.000
   .x8                0.488    0.074    6.573    0.000
   .x9                0.566    0.071    8.003    0.000
    visual            0.809    0.145    5.564    0.000
    textual           0.979    0.112    8.737    0.000
    speed             0.384    0.086    4.451    0.000

ところで、モデルの当てはまりの良さを表す指標を得たい場合には次のように引数を指定する。

summary(fit1, fit.measures=TRUE)

そうすると、出力が表示される次のようなブロックが表示される。

Model test baseline model:

  Minimum Function Test Statistic              918.852
  Degrees of freedom                                36
  P-value                                        0.000

User model versus baseline model:

  Comparative Fit Index (CFI)                    0.931
  Tucker-Lewis Index (TLI)                       0.896

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)              -3737.745
  Loglikelihood unrestricted model (H1)      -3695.092

  Number of free parameters                         21
  Akaike (AIC)                                7517.490
  Bayesian (BIC)                              7595.339
  Sample-size adjusted Bayesian (BIC)         7528.739

Root Mean Square Error of Approximation:

  RMSEA                                          0.092
  90 Percent Confidence Interval          0.071  0.114
  P-value RMSEA <= 0.05                          0.001

Standardized Root Mean Square Residual:

  SRMR                                           0.065

一番最初にbeseline modelというものについての情報が表示される。これは独立モデル(independence model)とも呼ばれるもので、変数と変数の間の関連を一切想定しないモデルのことを指す。当然そんなモデルは当てはまりがとっても悪いので、そことの比較においてモデルの当てはまりを評価しようというのが基本の考え方だ。

次のブロックには、指定したモデルとベースラインモデルを比較した指標が2つ出力される。どちらも、指定されたモデルとベースラインモデルのχ二乗値を用いて計算する。CFIもTLIも1に近いほど、当てはまりがよいとされる。

その次のブロックは対数尤度と情報量基準に基づく指標。対数尤度尤度の方はちょっとよく知らないけれど、情報量基準の方は2つ以上のモデルを比較するときに使われるようである。小さいほど当てはまりが良いとされ絶対値に意味はない。

さらに次のブロックにはRMSEAとSRMRは0.1以上だと当てはまりがよくないとされるとか。

χ二乗値が計算のもとになるのだけれど、これが自由度とかサンプルサイズとかの影響を受けるから、それをうけないように指標を改良してきているようである。

モデルの比較

せっかくなので、モデルを作って比較してみる。以前の記事で探索的因子分析の際に、そういえば観測変数x9が因子2と因子3のどちらにも影響されていてどっちつかずだったので、x9の観測変数を削除したモデルと比較してみることにしよう。モデルを次のように書いて、同様に結果を出力する。

model2 <- "
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8
"

fit2 <- cfa(model = model2,
            data = dat
            )
summary(fit2, fit.measure=TRUE)

これを、最初に行ったモデル1の適合度の指標を比べてみる。

指標 モデル1 モデル2
chi2 85.306 44.409
Comparative Fit Index (CFI) 0.931 0.964
Tucker-Lewis Index(TLI) 0.896 0.941
AIC 7517.490 6742.935
Bayesian (BIC) 7595.339 6813.370
Sample-size adjusted Bayesian (BIC) 7528.739 6753.113
RMSEA 0.092 0.073
SRMR 0.065 0.052

モデルのデータへの当てはまりのみを見れば項目9を削除した方が良いことがわかる。もっとも、モデルの選択であったり項目の削除であったりは、適合度のみをもとに判断するのでなく、理論的な仮説に基づいてなされるべきであるので、モデルが想定する変数間の関係を変えて適合度を上げつつも、理論的な仮説と整合性のあるモデルを組まなければならないようだ。(だから今回のように内容を吟味することもなく数値だけで判断するのは過ちということである)

いきなり因子分析(その5):高次因子分析モデルと双因子モデル)

因子分析についてのシリーズものの記事。

その1(とりあえず試した編)
その2(因子数の決定編)
その3(様々な推定法編)
その4(因子軸の回転編)

今回はもう少し複雑なモデルということで高次因子分析モデルと階層因子分析モデルのうち双因子モデルを試してみる。

高次因子分析って?

高次因子分析というのは、複数の因子をさらに少数の因子で説明するモデルである。たとえば、因子分析の結果因子間に強い相関が仮定される場合に、それらに影響を与える上位の因子を想定するモデルである。

こうした場合に、観測を変数を説明する因子(通常の因子分析で想定する因子)を1次因子と言い、それらの因子を説明する上位の因子を2次因子というそうだ。

パス図で表すとこんな感じ。

f:id:nekomosyakushimo:20180714232805p:plain

双因子モデルって?

双因子モデルというのは、階層因子分析の一部とされる。階層因子分析というのは、一般因子とグループ因子で階層を分けて観測変数を説明するモデルである。一般因子は、全ての観測変数に影響を与えている因子で、グループ因子は通常の因子分析と同じように、一部の観測辺境に影響を与えるモデルらしい。

パス図で表すとこんな感じ。

f:id:nekomosyakushimo:20180714232832p:plain

Rで高次因子分析を実践

psychを使って、高次因子分析モデルや双因子モデルで因子分析する方法は次のpdfを参考にした。 http://personality-project.org/r/psych/HowTo/factor.pdf

データは前の記事と変わらずHolzingerSwineford1939を使うことにする。

どちらのモデルもpsychに含まれるomega()という関数を使ってできるらしい。

ちなみに日本語でこの関数についてググると大体がω係数の算出に関連したものがヒットする。ω係数とは、信頼性の指標であるαと同じように、内的一貫性を表す指標だとされている。因子負荷と誤差分散を用いて算出される値らしく、この関数を使うとそれが算出されるのでそうした記事がヒットするようだ。

それはさておき、さっそく使ってみる。

library(psych)
model_higherorder <- omega(m=dat[7:15],
                           nfactors = 3,
                           sl=FALSE
                           )

引数のm=でデータの指定。ローデータでも相関行列でもどちらでも良いっぽい。

nfactorsには因子数を指定する。省略すると、何だかよく分からないメカニズムで因子数を推定して結果を返すけれどもきっと指定した方が良いでしょう。(因子数の決め方についてはその2 を参照)

slというのは、Schmid Leiman変換というもので、これをTRUEにすると双因子モデルになり、FALSEだと高次因子分析モデルになる。

他にもfmで、前記事で扱った推定法を選べたりもする。(デフォルトは"minres"なようだ)

さて、結果を見てみる。summary(model_higherorder)を打つと結果が表示される。

> summary(model_higherorder)
Omega
Alpha:                 0.76
G.6:                   0.81
Omega Hierarchical:    0.45
Omega H asymptotic:    0.53
Omega Total            0.85

まず最初のブロックはω係数とか信頼係数の値が表示されている。Aplhaはクロンバックのα係数、G.6はガットマンの提案した信頼性係数の6つ目の下界である(詳しく知りたい人は岡田先生のこの論文を読むと良い)。

そこから3つがω係数について。調べた感じだとOmega Hierarchicalは一般因子の因子負荷係数を基に算出したもので、Omega Totalの方は一般因子だけでなく、グループ因子も足して計算したものらしい。omega H asymptoticというのは、計算式をみる限りだと全ての因子負荷に占める一般因子部分の割合みたいに読めるのだけれどちょっとどう解釈して良いかは分からない。信頼性係数として報告する際に使うのはOmega Totalの値である。

With eigenvalues of:
   g  F1*  F2*  F3*
1.44 1.62 0.77 1.03

次には一般因子とグループ因子の固有値が出てくる。まぁこれはそのまま固有値ですね。

The degrees of freedom for the model is 12  and the fit was  0.08
The number of observations was  301  with Chi Square =  22.56  with prob <  0.03

The root mean square of the residuals is  0.02
The df corrected root mean square of the residuals is  0.05

RMSEA and the  0.9 confidence intervals are  0.055 0.016 0.088
BIC =  -45.93Explained Common Variance of the general factor =  0.3

そこから先は、モデルの評価に関するあれやこれやの指標。普通に因子分析をしたときと同様のものが出てくる。一つ違うのは、一般因子のECVというものが出力されていることでしょうか。これはomega()のhelpに次のように書いてあった。

Following suggestions by Steve Reise, the Explained Common Variance (ECV) is also reported. This is the ratio of the general factor eigen value to the sum of all of the eigen values. As such, it is a better indicator of unidimensionality than of the amount of test variance accounted for by a general factor.

Total, General and Subset omega for each subset
                                                g  F1*  F2*  F3*
Omega total for total scores and subscales    0.85 0.89 0.65 0.72
Omega general for total scores and subscales  0.45 0.24 0.27 0.19
Omega group for total scores and subscales    0.35 0.65 0.37 0.53

最後には、それぞれの因子が合計得点および下位尺度のばらつきをどの程度説明されるかが表示される。1行目はOmega totalとは、一般因子Omega generalと(2行目)とグループ因子Omega group(3行目)を足したものである。一番左の列は、結果の最初に表示されたω係数と同じである。右の3列のサブセットというのがそれぞれの一次因子が負荷していると想定される観測変数のω係数ということだろうか(よく分からん)。

summaryだと、これしか出ないのでが、高次の因子から下位の因子への負荷ぐらい知っておきたいだろう。そういう場合は、次のようにするとその値を取り出せる。

print(model_higherorder$schmid$gloading, digits=2)

Loadings:
   MR1  
F1 0.514
F2 0.628
F3 0.415

                 MR1
SS loadings    0.831
Proportion Var 0.277

双因子モデルを実践

今度は、双因子モデルを試してみる。といっても、引数sl=TRUEと指定するのみだけれど。

model_higherorder <- omega(m=dat[7:15],
                           nfactors = 3,
                           sl=TRUE
)

summary(model_bifactor)

結果は高次因子分析モデルと同様に出てくる。ただ、一般因子から観測因子への負荷が知りたいこともあるでしょう。その場合には、print(model_bifactor$schmid$sl, digit = 2)とすると値が取り出せる。

> print(model_bifactor$schmid$sl, digit = 2)
      g    F1*    F2*     F3*   h2   u2   p2
x1 0.49  0.169  0.460  0.0282 0.48 0.52 0.49
x2 0.29  0.037  0.396 -0.1113 0.26 0.74 0.33
x3 0.41 -0.053  0.534  0.0173 0.45 0.55 0.37
x4 0.45  0.726  0.013  0.0070 0.73 0.27 0.28
x5 0.42  0.760 -0.050  0.0061 0.75 0.25 0.23
x6 0.46  0.691  0.062 -0.0118 0.69 0.31 0.30
x7 0.23  0.038 -0.119  0.6701 0.52 0.48 0.10
x8 0.35 -0.029  0.097  0.6246 0.52 0.48 0.23
x9 0.45  0.027  0.297  0.4152 0.46 0.54 0.43

普通に因子分析をしたときと同じように見ると良いとが一番右にp2という列がある。Helpによれば、一般因子が共通性に占める割愛のようで、これが高いほど観測変数は一般因子で説明されるということみたいである。

おわり。

いきなり因子分析(その4):因子軸の様々な回転

まだまだ因子分析の続き。

その1(とりあえず試した編)
その2(因子数の決定編)
その3(様々な推定法編)

今回は因子の回転について。

様々な因子軸の回転法

このシリーズものの記事では、Rのpsychというパッケージのfa()という関数を使っているのだが、これは 実にたくさんの回転方法を指定できる。

関数のHelpを見てみると因子間の相関を仮定しない直交回転だけで次のものがある。

  • 初期解 "none"
  • リマックス "varimax"
  • クォーティマックス "quartimax"
  • ベントラー(直交) "bentlerT"
  • エカマックス "equamax"
  • バリミン? "varimin"
  • ジオミン(直交)"geominT"
  • バイファクター"bifactor"

因子間に相関を仮定する斜交回転では次のものがある。

  • プロマックス "Promax"
  • プロマックス"promax"
  • オブリミン "oblimin"
  • シンプリマックス"simplimax"
  • ベントラー(斜交)"bentlerQ
  • ジオミン(斜交)"geominQ"
  • バイクォーティミン"biquartimin"
  • 独立クラスタ"cluster"

日本語で調べても出てこないものについては、カタカタでこれがあっているのかも分からない。とにかくたくさん回転法があることは分かった。

ちなみにプロマックスが2つあるのは、SPSSのプロマックスは回転の前になんか処理を施しているので、その処理を施した版とそうでない版を用意しているとのことだ。Helpには以下のように書いてある。

SPSS seems to do a Kaiser normalization before doing Promax, this is done here by the call to "promax" which does the normalization before calling Promax in GPArotation.

回転方法の決め方

これだけあるとどれを選べばよいか分からないとの声が聞こえてきそうだ。少なくとも私はどれを選べばよいか分からない。

そういう困ったことがあるときは、清水先生のHPをみると良い。

因子分析における因子軸の回転法についてなんていう素敵な記事が見つかる。回転法についてそれぞれがどんな回転なのかを説明してくれている激アツ資料である。

そこには、次のような回転法の決め方の提案がなされている。

  1. 基本はプロマックス回転で問題ない。プロマックス回転で思うような解が得られない場合は他の回転を試す。
  2. より単純構造を目指したいなら独立クラスター回転がオススメ。
  3. そもそも単純構造にならないと思われるデータの場合、ジェオミンやオブリミン回転がオススメ。
  4. 先行研究を再現したい場合は、出てほしい因子負荷量を指定してプロクラステス回転を使う。
  5. 他にも確認的因子分析という選択肢もある。

Rで実験

早速Rで試してみる。プロマックスと単純構造度合いが強いとされる独立クラスタと複雑な構造を許容するというジオミンを試してみることにする。

データは前記事と同じようにlavaanHolzingerSwineford1939を使うことにする。

result_promax <- fa(r = dat[7:15],
                    nfactor = 3,
                    fm = "ml",
                    rotate= "promax",
                    use = "pairwise"
                    )
result_cluster <- fa(r = dat[7:15],
                    nfactor = 3,
                    fm = "ml",
                    rotate= "cluster",
                    use = "pairwise"
                    )
result_geominQ <- fa(r = dat[7:15],
                     nfactor = 3,
                     fm = "ml",
                     rotate= "geominQ",
                     use = "pairwise"
                     )

結果のうち因子負荷量の部分を示すとそれぞれ次の通り。

プロマックス

ML1 ML2 ML3 h2 u2 com
x1 0.153 0.036 0.609 0.487 0.513 1.13
x2 0.013 -0.116 0.525 0.251 0.749 1.1
x3 -0.115 0.029 0.703 0.457 0.543 1.06
x4 0.844 0.005 0.01 0.721 0.279 1
x5 0.898 0.007 -0.082 0.757 0.243 1.02
x6 0.807 -0.011 0.068 0.695 0.305 1.01
x7 0.044 0.743 -0.215 0.498 0.502 1.17
x8 -0.049 0.722 0.049 0.531 0.469 1.02
x9 0.005 0.479 0.335 0.457 0.543 1.79

独立クラスタ

ML1 ML2 ML3 h2 u2 com
x1 0.151 0.083 0.595 0.487 0.513 1.17
x2 0.012 -0.078 0.516 0.251 0.749 1.1
x3 -0.115 0.083 0.682 0.457 0.543 1.09
x4 0.839 0.005 0.02 0.721 0.279 1
x5 0.892 0.000 -0.065 0.757 0.243 1.01
x6 0.802 -0.007 0.081 0.695 0.305 1.02
x7 0.043 0.735 -0.238 0.498 0.502 1.21
x8 -0.049 0.735 0.019 0.531 0.469 1.01
x9 0.004 0.510 0.308 0.457 0.543 1.64

ジオミン

ML1 ML3 ML2 h2 u2 com
x1 0.188 0.604 0.029 0.487 0.513 1.20
x2 0.044 0.507 -0.119 0.251 0.749 1.1
x3 -0.073 0.691 0.020 0.457 0.543 1.02
x4 0.839 0.024 0.01 0.721 0.279 1
x5 0.887 -0.065 0.010 0.757 0.243 1.01
x6 0.806 0.080 -0.009 0.695 0.305 1.02
x7 0.031 -0.150 0.726 0.498 0.502 1.09
x8 -0.045 0.106 0.703 0.531 0.469 1.05
x9 0.025 0.368 0.463 0.457 0.543 1.91

結果を比べてみると独立クラスター回転では、観測変数x9のメリハリをよりつけるようにしたということで、言われているように単純構造を目指している感じがする。一方でジオミンはやや、もやっとした構造になったような印象を受ける。

ちなみに複雑性の平均を比較すると次の通り(複雑性は1に近いほど単純構造ということであった)。

回転法 複雑性の平均
プロマックス 1.1
独立クラスタ 1.1
ジオミン 1.2

やはりジオミンはやや複雑な構造を持つと想定される場合に使うと良いのかもしれない。

通常学級の特別支援本

通常学級の特別支援─今日からできる! 40の提案─

通常学級の特別支援─今日からできる! 40の提案─

通常学級の特別支援セカンドステージ―6つの提言と実践のアイデア50

通常学級の特別支援セカンドステージ―6つの提言と実践のアイデア50

最近読んだ実践的な本の紹介でも。通常学級の特別支援に焦点を当てた本。著者の現場での経験がコンパクトにまとまっており、通常学級の特別支援を考える人には参考になる点も多いだろう。

今だったらきっと「ユニバーサルデザイン」って書名のどこかについて売られるだろうなぁ、というような実践でのアイディアや考え方が多い。1冊目の方は2008年が初版なので特別支援教室が始まって1年ぐらいしかたっていない時期だが、古さは感じないし、UD界隈で言われるような考え方の多くが本の中でも紹介されている。

本全体の印象として、当事者や保護者へ敬意というか礼儀というかがあって、そこには大変好感が持てた。たとえば、1冊目の提案1は「子どもの本音の思いを大切に」というタイトルがついていて、当事者の声から始まっている。本では様々な提案がなされるのだが、出発点が「子どもの見え方・聞こえ方・感じ方」というのが、この本全体の持つ性格を象徴しているように思う。

2冊目の方に、家庭との連携の際のプライバシーへの配慮や親の生活への影響についてかかれており、次のような記述がある。

少なくとも次のことは言えるーー子供の家庭・地域生活づくりは親の生活の一部の変更を確実に迫る。子どもの生活づくりは親の生活づくりでもある。このことを心に留めて支援したい。(p.156)

忘れがちなことだけど、大変大切なことに思う。

本の内容には関係ないが、読んでいるうちに著者の研修を昔1度受けたことがあるのを思い出した。小ネタの多い研修で現場の先生方には大ウケだったと記憶している。(私はアクティブな研修があまり好きでないので微妙な気持ちで参加していたのだけれど・・・)

以下は自分のための雑多なメモ。

1冊目

-「 "できる""やれる"支援条件を見つける(p.23)」というアイディア(ワーディング)。

見方を変えれば、この子どもたちは、周りの環境要因を含む支援条件に大変敏感なのだ。私たちの日常を考えても、いつも"できる"こともー聞き慣れない大音量の音楽が流れている部屋ではー集中できず"できない"ことだってある。中略。子ども理解とは、子どもの周りの支援条件を理解することでもある。そして、支援条件に目を向け、改善の努力を図ることが子どもの味方になる第一歩と考えたい。(p.23-24)

  • 学年会を中心とした仕組みづくりで学級を支える話(p.33)
  • 授業は子どもたちの学校生活の中心、特別支援教育の根幹に位置づくべき(p.77)
  • 授業のユニット化の話(p.86)
  • 丸つけをされるのが嬉しくなるキャラクターのアイディア(p.88)

2冊目

  • 「ないと困る支援の専門性:学級経営、授業づくりの専門性(p.iii)」
  • 特別支援教育ユニバーサルデザイン)が「超教科的・超領域的(p.40)」という話は今までの特別支援の発展やSENSのセミナーの内容を見るとよく分かる。もっとも授業UDとかは超教科的な基盤の上に、教科としての専門性を成り立たせようという発想も感じるけど(あくまで外からみた感想)
  • 「〇〇したい。」で終わる分がやたら多い。なんかロマンチック。
  • 「不適応という適応の現実(p.147)」ここらへんの記述は特別支援の実践家であると同時に障害児の親という立場である筆者だからこそリアリティのある書き方ができるのかもしれない。
  • 家庭生活に関連することの聞き取りについてのプライバシーのくだり(p.153) とかこの感覚大事。

いきなり因子分析(その3):様々な推定法を試してみた

以前からの続きで因子分析を試してみるシリーズ。

その1

その2

今回はいろいろな推定方法について。あいもかわらず清水先生のスライドを頼りにしながら試してみる。

Rで因子分析 商用ソフトで実行できない因子分析のあれこれ

データセットは前回、前々回と同じでlavaanに入っているHolzingerSwineford1939を利用する。

psychパッケージ

一番最初の記事で、デフォルトで入っているfactanal()という関数を用いたが、この関数は最尤法しか使えないらしい。違った推定方法ができるものとしてpsychfa()という関数があるようだ。

あと、因子分析を実行する際にGPArotationというパッケージも必要になるようなのでその2つをインストールして実行してみる。

library(psych)
library(GPArotation)

result <- fa(dat[7:15],
             nfactor = 3,
             fm = "minres",
             rotate= "promax",
             use = "pairwise"
             )

rという引数にデータを指定する。今回のようにローデータを突っ込んでもよいし、相関行列を突っ込んでもよいらしい。

nfactorは因子数である。前回記事で3つがおすすめという結論が出たのでそれを使う。

fmの引数に推定法を指定する。最小残差法minres、最尤法ml、反復主因子法pa、一般化最小二乗法gls、重みつき最小二乗法wlsといったように色々指定できる。(知らないものが多い)

useの引数は欠損値の処理でデフォルトはペアワイズ。Helpをみると、cor()と同じ引数が指定できるよと書いてある。

結果を見てみよう。

print(result, digits = 3)

Factor Analysis using method =  minres
Call: fa(r = dat[7:15], nfactors = 3, rotate = "promax", fm = "minres",
    use = "pairwise")
Standardized loadings (pattern matrix) based upon correlation matrix
      MR1    MR2    MR3    h2    u2  com
x1  0.157  0.039  0.598 0.477 0.523 1.15
x2  0.010 -0.120  0.531 0.255 0.745 1.10
x3 -0.109  0.028  0.699 0.453 0.547 1.05
x4  0.849  0.007  0.005 0.728 0.272 1.00
x5  0.895  0.005 -0.078 0.754 0.246 1.02
x6  0.804 -0.013  0.072 0.691 0.309 1.02
x7  0.047  0.759 -0.229 0.519 0.481 1.19
x8 -0.050  0.710  0.060 0.520 0.480 1.02
x9  0.001  0.476  0.344 0.460 0.540 1.82

                        MR1   MR2   MR3
SS loadings           2.211 1.328 1.318
Proportion Var        0.246 0.148 0.146
Cumulative Var        0.246 0.393 0.540
Proportion Explained  0.455 0.273 0.271
Cumulative Proportion 0.455 0.729 1.000

 With factor correlations of
      MR1   MR2   MR3
MR1 1.000 0.257 0.391
MR2 0.257 1.000 0.350
MR3 0.391 0.350 1.000

Mean item complexity =  1.2
Test of the hypothesis that 3 factors are sufficient.

The degrees of freedom for the null model are  36  and the objective function was  3.053 with Chi Square of  904.097
The degrees of freedom for the model are 12  and the objective function was  0.077

The root mean square of the residuals (RMSR) is  0.019
The df corrected root mean square of the residuals is  0.033

The harmonic number of observations is  301 with the empirical chi square  7.867  with prob <  0.795
The total number of observations was  301  with Likelihood Chi Square =  22.555  with prob <  0.0317

Tucker Lewis Index of factoring reliability =  0.9633
RMSEA index =  0.0553  and the 90 % confidence intervals are  0.0157 0.0882
BIC =  -45.93
Fit based upon off diagonal values = 0.996
Measures of factor score adequacy             
                                                    MR1   MR2   MR3
Correlation of (regression) scores with factors   0.944 0.859 0.848
Multiple R square of scores with factors          0.892 0.738 0.719
Minimum correlation of possible factor scores     0.783 0.477 0.439

上のブロックには因子負荷量が表示されている。h2の列が共通性でu2が独自性。デフォルトの関数と違ってどっちも表示されるので便利。comというのは複雑性というもので1に近いほどその変数は単純構造に近いらしい。上の観測変数だとx9は因子2と因子3でどっちつかずなため複雑性が1から離れている。

その下のブロックには因子寄与や因子寄与率、累積因子寄与などが表示されている。説明率と累積説明率もその下に表示される。

次のブロックは因子間相関が表示される。ちなみに直行回転だとこれは表示されない。

そのあとは複雑性の平均が表示されている。きっと観測変数のx9を削除するともう少しあがるのかな?

下の方には、様々な適合度の指標が順番に表示される。χ2乗値、RMSR、Tucker Lewis Index、、RMSEA、BICなど。

RMSEAだと0.050以下が良好だとか、TuckerLewisIndexは1に近いほどよいだとか指標ごとにあてはまりの程度を教えてくれる。BICは同一データからのモデル間の比較などに使われて小さいほど良いとか言われているもの。(どうやって算出しているかは調べたことがないから知らない)

違った推定法を比較してみる

fa()の用法の説明を終えたところで、推定法だけを変えてみて結果がどう変わるのかを見てみる。

result_minres <- fa(r = dat[7:15],
                    nfactor = 3,
                    fm = "minres",
                    rotate= "promax",
                    use = "pairwise"
                    )
result_ml <- fa(r = dat[7:15],
                nfactor = 3,
                fm = "ml",
                rotate= "promax",
                use = "pairwise"
                )
result_gls <- fa(r = dat[7:15],
                 nfactor = 3,
                 fm = "gls",
                 rotate= "promax",
                 use = "pairwise"
                 )
print(result_minres, digits = 3)
print(result_ml, digits = 3)
print(result_gls, digits = 3)
print(result_wls, digits = 3)
print(result_pa, digits = 3)

最小残差法

MR1 MR2 MR3 h2 u2 com
x1 0.157 0.039 0.598 0.477 0.523 1.15
x2 0.010 -0.120 0.531 0.255 0.745 1.1
x3 -0.109 0.028 0.699 0.453 0.547 1.05
x4 0.849 0.007 0.005 0.728 0.272 1
x5 0.895 0.005 -0.078 0.754 0.246 1.02
x6 0.804 -0.013 0.072 0.691 0.309 1.02
x7 0.047 0.759 -0.229 0.519 0.481 1.19
x8 -0.050 0.710 0.060 0.520 0.480 1.02
x9 0 0.476 0.344 0.460 0.540 1.82

最尤法

ML1 ML2 ML3 h2 u2 com
x1 0.153 0.036 0.609 0.487 0.513 1.13
x2 0.013 -0.116 0.525 0.251 0.749 1.1
x3 -0.115 0.029 0.703 0.457 0.543 1.06
x4 0.844 0.005 0.010 0.721 0.279 1
x5 0.898 0.007 -0.082 0.757 0.243 1.02
x6 0.807 -0.011 0.068 0.695 0.305 1.01
x7 0.044 0.743 -0.215 0.498 0.502 1.17
x8 -0.049 0.722 0.049 0.531 0.469 1.02
x9 0 0.479 0.335 0.457 0.543 1.79

一般化最小二乗法

GLS1 GLS2 GLS3 h2 u2 com
x1 0.154 0.036 0.613 0.492 0.508 1.13
x2 0.011 -0.117 0.525 0.251 0.749 1.1
x3 -0.108 0.028 0.696 0.451 0.549 1.05
x4 0.845 0.007 0.009 0.723 0.277 1
x5 0.902 0.007 -0.081 0.766 0.234 1.02
x6 0.807 -0.013 0.073 0.697 0.303 1.02
x7 0.049 0.741 -0.219 0.497 0.503 1.18
x8 -0.051 0.728 0.054 0.542 0.458 1.02
x9 0 0.483 0.347 0.469 0.531 1.82

重み付き最小二乗法

WLS1 WLS2 WLS3 h2 u2 com
x1 0.156 0.038 0.603 0.482 0.518 1.14
x2 0.011 -0.118 0.528 0.253 0.747 1.1
x3 -0.110 0.028 0.701 0.456 0.544 1.05
x4 0.846 0.007 0.007 0.723 0.277 1
x5 0.896 0.006 -0.080 0.756 0.244 1.02
x6 0.806 -0.013 0.071 0.693 0.307 1.02
x7 0.047 0.746 -0.222 0.502 0.498 1.18
x8 -0.051 0.720 0.056 0.531 0.469 1.02
x9 0 0.477 0.342 0.459 0.541 1.81

反復主因子法

PA1 PA2 PA3 h2 u2 com
x1 0.157 0.039 0.598 0.477 0.523 1.15
x2 0.010 -0.120 0.531 0.256 0.744 1.1
x3 -0.109 0.028 0.699 0.453 0.547 1.05
x4 0.850 0.007 0.005 0.728 0.272 1
x5 0.894 0.005 -0.078 0.753 0.247 1.02
x6 0.805 -0.013 0.072 0.692 0.308 1.02
x7 0.047 0.754 -0.226 0.512 0.488 1.19
x8 -0.050 0.714 0.059 0.524 0.476 1.02
x9 0 0.477 0.344 0.461 0.539 1.82

どれもそう大きく数値的に違いはなさそうである。じゃあ、どれ使えばいいんだよという声があるかもしれないが、清水先生のHPによれば、ファーストチョイスは最尤法で良いらしい。サンプルサイズが小さくて最尤法をあきらめるとなったときに、最小二乗法、反復主因子法などが候補にあがってくるそうだ。ここらへんは、不適解が出るかであったり、収束するかであったりを見ながら決めていくと良いとのことである。

おまけ

調べ物をする最中で次の論文は良さそうだと感じた。

Current Methodological Considerations in Exploratory and Confirmatory Factor Analysis

Schmittという人が2011年に書いたもの。Abstractには、因子分析を行う際には色々決めなければあると紹介した後に次のように書いている。

Unfortunately, researchers continue to use outdated methods in each of these areas.The present article provides a current overview of these areas in an effort to provide researchers with up-to-date methods and considerations in both exploratory and confirmatory factor analysis.

サンプルサイズ・因子数の決定・推定法・回転など何をベースに考えて行けば良いのか注意書きのようなものがあり、因子分析を使う人にとっては役に立つ内容だろう。

おまけ2

Rで出力した表をMarkdownに表記にしたいと思ったときには次のサイトが便利だった。

Markdown Tables Generator

エクセルとかナンバーズにコピペした後に、このサイトに貼り付けるとMarkdown表記にしてくれる。

Rでの大規模データ処理:ffパッケージを試してみる

Rは大変便利な言語なのだけれど、大規模なデータの扱いは苦手と言われている。その理由は、データをメモリに全部載せて作業するので、メモリの量によって扱える操作が制限されるとういもの。最近やや大規模なデータを扱う必要があってその処理の仕方を勉強しているのでその覚書。

ffパッケージ

大規模データを扱うためのパッケージにffというものがある。このスライドThe ff packageとかこのスライド Managing large datasets in R – ff examples and conceptsを参考に試してみる。日本語で書かれているものでff packageの使い方も大変参考になった。

ffはメモリに全部保存するかわりに、データそのものはストレージにおいておき必要なときに必要な分だけ出せるようにする仕組みらしい。

ちなみに、Rは整数だと2^{31}-1までのデータしか表示できない。次のように打つと華麗にエラーを吐き出す。

> as.integer(2^31-1)
[1] 2147483647
> as.integer(2^31)
[1] NA
 警告メッセージ:
NAs introduced by coercion to integer range

ffはこの問題も解決してくれるらしい。

ffの仕組み

ffのパッケージは、データの入れものとして機能する新しいオブジェクトの型が導入されるようである。ffC++で書かれる層とRで書かれる層の2層からなっているそうだ。

Rの層では普段のベクトルやデータフレームと同じ要領で操作ができるようなので色々と操作をしてみる。ff() という関数にlengthを指定してやればベクトルが生成できるらしい。いくつか作ってみる。

# ベクトルの生成
library(ff)
a <- ff(0, length=1)  # 何も指定しないとデータの最初をみて型を判断
b <- ff(0, length=1, vmode="integer") # vmodeで型の明示も可能
c  <- ff(2^31, length=1)  #大きいデータも扱える
d <- ff(1:1000, length=1000) # ベクトルだってOK
e <- ff(1:1e+8, length=1e+8) # 長さ一億のベクトルだってOK

f:id:nekomosyakushimo:20180708151329p:plain

ベクトルへの数値の代入

数値の代入は普通のRと同じようにできる。例えば、先ほどのaというオブジェクトに数字を代入したければ次のようにすれば良い。

a[1] <- 9999

基本的にobject名[]を使って中の要素にアクセスできるのだけれど不用意に大規模データの中身にアクセスしようとする大量のベクトルデータが返ってきて酷い目に合う。

サイズの違い

ffのオブジェクトは、中に大きなデータが入っていてもRの上では必要な分だけデータをとってくるので、オブジェクトのサイズが大きく違う。簡単な数値で実験してみると違いが分かる。

# Rで一様分布
g <- double(10000000)
g[] <- runif(10000000)

# ffで一様分布
h<- ff(vmode="double", length=10000000)
for (i in chunk(h)) h[i] <- runif(sum(i))

print(object.size(g),units="auto")  # 76.3 Mb
print(object.size(h),units="auto")  # 3.1Kb

試しに一様乱数で1,000,000のデータを発生させてみたのだが、普通のRのオブジェクトだと76.3 Mbなのに対して、なんとffオブジェクトは3.1Kbである。ffで生成されたオブジェクトは基本的に必要なときに必要なデータを持ってくるので、実質的に入っているデータがこの10倍になろうがR上では変わらない。

 ファイルの場所

Rのオンメモリのファイルサイズが小さいといってもデータは存在している訳で、デフォルトだと/tmp/Rtmp******とテンポラリーフォルダ内に.ffという拡張子のファイルが生成されている。

オブジェクトの本物のファイルがどこに保存されているかをを確認したいときには、getOption("fftempdir")を用いるとファイルが生成されるフォルダを返す。ためしいにいくつかffオブジェクトを作成すると、その中に.ffの拡張子のファイルが生成されるのが確認できる。

フォルダ内のファイルを確認するとファイルの中身によって大きさが違っており、たとえばさきほど作った長さ1のベクトルが8KBなのに対して、長さ1億のベクトルは400MBのファイルサイズであった。

rm()を用いてr上でffオブジェクトを消してもこの.ffの本体のファイルが消える訳ではないらしい。本体のファルを消したいときにはdelete.ff()と言う関数を用いると関数が開けるかどうかを真偽値で返して、TRUEなのであればそのファイルを削除してくれる。ちなみに、生成した.ffのファイル本体をFinder上で消してからdelete.ff()を実行すると次のようなエラーを返すため、ファイルが開ける場合にのみ削除が実行されるようである。

f:id:nekomosyakushimo:20180708151410p:plain

ファイルの保存場所を指定したいときはoptions(fftempdir="")で指定する。たとえば、options(fftempdir=paste(getwd(),"/fftemp", sep=""))あたりで指定すると、現在の作業ディレクトリの下のフォルダにできたりもする。

行列で実験

ベクトルだけでなく行列も作ってみる。基本的にdimを指定してやると行列になるっぽい。

# 行列の作成
i <- ff(1:12, dim=c(3,4))
t(i) # 転置もできるぞ
j <- ff(0, dim=c(10000000,100)) # 大きい行列だってOK

データフレームの作成

データフレームはffdf()という関数でffのベクトルをくっつけるとできるようだ。

# データフレームの作成
m <- ff(0, vmode="integer",length=100)
n <- ff(0, vmode="double", length=100)
o <- ffdf(m,n)
colnames(o) <- c("A", "B")  # 列名をつけるのもOK

f:id:nekomosyakushimo:20180708151429p:plain

データの読み込みと書き出し

ffdfオブジェクトでデータフレームの読み書きにはwrite.csv.ffdf()read.csv.ffdf()という関数を使うと良いようだ(csvで扱う場合)。

# データフレームの読み書き
write.csv.ffdf(o, file="ff_sample.csv")
p <- read.csv.ffdf(file="ff_sample.csv", fileEncoding="utf-8", header=TRUE)

ここら辺もRで普段データフレームを扱っているのと同じ要領で処理できる。

とりあえずここら辺までで、また色々と調べたら書き足したい。

いきなり因子分析(その2):MAPやBICや平行分析による因子数の決定

以前に書いた記事の続き。

nekomosyakushimo.hatenablog.com

前回は固有値を求めスクリープロットを書いて因子数を決定したが、因子数の決め方にも色々あるらしいので、それを試してみる。参考にしたのは関西学院大の清水先生の次のスライド。この記事の内容はほぼスライドにしたがっただけなので、詳しい解説が読みたい方はスライドをみると良い。

Rで因子分析 商用ソフトで実行できない因子分析のあれこれ

あと、これらの方法の説明については堀先生による次の論文も参考にした。

因子分析における因子数決定法

因子数の決め方いろいろ

因子数を決める基準も様々にあるそうだ。前回の記事では、ガットマン基準とスクリープロット基準の合わせ技のような感じで決めたと思う。

カイザー・ガットマン基準

相関行列の固有値が1を超えるものを選ぶ方法。

スクリープロット基準

固有値をプロットして、大きさの変化がなだらかになる直前の因子数を選ぶもの。

MAP基準

MAPとはMinimum Average Partialの略で日本語にすると最小偏相関平均らしい。主成分を統制変数におき、観測変数間の偏相関係数を求めて、その2乗平均を最小にする主成分数を因子数にするものらしい。

情報量基準

ざっくりと言えばモデルのデータへのあてはまりとモデルの単純さを組み合わせたものが情報量基準。データにほどよくあてはまってほしいけれど、複雑になりすぎてほしくないみたいな話。AICだのBICだの色々あるのだけれど、スライドだとBICをオススメしている。(BICの方がパラメータ数のペナルティが大きいようだ)

平行分析

アイディアとしてはカイザー・ガットマン基準の改定。これはデータと同じサンプルサイズと変数の数で乱数を発生させて、その相関行列の固有値の推定を行い、乱数データの固有値の方が大きくなる一つ前の因子までとる方法らしい。乱数データを扱うことで、標本誤差の影響を考慮できるそうだ。

MAPとBICを算出

という訳で実際にやってみる。使用するデータは前回と同じでHolzingerSwineford1939を使用する。

MAPで求めるときにはpsychパッケージのvss()という関数を使うらしい。引数にはデータと、推定する因子数を予想より多めに入力するそうだ。

library(psych)
vss(dat[,7:15], n=8)

f:id:nekomosyakushimo:20180705082403p:plain

まずMAPについてだが、中断のあたりにThe Velicer MAP achieves a minimum of 0.06 with 2 factors と書いてありMAP基準だと2因子を提案している。

そのすぐ下にBIC achieves a minimum of NA with 3 factorsと書いてありBICによる情報量基準だと3因子が提案されている。

下の方には、因子数ごとのガチャガチャとした統計量が出ているが、MAPという列に注目すると、たしかに因子数2の行が最小の値になっている。また、BICという列に注目すると因子数3の行がやはり最小になっている。

平行分析

こんどは平行分析の方をやってみる。fa.parallel()という関数を使う。スライドをみると、引数にリストワイズしたサンプルサイズを指定する必要があるそうだが、使っているデータは欠損を含まないので、そのままサンプルサイズの301を入力する。n.iterというところはとりあえずスライドと同じ50にしておいた。関数のhelpをみるとデフォルトは20。

fa.parallel(cor(dat[,7:15]),
            n.obs=301,
            n.iter=50
            )

f:id:nekomosyakushimo:20180705082420p:plain

コンソール画面にはシンプルに推奨する因子数がでてくる。

f:id:nekomosyakushimo:20180705082459p:plain

そして、前回の記事に書いたスクリープロットにいくらか線が加わったものが作図される。PCというのが主成分分析でFAというのが因子分析の結果である。赤い色の点線と交差した、1つ前の因子数/主成分数が提案される数ということで、今回の場合にはどちらも3つを提案している。

ちなみに、主成分のACUTUAL DATAの線は前回、固有値をプロットしたスクリープロットと全く同じである。

因子数の決定基準

様々な基準があったが推奨する因子数が全て一致すればそれを採用すれば良いらしい。違う場合には、MAPはマイナー因子を拾わない傾向があり、平行分析はやや多めに因子数を提案するようなので、理論や解釈可能性をベースにその間で決めていくと良いらしい。

ちなみに今回の分析をまとめると次の表のようになるので、このデータについては因子数3を採用するのが良いと言えそうである。

基準 因子数
ガットマン 3
スクリー 3
MAP 2
BIC 3
平行分析 3