ねこすたっと

先ほど作図したプロットqに、当てはめられた直線*1を重ねてみます。

q + geom_abline(intercept = fit$coefficients[1],
                slope = fit$coefficients[2],
                colour = "indianred1",
                lwd = 1)

見るからに当てはまってないですね。

segmented modelとは

ここでは下の数式で表されるモデルを当てはめます。

$\begin{aligned} Y_i &\sim Binomial(1, \theta_i) \\ logit(\theta_i) &= \beta_0 + \beta_1*age_i + \beta_2*(age_i - \psi)*I(age_i > \psi) \end{aligned}$

ここで $\psi$ は変化点（breakpoint）、 $I(age_i > \psi)$ は年齢が変化点よりも大きいときは1を、大きくないときは0を返す関数です。係数 $\beta_2$ は、傾きの変化量（difference-in-slope）を表していて、年齢が変化点の手前のときは傾きが $\beta_1$ 、変化点よりも大きいときの傾きは $\beta_1 + \beta_2$ になるように推定されます。

Muggeo VM (2008) に推定方法について色々書かれていました。私が押さえ（られ）たのは次の2点です。

変化点においては2つの直線間の距離（gap）が0になるはず。このgapが小さくなるように変化点を探していく。
変化点があるということは、係数 $\beta_2$ が0ではないということ。帰無仮説 $H_0: \beta_2 =0$ が棄却されなかったら、そこは統計学的に有意な変化点がないと判断するということ。

変化点は必ずしも1つとは限りませんし、複数の変数について探索することも可能です。

segmented( )で変化点を探す

segmentedパッケージのsegmented( )で変化点を探す方法は意外に簡単でした。さっき当てはめたモデル（のオブジェクト）をsegmented( )に渡すだけです。

指定する代表的な引数は以下のとおり。

seg.Z：変化点を探したい変数名を~xのようにして指定する。複数指定したいときは~x+yのように追加する。
psi：変化点を探索する際の初期値

説明変数が1つしかなければseg.Zは明示しなくてもいいですし、初期値psiも必須ではありません。

fit.seg <- segmented(fit, seg.Z = ~age, psi = 35)

結果はsummary( )で取り出します。

> summary(fit.seg)

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.glm(obj = fit)

Estimated Break-Point(s):
            Est. St.Err
psi1.age 31.081  0.724

Meaningful coefficients of the linear terms:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -6.78244    0.43141 -15.722   <2e-16 ***
age         -0.01341    0.01795  -0.747    0.455    
U1.age       0.27422    0.02324  11.800       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null     deviance: 625.210  on 29  degrees of freedom
Residual deviance:  43.939  on 26  degrees of freedom
AIC: 190.82

Boot restarting based on 6 samples. Last fit:
Convergence attained in 5 iterations (rel. change 4.1199e-11)

変化点

Estimated Break-Point(s)の後には、検出された変化点とその推定誤差が書かれています。このデータだと31歳で傾きが変わっているようです。

confint.segmented( )で変化点の信頼区間を表示することができます。

> confint.segmented(fit.seg)
            Est. CI(95%).low CI(95%).up
psi1.age 31.0811     29.5925    32.5698

傾き

Meaningful coefficients of the linear termsのところには、当てはめられた直線の切片・傾きが書かれています。ageは変化点の前、U1.ageは変化点の後で傾きがどれくらい変化したか、を表しています。つまり先程の式で、 $\beta_1 = -0.01341$ , $\beta_2 = 0.27422$ ということです。

変化点でつながる2つの直線の切片・傾きは、intercept( ), slope( )で得ることができます。手書きで作図したいときはこれらを使います。

> intercept(fit.seg)
$age
               Est.
intercept1  -6.7824
intercept2 -15.3060

> slope(fit.seg)
$age
           Est.  St.Err.  t value CI(95%).l CI(95%).u
slope1 -0.01341 0.017947 -0.74722 -0.048586  0.021765
slope2  0.26081 0.014764 17.66500  0.231870  0.289750

slope2は $\beta_1 + \beta_2 = -0.01341 + 0.27422 = 0.26081$ となっていることを確認しました。

Davies検定

先程の結果サマリーで、U1.ageのP値は出ていませんでした。 U1.ageの係数 $\beta_2$ が有意に0と異なるか、つまり変化点があるかどうかを検定する方法がDavies検定（Davies' test）です。Muggeo VM (2008) によると、K個の変化点を設定したときに計算される統計量をもとに検定するようです。davies.test( )のエラーメッセージによるとKは10以上が推奨と出ますが、Muggeo VM (2008)ではKは5-10程度で十分と書かれています。

davies.test(fit, seg.Z=~age, k=10)

出力結果は以下のとおりです。

 Davies' test for a change in the slope

data:  formula = cases/births ~ age ,   method = glm 
model = binomial , link = logit , statist = lrt 
segmented variable = age
'best' at = 30.5, n.points = 10, p-value < 2.2e-16
alternative hypothesis: two.sided

p-value < 2.2e-16なので統計学的に有意な変化点があって、最良な点は30.5歳という結果でした。Kをめちゃくちゃ大きくすると、前述の推定値31.08に近くなっていくようです。

