確証的因子分析の際のカイ二乗値はどこからくるの？

確証的因子分析を行うと、モデルがどの程度データに当てはまっているかの指標が複数算出される。その解釈の仕方については検索するとよく出てくるが、その算出については調べてもあまり出てこないので手元で色々と試してた記録を残しておく。

確証的因子分析を実行

何はともあれ、確証的因子分析をとりあえず実行する。データはお得意のlavaanのパッケージに入っているHolzingerSwineford1939を使う。

dat <- HolzingerSwineford1939
model1 <- "
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
"
fit <- cfa(model = model1,
            data = dat
           )
summary(fit, fit.measures=T)

こんな感じでやると、適合度の指標がたくさん出てくるのだが、最初の方に次のような値が出力される。

Minimum Function Test Statistic               85.306
Degrees of freedom                                24

これが指標のうちの一つであるカイ二乗値である。カイ二乗値はモデル評価の指標として問題もあるので、これ単体にでモデルの評価をしようと勧める文献は見ないのだが、その他の指標を出す際に使われたりもするので、これがどう算出されているかを今回の記事では扱う。

カイ二乗値はどこから？

カイ二乗値がどのように算出されているかというと『因子分析入門』（東京図書）という本のp.105には次のような式が乗っている。

f:id:nekomosyakushimo:20180722151440p:plain

ここで、 $S$ は標本共分散行列、 $\hat{\Sigma}$ は母数で再現された共分散行列、 $n$ は観測対象の数、 $p$ は観測変数の数である。

これらの行列を得ることができれば、この式通りに計算して、先ほど算出されたカイ二乗値と同じ値になるはずである。さっそくやってみよう。

lavaanの結果から必要な行列を取り出す

必要な行列を得る。標本の共分散行列は普通にcov()の関数で求めれば問題ない。問題は、母数で再現された共分散行列であるが、cfa()を実行した結果オブジェクトを、fitted.values()という関数に与えると、推定値による共分散行列と推定値の平均ベクトルを得ることができるのでそれを用いて推定値による共分散行列を取得する。

sample_cov <- cov(dat[7:15])  #標本共分散分散行列
sigma_hat <- fitted.values(fit)$cov 　#推定された共分散行列

計算式にあてはめ

これらの行列を用いて式通りに計算していく。部分に分けて計算していこう。

f:id:nekomosyakushimo:20180722151505p:plain

まず、Aの部分から。Rでは逆行列を求めるにはsolve()を用いると良い。トレースとは、行列の対角成分の和のことであるが、Rだとトレースを一発で求める関数はないのでdiag()で対角成分を抜き出したあとにsum()でそれらを足すという手順をとる。

A = sum(diag(sample_cov %*% solve(sigma_hat)))

次はBの部分。行列式はdet()で求める。自然対数をとるにはlog()を用いる。

B = log(det(solve(sigma_hat) %*% sample_cov))

あとは、これらを組み合わせて計算するのみである。

n.obs <- 301  #観察対象の数
n.variables <- 9  # 観測変数の数
A <- sum(diag(sample_cov %*% solve(sigma_hat)))
B <- log(det(solve(sigma_hat) %*% sample_cov))

chisquare <- (n.obs-1) * (A - B - n.variables)
print(chisquare)

[1] 85.03708

微妙にずれている。そこで、計算式の(n.obs - 1)をn.obsに直して計算すると一致した。

chisquare1 <- (n.obs) * (A - B - n.variables)
print(chisquare1)
[1] 85.32054

『因子分析入門』のp.218には、カイ二乗値の定義にもいくつか流儀があると書いてあり、n-1をかける場合もnをかける場合も紹介されている。lavaanのパッケージでは、nをかけて算出しているようである。

独立モデルの場合

モデルの適合度を調べるときに、あてはまりの最も悪いモデルの値を算出して、そのモデルからどれだけ改善が見られるかという考え方をする指標がある。このときに、使われる観測変数間に一切パスを引かないモデルを独立モデルという。最初に実施したcfa()の結果だと次のように表示されている部分である。

Model test baseline model:

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

ついでなので、独立モデルのカイ二乗値を計算値に当てはめて求めてみる。独立モデルは対角行列であり、その最尤推定量は常に $\Sigma(\hat{\theta})$ であるので（『共分散構造分析[入門編]』（朝倉書店）p.176）、それを求めて計算してみよう。

標本共分散行列の非対角成分を0にするにはupper.tri()をlower.tri()用いて、そこの部分に0を代入してやると良い。

baseline <- sample_cov
baseline[upper.tri(baseline)] <- 0
baseline[lower.tri(baseline)] <- 0

そしてそれを先ほどの計算式に代入して計算する。

n.obs <- 301
n.variables <- 9
A <- sum(diag(sample_cov %*% solve(baseline)))
B <- log(det(solve(baseline) %*% sample_cov))
chisquare_bl <- (n.obs) * (A - B - n.variables)

print(chisquare_bl)
[1] 918.8516