きわめて個人的なR言語のメモ

2005-01-22Rで対数線形モデル

近日中に全面修正します。

家事に対する意見

 放送大学の教材「社会調査法」(1992年初版)に載っている対数線形モデル?の例題を、に搭載されている対数線形モデルのための関数loglinを使って分析してみる。

 「「家事は女の仕事である」という意見に対する賛否」を「結婚時の就職状況」と「現在の就職状況」から分類したもの。

 既存の多重クロス集計表をオブジェクト"shigoto"に入力する。オブジェクトの名前は何でもいい。

shigoto <- array(data=c(7,5,14,23,10,1,16,11),dim=c(2,2,2))

 各カテゴリに名前を付ける。日本語対応版のRなら日本語のカテゴリ名を付けられます。

dimnames(shigoto) <- list(c("結婚時有職","結婚時無職"),c("現在有職","現在無職"),c("賛成","反対"))

 "shigoto"の内容を表示させる。

shigoto

 見やすい形の多重クロス集計表に整形して表示させる。

ftable(shigoto)

いよいよloglin

さまざまなモデルでログリニア分析を行い、モデルを選択する。

ここで、liglin関数の二つ目の引数でモデルを指定するのですが、この例題の場合、1=「結婚時の就職状況(有職/無職)」、2=「現在の就職状況(有職/無職)」、3=「「家事は女の仕事」に対する賛否(賛成/反対)」を示しています。また、loglin関数では階層的なモデルを想定しており、高次のモデルには下位のモデルが含まれます。つまり、 list(c(A,B)) は A*B という二次の交互作用と、AとBそれぞれの主効果の両方を含んでおり、模式的に表すと【A+B+A*B+e】という意味になります。

飽和モデル

 三変数が全て関連しているモデル。

 【A+B+C+A*B+B*C+C*A+A*B*C+e】

model0 <- loglin(shigoto, list(c(1,2,3)), param=TRUE)

対連関モデル(均一連関モデル)

 飽和モデルから三次の交互作用を取り除いたモデル。

 【A+B+C+A*B+B*C+C*A+e】

model1 <- loglin(shigoto, list(c(1, 2), c(1, 3), c(2, 3)), param=TRUE)

条件付独立モデル

 対連関モデルから特定の二変数の交互作用を取り除いたモデル。

 【A+B+C+A*B+B*C+e】

model2 <- loglin(shigoto, list(c(1, 2), c(1, 3)), param=TRUE)
model3 <- loglin(shigoto, list(c(1, 2), c(2, 3)), param=TRUE)
model4 <- loglin(shigoto, list(c(1, 3), c(2, 3)), param=TRUE)

一変数独立モデル

 変数Cが変数Aと変数Bの組み合わせから独立であるモデル。

 【A+B+C+A*B+e】

model5 <- loglin(shigoto, list(c(1, 2),3), param=TRUE)
model6 <- loglin(shigoto, list(c(1, 3),2), param=TRUE)
model7 <- loglin(shigoto, list(c(2, 3),1), param=TRUE)

三変数独立モデル

 三つの変数がそれぞれ独立しているモデル。(χ二乗検定が想定しているモデル。)

 【A+B+C+e】

model8 <- loglin(shigoto, list(1, 2, 3), param=TRUE)

分析結果を表示する。

 出力のうち、$lrtは尤度比、$pearsonはχ二乗値、$dfは自由度、$marginは出力されたモデル。

model0
model1
model2
model3
model4
model5
model6
model7
model8

 各モデルにおける観測度数と推定度数の適合度の検定の有意水準(p値)を計算する。

p0 <- 1-pchisq(model0$lrt, model0$df)
p1 <- 1-pchisq(model1$lrt, model1$df)
p2 <- 1-pchisq(model2$lrt, model2$df)
p3 <- 1-pchisq(model3$lrt, model3$df)
p4 <- 1-pchisq(model4$lrt, model4$df)
p5 <- 1-pchisq(model5$lrt, model5$df)
p6 <- 1-pchisq(model6$lrt, model6$df)
p7 <- 1-pchisq(model7$lrt, model7$df)
p8 <- 1-pchisq(model8$lrt, model8$df)

 計算したp値を表示させます。

 0.05以下のモデルは当てはまりが悪いということになる(んだろうと思う)。

p0
p1
p2
p3
p4
p5
p6
p7
p8

 堀先生のページを読むと、AIC=χ2-2*df ということなので、これも算出してみます。AICが最も小さいモデルを選択することになります。ただし、「2.0以内の違いは意味のある違いと解釈しない方がいい」ということです。

AIC0 <- model0$pearson-2*model0$df
AIC1 <- model1$pearson-2*model1$df
AIC2 <- model2$pearson-2*model2$df
AIC3 <- model3$pearson-2*model3$df
AIC4 <- model4$pearson-2*model4$df
AIC5 <- model5$pearson-2*model5$df
AIC6 <- model6$pearson-2*model6$df
AIC7 <- model7$pearson-2*model7$df
AIC8 <- model8$pearson-2*model8$df

 計算したAICを表示させます。

AIC0
AIC1
AIC2
AIC3
AIC4
AIC5
AIC6
AIC7
AIC8

 ついでにモザイク図を作成し、指定したモデルでの残差で色を付ける。

 コピペすると一気にグラフが表示されますが、「PageUp」と「PageDown」キーで前後のグラフを見ることが出来ます。見られない場合は、グラフのウィンドウのメニューからHistory→Recordingと選択して、もう一度コマンドを入れてください。

mosaicplot(shigoto, main="model0飽和モデル", margin=list(c(1,2,3)), shade=TRUE)
mosaicplot(shigoto, main="model1対連関モデル", margin=list(c(1, 2), c(2, 3), c(1, 3)), shade=TRUE)
mosaicplot(shigoto, main="model2条件付独立モデル", margin=list(c(1, 2), c(1, 3)), shade=TRUE)
mosaicplot(shigoto, main="model3条件付独立モデル", margin=list(c(1, 2), c(2, 3)), shade=TRUE)
mosaicplot(shigoto, main="model4条件付独立モデル", margin=list(c(1, 3), c(2, 3)), shade=TRUE)
mosaicplot(shigoto, main="model5一変数独立モデル", margin=list(c(1, 2),3), shade=TRUE)
mosaicplot(shigoto, main="model6一変数独立モデル", margin=list(c(1, 3),2), shade=TRUE)
mosaicplot(shigoto, main="model7一変数独立モデル", margin=list(c(2, 3),1), shade=TRUE)
mosaicplot(shigoto, main="model8三変数独立モデル", margin=list(1, 2, 3), shade=TRUE)

最も当てはまりが良さそうなmodel2のモザイク図がこれ。残差が小さいので白いです。

f:id:bob3:20050123222855:image

#ここまでは出来ましたが、結果をどう解釈すればよいのか解りません(^^;;;

トラックバック - http://r-user.g.hatena.ne.jp/bob3/20050122