多変量解析_07_構造方程式モデリング

Chapter 7 構造方程式モデリング Chapter 7 構造方程式モデリング 本資料は，神戸大学経営学研究科で 2022 年度より担当している「統計的方法論特殊研究（多 変量解析）」の講義資料です。CC BY-NC 4.0 ライセンスの下に提供されています。 作成者連絡先 神戸大学大学院経営学研究科 分寺杏介（ぶんじ・きょうすけ） mail: [email protected] website: https://www2.kobe-u.ac.jp/~bunji/ 今回からは構造方程式モデリング (Structural Equation Modeling; SEM) の話です。はじ めに SEM の基本的な考え方を紹介し，R で実際に実行してみます。 Contents 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 7.10 7.11 7.12 7.13 7.14 因子分析の見方を変える . . . . . . . . . . . . . . . . . . . . . . . . . R でやってみる（検証的因子分析） . . . . . . . . . . . . . . . . . . . 回帰分析の見方を変える . . . . . . . . . . . . . . . . . . . . . . . . . R でやってみる（回帰分析・パス解析） . . . . . . . . . . . . . . . . . 回帰分析と因子分析を組み合わせる . . . . . . . . . . . . . . . . . . . . モデルの適合度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . モデルの修正 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . カテゴリカル変数の SEM . . . . . . . . . . . . . . . . . . . . . . . . 多母集団同時分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lavaan の記法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 同値モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 𝜔ℎ を求める . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . もう一つの SEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . （おまけ）lavaan で探索的因子分析ができるようになっていた件について . . . . . . . . . . . . . . 2 8 16 19 23 29 37 50 57 70 81 83 85 91 （いつにもましてページが多いですが，lavaan の出力が多いため…だと思います…） R での SEM は lavaan というパッケージを用いるのが最もメジャーです。事前にインスト ール＆読み込みをしておきましょう。合わせて，このチャプターで使用するいくつかのパッ ケージもインストールしておきます。 1

• https://creativecommons.org/licenses/by-nc/4.0/
• https://www2.kobe-u.ac.jp/~bunji/

7.1 因子分析の見方を変える 事前準備 1 # まだインストールしていない人は 2 # install.packages("lavaan") 3 # install.packages("semPlot") 4 # install.packages("semTools") 5 6 # 前半ではlavaanのみ使います。 7 library(lavaan) 8 dat <- readRDS("chapter04.rds") 7.1 因子分析の見方を変える 7.1.1 1 因子モデル 図7.1 因子分析モデル まずは図7.1のような１因子３項目の因子分析モデルを考えてみます。このモデルでは，各観 測変数はそれぞれ ⎧𝐲1 = 𝐟 𝑏1 + 𝛆1 { 𝐲 = 𝐟 𝑏 2 + 𝛆2 ⎨ 2 {𝐲 = 𝐟 𝑏 + 𝛆 ⎩ 3 3 3 (7.1) と表されました。したがって，３つの観測変数のモデル上の共分散行列 𝚺 は 𝑏12 + 𝜎𝛆2 1 ⎡ 𝚺 = ⎢ 𝑏2 𝑏1 ⎣ 𝑏3 𝑏1 𝑏22 𝑏1 𝑏2 + 𝜎𝛆2 2 𝑏3 𝑏2 𝑏1 𝑏3 𝑏2 𝑏3 ⎤ ⎥ 2 𝑏3 + 𝜎𝛆2 3 ⎦ (7.2) となります。また，３つの観測変数のデータでの共分散行列も 𝜎𝑦21 ⎡ 𝐒 = ⎢𝜎𝑦2 ,𝑦1 ⎣𝜎𝑦3 ,𝑦1 𝜎𝑦1 ,𝑦2 𝜎𝑦22 𝜎𝑦3 ,𝑦2 𝜎𝑦1 ,𝑦3 𝜎𝑦2 ,𝑦3 ⎤ ⎥ 𝜎𝑦23 ⎦ (7.3) と表すことができ，これは実際にデータから計算が可能です。例えば Q1_1 から Q1_3 の３ 2