作図

推定したsegmented modelをグラフにしてみます。fit.segをplot( )に渡したら何となく描いてくれますが、以下の引数を知っておくと良さそうです。

conf.level：信頼区間の幅
shade：信頼区間の帯をつけるかどうか
link：リンク関数のスケールで描くかどうか。デフォルトはTRUEなので、前述の例で言うとy軸がロジット変換されたものになっている。

その他に線幅・線色など指定可能。

変化点に関しては、points( )で点推定値を、lines( )で信頼区間を追加することができるみたいです。

plot(fit.seg, 
     conf.level=0.95, 
     shade=TRUE,
     link=FALSE,
     lwd=2,col="red")
points(fit.seg,
       link=FALSE,
       pch=5,col="darkgreen")
lines(fit.seg,
      lwd=2,col="blue")

指定できる引数は、?plot.segmented, ?point.segmented, ?lines.segmentedで見ることができます。points( )ではlinkを指定しますが、lines( )では指定できません。

前述のコードの出力は次のとおり。

おわりに

混合効果モデルはnlmeパッケージのlme( )には対応しているようですが、lme4パッケージには対応していないようです（多分）。glmer( )で当てはめた一般線形モデルでsegmented regressionができるといいんですが...

参考資料

Muggeo VM (2008). “segmented: an R Package to Fit Regression Models with Broken-Line Relationships.” R News, 8(1), 20–25. https://cran.r-project.org/doc/Rnews/.

*1:y軸がロジット変換されているため

2023-03-24

「割合」や「率」に対してオフセット項付きポアソン回帰モデルを当てはめる [R]

dichotomous R neco

アウトカムが「割合」、というか「試行数と成功数」として与えられているときに二項回帰モデルを当てはめる方法は以前まとめました。

今回は「オフセット項」を使った回帰モデルに当てはめる方法をまとめてみようと思います。

オフセット項（offset）とは
使用するデータ
glm( )でオフセット項付きポアソン回帰モデルを当てはめる
おわりに
参考資料

オフセット項（offset）とは

オフセット項とは回帰係数か1に固定される説明変数のことです。 ポアソン回帰モデルはアウトカムが発生数などのカウント変数のときに使われますが、観察された総人数や期間が違えば発生数も違って当然ですよね。オフセット項を使えば、この「分母」をモデルに組み入れることができます。

通常の対数-ポアソン回帰モデルは、次の式で表されます。

$\begin{aligned} Y_i &\sim Poisson(\lambda_i) \\ \log(\lambda_i) &= X_i \beta \end{aligned}$

ポアソン分布は強度（intensity）λというパラメータで決まり、分布の平均・分散はともにλです。対数-ポアソンモデルでは、λの対数が説明変数の線形和（=Xβ）で決まる、として係数βが推定されます。つまり、発生数の期待値が説明変数で説明されています。

これに対して、オフセット項が入った対数-ポアソン回帰モデルは、次のようになります。

$\begin{aligned} Y_i &\sim Poisson(\lambda_i) \\ \log(\lambda_i) &= X_i \beta + \log(offset_i) \end{aligned}$

ここで、オフセット変数（offset）は、発生数の分母にくる観察された総人数や期間で、対数変換して加えます。第2式を変形すると、

$\begin{aligned} \log \left( \frac{\lambda_i}{offset_i} \right) &= X_i \beta \\ \end{aligned}$

となり、このモデルでは「発生数/オフセット」の期待値が説明変数で説明されていることになります。

オフセット変数が観察総数のときは「発生数/オフセット = 発生割合」になりますし、オフセット変数が期間のときは「発生数/オフセット = 単位期間あたりの発生数」、つまり発生率になります。

使用するデータ

今回はsegmentedパッケージのdownデータを使います。

library(segmented)
data(down)

このデータは母体年齢と出生児のダウン症（21トリソミー）の発生の関連を示したデータで、以下の3つの変数を含みます。

age：母体の平均年齢
births：出生数
cases：ダウン症を有する児の数

実際のデータはこんな感じ。

> head(down)
   age births cases
1 17.0  13555    16
2 18.5  13675    15
3 19.5  18752    16
4 20.5  22005    22
5 21.5  23896    16
6 22.5  24667    12

このデータのように、「試行数」と「成功数」が与えられている場合に二項回帰モデルを当てはめる方法は以前まとめました。

glm( )でオフセット項付きポアソン回帰モデルを当てはめる

分母に来てほしいbirthsをオフセット項に指定します（logを忘れないように）。書き方は2通り。

(1) offset( )を使ってモデル式の中に書く

fit.pois <- glm(cases ~ age + offset(log(births)),
                family=poisson(link=log), data=down)

(2) 引数offsetとして指定する

fit.pois <- glm(cases ~ age, offset=log(births),
                family=poisson(link=log), data=down)

結果はいずれも以下のとおり。

> summary(fit.pois)

