読者です 読者をやめる 読者になる 読者になる

カテゴリカルデータ解析読書メモ(第8章)

対数線形モデル

対数線形モデルにおいては、分割表における各セルの期待度数を予測する。すなわち、複数のカテゴリカル変数の中から特定の目的変数を設定するのでなく、全ての変数について相互の関係を調べることを目的としている。

lonlin()とloglm()

対数線形モデルによる解析を行う関数としてloglin()とloglm()が用意されている。
loglin()は以下のように使用する。

> library(vcd)
> tabledata <- xtabs(Freq ~ Status + Urban, data = DanishWelfare)
> result <- loglin(tabledata, margin = list(1, 2), param = TRUE, fit = TRUE)
2 iterations: deviation 2.273737e-13 
> result$para
$`(Intercept)`
[1] 5.383959

$Status
     Widow    Married  Unmarried 
-0.7693390  1.0654469 -0.2961079 

$Urban
   Copenhagen SubCopenhagen     LargeCity          City       Country 
   -0.4836304    -0.3771835    -0.4102991     0.6787275     0.5923855 

対してloglm()は次のようにモデル式で入力する。

> result <- loglm(~ Status + Urban, data = tabledata)
> result
Call:
loglm(formula = ~Status + Urban, data = tabledata)

Statistics:
                      X^2 df P(> X^2)
Likelihood Ratio 153.9505  8        0
Pearson          158.1145  8        0
> coef(result)
$`(Intercept)`
[1] 5.383959

$Status
     Widow    Married  Unmarried 
-0.7693390  1.0654469 -0.2961079 

$Urban
   Copenhagen SubCopenhagen     LargeCity          City       Country 
   -0.4836304    -0.3771835    -0.4102991     0.6787275     0.5923855 

計算される値はどちらも同じだが、入出力の理解のしやすさはloglm()の方が上の様に思う。以降、loglm()を使用する。

期待度数と残差を求める

期待度数はfitted()、残差はresiduals()を用いて計算できる。

> round(fitted(result), 1)
Re-fitting to get fitted values
           Urban
Status      Copenhagen SubCopenhagen LargeCity   City Country
  Widow           62.2          69.2      67.0  199.0   182.5
  Married        389.9         433.6     419.5 1246.5  1143.4
  Unmarried       99.9         111.1     107.5  319.4   293.0
> round(residuals(result), 1)
Re-fitting to get frequencies and fitted values
           Urban
Status      Copenhagen SubCopenhagen LargeCity City Country
  Widow            6.1           2.5      -1.5  1.0    -6.3
  Married         -5.4          -0.2       0.0  0.0     3.1
  Unmarried        4.5          -1.8       1.2 -0.8    -1.7

三元表での対数線形モデル

変数が3つ以上になると作成できるモデルの数がかなり多くなる。
交互作用項について細かく検討するのであればupdate()を利用すると入力の手間を節約できる。

> Modelindependent <- loglm(~ Gender + Admit + Dept, data = UCBAdmissions)
> Modelindependent
Call:
loglm(formula = ~Gender + Admit + Dept, data = UCBAdmissions)

Statistics:
                      X^2 df P(> X^2)
Likelihood Ratio 2097.671 16        0
Pearson          2000.328 16        0
> # GenderとAdmitの交互作用を加える
> ModelGA <- update(Modelindependent, .~. + Gender:Admit)
> ModelGA
Call:
loglm(formula = . ~ Gender + Admit + Dept + Gender:Admit, data = UCBAdmissions)

Statistics:
                      X^2 df P(> X^2)
Likelihood Ratio 2004.222 15        0
Pearson          1748.160 15        0

テキストでは全てのモデルを作成して比較しているが、面倒なのでstep()を使用してモデル選択をしてみよう。
まずフルモデルを作成する。

> ModelFull <- loglm(~ Gender * Admit * Dept, data = UCBAdmissions)
> ModelFull
Call:
loglm(formula = ~Gender * Admit * Dept, data = UCBAdmissions)

Statistics:
                 X^2 df P(> X^2)
Likelihood Ratio   0  0        1
Pearson            0  0        1

フルモデルでは期待度数と観測度数が一致する。
では、AICによるモデル選択を試みる。

> step(ModelFull)
Start:  AIC=48
~Gender * Admit * Dept

                    Df    AIC
<none>                 48.000
- Gender:Admit:Dept  5 58.204
Call:
loglm(formula = ~Gender * Admit * Dept, data = UCBAdmissions)

Statistics:
                 X^2 df P(> X^2)
Likelihood Ratio   0  0        1
Pearson            0  0        1

AICではフルモデルが選択された。
次にBICによるモデル選択を試みる。

> step(ModelFull, k = log(sum(UCBAdmissions)))
Start:  AIC=202.02
~Gender * Admit * Dept

                    Df    AIC
- Gender:Admit:Dept  5 180.14
<none>                 202.02

Step:  AIC=180.14
~Gender + Admit + Dept + Gender:Admit + Gender:Dept + Admit:Dept

               Df     AIC
- Gender:Admit  1  173.25
<none>             180.14
- Admit:Dept    5  901.45
- Gender:Dept   5 1266.75

Step:  AIC=173.25
~Gender + Admit + Dept + Gender:Dept + Admit:Dept

              Df     AIC
<none>            173.25
- Admit:Dept   5  986.49
- Gender:Dept  5 1351.78
Call:
loglm(formula = ~Gender + Admit + Dept + Gender:Dept + Admit:Dept, 
    data = UCBAdmissions, evaluate = FALSE)

Statistics:
                      X^2 df    P(> X^2)
Likelihood Ratio 21.73551  6 0.001351993
Pearson          19.93841  6 0.002840164

BICを基準にした場合は、Admitとgenderの相互作用がモデルに含まれていない。これは第5章での解析結果とも一致する。
AICBICで結果が異なったが、データ数が多い場合にAICを基準とすると、フルモデルを選択しやすい傾向があり、BICをモデル選択に使用する場合が多い、とのことである。