7.1 因子分析の見方を変える 項目のみを考えれば cov(dat[,paste0("Q1_",1:3)]) によって*1 具体的に 1.976 0.581 𝐒=⎡ ⎢0.581 1.375 ⎣0.504 0.762 0.504 0.762⎤ ⎥ 1.706⎦ という値が得られます。構造方程式モデリングでは，𝚺 の値が 𝐒 にできるだけ近づくように 母数（パラメータ）を推定します。実際に 𝚺 と 𝐒 の同じ項をイコールで結ぶと， 2 ⎧ 𝜎 𝑦1 { { 𝜎𝑦22 { { 𝜎2 𝑦3 ⎨𝜎 𝑦 ,𝑦 { 1 2 {𝜎 { 𝑦1 ,𝑦3 { ⎩𝜎𝑦2 ,𝑦3 = 𝑏12 + 𝜎𝛆2 1 = 𝑏22 + 𝜎𝛆2 2 = 𝑏32 + 𝜎𝛆2 3 = 𝑏 1 𝑏2 (7.4) = 𝑏 1 𝑏3 = 𝑏 2 𝑏3 と い う 連 立 方 程 式 を 得 る こ と が 出 来 ま す。 こ の モ デ ル で 推 定 す る べ き パ ラ メ ー タ は (𝑏1 , 𝑏2 , 𝑏3 , 𝜎𝛆2 1 , 𝜎𝛆2 2 , 𝜎𝛆2 3 ) の６つです。そして式も６つあるため，このケースでは明確な 解を求めることが出来ます。例えば Q1_1 から Q1_3 の３項目であれば 2 2 ⎧1.976 = 𝑏1 + 𝜎𝛆1 { {1.375 = 𝑏22 + 𝜎𝛆2 2 { {1.706 = 𝑏2 + 𝜎2 𝛆3 3 ⎨0.581 = 𝑏 𝑏 1 2 { {0.504 = 𝑏 𝑏 { 1 3 { ⎩0.762 = 𝑏2 𝑏3 となり，これを解くと ⎧ 𝑏1 { 𝑏 { 2 { { 𝑏3 ⎨𝜎𝛆2 1 { {𝜎2 { 𝛆2 { 2 ⎩𝜎𝛆3 = 0.620 = 0.937 = 0.813 = 1.591 = 0.497 = 1.045 という解が得られます。図7.1に書き加えたならば図7.2のようになります。 グラフィカルモデルでは，観測変数を四角形で，潜在変数を楕円で表していました。図7.2と パラメータの関係を見てみると，2 つの観測変数の共分散は，2 つの四角形を結ぶ矢印にか かる係数の積になっていることがわかります。つまり，グラフィカルモデル上で経路が存在 する観測変数間には非ゼロの共分散があることを仮定していたわけです。また，1 つの観測 変数の分散は，その四角形に向かっている矢印の係数の二乗和になっています。なお，独自 因子には自分自身に対する両方向の矢印が引かれています。グラフィカルモデルでは両方向 の矢印で共分散を表すので，これは自分自身との共分散＝分散を表していることになります。 先程の結果は観測変数を標準化していない状態での解でした。そのため，独自因子の分散が 1 を超えたり，因子負荷の二乗との合計が 1 にならなかったりしています。例えば項目 𝐲1 に *1 3 実際には cov() 関数で得られるのは不偏共分散行列です。𝐒 は標本共分散行列なので，𝑃 /𝑃 − 1 倍する必 要があったりするのですが，今回は 2000 人以上いるので気にしないことにします。

7.1 因子分析の見方を変える 図7.2 因子分析モデルの解 ついて見てみると，因子負荷の二乗と独自因子の分散の和は 0.6202 + 1.591 ≃ 1.976 となり ます。これは，𝐲1 の分散 𝜎𝑦21 の値 1.976 に対応しているわけです。ここで前回までの因子分 析と同じように全ての観測変数を標準化した状態，つまり 𝐒 に共分散行列ではなく相関行列 を与えて推定をした場合に得られる解のことを標準化解や標準化推定値などと呼びます。そ してSEMでは多くの場合，結果には標準化解を載せています。絶対のルールでは無いのです が，特に潜在変数がある場合には単位の影響を受けない標準化解を用いたほうがしっくり来 る，ということなのでしょう。ということでこの先も，潜在変数が含まれるモデルを扱う場 合には 𝐒 には相関行列を想定した，標準化解を考えることにします。 実 際 に 先 程 の 例 で 標 準 化 解 を 計 算 し て み る と， デ ー タ に お け る 相 関 行 列 𝐒 は cor(dat[,paste0("Q1_",1:3)]) によって具体的に*2 1 0.353 𝐒=⎡ ⎢ ⎣0.274 0.353 1 0.497 0.274 0.497⎤ ⎥ 1 ⎦ という値が得られるため，連立方程式は 1 = 𝑏12 + 𝜎𝛆2 1 ⎧ { { 1 = 𝑏22 + 𝜎𝛆2 2 { { 1 = 𝑏32 + 𝜎𝛆2 3 ⎨0.353 = 𝑏 𝑏 1 2 { {0.274 = 𝑏 𝑏 { 1 3 { 0.497 = 𝑏 ⎩ 2 𝑏3 となり，これを解くと ⎧ 𝑏1 { 𝑏 { 2 { { 𝑏3 ⎨𝜎𝛆2 1 { {𝜎 2 { 𝛆2 { 2 ⎩𝜎𝛆3 = 0.441 = 0.800 = 0.622 = 0.805 = 0.361 = 0.613 という標準化解を得ることができます。 *2 4 もちろんこれも厳密にはポリコリック相関のほうが良い可能性が高いですが，ここでは深く考えずにとりあ えずピアソンの積率相関係数で話を進めます。

7.1 因子分析の見方を変える 7.1.2 2 因子モデル 図7.3 2 因子モデル 続いて図7.3のような 2 因子 6 項目の因子分析モデルで同じように見てみましょう。ここで は，因子 1 は最初の 3 項目にのみ，因子 2 は後ろの 3 項目にのみ影響を与えているという仮 定を置いています。また，因子間相関 𝑟𝑓1,𝑓2 の間には両矢印を描いています。このように， （共）分散を表す時には両矢印を使います。つまりここでは因子間の共分散（＝相関）が 0 で なくても良いという仮定（斜交回転）を置いているのです。 このモデルでは，因子負荷行列 𝐁⊤ は 𝑏11 ⎡𝑏 ⎢ 21 𝑏 𝐁⊤ = ⎢ 31 ⎢0 ⎢0 ⎣0 0 0⎤ ⎥ 0⎥ 𝑏42 ⎥ 𝑏52 ⎥ 𝑏62 ⎦ (7.5) となります。ちなみに SEM で多因子モデルを扱う場合の多くは，このように 1 項目が 1 因 子からのみ影響を受けていると考えると思います*3 。 *3 5 理論的には 𝑇 因子構造の場合，因子負荷行列のうち最低でも𝑇 (𝑇 − 1) 個の値を固定したら推定は可能です。

7.1 因子分析の見方を変える すると，6 つの観測変数のモデル上の相関行列 𝚺 は 𝚺 = 𝐁⊤ 𝚽𝐁 + 𝚿 𝑏11 ⎡𝑏 ⎢ 21 𝑏 = ⎢ 31 ⎢0 ⎢0 ⎣0 0 0⎤ ⎥ 1 𝑟𝑓1,𝑓2 𝑏11 𝑏21 𝑏31 0 0⎥ 0 0 [ ][ ]+𝚿 𝑏42 ⎥ 𝑟𝑓1,𝑓2 1 0 0 0 𝑏42 𝑏52 𝑏62 𝑏52 ⎥ 𝑏62 ⎦ 2 𝑏11 𝑏21 𝑏11 𝑏31 𝑟𝑓1,𝑓2 𝑏11 𝑏42 𝑟𝑓1,𝑓2 𝑏11 𝑏52 𝑏11 2 ⎡ 𝑏 𝑏 𝑏 𝑏 𝑏 𝑟𝑓1,𝑓2 𝑏21 𝑏42 𝑟𝑓1,𝑓2 𝑏21 𝑏52 21 11 21 31 21 ⎢ 2 𝑏 𝑏 𝑏 𝑏 𝑏 𝑟𝑓1,𝑓2 𝑏31 𝑏42 𝑟𝑓1,𝑓2 𝑏31 𝑏52 31 11 31 21 31 =⎢ 2 𝑟 𝑏 𝑏 𝑟 𝑏 𝑏 𝑟 𝑏 𝑏 𝑏42 𝑏42 𝑏52 ⎢ 𝑓1,𝑓2 42 11 𝑓1,𝑓2 42 21 𝑓1,𝑓2 42 31 2 ⎢𝑟𝑓1,𝑓2 𝑏52 𝑏11 𝑟𝑓1,𝑓2 𝑏52 𝑏21 𝑟𝑓1,𝑓2 𝑏52 𝑏31 𝑏52 𝑏42 𝑏52 𝑏62 𝑏42 𝑏62 𝑏52 ⎣𝑟𝑓1,𝑓2 𝑏62 𝑏11 𝑟𝑓1,𝑓2 𝑏62 𝑏21 𝑟𝑓1,𝑓2 𝑏62 𝑏31 𝜎𝛆2 0 0 0 0 0 ⎤ ⎡ 0 1 𝜎2 0 0 0 0 𝛆 2 ⎥ ⎢ 2 0 𝜎 𝛆3 0 0 0 ⎥ ⎢ 0 +⎢ 0 0 𝜎𝛆2 4 0 0 ⎥ ⎥ ⎢ 0 2 0 ⎥ 0 0 0 𝜎 𝛆5 ⎢ 0 0 0 0 0 𝜎𝛆2 6 ⎦ ⎣ 0 2 𝑏11 + 𝜎𝛆2 𝑏11 𝑏21 𝑏11 𝑏31 2 2 ⎡ 𝑏 𝑏 1 𝑏21 𝑏31 𝑏21 + 𝜎𝛆2 21 11 ⎢ 2 𝑏31 𝑏21 𝑏31 + 𝜎𝛆2 3 ⎢ 𝑏31 𝑏11 =⎢ ⎢𝑟𝑓1,𝑓2 𝑏42 𝑏11 𝑟𝑓1,𝑓2 𝑏42 𝑏21 𝑟𝑓1,𝑓2 𝑏42 𝑏31 ⎢𝑟𝑓1,𝑓2 𝑏52 𝑏11 𝑟𝑓1,𝑓2 𝑏52 𝑏21 𝑟𝑓1,𝑓2 𝑏52 𝑏31 ⎣𝑟𝑓1,𝑓2 𝑏62 𝑏11 𝑟𝑓1,𝑓2 𝑏62 𝑏21 𝑟𝑓1,𝑓2 𝑏62 𝑏31 𝑟𝑓1,𝑓2 𝑏11 𝑏42 𝑟𝑓1,𝑓2 𝑏21 𝑏42 𝑟𝑓1,𝑓2 𝑏31 𝑏42 2 𝑏42 + 𝜎𝛆2 4 𝑏52 𝑏42 𝑏62 𝑏42 𝑟𝑓1,𝑓2 𝑏11 𝑏52 𝑟𝑓1,𝑓2 𝑏21 𝑏52 𝑟𝑓1,𝑓2 𝑏31 𝑏52 𝑏42 𝑏52 2 𝑏52 + 𝜎𝛆2 5 𝑏62 𝑏52 𝑟𝑓1,𝑓2 𝑏11 𝑏62 𝑟𝑓1,𝑓2 𝑏21 𝑏62 ⎤ ⎥ 𝑟𝑓1,𝑓2 𝑏31 𝑏62 ⎥ 𝑏42 𝑏62 ⎥ 𝑏52 𝑏62 ⎥ 2 𝑏62 ⎦ 𝑟𝑓1,𝑓2 𝑏11 𝑏62 𝑟𝑓1,𝑓2 𝑏21 𝑏62 ⎤ ⎥ 𝑟𝑓1,𝑓2 𝑏31 𝑏62 ⎥ 𝑏42 𝑏62 ⎥ ⎥ 𝑏52 𝑏62 ⎥ 2 𝑏62 + 𝜎𝛆2 6 ⎦ (7.6) となります。例えば y_1 と y_4 の繋がりを図7.3で見てみると，y_1 → f_1 → f_2 → y_4 という経路になっており，この 2 変数の共分散はそれぞれの矢印に係る係数の積である 𝑟𝑓1,𝑓2 𝑏42 𝑏11 となっていることがわかります。もし因子間に相関がなければ，y_1 と y_4 を 結ぶ経路は存在せず，数式的には 𝑟𝑓1,𝑓2 = 0 となることからも，この 2 変数が完全に無相関 であることを意味するようになります。 このモデルの重要なポイントは，分散共分散行列が 𝐁⊤ 𝚽𝐁 + 𝚿 という形で分解可能だとい う点は前回までの因子分析と全く同じだという点です。SEM の枠組みでは，因子負荷行列 𝐁 や因子間相関 𝚽 に関して要素の値を固定（通常は 0 に固定）するため，それに応じて観測変 数間の共分散の形は変わりますが，違いといえばそれだけです*4 。特定の因子が特定の項目 にのみ負荷を持っているという仮定は，分析者が仮説などに基づいて行うものです。そして 得られた結果（𝚺 と 𝐒 の乖離度）は，その仮説（に基づく共分散構造）がどの程度妥当そう かを表したものといえます。そのような意味で，SEM における因子分析を検証的（確認的） 因子分析 (confirmatory factor analysis; CFA) と呼びます。一方，前回まで行っていた因 子分析は，因子構造を探し当てるのが目的だったため探索的因子分析 (exploratory factor analysis; EFA) と呼んでいます。 *4 6 残念ながら固有値分解によって因子負荷行列を求めることはできなくなりますが，コンピュータによる反復 計算ならば力技で近い値を求めることができるでしょう。

7.1 因子分析の見方を変える 上述のモデルでは，推定するパラメータは全部で 13 個（6 個の因子負荷＋ 6 個の誤差分散＋ 1 個の因子間相関）あり，対する観測変数間の分散共分散の数が全部で 21 個（分散が 6 個＋ 共分散が 15 個）あります。ということで 1 因子 3 項目の場合のように解を一意に求めるこ とは出来ないのですが，なるべく誤差が小さくなるように全てのパラメータをまとめて推定 することは可能なわけです。 7.1.3 もう少しだけ複雑なモデル 図7.4 独立ではない 2 因子モデル CFA では，もちろん「ある項目は 2 つの因子の影響を受けている」というモデルを作ること も可能です。図7.4は，項目 3 だけが 2 つの因子の影響を受けていると仮定しています。この 時，因子負荷行列は 𝑏11 ⎡𝑏 ⎢ 21 𝑏 ⊤ 𝐁 = ⎢ 31 ⎢0 ⎢0 ⎣0 0 0⎤ ⎥ 𝑏32 ⎥ 𝑏42 ⎥ 𝑏52 ⎥ 𝑏62 ⎦ (7.7) となります。このとき共分散行列はどうなるでしょうか。すべて書くのは面倒なので項目 1 と 3 の共分散だけ抜き出してみると 𝜎𝑦1,𝑦3 = 𝑏31 𝑏11 + 𝑟𝑓1,𝑓2 𝑏32 𝑏11 (7.8) という形になっています。図7.4で y_1 と y_3 の経路を見てみると， • y_1 → f_1 → y_3 (𝑏31 𝑏11 ) • y_1 → f_1 → f_2 → y_3(𝑟𝑓1,𝑓2 𝑏32 𝑏11 ) という 2 つの経路が存在しており，それに対応した係数の積の合計がこの 2 変数の共分散を 表すようになっています。また，項目 3 の分散は 2 2 2 𝜎𝑦3 = 𝑏31 + 𝑏32 + 2𝑟𝑓1,𝑓2 𝑏31 𝑏32 + 𝜎𝛆2 3 (7.9) となります。1 因子モデルのときも分散は係数の二乗和になっていましたが，より正確に 言えば 2 • y_3 → f_1 → y_3 (𝑏31 ) 7

7.2 R でやってみる（検証的因子分析） 2 • y_3 → f_2 → y_3 (𝑏32 ) • y_3 → e_3 → y_3 (𝜎𝛆2 3 ) • y_3 → f_1 → f_2 → y_3 (𝑟𝑓1,𝑓2 𝑏31 𝑏32 ) • y_3 → f_2 → f_1 → y_3 (𝑟𝑓1,𝑓2 𝑏31 𝑏32 ) という要領で，その変数を出発してから戻るまでの経路の全パターンに関して同様の計算を していたわけです。 7.2 R でやってみる（検証的因子分析） それでは，図7.3のモデルを実際に SEM の枠組み（検証的因子分析）で推定してみましょう。 lavaan では，lm() などと同じように第一引数にモデル式を与えます。複雑なモデルになる と結構な長さになるので，以下のように事前にモデル式をオブジェクトに入れておくのがお すすめです。 モデル式を用意する 1 model_cfa <- " 2 # 改行とコメントアウトはご自由に 3 4 f_1 =~ Q1_1 + Q1_2 + Q1_3 5 f_2 =~ Q1_6 + Q1_7 + Q1_8 6 " モデル式は character 型で用意するため，ダブルクオート（" "）かシングルクオート（' '） で全体をくくります。その中に，回帰式のようなものを書いていきます。lavaan の構文で は，潜在変数を規定する場合にはイコールチルダ（=~）を使用します。つまり，上のコード の 4 行目では「因子 f_1 は Q1_1 と Q1_2 と Q1_3 に影響を与えているもの」であることを 表しているわけです。 ちなみに観測変数に関してはデータフレームの列名で指定しますが，潜在変数はデータフレ ーム内に存在しないので，自分でわかりやすい名前をつけても構いません。例えば Q1_1 と Q1_2 と Q1_3 は”Agreeableness” を測定する項目だとわかっているので，1 行目は f_A =~ Q1_1 + Q1_2 + Q1_3 としても良いわけです。lavaan では，データフレームに存在しない 変数は自動的に潜在変数であるとみなして分析を行ってくれます。 モデル式が用意できたら，あとは分析を実行するだけです。検証的因子分析の場合には cfa() という関数を利用しましょう*5 。 lavaan で CFA 1 result_cfa <- cfa(model_cfa, data = dat) 2 # 結果を見るときにはlm()と同じようにsummary()にかける *5 8 本当は lavaan パッケージの基本関数は同名の lavaan() というものです。この関数ではモデルの推定のた めの設定をいろいろと細かく変更できるのですが，その設定が割と細い上にデフォルトでは思っていたよう な設定になかなかたどり着けません。いわば上級者向けの関数です。

7.2 R でやってみる（検証的因子分析） 3 summary(result_cfa) 1 lavaan 0.6.15 ended normally after 38 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters ML NLMINB 13 6 7 Number of observations 2432 8 9 Model Test User Model: 10 11 Test statistic 12 Degrees of freedom 13 P-value (Chi-square) 62.061 8 0.000 14 15 Parameter Estimates: 16 17 Standard errors Standard 18 Information Expected 19 Information saturated (h1) model Structured 20 21 Latent Variables: 22 Estimate Std.Err z-value P(>|z|) 23 f_1 =~ 24 Q1_1 1.000 25 Q1_2 1.548 0.106 14.538 0.000 26 Q1_3 1.348 0.082 16.362 0.000 27 f_2 =~ 28 Q1_6 1.000 29 Q1_7 1.220 0.075 16.316 0.000 30 Q1_8 0.894 0.053 16.756 0.000 Estimate Std.Err z-value P(>|z|) 0.115 0.016 7.349 0.000 31 32 Covariances: 33 34 f_1 ~~ 35 f_2 36 37 Variances: 38 9 Estimate Std.Err z-value P(>|z|) 39 .Q1_1 1.607 0.051 31.229 0.000 40 .Q1_2 0.492 0.054 9.151 0.000 41 .Q1_3 1.036 0.050 20.917 0.000 42 .Q1_6 0.950 0.043 22.233 0.000 43 .Q1_7 0.888 0.055 16.193 0.000 44 .Q1_8 1.205 0.044 27.532 0.000

7.2 R でやってみる（検証的因子分析） 45 f_1 0.368 0.040 9.180 0.000 46 f_2 0.567 0.047 11.995 0.000 結果を上から見ていきましょう。 3 Estimator 4 Optimization method 5 Number of model parameters ML NLMINB 13 Estimator のところで，最尤推定法（ML）によってパラメータを求めていることがわかりま す。パラメータは反復計算で求めているのですが，Optimization method はその計算アル ゴリズムを示しています。基本的には触れなくて OK です。Number of model parameters はモデルのパラメータ数です。6 個の因子負荷＋ 6 個の誤差分散＋ 1 個の因子間相関で 13 個ということは先程説明した通りです。 9 Model Test User Model: 10 11 Test statistic 12 Degrees of freedom 13 P-value (Chi-square) 62.061 8 0.000 モデルに対して 𝜒2 検定を行っています。帰無仮説は「モデルの共分散行列とデータの共分 散行列が同じ」ということで，𝑝 値が大きい（帰無仮説が棄却されない）ほど当てはまりが 良いことを意味します。つまり，できればこの検定は棄却される（𝑝 > .05）と嬉しいのです が，統計的仮説検定の仕組み上サンプルサイズが多いほど帰無仮説が棄却されやすかったり となかなか使い勝手の良いものではないので，基本的には無視でも良いでしょう。 21 Latent Variables: 22 Estimate Std.Err z-value P(>|z|) 23 f_1 =~ 24 Q1_1 1.000 25 Q1_2 1.548 0.106 14.538 0.000 26 Q1_3 1.348 0.082 16.362 0.000 27 f_2 =~ 28 Q1_6 1.000 29 Q1_7 1.220 0.075 16.316 0.000 30 Q1_8 0.894 0.053 16.756 0.000 因子負荷です。lavaan はデフォルトでは，パラメータの制約として「因子の分散が 1」の代 わりに「各因子の第 1 項目に対する因子負荷が 1」という制約を置いています。因子分析の 資料の 3 ページで潜在変数のスケールの不定性を説明しましたが，lavaan のデフォルトで は因子ごとに最初の項目の 𝑏 の値を一つ固定することによってその不定性を封じているわけ です。ですが一般的には「因子の分散が 1」の方がよく用いられる制約なので，そのように した場合の推定値も出力してみます。その方法は 2 種類あり，1 つ目の方法は「推定時に引 10

7.2 R でやってみる（検証的因子分析） 数std.lv = TRUEを与える」というものです。 推定時に潜在変数の分散を 1 にする 1 result_cfa <- cfa(model_cfa, data = dat, std.lv = TRUE) 2 summary(result_cfa) 1 lavaan 0.6.15 ended normally after 27 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters ML NLMINB 13 6 7 Number of observations 2432 8 9 Model Test User Model: 10 11 Test statistic 12 Degrees of freedom 13 P-value (Chi-square) 62.061 8 0.000 14 15 Parameter Estimates: 16 17 Standard errors Standard 18 Information Expected 19 Information saturated (h1) model Structured 20 21 Latent Variables: 22 Estimate Std.Err z-value P(>|z|) 23 f_1 =~ 24 Q1_1 0.607 0.033 18.359 0.000 25 Q1_2 0.939 0.034 27.729 0.000 26 Q1_3 0.818 0.034 24.120 0.000 27 f_2 =~ 28 Q1_6 0.753 0.031 23.990 0.000 29 Q1_7 0.919 0.035 26.112 0.000 30 Q1_8 0.673 0.031 21.419 0.000 Estimate Std.Err z-value P(>|z|) 0.251 0.028 8.852 0.000 Estimate Std.Err z-value P(>|z|) 31 32 Covariances: 33 34 f_1 ~~ 35 f_2 36 37 Variances: 38 11 39 .Q1_1 1.607 0.051 31.229 0.000 40 .Q1_2 0.492 0.054 9.151 0.000

7.2 R でやってみる（検証的因子分析） 41 .Q1_3 1.036 0.050 20.917 0.000 42 .Q1_6 0.950 0.043 22.233 0.000 43 .Q1_7 0.888 0.055 16.193 0.000 44 .Q1_8 1.205 0.044 27.532 0.000 45 f_1 1.000 46 f_2 1.000 std.lv = TRUE とすることで，潜在変数 (latent variable) を標準化 (standardized) するこ とができます。また，もう一つの方法は「summary()時に引数standardized = TRUEを与 える」というものです。この方法では因子の分散を 1 に調整した場合の推定値を事後的に変 換して出力してくれます。 事後的に推定値を変換する 1 # 説明のためstd.lv = TRUEを消して再推定 2 result_cfa <- cfa(model_cfa, data = dat) 3 summary(result_cfa, standardized = TRUE) 1 lavaan 0.6.15 ended normally after 38 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters ML NLMINB 13 6 7 Number of observations 2432 8 9 Model Test User Model: 10 11 Test statistic 12 Degrees of freedom 13 P-value (Chi-square) 62.061 8 0.000 14 15 Parameter Estimates: 16 17 Standard errors Standard 18 Information Expected 19 Information saturated (h1) model Structured 20 21 Latent Variables: 22 12 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.607 0.432 23 f_1 =~ 24 Q1_1 1.000 25 Q1_2 1.548 0.106 14.538 0.000 0.939 0.801 26 Q1_3 1.348 0.082 16.362 0.000 0.818 0.626 27 f_2 =~ 28 Q1_6 0.753 0.611 1.000

7.2 R でやってみる（検証的因子分析） 29 Q1_7 1.220 0.075 16.316 0.000 0.919 0.698 30 Q1_8 0.894 0.053 16.756 0.000 0.673 0.523 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.115 0.016 7.349 0.000 0.251 0.251 31 32 Covariances: 33 34 f_1 ~~ 35 f_2 36 37 Variances: 38 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 39 .Q1_1 1.607 0.051 31.229 0.000 1.607 0.814 40 .Q1_2 0.492 0.054 9.151 0.000 0.492 0.358 41 .Q1_3 1.036 0.050 20.917 0.000 1.036 0.608 42 .Q1_6 0.950 0.043 22.233 0.000 0.950 0.626 43 .Q1_7 0.888 0.055 16.193 0.000 0.888 0.513 44 .Q1_8 1.205 0.044 27.532 0.000 1.205 0.727 45 f_1 0.368 0.040 9.180 0.000 1.000 1.000 46 f_2 0.567 0.047 11.995 0.000 1.000 1.000 引数 standardized = TRUE を与えると，係数の出力の右に Std.lv と Std.all という 2 つの列が追加されます。Std.lv 列は，推定時に引数 std.lv = TRUE を与えた時と同じ「潜 在変数のみ標準化された場合」の推定値です。一方 Std.all 列はその名の通り「観測変数も 含めて全て標準化された場合」の推定値，すなわち標準化解を出しています。 cfa(std.lv = TRUE) のみだと観測変数まで標準化した標準化解が出力されません。一方 summary(standardized = TRUE) のみだと各因子の第 1 項目の因子負荷が推定されなくな るため，仮説検定を行ってくれません。ということで，推定時には std.lv = TRUE をいれ ると同時に，summary() 時に standardized = TRUE を入れるようにすると良いのではな いかと思います。 標準化解の出し方（おすすめ） 1 result_cfa <- cfa(model_cfa, data = dat, std.lv = TRUE) 2 summary(result_cfa, standardized = TRUE) 1 lavaan 0.6.15 ended normally after 27 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters ML NLMINB 13 6 7 Number of observations 2432 8 9 Model Test User Model: 10 11 13 Test statistic 62.061

7.2 R でやってみる（検証的因子分析） 12 Degrees of freedom 13 P-value (Chi-square) 8 0.000 14 15 Parameter Estimates: 16 17 Standard errors Standard 18 Information Expected 19 Information saturated (h1) model Structured 20 21 Latent Variables: 22 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 23 f_1 =~ 24 Q1_1 0.607 0.033 18.359 0.000 0.607 0.432 25 Q1_2 0.939 0.034 27.729 0.000 0.939 0.801 26 Q1_3 0.818 0.034 24.120 0.000 0.818 0.626 27 f_2 =~ 28 Q1_6 0.753 0.031 23.990 0.000 0.753 0.611 29 Q1_7 0.919 0.035 26.112 0.000 0.919 0.698 30 Q1_8 0.673 0.031 21.419 0.000 0.673 0.523 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.251 0.028 8.852 0.000 0.251 0.251 31 32 Covariances: 33 34 f_1 ~~ 35 f_2 36 37 Variances: 38 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 39 .Q1_1 1.607 0.051 31.229 0.000 1.607 0.814 40 .Q1_2 0.492 0.054 9.151 0.000 0.492 0.358 41 .Q1_3 1.036 0.050 20.917 0.000 1.036 0.608 42 .Q1_6 0.950 0.043 22.233 0.000 0.950 0.626 43 .Q1_7 0.888 0.055 16.193 0.000 0.888 0.513 44 .Q1_8 1.205 0.044 27.532 0.000 1.205 0.727 45 f_1 1.000 1.000 1.000 46 f_2 1.000 1.000 1.000 ということで改めて因子負荷です。 21 Latent Variables: 22 14 Estimate Std.Err z-value P(>|z|) Std.lv Std.all Q1_1 0.607 0.033 18.359 0.000 0.607 0.432 Q1_2 0.939 0.034 27.729 0.000 0.939 0.801 26 Q1_3 0.818 0.034 24.120 0.000 0.818 0.626 27 f_2 =~ 28 Q1_6 0.753 0.031 23.990 0.000 0.753 0.611 23 f_1 =~ 24 25

7.2 R でやってみる（検証的因子分析） 29 Q1_7 0.919 0.035 26.112 0.000 0.919 0.698 30 Q1_8 0.673 0.031 21.419 0.000 0.673 0.523 左から「推定値」 「標準誤差」 「検定統計量 𝑧 」 「𝑝 値」 「Std.lv」 「Std.all」と並んでいます。 回帰分析のときと同じように「因子負荷が 0」という帰無仮説に対して仮説検定を行ってい ます。検定の結果は標準化してもしなくても変わらないのでご安心ください。 32 Covariances: 33 34 f_1 ~~ 35 f_2 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.251 0.028 8.852 0.000 0.251 0.251 共分散です。Std.lv と Std.all では潜在変数の分散は 1 になっているため，0.251 が 2 因 子の因子間相関です。 37 Variances: 38 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 39 .Q1_1 1.607 0.051 31.229 0.000 1.607 0.814 40 .Q1_2 0.492 0.054 9.151 0.000 0.492 0.358 41 .Q1_3 1.036 0.050 20.917 0.000 1.036 0.608 42 .Q1_6 0.950 0.043 22.233 0.000 0.950 0.626 43 .Q1_7 0.888 0.055 16.193 0.000 0.888 0.513 44 .Q1_8 1.205 0.044 27.532 0.000 1.205 0.727 45 f_1 1.000 1.000 1.000 46 f_2 1.000 1.000 1.000 各変数の分散です。Std.lv と Std.all では一番下の 2 つの分散は 1 に固定されています。 変数名の左に. がついているものに関しては独自因子の分散（または誤差分散）を表してい ます。 独自因子の分散のなかでも Std.all 列に表示された値は，探索的因子分析と同じように独自 性の値として考えることが出来ます。例えば項目 Q1_1 では，Std.all の値＝独自因子の分 散が 0.814 です。これに対して因子負荷は 0.432 でした。実際に因子負荷の二乗と独自因子 の分散の和は 0.4322 + 0.814 ≃ 1 となります。 ちなみに因子得点は，predict() 関数によって出すことが出来ます。 1 15 head(predict(result_cfa)) 1 f_1 2 [1,] -0.8119757 -1.46119237 3 [2,] -0.2904413 -0.08641555 4 [3,] -0.5899198 -0.03237852 5 [4,] -0.3137981 -0.54036752 f_2

7.3 回帰分析の見方を変える 6 [5,] -1.2359450 -0.24949133 7 [6,] 0.4269696 1.30844315 探索的因子分析のときと同様に，この因子得点の相関はモデル上の因子間相関（今回の場合 0.267）とは異なります。利用する際は気をつけてください。 因子得点の相関 1 cor(predict(result_cfa)) 1 f_1 2 f_1 1.0000000 0.3259455 3 f_2 0.3259455 1.0000000 f_2 7.3 回帰分析の見方を変える SEM では回帰分析も共通の枠組みで取り扱うことが出来ます。つまり観測変数間の共分 散構造をモデルのパラメータで表し，推定を行います。最もシンプルな回帰モデルとして， 図7.5のように協調性の得点 Q1_A を年齢 age で回帰してみましょう。回帰式は (7.10) 𝑥𝐴 = 𝑏0 + 𝑏1 𝑥age + 𝜀𝐴 となります。したがって，この 2 変数の共分散行列は 𝚺=[ 𝜎𝑥2 𝐴 ] 𝜎𝑥2 age 𝜎𝑥𝐴 ,𝑥age =[ 𝑏12 𝜎𝑥2 age + 𝜎𝜀2𝐴 𝑏1 𝜎𝑥2 age 𝜎𝑥2 age ] (7.11) と表せます。回帰モデルなので，変数 𝑥𝐴 は 𝑥age の値に応じて変化していると見ることも 出来ます。言い換えると，𝑥𝐴 の値はモデルの中で決まっているということです。このよう に，モデル式の中で他の変数の関数によって表される変数のことを内生変数 (endogenous variable) と呼びます。一方，変数 𝑥age は他の変数の値によって変化することはありません。 こうした変数を，内生変数に対して外生変数 (exogenous variable) と呼びます。グラフィカ ルモデル的には，他の変数からの一方向矢印が一つでも刺さっているものは全て内生変数と 呼びます。SEMではこの「内生変数」と「外生変数」の分類がちょくちょく出てくるので， 覚えておいてください。 図7.5 単回帰モデル さて，先程の回帰モデルを SEM で見た場合，推定しないといけないパラメータは 𝑏1 と 𝜎𝜀2𝐴 の 2 つです（𝜎𝑥2 age はデータから直接求めたら良いので） 。そしてデータからは（cov() 関数 を使えば）2 つの変数の分散および共分散 𝐒=[ 16 𝑠2𝑥𝐴 𝑠𝑥𝐴 ,𝑥age 𝑠2𝑥age ]=[ 20.667 8.906 120.303 ]

7.3 回帰分析の見方を変える の計 3 つが計算できるため，これでモデルのパラメータが推定できます。一応連立方程式を 立ててみると， 2 2 2 ⎧ 20.667 = 𝑏1 𝜎𝑥age + 𝜎𝜀𝐴 ⎫ { } 20.667 = 120.303𝑏12 + 𝜎𝜀2𝐴 8.906 = 𝑏1 𝜎𝑥2 age = { } ⎨ ⎬ 8.906 = 120.303𝑏1 } {120.303 = 𝜎2 𝑥age ⎭ ⎩ となります。このように，外生変数の分散が連立方程式に含まれる場合には，その外生変数 のデータでの分散をそのまま入れてあげることが出来ます。なお，共分散構造のみを扱うな らば，切片𝑏0 は推定の対象になっていないという点には気をつけてください。切片もモデル の推定対象にしたい場合には，平均構造も含めた SEM を実行する必要があります*6 。 実際に計算してみると，(𝑏1 , 𝜎𝜀2𝐴 ) = (0.074, 20.007) という値が得られます。これを 𝚺 にあ てはめてみると， 𝚺=[ 0.0742 × 120.303 + 20.007 0.074 × 120.303 20.667 ]≃[ 120.303 8.906 120.303 ] ということで，データの共分散行列 𝐒 に一致する値が得られました。 7.3.1 もう少し複雑なモデル（パス解析） モデルが複雑になってもやることは同じです。回帰式を順に展開していき，全ての変数の分 散共分散を「外生変数の分散とパラメータ」で表してあげます。例えば図7.6は重回帰分析に 単回帰分析をくっつけたようなモデルです。ちなみにこのように，単一の回帰式で表せない ような，回帰分析をくっつけたモデルは，SEM の中でも特にパス解析 (path analysis) と呼 ばれることがあります。 図7.6 パス解析モデル 回帰式はそれぞれ { 𝑥𝐴 = 𝑏0𝐴 + 𝑏gender 𝑥gender + 𝑏age 𝑥age + 𝜀𝐴 𝑥𝐸 = 𝑏0𝐸 + 𝑏𝐴 𝑥𝐴 + 𝜀𝐸 (7.12) となりますが，さらに下の式の 𝑥𝐴 を展開することで ⎧𝑥𝐴 = 𝑏gender 𝑥gender + 𝑏age 𝑥age + 𝑏0𝐴 + 𝜀𝐴 { 𝑥 = 𝑏𝐴 (𝑏0𝐴 + 𝑏gender 𝑥gender + 𝑏age 𝑥age + 𝜀𝐴 ) + 𝑏0𝐸 + 𝜀𝐸 ⎨ 𝐸 { = 𝑏𝐴 𝑏0𝐴 + 𝑏𝐴 𝑏gender 𝑥gender + 𝑏𝐴 𝑏age 𝑥age + 𝑏𝐴 𝜀𝐴 + 𝑏0𝐸 + 𝜀𝐸 ⎩ *6 17 (7.13) 平均構造も統一的な枠組みで扱うために，近年では「共分散構造分析」ではなく「構造方程式モデリング」と いう呼び方が一般的なのだと思ってます。

7.3 回帰分析の見方を変える となります。なお，検証的因子分析での独自因子と同じように，パス解析でも内生変数に誤 差項を明示することがあります（図7.7）*7 。こうすることによって，各変数の分散および共 分散を考えるのが多少楽になります（後述）。 図7.7 パス解析モデル（誤差項を追加） 4 つの変数のモデル上の共分散行列 𝚺 の各成分は，頑張って展開すると 𝜎𝑥2 age ⎧ { { 𝜎𝑥2 gender { {𝜎 { 𝑥gender ,𝑥age { { 𝜎𝑥2 𝐴 { { 𝜎𝑥2 𝐸 { = 𝜎𝑥2 age ⎨ 𝜎𝑥 ,𝑥 𝐴 age { { { 𝜎𝑥𝐸 ,𝑥age { { 𝜎 { 𝑥𝐴 ,𝑥gender { { 𝜎𝑥𝐸 ,𝑥gender { { 𝜎𝑥𝐸 ,𝑥𝐴 ⎩ = 𝑏gender 𝜎𝑥gender ,𝑥age + 𝑏age 𝜎𝑥2 age = 𝜎𝑥2 gender = 𝜎𝑥gender ,𝑥age 2 2 𝜎𝑥2 age + 2𝑏gender 𝑏age 𝜎𝑥gender ,𝑥age + 𝜎𝜀2𝐴 = 𝑏gender 𝜎𝑥2 gender + 𝑏age 2 2 2 2 2 2 2 𝑏gender 𝑏age 𝜎𝑥gender ,𝑥age + 𝑏𝐴 𝜎𝜀𝐴 + 𝜎𝜀2𝐸 𝑏age 𝜎𝑥2 age + 2𝑏𝐴 = 𝑏𝐴 𝑏gender 𝜎𝑥2 gender + 𝑏𝐴 = 𝑏𝐴 𝑏gender 𝜎𝑥age ,𝑥gender + 𝑏𝐴 𝑏age 𝜎𝑥2 age = 𝑏gender 𝜎𝑥2 gender + 𝑏age 𝜎𝑥age ,𝑥gender = 𝑏𝐴 𝑏gender 𝜎𝑥2 gender + 𝑏𝐴 𝑏age 𝜎𝑥age ,𝑥gender 2 2 𝜎𝑥2 age + 2𝑏𝐴 𝑏gender 𝑏age 𝜎𝑥gender ,𝑥age + 𝑏𝐴 𝜎𝜀2𝐴 = 𝑏𝐴 𝑏gender 𝜎𝑥2 gender + 𝑏𝐴 𝑏age (7.14) となります。これを図7.7に照らし合わせると，例えば 𝜎𝑥2 𝐴 は 2 • Q1_A → gender → Q1_A (𝑏gender 𝜎𝑥2 gender ) 2 • Q1_A → age → Q1_A (𝑏age 𝜎𝑥2 age ) • Q1_A → gender → age → Q1_A (𝑏gender 𝑏age 𝜎𝑥gender ,𝑥age ) • Q1_A → age → gender → Q1_A (𝑏gender 𝑏age 𝜎𝑥gender ,𝑥age ) • Q1_A → e_A → Q1_A (𝜎𝜀2𝐴 ) の和になっています。検証的因子分析のときと同じように，自分自身を出発してから戻るま での全経路の係数の和になっているわけです。同様に，共分散 𝜎𝑥𝐸 ,𝑥𝐴 では 2 • Q1_E → Q1_A → gender → Q1_A (𝑏𝐴 𝑏gender 𝜎𝑥2 gender ) 2 • Q1_E → Q1_A → age → Q1_A (𝑏𝐴 𝑏age 𝜎𝑥2 age ) • Q1_E → Q1_A → gender → age → Q1_A (𝑏𝐴 𝑏gender 𝑏age 𝜎𝑥gender ,𝑥age ) • Q1_E → Q1_A → age → gender → Q1_A (𝑏𝐴 𝑏gender 𝑏age 𝜎𝑥gender ,𝑥age ) • Q1_E → Q1_A → e_A → Q1_A (𝑏𝐴 𝜎𝜀2𝐴 ) *7 18 パス解析での誤差項は潜在変数としても扱われないため，楕円で囲わないのが一般的です。

7.4 R でやってみる（回帰分析・パス解析） という感じで，Q1_E から Q1_A への全経路での係数の和の形になっていると見ることができ ます*8 。 というわけで，モデル上の共分散行列 𝚺 を求めることが出来ました。あとはこれに対応するデ ータ上の分散（4 つ）と共分散（6 つ）を用いて，5 つのパラメータ (𝑏𝐴 , 𝑏age , 𝑏gender , 𝜎𝜀2𝐴 , 𝜎𝜀2𝐸 ) を推定するだけです。 実際に計算してみると， 𝜎𝑥2 age ⎧ { { 𝜎𝑥2 gender { { {𝜎𝑥gender ,𝑥age { { 𝑏𝐴 ⎨ 𝑏 age { { 𝑏 { gender { { 𝜎𝜀2𝐴 { { 𝜎𝜀2𝐸 ⎩ = 120.254 = 0.221 = 0.204 = 0.550 = 0.071 = 1.891 = 19.210 = 22.154 という解が得られます（上 3 つはデータから直接求めた値） 。これを(7.14)式に当てはめてみ ると， 𝜎𝑥2 age ⎡𝜎 𝑥gender ,𝑥age 𝚺=⎢ ⎢ 𝜎 ⎢ 𝑥𝐴 ,𝑥age ⎣ 𝜎𝑥𝐸 ,𝑥age 𝜎𝑥2 gender 𝜎𝑥𝐴 ,𝑥gender 𝜎𝑥𝐸 ,𝑥gender 𝜎𝑥2 𝐴 𝜎𝑥𝐸 ,𝑥𝐴 120.254 ⎤ ⎡ 0.221 ⎥ = ⎢ 0.204 ⎥ ⎢ 8.902 0.432 ⎥ 2 4.892 0.238 𝜎𝑥𝐸 ⎦ ⎣ 20.685 11.354 ⎤ ⎥ ⎥ 28.395⎦ となります。また，データ上の共分散行列の値は 120.254 ⎡ 0.204 0.221 𝐒=⎢ ⎢ 8.906 0.433 0.260 ⎣ 4.198 20.667 11.359 ⎤ ⎥ ⎥ 28.407⎦ となり，だいぶ近い値が得られました。 7.4 R でやってみる（回帰分析・パス解析） それでは，lavaan を使って回帰分析をやってみましょう。sem() 関数を使うと，回帰分析 を含むパス解析に適した設定で自動的に推定を行ってくれます。 モデル構文では，lm() や glm() のときと同じように，両辺をチルダ（~）でつなげることで 観測変数どうしの回帰を表します。ということで，単回帰分析であれば以下のように 1 文だ けのモデル構文を書けばよいわけです。 *8 「途中で Q1_A を通っているじゃないか」という点は大目に見てください。SEM では全ての内生変数の分散・ 共分散を外生変数の分散の関数で表す必要があるため，パスの経路で考える場合には必ず外生変数を通過す る必要があるのです。 19

7.4 R でやってみる（回帰分析・パス解析）
lavaan で単回帰分析
1

model <- "Q1_A ~ age"

2
3

result_lm <- sem(model, data = dat)

4

summary(result_lm)

1

lavaan 0.6.15 ended normally after 1 iteration

2
3

Estimator

4

Optimization method

5

Number of model parameters

ML
NLMINB
2

6
7

Number of observations

2432

8
9

Model Test User Model:

10
11

Test statistic

12

Degrees of freedom

0.000
0

13
14

Parameter Estimates:

15
16

Standard errors

Standard

17

Information

Expected

18

Information saturated (h1) model

Structured

19
20

Regressions:

21
22

Q1_A ~

23

age

Estimate

Std.Err

z-value

P(>|z|)

0.074

0.008

8.952

0.000

Estimate

Std.Err

z-value

P(>|z|)

19.999

0.574

34.871

0.000

24
25

Variances:

26
27

.Q1_A

基本的な結果の見方は因子分析の時と同じです。出力の Regressions: の部分が回帰係数を，

Variances: の部分が誤差分散を表しています。先程説明したように，SEM では共分散構
造のみを扱うならば，切片𝑏0 は推定の対象にならないため，切片項 𝑏0 の推定値は出力され
ていません。これを出すためには，モデル式に少し手を加える必要があります。これも lm()
や glm() の時と同じなのですが，回帰式に"+ 1" を加えると，切片項を明示的に推定するよ
うに指示することが出来ます。
切片項も出してください

20

1

model <- "

2

Q1_A ~ age + 1

7.4 R でやってみる（回帰分析・パス解析）

3

"

4
5

result_lm <- sem(model, data = dat)

6

summary(result_lm)

1

lavaan 0.6.15 ended normally after 8 iterations

2
3

Estimator

4

Optimization method

5

Number of model parameters

ML
NLMINB
3

6
7

Number of observations

2432

8
9

Model Test User Model:

10
11

Test statistic

12

Degrees of freedom

0.000
0

13
14

Parameter Estimates:

15
16

Standard errors

Standard

17

Information

Expected

18

Information saturated (h1) model

Structured

19
20

Regressions:

21
22

Q1_A ~

23

age

Estimate

Std.Err

z-value

P(>|z|)

0.074

0.008

8.952

0.000

Estimate

Std.Err

z-value

P(>|z|)

21.122

0.253

83.433

0.000

Estimate

Std.Err

z-value

P(>|z|)

19.999

0.574

34.871

0.000

24
25

Intercepts:

26
27

.Q1_A

28
29

Variances:

30
31

.Q1_A

出力を見ると，新たに Intercepts: という欄ができています。これが切片を表している部
分です。
一応，ここまでの結果を，lm() で普通に回帰分析を行った場合と比較してみましょう。

lm() で回帰分析
1

21

summary(lm(Q1_A ~ age, data = dat))

7.4 R でやってみる（回帰分析・パス解析）

1
2

Call:

3

lm(formula = Q1_A ~ age, data = dat)

4
5

Residuals:

6

Min

1Q

Median

3Q

Max

7

-18.565

-2.602

0.621

3.361

8.064

8
9

Coefficients:

10

Estimate Std. Error t value Pr(>|t|)

11

(Intercept) 21.121738

0.253261

83.399

<2e-16 ***

12

age

0.008273

8.949

<2e-16 ***

13

---

14

Signif. codes:

0.074029

0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

15
16

Residual standard error: 4.474 on 2430 degrees of freedom

17

Multiple R-squared:

18

F-statistic: 80.08 on 1 and 2430 DF,

0.0319,

Adjusted R-squared:

0.0315

p-value: < 2.2e-16

lm() では最小二乗法で解を求めている一方，lavaan はデフォルトでは最尤法で解を求めて
います。推定法の違いはありますが，回帰係数および切片はほぼ一致する値が得られました。
続いてパス解析もやってみましょう。基本的には一つ一つの回帰式を順番に書いていけば良
いだけです。簡単ですね。
1

model <- "

2

Q1_A ~ gender + age

3

Q1_E ~ Q1_A

4

"

5
6

result_sem <- sem(model, data = dat)

7

summary(result_sem)

1

lavaan 0.6.15 ended normally after 1 iteration

2
3

Estimator

4

Optimization method

5

Number of model parameters

ML
NLMINB
5

6
7

Number of observations

2432

8
9

Model Test User Model:

10

22

11

Test statistic

12

Degrees of freedom

0.709
2

7.5 回帰分析と因子分析を組み合わせる 13 P-value (Chi-square) 0.702 14 15 Parameter Estimates: 16 17 Standard errors 18 Information 19 Information saturated (h1) model Standard Expected Structured 20 21 Regressions: 22 23 Estimate Std.Err z-value P(>|z|) Q1_A ~ 24 gender 1.891 0.189 9.996 0.000 25 age 0.071 0.008 8.733 0.000 26 Q1_E ~ 27 Q1_A 0.550 0.021 26.173 0.000 Estimate Std.Err z-value P(>|z|) 28 29 Variances: 30 31 .Q1_A 19.210 0.551 34.871 0.000 32 .Q1_E 22.154 0.635 34.871 0.000 7.5 回帰分析と因子分析を組み合わせる 図7.6に基づくパス解析では，因子 A と因子 E に関する変数として和得点（Q1_A と Q1_E） を使用しました。ですが和得点には独自因子（構成概念からすると誤差）の影響が残ってい るため，可能であれば因子得点を用いた方が良い場合があります。あるいは（構成概念）妥 当性の考え方の根底にあった「観測得点はあくまでも構成概念の顕在化」という考えに従え ば，仮説モデルでは構成概念そのもの同士の関係（法則定立ネットワーク）を見るべきだと 言えるでしょう。いま，手元には構成概念（因子）を構成する各項目の回答そのものがあり ます。そこで，図7.3の因子分析と図7.6の回帰分析を組み合わせて，図7.8のようなモデルを 考えてみましょう。 図7.8 回帰分析と因子分析を組み合わせてみた 23

7.5 回帰分析と因子分析を組み合わせる モデルが複雑になってもやることは同じです。全ての（標準化された）観測変数間の相関行列 1 ⎡𝑟 ⎢ gender,age ⎢ 𝑟𝑦1,age ⎢ 𝑟 𝚺 = ⎢ 𝑦2,age 𝑟 ⎢ 𝑦3,age ⎢ 𝑟𝑦6,age ⎢ 𝑟𝑦7,age ⎣ 𝑟𝑦8,age 1 𝑟𝑦1,gender 𝑟𝑦2,gender 𝑟𝑦3,gender 𝑟𝑦6,gender 𝑟𝑦7,gender 𝑟𝑦8,gender 1 𝑟𝑦2,𝑦1 𝑟𝑦3,𝑦1 𝑟𝑦6,𝑦1 𝑟𝑦7,𝑦1 𝑟𝑦8,𝑦1 1 𝑟𝑦3,𝑦2 𝑟𝑦6,𝑦2 𝑟𝑦7,𝑦2 𝑟𝑦8,𝑦2 1 𝑟𝑦6,𝑦3 𝑟𝑦7,𝑦3 𝑟𝑦8,𝑦3 1 𝑟𝑦7,𝑦6 𝑟𝑦8,𝑦6 1 𝑟𝑦8,𝑦7 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 1⎦ と，データ上の相関行列 1 ⎡0.039 ⎢ ⎢0.156 ⎢0.113 𝐒=⎢ 0.064 ⎢ ⎢0.095 ⎢0.018 ⎣0.065 1 0.157 0.184 0.137 0.008 0.060 0.048 1 0.353 0.274 −0.011 −0.011 0.024 1 0.497 0.098 0.124 0.185 1 0.110 0.142 0.126 1 0.435 0.313 1 0.359 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 1⎦ がなるべく近くなるように各種パラメータを推定してあげます。 ここでいくつかの変数間のモデル上の相関係数を具体的に見てみましょう。同じ因子に属す る Q1_1 と Q1_2 の相関は， （検証的）因子分析のときには 2 項目の因子負荷の積 (𝑏11 𝑏21 ) と いう簡単な形で表すことができていました。今回のモデルではどうでしょうか。因子得点は 内生変数として扱われており，𝑓𝐴 = 𝑏age 𝑥age + 𝑏gender 𝑥gender + 𝜀𝐴 という回帰式になって います。これまでと同じように Q1_1 と Q1_2 をそれぞれ外生変数（age と gender）の関数 まで展開してあげるならば， 𝑦1 = 𝑓𝐴 𝑏1 + 𝜀1 = (𝑏age 𝑥age + 𝑏gender 𝑥gender + 𝜀𝐴 )𝑏1 + 𝜀1 𝑦2 = 𝑓𝐴 𝑏2 + 𝜀2 (7.15) = (𝑏age 𝑥age + 𝑏gender 𝑥gender + 𝜀𝐴 )𝑏2 + 𝜀2 となり，これを頑張って展開して相関を求める必要がありそうです。しかしよく考えてみる と，因子分析モデルでは潜在変数の分散を 1 に固定する，という制約がありました。つまり 𝑓𝐴 = (𝑏age 𝑥age + 𝑏gender 𝑥gender + 𝜀𝐴 ) の部分は，これをすべて合わせて分散 1 になるよう に調整される，ということです。したがってこの 2 項目の相関を求めると，結局前回資料の 6.3.1 項の議論に帰着して*9 𝑟𝑦1,𝑦2 = 𝑟(𝑓𝐴 𝑏1 +𝜀1 ),(𝑓𝐴 𝑏2 +𝜀2 ) = 𝑟(𝑓𝐴 𝑏1 ),(𝑓𝐴 𝑏2 ) = 𝑏1 𝑏2 (𝜎𝑓𝐴 ) (7.16) = 𝑏 1 𝑏2 となるのです。また，誤差分散 𝜎𝜀2𝐴 の値は 1 から回帰係数の二乗和をひいた値になっていま す。このあたりの議論は，矢印の向きこそ反対ですが因子分析の時に見た「観測変数の分散 は因子負荷の二乗和と独自因子の分散の和になっている」というのと同じですね。 *9 24 共通因子と独自因子の間は無相関であり，またモデル上独自因子間も無相関と仮定しているため，これらの 項が全て消えています。

7.5 回帰分析と因子分析を組み合わせる 続いて最も距離の遠い age と Q1_8 の相関を見てみます。Q1_8 を外生変数（age と gender） の関数に戻してあげると， 𝑦8 = 𝑓𝐸 𝑏8 + 𝜀8 = (𝑏𝐴 𝑓𝐴 + 𝜀𝐸 )𝑏8 + 𝜀8 = (𝑏𝐴 [𝑏age 𝑥age + 𝑏gender 𝑥gender + 𝜀𝐴 ] + 𝜀𝐸 ) 𝑏8 + 𝜀8 (7.17) error ⏞⏞⏞ = 𝑏𝐴 𝑏8 (𝑏age 𝑥age + 𝑏gender 𝑥gender ) + ⏞⏞ 𝑏𝐴 𝑏8⏞ 𝜀⏞⏞⏞ 𝐴 + 𝑏 8 𝜀𝐸 + 𝜀 8 となります。後ろの 𝜀 に関する項 (error) は，観測変数 age との相関はゼロなはずなので無 視して分散を求めると， 2 𝑟age,𝑦8 = 𝑏𝐴 𝑏8 𝑏age 𝜎age + 𝑏𝐴 𝑏8 𝑏gender 𝜎age,gender = 𝑏𝐴 𝑏8 𝑏age + 𝑏𝐴 𝑏8 𝑏gender 𝜎age,gender (7.18) となります。ここでもやはり • Q1_8 → f_E → f_A → age (𝑏𝐴 𝑏8 𝑏age ) • Q1_8 → f_E → f_A → gender → age (𝑏𝐴 𝑏8 𝑏gender 𝜎age,gender ) という要領で，2 変数間の全経路の係数の積の和の形になっていることがわかります。 モデルが複雑になっても，内部でやっていることは今までと同じであることを確認したので， 実際に推定を行ってみましょう。 因子分析と回帰分析を同時にやる 1 model <- " 2 f_A =~ Q1_1 + Q1_2 + Q1_3 3 f_E =~ Q1_6 + Q1_7 + Q1_8 4 f_A ~ age + gender 5 f_E ~ f_A 6 age ~~ gender 7 " 6 行目（age ~~ gender）は，age と gender の共分散を明示的に推定するように指示して います。書かなくても勝手にデータから計算した値を使ってくれるので，この行の有無でパ ラメータの推定値が変わることはないのですが，この行を書くことによって，結果の出力に データの共分散の値を表示してくれるようになります。 1 result_sem <- sem(model, data = dat, std.lv = TRUE) 2 summary(result_sem, standardized = TRUE) 1 lavaan 0.6.15 ended normally after 36 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters 6 25 ML NLMINB 18

7.5 回帰分析と因子分析を組み合わせる 7 Number of observations 2432 8 9 Model Test User Model: 10 11 Test statistic 12 Degrees of freedom 13 P-value (Chi-square) 129.227 18 0.000 14 15 Parameter Estimates: 16 17 Standard errors 18 Information 19 Information saturated (h1) model Standard Expected Structured 20 21 Latent Variables: 22 Estimate Std.Err z-value P(>|z|) Std.lv Std.all Q1_1 0.602 0.031 19.255 0.000 0.628 0.447 25 Q1_2 0.888 0.030 29.426 0.000 0.926 0.790 26 Q1_3 0.786 0.030 25.796 0.000 0.819 0.627 27 f_E =~ 28 Q1_6 0.728 0.031 23.767 0.000 0.753 0.611 29 Q1_7 0.888 0.035 25.607 0.000 0.918 0.698 30 Q1_8 0.651 0.031 21.262 0.000 0.673 0.523 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 23 f_A =~ 24 31 32 Regressions: 33 34 f_A ~ 35 age 0.014 0.002 6.275 0.000 0.013 0.147 36 gender 0.522 0.053 9.894 0.000 0.501 0.236 0.252 0.030 8.402 0.000 0.254 0.254 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.204 0.105 1.945 0.052 0.204 0.039 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 37 f_E ~ 38 f_A 39 40 Covariances: 41 42 43 age ~~ gender 44 45 Variances: 46 26 47 .Q1_1 1.581 0.051 31.144 0.000 1.581 0.801 48 .Q1_2 0.517 0.046 11.273 0.000 0.517 0.376 49 .Q1_3 1.034 0.046 22.666 0.000 1.034 0.607 50 .Q1_6 0.950 0.043 22.240 0.000 0.950 0.626 51 .Q1_7 0.889 0.055 16.228 0.000 0.889 0.513 52 .Q1_8 1.205 0.044 27.533 0.000 1.205 0.727

7.5 回帰分析と因子分析を組み合わせる 53 54 age 120.254 3.449 34.871 0.000 120.254 0.221 0.006 34.871 0.000 gender 1.000 0.221 1.000 55 .f_A 1.000 0.920 0.920 56 .f_E 1.000 0.935 0.935 結果の見方も同じなので省略します。一応全ての標準化解を図7.8に書き加えると，図7.9の ようになるわけです。実際に論文などで報告する際には，全ての値を載せるとごちゃごちゃ してしまう可能性があります。そこで，以下のような省略を行い簡潔化することがあります。 このあたりは自由というか，どういう基準で簡略化したかがハッキリわかればOKです。 • 因子分析の部分は書かない（因子だけ残す） • 係数が一定以下のところには矢印をひかない • 統計的に有意ではない係数のところは矢印をひかない 図7.9 結果 semPlot パッケージにある semPaths() という関数は，結果のオブジェクトから自動的にパ ス図を作成してくれます。残念ながらそのまま論文に載せられるほど高性能ではないのです が，自分の書いたモデル式に間違いが無いかを確認したり，係数が一定以下のところを除外 したときの図をさっと確認するのには使えると思います。 0.80 0.38 0.61 1.00 1.00 Q1_1 Q1_2 Q1_3 age gnd 0.45 0.79 0.63 0.15 0.24 f_A 0.94 0.25 0.92 f_E 0.61 0.52 Q1_6 Q1_7 Q1_8 0.63 0.51 0.73 図7.10 27 0.70 semPaths()

7.5 回帰分析と因子分析を組み合わせる 7.5.1 （おまけ）SEM モデルを一つの式で表す ここまでに見てきたように，SEM は回帰分析と因子分析を内包する上位の概念です。回帰分 析も因子分析も，「被説明変数を説明変数の線形和で表現する」という点は変わらないので， これらがくっついたモデル（図7.8）であっても，これまでと同じように線形代数を使えば一 つの式で表すことができます (Adachi, 2020)。 まず，因子分析の部分（潜在変数から観測変数に矢印が引かれている部分）だけ見てみると， これは単純な 2 因子モデルなので， (1) 𝐲𝑝 = 𝐟𝑝 𝐁(1) + 𝛆𝑝 (7.19) と表すことができます。上付きの添字 (1) は，因子分析パートと回帰分析パートの 𝐁 や 𝛆 を 区別するために便宜的につけているものです。それぞれを展開すると， [𝑦𝑝1 𝑦𝑝2 𝑦𝑝3 𝑦𝑝6 𝑦𝑝7 𝑦𝑝8 ] = [𝑓𝑝A 𝑓𝑝E ] [ 𝜀𝑝2 + [𝜀𝑝1 𝑏11 0 𝜀𝑝3 𝑏21 0 𝑏31 0 𝜀𝑝4 0 𝑏42 𝜀𝑝5 0 𝑏52 0 ] 𝑏62 (7.20) 𝜀𝑝6 ] となります。同様に，回帰分析の部分は， (2) (7.21) 𝑓𝑝E = 𝑓𝑝A 𝑏A + 𝛆𝑝 となります。ここで，因子分析の表現に合わせるため，式の左辺にすべての潜在変数が登場 するように少し式を書き換えます。具体的には { 𝑓𝑝A = 𝑓𝑝A (7.22) 𝑓𝑝E = 𝑓𝑝A 𝑏A +𝜀𝑝A という連立方程式を用意し，これを行列によって [𝑓𝑝A 𝑓𝑝E ] = [𝑓𝑝A 𝑓𝑝E ] [ 0 𝑏A ] + [𝑓𝑝A 0 0 𝜀𝑝A ] (2) ⟶ 𝐟𝑝 = 𝐟𝑝 𝐁(2) + 𝛆𝑝 (7.23) (7.24) と表現します。 これで準備は完了です。あとは(7.20)式と(7.24)式をうまく組み合わせると， [𝐟𝑝 𝐲𝑝 ] = [𝐟𝑝 𝐲𝑝 ] [ 𝐁(1) 𝟎 𝐁(2) ] + [𝛆(1) 𝑝 𝟎 ⟶ 𝐯 𝑝 = 𝐯 𝑝 𝐁 + 𝛆𝑝 (2) 𝛆𝑝 ] (7.25) (7.26) となります。回帰係数行列の 𝟎 は，今回のモデルでは 𝑦 から 𝑓 への回帰や 𝑦 同士の回帰は一 つも無いため，ゼロ行列を置いています。もしも 𝑦 同士の回帰などがあれば，それに応じて 係数を入れる必要があります。とにかく，これによってすべての（観測・潜在）変数同士の 回帰分析のモデルを統一的に表すことができました。 28

7.6 モデルの適合度 ちなみに展開すると横長になりすぎてかなり見にくいですが，一応 [𝑓𝑝A [𝑓𝑝A 𝑓𝑝E 𝑓𝑝E + [𝑓𝑝A 𝑦𝑝1 𝑦𝑝1 𝜀𝑝A 𝑦𝑝2 𝑦𝑝2 𝜀𝑝1 𝜀𝑝2 𝑦𝑝3 𝑦𝑝3 𝑦𝑝6 𝑦𝑝6 𝜀𝑝3 𝜀𝑝4 𝑦𝑝7 𝑦𝑝8 ] = 𝑦𝑝7 0 𝑏A ⎡0 0 ⎢ ⎢0 0 ⎢0 0 𝑦𝑝8 ] ⎢ 0 0 ⎢ ⎢0 0 ⎢0 0 ⎣0 0 𝜀𝑝5 𝑏11 0 0 0 0 0 0 0 𝑏21 0 0 0 0 0 0 0 𝑏31 0 0 0 0 0 0 0 0 𝑏42 0 0 0 0 0 0 0 𝑏52 0 0 0 0 0 0 0 𝑏62 ⎤ ⎥ 0⎥ 0⎥ 0⎥ ⎥ 0⎥ 0⎥ 0⎦ 𝜀𝑝6 ] (7.27) という形をしています。 あとは(7.26)式を，左辺が 𝐲 のみになるように変形していくだけです。(𝑇 + 𝐼) 単位行列 𝐈(𝑇 +𝐼) を用いると 𝐯𝑝 𝐈(𝑇 +𝐼) = 𝐯𝑝 𝐁 + 𝛆𝑝 (7.28) と表せるので，𝐯𝑝 𝐁 を移項させてから両辺に逆行列をかけると 𝐯𝑝 = 𝛆𝑝 (𝐈(𝑇 +𝐼) − 𝐁)−1 (7.29) と変形させることができます。ここで， 𝐯𝑝 = [𝐟𝑝 𝐲𝑝 ] (7.30) であることを踏まえると，𝐯𝑝 のうち 𝐟𝑝 に当たる部分が 0 になったものが 𝐲𝑝 なので，(𝑇 × 𝑇 ) のゼロ行列 𝟎(𝑇 ) と，(𝐼 × 𝐼) の単位行列 𝐈(𝐼) を縦に並べた行列 𝐇 を使って， 𝐲𝑝 = [𝐟𝑝 𝐲𝑝 ] [ 𝟎(𝑇 ) ] 𝐈(𝐼) (7.31) = 𝐯𝑝 𝐇 と表せます。あとは，これに(7.29)式を代入することで， 𝐲𝑝 = 𝛆𝑝 (𝐈(𝑇 +𝐼) − 𝐁)−1 𝐇 (7.32) という形で，どんな SEM モデルであっても，すべての観測変数をすべてのパラメータの関 数で表す一般化が完成しました。 7.6 モデルの適合度 ここまで，いくつかの簡単なモデルで SEM の基本的な仕組みを紹介しました。SEM の目的 はなるべく少ないモデルパラメータによってデータの共分散行列に近い値を再現することで す。データ（の共分散行列）に対するモデル（の共分散行列）の当てはまりの程度を表す指 標を（モデルの）適合度と呼び，その指標には様々な種類があります。以下でいくつかの代 表的な指標を紹介します。どれか一つ・一種類の指標を使えば OK というものではありませ ん。様々な視点からの評価を行い，総合的にモデルの適合度を判定するようにしましょう。 適合度指標に関する議論は，星野他 (2005) あたりが詳しいです。 29

7.6 モデルの適合度 7.6.1 相対的指標 相対的指標では， 「最悪のモデル」と比べてどの程度当てはまりが改善したかを考えます。こ こでの「最悪のモデル」とは，全ての変数が完全に無相関だと仮定したモデルです。例えば 図7.7のパス解析で考えると，図7.11のような状態だということです。 図7.11 なにもないモデル 何もないということは，モデル上の相関行列 𝚺0 は 𝜎𝑥2 ⎡ 0age 𝚺0 = ⎢ ⎢ 0 ⎣ 0 𝜎𝑥2 gender 0 0 𝜎𝑥2 𝐴 0 ⎤ ⎥ ⎥ 𝜎𝑥2 𝐸 ⎦ (7.33) ということです。もちろんこれ以上当てはまりの悪いモデルは考えられません。ということ で，相対的指標ではこの独立モデル（ヌルモデル）と比べて，設定したモデルではどの程度 当てはまりが改善したかを確認します*10 。ここで用いる「当てはまり」の指標は，最尤推 定の場合は尤度を元に設定します。分析者が設定した仮説モデルに関する当てはまりの指標 𝐿ℎ は 𝐿ℎ = (𝑃 − 1) {tr(𝚺−1 𝐒) − log |𝚺−1 𝐒| − 𝐼} (7.34) となります。この式は，最尤推定量におけるモデルの尤度（波括弧の中身）を 𝑃 − 1 倍した ものなのですが，細かい式の意味や導出はともかく，この値が「データにおける観測された 共分散行列𝐒」と「仮説モデルの共分散行列𝚺」の比を調整した値だということだけ理解し ておいてください。つまりこの値は 𝚺 と 𝐒 の値が近くなるほど小さな値を取る，いわば乖 離度を表しています*11 。そして同様に，独立モデルでの当てはまりの指標 𝐿0 を −1 𝐿0 = (𝑃 − 1) {tr(𝚺−1 0 𝐒) − log |𝚺0 𝐒| − 𝐼} (7.35) とします。相対的指標では，この 𝐿ℎ と 𝐿0 の比較を行います。 相対的指標でよく用いられている指標の一つは，TLI (Tucker-Lewis index; Tucker and Lewis, 1973) でしょう。 TLI = *10 𝐿0 𝐿ℎ 𝑑𝑓0 − 𝑑𝑓ℎ 𝐿0 𝑑𝑓0 − 1 (7.36) 独立モデルとは反対に，全ての観測変数間に両方向の矢印を引いたモデル（共分散を独立に推定するモデル） のことを飽和モデル（フルモデル）と呼びます。わざわざ飽和モデルを SEM で実行する必要はまったくな いのですが「全てのモデルは独立モデルと飽和モデルの間に位置している」という感覚は持っておくと良い かもしれません。 *11 もう少し詳しく説明すると，𝐿 は「対立仮説を飽和モデルとしたときの 𝚺 の尤度比検定統計量」になって ℎ います。尤度比検定では，𝐿ℎ つまり乖離度が大きい場合には対立仮説が採択されるため，「いま考えている モデル 𝚺 よりは飽和モデルが真のモデル構造と考えられる」となります。 30

7.6 モデルの適合度 𝑑𝑓ℎ と 𝑑𝑓0 はそれぞれ仮説モデルと独立モデルでの自由度です。自由度はデータの分散共分 散の総数（ 𝐼(𝐼+1) ）からパラメータ数を引いた数になります。例えば図7.7のモデルでは，分 2 散共分散の総数 10 個を 8 つのパラメータで表現しました。そのため自由度は 𝑑𝑓ℎ = 2 とな ります。本来データの分散共分散は 10 個あったのに，8 個だけでそれを再現したという意 味で，自由度は「何個のパラメータを削減しているか」という見方もできそうです。あらた 𝐿 めて TLI の式を見てみます。 𝑑𝑓0 は独立モデルでの乖離度です。削減したパラメータ 1 個に 0 𝐿 つきどの程度のズレが発生したかを表しています*12 。これに対して， 𝑑𝑓ℎ は仮説モデルでの ℎ 乖離度です。あってもなくても良いパラメータは削減しても大きなズレは発生しないですが， 重要なパラメータを削減すると大きなズレが発生します。ズレの大きさをパラメータ数で割 ることで，「追加したパラメータ数の割に当てはまりが改善したか」を評価しているわけで すね。 似たような指標に CFI (Comparative fit index; Bentler, 1990) があり，これもよく用いら れています。 CFI = 1 − max(𝐿ℎ − 𝑑𝑓ℎ , 0) max(𝐿ℎ − 𝑑𝑓ℎ , 𝐿0 − 𝑑𝑓0 , 0) (7.37) CFI では乖離度 𝐿ℎ から自由度 𝑑𝑓ℎ を引いています。が，あとは大体同じです。分母の方の max 関数では，普通は最も乖離度の大きい 𝐿0 − 𝑑𝑓0 が選ばれるはずです。そうなると CFI = 1 − 𝐿ℎ − 𝑑𝑓ℎ 𝐿0 − 𝑑𝑓0 (7.38) となり，仮説モデルにおける当てはまりの改善度を表していることがよくわかります。 TLI や CFI は 1 に近いほど良い指標です。一般的には 0.9 くらいはほしいだの 0.95 あると 嬉しいだの言われています。しかし 𝐿ℎ は 𝚺 と 𝐒 のズレの大きさに規定されている以上，基 本的にはパラメータ数を増やすほど値が大きくなりやすいという点には注意が必要です。 7.6.2 絶対的指標 データの相関行列 𝐒 とモデルの相関行列 𝚺 の間に 𝐒=𝚺+𝐄 (7.39) という関係があると考えれば，𝐄 の部分がどれだけ小さいか，という指標を考えることがで きます。例えば図7.7のパス解析では， 𝐒=𝚺+𝐄 120.254 ⎡ 0.204 0.221 =⎢ ⎢ 8.902 0.432 0.238 ⎣ 4.892 120.254 ⎡ 0.204 0.221 =⎢ ⎢ 8.906 0.433 0.260 ⎣ 4.198 *12 31 20.685 11.354 0 ⎤ ⎡ 0 0 ⎥+⎢ ⎥ ⎢ 0.004 0.001 −0.018 28.395⎦ ⎣−0.694 0.022 0.005 20.667 11.359 ⎤ ⎥ ⎥ 28.407⎦ ⎤ ⎥ ⎥ 0.012⎦ ちなみに 𝐿ℎ や 𝐿0 を 𝑃 − 1 で割ったものの期待値は，対応する自由度に一致することが知られています。 つまり，自由度が小さいほど＝パラメータ数が多いほどズレの期待値は当然のように小さくなるわけです。

7.6 モデルの適合度 という関係になっています。ここから適合度を単純に考えると，𝐄 の平均のようなものを 計算すると良さそうです。この指標を RMR (root mean-squared residual) と呼び，式と しては √ √ RMR = √ ⎷ 𝐼 𝐼 2 ∑ ∑(𝜎𝑖𝑗 − 𝑠𝑖𝑗 )2 𝐼(𝐼 − 1) 𝑖=1 𝑗=1 (7.40) と表されます。ここで 𝜎𝑖𝑗 , 𝑠𝑖𝑗 はそれぞれ 𝚺, 𝐒 の (𝑖, 𝑗) 成分の値を意味しています。しかし， RMR は変数のスケールの影響を受けてしまいます。例えば全ての変数の値を 10 倍（＝分散 は 100 倍）にしたとしたら，𝐒 と 𝚺 の各要素の値は単純に 100 倍されてしまうため，𝐄 の 平均である RMR も 100 倍になってしまいます。 こ れ で は 基 準 を 決 め ら れ な い の で，RMR を 標 準 化 す る こ と で 生 ま れ た の が SRMR(standardized RMR; Hu and Bentler, 1999) で す。 標 準 化 の 仕 方 は， 「共 分 散から相関係数を導出する時」と同様にすることで， √ √ SRMR = √ ⎷ 𝐼 𝐼 𝜎𝑖𝑗 − 𝑠𝑖𝑗 2 ∑∑( ) 𝐼(𝐼 − 1) 𝑖=1 𝑗=1 𝑠𝑖𝑖 𝑠𝑗𝑗 2 (7.41) と定式化されます。 SRMR の他にもよく使われる指標としては，GFI(Goodness of fit index; Bentler, 1983) が あります。これも基本的な考え方は同じで，𝐒 と 𝚺 の各要素の値が近いほど良くなります。 𝑖 𝐼 GFI = 1 − ∑𝑖=1 ∑𝑗=1 (𝑠𝑖𝑗 − 𝜎𝑖𝑗 ) 𝐼 𝑖 ∑𝑖=1 ∑𝑗=1 𝑠𝑖𝑗 (7.42) SRMR や GFI は，回帰分析の決定係数や因子分析の分散説明率のようなもので，パラメー タの数が多いほど値は改善していきます。因子分析の分散説明率では「最低限これくらいは 欲しい」といった考え方を紹介しました。その観点ではこれらの指標は役に立ちます。例え ば SRMR では 0.08 や 0.05 以下くらいだと良いと言われていたり，GFI や後述する AGFI は 0.95 以上だと良いと言われていたりします (Hancock et al., 2019; 豊田, 2014)。ですが SEM の目的は「なるべく少ないパラメータで」達成したいところです。そう考えると，これ らの指標だけを用いていると少し問題があるわけです。 ちなみに，GFI には自由度で調整した AGFI(Adjusted GFI) というものがあります。 AGFI = 1 − 𝐼(𝐼 + 1) 1 (1 − GFI) 2 𝑑𝑓 (7.43) 推定するパラメータ数が多くなるほど相関行列の乖離度 (1 − GFI) を過大評価するような調 整がされているわけです。 7.6.3 倹約的指標 データへの当てはまりは良いほど嬉しいのですが，「なるべく少ないパラメータで」を満た すために，先ほど紹介した AGFI のように推定するパラメータ数を考慮した適合度指標を考 32

7.6 モデルの適合度 えます。このクラスで多分最もよく用いられているのが RMSEA (Root Mean Square error of approximation; Steiger and Lind, 1980) です。 RMSEA = √max ( 1 𝐿ℎ − , 0) 𝑑𝑓ℎ 𝑃 −1 (7.44) TLI のところで見た「削減したパラメータ 1 個あたりの乖離度の大きさ」の指標です。同じ 乖離度であれば，パラメータ数が少ない＝自由度が大きいほど良い値になります。RMSEA に関しては，0.05 を下回っていると良いとか，90% 信頼区間が 0.1 を含まないと良いなどと 言われています。 倹約的指標には，他にも PGFI (parsimonious GFI; Mulaik et al., 1989) といったものもあ ります。 PGFI = 𝑑𝑓ℎ 𝐼(𝐼−1) 2 GFI (7.45) ただ，PGFI はシンプルなモデルでは変動が大きすぎることがあります。例えば図7.7のモデル では， 10 個の分散共分散を 8 個のパラメータで説明していますが，この場合 𝑃 𝐺𝐹 𝐼 = 2 10 𝐺𝐹 𝐼 となり，GFI の 5 分の 1 になってしまいます。独立モデルだったとしても各変数の分散パラ メータは必ず存在しているため，𝑃 𝐺𝐹 𝐼 = 6 10 𝐺𝐹 𝐼 となります。ということで，PGFI と同 様の補正がかかる倹約的指標では，どの程度の値なら十分かを明確に決めづらいためあまり 使われることはありません。 7.6.4 モデル比較の指標 適合度指標とは少し異なりますが，パラメータ数の割に当てはまりの良いモデルを選ぶため の指標として情報量規準というものがあります。代表的な情報量規準には AIC = 𝐿ℎ − 2𝑑𝑓ℎ (7.46) BIC = 𝐿ℎ − log(𝑃 )𝑑𝑓ℎ (7.47) などがあります (Akaike, 1998; Schwarz, 1978)。式からわかるように，パラメータ数が多く なるほど第 2 項の値が小さくなり，結果的に情報量規準の値が大きくなってしまいます。 情報量規準には絶対的な基準はありません。単体のモデルに対する値を見て一喜一憂するも のではなく，複数のモデルにおける値を比較してより小さい方を選ぶためのものです。 一般的に，乖離度 𝐿ℎ はサンプルサイズが大きくなるほど大きな値になりやすいものです。 AIC ではパラメータを 1 つ追加した際のペナルティがサンプルサイズによらず 2 のため，サ ンプルサイズが大きくなるほど「パラメータを追加したことによる当てはまりの改善」のほ うが「パラメータを追加したことによるペナルティ」よりも大きくなりやすく，結果的にパ ラメータ数の多いモデルを選びやすくなる，という性質があったりします。また 2 つのモデ ルでの値の差がわずかな場合には，標本誤差によって大小関係が変わりうるため，単一の情 報量規準だけを用いてモデルを選択するのは危険なときがあるかもしれません。少しだけ注 意しましょう。 33

7.6 モデルの適合度 7.6.5 確認してみよう それでは，これまでに紹介したものを含めた様々な適合度指標を lavaan で出してみましょ う。fitmeasures() 関数を使うと簡単に出すことが出来ます*13 。 適合度指標を出す 1 fitmeasures(result_sem) 1 npar fmin chisq 2 18.000 0.027 129.227 3 df pvalue baseline.chisq 4 18.000 0.000 2331.678 5 baseline.df baseline.pvalue cfi 6 28.000 0.000 0.952 7 tli nnfi rfi 8 0.925 0.925 0.914 9 nfi pnfi ifi 10 0.945 0.607 0.952 11 rni logl unrestricted.logl 12 0.952 -34146.289 -34081.676 13 aic bic ntotal 14 68328.579 68432.915 2432.000 15 bic2 rmsea rmsea.ci.lower 16 68375.725 0.050 0.042 17 rmsea.ci.upper rmsea.ci.level rmsea.pvalue 18 0.059 0.900 0.450 rmsea.close.h0 rmsea.notclose.pvalue rmsea.notclose.h0 19 20 0.050 0.000 0.080 21 rmr rmr_nomean srmr 22 0.308 0.308 0.036 23 srmr_bentler srmr_bentler_nomean crmr 24 0.036 0.036 0.041 25 crmr_nomean srmr_mplus srmr_mplus_nomean 26 0.041 0.036 0.036 27 cn_05 cn_01 gfi 28 544.308 656.021 0.987 29 agfi pgfi mfi 30 0.974 0.494 0.977 31 ecvi 32 0.068 これまでに紹介したもの以外にも色々出ています。紹介したものに関する指標だけ説明 すると *13 summary() に引数 fit = TRUE を与えても良いのですが，この方法だと一部の適合度指標しか出してくれま せん。 34

7.6 モデルの適合度 指標名 解説 chisq 𝜒2 統計量（𝐿ℎ ） df 自由度（𝑑𝑓ℎ ） 𝜒2 検定の 𝑝 値 pvalue 独立モデルでの値（𝐿0 , 𝑑𝑓0 ） baseline.*** cfi CFI aic AIC bic BIC rmsea RMSEA rmsea.ci.lower RMSEA の 90% 信頼区間の下限 rmsea.ci.upper RMSEA の 90% 信頼区間の上限 rmsea.pvalue 「RMSEA が 0.05 以下」という帰無仮説に対する 𝑝 値 srmr SRMR gfi GFI agfi AGFI pgfi PGFI といった感じです。 実際に図7.3の「本来の 2 因子モデル」と，これに「余計な矢印が入った 2 因子モデル」の 2 つを比較してみましょう。 2 つのモデルを比較 1 model1 <- " 2 # 本来の2因子モデル 3 f_1 =~ Q1_1 + Q1_2 + Q1_3 4 f_2 =~ Q1_6 + Q1_7 + Q1_8 5 " 6 model2 <- " 7 # 余計なものが入っている 8 f_1 =~ Q1_1 + Q1_2 + Q1_3 + Q1_6 9 f_2 =~ Q1_6 + Q1_7 + Q1_8 10 " 11 12 result1 <- cfa(model1, data = dat) 13 result2 <- cfa(model2, data = dat) 14 15 # 見やすいように横に並べて丸めました 16 round(data.frame(fitmeasures(result1), fitmeasures(result2)), 4) 1 2 35 fitmeasures.result1. fitmeasures.result2. npar 13.0000 14.0000

7.6 モデルの適合度 3 fmin 0.0128 0.0118 4 chisq 62.0615 57.1972 5 df 8.0000 7.0000 6 pvalue 7 baseline.chisq 8 baseline.df 9 baseline.pvalue 10 11 0.0000 0.0000 2117.1086 2117.1086 15.0000 15.0000 0.0000 0.0000 cfi 0.9743 0.9761 tli 0.9518 0.9488 12 nnfi 0.9518 0.9488 13 rfi 0.9450 0.9421 14 nfi 0.9707 0.9730 15 pnfi 0.5177 0.4541 16 ifi 0.9744 0.9762 17 rni 0.9743 0.9761 18 logl -23329.6593 -23327.2272 19 unrestricted.logl -23298.6286 -23298.6286 20 aic 46685.3187 46682.4544 21 bic 46760.6728 46763.6050 22 ntotal 2432.0000 2432.0000 23 bic2 46719.3688 46719.1238 24 rmsea 0.0527 0.0543 25 rmsea.ci.lower 0.0410 0.0418 26 rmsea.ci.upper 0.0653 0.0678 27 rmsea.ci.level 0.9000 0.9000 28 rmsea.pvalue 0.3334 0.2706 29 rmsea.close.h0 0.0500 0.0500 30 rmsea.notclose.pvalue 0.0001 0.0007 31 rmsea.notclose.h0 0.0800 0.0800 32 rmr 0.0590 0.0569 33 rmr_nomean 0.0590 0.0569 34 srmr 0.0347 0.0330 35 srmr_bentler 0.0347 0.0330 36 srmr_bentler_nomean 0.0347 0.0330 37 crmr 0.0410 0.0391 38 crmr_nomean 0.0410 0.0391 39 srmr_mplus 0.0347 0.0330 40 srmr_mplus_nomean 0.0347 0.0330 41 cn_05 608.6845 599.1284 42 cn_01 788.2753 786.5616 43 gfi 0.9916 0.9923 44 agfi 0.9779 0.9768 45 pgfi 0.3777 0.3308 46 mfi 0.9889 0.9897 47 ecvi 0.0362 0.0350 （あまり大きな差は無いですが…）絶対的指標である GFI や SRMR はモデル 2 の方が良い値 36

7.7 モデルの修正 を示していますが，自由度を考慮した AGFI や RMSEA，TLI などはモデル 1 の方が良い 値を示しています。一方 CFI などはモデル 2 の方が良い値を示しています。このあたりに， 各指標の観点の違いが現れています。そのため，結果を報告する際には複数の指標を提示し， 総合的に判断する必要があるわけです。 情報量規準を見ると，AIC はモデル 2，BIC はモデル 1 のほうが大きな値になっています。 先程説明したように，AIC のほうがサンプルサイズが大きくなるほど複雑な（パラメータ数 の多い）モデルを好む傾向があり，今回の比較でもその傾向が表れたということです。 （おまけ）適合度指標の関係 図7.12に，各種適合度指標の関係を表してみました。 • 相対的指標は「独立モデルを０，飽和モデルを１」としたときの仮説モデルの位置づ けを表しています。 • 絶対的指標は「飽和モデル＝データ」からのズレを表す指標です。理論上はたぶん下 限は決まっていませんが，基本的に０から１の範囲内に収まるはずです。 • 倹約的指標はものによります。自由度で調整した値です。 • 情報量規準は独立モデル・飽和モデルと何の関係もありません。尤度を自由度で調整 しただけです。図の例では，モデル B のほうが絶対的な当てはまりは良いのですが， 自由度による調整の結果情報量規準の値はモデル A のほうが良い，という状態を表し ています。 （ここからはプラクティカルな話題です。 ） 7.7 モデルの修正 SEM では，仮説に基づいてモデルを作成（パス図に線を引く作業を）します。そして，実際 のデータと仮説モデルとで共分散行列がどの程度ずれているかを評価し，ズレが小さければ 「問題ないモデル」と判断します。ズレの大きさを評価する適合度指標には，大きく分けると 「絶対的なズレの大きさを評価するもの（SRMR や GFI など）」と「パラメータ数（自由度） を考慮して評価するもの（AGFI や RMSEA）」の 2 種類がありました。 仮説モデルにおいて適合度指標の値が悪ければ，モデルを修正してデータとの整合性 and/or モデルの倹約度を高める必要があります。データとの整合性を高めるためには，矢印を追加 してあげれば良いです。矢印を追加するということは推定するパラメータが増えるため，回 帰分析の説明変数や因子分析の因子数を増やした時と同じように，パラメータを増やした分 だけデータとモデルのズレは小さくできるはずです。一方，モデルの倹約度を高めるために は矢印を消去する必要もあります。ある矢印が，あってもなくてもデータとモデルのズレに 大きな影響を与えないのであれば，そのような矢印はなくなったほうが RMSEA などは向上 するはずです。 ということで，まずはモデルを修正してデータとの整合性を高めるための方法について説明 します。ひとつひとつ矢印を足したり引いたりしたモデルを試しまくるのはさすがに効率が 悪いので，もう少しスマートに修正案にアタリをつけたいと思います。ここでは，図7.8の 37

7.7 モデルの修正 図7.12 適合度指標の関係 SEM モデルを題材にモデルの修正を試してみます。 1 model1 <- " 2 f_A =~ Q1_1 + Q1_2 + Q1_3 3 f_E =~ Q1_6 + Q1_7 + Q1_8 4 f_A ~ age + gender 5 f_E ~ f_A 6 age ~~ gender 7 " 8 9 10 1 result1 <- sem(model1, data = dat, std.lv = TRUE) summary(result1, standardized = TRUE) lavaan 0.6.15 ended normally after 36 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters 6 38 ML NLMINB 18

7.7 モデルの修正 7 Number of observations 2432 8 9 Model Test User Model: 10 11 Test statistic 12 Degrees of freedom 13 P-value (Chi-square) 129.227 18 0.000 14 15 Parameter Estimates: 16 17 Standard errors 18 Information 19 Information saturated (h1) model Standard Expected Structured 20 21 Latent Variables: 22 Estimate Std.Err z-value P(>|z|) Std.lv Std.all Q1_1 0.602 0.031 19.255 0.000 0.628 0.447 25 Q1_2 0.888 0.030 29.426 0.000 0.926 0.790 26 Q1_3 0.786 0.030 25.796 0.000 0.819 0.627 27 f_E =~ 28 Q1_6 0.728 0.031 23.767 0.000 0.753 0.611 29 Q1_7 0.888 0.035 25.607 0.000 0.918 0.698 30 Q1_8 0.651 0.031 21.262 0.000 0.673 0.523 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 23 f_A =~ 24 31 32 Regressions: 33 34 f_A ~ 35 age 0.014 0.002 6.275 0.000 0.013 0.147 36 gender 0.522 0.053 9.894 0.000 0.501 0.236 0.252 0.030 8.402 0.000 0.254 0.254 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.204 0.105 1.945 0.052 0.204 0.039 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 37 f_E ~ 38 f_A 39 40 Covariances: 41 42 43 age ~~ gender 44 45 Variances: 46 39 47 .Q1_1 1.581 0.051 31.144 0.000 1.581 0.801 48 .Q1_2 0.517 0.046 11.273 0.000 0.517 0.376 49 .Q1_3 1.034 0.046 22.666 0.000 1.034 0.607 50 .Q1_6 0.950 0.043 22.240 0.000 0.950 0.626 51 .Q1_7 0.889 0.055 16.228 0.000 0.889 0.513 52 .Q1_8 1.205 0.044 27.533 0.000 1.205 0.727

7.7 モデルの修正 53 54 age 120.254 3.449 34.871 0.000 120.254 0.221 0.006 34.871 0.000 gender 1.000 0.221 1.000 55 .f_A 1.000 0.920 0.920 56 .f_E 1.000 0.935 0.935 適合度はすでにかなり良いんですが… 1 fitmeasures(result1) 1 npar 2 18.000 0.027 129.227 3 df pvalue baseline.chisq fmin chisq 4 18.000 0.000 2331.678 5 baseline.df baseline.pvalue cfi 6 28.000 0.000 0.952 7 tli nnfi rfi 8 0.925 0.925 0.914 9 nfi pnfi ifi 10 0.945 0.607 0.952 11 rni logl unrestricted.logl 12 0.952 -34146.289 -34081.676 13 aic bic ntotal 14 68328.579 68432.915 2432.000 15 bic2 rmsea rmsea.ci.lower 16 68375.725 0.050 0.042 17 rmsea.ci.upper rmsea.ci.level rmsea.pvalue 18 0.059 0.900 0.450 rmsea.close.h0 rmsea.notclose.pvalue rmsea.notclose.h0 19 20 0.050 0.000 0.080 21 rmr rmr_nomean srmr 22 0.308 0.308 0.036 23 srmr_bentler srmr_bentler_nomean crmr 24 0.036 0.036 0.041 25 crmr_nomean srmr_mplus srmr_mplus_nomean 26 0.041 0.036 0.036 27 cn_05 cn_01 gfi 28 544.308 656.021 0.987 29 agfi pgfi mfi 30 0.974 0.494 0.977 31 ecvi 32 0.068 7.7.1 不要なパスを消す 普通に考えると年齢 age と性別 gender は無相関だと考えるのが自然ですね。 先程の分析の結 果を見ても， 共分散 age ~~ gender の検定結果は有意にはなっていません。 summary(sem()) 40

7.7 モデルの修正
で出力される検定では，帰無仮説に「そのパラメータの母数が 0 に等しい」を設定した Wald
検定の結果を出しています。Wald 検定の中身はともかく，とりあえずこの P(>|z|) 列を見
れば，そのパラメータを削除してよいかは判断できるということです。
ということで，検定結果が有意ではなかった age ~~ gender の共分散をゼロに固定してみ
ましょう。
パラメータを制約する
1

model2 <- "

2

f_A =~ Q1_1 + Q1_2 + Q1_3

3

f_E =~ Q1_6 + Q1_7 + Q1_8

4

f_A ~ age + gender

5

f_E ~ f_A

6

age ~~ 0 * gender

7

"

lavaan の記法では，このように右辺の変数名の前に 0* を付けると，そのパラメータの値を
0 に固定することが出来ます*14 。
1

result2 <- sem(model2, data = dat, std.lv = TRUE)

2

summary(result2, standardized = TRUE)

1

lavaan 0.6.15 ended normally after 30 iterations

2
3

Estimator

4

Optimization method

5

Number of model parameters

ML
NLMINB
17

6
7

Number of observations

2432

8
9

Model Test User Model:

10
11

Test statistic

12

Degrees of freedom

13

P-value (Chi-square)

133.020
19
0.000

14
15

Parameter Estimates:

16
17

Standard errors

18

Information

19

Information saturated (h1) model

Standard
Expected
Structured

20
21
22

*14

41

Latent Variables:
Estimate

Std.Err

z-value

P(>|z|)

Std.lv

Std.all

同じように，1* としたら 1 に固定したりもできます。0 以外はあまり使い道ないかもしれませんが…

7.7 モデルの修正 23 f_A =~ 24 Q1_1 0.602 0.031 19.216 0.000 0.627 0.446 25 Q1_2 0.888 0.030 29.335 0.000 0.925 0.789 26 Q1_3 0.786 0.031 25.715 0.000 0.818 0.627 27 f_E =~ 28 Q1_6 0.728 0.031 23.763 0.000 0.753 0.611 29 Q1_7 0.888 0.035 25.602 0.000 0.918 0.698 30 Q1_8 0.651 0.031 21.259 0.000 0.673 0.523 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 31 32 Regressions: 33 34 f_A ~ 35 age 0.014 0.002 6.280 0.000 0.013 0.147 gender 0.522 0.053 9.901 0.000 0.502 0.236 0.252 0.030 8.387 0.000 0.254 0.254 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.000 0.000 36 37 f_E ~ 38 f_A 39 40 Covariances: 41 42 43 age ~~ gender 0.000 44 45 Variances: 46 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 47 .Q1_1 1.581 0.051 31.136 0.000 1.581 0.801 48 .Q1_2 0.517 0.046 11.237 0.000 0.517 0.377 49 .Q1_3 1.034 0.046 22.625 0.000 1.034 0.607 50 .Q1_6 0.950 0.043 22.238 0.000 0.950 0.626 51 .Q1_7 0.889 0.055 16.225 0.000 0.889 0.513 52 .Q1_8 1.205 0.044 27.532 0.000 1.205 0.727 53 age 120.254 3.449 34.871 0.000 120.254 1.000 54 gender 0.221 0.006 34.871 0.000 0.221 1.000 55 .f_A 1.000 0.923 0.923 56 .f_E 1.000 0.936 0.936 確かに共分散 age ~~ gender の値が推定されず，0 に固定されています。それでは，パス 削除前（result1）と削除後（result2）の適合度を比較してみましょう。ここでは，複数 のモデルの適合度を並べてくれる semTools::compareFit() という関数を使ってみます。 compareFit() 42 1 # インストールがまだの人は 2 # install.packages("semTools") 3 library(semTools) 4 # summary()にかけないと中身がみえない 5 summary(compareFit(result1, result2))

7.7 モデルの修正 1 ################### Nested Model Comparison ######################### 2 3 Chi-Squared Difference Test 4 5 Df AIC BIC Chisq Chisq diff 6 result1 18 68329 68433 129.23 7 result2 19 68330 68429 133.02 8 --- 9 Signif. codes: RMSEA Df diff Pr(>Chisq) 3.7929 0.033888 1 0.05147 . 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 10 11 ####################### Model Fit Indices ########################### 12 chisq df pvalue 13 result1 129.227† 18 14 result2 15 133.020 .000 19 rmsea .050 .000 .050† cfi tli .952† .925 .951 .927† srmr aic .036† 68328.579† .037 68330.371 bic 16 result1 17 result2 68428.911† 68432.915 18 19 ################## Differences in Fit Indices ####################### 20 21 df result2 - result1 rmsea cfi tli srmr aic bic 1 -0.001 -0.001 0.002 0.001 1.793 -4.004 最初の Nested Model Comparison というところは，ネストされたモデル間での比較です。 result2 のモデルは，result1 のモデルで age ~~ gender の値が 0 になった場合という意 味で，result1 にネストされている状態です。こうしたモデル間の比較では，尤度比検定を 行うことができます。ネストされたモデルでは，パラメータ数の減少に伴い少なからずデー タに対する当てはまりが悪化（尤度が減少）しているはずですが，もしも尤度比検定が有意 にならない場合，その当てはまりの悪化が「減らしたパラメータ数的に許容できる」という 感じになり，それならばより節約的なモデルを選ぼう，という流れになります。今回の場合， Pr(>Chisq) が 0.05 より大きいため検定は（ギリギリ）有意になっていません。したがって， 尤度比検定的には result2 を選ぼう，ということになるわけです*15 。 その下の Model Fit Indices では各種適合度指標を並べてくれています。ダガー（†）の付い ているものは，その適合度指標においてベストなモデル，という意味です。また Differences in Fit Indices は読んで字のごとく適合度指標の差分です。僅かな差ですが，RMSEA や TLI の値は改善しているようです。一方で絶対的指標である SRMR は悪化しているようで すね。 このように，単一のパスを削除する場合には，cfa() や sem() の結果の検定のところを見る のが第一歩です。ですが，統計的仮説検定の枠組みの中にあるため，サンプルサイズの影響 を受けてしまう点は注意が必要です。例えば先程の分析を，サンプルサイズを 100 人にして 行うと，さっきは有意だった f_A ~ age の回帰係数も有意ではなくなってしまいます。検 定の結果だけを鵜呑みにするのではなく，例えば実際の推定値の大きさを踏まえて決定する *15 43 例によって，統計的仮説検定なのでサンプルサイズの影響を受ける点に注意してください。

7.7 モデルの修正 などの対応が必要でしょう。 サンプルサイズを減らしてみると 1 result1_100 <- sem(model1, data = dat[1:100, ], std.lv = TRUE) 2 summary(result1_100, standardized = TRUE) 1 lavaan 0.6.15 ended normally after 35 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters ML NLMINB 18 6 7 Number of observations 100 8 9 Model Test User Model: 10 11 Test statistic 12 Degrees of freedom 13 P-value (Chi-square) 29.972 18 0.038 14 15 Parameter Estimates: 16 17 Standard errors Standard 18 Information Expected 19 Information saturated (h1) model Structured 20 21 Latent Variables: 22 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 23 f_A =~ 24 Q1_1 0.535 0.151 3.540 0.000 0.576 0.394 25 Q1_2 0.926 0.147 6.301 0.000 0.996 0.857 26 Q1_3 0.674 0.141 4.781 0.000 0.725 0.537 27 f_E =~ 28 Q1_6 0.764 0.126 6.085 0.000 0.840 0.679 29 Q1_7 0.829 0.135 6.163 0.000 0.911 0.690 30 Q1_8 0.822 0.132 6.222 0.000 0.904 0.698 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 31 32 Regressions: 33 34 f_A ~ 35 age 0.011 0.010 1.050 0.294 0.010 0.115 36 gender 0.720 0.254 2.830 0.005 0.669 0.330 0.424 0.148 2.874 0.004 0.415 0.415 37 f_E ~ 38 f_A 39 40 44 Covariances:

7.7 モデルの修正 41 42 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 1.028 0.580 1.772 0.076 1.028 0.180 age ~~ 43 gender 44 45 Variances: 46 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 47 .Q1_1 1.801 0.273 6.606 0.000 1.801 0.844 48 .Q1_2 0.360 0.234 1.535 0.125 0.360 0.266 49 .Q1_3 1.296 0.225 5.769 0.000 1.296 0.711 50 .Q1_6 0.823 0.172 4.788 0.000 0.823 0.539 51 .Q1_7 0.915 0.197 4.647 0.000 0.915 0.524 52 .Q1_8 53 age 54 gender 0.860 0.190 4.533 0.000 0.860 0.513 133.886 18.934 7.071 0.000 133.886 1.000 0.244 0.034 7.071 0.000 0.244 1.000 55 .f_A 1.000 0.864 0.864 56 .f_E 1.000 0.828 0.828 7.7.2 新しいパスを増やす 続いて，新しいパスを追加することを考えたいと思います。パスを追加する場合，基本的に 絶対的な適合度は上昇するため，先程の尤度比検定の考え方とは反対に「パラメータを 1 個 追加する価値があるだけの適合度の改善が見られる」ようなパラメータは追加してあげれば よいわけです。modificationIndices() という関数は，現状モデル内でパラメータが自由 推定になっていない全てのパス*16 について，「そこにパスを追加するとどれだけ適合度が改 善するか」をまとめて表してくれます。 修正指標 1 # 実はmodindices()でも良い 2 # 引数sort = TRUEを与えると，改善度の大きい順に並べてくれる 3 modificationIndices(result2, sort = TRUE) 1 rhs mi epc sepc.lv sepc.all sepc.nox 2 24 -0.203 -0.144 -0.144 3 35 Q1_2 ~~ Q1_8 23.532 0.114 0.114 0.145 0.145 4 39 Q1_6 ~~ Q1_7 19.848 0.527 0.527 0.574 0.574 5 23 f_A =~ Q1_8 19.845 0.134 0.139 0.108 0.108 6 30 Q1_1 ~~ Q1_7 9.220 -0.092 -0.092 -0.078 -0.078 7 26 f_E =~ Q1_3 7.717 0.087 0.089 0.069 0.069 8 37 Q1_3 ~~ Q1_7 5.463 0.062 0.062 0.065 0.065 f_E =~ Q1_1 33.281 -0.196 9 41 Q1_7 ~~ Q1_8 4.451 -0.196 -0.196 -0.189 -0.189 10 21 f_A =~ Q1_6 4.451 -0.062 -0.065 -0.053 -0.053 11 47 age 4.344 0.611 0.056 0.056 *16 45 lhs op ~ f_E 0.591 「自由推定になっていない」ということは，例えば 1* によって係数を 0 以外の定数に固定している場所も検 討対象にしてくれる，ということです。

7.7 モデルの修正

12

29

13

49 gender

14

Q1_1 ~~ Q1_6
~

f_A

3.977 -0.058
3.790

0.121

-0.058

-0.047

-0.047

0.126

0.268

0.268

[ reached 'max' / getOption("max.print") -- omitted 20 rows ]

引数 sort = TRUE を与えると，改善度（出力の一番左，mi 列）の大きい順に並べてくれま
す。したがって今回の場合，一つだけパスを追加するならば f_E =~ Q1_1 というパスを追
加すると適合度が最も改善する，ということです。出力では mi 列に続いて，
そのパスを自由推定にした際のパラメータの変化量の予測値。パスが引かれていない

epc

場合は「0 からの変化量」なので推定値そのものの予測値を意味する。

sepc.lv 潜在変数 lv を標準化した際の EPC。
sepc.all 全ての変数を標準化した際（標準化解）の EPC。
sepc.nox 観測変数の外生変数以外を全て標準化した際の EPC。
が表示されています。きっと多くの場合では mi が大きいところは sepc.all が大きくなっ
ていることでしょう。ということで，パスを追加する際には mi が大きいほうから試してみ
ると良いでしょう。さっそく f_E =~ Q1_1 を追加して推定してみます。

modindices の仰せのままに
1

model3 <- "

2

f_A =~ Q1_1 + Q1_2 + Q1_3

3

f_E =~ Q1_6 + Q1_7 + Q1_8 + Q1_1

4

f_A ~ age + gender

5

f_E ~ f_A

6

age ~~ 0*gender

7

"

8
9
10

1

result3 <- sem(model3, data = dat, std.lv = TRUE)
summary(result3, standardized = TRUE)

lavaan 0.6.15 ended normally after 31 iterations

2
3

Estimator

4

Optimization method

5

Number of model parameters

ML
NLMINB
18

6
7

Number of observations

2432

8
9

Model Test User Model:

10

46

11

Test statistic

12

Degrees of freedom

13

P-value (Chi-square)

98.800
18
0.000

7.7 モデルの修正 14 15 Parameter Estimates: 16 17 Standard errors Standard 18 Information Expected 19 Information saturated (h1) model Structured 20 21 Latent Variables: 22 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 23 f_A =~ 24 Q1_1 0.678 0.034 19.682 0.000 0.707 0.503 25 Q1_2 0.867 0.029 29.957 0.000 0.904 0.772 26 Q1_3 0.792 0.030 26.483 0.000 0.826 0.633 27 f_E =~ 28 Q1_6 0.720 0.030 23.860 0.000 0.752 0.610 29 Q1_7 0.879 0.034 25.758 0.000 0.918 0.698 30 Q1_8 0.647 0.030 21.333 0.000 0.675 0.524 31 Q1_1 -0.201 0.034 -5.875 0.000 -0.210 -0.149 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 32 33 Regressions: 34 35 f_A ~ 36 age 0.014 0.002 6.466 0.000 0.014 0.152 37 gender 0.529 0.053 9.982 0.000 0.507 0.238 0.287 0.031 9.131 0.000 0.287 0.287 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.000 0.000 38 f_E ~ 39 f_A 40 41 Covariances: 42 43 44 age ~~ gender 0.000 45 46 Variances: 47 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 48 .Q1_1 1.514 0.052 29.360 0.000 1.514 0.767 49 .Q1_2 0.555 0.042 13.069 0.000 0.555 0.404 50 .Q1_3 1.021 0.044 23.011 0.000 1.021 0.600 51 .Q1_6 0.952 0.042 22.594 0.000 0.952 0.628 52 .Q1_7 0.889 0.054 16.558 0.000 0.889 0.513 53 .Q1_8 1.202 0.044 27.589 0.000 1.202 0.725 54 age 120.254 3.449 34.871 0.000 120.254 1.000 55 gender 0.221 0.006 34.871 0.000 0.221 1.000 56 .f_A 1.000 0.920 0.920 57 .f_E 1.000 0.918 0.918 パス追加前（result2）と追加後（result3）の適合度を比較してみましょう。 47

7.7 モデルの修正

1

# summary()のほうで引数fit.measuresを指定すると，出力する適合度指標を選べます

2

# 試しにAGFIを追加してみます

3

summary(compareFit(result2, result3), fit.measures = c("chisq", "df",
"pvalue", "cfi", "tli", "rmsea", "srmr", "agfi", "aic", "bic"))

1

################### Nested Model Comparison #########################

2
3

Chi-Squared Difference Test

4
5

Df

AIC

BIC

Chisq Chisq diff

6

result3 18 68298 68402

7

result2 19 68330 68429 133.02

8

---

9

Signif. codes:

RMSEA Df diff Pr(>Chisq)

98.80
34.22 0.11688

1

4.921e-09 ***

0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

10
11

####################### Model Fit Indices ###########################

12

chisq df pvalue

13

result3 98.800† 18

14

result2 133.020

15

cfi

tli

rmsea

srmr

agfi

.000 .965† .945† .043† .029† .980†

19

.000

aic

.951

16

result3 68298.151† 68402.488†

17

result2

68330.371

.927

.050

.037

.975

bic
68428.911

18
19

################## Differences in Fit Indices #######################

20
21

df
result2 - result3

cfi

tli rmsea

srmr

agfi

aic

bic

1 -0.014 -0.018 0.007 0.008 -0.005 32.22 26.424

尤度比検定の結果は有意になっているため，パス追加前（result2）の方が当てはまりが有
意に悪い，ということになりました。ほかの適合度指標を見ても，どうやら確かに追加した
あとの方が改善しているようです。

7.7.3 実際のところ
ここまで，モデルを修正するための候補を統計的な視点から紹介しました。しかし SEM で
重要なのは，因子分析と同じようにその結果をきちんと解釈できるか（解釈可能性）です。
実際にモデルを修正するかどうかは，そのパスが持つ意味をよく考えて行う必要があります。
先程の修正では，f_E =~ Q1_1 を追加したことによって図7.13のような状態になりました。
つまり完全な単純構造ではないと考えたほうがデータ的には当てはまりが良い，ということ
です。ですが修正指標（mi）も，尤度に基づく指標です。つまりサンプルサイズが大きいほ
ど当てはまりは大幅に改善するように見えてしまいます。同様に，モデルが複雑になるほど
自由度も大きくなり共分散構造の乖離度も大きくなりやすいと考えられます。したがって，
一見無意味なパスでもそれなりの改善を示す可能性が生じます。パスを追加・削除すること
が自分の仮説と照らし合わせて受け入れられるかをよく考えた上で修正を行うようにしまし
ょう。

48

7.7 モデルの修正 図7.13 より当てはまりの良いモデル。でもいいの？ホントにそれで 実践場面でよく見られる修正としては，図7.14のように独自因子の間の相関を追加するパタ ーンです。例えば時系列モデルにおいて，異なる時期に聞かれた同一の項目間においたり，あ るいは DIF のように特定の属性などが影響していると言えそうな場合には妥当なパスの追 加だと思われます。このような場合にも，パスを追加する場合にはなぜそのパスを追加する べき（しても問題ない）と考えたかをきちんと説明するようにしましょう。裏を返せば，そ のパスが意味するところが説明できないような場所には，いくら適合度が改善するにしても パスを追加しないほうが良いだろう，ということです。 図7.14 独自因子間の相関 また，モデルの修正は割と終わりなき戦いです。というのも，modificationIndices() で 出力される mi は必ず正の値を取ります。絶対的な適合度指標と同じように，パスを追加す ればするほど当てはまりは良くなるということです。モデルの修正を行う場合には， 1. まず仮説で立てたモデルを基準として 2. もしそのままでは当てはまりが良くない（適合度指標が許容範囲を超えている）とい う場合にのみ 3. もともとの仮説で受け入れられるパスの追加・削除を少しずつ行い 4. 適合度指標が許容範囲になったらその時点でストップ くらいに考えておくのが良いのではないかと思います。きっと仮説で立てたモデルが理論的 には最も整合的なはずなので，そこからあまり変えないようにできればそれに越したことは 無いわけです。適合度指標の許容範囲についてはゼロイチで判断できるものではないですが， 自分が納得できる＆査読者などを納得させられるくらいの値であれば良いでしょう。 49

7.8 カテゴリカル変数の SEM 7.8 カテゴリカル変数の SEM 資料の前半で見てきたように，SEM は因子分析や回帰分析を内包する分析手法です。したが って，カテゴリカル変数を扱う場合も因子分析や回帰分析と同じように考えれば OK です。 7.8.1 内生変数がカテゴリカルな回帰分析・パス解析 内生変数，つまり結果変数がカテゴリカルな場合，ポリコリック相関のように，本来連続的 なものが閾値で切り分けられてカテゴリになっている可能性があります。このような場合は， ダミー変数化しても意味がありません。図7.15のように，性別と勤勉性（Q1_C）が学歴に与 える影響を考える回帰分析モデルを実行してみましょう。 図7.15 内生変数がカテゴリカルな場合 lavaan には，内生観測変数がカテゴリカルな場合のための引数が用意されています。 内生変数がカテゴリカル 1 # まずは普通にモデル式を作成 2 model <- " 3 education ~ gender + Q1_C 4 " 5 # 引数orderedに指定する 6 result <- sem(model, data = dat, ordered = "education") 7 summary(result) 1 lavaan 0.6.15 ended normally after 3 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters DWLS NLMINB 6 6 7 8 Number of observations Used Total 2232 2432 Standard Scaled 0.000 0.000 0 0 9 10 Model Test User Model: 11 12 Test Statistic 13 Degrees of freedom 14 15 50 Parameter Estimates:

7.8 カテゴリカル変数の SEM 16 17 Standard errors 18 Information 19 Information saturated (h1) model Robust.sem Expected Unstructured 20 21 Regressions: 22 Estimate Std.Err z-value P(>|z|) 23 education ~ 24 gender 0.007 0.047 0.150 0.881 25 Q1_C 0.006 0.005 1.233 0.217 Estimate Std.Err z-value P(>|z|) 26 27 Intercepts: 28 29 .education 0.000 30 31 Thresholds: 32 Estimate Std.Err z-value P(>|z|) 33 education|t1 -1.217 0.123 -9.858 0.000 34 education|t2 -0.706 0.123 -5.757 0.000 35 education|t3 0.608 0.124 4.914 0.000 36 education|t4 1.117 0.125 8.961 0.000 Estimate Std.Err z-value P(>|z|) Std.Err z-value P(>|z|) 37 38 Variances: 39 40 .education 1.000 41 42 Scales y*: 43 44 Estimate education 1.000 結果を見ると，Intercepts: と Variances: のところにある education に，潜在変数であ ることを示す. がついて.education になっています。そしてこれらの推定値がそれぞれ 0, 1 になっています。つまりこれは，education という変数の背後にある連続的なものが平均 0，分散 1 に標準化された値として考慮されている，ということです。 そのもとで，Thresholds: のところには各カテゴリの閾値が表示されています。図 4.8 に当 てはめると図7.16のような状態だということです。 1 Warning: Using `size` aesthetic for lines was deprecated in ggplot2 2 i Please use `linewidth` instead. 3 This warning is displayed once every 8 hours. 4 Call `lifecycle::last_lifecycle_warnings()` to see where this warning 5 was generated. 3.4.0. 51

7.8 カテゴリカル変数の SEM 1 2 -1.208 3 -0.695 図7.16 4 5 0.628 1.137 education の閾値 7.8.2 （おまけ）外生変数がカテゴリカルな回帰分析・パス解析 「外生変数がカテゴリカル」ということは，回帰分析の説明変数がカテゴリカル，ということ です。例えば education を説明変数に加えた，図7.17の回帰分析モデルを考えます。 図7.17 外生変数がカテゴリカルな場合 education は順序カテゴリカル変数ですが，順序を無視して考えるならば，ダミー変数 へ置き換えたら良いですね。ダミー変数をデータフレームに追加する簡単な方法として， contr.treatment() 関数を利用する方法があります。 一気にダミー変数を作る 1 # ダミーコーディングの形を作る関数 2 # 数字はカテゴリ数 3 contr.treatment(5) 1 2 3 4 5 2 1 0 0 0 0 3 2 1 0 0 0 4 3 0 1 0 0 5 4 0 0 1 0 6 5 0 0 0 1 contr.treatment() 関数では第一カテゴリを基準カテゴリとしたダミー変数を作成してく れます。 1 1 2 52 contr.treatment(5)[dat[, "education"], ] 2 3 4 5 <NA> NA NA NA NA

7.8 カテゴリカル変数の SEM

3

<NA> NA NA NA NA

4

<NA> NA NA NA NA

5

<NA> NA NA NA NA

6

<NA> NA NA NA NA

7

3

8

<NA> NA NA NA NA

9

2

0
1

1
0

0
0

0
0

10

<NA> NA NA NA NA

11

1

12

<NA> NA NA NA NA

13

<NA> NA NA NA NA

14

1

15

<NA> NA NA NA NA

16

<NA> NA NA NA NA

17

<NA> NA NA NA NA

18

<NA> NA NA NA NA

19

<NA> NA NA NA NA

20

<NA> NA NA NA NA

21

<NA> NA NA NA NA

22

5

0

0

0

1

23

2

1

0

0

0

24

1

0

0

0

0

25

3

0

1

0

0

26

5

0

0

0

1

27

0

0

0

0

0

0

0

0

[ reached getOption("max.print") -- omitted 2407 rows ]

contr.treatment() で作成したダミー行列の dat[,"education"] 行目を取る，というこ
とは，

• education == 1 の人なら 0 0 0 0 の行
• education == 2 の人なら 1 0 0 0 の行
• education == 3 の人なら 0 1 0 0 の行
• education == 4 の人なら 0 0 1 0 の行
• education == 5 の人なら 0 0 0 1 の行
• education == NA の人なら NA NA NA NA の行
が返ってくるわけです。これをもとの dat の存在しない列名に代入してあげれば OK ですね。
1

dat[, paste0("edu_", 2:5)] <- contr.treatment(5)[dat[, "education"], ]

こうして作成したダミー変数を入れたモデルを作成して実行しましょう。

53

1

model <- "

2

Q1_A ~ age + edu_2 + edu_3 + edu_4 + edu_5 + 1

3

"

7.8 カテゴリカル変数の SEM 4 5 summary(sem(model, data = dat)) 1 lavaan 0.6.15 ended normally after 40 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters ML NLMINB 7 6 7 8 Number of observations Used Total 2232 2432 9 10 Model Test User Model: 11 12 Test statistic 13 Degrees of freedom 0.000 0 14 15 Parameter Estimates: 16 17 Standard errors 18 Information 19 Information saturated (h1) model Standard Expected Structured 20 21 Regressions: 22 Estimate Std.Err z-value P(>|z|) 23 Q1_A ~ 24 age 25 edu_2 26 edu_3 1.097 0.336 3.263 0.001 27 edu_4 -0.132 0.393 -0.337 0.736 28 edu_5 0.413 0.394 1.046 0.295 Estimate Std.Err z-value P(>|z|) 20.897 0.382 54.641 0.000 Estimate Std.Err z-value P(>|z|) 18.842 0.564 33.407 0.000 0.068 0.009 7.470 0.000 -0.207 0.417 -0.497 0.619 29 30 Intercepts: 31 32 .Q1_A 33 34 Variances: 35 36 .Q1_A Number of observations のところに Total という値が追加されています。当然ですが欠 測値がある場合にはそのままだと分析が出来ないため，デフォルトでは「分析に使用する変 数に欠測が一つもない人」だけを使って分析をしているのです。 54

7.8 カテゴリカル変数の SEM 順序カテゴリをダミー化するなら 順序カテゴリの順序性を保ったままダミー化する場合は，下の表のようにダミーの効 果を加算していく，という方法が考えられます。カテゴリ 2 では edu_2 の効果のみ， カテゴリ 3 では edu_2+edu_3 の効果，…というようになります。 表7.2 順序カテゴリのコーディング例 ダミー カテゴリ 2 3 4 5 1 0 0 0 0 2 1 0 0 0 3 1 1 0 0 4 1 1 1 0 5 1 1 1 1 7.8.3 カテゴリカルな確認的因子分析 確 認 的 因 子 分 析 の 場 合 も， 引 数 ordered を 使 え ば OK で す。 全 て の 変 数 を 名 指 し で ordered=paste0("Q1_",c(1,2,3,6,7,8)) という感じで指定しても良いのですが，も し全ての内生変数である観測変数がカテゴリの場合*17 ，ordered = TRUE とすることで全 てカテゴリカルになります。 ということで，前回行った図7.3の単純構造 2 因子モデルで試してみましょう。 カテゴリカル検証的因子分析 1 # まずは普通にモデル式を作成 2 model <- " 3 f_1 =~ Q1_1 + Q1_2 + Q1_3 4 f_2 =~ Q1_6 + Q1_7 + Q1_8 5 " 6 # 引数orderedに指定する 7 result <- cfa(model, data = dat, std.lv = TRUE, ordered = TRUE) 8 summary(result, standardized = TRUE) 1 lavaan 0.6.15 ended normally after 14 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters DWLS NLMINB 37 6 *17 55 図7.15の Q1_C を和得点ではなく因子得点にする場合，ordered = TRUE にすると，自動的に education も カテゴリ変数として扱われる，ということです。このように観測変数の回帰と因子分析を組み合わせた場合 には意図せぬカテゴリ化を行ってしまう可能性があるので気をつけましょう。

7.8 カテゴリカル変数の SEM 7 Number of observations 2432 8 9 Model Test User Model: 10 Standard Scaled 72.885 100.302 11 Test Statistic 12 Degrees of freedom 13 P-value (Chi-square) 14 Scaling correction factor 0.733 15 Shift parameter 0.815 16 8 8 0.000 0.000 simple second-order correction 17 18 Parameter Estimates: 19 20 Standard errors 21 Information 22 Information saturated (h1) model Robust.sem Expected Unstructured 23 24 Latent Variables: 25 Estimate Std.Err z-value P(>|z|) Std.lv Std.all Q1_1 0.472 0.020 24.163 0.000 0.472 0.472 Q1_2 0.854 0.021 41.558 0.000 0.854 0.854 29 Q1_3 0.681 0.019 36.343 0.000 0.681 0.681 30 f_2 =~ 31 Q1_6 0.649 0.020 32.847 0.000 0.649 0.649 32 Q1_7 0.729 0.019 37.962 0.000 0.729 0.729 33 Q1_8 0.561 0.018 30.590 0.000 0.561 0.561 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.271 0.027 10.191 0.000 0.271 0.271 Estimate Std.Err z-value P(>|z|) 26 f_1 =~ 27 28 34 35 Covariances: 36 37 f_1 ~~ 38 f_2 39 40 Intercepts: 41 Std.lv Std.all 42 .Q1_1 0.000 0.000 0.000 43 .Q1_2 0.000 0.000 0.000 44 .Q1_3 0.000 0.000 0.000 45 .Q1_6 0.000 0.000 0.000 46 .Q1_7 0.000 0.000 0.000 47 .Q1_8 0.000 0.000 0.000 48 f_1 0.000 0.000 0.000 49 f_2 0.000 0.000 0.000 Std.lv Std.all 50 51 52 56 Thresholds: Estimate Std.Err z-value P(>|z|)

7.9 多母集団同時分析 53 Q1_1|t1 -1.887 0.051 -36.934 0.000 -1.887 -1.887 54 Q1_1|t2 -1.234 0.034 -36.437 0.000 -1.234 -1.234 55 Q1_1|t3 -0.743 0.028 -26.414 0.000 -0.743 -0.743 56 Q1_1|t4 -0.326 0.026 -12.588 0.000 -0.326 -0.326 57 Q1_1|t5 0.434 0.026 16.487 0.000 0.434 0.434 58 Q1_2|t1 -2.144 0.064 -33.743 0.000 -2.144 -2.144 59 Q1_2|t2 -1.534 0.040 -38.430 0.000 -1.534 -1.534 60 Q1_2|t3 -1.179 0.033 -35.716 0.000 -1.179 -1.179 61 Q1_2|t4 -0.478 0.027 -18.047 0.000 -0.478 -0.478 62 Q1_2|t5 0.481 0.027 18.127 0.000 0.481 0.481 63 Q1_3|t1 -1.829 0.049 -37.433 0.000 -1.829 -1.829 64 Q1_3|t2 -1.303 0.035 -37.181 0.000 -1.303 -1.303 65 Q1_3|t3 -0.963 0.030 -31.885 0.000 -0.963 -0.963 Q1_3|t4 -0.328 0.026 -12.668 0.000 -0.328 -0.328 Std.lv Std.all 66 67 [ reached getOption("max.print") -- omitted 16 rows ] 68 69 Variances: 70 Estimate Std.Err z-value P(>|z|) 71 .Q1_1 0.777 0.777 0.777 72 .Q1_2 0.270 0.270 0.270 73 .Q1_3 0.536 0.536 0.536 74 .Q1_6 0.578 0.578 0.578 75 .Q1_7 0.469 0.469 0.469 76 .Q1_8 0.686 0.686 0.686 77 f_1 1.000 1.000 1.000 78 f_2 1.000 1.000 1.000 79 80 Scales y*: 81 Std.lv Std.all 82 Q1_1 Estimate 1.000 Std.Err z-value P(>|z|) 1.000 1.000 83 Q1_2 1.000 1.000 1.000 84 Q1_3 1.000 1.000 1.000 85 Q1_6 1.000 1.000 1.000 86 Q1_7 1.000 1.000 1.000 87 Q1_8 1.000 1.000 1.000 ordered を指定せずにやった結果と比べると，因子負荷の値（Std.all 列）が少し大きくな っていることがわかります。このあたりは探索的因子分析のときと同じですね。 7.9 多母集団同時分析 分析結果の一般化可能性を高めるためには「異なる属性間で一貫して同じ結果が得られるか」 を検討するのが一つの手段となります。あるいは異なる属性・集団の間での差異に関心があ る場合にも，グループごとに分析を行い結果の一貫性を検討することが多々あります。例え ば因子分析によって得られる因子得点に関して，男性と女性で平均点に差があるか・どちら のほうが平均点が高いかを検討したいとしましょう。SEM では多母集団同時分析という枠 57

7.9 多母集団同時分析 組みによってこうした分析を行うことができます。 7.9.1 多母集団同時分析の基本的な考え方 基本的な考え方はシンプルです。図7.18のように，同じモデルをグループごとに適用し，対 応するパラメータの推定値が同じかどうかを評価していきます。なおここでは，上付き添字 (𝑀) はグループごとに異なる係数を表すと思ってください。つまり 𝑏1 (𝐹 ) と 𝑏1 は同じ因子から 同じ項目への因子負荷ですが，異なるパラメータとして別々に推定される，ということです。 図7.18 多母集団同時分析 lavaan では，引数 group に，グループを表す変数の名前を指示してあげればすぐに分析が 可能です。 とりあえずグループごとに分析してみる 1 model <- "f_A =~ Q1_1 + Q1_2 + Q1_3 + Q1_4" 2 3 # もちろんordered=TRUEにしてもOKです。 4 # ここでは簡略化のためorderedは指定しません。 5 result_base <- cfa(model, data = dat, std.lv = TRUE, group = "gender") 6 summary(result_base, standardized = TRUE) 1 lavaan 0.6.15 ended normally after 43 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters ML NLMINB 24 6 7 Number of observations per group: 8 1 802 9 2 1630 10 11 58 Model Test User Model:

7.9 多母集団同時分析 12 13 Test statistic 14 Degrees of freedom 15 P-value (Chi-square) 16 Test statistic for each group: 34.458 4 0.000 17 1 7.472 18 2 26.986 19 20 Parameter Estimates: 21 22 Standard errors 23 Information 24 Information saturated (h1) model Standard Expected Structured 25 26 27 Group 1 [1]: 28 29 Latent Variables: 30 Estimate Std.Err Q1_1 0.549 Q1_2 0.857 34 Q1_3 35 Q1_4 31 f_A =~ 32 33 z-value P(>|z|) Std.lv Std.all 0.058 9.544 0.000 0.549 0.385 0.051 16.813 0.000 0.857 0.682 0.965 0.054 17.784 0.000 0.965 0.731 0.761 0.059 13.003 0.000 0.761 0.515 36 37 Intercepts: 38 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 39 .Q1_1 4.278 0.050 85.060 0.000 4.278 3.004 40 .Q1_2 4.494 0.044 101.352 0.000 4.494 3.579 41 .Q1_3 4.347 0.047 93.246 0.000 4.347 3.293 42 .Q1_4 4.434 0.052 84.932 0.000 4.434 2.999 43 f_A 0.000 0.000 0.000 44 45 Variances: 46 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 47 .Q1_1 1.728 0.093 18.624 0.000 1.728 0.852 48 .Q1_2 0.842 0.070 11.949 0.000 0.842 0.534 49 .Q1_3 0.811 0.082 9.916 0.000 0.811 0.465 50 .Q1_4 1.607 0.094 17.037 0.000 1.607 0.735 51 f_A 1.000 1.000 1.000 Std.lv Std.all 52 53 54 Group 2 [2]: 55 56 57 59 Latent Variables: Estimate Std.Err z-value P(>|z|)

7.9 多母集団同時分析 58 f_A =~ 59 Q1_1 0.575 0.039 14.834 0.000 0.575 0.420 60 Q1_2 0.811 0.032 25.073 0.000 0.811 0.738 61 Q1_3 0.852 0.037 23.087 0.000 0.852 0.665 62 Q1_4 0.695 0.041 16.883 0.000 0.695 0.474 63 64 Intercepts: 65 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 66 .Q1_1 4.748 0.034 139.992 0.000 4.748 3.467 67 .Q1_2 4.953 0.027 182.102 0.000 4.953 4.510 68 .Q1_3 4.728 0.032 149.135 0.000 4.728 3.694 69 .Q1_4 4.818 0.036 132.828 0.000 4.818 3.290 70 f_A 0.000 0.000 0.000 71 72 Variances: 73 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 74 .Q1_1 1.544 0.059 26.000 0.000 1.544 0.824 75 .Q1_2 0.549 0.041 13.298 0.000 0.549 0.455 76 .Q1_3 0.913 0.052 17.571 0.000 0.913 0.557 77 .Q1_4 1.662 0.066 25.038 0.000 1.662 0.775 78 f_A 1.000 1.000 1.000 結果を見ると，Number of observations per group:（7 行目）のところで，各グループ での人数を表示しています。 そして Group 1 [1]:（27 行目から）と Group 2 [2]:（54 行目から）という形で，それぞ れのグループに対する推定結果が表示されています。この段階では，グループごとに因子得 点が標準化されている（52, 78 行目）ため，当然ながら平均値などの比較はできません。そ して因子得点をグループ間で比較可能にするためには，因子負荷と切片がともにグループ間 で等しいという制約を置く必要があります。多母集団同時分析の目的の一つは，因子構造が 同じかどうかを検証することですが，同一性が満たされて初めて，異なるグループ間での因 子得点の平均が比較可能になります。 ここでいう「切片」は，回帰分析の「切片」と同じようなものです。因子分析では説明変数 にあたる因子得点が潜在変数のため，少し厄介なことになっているのです。これまで 1 因子 の因子分析モデルは 𝐲𝑖 = 𝐟 𝑏𝑖 + 𝛆𝑖 (7.48) という形を取っていました。ですが本当は，回帰分析での 𝑏0 にあたる切片項が存在しており， 𝐲𝑖 = 𝜏 + 𝐟 𝑏 𝑖 + 𝛆 𝑖 (7.49) という式になっています。とはいえ通常，観測変数 𝐲𝑖 は標準化されているため，暗黙に 𝜏 = 0 となります。そのためにこれまでは 𝜏 を無視してきました。ですが多母集団同時分析では， この切片項まで考慮しないと因子得点の比較ができなくなってしまいます*18 。例えば，デー *18 60 多母集団同時分析以外でも切片を考えても良いのですが，グループ間での比較をしない限りはさほど重要で

7.9 多母集団同時分析 タ全体で標準化された（全体平均が 0 の）観測変数 𝐲𝑖 の平均値が，男性のみで計算すると 0.3, 女性のみで計算すると-0.3 だった，つまり男性のほうが平均的に高い得点をつけていた としましょう。この 0.6 の差は「男性のほうが関連する因子得点 𝐟 が高い」ということを意 味しているように見えるかもしれませんが，まだそうではない可能性を持っています。グル ープごとに異なる切片 𝜏 (𝑔) があるとすると，先程の式は 𝐲𝑖 = 𝜏 (𝑔) + 𝐟 𝑏𝑖 + 𝛆𝑖 (7.50) となります。ここで，もし男性の切片 𝜏 (M) = 0.3 で女性の切片 𝜏 (F) = −0.3 だとすると，観 測変数𝐲𝑖 の平均値の差0.6は切片の差に起因するものになってしまい，因子得点 𝐟 の部分では 平均値には差がなくなります。回帰分析で言えば gender というダミー変数に対する回帰係 数が 0.6 だった，ということです。例えばその項目が因子とは別の理由で男性の方が「あて はまる」と答えやすい項目などのケースでは，こういったことが起こると考えられます。そ こで，もし 𝜏 (M) = 𝜏 (F) という制約をおけば，観測変数 𝐲𝑖 の平均値の差 0.6 はどうやら因子 得点 𝐟 の部分の違いによるものだ，ということができるわけです。こういった理由から，グ ループ間で因子得点を比較可能にするためには，因子負荷だけでなく切片にも等値制約を置 く必要があるのです。 7.9.2 測定不変性の検証 多母集団同時分析では，先程実行したモデル（配置不変モデル：configural invariance など と呼ばれます）をベースラインとして，グループ間で等値制約をどこまで追加しても当ては まりが悪化しないかを検証していきます。ということで，もし配置不変モデルの段階で適合 度が悪かったら終了です。そのような場合は，少なくとも group で分ける前から適合度が悪 いはず… 配置不変モデルの適合度が問題なければ制約を追加していきましょう。ここからは 弱測定不変モデル：weak invariance 因子負荷が同じ 強測定不変モデル：strong invariance さらに切片が同じ 厳密な測定不変モデル：strict invariance （ほぼ完全な測定不変モデル） （完全な測定不変モデル） さらに独自因子の分散（共分散）が同じ さらに共通因子の分散共分散が同じ さらに因子得点の平均が同じ という順番で徐々に制約を追加していきます*19 。つまり最低でも強測定不変モデルが許容さ れるならば，因子得点の平均の比較が可能になるというわけです。また，もしも一番下まで クリア出来てしまった場合には，多母集団同時分析を行う必要はありません。 胸を張って， 単一母集団モデルとしての結果を報告しましょう。 lavaan でこれらの等値制約をかける場合には，引数 group.equal に制約をかけたいパラメ ータの種類を指定します。指定は何種類でも OK です。 *19 61 はないのでこれまで無視してきました。 下 2 つのモデルには決まった名称が無いため，カッコをつけています。

7.9 多母集団同時分析 指定 制約されるもの loadings 因子負荷 intercepts 切片 means 因子得点の平均 thresholds 閾値（ordered のときの） regressions 回帰係数 residuals 独自因子の分散 residual.covariances 独自因子の共分散 lv.variance 共通因子の分散 lv.covariance 共通因子の共分散 ということで，まずは弱測定不変モデル（因子負荷が同じ）を実行してみます。 弱測定不変モデル 1 # 因子負荷に等値制約を置く 2 3 result_weak <- cfa(model, data = dat, group = "gender", std.lv = TRUE, g roup.equal = c("loadings")) summary(result_weak, standardized = TRUE) 1 lavaan 0.6.15 ended normally after 43 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters 6 Number of equality constraints ML NLMINB 25 4 7 8 Number of observations per group: 9 1 802 10 2 1630 11 12 Model Test User Model: 13 14 Test statistic 15 Degrees of freedom 16 P-value (Chi-square) 17 Test statistic for each group: 35.926 7 0.000 18 1 8.416 19 2 27.509 20 21 Parameter Estimates: 22 62 23 Standard errors Standard 24 Information Expected

7.9 多母集団同時分析 25 Information saturated (h1) model Structured 26 27 28 Group 1 [1]: 29 30 Latent Variables: 31 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 32 f_A =~ 33 Q1_1 (.p1.) 0.596 0.038 15.684 0.000 0.596 0.415 34 Q1_2 (.p2.) 0.869 0.039 22.169 0.000 0.869 0.692 35 Q1_3 (.p3.) 0.937 0.042 22.576 0.000 0.937 0.714 36 Q1_4 (.p4.) 0.756 0.041 18.308 0.000 0.756 0.512 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 37 38 Intercepts: 39 40 .Q1_1 4.278 0.051 84.317 0.000 4.278 2.977 41 .Q1_2 4.494 0.044 101.301 0.000 4.494 3.577 42 .Q1_3 4.347 0.046 93.738 0.000 4.347 3.310 43 .Q1_4 4.434 0.052 84.975 0.000 4.434 3.001 44 f_A 0.000 0.000 0.000 45 46 Variances: 47 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 48 .Q1_1 1.709 0.092 18.650 0.000 1.709 0.828 49 .Q1_2 0.823 0.062 13.282 0.000 0.823 0.522 50 .Q1_3 0.846 0.068 12.457 0.000 0.846 0.491 51 .Q1_4 1.612 0.091 17.623 0.000 1.612 0.738 52 f_A 1.000 1.000 1.000 53 54 55 Group 2 [2]: 56 57 Latent Variables: 58 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 59 f_A =~ 60 Q1_1 (.p1.) 0.596 0.038 15.684 0.000 0.552 0.405 61 Q1_2 (.p2.) 0.869 0.039 22.169 0.000 0.804 0.732 62 Q1_3 (.p3.) 0.937 0.042 22.576 0.000 0.867 0.675 63 Q1_4 (.p4.) 0.756 0.041 18.308 0.000 0.700 0.478 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 64 65 Intercepts: 66 63 67 .Q1_1 4.748 0.034 140.596 0.000 4.748 3.482 68 .Q1_2 4.953 0.027 182.141 0.000 4.953 4.511 69 .Q1_3 4.728 0.032 148.704 0.000 4.728 3.683 70 .Q1_4 4.818 0.036 132.793 0.000 4.818 3.289

7.9 多母集団同時分析 71 f_A 0.000 0.000 0.000 72 73 Variances: 74 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 75 .Q1_1 1.555 0.059 26.460 0.000 1.555 0.836 76 .Q1_2 0.559 0.038 14.861 0.000 0.559 0.464 77 .Q1_3 0.896 0.049 18.158 0.000 0.896 0.544 78 .Q1_4 1.656 0.065 25.361 0.000 1.656 0.772 79 f_A 0.856 0.075 11.434 0.000 1.000 1.000 因子負荷のところに (.p1.) などといったものが表示されるようになりました。これは lavaan における等値制約の表し方です。Group 1 [1]: と Group 2 [2]: のそれぞれ対応 する因子負荷には同じ名前がついているため，等値制約が置かれたことがわかります。 また Intercepts: 欄には各観測変数の切片（＝グループ平均）が表示されています。これを 見る限りでは Group 2 [2]:，つまり女性の方が平均点が高そうですが，先ほど説明した理 由よりこれが「因子得点の差」なのか「それ以外の要因の差」なのかはまだわかりません。 そしてもう一点。Group 2 [2]: の因子の分散 (f_A) の値が 1 では無くなっています。一方 Group 1 [1]: のほうでは相変わらず因子の分散は 1 に固定されています，これは，因子負 荷に等値制約を置いたことで，共通因子のスケールは Group 1 [1]: 側で決まり，Group 2 [2]: 側では制約を置く必要がなくなったためです。その結果，因子負荷の Std.lv 列は「全 ての潜在変数の分散を 1 にした」時の推定値なので，グループ間で値が変わってしまいます。 せっかく等値制約を置いたのに，これでは多少違和感がありますが，そういうものなのであ まり気にしないでください。実際に弱測定不変モデルの結果を報告する際には，たぶん因子 負荷は Estimate 列の値を載せて，各グループの因子得点の分散を記入することが多いと思 います。 続いて，弱測定不変モデルと配置不変モデルの間で適合度が著しく悪化していないかを検証 します。比較の方法は7.7節と同じです。多母集団モデルでは，各パラメータをグループごと に異なるものとしてそれぞれ自由に推定していました。ここに等値制約を置くということは， 自由に推定できるパラメータの数が減少していることを意味するので，それだけ適合度は 悪くなるはずです。モデル修正におけるパスの削除が「そのパスの係数を0に固定する」だ (F) (M) とすると，多母集団モデルにおける等値制約は例えば「因子負荷𝑏𝑖 の値を𝑏𝑖 に固定する 」としているようなもので，やっていることは同じなわけですね。 1 summary(compareFit(result_base, result_weak)) 1 ################### Nested Model Comparison ######################### 2 3 Chi-Squared Difference Test 4 5 6 64 Df result_base AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq) 4 31636 31775 34.458

7.9 多母集団同時分析 7 result_weak 7 31631 31753 35.925 1.4678 0 3 0.6897 8 9 ####################### Model Fit Indices ########################### 10 chisq df pvalue 11 result_base 34.458† 12 result_weak 4 35.926 .000 13 7 rmsea .079 cfi .979 tli .936 .000 .058† .980† .965† srmr .022† .023 aic 31635.698 31631.166† bic 14 result_base 15 result_weak 31752.891† 31774.813 16 17 ################## Differences in Fit Indices ####################### 18 19 df result_weak - result_base rmsea cfi tli srmr aic bic 3 -0.021 0.001 0.029 0.001 -4.532 -21.922 今回の等値制約では，4 つの項目の因子負荷にそれぞれ等値制約を置いたため，パラメータ が 4 つ減少しています。一方で，先程説明したように Group 2 [2]: の因子の分散が自由推 定可能になったことでパラメータが実は 1 つ増えており，弱測定不変モデルでは差し引き 3 つのパラメータが減少しています。これが自由度 Df の違いに表れているのです。 結果を見ると，どうやら等値制約による当てはまりの悪化よりもパラメータを減らした（倹 約度の向上）効果の方が大きいらしく，RMSEA や CFI は改善しました。尤度比検定も有意 ではなかったということで，弱測定不変性が確認されたと言って良いでしょう。 続いて，強測定不変モデル（因子負荷＋切片）の検討に移ります。 強測定不変モデル 1 # group.equalに"intercepts"を追加 2 result_strong <- cfa(model, data = dat, group = "gender", std.lv = 3 summary(result_strong, standardized = TRUE) 1 lavaan 0.6.15 ended normally after 37 iterations TRUE, group.equal = c("loadings", "intercepts")) 2 3 Estimator 4 Optimization method 5 Number of model parameters 6 Number of equality constraints ML NLMINB 26 8 7 8 Number of observations per group: 9 1 802 10 2 1630 11 12 Model Test User Model: 13 14 65 Test statistic 48.938

7.9 多母集団同時分析 15 Degrees of freedom 16 P-value (Chi-square) 17 Test statistic for each group: 10 0.000 18 1 18.205 19 2 30.733 20 21 Parameter Estimates: 22 23 Standard errors 24 Information 25 Information saturated (h1) model Standard Expected Structured 26 27 28 Group 1 [1]: 29 30 Latent Variables: 31 Estimate Std.Err z-value P(>|z|) Std.lv Std.all (.p1.) 0.626 0.037 16.693 0.000 0.626 0.432 (.p2.) 0.882 0.039 22.885 0.000 0.882 0.701 Q1_3 (.p3.) 0.910 0.040 22.835 0.000 0.910 0.696 Q1_4 (.p4.) 0.752 0.040 18.705 0.000 0.752 0.509 32 f_A =~ 33 Q1_1 34 Q1_2 35 36 37 38 Intercepts: 39 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 40 .Q1_1 (.10.) 4.385 0.037 117.490 0.000 4.385 3.025 41 .Q1_2 (.11.) 4.504 0.040 112.782 0.000 4.504 3.582 42 .Q1_3 (.12.) 4.296 0.042 103.377 0.000 4.296 3.285 43 .Q1_4 (.13.) 4.437 0.041 108.640 0.000 4.437 3.004 44 f_A 0.000 0.000 0.000 45 46 Variances: 47 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 48 .Q1_1 1.710 0.092 18.510 0.000 1.710 0.814 49 .Q1_2 0.804 0.062 13.057 0.000 0.804 0.508 50 .Q1_3 0.883 0.066 13.277 0.000 0.883 0.516 51 .Q1_4 1.616 0.092 17.663 0.000 1.616 0.741 52 f_A 1.000 1.000 1.000 53 54 55 Group 2 [2]: 56 57 Latent Variables: 58 66 59 f_A =~ 60 Q1_1 (.p1.) Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.626 0.037 16.693 0.000 0.579 0.422

7.9 多母集団同時分析 61 Q1_2 (.p2.) 0.882 0.039 22.885 0.000 0.816 0.743 62 Q1_3 (.p3.) 0.910 0.040 22.835 0.000 0.842 0.659 63 Q1_4 (.p4.) 0.752 0.040 18.705 0.000 0.696 0.475 64 65 Intercepts: 66 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 67 .Q1_1 (.10.) 4.385 0.037 117.490 0.000 4.385 3.198 68 .Q1_2 (.11.) 4.504 0.040 112.782 0.000 4.504 4.099 69 .Q1_3 (.12.) 4.296 0.042 103.377 0.000 4.296 3.362 70 .Q1_4 (.13.) 4.437 0.041 108.640 0.000 4.437 3.028 71 f_A 0.504 0.053 9.503 0.000 0.545 0.545 72 73 Variances: 74 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 75 .Q1_1 1.545 0.059 26.271 0.000 1.545 0.822 76 .Q1_2 0.541 0.037 14.646 0.000 0.541 0.448 77 .Q1_3 0.924 0.048 19.413 0.000 0.924 0.566 78 .Q1_4 1.663 0.065 25.474 0.000 1.663 0.774 79 f_A 0.856 0.075 11.413 0.000 1.000 1.000 Group 2 [2]: の f_A の Intercepts: の値がゼロではなくなりました。切片まで固定した ことによって，因子得点の平均値が自由推定になった状態です。推定値が 0.504 ということ は，女性の方が各観測変数のグループ平均が高かったのは確かに因子得点が平均的に高いた めだ，といえそうです。一方で，もしも適合度がひどく悪化していたら「切片がグループ間 で同じというのはおかしい」ということになり，観測変数のグループ間の差は，むしろ切片 項に起因するものだと言えそうです。 ということで適合度を比較 1 # 3モデル以上でも一気に出来ます 2 summary(compareFit(result_base, result_weak, result_strong)) 1 ################### Nested Model Comparison ######################### 2 3 Chi-Squared Difference Test 4 5 Df BIC Chisq Chisq diff RMSEA Df diff result_base 4 31636 31775 34.458 7 result_weak 7 31631 31753 35.925 1.4678 0.000000 3 8 result_strong 10 31638 31743 48.938 13.0123 0.052389 3 9 67 AIC 6 Pr(>Chisq) 10 result_base 11 result_weak 0.68973 12 result_strong 0.00461 ** 13 --- 14 Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

7.9 多母集団同時分析 15 16 ####################### Model Fit Indices ########################### 17 chisq df pvalue 18 result_base 34.458† 19 result_weak 35.926 20 result_strong 48.938 21 4 rmsea .000 7 10 .079 .000 .058 .000 .057† cfi tli srmr aic .979 .936 .022† 31635.698 .980† .965 .023 31631.166† .973 .967† .030 31638.178 bic 22 result_base 31774.813 23 result_weak 31752.891 24 result_strong 31742.514† 25 26 ################## Differences in Fit Indices ####################### 27 df rmsea cfi tli srmr 28 result_weak - result_base 3 -0.021 29 result_strong - result_weak 3 -0.002 -0.007 0.002 0.006 aic bic 0.001 0.029 0.001 -4.532 -21.922 7.012 -10.377 尤度比検定は有意になった（＝当てはまりが有意に悪くなった）のですが，適合度指標はさ ほど悪化していません。というか RMSEA や TLI は更に改善しました。尤度比検定はサン プルサイズが多いとすぐ有意になってしまうこともありますし，BIC も強測定不変モデルを 支持しているため，ここは強測定不変性が認められたと結論付けて先に進んでみましょう。 あとはどこまで測定不変性が認められるかを試すだけです。とりあえず一気に残りの全ての モデルを推定してみます。 残りのモデルも一気にやってしまう 1 # 厳密な測定不変モデル 2 # 独自因子の分散共分散を制約に追加 3 # 独自因子の共分散を追加している場合は"residual.covariances"も入れる 4 result_strict <- cfa(model, data = dat, group = "gender", std.lv = TRUE, group.equal = c("loadings", "intercepts", "residuals")) 5 6 # （ほぼ完全な測定不変モデル） 7 # 共通因子の分散共分散を制約に追加 8 # 2因子以上の場合は"lv.covariances"も入れる 9 result_almost <- cfa(model, data = dat, group = "gender", std.lv = TRUE, group.equal = c("loadings", "intercepts", "residuals", "lv.variances")) 10 11 # （完全な測定不変モデル） 12 # 因子得点の平均を制約に追加 13 result_complete <- cfa(model, data = dat, group = "gender", std.lv = TRUE, group.equal = c("loadings", "intercepts", "residuals", "lv.variances", "means")) そして全部一気に適合度を比較してみます。 68

7.9 多母集団同時分析

1

summary(compareFit(result_base, result_weak, result_strong,
result_strict, result_almost, result_complete))

1

################### Nested Model Comparison #########################

2
3

Chi-Squared Difference Test

4
5

Df

AIC

BIC

Chisq Chisq diff

RMSEA Df diff

6

result_base

4 31636 31775

34.458

7

result_weak

7 31631 31753

35.925

1.468 0.000000

3

8

result_strong

10 31638 31743

48.938

13.012 0.052389

3

9

result_strict

14 31652 31733

70.463

21.526 0.060026

4

10

result_almost

15 31657 31732

77.395

6.931 0.069841

1

11

result_complete 16 31760 31830 183.159

105.765 0.293522

1

12

Pr(>Chisq)

13

result_base

14

result_weak

0.6897314

15

result_strong

0.0046101 **

16

result_strict

0.0002491 ***

17

result_almost

0.0084694 **

18

result_complete

< 2.2e-16 ***

19

---

20

Signif. codes:

0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

21
22

####################### Model Fit Indices ###########################

23

chisq df pvalue

rmsea

cfi

tli

srmr

24

result_base

34.458†

25

result_weak

35.926

7

26

result_strong

48.938

10

.000 .057†

27

result_strict

70.463

14

.000

.058

.960

.966

.036

28

result_almost

77.395

15

.000

.058

.956

.965

.048

29

result_complete 183.159

16

.000

.093

.883

.912

.094

30

4

.000

.079

.000

aic

.058

.979

.936

.022†

.980†

.965

.023

.973

.967†

.030

bic

31

result_base

31635.698

32

result_weak

31631.166†

33

result_strong

31638.178

34

result_strict

31651.703

31732.854

35

result_almost

31656.635

31731.989†

36

result_complete

31760.399

31829.957

31774.813
31752.891
31742.514

37
38

################## Differences in Fit Indices #######################

39

69

df

rmsea

cfi

srmr

aic

3 -0.021

0.001

0.029 0.001

-4.532

result_strong - result_weak

3 -0.002 -0.007

0.002 0.006

7.012

result_strict - result_strong

4

0.001 -0.012 -0.001 0.007

13.526

result_almost - result_strict

1

0.001 -0.004 -0.001 0.011

4.931

40

result_weak - result_base

41
42
43

tli

7.10 lavaan の記法

44

result_complete - result_almost

1

45

0.034 -0.074 -0.053 0.046 103.765
bic

46

result_weak - result_base

-21.922

47

result_strong - result_weak

-10.377

48

result_strict - result_strong

49

result_almost - result_strict

-0.865

50

result_complete - result_almost

97.968

-9.660

この結果を踏まえてどこまでの測定不変性を認めるかについては諸説ありそうです。基本的
には適合度指標が許容範囲内にある中で最も制約の強いモデルを選べば良いのではないか，
と思いますが，人によって許容範囲内であっても前のモデルから適合度指標が著しく低下し
ている場合にはアウトとする人もいます。

7.10 lavaan の記法
知っていると役に立つ日が来るかもしれない，lavaan のより高度な記法について少しだけ
紹介します。

7.10.1 パラメータの制約
尺度を作成する際には，探索的因子分析の結果 2 項目だけの因子が誕生してしまうことがあ
ります。2 項目だけの因子は，単体では結構危うい存在です。というのも，2 項目の場合推
定するパラメータは因子負荷と独自因子の分散がそれぞれ 2 つずつで計 4 つあるにも関わら
ず，データの分散共分散は 3 つしかありません。つまりそのままでは推定に失敗する，とい
うことです*20 。試しに 2 項目 1 因子のモデルをやってみましょう。
うまくいかない 2 項目 1 因子モデル
1

model <- "f_A =~ Q1_1 + Q1_2"

2

summary(cfa(model, data = dat, std.lv = TRUE), standardized = TRUE)

1

Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats =
lavsamplestats, : lavaan WARNING:

2

Could not compute standard errors! The information matrix could

3

not be inverted. This may be a symptom that the model is not

4

identified.

1

lavaan 0.6.15 ended normally after 9 iterations

2
3

Estimator

4

Optimization method

5

Number of model parameters

ML
NLMINB
4

6

*20

70

他の因子との因子間相関（両方向の矢印）があれば，2 項目でも問題なく推定できるようになります。

7.10 lavaan の記法

7

Number of observations

2432

8
9

Model Test User Model:

10
11

Test statistic

NA

12

Degrees of freedom

-1

13

P-value (Unknown)

NA

14
15

Parameter Estimates:

16
17

Standard errors

18

Information

19

Information saturated (h1) model

Standard
Expected
Structured

20
21

Latent Variables:

22

Estimate

Std.Err

Q1_1

1.010

Q1_2

23

f_A =~

24
25

z-value

P(>|z|)

Std.lv

Std.all

NA

1.010

0.719

0.575

NA

0.575

0.491

Estimate

Std.Err

Std.lv

Std.all

26
27

Variances:

28

z-value

P(>|z|)

29

.Q1_1

0.954

NA

0.954

0.483

30

.Q1_2

1.044

NA

1.044

0.759

31

f_A

1.000

1.000

1.000

出力を見ると Degrees of freedom が-1 になっており，検定も行われていません。一応因
子負荷は出ているのですが，実際のところこれが正しいという保証は全くありません。この
ような場合には，2 項目の因子負荷に等値制約をかけることで推定するパラメータ数を減ら
す事ができます。
等値制約をかける
1

model <- "f_A =~ b1 * Q1_1 + b1 * Q1_2"

モデルの修正のところで，右辺の変数名の前に 0* をかけることで共分散をゼロに固定する
方法を紹介しました。それと同じように，ここでは右辺の変数名の前に共通の b1 をかけて
います。
1

summary(cfa(model, data = dat, std.lv = TRUE), standardized = TRUE)

1

lavaan 0.6.15 ended normally after 11 iterations

2

71

3

Estimator

4

Optimization method

ML
NLMINB

7.10 lavaan の記法 5 Number of model parameters 4 6 Number of equality constraints 1 7 8 Number of observations 2432 9 10 Model Test User Model: 11 12 Test statistic 13 Degrees of freedom 0.000 0 14 15 Parameter Estimates: 16 17 Standard errors 18 Information 19 Information saturated (h1) model Standard Expected Structured 20 21 Latent Variables: 22 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 23 f_A =~ 24 Q1_1 (b1) 0.762 0.023 32.813 0.000 0.762 0.543 25 Q1_2 (b1) 0.762 0.023 32.813 0.000 0.762 0.650 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 26 27 Variances: 28 29 .Q1_1 1.394 0.051 27.466 0.000 1.394 0.706 30 .Q1_2 0.793 0.039 20.520 0.000 0.793 0.577 31 f_A 1.000 1.000 1.000 結果を見ると，因子負荷の欄に測定不変性の検証の時と同じように (b1) という表記が追加 されています。パラメータ数を減らしたことによって，Degrees of freedom が 0（つまり 飽和モデル）になり，晴れてパラメータの検定が行われるようになりました。 （意見）尺度得点を和得点にするなら因子負荷の等値制約が必要じゃないか 尺度得点として，因子得点ではなく和得点を使うケースが多いことはすでに説明した とおりです。このような場合でも，尺度作成論文では通常「1 項目が 1 因子にのみ負 荷している状態（完全な単純構造） 」で CFA を行い，モデル適合度を確認しています。 ですがよく考えると，尺度得点が和得点であるならば，CFA を行う際に因子負荷に は等値制約をかけたほうがよいのではないかという気がしています。 因子負荷の大きさは，その項目が因子得点の推定に与える影響の強さを表していま す。言い換えると，因子負荷はその項目の重みという見方ができるということです。 そう考えると，等値制約を置かない CFA によるモデル適合度の検証は「重み付け得 点を使用した場合の適合度」であり，もしも項目の重みが極端に異なる場合は，重み 付けない単純な和得点では著しく適合度が悪化する可能性があります。 したがって，「和得点でのモデル適合度」を検証するために，因子ごとに全ての因子 72

7.10 lavaan の記法 負荷に等値制約をおいた状態でのモデル適合度を評価すべきではないか，という気が するのですが実際にそのようにしている研究は見たことが無い気が…と思っていたら こんな論文がありました。→(McNeish and Wolf, 2020) 等値制約は他にも，多母集団同時分析で特定の回帰係数が同じかどうかを評価する際や，複 数時点で同一のモデルを測定する場合に共通のパスが同じ値かを検討する際にも使ったりし ます。多母集団同時分析に関して，例えば図7.17の回帰分析モデルを性別ごとに行い，「age が Q1_A におよぼす影響」が男女間で同じかを検討したい場合には，以下のようにすると良 いです。 部分的な等値制約 1 model <- " 2 Q1_A ~ c(b,b) * age + education 3 " 4 5 result <- sem(model, data = dat, group = "gender") 6 summary(result) 1 lavaan 0.6.15 ended normally after 29 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters 8 6 Number of equality constraints 1 ML NLMINB 7 8 Used Total 9 Number of observations per group: 1 732 802 10 2 1500 1630 11 12 Model Test User Model: 13 14 Test statistic 15 Degrees of freedom 16 P-value (Chi-square) 17 Test statistic for each group: 2.008 1 0.157 18 1 1.472 19 2 0.536 20 21 Parameter Estimates: 22 23 Standard errors Standard 24 Information Expected 25 Information saturated (h1) model 26 73 Structured

7.10 lavaan の記法

27
28

Group 1 [1]:

29
30

Regressions:

31
32

Q1_A ~

33

age

34

education

Estimate

Std.Err

z-value

P(>|z|)

0.055

0.009

6.316

0.000

-0.097

0.144

-0.671

0.502

Estimate

Std.Err

z-value

P(>|z|)

20.779

0.522

39.841

0.000

Estimate

Std.Err

z-value

P(>|z|)

21.509

1.124

19.131

0.000

Estimate

Std.Err

z-value

P(>|z|)

0.055

0.009

6.316

0.000

0.145

0.101

1.430

0.153

Estimate

Std.Err

z-value

P(>|z|)

22.035

0.383

57.560

0.000

Estimate

Std.Err

z-value

P(>|z|)

16.617

0.607

27.386

0.000

(b)

35
36

Intercepts:

37
38

.Q1_A

39
40

Variances:

41
42

.Q1_A

43
44
45

Group 2 [2]:

46
47

Regressions:

48
49

Q1_A ~

50

age

51

education

(b)

52
53

Intercepts:

54
55

.Q1_A

56
57

Variances:

58
59

.Q1_A

このように掛け算する時に group の数と同じ長さの c() を与えると，グループの順にそれ
ぞれパラメータに名前を付けることができます。もちろん c(0.5, 0.6)* というように直接
数字を入力することによって，グループごとに異なる値に固定することもできます。

7.10.2 等式・不等式制約
図7.19のように，2 番目の項目を逆の因子に入れてしまった因子分析モデルを考えます。
このモデルを推定すると，（設定をミスっているので）推定値が変なことになります。

74

1

# 間違えちゃった

2

model <- "

7.10 lavaan の記法

図7.19

間違えちゃった

3

f_A =~ Q1_1 + Q1_7 + Q1_3

4

f_E =~ Q1_6 + Q1_2 + Q1_8

5

"

6

result <- cfa(model, data = dat, std.lv = TRUE)

1

Warning in lav_object_post_check(object): lavaan WARNING: covariance
matrix of latent variables

2

is not positive definite;

3

use lavInspect(fit, "cov.lv") to investigate.

なにやら Warning が出てきました。どうやらメッセージを読む限りは「潜在変数の共分散
行列が正定値 (positive definite) ではない」ということですが，つまりどういうことでしょ
うか。
1

summary(result, standardized = TRUE)

1

lavaan 0.6.15 ended normally after 33 iterations

2
3

Estimator

4

Optimization method

5

Number of model parameters

ML
NLMINB
13

6
7

Number of observations

2432

8
9

Model Test User Model:

10
11

Test statistic

12

Degrees of freedom

13

P-value (Chi-square)

844.850
8
0.000

14
15

Parameter Estimates:

16
17

75

Standard errors

Standard

7.10 lavaan の記法 18 Information 19 Information saturated (h1) model Expected Structured 20 21 Latent Variables: 22 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 23 f_A =~ 24 Q1_1 0.500 0.033 15.207 0.000 0.500 0.356 25 Q1_7 0.407 0.030 13.355 0.000 0.407 0.309 26 Q1_3 0.726 0.034 21.515 0.000 0.726 0.556 27 f_E =~ 28 Q1_6 0.318 0.029 11.027 0.000 0.318 0.258 29 Q1_2 0.753 0.034 22.179 0.000 0.753 0.642 30 Q1_8 0.372 0.030 12.325 0.000 0.372 0.289 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 1.290 0.057 22.693 0.000 1.290 1.290 31 32 Covariances: 33 34 f_A ~~ 35 f_E 36 37 Variances: 38 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 39 .Q1_1 1.725 0.053 32.801 0.000 1.725 0.874 40 .Q1_7 1.566 0.047 33.520 0.000 1.566 0.904 41 .Q1_3 1.177 0.048 24.574 0.000 1.177 0.691 42 .Q1_6 1.415 0.042 33.927 0.000 1.415 0.933 43 .Q1_2 0.808 0.046 17.491 0.000 0.808 0.588 44 .Q1_8 1.519 0.045 33.623 0.000 1.519 0.916 45 f_A 1.000 1.000 1.000 46 f_E 1.000 1.000 1.000 Covariances: の推定値（35 行目）が 1.290 になっています。この値は，2 つの因子 f_A と f_E の相関係数を表しているはずなので 1 を超えるのは明らかにおかしな話です。このよう に，モデルの設定をめちゃめちゃ失敗すると不適解が発生してしまいます。不適解が生じた 場合には，まずモデルの誤設定を修正する（不要なパスを削除したり，変数そのものを削除 したり…）べきです。ただ場合によってはモデルの設定は問題ないという自負があるにも関 わらず不適解が生じることがあります。よくあるのはデータが悪い状況です。SEM では，デ ータにおいて観測された分散共分散をもとにして，モデル上の分散共分散を推定していきま す。したがって，データが酷い（サンプルサイズが小さい・欠測値が多すぎる・適当な回答が 見られる）ときには，モデルが正しくてもうまく推定できないことがあります*21 。Chapter 3 において適当な回答をしているデータを除外する方法を紹介したのはこのような事態を回 避するためでもあったのです。 そして特にモデルが複雑な場合やデータがひどい場合には，独自因子の分散がマイナスにな *21 R に慣れた方は，試しに「変数間の相関関係を完全に無視してデータのサイズだけ合わせたデータ」を乱数 生成して試してみてください。ほぼ確実に推定失敗すると思います。 76

7.10 lavaan の記法
るなどの事態（Heywood case）が起こることがあります。基本的には不適解はそのまま報
告してはいけないのですが，

• どう考えてもモデルは正しいし，
• データには問題ない（と思っている）し，
• 不適解の部分以外はいい感じの推定値だし，
• 不適解といってもギリギリアウト（独自因子の分散が-0.005 とか）
という感じで，自分は悪くないと言い切れる自信がある状況であれば，多少の不適解には目
をつぶって*22 不等式制約をかけてとりあえず結果を報告することがあります。そういう場合
には，パラメータに上限・下限を定めて推定するという方法をとります。
パラメータの上限・下限を決める
1

model <- "

2

f_A =~ Q1_1 + Q1_7 + Q1_3

3

f_E =~ Q1_6 + Q1_2 + Q1_8

4
5

# 制約をかけたいパラメータにラベルを付ける

6

f_A ~~ cor * f_E

7
8

# パラメータの制約をふつうに記述する

9

cor < 1

10
11

# 同様に，独自因子の分散を0以上にしたいなら

12

Q1_1 ~~ var_e1 * Q1_1

13

var_e1 > 0

14

"

15

result <- cfa(model, data = dat, std.lv = TRUE)

16

summary(result, standardized = TRUE)

1

lavaan 0.6.15 ended normally after 62 iterations

2
3

Estimator

4

Optimization method

5

Number of model parameters

6

Number of inequality constraints

ML
NLMINB
13
2

7
8

Number of observations

2432

9
10

Model Test User Model:

11
12

Test statistic

13

Degrees of freedom

*22

77

862.410
8

独自因子の分散の真値がほぼゼロの場合，標本誤差の影響などで 0 を下回ってしまうことが実際にあるので
す。そうだと言い切れるなら仕方ないのです。

7.10 lavaan の記法 14 P-value (Chi-square) 0.000 15 16 Parameter Estimates: 17 18 Standard errors 19 Information 20 Information saturated (h1) model Standard Expected Structured 21 22 Latent Variables: 23 Estimate Std.Err z-value P(>|z|) Std.lv Std.all Q1_1 0.585 0.033 17.759 0.000 0.585 0.416 Q1_7 0.315 0.031 10.041 0.000 0.315 0.239 27 Q1_3 0.843 0.032 26.209 0.000 0.843 0.645 28 f_E =~ 29 Q1_6 0.253 0.029 8.579 0.000 0.253 0.205 30 Q1_2 0.884 0.030 29.289 0.000 0.884 0.754 31 Q1_8 0.347 0.031 11.336 0.000 0.347 0.270 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 1.000 NA 1.000 1.000 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 24 f_A =~ 25 26 32 33 Covariances: 34 35 f_A ~~ 36 f_E (cor) 37 38 Variances: 39 40 .Q1_1 1.633 0.052 31.646 0.000 1.633 0.827 41 .Q1_7 1.632 0.048 33.979 0.000 1.632 0.943 42 .Q1_3 0.995 0.047 21.346 0.000 0.995 0.584 43 .Q1_6 1.453 0.042 34.230 0.000 1.453 0.958 44 .Q1_2 0.594 0.043 13.723 0.000 0.594 0.432 45 .Q1_8 1.537 0.046 33.715 0.000 1.537 0.927 46 f_A 1.000 1.000 1.000 47 f_E 1.000 1.000 1.000 (vr_1) 48 49 Constraints: 50 |Slack| 51 1 - (cor) 0.000 52 var_e1 - 0 1.633 ということで，Warning も出ずに無事推定できました。ただ，この方法はあくまでも対処療 法です。まずはモデルの誤設定の可能性を完全に消し去るのが第一です。 不等式制約と同じように等式制約や四則演算による制約を置くことも可能です。例えば（意 味があるかはともかく） 「項目 1 の因子負荷」が「項目 2 と項目 3 の因子負荷の平均」よりも 大きいということが仮説から言えるとしたら，次のような書き方が出来ます。 78

7.10 lavaan の記法 1 # 意味があるかはわからないが 2 model <- " 3 f_A =~ b1*Q1_1 + b2*Q1_2 + b3*Q1_3 + Q1_4 + Q1_5 4 b1 > (b2 + b3)/2 5 " 7.10.3 推定しないパラメータの定義（媒介分析） SEM でよくあるのが，ある変数が別の変数に及ぼす効果を「直接的な効果（直接効果）」と 「別の変数を介して及ぼす効果（間接効果）」に分解して考えるモデルです。図7.20では，学 歴に性別が及ぼす効果について，「性別が勤勉性 (Q1_C) を介して学歴に影響している」とい う間接効果を考えています。 図7.20 媒介モデル このように，ある変数間の関係に第 3 の変数を考える分析のことを媒介分析 (mediation analysis) と呼びます。Baron and Kenny (1986) の枠組みでは，gender から education に対する Q1_C の媒介効果を確認するためには • Q1_C をモデルから消したときには gender から education への単回帰が有意である • その上で gender から Q1_C への単回帰も Q1_C から age への単回帰も有意になる • Q1_C をモデルに含めた場合には，gender から education への直接効果が有意では なくなる ことを示せば良いとしています。理想的には 3 点目に関して gender から education への 直接効果が完全にゼロになってくれれば，「gender は Q1_C を通してのみ education に影 響を与える」ということがハッキリ言えるわけですが，実際にはなかなかそうもいきません。 とりあえず gender から education への直接効果が有意に減少してくれれば万々歳です。 図7.20のモデルでは，gender から education への経路は 2 つあります。こうしたケースで は 2 変数の共分散は「全ての経路のパス係数の積の和」でした。媒介分析でも同じように， 直接効果が 𝑐，間接効果が 𝑎𝑏 となり，これらを合わせた 𝑎𝑏 + 𝑐 を gender から education への総合効果と呼びます。先程の 1 点目は，言い換えると「gender から education への総 合効果が有意である」ということです。 したがって，先程の Baron and Kenny (1986) の枠組みを SEM のパラメータで表すと • 𝐻1 ∶ 𝑎𝑏 + 𝑐 ≠ 0 • 𝐻2 ∶ 𝑎 ≠ 0 79

7.10 lavaan の記法 • 𝐻3 ∶ 𝑏 ≠ 0 • 𝐻4 ∶ 𝑎𝑏 ≠ 0 を全て示すことができれば，Q1_C による媒介効果を認めることができるわけです*23 。 それでは lavaan で媒介分析を行ってみましょう。 媒介分析 1 model <- " 2 education ~ b*Q1_C + c*gender 3 Q1_C ~ a*gender 4 5 # 制約ではなく事後的な計算をして欲しい 6 ab := a*b 7 total := ab + c 8 " ":=" という記号を使うと，ラベル付けしたパラメータの推定値に基づいて計算をしてくれ ます。"ab := a*b" と書けば， 「a*b を ab と定義します」という意味であり，この記述によ って ab というパラメータを推定するということはありません。ただ単に推定された a と b の値をもとに，事後的に ab を計算して出力してくれるようになります。一方，ab = a*b と してしまうと，ab もパラメータとして推定されてしまいます。 1 result <- sem(model, data = dat) 2 summary(result, standardized = TRUE) 1 lavaan 0.6.15 ended normally after 14 iterations 2 3 Estimator 4 Optimization method 5 Number of model parameters ML NLMINB 5 6 7 8 Number of observations Used Total 2232 2432 9 10 Model Test User Model: 11 12 Test statistic 13 Degrees of freedom 0.000 0 14 15 Parameter Estimates: 16 17 *23 80 Standard errors Standard もしこれに加えて 𝐻5 ∶ 𝑐 = 0 が示せれば，完全な媒介効果だといえます。因果関係における中継点を見つけ た，という意味で因果関係のメカニズム解明にかなり役立つ情報になりえます。

7.11 同値モデル 18 Information 19 Information saturated (h1) model Expected Structured 20 21 Regressions: 22 23 24 25 26 27 Estimate Std.Err z-value P(>|z|) Std.lv Std.all education ~ Q1_C (b) 0.006 0.005 1.187 0.235 0.006 0.025 gender (c) 0.007 0.050 0.134 0.893 0.007 0.003 (a) 0.822 0.213 3.852 0.000 0.822 0.081 Estimate Std.Err z-value P(>|z|) Std.lv Std.all Q1_C ~ gender 28 29 Variances: 30 31 .education 32 .Q1_C 1.235 0.037 33.407 0.000 1.235 0.999 22.423 0.671 33.407 0.000 22.423 0.993 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 33 34 Defined Parameters: 35 36 ab 0.005 0.004 1.134 0.257 0.005 0.002 37 total 0.012 0.050 0.232 0.817 0.012 0.005 Defined Parameters: という箇所が増えました。定義されたパラメータについても，自由推 定されたパラメータと同じように仮説検定まで行ってくれます。ということで，これを使えば 先程の仮説群をまとめて検証することが出来ます。今回の場合では，そもそも 𝐻1 ∶ 𝑎𝑏 + 𝑐 ≠ 0 が有意ではないのでおしまいです*24 。 7.11 同値モデル 図7.21に 2 つの異なるモデルを考えました。上は「gender が education に及ぼす影響」を考 えるモデルです。これに対して下のモデルでは矢印の向きを全てひっくり返して「education が gender に及ぼす影響」を推定するモデルになっています。 図7.21 矢印の向きを変えてみた 普通に考えたら，少なくとも「学歴によって性別が影響を受ける」よりは「性別によって学 歴が影響を受ける」のほうが自然な気がしますが，とりあえずこれらのモデルを分析してみ ましょう。 *24 81 当然ですが，total の値は lm(education ~ gender,data=dat) の結果と一致します。

7.11 同値モデル 1 model1 <- " 2 education ~ Q1_C 3 Q1_C ~ Q1_A 4 Q1_A ~ gender 5 " 6 7 model2 <- " 8 gender ~ Q1_A 9 Q1_A ~ Q1_C 10 Q1_C ~ education 11 " 12 13 result1 <- sem(model1, data = dat) 14 result2 <- sem(model2, data = dat) 15 # 適合度指標を比較，ネストしていないので nested = FALSE をつけて 16 summary(compareFit(result1, result2, nested = FALSE)) 1 ####################### Model Fit Indices ########################### 2 chisq df pvalue 3 result1 5.916† 4 result2 5 5.916 3 rmsea cfi tli .116 .021† .988† .977† 3 .116 .021 .988 srmr .016 .977 aic 32827.430 .016† 28979.094† bic 6 result1 7 result2 29013.358† 32861.694 この 2 つのモデルでは，RMSEA, CFI, TLI, SRMR といった適合度指標は全て完全に同じ 値になります*25 。ダガー（†）がついているのは，数値計算上の細かいズレのためであり，こ れも小数点第 9 位くらいでやっと差が出る程度の違いしかありません。適合度指標が完全に 一致しているということは，モデル上の共分散行列（𝚺）に関しても全く同じ値になってい るため，これらの 2 つのモデルは統計的には優劣をつけることが出来ません。 こうしたモデルのことを同値モデル (equivalent model) と呼びます。基本的には，回帰を含 む SEM のモデルには多少なりとも同値モデルが存在しており，可能であればこうした同値 モデルに関しても議論をしておくべきだという意見もあります (Kang and Ahn, 2021) *26 。 多くの場合，同値モデルは図7.20のように，どこかしらの矢印の向きを反対にしたモデルで す。一部の人は SEM の矢印を「因果関係」かのように扱ってしまいますが，これは誤った 考え方です。統計的にはデータに対する当てはまり（＝モデルの「良さ」 ）が完全に同じモデ ルがあるということは，分析上は矢印の向きはどちらでも良いといえます。これは回帰分析 の 𝑥 と 𝑦 は入れ替えても問題ないのと同じです。したがって，分析者は統計以外の側面から， つまり先行研究などを踏まえてなるべく理論的に矢印の向きを指定する必要があります。 「同 値モデルに関しても議論をしておくべき」というのは「なぜ矢印の向きが逆ではいけないの *25 *26 82 AIC や BIC がぜんぜん違う値になっているのですが，この理由は分かりません…何か計算方法のせいだとは 思うのですが，とりあえず今回のケースでは AIC や BIC はあてにならないと思ってください。 ちなみに同値モデルの作り方に関しては Mayekawa (1994) などにあります。

7.12 𝜔ℎ を求める か」をきちんと説明しておく必要があるということです。 7.12 𝜔ℎ を求める Chapter 6 で紹介した 𝜔 係数のお話の続きです。今回は，分析者が事前に設定した因子構造 における 𝜔ℎ を求める方法を紹介します。 7.12.1 双因子モデル 図7.22の 10 項目２因子のモデルに関してやってみましょう。 図7.22 𝜔 係数を求める まずは今までと同じように cfa() でパラメータの推定を行います。 omega_h を求めるために 1 model <- " 2 g =~ Q1_1 + Q1_2 + Q1_3 + Q1_4 + Q1_5 + Q1_6 + Q1_7 + Q1_8 + Q1_9 + 3 f_A =~ Q1_1 + Q1_2 + Q1_3 + Q1_4 + Q1_5 4 f_C =~ Q1_6 + Q1_7 + Q1_8 + Q1_9 + Q1_10 Q1_10 5 6 # G因子と特殊因子は無相関にしないと推定できない 7 g ~~ 0*f_A + 0*f_C 8 " 9 10 result <- cfa(model, data = dat, ordered = TRUE) 𝜔 係数を算出するためには semTools::compRelSEM() という関数を使います*27 。 omega_h を求める 83 1 compRelSEM(result) *27 少し前までは semTools:reliability() という関数もあったのですが，すでに deprecated になっている ので semTools::compRelSEM() を使います。

7.12 𝜔ℎ を求める 1 2 g f_A f_C 0.566 0.633 0.034 出てくる値は，それぞれの因子のみの 𝜔 係数です。したがって，g のところにある 0.566 とい う値が，𝜔ℎ を表しています。第 10 回で出た値（0.42）と比べるとそこそこ大きな値になって います。これは「カテゴリカル因子分析にした」ことと「もう一方の因子に対する因子負荷を 完全にゼロにした」ことが理由だと思われます。g 以外の特殊因子の値は， 「G 因子の影響を 取り除いた 𝜔」という位置づけなので，使う機会は無いでしょう。また，引数 return.total = TRUE とすると，𝜔𝑡 も出してくれます。ただ第 10 回でお話したように 𝜔𝑡 を優先的に使う ことは無いと思うので，基本的にはそのままで OK です。 7.12.2 高次因子モデル G 因子の位置づけに関して，場合によっては図7.23のように各特殊因子の上位にあると考え るほうが妥当なケースがあるかもしれません（総合的な学力に対する「文系」 「理系」なんか はまさにそれ）。 図7.23 高次因子モデル このような場合に，全ての項目をまとめて一つの得点と考えた場合の内的一貫性は，𝜔ℎ とは 異なる形になります。項目 1 について式を展開すると， 𝐲1 = 𝐟𝐴 𝑏1 + 𝛆1 𝑔 = (𝐆𝑏𝐴 + 𝛆𝐴 )𝑏1 + 𝛆1 (7.51) error = 𝑔 𝐆𝑏𝐴 𝑏1 ⏞⏞ +⏞ 𝑏1⏞ 𝛆⏞ 𝐴 + 𝛆1 𝑔 𝑔 となるため，真値 𝑡 に相当する部分は 𝐆𝑏𝐴 𝑏1 となります。したがって，その係数 𝑏𝐴 𝑏1 を使 って，高次因子モデルにおける 𝜔ℎ𝑜 は 𝐼 𝜔ℎ𝑜 = 𝑔 (∑𝑖=1 𝑏⋅𝑔 𝑏𝑖 ) 2 𝜎𝑦2𝑡𝑜𝑡𝑎𝑙 (7.52) 𝑔 となります。𝑏⋅𝑔 は項目によって 𝑏𝐴 か 𝑏𝐶 になるものです。 数理的にはそういう扱いですが，lavaan で行う場合にはそれほど難しくはありません。ま 84

7.13 もう一つの SEM
ずはいつも通りモデルを組みます。
（2 因子だと解が求められなかったので 3 因子分用意して
います。）

omega_ho を求めるために
1

model <- "

2

f_A =~ Q1_1 + Q1_2 + Q1_3 + Q1_4 + Q1_5

3

f_C =~ Q1_6 + Q1_7 + Q1_8 + Q1_9 + Q1_10

4

f_E =~ Q1_11 + Q1_12 + Q1_13 + Q1_14 + Q1_15

5
6

# 高次因子の測定方程式

7

g =~ f_A + f_C + f_E

8

"

9
10

# ordered = TRUEにするとcompRelSEM()でエラーが出る

11

result <- cfa(model, data = dat)

あとは semTools::compRelSEM() に入れるだけです。

omega_ho を求める
1

1
2

compRelSEM(result, higher = "g")

f_A

f_C

f_E

g

0.714 0.745 0.767 0.638

引数 higher に高次因子の名前を入れてあげるだけです。先ほどと同じように g のところに
ある 0.638 という値が，𝜔ℎ𝑜 を表しているわけです。
このように，semTools::compRelSEM() を利用することで，
「完全な単純構造ではない場合
の 𝜔」や「独自因子の間に相関を認めた場合の 𝜔」などもフレキシブルに計算可能です。

7.13 もう一つの SEM
これまでに紹介してきた SEM は，共分散構造分析という別名にも表れているように，観測
変数間の共分散（相関係数）をうまく復元できるようにモデルのパラメータを推定する方法
です。この共分散に基づく SEM (Covariance Based SEM; CB-SEM) に対して，主にマー
ケティングの分野として（とはいえ今では様々な分野で使われている）
，異なる SEM が独自
に発展しているようです。ということで最後に，もう一つの SEM こと PLS-SEM (Partial

Least Squares SEM) を簡単に紹介します。
7.13.1 PLS-SEM の考え方
SEM は多変量解析の手法ですが，もともと多変量解析の手法としては重回帰分析や分散分
析，クラスター分析などのシンプルな手法が存在しています。20 世紀にはこうした手法によ
って（あるいはこうした手法に当てはめられるように RQ やデータを設計することで）多変
量解析の研究が行われていました。しかしこうした手法には，(a) 複雑なモデルを考えること

85

7.13 もう一つの SEM ができない，(b) 潜在変数を考えることができない（すべての変数は直接観測可能だとみな す），(c) 測定誤差というものを明示的に考えることができない，といった限界がありました (Haenlein and Kaplan, 2004)。これらの問題を解決するために考案されたのが (CB-)SEM ということのようです。そもそも SEM というものは「構造方程式モデリング」と訳される名 が示すように，潜在変数を含む複雑な変数間の関係性を（連立）方程式によって表現し，そ のパラメータをデータから推定する方法でした。PLS-SEM とは，その定式化および推定の 方法として，CB-SEM とは別の文脈から生まれてきたもののようです。 モデルの設計 したがって，PLS-SEM においても基本的なコンセプトは CB-SEM と大きく変わりません。 図7.24は，ホテル利用者が行う評価に及ぼす影響をイメージした仮想的なグラフィカルモデ ルです*28 。このモデルには，以下の 3 つの潜在変数およびそれらを測定する計 7 項目の観測 変数があるとします。 図7.24 満足度 仮想的なグラフィカルモデル ホテルのサービスに対する満足度。以下の 3 項目によって，それぞれ 7 段階評価で 測定される。 1. スタッフの対応は良かった。(sat_1) 2. ホテルは綺麗だった。(sat_2) 3. 食事は美味しかった。(sat_3) 価格 （主観的な）価格。 「1 泊あたりの価格はどう感じましたか。(price)」という 1 項目に 対して「とても安かった」から「とても高かった」までの 7 段階評価。 評価 ホテルの評価。以下の 3 項目によって，それぞれ 7 段階評価で測定される。 1. このホテルのファンになった。(rate_1) 2. 機会があればまたこのホテルに泊まりたい。(rate_2) *28 86 このモデルはあくまでも仮想です。一応元ネタはAlbers (2010) ですが，皆さんのお持ちの知見に照らし合 わせたときに「そんな因果関係なわけないだろ！」といった批判は甘んじて受け入れます。

7.13 もう一つの SEM 3. このホテルを他の人に薦めたい。(rate_3) このような図は，CB-SEM でも見ることがあるかと思いますが，いくつか気になる点があり ます。 1 点目は，3 つの観測変数から「満足度」という潜在変数に対して矢印が刺さっている点で す。これまでの CB-SEM の説明では，矢印は必ず潜在変数から観測変数に向かって刺さっ ていました。つまり逆向きです。実は PLS-SEM では，このような逆向きの矢印もモデルに 含めることができると言われています。 では，具体的にこの矢印の向きの違いは何を表しているのでしょうか。CB-SEM で用いられ てきた矢印の向きは， 「複数の観測変数の共通原因としての潜在変数が存在する」ということ を考えていました。ここでは，潜在変数を測定するために「観測変数に反映された共通部分」 を見出そうとします。そういった意味で，reflective measurement と呼ばれることがあり ます。rate_1-3 の項目を見るとわかりますが，多くの人では，これらの項目への回答は互 いに似たものになり，それは「ホテルに高い評価をした人ほどいずれも『当てはまる』寄り の回答をするだろう」という推測が立ちます。 一方で，図7.24に登場した逆向きの矢印は， 「複数の観測変数の総体として潜在変数が形作ら れる*29 」と考えます。そのため，こちらの矢印は formative measurement を表したものだ と呼ばれます。sat_1-3 の項目は，一部共通する側面はあるものの，共通の「原因」を持っ ていないと考えられます。つまり，「従業員の対応は良いがあまり綺麗ではないホテル」や 「ご飯は美味しいが従業員が無愛想なホテル」などが存在することは容易に想像でき，また 「（全体的な）満足度」というものはこうした個別の要素を組み合わせて決まってくるはずで す。「総合的に見て満足だと感じたから，ご飯は美味しかったし従業員は良くしてくれた…気 がする」というようにレビューをする人はいないでしょう。 reflective と formative の違いを表した図が，図7.25です。reflective measurement では，各 観測変数が尋ねている内容（構成概念の一側面）の共通部分（積集合）が潜在変数として表 れている，と考えます。そして各観測変数 (rate_1-3) は共通の原因変数（ 「評価」 ）の影響を 受けているため，基本的に相関係数は高いはずです。一方で formative measurement では， そもそも各観測変数が高い相関を持っている必要はありません。そして潜在変数は，観測変 数の和集合として規定されている，と考えているわけです。 ここまでの話は，どちらが良いとか悪いとかの話ではなく，あくまでも潜在変数の考え方の 違いの話です。どちらの考え方も本質的にはあり得るわけですが，CB-SEM では formative measurement は表現できない，という点において PLS-SEM を積極的に使う理由が生まれ るかもしれませんね。 図7.24のモデルの違和感の 2 点目は，「価格」という潜在変数が 1 つの観測変数のみと，し かも矢印ではなくただの線で繋がっている点です。これも CB-SEM では，モデルの識別の ために，各潜在変数に対して少なくとも 3 つの観測変数が必要でした。一方 PLS-SEM では， このように単一の観測変数をそのまま潜在変数として利用することができると言われていま *29 87 その点で，PLS-SEM は composite-based SEM と呼ばれることもあるようです (Hwang et al., 2020)。

7.13 もう一つの SEM 図7.25 ２つの測定方法の違い (Hair, 2017, Exhibit 2.7 をもとに作成) す。この点に関しては，どうやら計算上の都合によるところが大きいです。というのも後ほ ど説明しますが，PLS-SEM では「観測変数から潜在変数の得点を推定する（測定モデル）」 「潜在変数同士の関係（パス係数）を推定する（構造モデル）」という 2 つの計算を別々に行 い，これらを交互に繰り返します。PLS-SEM ではこのうち構造モデルにおいて，観測変数 と潜在変数の間のパス解析（回帰）を考えないために，観測変数を便宜的にそのまま潜在変 数として扱っているような気がします。したがって，単一観測変数から潜在変数の得点を推 定するといっても結局やっていることは恒等変換にすぎません。ちなみに CB-SEM の推定 では構造モデルと測定モデルを区別せずにまとめて同時に推定しているので，このような場 合はわざわざ「潜在変数に置き換える」ことはせずに，潜在変数のまま扱うことができます。 推定 続いて PLS-SEM のパラメータ推定のイメージを紹介します。その前に CB-SEM における パラメータ推定の基本理念をおさらいすると，CB-SEM では「推定されたパラメータから復 元された共分散行列」と「データにおける観測された共分散行列」のズレが最も小さくなる ようにパラメータを推定していました。これはいわば，データを最も良く説明できるパラメ ータの組を求める行為です。 これに対して PLS-SEM の基本理念は，内生潜在変数 (endogenous latent variable) を予 測した際の分散説明率を最大化するようにパラメータを推定していきます*30 。内生潜在変数 ）です。 とは，他の（外生）潜在変数からの矢印が刺さっている潜在変数（図7.24では「評価」 PLS-SEM では，最終的に内生潜在変数を被説明変数とした（重）回帰分析における分散説 明率＝決定係数（𝑅2 ）の値が最大になるように，以下に示した 4 つのステップを繰り返して 推定を行います。具体的な各ステップを見てみましょう (Henseler et al., 2012)。 【ステップ 1】 とりあえず各潜在変数の得点を観測変数から求める（図7.26）。この時点で は矢印の向きによらず，すべて観測変数の重み付け和として求める。初期値は全ての重みを 1 とする（＝単純な和得点から始める）。 *30 88 その点で，PLS-SEM は variance-based SEM と呼ばれることもあるようです。

7.13 もう一つの SEM 図7.26 推定のステップ 1 【ステップ 2】 各内生潜在変数を対応する外生潜在変数によって回帰してパス係数を求める （図7.27）。普通に最小二乗法で回帰してあげれば自ずと決定係数（𝑅2 ）が最大になるパス係 数を求めることができる。 図7.27 推定のステップ 2 【ステップ 3】 得られたパス係数を使って各内生潜在変数の値を再計算する（図7.28）。こ こでは，内生潜在変数の値を外生潜在変数の重み付け和として計算し，ステップ 1 で得られ た値から置き換える。外生潜在変数の値については多分ノータッチ。 図7.28 推定のステップ 3 【ステップ 4】 ステップ 1 で用いる潜在変数得点計算時の各観測変数の重みを以下の方法で 。 置き換える（図7.29） 89

7.13 もう一つの SEM a. reflective measurement の箇所は潜在変数得点と観測変数との相関係数に置き換える （mode A）。 b. formative measurement の箇所は観測変数一つ一つから潜在変数への単回帰を行った 回帰係数に置き換える（mode B）。 図7.29 推定のステップ 4 以上 4 つのステップを解が変わらなくなるまで繰り返すことで，最終的な推定値を得ます。 細かく見ると，PLS-SEM では図7.26‐7.29の中の点線で囲った部分ごとに計算をしているこ とになります。 7.13.2 CB-SEM と PLS-SEM の使い分け これまでに紹介してきた PLS-SEM の特徴を踏まえると，CB-SEM との使い分けのポイン トが見えてきます。ここでは，Hair et al. (2021, Table 1.6) をもとに，使い分けのポイント をいくつか見てみましょう。 • 分析の目的が説明（観測されたデータ全体に最もよく当てはまるモデルパラメータを 推定すること）よりも予測（外生潜在変数から内生潜在変数をできるだけ精度良く予 測・説明すること）に関心がある場合は PLS-SEM が良い。 • formative measure な潜在変数が含まれる場合は，CB-SEM では難しいので PLSSEM を使うしかない。 • サンプルサイズが小さい場合は PLS-SEM のほうが良い。というのも PLS-SEM は基 本的に回帰分析の組み合わせなので，CB-SEM と比べると少ないサンプルサイズで もそれなりに安定した解が出やすい，らしい。 • 観測変数の（多変量）正規性が怪しい場合は PLS-SEM のほうが良い。というのも CB-SEM は変数の分散共分散行列をもとに推定を行うが，その方法として，多変量正 規性を仮定した最尤法などが用いられるため，仮定が満たされない場合は厳しい。 PLS-SEM の利用に関するいくつかのレビュー論文によれば，PLS-SEM を使う理由として 最も多く挙げられているのが「サンプルサイズが小さい」 ，そして「データの非正規性」 「予測 が目的」 「モデルが複雑」などが続いているようです (e.g., Bayonne et al., 2020; Lin et al., 2020; Magno et al., 2022)。 しかし PLS-SEM にはそれなりの反対意見もあるようで，一部の人は「PLS-SEM はくそく 90

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について
らえだ」という感じもあります (e.g., Rönkkö et al., 2015)。ということで，理論（研究仮
説）的に PLS-SEM のほうが適切だと言える場合には PLS-SEM を使えばよいのですが，ど
ちらを使えば良いかわからない（どちらでもいける）場合には，とりあえず CB-SEM を試
してみて，解が得られないなど特殊な事情がある場合に初めて PLS-SEM の利用を考えるく
らいで良いのかもしれません*31 。

7.13.3 R で PLS-SEM
PLS-SEM を実行するためのソフトウェアとして，SmartPLSというものがよく用いられて
いるようです。一応 R でも PLS-SEM を行うためのパッケージとして plssem や seminr

(SEM in R) パッケージというものが作成されているので，関心がある方はそちらのパッケ
ージを見てみてください。ちなみにHair et al. (2021) は seminr の使い方を紹介している

OA 本です。
7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について
2023 年 1 月にリリースされたバージョン 0.6-13 以降では，lavaan パッケージで探索的因子
分析が可能になっていることが判明しました。ただの探索的因子分析なら psych::fa() で
十分なのですが，lavaan では SEM の枠組みと組み合わせることで EFA と CFA や回帰分
析を組み合わせることができる点が，もしかしたら便利かもしれません。こうしたモデルは
特に Exploratory SEM (ESEM; Asparouhov and Muthén, 2009) と呼ばれている比較的新
しい手法です。世間の注目を集めているかはわかりません。ということでこのおまけセクシ
ョンでは，lavaan で探索的因子分析および ESEM を行う方法を紹介していきます。

7.14.1 lavaan で EFA を試してみる
まずは普通に探索的因子分析を試してみます。6.4 節の 10 項目 2 因子モデル（図 6.10）を再
現してみましょう。lavaan で EFA を指定する方法は大きく分けると 2 種類です。

efa() 関数
シンプルな方法は，efa() 関数を利用するものです。

efa() 関数
1

# 使い方はpsych::fa()とほぼ同じ

2

res_efa <- efa(dat[, paste0("Q1_", 1:10)], nfactors = 2, rotation =
"oblimin")

3

summary(res_efa)

1

This is lavaan 0.6.15 -- running exploratory factor analysis

2

*31

91

とはいえ解が得られない理由がサンプルサイズが小さいせいで，努力で増やせるようなもの（企業データな
どそもそも集めにくいものでないの）であれば，PLS-SEM を使うことよりもサンプルサイズを増やすこと
を先に検討すべきだと思いますよ。

• https://www.smartpls.com

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について 3 Estimator 4 Rotation method 5 Oblimin gamma 6 Rotation algorithm (rstarts) 7 Standardized metric TRUE 8 Row weights None Number of observations 2432 ML OBLIMIN OBLIQUE 0 GPA (30) 9 10 11 12 Fit measures: 13 aic 14 bic sabic chisq df pvalue nfactors = 2 78530.51 78698.61 78606.47 375.946 26 cfi rmsea 0 0.932 0.074 15 16 Eigenvalues correlation matrix: 17 18 ev1 ev2 ev3 ev4 ev5 ev6 ev7 ev8 ev9 19 3.078 1.834 0.896 0.823 0.722 0.684 0.530 0.523 0.469 20 ev10 21 0.440 22 23 Standardized loadings: (* = significant at 1% level) 24 25 f1 f2 26 Q1_1 0.408* 27 Q1_2 0.661* 28 Q1_3 0.779* * 29 Q1_4 0.464* .* 30 Q1_5 0.624* 31 Q1_6 32 Q1_7 33 34 35 unique.var communalities 0.846 0.154 0.558 0.442 0.408 0.592 * 0.723 0.277 0.603 0.397 0.567* 0.690 0.310 0.623* 0.615 0.385 Q1_8 0.551* 0.689 0.311 Q1_9 0.691* 0.527 0.473 Q1_10 0.587* 0.634 0.366 36 37 f2 f1 total 38 Sum of sq (obliq) loadings 1.876 1.832 3.707 39 Proportion of total 0.506 0.494 1.000 40 Proportion var 0.188 0.183 0.371 41 Cumulative var 0.188 0.371 0.371 42 43 Factor correlations: (* = significant at 1% level) 44 45 f1 46 f1 1.000 47 f2 0.317* f2 1.000 回転方法を指定する引数は rotation です。psych::fa() では rotate だったのでご注意 92

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について ください。ということで，かなり簡単に lavaan の出力形式で EFA の結果を出すことがで きました。lavaan::efa() の利点の一つは，自動的に因子負荷の検定を行ってくれるところ です。SEM の枠組みにしたがって，回転後の因子負荷が0であるという帰無仮説のもとで検 定をしています。デフォルトでは 1% 水準で有意のところに * がつくようになっていますが， これも引数で調整可能です。ついでに因子負荷が小さいところも全て表示されるように書き 換えてみましょう。 1 # cutoff: 指定した値より小さい因子負荷を表示しない 2 # alpha.level: (*)をつける有意水準を指定する 3 summary(res_efa, cutoff = 0, alpha.level = 0.05) 1 This is lavaan 0.6.15 -- running exploratory factor analysis 2 3 Estimator 4 Rotation method 5 Oblimin gamma 6 Rotation algorithm (rstarts) 7 Standardized metric TRUE 8 Row weights None Number of observations 2432 ML OBLIMIN OBLIQUE 0 GPA (30) 9 10 11 12 Fit measures: 13 14 aic bic sabic chisq df pvalue nfactors = 2 78530.51 78698.61 78606.47 375.946 26 cfi rmsea 0 0.932 0.074 15 16 Eigenvalues correlation matrix: 17 18 ev1 ev2 ev3 ev4 ev5 ev6 ev7 ev8 ev9 19 3.078 1.834 0.896 0.823 0.722 0.684 0.530 0.523 0.469 20 ev10 21 0.440 22 23 Standardized loadings: (* = significant at 5% level) 24 25 93 f1 f2 unique.var communalities 26 Q1_1 0.408* -0.069* 0.846 0.154 27 Q1_2 0.661* 0.558 0.442 28 Q1_3 0.779* -0.032* 0.408 0.592 29 Q1_4 0.464* 0.142* 0.723 0.277 30 Q1_5 0.624* 0.019 0.603 0.397 31 Q1_6 -0.036 0.567* 0.690 0.310 32 Q1_7 -0.008 0.623* 0.615 0.385 33 Q1_8 0.022 0.551* 0.689 0.311 34 Q1_9 -0.012 0.691* 0.527 0.473 0.011

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について 35 Q1_10 0.049* 0.587* 0.634 0.366 36 37 f2 f1 total 38 Sum of sq (obliq) loadings 1.876 1.832 3.707 39 Proportion of total 0.506 0.494 1.000 40 Proportion var 0.188 0.183 0.371 41 Cumulative var 0.188 0.371 0.371 42 43 Factor correlations: (* = significant at 5% level) 44 45 f1 46 f1 1.000 47 f2 0.317* f2 1.000 ちなみに lavaan のデフォルトでは，最尤法によってパラメータの推定を行っています。した がって，上記の結果は psych::fa() で fm = "ml" を指定したときに同じになるはずです。 psych::fa() 関数と比べる 1 res_fa <- psych::fa(dat[, paste0("Q1_", 1:10)], nfactors = 2, fm = "ml") 2 # 因子負荷だけ表示する 3 res_fa$loadings 1 2 Loadings: 3 ML2 ML1 4 Q1_1 0.408 5 Q1_2 0.661 6 Q1_3 0.779 7 Q1_4 8 Q1_5 0.142 0.464 0.624 9 Q1_6 0.567 10 Q1_7 0.623 11 Q1_8 0.551 12 Q1_9 0.691 13 Q1_10 0.587 14 15 ML2 ML1 16 SS loadings 17 Proportion Var 0.186 0.182 18 Cumulative Var 0.186 0.368 1.863 1.819 分散説明率のところには細かい違いはあるのですが，因子負荷に関しては確かに同じ結果が 得られました。 lavaan::efa() 関 数 の 別 の 利 点 は，複数の因子数を一気に実行できるこ と で す。 引 数 nfactors に整数のベクトルを与えてあげると，すべての因子数で実行した結果をまと 94

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について めて一つのオブジェクトで返してくれます。 色々な因子数を試す 1 # 一個ずつ順番にやっているだけなので指定した数が多いほど時間はかかる 2 # ちなみにrotationのデフォルトは"geomin" 3 res_efa <- efa(dat[, paste0("Q1_", 1:10)], nfactors = 1:3) 4 summary(res_efa) 1 This is lavaan 0.6.15 -- running exploratory factor analysis 2 3 Estimator 4 Rotation method 5 Geomin epsilon 6 Rotation algorithm (rstarts) 7 Standardized metric TRUE 8 Row weights None Number of observations 2432 ML GEOMIN OBLIQUE 0.001 GPA (30) 9 10 11 12 Overview models: 13 aic bic sabic chisq df pvalue cfi rmsea 14 nfactors = 1 80259.06 80374.99 80311.44 2122.495 35 0 0.595 0.157 15 nfactors = 2 78530.51 78698.61 78606.47 375.946 26 0 0.932 0.074 16 nfactors = 3 78387.85 78602.32 78484.76 217.286 18 0 0.961 0.067 17 18 Eigenvalues correlation matrix: 19 20 ev1 ev2 ev3 ev4 ev5 ev6 ev7 ev8 ev9 21 3.078 1.834 0.896 0.823 0.722 0.684 0.530 0.523 0.469 22 ev10 23 0.440 24 25 Number of factors: 1 26 27 Standardized loadings: (* = significant at 1% level) 28 29 95 unique.var communalities 30 Q1_1 f1 .* 0.934 0.066 31 Q1_2 0.502* 0.748 0.252 32 Q1_3 0.523* 0.727 0.273 33 Q1_4 0.487* 0.763 0.237 34 Q1_5 0.480* 0.769 0.231 35 Q1_6 0.439* 0.808 0.192 36 Q1_7 0.507* 0.742 0.258 37 Q1_8 0.480* 0.769 0.231 38 Q1_9 0.552* 0.696 0.304

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について 39 Q1_10 0.531* 0.718 0.282 40 41 f1 42 Sum of squared loadings 2.327 43 Proportion of total 1.000 44 Proportion var 0.233 45 Cumulative var 0.233 46 47 Number of factors: 2 48 49 Standardized loadings: (* = significant at 1% level) 50 51 f1 52 Q1_1 0.407* 53 Q1_2 0.660* 54 Q1_3 0.778* 55 Q1_4 0.464* 56 Q1_5 0.623* 57 Q1_6 58 Q1_7 59 60 61 f2 unique.var communalities * 0.846 0.154 0.558 0.442 0.408 0.592 0.723 0.277 0.603 0.397 0.567* 0.690 0.310 0.622* 0.615 0.385 Q1_8 0.550* 0.689 0.311 Q1_9 0.690* 0.527 0.473 Q1_10 0.587* 0.634 0.366 .* 62 63 f2 f1 total 64 Sum of sq (obliq) loadings 1.876 1.831 3.707 65 Proportion of total 0.506 0.494 1.000 66 Proportion var 0.188 0.183 0.371 67 Cumulative var 0.188 0.371 0.371 68 69 Factor correlations: (* = significant at 1% level) 70 71 f1 72 f1 1.000 73 f2 0.309* f2 1.000 74 75 Number of factors: 3 76 77 Standardized loadings: (* = significant at 1% level) 78 79 96 f1 f2 f3 80 Q1_1 0.408* .* 81 Q1_2 0.671* * 82 Q1_3 0.796* * 83 Q1_4 0.475* .* 84 Q1_5 0.636* unique.var communalities 0.798 0.202 0.555 0.445 0.400 0.600 0.725 0.275 0.601 0.399

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について 85 Q1_6 0.379* 0.448* 0.637 0.363 86 Q1_7 0.472* 0.492* 0.496 0.504 87 Q1_8 .* 0.472* 0.709 0.291 88 Q1_9 0.793* 0.386 0.614 89 Q1_10 0.600* 0.614 0.386 90 91 f1 f3 f2 total 92 Sum of sq (obliq) loadings 1.873 1.714 0.491 4.079 93 Proportion of total 0.459 0.420 0.120 1.000 94 Proportion var 0.187 0.171 0.049 0.408 95 Cumulative var 0.187 0.359 0.408 0.408 96 97 Factor correlations: (* = significant at 1% level) 98 99 f1 f2 100 f1 1.000 101 f2 0.142 1.000 102 f3 0.325* 0.076 f3 1.000 出 力 さ れ る オ ブ ジ ェ ク ト は， あ く ま で も lavaan の 形 式 に 沿 っ て い る た め， 例 え ば fitMeasures() 関数を利用して（SEM 的な）適合度指標を出して，因子数を決める参考に することも可能だったりします。 EFA で SEM 的適合度指標 1 fitMeasures(res_efa) 1 97 2 npar 3 fmin 4 chisq 5 df 6 pvalue 7 baseline.chisq 8 baseline.df 9 nfct=1 nfct=2 nfct=3 20.000 29.000 37.000 0.436 0.077 0.045 2122.495 375.946 217.286 35.000 26.000 18.000 0.000 0.000 0.000 5201.585 5201.585 5201.585 45.000 45.000 45.000 baseline.pvalue 0.000 0.000 0.000 10 cfi 0.595 0.932 0.961 11 tli 0.480 0.883 0.903 12 nnfi 0.480 0.883 0.903 13 rfi 0.475 0.875 0.896 14 nfi 0.592 0.928 0.958 15 pnfi 0.460 0.536 0.383 16 ifi 0.596 0.932 0.962 17 rni 0.595 0.932 0.961 18 logl -40109.529 -39236.255 -39156.925 19 unrestricted.logl -39048.282 -39048.282 -39048.282

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について 20 aic 80259.058 78530.510 78387.849 21 bic 80374.987 78698.607 78602.319 22 ntotal 2432.000 2432.000 2432.000 23 bic2 80311.443 78606.467 78484.761 24 rmsea 0.157 0.074 0.067 25 rmsea.ci.lower 0.151 0.068 0.060 26 rmsea.ci.upper 0.162 0.081 0.076 27 rmsea.ci.level 0.900 0.900 0.900 28 rmsea.pvalue 0.000 0.000 0.000 29 rmsea.close.h0 0.050 0.050 0.050 30 rmsea.notclose.pvalue 1.000 0.087 0.006 31 rmsea.notclose.h0 0.080 0.080 0.080 32 rmr 0.201 0.065 0.049 33 rmr_nomean 0.201 0.065 0.049 34 srmr 0.115 0.034 0.026 35 srmr_bentler 0.115 0.034 0.026 36 srmr_bentler_nomean 0.115 0.034 0.026 37 crmr 0.127 0.038 0.028 38 crmr_nomean 0.127 0.038 0.028 39 srmr_mplus 0.115 0.034 0.026 40 srmr_mplus_nomean 0.115 0.034 0.026 41 cn_05 58.064 252.548 324.123 42 cn_01 66.704 296.256 390.563 43 gfi 0.800 0.970 0.982 44 agfi 0.685 0.936 0.946 45 pgfi 0.509 0.458 0.322 46 mfi 0.651 0.931 0.960 47 ecvi 0.889 0.178 0.120 今回の結果では，どうやら 2 因子モデルよりも 3 因子モデルのほうがモデル適合度は高いよ うです。ただこのあたりは因子分析の Chapter でお話したように，最終的には解釈可能性ま でを視野に入れて決めるようにしましょう。 モデル式で指定する lavaan で EFA を実行するもう一つの方法は， 「モデル式で EFA であることを指定する」と いうものです。これまでに本章では，変数名の前に 0* をつけることでその値を固定したり （7.7.1項），名前をつけることで等値制約を置く方法（7.10節）を紹介しました。これと同じ 要領で，因子名の前に「この因子にどの項目がかかるかはEFA的に決めてください」という ことを指定する記号を書き足すことができるのです。 モデル式で EFA を指定する 98 1 # 変数名が多いとしんどいね… 2 model_efa <- ' 3 efa("f12")*f1 =~

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について

4
5

Q1_1 + Q1_2 + Q1_3 + Q1_4 + Q1_5 + Q1_6 + Q1_7 + Q1_8 + Q1_9 + Q1_10
efa("f12")*f2 =~

6
7

Q1_1 + Q1_2 + Q1_3 + Q1_4 + Q1_5 + Q1_6 + Q1_7 + Q1_8 + Q1_9 + Q1_10
'

EFA 的に推定＆回転をしてほしい因子名の前に efa() という記号をかけてあげます。カッ
コの中は特に指定はありません（後で説明）。とりあえずここでは F1 と F2 の間の EFA で
あることがわかりやすいように"f12" としましたが，例えば efa("xxx")*f1 =~ などとし
ても（f1 と f2 に同じ名前がついていれば）OK です。
1

result_efa <- cfa(model_efa, data = dat, rotation = "oblimin")

2

summary(result_efa, standardized = TRUE)

1

lavaan 0.6.15 ended normally after 1 iteration

2
3

Estimator

4

Optimization method

5

Number of model parameters

ML
NLMINB
29

6
7

Rotation method

8

Oblimin gamma

9

Rotation algorithm (rstarts)

OBLIMIN OBLIQUE
0
GPA (30)

10

Standardized metric

TRUE

11

Row weights

None

Number of observations

2432

12
13
14
15

Model Test User Model:

16
17

Test statistic

18

Degrees of freedom

19

P-value (Chi-square)

375.946
26
0.000

20
21

Parameter Estimates:

22
23

Standard errors

24

Information

25

Information saturated (h1) model

Standard
Expected
Structured

26
27

Latent Variables:

28

99

Estimate

Std.Err

z-value

P(>|z|)

Std.lv

Std.all

Q1_1

0.574

0.034

17.024

0.000

0.574

0.408

Q1_2

0.775

0.026

29.811

0.000

0.775

0.661

Q1_3

1.017

0.028

36.231

0.000

1.017

0.779

29

f1 =~ f12

30
31
32

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について 33 Q1_4 34 Q1_5 0.788 0.028 27.908 0.000 0.788 0.624 35 Q1_6 -0.044 0.024 -1.878 0.060 -0.044 -0.036 36 Q1_7 -0.011 0.023 -0.459 0.647 -0.011 -0.008 37 Q1_8 0.029 0.025 1.133 0.257 0.029 0.022 38 Q1_9 -0.016 0.020 -0.795 0.427 -0.016 -0.012 39 Q1_10 0.079 0.031 2.592 0.010 0.079 0.049 40 f2 =~ f12 41 Q1_1 -0.097 0.033 -2.981 0.003 -0.097 -0.069 Q1_2 0.013 0.020 0.682 0.495 0.013 0.011 42 43 0.686 0.034 20.081 0.000 0.686 0.464 [ reached getOption("max.print") -- omitted 8 rows ] 44 45 Covariances: 46 47 f1 ~~ 48 f2 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.317 0.025 12.921 0.000 0.317 0.317 49 50 Variances: 51 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 52 .Q1_1 1.672 0.051 32.747 0.000 1.672 0.846 53 .Q1_2 0.767 0.029 26.155 0.000 0.767 0.558 54 .Q1_3 0.695 0.037 18.788 0.000 0.695 0.408 55 .Q1_4 1.585 0.050 31.457 0.000 1.585 0.723 56 .Q1_5 0.962 0.035 27.799 0.000 0.962 0.603 57 .Q1_6 1.046 0.036 29.327 0.000 1.046 0.690 58 .Q1_7 1.064 0.039 27.338 0.000 1.064 0.615 59 .Q1_8 1.142 0.038 29.675 0.000 1.142 0.689 60 .Q1_9 0.998 0.041 24.060 0.000 0.998 0.527 61 .Q1_10 1.687 0.059 28.395 0.000 1.687 0.634 62 f1 1.000 1.000 1.000 63 f2 1.000 1.000 1.000 だいぶ lavaan の見慣れた形で出力されてくれました。結果はもちろん lavaan::efa() の ときと同じ値なので説明は省略します。 実際のところ，EFA をするだけであれば efa() 関数を使っておけば良いと思います*32 。モ デル式での指定方法は，このあとの ESEM で必要になってくるのです…。 7.14.2 Exploratory SEM (ESEM) ESEM は，SEM と EFA を組み合わせたようなものです。この背景には，そもそも CFA を 探索的に使っているケースが割と存在する，ということがあるようです (Browne, 2001)。例 えば7.7節で紹介したように，モデルの適合度を改善するためにある場所にパスを足したり引 いたりすることは，まさに explotaroty といって良いでしょう。また，自分で作成した尺度 *32 100 もっと言えば psych::fa() でも何一つ問題はないんですけどね。

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について や，手元のデータで因子構造が再現できる怪しいデータを用いて SEM を実行する場合には， 通常まず EFA で因子構造・因子得点を出した上で，これを SEM の中に投入することが多く あります。ならばもういっそ全部ひっくるめて一気にやってしまおう，というモチベーショ ンもあるのかもしれません。 ということで ESEM を試してみましょう。ここでは図7.30のようなモデルを考えてみます。 左側の 6 項目と右側の 6 項目は，それぞれ異なる尺度を表しています（以後，左側を尺度 A， 右側を尺度 B と呼びます）。想定する因子数だけは決まっていて，それぞれが 2 因子を測定 する尺度です。さらに，尺度 A の一方の因子が尺度 B の一方の因子に影響を与える，という 因子間の関係性（回帰）も仮説として検証したいとします。 図7.30 ESEM モデル このモデルを SEM で実行する場合には，尺度 A・B の中の各項目がそれぞれどちらの因子 に所属するかを明確に決める必要がありました。図7.30で言えば，青い点線のところにパスを 許容しない（＝因子負荷は厳密にゼロとする）のが SEM でした。一方 ESEM では，この青 い点線のパス (cross-loadings) を許容します*33 。つまり，尺度 A に含まれる 2 因子（𝑓1 , 𝑓2 ） は 6 項目から EFA によって推定し，尺度 B の 6 項目からも同様に 2 因子（𝑓3 , 𝑓4 ）を推定 する，ということです。そして，小さい因子負荷もそのまま残して SEM の適合度指標の計 算などを行います。言い換えると，ESEM は SEM の制約を緩めたモデルあるいは EFA に 部分的な制約を加えたモデルという位置づけになっているわけです。 では実際に図7.30の ESEM モデルを書き下してみましょう。先程出てきた efa("")* の記 法を使って書いていきます。 ESEM のモデル式 101 1 model_esem <- ' 2 # 尺度AのEFA *33 もちろん，そもそもどこのパスが黒か青かわからない（完全に探索的な）状態でも ESEM は実行可能です。 というか ESEM の推定ではパスの色は区別していない，はず。

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について

3

efa("f12")*f1 =~ Q1_1 + Q1_2 + Q1_3 + Q1_6 + Q1_7 + Q1_8

4

efa("f12")*f2 =~ Q1_1 + Q1_2 + Q1_3 + Q1_6 + Q1_7 + Q1_8

5

# 尺度BのEFA

6

efa("f34")*f3 =~ Q1_11 + Q1_12 + Q1_13 + Q1_16 + Q1_17 + Q1_18

7

efa("f34")*f4 =~ Q1_11 + Q1_12 + Q1_13 + Q1_16 + Q1_17 + Q1_18

8
9

# 因子間の回帰

10

f3 ~ f1

11

f4 ~ f2

12

'

尺度 A と尺度 B はそれぞれ異なる因子を測定しているため，例えば Q1_1 が f3 や f4 にか
かることはありません（因子負荷は 0）。そして因子の回転は，「f1 と f2 が単純構造に近づ
くように」「f3 と f4 が単純構造に近づくように」を独立に実行していきます。このことを
表すために，f1 と f2 の前には efa("f12")*，f3 と f4 の前には efa("f34")* というよう
に，異なる名前の efa をかけているのです。
あとは推定するだけです。
1

# EFA部分があるのでrotationを指定可能

2

result_esem <- sem(model_esem, data = dat, rotation = "oblimin")

3

summary(result_esem, standardized = TRUE)

1

lavaan 0.6.15 ended normally after 29 iterations

2
3

Estimator

4

Optimization method

5

Number of model parameters

ML
NLMINB
36

6
7

Rotation method

8

Oblimin gamma

9

Rotation algorithm (rstarts)

OBLIMIN OBLIQUE
0
GPA (30)

10

Standardized metric

TRUE

11

Row weights

None

Number of observations

2432

12
13
14
15

Model Test User Model:

16
17

Test statistic

18

Degrees of freedom

19

P-value (Chi-square)

20
21
22

102

Parameter Estimates:

422.462
42
0.000

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について 23 Standard errors 24 Information 25 Information saturated (h1) model Standard Expected Structured 26 27 Latent Variables: 28 Estimate Std.Err z-value P(>|z|) Std.lv Std.all Q1_1 0.607 0.034 17.755 0.000 0.607 0.432 Q1_2 0.856 0.029 29.726 0.000 0.856 0.730 29 f1 =~ f12 30 31 32 Q1_3 0.909 0.032 28.397 0.000 0.909 0.696 33 Q1_6 -0.021 0.021 -0.996 0.319 -0.021 -0.017 34 Q1_7 -0.017 0.016 -1.039 0.299 -0.017 -0.013 35 Q1_8 0.099 0.028 3.559 0.000 0.099 0.077 36 f2 =~ f12 37 Q1_1 -0.067 0.033 -2.028 0.043 -0.067 -0.048 38 Q1_2 0.003 0.015 0.179 0.858 0.003 0.002 39 Q1_3 0.013 0.019 0.660 0.509 0.013 0.010 40 Q1_6 0.763 0.034 22.180 0.000 0.763 0.620 41 Q1_7 0.930 0.038 24.354 0.000 0.930 0.707 42 Q1_8 0.634 0.033 19.016 0.000 0.634 0.493 43 [ reached getOption("max.print") -- omitted 14 rows ] 44 45 Regressions: 46 47 f3 ~ 48 f1 49 f4 ~ 50 f2 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.637 0.039 16.261 0.000 0.541 0.541 -0.053 0.027 -1.971 0.049 -0.052 -0.052 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 0.264 0.028 9.314 0.000 0.264 0.264 -0.220 0.028 -7.764 0.000 -0.220 -0.220 Estimate Std.Err z-value P(>|z|) Std.lv Std.all 51 52 Covariances: 53 54 55 f1 ~~ f2 56 .f3 ~~ 57 .f4 58 59 Variances: 60 61 .Q1_1 1.623 0.051 31.654 0.000 1.623 0.822 62 .Q1_2 0.641 0.037 17.164 0.000 0.641 0.466 63 .Q1_3 0.872 0.045 19.547 0.000 0.872 0.512 64 .Q1_6 0.942 0.045 20.888 0.000 0.942 0.621 65 .Q1_7 0.874 0.059 14.784 0.000 0.874 0.505 66 .Q1_8 1.212 0.043 28.146 0.000 1.212 0.731 67 .Q1_11 1.619 0.066 24.478 0.000 1.619 0.610 68 .Q1_12 1.190 0.066 18.083 0.000 1.190 0.459 103

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について

69

.Q1_13

1.252

0.045

27.828

0.000

1.252

0.689

70

.Q1_16

0.657

0.044

14.817

0.000

0.657

0.265

71

.Q1_17

0.705

0.041

17.221

0.000

0.705

0.301

72

.Q1_18

1.431

0.048

30.098

0.000

73
74
75

1.431

0.564

f1

1.000

1.000

1.000

f2

1.000

1.000

1.000

[ reached getOption("max.print") -- omitted 2 rows ]

無事，小さい因子負荷を残して推定が完了しました。ESEM の利点は，このように小さい因
子負荷を許容することでモデル適合度が SEM よりも良くなる可能性がある，という点にあ
ります。試しに図7.30から青い点線を取り除いた SEM モデルと適合度を比べてみましょう。
1

# SEMモデルを作る

2

model_sem <- "

3

# 尺度A

4

f1 =~ Q1_1 + Q1_2 + Q1_3

5

f2 =~ Q1_6 + Q1_7 + Q1_8

6

# 尺度B

7

f3 =~ Q1_11 + Q1_12 + Q1_13

8

f4 =~ Q1_16 + Q1_17 + Q1_18

9
10

# 因子間の回帰

11

f3 ~ f1

12

f4 ~ f2

13

"

14
15

# 推定

16

result_sem <- sem(model_sem, data = dat, std.lv = TRUE)

17
18

# 適合度の比較

19

# SEMはESEMの特殊なケース（一部の因子負荷を0に固定）なのでnested

20

summary(compareFit(result_esem, result_sem))

1

################### Nested Model Comparison #########################

2
3

Chi-Squared Difference Test

4
5

Df

AIC

BIC

Chisq Chisq diff

6

result_esem 42 96399 96607 422.46

7

result_sem

50 96485 96647 524.93

RMSEA Df diff Pr(>Chisq)

102.47 0.069681

8

8
9

result_esem

10

result_sem

11

---

12

Signif. codes:

104

***
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

< 2.2e-16

7.14 （おまけ）lavaan で探索的因子分析ができるようになっていた件について 13 14 ####################### Model Fit Indices ########################### 15 chisq df pvalue 16 result_esem 422.462† 42 17 result_sem 18 524.928 50 rmsea cfi tli srmr aic .000 .061† .944† .913† .046† 96398.512† .000 .062 .931 .908 .054 96484.978 bic 19 result_esem 96607.185† 20 result_sem 96647.279 21 22 ################## Differences in Fit Indices ####################### 23 24 df rmsea result_sem - result_esem cfi tli srmr aic bic 8 0.001 -0.014 -0.004 0.007 86.466 40.094 今回の比較では，ESEM の結果の方が適合度が良いことがわかります。特に青い点線にあた るところに 0 ではない因子負荷があるような場合には，SEM で無理に 0 に制約をかけてし まうと適合度が悪化してしまう可能性がある，ということですね。 105

参考文献 参考文献 Adachi, K. (2020). Matrix-based introduction to multivariate data analysis (2nd ed.). Springer. https://doi.org/10.1007/978-981-15-4103-2 Akaike, H. (1998). Information theory and an extension of the maximum likelihood principle. In Parzen, E., Tanabe, K., and Kitagawa, G. (Eds.) Selected papers of Hirotugu Akaike (pp. 199–213.). Springer. https://doi.org/10.1007/978-1-4612-16940_15 (Original work published 1969) Albers, S. (2010). PLS and success factor studies in marketing. In Esposito Vinzi, V., Chin, W. W., Henseler, J., and Wang, H. (Eds.) Handbook of partial least squares: Concepts, methods and applications (pp. 409–425.). Springer Berlin Heidelberg. https: //doi.org/10.1007/978-3-540-32827-8 Asparouhov, T. and Muthén, B. (2009). Exploratory structural equation modeling. Structural Equation Modeling: A Multidisciplinary Journal, 16(3), 397–438. https: //doi.org/10.1080/10705510903008204 Baron, R. M. and Kenny, D. A. (1986). The moderator‒mediator variable distinction in social psychological research: Conceptual, strategic, and statistical considerations.. Journal of Personality and Social Psychology, 51(6), 1173–1182. https://doi.org/10 .1037/0022-3514.51.6.1173 Bayonne, E., Marin-Garcia, J. A., and Alfalla-Luque, R. (2020). Partial least squares (PLS) in operations management research: Insights from a systematic literature review. Journal of Industrial Engineering and Management, 13(3), 565–597. https: //doi.org/10.3926/jiem.3416 Bentler, P. M. (1983). Some contributions to efficient statistics in structural models: Specification and estimation of moment structures. Psychometrika, 48(4), 493–517. https://doi.org/10.1007/BF02293875 Bentler, P. M. (1990). Comparative fit indexes in structural models.. Psychological Bulletin, 107 (2), 238–246. https://doi.org/10.1037/0033-2909.107.2.238 Browne, M. W. (2001). An overview of analytic rotation in exploratory factor analysis. Multivariate Behavioral Research, 36(1), 111–150. https://doi.org/10.1207/S15327 906MBR3601_05 106

• https://doi.org/10.1007/978-981-15-4103-2
• https://doi.org/10.1007/978-1-4612-1694-0_15
• https://doi.org/10.1007/978-3-540-32827-8
• https://doi.org/10.1080/10705510903008204
• https://doi.org/10.1037/0022-3514.51.6.1173
• https://doi.org/10.3926/jiem.3416
• https://doi.org/10.1007/BF02293875
• https://doi.org/10.1037/0033-2909.107.2.238
• https://doi.org/10.1207/S15327906MBR3601_05

参考文献 Haenlein, M. and Kaplan, A. M. (2004). A beginner’s guide to partial least squares analysis. Understanding Statistics, 3(4), 283–297. https://doi.org/10.1207/s1532803 1us0304_4 Hair, J. F. (Ed.) (2017). A primer on partial least squares structural equation modeling (PLS-SEM) (2nd ed.). Sage. Hair, J. F., Hult, T. M., Ringle, C. M., Sarstedt, M., Danks, N. P., and Ray, S. (2021). Partial least squares structural equation modeling (PLS-SEM) using R: A workbook. Springer. Hancock, G. R., Stapleton, L. M., and Mueller, R. O. (2019). Structural equation modeling. In The reviewerʼs guide to quantitative methods in the social sciences (2nd ed., pp. 445–456.). Routledge. Henseler, J., Ringle, C. M., and Sarstedt, M. (2012). Using partial least squares path modeling in advertising research: Basic concepts and recent issues. In Okazaki, S. (Ed.) Handbook of research on international advertising (pp. 252–276.). E. Elgar. 星野崇宏・岡田謙介・前田忠彦 (2005). 構造方程式モデリングにおける適合度指標とモデ ル改善について: 展望とシミュレーション研究による新たな知見 行動計量学，32(2), 209–235. https://doi.org/10.2333/jbhmk.32.209 Hu, L.-t. and Bentler, P. M. (1999). Cutoff criteria for fit indexes in covariance structure analysis: Conventional criteria versus new alternatives. Structural Equation Modeling: A Multidisciplinary Journal, 6(1), 1–55. https://doi.org/10.1080/10705519909540118 Hwang, H., Sarstedt, M., Cheah, J. H., and Ringle, C. M. (2020). A concept analysis of methodological research on composite-based structural equation modeling: Bridging PLSPM and GSCA. Behaviormetrika, 47 (1), 219–241. https://doi.org/10.1007/s412 37-019-00085-5 Kang, H. and Ahn, J.-W. (2021). Model setting and interpretation of results in research using structural equation modeling: A checklist with guiding questions for reporting. Asian Nursing Research, 15(3), 157–162. https://doi.org/10.1016/j.anr.2021.06.001 Lin, H.-M., Lee, M.-H., Liang, J.-C., Chang, H.-Y., Huang, P., and Tsai, C.-C. (2020). A review of using partial least square structural equation modeling in e-learning research. British Journal of Educational Technology, 51(4), 1354–1372. https://doi.org/10.111 1/bjet.12890 Magno, F., Cassia, F., and Ringle, C. M. (2022). A brief review of partial least squares structural equation modeling (PLS-SEM) use in quality management studies. The TQM Journal, ahead-of-print, https://doi.org/10.1108/TQM-06-2022-0197 Mayekawa, S.-i. (1994). Equivalent path models in linear structural equation models. Behaviormetrika, 21(1), 79–96. https://doi.org/10.2333/bhmk.21.79 107

• https://doi.org/10.1207/s15328031us0304_4
• https://doi.org/10.2333/jbhmk.32.209
• https://doi.org/10.1080/10705519909540118
• https://doi.org/10.1007/s41237-019-00085-5
• https://doi.org/10.1016/j.anr.2021.06.001
• https://doi.org/10.1111/bjet.12890
• https://doi.org/10.1108/TQM-06-2022-0197
• https://doi.org/10.2333/bhmk.21.79

参考文献 McNeish, D. and Wolf, M. G. (2020). Thinking twice about sum scores. Behavior Research Methods, 52(6), 2287–2305. https://doi.org/10.3758/s13428-020-01398-0 Mulaik, S. A., James, L. R., Van Alstine, J., Bennett, N., Lind, S., and Stilwell, C. D. (1989). Evaluation of goodness-of-fit indices for structural equation models.. Psychological Bulletin, 105(3), 430–445. https://doi.org/10.1037/0033-2909.105.3.430 Rönkkö, M., McIntosh, C. N., and Antonakis, J. (2015). On the adoption of partial least squares in psychological research: Caveat emptor. Personality and Individual Differences, 87, 76–84. https://doi.org/10.1016/j.paid.2015.07.019 Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464. https://doi.org/10.1214/aos/1176344136 Steiger, J. H. and Lind, J. C. (1980). Statistically based tests for the number of common factors. Paper Presented at the Psychometric Society Annual Meeting 豊田秀樹 (2014). 共分散構造分析 [R 編] 東京図書 Tucker, L. R. and Lewis, C. (1973). A reliability coefficient for maximum likelihood factor analysis. Psychometrika, 38(1), 1–10. https://doi.org/10.1007/BF02291170 108

• https://doi.org/10.3758/s13428-020-01398-0
• https://doi.org/10.1037/0033-2909.105.3.430
• https://doi.org/10.1016/j.paid.2015.07.019
• https://doi.org/10.1214/aos/1176344136
• https://doi.org/10.1007/BF02291170