Call:
glm(formula = cases ~ age + offset(log(births)), family = poisson(link = log), 
    data = down)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4070  -1.9389   0.5429   2.1276   4.7578  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -10.55277    0.21385  -49.35   <2e-16 ***
age           0.13713    0.00645   21.26   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 622.56  on 29  degrees of freedom
Residual deviance: 182.19  on 28  degrees of freedom
AIC: 325.26

Number of Fisher Scoring iterations: 5

おわりに

segmentedパッケージの勉強をしたついでに書いた記事なので、モデルの当てはまりは適切ではないかもしれません。
飼い猫の予防接種に初めて付き添ってきました。「ツメ切り希望」と書きたかったのですが、「瓜切り」と書いていたことに後で気づきました。

Statistical Methods In Medical Research, 4Ed [Paperback] [Jan 01, 2017] Armitage P.

参考資料

p.688にオフセットの記載があります。

作者:Books Wagon

2023-03-22

試行数と成功数が与えられたデータに対して二項回帰モデルを当てはめる [R]

dichotomous R

アウトカムが2値変数のときに使う回帰モデルについては、以前まとめたことがあります。

今回はアウトカムが「割合」、というか「試行数と成功数」として与えられているときの回帰モデルについてまとめてみます。

使用するデータ
cbind(成功数, 試行数)で二項モデルを当てはめる
割合を応答変数として試行数で重み付けして二項モデルを当てはめる
おわりに
参考資料

使用するデータ

今回はsegmentedパッケージのdownデータを使います。

library(segmented)
data(down)

このデータは母体年齢と出生児のダウン症（21トリソミー）の発生の関連を示したデータで、以下の3つの変数を含みます。

age：母体の平均年齢
births：出生数
cases：ダウン症を有する児の数

実際のデータはこんな感じ。

> head(down)
   age births cases
1 17.0  13555    16
2 18.5  13675    15
3 19.5  18752    16
4 20.5  22005    22
5 21.5  23896    16
6 22.5  24667    12

cbind(成功数, 試行数)で二項モデルを当てはめる

発生割合が（アウトカムの確率分布の）パラメータになっているモデルといえば、二項-ロジットモデル（ロジスティック回帰モデル）ですね。glm( )を使って、

fit <- glm(cases ~ age, family = binomial(link=logit), data = down)

と書きたいところですが、casesは2値変数ではないのでエラーになります。

こういうときは、cbind( )を使って次のように書きます。

fit.bin <- glm(cbind(cases, births) ~ age, 
               family = binomial(link=logit), data = down)

> summary(fit.bin)

Call:
glm(formula = cbind(cases, births) ~ age, family = binomial(link = logit), 
    data = down)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4053  -1.9359   0.5364   2.1314   4.7481  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -10.545157   0.213919   -49.3   <2e-16 ***
age           0.136872   0.006456    21.2   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 619.96  on 29  degrees of freedom
Residual deviance: 181.27  on 28  degrees of freedom
AIC: 324.15

Number of Fisher Scoring iterations: 5

リンク関数がロジットなので、係数を指数変換すればオッズ比として解釈できます。

割合を応答変数として試行数で重み付けして二項モデルを当てはめる

発生数を出生数で割って、発生割合にしてアウトカム変数とする方法です。

fit.prop <- glm(cases/births ~ age, 
               family = binomial(link=logit), data = down)

とすると「成功数が整数になっていない！」と怒られます。

次のように、引数weightに試行数（ここでは出生数）を指定します。試行数が大きいところの方が信頼度が高いので、全体の結果に寄与する程度も大きいからです。

fit.prop <- glm(cases/births ~ age, weight = births, 
               family = binomial(link=logit), data = down)

> summary(fit.prop)

Call:
glm(formula = cases/births ~ age, family = binomial(link = logit), 
    data = down, weights = births)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-3.412  -1.943   0.546   2.140   4.767  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -10.562280   0.214342  -49.28   <2e-16 ***
age           0.137524   0.006469   21.26   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 625.21  on 29  degrees of freedom
Residual deviance: 183.86  on 28  degrees of freedom
AIC: 326.74

Number of Fisher Scoring iterations: 5

だいたい似たような推定値になりました。

おわりに

(Dispersion parameter for binomial family taken to be 1)というメッセージが出てきているので、過剰分散（overdispersion）の話とか、擬似二項分布（quasibinomial）の話とかをまとめようと思っていますが、それはまた別の機会に。
「いやいや、割り算なんかせずにオフセット項のあるポアソン回帰だろ！」というご意見もあると思いますので、これも別の機会に。

参考資料

glm( )のヘルプに書いてあります。

2023-02-12

Hoyle拾い読み：検証的因子分析（CFA）

SEM latent variable neco

構造方程式モデリング（Structural Equation Modeling, SEM）の初学者が、タイトルの "Handbook" に誘われて買ってしまったHoyle先生の分厚い本を拾い読みしたメモです。といっても、私にとって重すぎる内容は拾い上げられていません。

今回は "Chapter 22：Confirmatory Factor Analysis" から。これまでもCFAについてまとめたことがありますが、実装方法が中心だったので、今回は教科書で学び直しです。

