対数線形モデル
対数線形モデルにおいては、分割表における各セルの期待度数を予測する。すなわち、複数のカテゴリカル変数の中から特定の目的変数を設定するのでなく、全ての変数について相互の関係を調べることを目的としている。
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章での解析結果とも一致する。
AICとBICで結果が異なったが、データ数が多い場合にAICを基準とすると、フルモデルを選択しやすい傾向があり、BICをモデル選択に使用する場合が多い、とのことである。
