R_対数線形モデル のバックアップ(No.2) - アールメカブ

アールメカブ


R_対数線形モデル のバックアップ(No.2)


_ 対数線形モデル

データは東大「人文社会科学の統計学」 p.277 飽和モデル,独立モデル

social.mobility2 <- 
 array(c(
   517, 244, 250,
   152, 351, 328,
   7, 13, 140,
   441, 159, 272,
   130, 227, 351,
   24, 21, 268
   ),
   dim = c(3,3, 2),
   dimnames =
   list(
    father = c("white","blue","agri" ), # = c(1) 行   
    child = c( "white","blue","agri"),  # = c(2) 列
    year= c("1985", "1975")             # = c(3) 次元
    )
   )
 class(social.mobility2 ) <- "table"
 ftable(social.mobility2)
 library(MASS)

# 3次の効果を含めた飽和モデルは,自由度ゼロなので必ず適合する。p. 271

social.mobility2.loglm <-
   loglm(~father + child + year +
      father:child + father:year + child:year +
      child:father:year,
    data=social.mobility2)
anova(social.mobility2.loglm, test="Chisq")

# 以下は,3次の効果を外した同時独立モデルだが,これに有意差は無い。つまり,この項は「効いていない」.そして,この項以外は削除できない.

stepAIC(social.mobility2.loglm)

#(テキストによれば) child と father の各組ごとの組み合わせに年代による差は無い。

 social.mobility2.loglm3 <-
  loglm(~father + child + year
   + father:child + father:year + child:year, data=social.mobility2)
 anova(social.mobility2.loglm3, test="Chisq")

# さらに同じことだが,loglin() 関数を利用してもよい

social.mobility2.loglin3 <- loglin(social.mobility2 ,
  list(c(1), c(2),c(3),  c(1,2), c(2,3), c(1,3)))

# あるいは交互作用に関わる項は暗黙のうちに含まれるので

social.mobility2.loglm3 <- loglm(social.mobility2,
   list(c(1,2), c(2,3), c(1,3)))    
# anova(social.mobility2.loglin3, test="Chisq")
# summary(social.mobility2.loglin3)
1 - pchisq(social.mobility2.loglin3$lrt,
  social.mobility2.loglin3$df)
1 - pchisq(social.mobility2.loglin3$pearson,
   social.mobility2.loglin3$df)#
# 東大「人文社会の統計学」 p.276
# log m = u... + ui.. + u.j. + u..k + uij. + ui.k + u.jk (9.33)
# ここで uij.  ui.k  u.jk がそれぞれ,father:child, father:year, child:year
# を表しているが
# uijk (p.276下から12行目), すなわち father:child:year
# は含まれていないモデル
# なお u...  は全体平均 grandmean は暗黙のうちに含まれる
# ui.. 行平均  u.j. 列平均 u..k  軸平均は それぞれ要因 father, child, year
# を指定して含まれる

対数線形モデルは対応のないことが前提とされている. 参考

http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc036/08918.html

2*2*5の3要因の度数の比較を行いたいのですが,一つの要因に対応があるような計画です。要因が多いので,対数線形モデル分析のように思えるのですが,対数線形モデル分析はすべての要因に対応がないことが前提ですよね。このような場合,どのような分析が適切なのでしょうか・・・。

SAS なら Proc catmod, proc genmod でできるようです。例えばgenmod の次の解説参照

http://www.id.unizh.ch/software/unix/statmath/sas/sasdoc/stat/chap29/sect1.htm

ここでGEEについて触れてますが,そうするとGEEのできるソフトならできるでしょう。例えばSTATA など。catmod の説明は

http://www.id.unizh.ch/software/unix/statmath/sas/sasdoc/stat/chap22/sect1.htm