因子分析（factor analysis, FA）とは
CFAモデルに含まれるパラメータ
パラメータ推定値の測定単位系
モデルから計算される分散・共分散の推定値
推定結果で確認してみる
モデルの評価
おわりに

因子分析（factor analysis, FA）とは

複数の観測変数（indicator）から潜在的な共通性を因子（factor）として抽出することを目的としています。観測変数の分散を、共通因子によって説明される部分（共通分散, common variance）と、個々の観測変数に固有の部分（独自分散, unique variance）に分けることで、因子-観測変数関係（factor-indicator relationship）を見ようとする解析です。

FAには次の2種類があります。

探索的因子分析（exploratory FA, EFA）：

因子の数や因子-観測変数間の関係性を事前に設定しない
探索的・記述的
測定誤差同士は独立しているという仮定が必要
最大の因子負荷を1、最小の因子負荷を0に近づけるように回転させる

検証的因子分析（confirmatory FA, CFA）：

研究者が先行研究やその領域の知見・経験をもとにして、因子の数や因子-観測変数間の関係性を事前に設定する
事前の経験的・理論的背景を確立するために行われる
測定誤差同士に相関関係があっても良い
因子負荷は回転させない

CFAは理論的に構築されたモデルの妥当性検証に必要不可欠な解析ツールです。理論的に類似していると考えられる観測変数の相関が強いことを確かめたり（convergent validity）、逆に理論的に異なっていると考えられる観測変数間の相関がそれほど高くないことを確かめたり（discriminant validity）するのに使われます。

他にも、ある測定モデルが地理的・文化的に異なった集団においても同じように用いることができるかを評価する際にも用いられます。

SEMは事前にモデルを設定するのでEFAとは別物です。ここではCFAのみ扱います。

CFAモデルに含まれるパラメータ

下のパス図は2因子6項目のCFAの例です。

CFAに含まれるパラメータは、

因子負荷量（factor loading）：因子が観測変数に与える影響の度合い（ $\lambda_{i}$ ）
独自分散（unique variance）：各観測変数に固有に与えられた独自因子の分散（ $\varepsilon_{i,i}$ ）
共通分散（common variance）：因子の分散（ $\phi_{i,i}$ ）

があります。これらに加えて、

因子間共分散 $\phi_{i,j}$ （covariance between factors）
誤差間共分散 $\varepsilon_{i,j}$ （error covariance）*1

を設定することができます。

パラメータ推定値の測定単位系

潜在変数には測定単位がないので、何らかの制約（単位基準）を設ける必要があります。

方法1：変数を標準化する方法
変数の分散を1に固定する方法です。これを標準化（standardized）*2と言います。

潜在変数・観測変数の両方を標準化する場合をcompletely standardization、片方だけを標準化する場合をpartially standardizationと呼びますが、とりあえず前者を知っていれば事足りそうです。

方法2：マーカー変数法
それぞれの因子に属する観測変数のうち、1つの変数の因子負荷量を1に固定する方法です。方法1に対して、こちらでは変数は必ずしも標準化されていません（unstandardized）。

モデルから計算される分散・共分散の推定値

モデルが正しいという仮定のもとでは、目的の観測変数間の経路で因子負荷量と因子の分散を掛け合わせていくことで求められます。例えば、

同じ因子内の観測変数の共分散は、

$\begin{aligned} Cov(X4, X5) = \lambda_4 \phi_{2,2} \lambda_5 \end{aligned}$

異なった因子内の観測変数の共分散は、

$\begin{aligned} Cov(X3, X4) = \lambda_3 \phi_{1,2} \lambda_4 \end{aligned}$

となります。

観測変数の誤差間に相関がある場合は、誤差間の共分散を足します。

$\begin{aligned} Cov(X1, X2) = \lambda_1 \phi_{1,1} \lambda_2 + \varepsilon_{1,2} \end{aligned}$

この式で2つの変数を同じにすると、観測変数の分散が属する因子の共通分散と独自分散の和になることが導かれます。

$\begin{aligned} Var(X1) &= Cov(X1, X1) \\ &= \lambda_1 \phi_{1,1} \lambda_1 + \varepsilon_{1,1} \\ &= \lambda_1 \phi_{1,1} \lambda_1 + \varepsilon_{1} \end{aligned}$

一方、観測されたデータからは、

個々の観測変数の分散
観測変数間の共分散

が得られます。モデルが正しいという仮定のもとでパラメータを使って表した分散・共分散と、実際に観測された分散・共分散を比べることでパラメータを推定していきます。

推定結果で確認してみる

以前lavaanパッケージでCFAを実行する方法をまとめたときと、同じ例を使って推定結果を確認してみます。

Handbook of Structural Equation Modeling

# パッケージとデータの準備
library(lavaan)
data <- HolzingerSwineford1939[-c(1:6)]

まずは使用するデータの共分散行列を見てみます。

> round(cov(data),digits=3)
      x1     x2    x3    x4    x5    x6     x7    x8    x9
