{gnm}を使ったクロス表の分析
【2021/03/12更新：p値をカイ2乗値から計算していましたがdeviance由来に修正しました】
- 最近は対数線形・対数乗法モデルの復習をしています。
- 色々と読み返してはいるのですが、以下の本は簡潔にまとまっており、非常に勉強になりました。
- オッズ比から始まって、RC(M)-L modelまで取り扱っています。
- 今回は備忘録も兼ねて、同書のExample2.2をRの{gnm}パッケージを使って一部再現しようと思います。
- Rでの対数線形・対数乗法モデルの実装方法に関しては、あまりまとまった日本語の資料がないようでした。
- モデル自体の解説は日本語でもすぐに資料が見つかると思います。
データの用意
- 同書に記載ある通り、学歴(EDUC)と現職(OCC)のクロス表
ex.2.2.txtを用意します。 - 元ネタはGSS(1985-1990 cumulative, 女性のみ)のようです。
,上層ノンマニュアル,下層ノンマニュアル,上層マニュアル,下層マニュアル,農業 大卒以上,518,95,6,35,5 短大,81,67,4,49,2 高校,452,1003,67,630,5 高校未満,71,157,37,562,12
実行環境
> sessionInfo() R version 3.4.4 (2018-03-15) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 18.04.2 LTS Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1 locale: [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8 [5] LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 LC_PAPER=C.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] purrr_0.3.3 dplyr_0.8.3 gnm_1.1-0 loaded via a namespace (and not attached): [1] Rcpp_1.0.1 lattice_0.20-35 crayon_1.3.4 relimp_1.0-5 assertthat_0.2.1 MASS_7.3-51.5 [7] grid_3.4.4 R6_2.4.0 magrittr_1.5 pillar_1.4.2 rlang_0.4.2 qvcalc_1.0.1 [13] rstudioapi_0.10 Matrix_1.2-18 tools_3.4.4 glue_1.3.1 compiler_3.4.4 pkgconfig_2.0.2 [19] tidyselect_0.2.5 nnet_7.3-12 tibble_2.1.3
ライブラリ読み込み
library(gnm) library(dplyr) library(purrr)
前処理
- 何はともあれ、行ラベル・列ラベル・頻度からなるtidyなデータフレームにします。
tab <- read.table("ex.2.2.txt",header = TRUE,sep = ",",row.names = 1) mat <- as.matrix(tab) df <- as.data.frame.table(mat) colnames(df) <- c("EDUC","OCC","Freq")
> df EDUC OCC Freq 1 大卒以上 上層ノンマニュアル 518 2 短大 上層ノンマニュアル 81 3 高校 上層ノンマニュアル 452 4 高校未満 上層ノンマニュアル 71 5 大卒以上 下層ノンマニュアル 95 6 短大 下層ノンマニュアル 67 7 高校 下層ノンマニュアル 1003 8 高校未満 下層ノンマニュアル 157 9 大卒以上 上層マニュアル 6 10 短大 上層マニュアル 4 11 高校 上層マニュアル 67 12 高校未満 上層マニュアル 37 13 大卒以上 下層マニュアル 35 14 短大 下層マニュアル 49 15 高校 下層マニュアル 630 16 高校未満 下層マニュアル 562 17 大卒以上 農業 5 18 短大 農業 2 19 高校 農業 5 20 高校未満 農業 12
分析
- まずはクロス表全体としてオッズ比が1になる independent model を作ります。
- 一般化非線形モデルを推定するためのパッケージ{gnm}を使います。
> gnm(formula = Freq ~ EDUC + OCC, data = df, family = poisson) Call: gnm(formula = Freq ~ EDUC + OCC, family = poisson, data = df) Coefficients: (Intercept) EDUC短大 EDUC高校 EDUC高校未満 5.2557 -1.1775 1.1858 0.2415 OCC下層ノンマニュアル OCC上層マニュアル OCC下層マニュアル OCC農業 0.1640 -2.2867 0.1286 -3.8448 Deviance: 1373.176 Pearson chi-squared: 1455.658 Residual df: 12
- ここから階層的にモデルを作っていきます。
- RモデルやCモデルを作るためには、推定の問題上、パラメータを標準化し、足して0になるように制約をかける必要があります。
Row <- scale(as.numeric(row(mat)), scale = TRUE) Col <- scale(as.numeric(col(mat)), scale = TRUE)
- 本来であればR+C+RCモデルなども記載があるのですが、パラメータ制約がややこしいので省略しています。
- あとで
map()するため、結果をlistに放り込んでいきます。
set.seed(1024) gnm_list <- list() gnm_list$O <- gnm(formula = Freq ~ EDUC + OCC, data = df, family = poisson) gnm_list$U <- update(gnm_list$O, ~ . + Row:Col) gnm_list$R <- update(gnm_list$O, ~ . + EDUC:Col) gnm_list$C <- update(gnm_list$O, ~ . + OCC:Row) gnm_list$RpC <- update(gnm_list$O, ~ . + EDUC:Col + OCC:Row) #R+Cモデル gnm_list$RC <- update(gnm_list$O, ~ . + Mult(EDUC,OCC)) gnm_list$RC2 <- update(gnm_list$RC, ~ . + Mult(EDUC,OCC,inst = 2))
- さて、無事に各モデルが推定できたはずです。
モデルの適合度
- 次は各モデルの適合度を見ていきます。どのモデルが最もクロス表をよく説明するでしょうか。
- 自由度df、逸脱度G2(-2log尤度)、p値、BICの4つを出力します。
fit_df <- function(gnm_object){ df <- gnm_object$df.residual g_sq <- gnm_object$deviance p_value <- 1 - pchisq(q = gnm_object$deviance ,df = gnm_object$df.residual) bic <- g_sq - gnm_object$df.residual * log(sum(gnm_object$data$Freq)) fit_df <- data.frame(df,g_sq,p_value,bic) return(fit_df) } gof <- bind_rows(map(gnm_list,fit_df),.id = "model")
- 結果を見てみます。
> gof model df g_sq p_value bic 1 O 12 1373.1757657 0.0000000 1274.08092 2 U 11 244.0177150 0.0000000 153.18077 3 R 9 205.9735836 0.0000000 131.65245 4 C 8 155.3722685 0.0000000 89.30903 5 RpC 6 91.6077617 0.0000000 42.06034 6 RC 6 125.0597486 0.0000000 75.51232 7 RC2 2 0.6001022 0.7407804 -15.91571
最後のRC(2)モデルが最もBICが低く、かつp値も大きいため、このモデルが選択されます。
- p値が大きいというのは、Χ2値が小さい、つまり、RC2モデルによって求められる各セルの期待度数が、元データからそれほど乖離していないことを示しています。
- 本来は、自由度をいくつ使ってG2が改善したかなども併せて見ていく必要がありますが…。
RC(2)モデルのパラメータはこのような感じになります。
> gnm_list$RC2 Call: gnm(formula = Freq ~ EDUC + OCC + Mult(EDUC, OCC) + Mult(EDUC, OCC, inst = 2), family = poisson, data = df) Coefficients: (Intercept) EDUC短大 5.066410 -0.916769 EDUC高校 EDUC高校未満 1.067655 0.217877 OCC下層ノンマニュアル OCC上層マニュアル 0.353842 -2.143860 OCC下層マニュアル OCC農業 0.129995 -3.973502 Mult(., OCC).EDUC大卒以上 Mult(., OCC).EDUC短大 -5.581051 -1.462538 Mult(., OCC).EDUC高校 Mult(., OCC).EDUC高校未満 0.975465 1.973569 Mult(EDUC, .).OCC上層ノンマニュアル Mult(EDUC, .).OCC下層ノンマニュアル -0.249405 0.099685 Mult(EDUC, .).OCC上層マニュアル Mult(EDUC, .).OCC下層マニュアル 0.222067 0.314239 Mult(EDUC, .).OCC農業 Mult(., OCC, inst = 2).EDUC大卒以上 -0.013403 1.240632 Mult(., OCC, inst = 2).EDUC短大 Mult(., OCC, inst = 2).EDUC高校 0.683454 -1.323081 Mult(., OCC, inst = 2).EDUC高校未満 Mult(EDUC, ., inst = 2).OCC上層ノンマニュアル 3.154406 -0.168216 Mult(EDUC, ., inst = 2).OCC下層ノンマニュアル Mult(EDUC, ., inst = 2).OCC上層マニュアル -0.246388 0.007613 Mult(EDUC, ., inst = 2).OCC下層マニュアル Mult(EDUC, ., inst = 2).OCC農業 0.094214 0.382994 Deviance: 0.6001022 Pearson chi-squared: 0.5746549 Residual df: 2
