SNAGeek

Sociologically Technological, and Technologically Sociological

{gnm}を使ったクロス表の分析

  • 最近は対数線形・対数乗法モデルの復習をしています。
  • 色々と読み返してはいるのですが、以下の本は簡潔にまとまっており、非常に勉強になりました。
  • オッズ比から始まって、RC(M)-L modelまで取り扱っています。

Association Models (Quantitative Applications in the Social Sciences) by Raymond Sin-Kwok Wong(2010-02-09)

Association Models (Quantitative Applications in the Social Sciences) by Raymond Sin-Kwok Wong(2010-02-09)

  • 作者:Raymond Sin-Kwok Wong
  • 出版社/メーカー: SAGE Publications, Inc
  • 発売日: 1787
  • メディア: ペーパーバック

  • 今回は備忘録も兼ねて、同書のExample2.2をRの{gnm}パッケージを使って一部再現しようと思います。
    • Rでの対数線形・対数乗法モデルの実装方法に関しては、あまりまとまった日本語の資料がないようでした。
    • モデル自体の解説は日本語でもすぐに資料が見つかると思います。

github.com

データの用意

  • 同書に記載ある通り、学歴(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値(Χ2検定)、BICの4つを出力します。
fit_df <- function(gnm_object){
  df <- gnm_object$df.residual
  g_sq <- gnm_object$deviance
  p_value <- chi2.p.value(gnm_object)
  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.000000 1274.08092
2     U 11  244.0177150 0.000000  153.18077
3     R  9  205.9735836 0.000000  131.65245
4     C  8  155.3722685 0.000000   89.30903
5   RpC  6   91.6077617 0.000000   42.06034
6    RC  6  125.0597486 0.000000   75.51232
7   RC2  2    0.6001022 0.750266  -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 
  • RC effectの1次元目は学歴・職業階層の一次元的な”序列”を表していそうです。
    • つまり、高い(低い)学歴が高い(低い)職業的地位と結びついている
  • 一方で2次元目は解釈が難しく、筆者によればこれはメリトクラシーでは説明できない、労働市場の参入障壁を表しているといいます。
    • RC(M)モデルの Mを増やせば説明力は上がりますが、その解釈は分析者の腕の見せどころですね…。