x1 1.363  0.409 0.582 0.507 0.442 0.456  0.085 0.265 0.460
x2 0.409  1.386 0.453 0.210 0.212 0.248 -0.097 0.110 0.245
x3 0.582  0.453 1.279 0.209 0.113 0.245  0.089 0.213 0.375
x4 0.507  0.210 0.209 1.355 1.101 0.899  0.220 0.126 0.244
x5 0.442  0.212 0.113 1.101 1.665 1.018  0.143 0.181 0.296
x6 0.456  0.248 0.245 0.899 1.018 1.200  0.145 0.166 0.237
x7 0.085 -0.097 0.089 0.220 0.143 0.145  1.187 0.537 0.375
x8 0.265  0.110 0.213 0.126 0.181 0.166  0.537 1.025 0.459
x9 0.460  0.245 0.375 0.244 0.296 0.237  0.375 0.459 1.018

相関係数rと共分散Covの関係は、

$\begin{aligned} r_{x,y} = \frac{Cov(x,y)}{SD(x)SD(y)} \end{aligned}$

なので、共分散行列で変数を標準化したものは相関行列になります。

> round(cor(data),digits=3)
      x1     x2    x3    x4    x5    x6     x7    x8    x9
x1 1.000  0.297 0.441 0.373 0.293 0.357  0.067 0.224 0.390
x2 0.297  1.000 0.340 0.153 0.139 0.193 -0.076 0.092 0.206
x3 0.441  0.340 1.000 0.159 0.077 0.198  0.072 0.186 0.329
x4 0.373  0.153 0.159 1.000 0.733 0.704  0.174 0.107 0.208
x5 0.293  0.139 0.077 0.733 1.000 0.720  0.102 0.139 0.227
x6 0.357  0.193 0.198 0.704 0.720 1.000  0.121 0.150 0.214
x7 0.067 -0.076 0.072 0.174 0.102 0.121  1.000 0.487 0.341
x8 0.224  0.092 0.186 0.107 0.139 0.150  0.487 1.000 0.449
x9 0.390  0.206 0.329 0.208 0.227 0.214  0.341 0.449 1.000

次にCFAモデルを設定してパラメータを推定します。

# 使用するCFAモデル
model = "
F1 =~ x1 + x2 + x3
F2 =~ x4 + x5 + x6
F3 =~ x7 + x8 + x9
"

# 標準化した場合
fit_std <- cfa(model = model, data=data, estimator = "MLR", std.lv=TRUE)
# 標準化しない場合
fit_unstd <- cfa(model = model, data=data, estimator = "MLR", std.lv=FALSE, auto.fix.first = TRUE)

結果1：潜在変数を標準化した場合

確認に使うところ以外は省略してます。

summary(fit_std, fit.measures=TRUE, standardized=TRUE)

~~~~~（省略）~~~~~
Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  F1 =~                                                                 
    x1                0.900    0.100    8.973    0.000    0.900    0.772
    x2                0.498    0.088    5.681    0.000    0.498    0.424
    x3                0.656    0.080    8.151    0.000    0.656    0.581
  F2 =~                                                                 
    x4                0.990    0.061   16.150    0.000    0.990    0.852
~~~~~（省略）~~~~~

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  F1 ~~                                                                 
    F2                0.459    0.073    6.258    0.000    0.459    0.459
~~~~~（省略）~~~~~

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .x1                0.549    0.156    3.509    0.000    0.549    0.404
   .x2                1.134    0.112   10.135    0.000    1.134    0.821
   .x3                0.844    0.100    8.419    0.000    0.844    0.662
   .x4                0.371    0.050    7.382    0.000    0.371    0.275
~~~~~（省略）~~~~~
    F1                1.000                               1.000    1.000
    F2                1.000                               1.000    1.000
    F3                1.000                               1.000    1.000

(1) まずは観測変数の分散についてです。

F1-x1の因子負荷量 = 0.900（← Latent Variablesに記載）
F1の分散 = 1.000（← Variancesに記載）

から、x1の共通分散の推定値は0.900 × 1.000 × 0.900 = 0.810となります。

また、Variancesに記載されている.x1はx1の独自分散の推定値で0.549です。両者を足すと0.810 + 0.549 = 1.359となり、x1の分散1.363に近い数値になっているので、推定値としてはまずまず良さそうです（過剰識別のモデルで再現性が悪ければかけ離れた数値になる）。

(2) 次に同じ因子に属する観測変数（x1, x2）の共分散です。

F1-x1の因子負荷量 = 0.900
F1の分散 = 1.000
F1-x2の因子負荷量 = 0.498

から、Cov(x1,x2)の推定値 = 0.900 × 1.000 × 0.498 = 0.448となります。実際の値は0.409ですのでまずまず良い推定値に見えます。

(3) 最後に異なった因子に属する観測変数（x1, x4）の共分散です。

F1-x1の因子負荷量 = 0.900
F1-F2の共分散 = 0.459
F2-x4の因子負荷量 = 0.990

なので、Cov(x1,x4)の推定値 = 0.900 × 0.459 × 0.990 = 0.409となります（実際の値は0.507でした）。

結果2：マーカー変数の因子負荷量を1に固定した場合

次にマーカー変数法（方法2）の場合です。

