重み付き最小二乗法についての調べ物

SEMなどの際の推定法の中にWeighted Least Square（重み付き最小2乗法）というものがある。通常の（重みをつけない）最小2乗法は標本共分散行列とモデルの共分散行列の残差が最も小さくなるようにモデルの母数を推定するが, WLSでは残差に対して異なる重みを与える。

Weighted Least Squareとは

SEMのオプションで一般にWLSというとBrown(1984)によるAsymptotically Distribution-Free法のことを指すようで, 分布によらない推定ができるらしい。データとモデルの不一致度の関数は次のように表される。

$\displaystyle F(\boldsymbol{S, \Sigma}) = (\boldsymbol{s - \sigma})' \boldsymbol{W}^{-1} (\boldsymbol{s-\sigma})$

ここで, $\boldsymbol{s}$ は標本共分散行列の対角要素を含む下三角要素を縦ベクトルにしたもので $\boldsymbol{\sigma}$ は, 同じ操作をモデルの共分散行列にしたものである。 $\boldsymbol{W}$ が重みの行列である。

重みの行列はサイズが, 観測変数の数を $p$ としたときに $p(p+1)/2 \times p(p+1)/2$ になる。例えば観測変数が2個だとしたら $3 \times 3$ の行列でそれぞれの成分は次のようになる。

$\displaystyle \begin{pmatrix} w_{11,11} & &sym \\\ w_{21,11} & w_{21,21} & \\\ w_{22,11} & w_{22,21} & w_{22,22} \end{pmatrix}$

ここで, $w_{ij,kl}$ は標本数を $N$ , 観測変数を $x$ としたときに次のようになる。

$\displaystyle w_{ij,kl} = s_{ij,kl} - s_{ij} s_{il} \\ s_{ij,kl} = \frac{1}{N-1} \Sigma (x_{in} - \bar{x_i})(x_{jn} - \bar{x_j})(x_{kn} - \bar{x_k})(x_{ln} - \bar{x_l})$

である。

lavaanで試してみる

具体的な数字を入れてみると分かりやすいだろうということでみんな大好きlavaanで試してみる。まず, データの変数の部分だけ取り出してモデルを指定。

library(lavaan)
dat <- HolzingerSwineford1939[,7:15]

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

推定法にestimator="WLS"とつけるとlavaanはWLSでやってくれる。

res.wls <- cfa(model=model1, data=dat, estimator="WLS")
summary(res.wls)

f:id:nekomosyakushimo:20190511234844p:plain

さて, このモデルで重み行列はどこに保存されているかというとres.wls@SampleStats@NACOVというところにある。9変数あるので45×45の行列である。

f:id:nekomosyakushimo:20190511234850p:plain

表示しきれていないだけで, 下にまだまだ続いている。

さて, 実際にさっきの計算式通りに計算して確かめてみよう。たとえば $w_{21,11}$ を計算してみる。lavaanでは, N-1でなくNで計算しているようなので出力結果に合うように300/301をかけている。

W <- res.wls@SampleStats@NACOV[[1]] # lavaanの出力結果のW
mean <- res.wls@SampleStats@mean[[1]] # 標本の各変数の平均
s <- res.wls@SampleStats@cov[[1]] # 標本の分散共分散行列


s_2111 <- mean((dat[,2]- mean[2])*(dat[,1] -mean[1])^3)
W_2111 <- s_2111 - (s[1,1]*300/301 * s[2,1]*300/301)