1. 導入
対数線形モデルは、分割表の各セルの値がどのようなモデルで生起したかを説明するモデルです。分割表が絡んだ分析では、一般的に利用される手法のようですが、恥ずかしながら、最近知りました。今回もお勉強のまとめですが、いまいちピンと来ていないので、間違いがありましたら、指摘していただきたいです。
2. 対数線形モデル
観測度数について、以下のような2元分割表が与えられているとします。
| $B_1$ | $B_2$ | $\cdots$ | $B_c$ | 計 | |
|---|---|---|---|---|---|
| $A_1$ | $n_{11}$ | $n_{12}$ | $\cdots$ | $n_{1c}$ | $\sum_{j=1}^cn_{1j}$ |
| $A_2$ | $n_{21}$ | $n_{22}$ | $\cdots$ | $n_{2c}$ | $\sum_{j=1}^cn_{2j}$ |
| $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | |
| $A_r$ | $n_{r1}$ | $n_{r2}$ | $\cdots$ | $n_{rc}$ | $\sum_{j=1}^cn_{rj}$ |
| 計 | $\sum_{i=1}^rn_{i1}$ | $\sum_{i=1}^rn_{i2}$ | $\cdots$ | $\sum_{i=1}^rn_{ic}$ | n |
この分割表に対する、確率分布表は
| $B_1$ | $B_2$ | $\cdots$ | $B_c$ | 計 | |
|---|---|---|---|---|---|
| $A_1$ | $p_{11}$ | $p_{12}$ | $\cdots$ | $p_{1c}$ | $\sum_{j=1}^cp_{1j}$ |
| $A_2$ | $p_{21}$ | $p_{22}$ | $\cdots$ | $p_{2c}$ | $\sum_{j=1}^cp_{2j}$ |
| $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | |
| $A_r$ | $p_{r1}$ | $p_{r2}$ | $\cdots$ | $p_{rc}$ | $\sum_{j=1}^cp_{rj}$ |
| 計 | $\sum_{i=1}^rp_{i1}$ | $\sum_{i=1}^rp_{i2}$ | $\cdots$ | $\sum_{i=1}^rp_{ic}$ | 1 |
となります。独立性の仮定の下、各セルの生起確率は
$$p_{ij} = P(A_i\cap B_j) = P(A_i)P(B_j) = p_{i.}p_{.j}$$
であり、確率の推定量と期待度数は
\begin{align}
&\hat{p}_{ij} = \hat{p}_{i.} \hat{p}_{.j} = \frac{n_{i.}}{n}\frac{n_{.j}}{n}\tag{1}\\
&\hat{n}_{ij} = n\hat{p}_{ij} = \frac{n_{i.}n_{.j}}{n}\tag{2}
\end{align}
であたえられます。期待度数の対数をとると
\begin{align}
\log(\hat{n_{ij}}) &= -\log n + \log n_{i.} + \log n_{.j}\\
&:= \beta + \{A\}_{i} + \{B_{j}\}
\end{align}
と表されます。このことから、期待度数は、分割表全体に共通する効果$\beta$と、要因$A_{i},B_{j}$の主効果の和としてあらわされることがわかります。対数線形モデルでは、これに要因$A_{i},B_{j}$の交互作用を表す項${A*B}_{ij}$を加えた、
\begin{align}
\log(\hat{n_{ij}}) := \beta + \{A\}_{i} + \{B_{j}\} + \{A*B\}_{ij}
\end{align}
を用いて、期待度数をモデリングします。このモデルは飽和モデルと呼ばれます。交互作用を含めるかどうかは、先行研究に従ったり、一旦交互作用項を入れたモデルで係数を推定し、係数が0であるとする帰無仮説の下での検定結果に基づいて判断します。係数の推定のためには
\begin{align}
\sum_{i}\{A\}_{i} = \sum_{j}\{B\}_{j} = \sum_{i}\{A*B\}_{ij} = \sum_{j}\{A*B\}_{ij} = 0
\end{align}
という制約を課します。係数の推定については、このページをご覧ください。
3. Rでの実行例
このままだと抽象的すぎるので、もう少し具体的な例を考えてみましょう。ある町で、300名を抽出し、風邪を予防するための怪しい民間療法を行うことの効果を調査したとします(適当に発生させました)。調査の結果、以下の表が得られました。
| 風邪を引いた | 風邪を引かなかった | |
|---|---|---|
| 民間療法あり | 76 | 104 |
| 民間療法なし | 47 | 73 |
この表の各セルの観測度数は、「怪しい民間療法」と「風邪」からどのような影響をうけて生成されているのでしょうか。対数線形モデルに当てはめてみます。Rでは一般化線形モデルの枠組みで対数線形モデルの係数の推定を行うことができます。そのため、ここではglm()を用いた解析例を載せます。loglin()やloglm()に関してはこちらのqiita記事(対数線形モデルですべての係数を表示させる方法 1)がとても参考になります。(というか、このリンク先とさきほどのリンク先を見ながらまとめています。このリンク先のほうが詳しいです。)
# データの用意
data = data.frame(民間療法 = c(rep("あり",2),rep("なし",2)),
風邪 = c(rep(c("ひいた","ひいてない"),2)),
freq = c(76,104,47,73))
data #データの表示
res = glm(freq ~ 民間療法 + 風邪 + 民間療法:風邪, data = data, family = "poisson")
# 被説明変数の分布をpoisson分布とすると対数線形モデルとなる。
# 被説明変数が二値データ(=ロジスティック回帰)だと、family = "binomial"と指定します。
summary(res)
実行結果は以下のようになりました。
> summary(res)
Call:
glm(formula = freq ~ 民間療法 + 風邪 + 民間療法:風邪, family = "poisson",
data = data)
Deviance Residuals:
[1] 0 0 0 0
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.3307 0.1147 37.754 <2e-16 ***
民間療法なし -0.4806 0.1856 -2.590 0.0096 **
風邪ひいてない 0.3137 0.1509 2.078 0.0377 *
民間療法なし:風邪ひいてない 0.1267 0.2403 0.527 0.5982
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 2.2133e+01 on 3 degrees of freedom
Residual deviance: 2.0650e-14 on 0 degrees of freedom
AIC: 32.477
Number of Fisher Scoring iterations: 2
推定結果は、「民間療法なし」かつ「風邪を引いていない」セルの係数を表示しています。この推定結果だと、民間療法と風邪の交互作用項の係数は有意ではないということが分かります。このセルに対するモデル式をかくと、
\begin{align}
\log(n_{22}) &= \mbox{切片} + \mbox{民間療法なしの係数} + \mbox{風邪を引いていないの係数} + \mbox{なしなしの交互作用}\\
&= 4.3307 -0.4806 + 0.3137 + 0.1267 = 4.2905\\
e^{4.2905} &\approx 73
\end{align}
となります。しかし、ここで、交互作用項は有意でないので、これを外して再び推定してみると、
res2 = glm(freq ~ 民間療法 + 風邪, data = data, family = "poisson")
summary(res2)
> summary(res2)
Call:
glm(formula = freq ~ 民間療法 + 風邪, family = "poisson", data = data)
Deviance Residuals:
1 2 3 4
0.2548 -0.2142 -0.3160 0.2601
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.3014 0.1017 42.275 < 2e-16 ***
民間療法なし -0.4055 0.1179 -3.440 0.000581 ***
風邪ひいてない 0.3640 0.1174 3.101 0.001932 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 22.13286 on 3 degrees of freedom
Residual deviance: 0.27837 on 1 degrees of freedom
AIC: 30.755
Number of Fisher Scoring iterations: 3
となり、すべての係数が有意水準1％で0でないと判断されます。
\begin{align}
\log(n_{22}) &= \mbox{切片} + \mbox{民間療法なしの係数} + \mbox{風邪を引いていないの係数} \\
&= 4.3014-0.4055+0.3640 = 4.2599\\
e^{4.2599} &\approx 71
\end{align}
となります。結局、この観測度数の各セルは、
\begin{align}
\log(n_{ij}) &= \mbox{切片} + \mbox{民間療法}_{i} + \mbox{風邪}_{j}\ (i = 1,2\ j = 1,2)
\end{align}
のように民間療法の有無、風邪を引いているかどうかと切片の係数のみのモデルで表現できることがわかりました。つまり、交互作用は考慮しなくてもよいということです。
今回は2要因の分割表に対して対数線形モデルを適用させましたが、要因が増えた場合は、すべての組み合わせをモデル式に入れなければなりません。例えば、先ほどのデータに性別要因が加わった場合、分割表は、
| 風邪を引いた | 風邪を引かなかった | ||
|---|---|---|---|
| 女性 | 民間療法あり | 76 | 104 |
| 民間療法なし | 30 | 90 | |
| 男性 | 民間療法あり | 56 | 124 |
| 民間療法なし | 47 | 73 |
のようになり、飽和モデルを考えるなら、モデル式は
\begin{align}
\log(n_{ij}) &= \mbox{切片} + \mbox{民間療法}_{i} + \mbox{風邪}_{j} + \mbox{性別}_{k}\\
&+ \mbox{民間療法*風邪}_{ij} + \mbox{風邪*性別}_{jk} + \mbox{性別*民間療法}_{ki}\\
&+ \mbox{民間療法*風邪*性別}_{ijk}
\end{align}
となります。要因数をnとすると、切片を含めた飽和モデルの項数は
$$2^n$$
となります。指数のオーダーで増えていくという点に注意しましょう。
References
同志社大学のページ(対数線形モデル以外にもRでの実行例が数多く載っており参考になります。)
対数線形モデルですべての係数を表示させる方法 1