summary(fit_unstd, fit.measures=TRUE, standardized=TRUE)

~~~~~（省略）~~~~~
Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  F1 =~                                                                 
    x1                1.000                               0.900    0.772
    x2                0.554    0.132    4.191    0.000    0.498    0.424
    x3                0.729    0.141    5.170    0.000    0.656    0.581
  F2 =~                                                                 
    x4                1.000                               0.990    0.852
~~~~~（省略）~~~~~

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  F1 ~~                                                                 
    F2                0.408    0.099    4.110    0.000    0.459    0.459
~~~~~（省略）~~~~~

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .x1                0.549    0.156    3.509    0.000    0.549    0.404
   .x2                1.134    0.112   10.135    0.000    1.134    0.821
   .x3                0.844    0.100    8.419    0.000    0.844    0.662
   .x4                0.371    0.050    7.382    0.000    0.371    0.275
~~~~~（省略）~~~~~
    F1                0.809    0.180    4.486    0.000    1.000    1.000
    F2                0.979    0.121    8.075    0.000    1.000    1.000
    F3                0.384    0.107    3.596    0.000    1.000    1.000

先程と同じ内容を確認してみます。

(1) まずは観測変数の分散についてです。

F1-x1の因子負荷量 = 1.000（← Latent Variablesに記載）
F1の分散 = 0.809（← Variancesに記載）

なので、x1の共通分散の推定値は1.000 × 0.809 × 1.000 = 0.809となります。標準化したときの値0.810とほぼ同じです（ズレは丸めの影響でしょう）。

(2) 次に同じ因子に属する観測変数（x1, x2）の共分散です。

F1-x1の因子負荷量 = 1.000
F1の分散 = 0.809
F1-x2の因子負荷量 = 0.554

から、Cov(x1,x2)の推定値 = 1.000 × 0.809 × 0.554 = 0.448となり、標準化したときと同じ値になりました。

(3) 最後に異なった因子に属する観測変数（x1, x4）の共分散です。

F1-x1の因子負荷量 = 1.000
F1-F2の共分散 = 0.408
F2-x4の因子負荷量 = 1.000

なので、Cov(x1,x4)の推定値 = 1.000 × 0.408 × 1.000 = 0.408で、標準化したときの値0.409と同じ値です。

寄与率

Std.lvは潜在変数（lv = latent variable）のみ標準化した結果、Std.allは観測変数も標準化した結果です。

ここで再々度、観測変数の分散についてです。

結果1でStd.allの値を使うと、F1-x1の因子負荷量 = 0.772なので、 x1の共通分散は0.772 × 1.000 × 0.772 = 0.596です。 x1の独自分散は0.404なので、両者の和は0.596 + 0.404 = 1となり、観測変数についても標準化されていることが確認できました。

ここで、分散全体に対する共通分散の割合（=0.596）は寄与率と呼ばれ、x1の変動のうち59.6%が因子F1によって説明されることを示しています。

モデルの評価

設定したモデルを使ってパラメータが推定できたら、以下の3点について評価します。項目を挙げるだけなので、詳細は成書で確認してください。

モデル全体の当てはまり
モデル局所の当てはまり
パラメータ推定値の解釈

設定したモデルがイマイチだった場合、因子数や変数同士の関係性を見直すことがあります（respecification）。

モデル全体の当てはまり

$\chi^2$ 値で評価すると、サンプル数が多いときにほとんど必ず有意差が出てしまう（つまり、モデルのデータがずれている）という問題があるそうです。なので、他にも次のような指標を使って評価します。

standardized root mean square residual（SRMR）：<0.08ならOK
root mean square error of approximation（RMSEA）：<0.06ならOK
Tucker-Lewis index（TLI）：>0.95ならOK
comparative fit index（CFI）：>0.95ならOK

モデル局所の当てはまり

全体の当てはまりが良くても、ある観測変数については再現性がよくないことがあります。次の指標を追加います。

standardized residuals：>1.96だとズレが大きい
modification index：パラメータに関する制約を外したときに $\chi^2$ 値がどれくらい変化すうるか。df=1の $\chi^2$ 値が3.84なので、これを超えて低下する場合はパラメータに関する制約が誤っているかもしれない。

パラメータ推定値の解釈

あくまでモデルの当てはまりが良いことが前提条件です。次の項目を確認します。

パラメータが取りうる範囲に収まっているか
大きさ・方向が概念上・経験上の観点から説明可能か
因子間の相関が大きすぎないか。1に近いときは、因子で重複している要素が多いことになる。

おわりに

CFAの章を読んで、少しSEMに対する苦手意識が和らぎました。最初にこの章を読めばよかったです。
なぜか数式のフォントが変です。どうしたら直るのでしょうか。（→追記：公開したら直りました）
猫が布団に潜り込んでこなくなって春の訪れを感じました。

Guilford Press

*1:他にもcorrelated error, correlated residual, correlated uniqueness

*2:平均=0, 分散=1に固定すること

2023-02-08

Hoyle拾い読み：SEMにおけるモデル設定と識別条件

