データは東大「人文社会科学の統計学」 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次の効果を外した同時独立モデルだが,これに有意差は無い。つまり,この項は「効いていない」 #(テキストでは) year の効果は無い。
social.mobility2.loglm3 <- loglm(~father + child + year + father:child + father:year + child:year, data=social.mobility2) anova(social.mobility2.loglm3, test="Chisq")
# 同じことだが,以下でもよい
social.mobility2.loglm3 <- loglm(social.mobility2, list(c(1), c(2), c(1,2), c(2,3), c(1,3)))
# さらに同じことだが,loglin() 関数を利用してもよい
social.mobility2.loglin3 <- loglin(social.mobility2 , list(c(1), c(2), 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 # を指定して含まれる