SEM latent variable testing terminology

今回は "Chapter 8：Model Specification in SEM" と "Chapter 9：Identification: A Nontechnical Discussion of a Technical Issue" から。

モデル設定（specification）のアプローチ
モデル識別（identification）の条件
- 自由度
- 過小識別・適度識別・過剰識別
識別不能への対処方法
おわりに
参考資料

モデル設定（specification）のアプローチ

モデルを設定する際のアプローチとして次の4つが紹介されています。

strictly confirmatory approach：
研究者がその分野の背景から想定した1つのモデルについてのみ、データとの適合性を評価するアプローチ。あまり使われない。

alternative models approach：
手元のデータを見ずに設定した代替モデルと元のモデルと比較するアプローチ。

model-generating approach：
設定したモデルを解析し修正していくアプローチ。別のデータで当てはまりを検証することを念頭においている。

model discovery approach：
変数から想定される全てのモデルを自動的に評価して、相応しいモデルを探索するアプローチ。

2番目と3番目が多そうです。

モデル識別（identification）の条件

パラメータの値が1つに定まるとき、モデルは識別可能であると言います。先に言ってしまうと、モデルが式べえつ可能な条件を数式で理解するのは（私には）非常に難しそうです。実際にPCで計算してみて、識別できるかどうかで判断することになりそう...。

自由度

モデルが識別されるためには、「データから与えられる既知の情報」が「推定すべき未知の情報」よりも多いことが必要条件であることは直感的に理解できます（「既知の情報>未知の情報」だからといって必ずしも同定可能とは限らない）。 既知の情報個数と未知の情報個数の差を、自由度（degree of freedom, df）と呼びます。

k個の観測変数がある場合の既知情報は以下のとおりです。

k個の分散（variance）
k(k-1)/2個（←2組の異なる変数の組み合わせの数）の共分散（covariance）
必要ならk個の平均（mean）

これに対し、未知の情報はモデルの中で推定が必要なパラメータ（分散・共分散、係数、因子負荷量）の個数です。

過小識別・適度識別・過剰識別

既知の情報が未知の情報よりも少ない状態（df<0）を過小識別（underidentified）と言います。既知の情報が不足しているので、モデルを識別することができません。一見、既知の情報が足りてそうでも、非常に類似した変数が含まれていて実際には解が得られないことを、経験的過小識別*1（empirically underidentification）と言います。回帰モデルでいうところの多重共線性と似たような話だと思います。

既知の情報と未知の情報が等しい状態（df=0）を適度識別（just-identified）と言います。全てのパラメータに対して、ちょうど1つの推定値が得られます。1つしか得られないということはバラツキがないということなので、推定誤差は分かりませんし、モデルの当てはまり具合を評価することはできません。

既知の情報が未知の情報よりも多い状態（df>0）を過剰識別（overidentified）と言います。全ての条件をちょうど満たすことができる解が存在せず（解なし）、そのため得られた解は多少なり誤りを含んでいることになります。この誤差を評価することでモデルへの当てはまり指標が得られます。

識別不能への対処方法

モデルが識別できないときは、(1)モデルが複雑すぎる、(2)少なくとも1つのパラメータが固定されていない、が原因として考えられます。モデルの構造については、特に形成的因子（formative factor）で識別不能が問題になりやすいとのこと。1つの因子に対して3-4個のindicatorが望ましいと書かれています。

本の中ではモデルのパターンごとに識別に必要な条件が書かれています。例えば、フィードバックループのないパス解析モデルでは、「各々の内生変数のdisturbanceが、その内生変数の原因となっている変数と関連がないこと」が必要です（regression rule）。通常の回帰モデルで言えば、「従属変数の誤差項が説明変数と独立である」ということにあたります。

検証的因子分析（CFA）については、まずlatent variableの測定単位を固定する必要があります。このためには、

1つの因子負荷量を1に固定する
外生的潜在変数の分散を1に固定する

などが必要になります。CFAについてはこのほかにも、

同一の因子に属する観測変数同士の制約
異なった因子に属する観測変数同士の制約
同一の因子の観測変数間と他の因子の観測変数との間に関する制約

について説明されています（vanishing tetrads）。

フィードバックループのあるパス解析モデル、欠落変数があるパス解析モデル*2に関して操作変数の話も書かれていますが、詳細はまた別の機会に勉強します。

おわりに

非常に表面的になぞっただけで終わってしまいました。

How to Calculate Degrees of Freedom in a SEM Model - YouTube

参考資料

欠落変数バイアスについて参考にしました。

yukiyanai.github.io

SEMにおける自由度の計算方法についての解説動画。

Handbook of Structural Equation Modeling

Guilford Press

Handbook of Structural Equation Modeling

*1:訳語があっているか分かりません

*2:疫学でいうところの未測定交絡。計量経済学の領域と疫学の領域ではここらへんの言葉の使い方が微妙に違うので注意が必要そう

2023-02-07

Hoyle拾い読み：SEMのパス図と構成要素

SEM latent variable terminology neco

今回は "Chapter 3：Graphical Representation of SEM using Path Diagram" と "Chapter 4：Latent Variables in SEM" から。

ノードの種類
矢印の種類
外生変数と内生変数
モデルの分類
- 観測モデル
- 構造モデル
おわりに

SEMを表す方法として、

パス図を使って表す方法
個々の関係性を方程式を使って表す方法
方程式を行列を使ってまとめて表す方法

が紹介されていますが、直感的に1番理解しやすい「パス図を使った表記方法」をまとめてみます。

下の例のように、パス図は丸や四角で囲った変数（定数もあるのでまとめてノードと呼びます）と、それらを結ぶ矢印で描かれています。

ノードの種類

パス図に含まれるノードには以下の3種類あります。

観測変数（observed variable）
潜在変数（latent variable）
切片/定数項（intercept/constant term）

観測変数はデータとして直接測定される変数です。後述の因子（factor）と繋がっているものは、"Indicator" とも呼ばれます。パス図の中では四角囲みで表します。

潜在変数は概念的で、データとして測定されない変数です。潜在変数には、研究において興味がある理論的な要素である共通因子（common factor）あるいは単に因子（factor）と、測定や関係性の揺らぎに相当する誤差項（error term）があります。潜在変数はパス図の中では丸囲みで表します。

切片/定数項は、変数ではなく固定された数字で、大抵の場合は1に固定されています。パス図の中では三角囲みで表します（冒頭のパス図例2参照）。

矢印の種類

Nodeをつなぐ矢印には、

単方向矢印（→）
双方向矢印（↔︎）

があります。

単方向矢印（→）は、原因X→結果Yのように、XからYへの直接効果（direct effect）を表します。
ここで、変数Xは独立変数（independent variable）や予測変数（predictor）などと、変数Yは従属変数（dependent variable）, アウトカム（outcome）, 効果（effect）*1などと呼ばれます。

双方向矢印（↔︎）は、方向性のない2変数の関連（nondirectional association）を表すほか、外生変数（後述）の分散を表すときにも用いられます。

矢印の隣には次のようなパラメータが書き添えられることがあります。

因子負荷量
重み付け係数（回帰係数）
分散・共分散

データをもとに推定したいパラメータ（free parameter）はギリシャ文字や "＊" などが添えられます。ある値に固定されているパラメータ（fixed parameter）の場合はその数値を書きます（ほとんどの場合1）。それぞれの対象者の変数の値によって係数（因子負荷量など）が変わる場合は、その変数（definition variableと呼びます）を菱形で書き添えることもあります（冒頭のパス図例2参照）。

外生変数と内生変数

モデルの中で「親」となる原因を有するかどうかによる分類です。

モデルの中で原因となる変数がないものを外生変数（exogeneous variable）と呼びます。外生変数の平均・分散・共分散は、外から値が与えられることになります。

外生変数ではないものは全て内生変数（endogeneous variable）です。内生変数の平均・分散・共分散は、外生変数の関数として与えられます。

モデルの分類

モデルの分類といっても自由度が高いので、ここでは全体を分類するのではなく、モデルの中に見られる「パターン」を見ていきます。矢印が刺さる変数を左側にして数式を書くと、右辺はいずれの場合も線型結合（つまり、βXの和）になっているもののみ扱います（非線形の扱いはまだコンセンサスなし）。

観測モデル

観測モデル（measurement model）とは、潜在変数とそれに繋がる観測変数の関係性を示したモデルです。因子とindicatorの関係性によって、さらに次の2つに分類されます。

形成的モデル（formative model）：
観測変数Xによって因子Fが形成されるモデル、つまりパス図で "X→F" と示される部分です。

観測変数（formative indicator, cause indicator）の係数は回帰の重み付け係数（regression weight）
因子の測定誤差（= 撹乱項, disturbance）を設定しない場合は、因子ではなく単なる重み付け合計スコアになる
観測変数同士には相関があってもなくてもよい

反映的モデル（reflective model）：
観測変数Xを特徴づける因子Fが上流にあるモデル、つまりパス図で "F→X" と示される部分です。子にあたる観測変数の共通性（commonality）が因子Fに反映されることになるのでこのような名前になっています。

観測変数（reflective indicator, effect indicator）の係数は因子負荷量（factor loading）
「観測変数の測定誤差（= 測定誤差, measurement error）と因子の間には相関関係がない」と設定される（多分モデル識別に必要）
観測変数同士は因子を介して相関関係を持つ。因子の影響が強いほど、観測変数間の関連は強くなる

構造モデル

これに対し、構造モデル（structural model）は潜在変数同士の関係性を示したモデルです。因子のindicatorになっていない観測変数同士の関係性も構造モデルに分類されます。*2

構造モデルの頻出パターンはまた別の機会にまとめます。

閾値モデル（threshold model）
潜在成長モデル（latent growth model）
成長混合モデル（growth mixture model）
2次因子モデル（second-order model）, 高次因子モデル（high-order model）
双因子モデル（bifactor model）

おわりに

テキストを読むにはまず用語から。
iRobot（ルンバ j7+）が来て1週間、ネコたちも少し慣れてきました。

Guilford Press