Chapter 4 多元クロス集計表の処理
つぎに、変数が3つ以上の多元クロス集計表の処理について見ていく。例えば、夫婦の学歴組み合わせのクロス集計表が複数の調査年について得られている状況を想定してみよう。この場合、
妻の学歴 × 夫の学歴 × 調査年
の3変数からなる多元クロス集計表が得られることになる。こうした多元クロス集計表をRの中で表現・処理する方法はいくつかあるが、以下では配列（array）を用いた方法について紹介する。
4.1 配列（array）
配列（array）とは行列（matrix）を三次元以上に一般化したものである。言い換えると、二次元の配列は行列として処理することができる。Rで配列を作成する方法は主に2つある。
- 関数
array()を用いて配列を直接作成する - 複数の行列を
abind()で結合する
以下、順に見ていこう。まずは関数array()の使用方法についてである。
# 2×2×2の3次元からなり、要素がすべてNAの空の配列
array(data = NA, dim = c(2,2,2))## , , 1
##
## [,1] [,2]
## [1,] NA NA
## [2,] NA NA
##
## , , 2
##
## [,1] [,2]
## [1,] NA NA
## [2,] NA NA# 任意のベクトルを引き渡してもよい。列方向に順番に格納されていく。
array(data = 1:8, dim = c(2,2,2))## , , 1
##
## [,1] [,2]
## [1,] 1 3
## [2,] 2 4
##
## , , 2
##
## [,1] [,2]
## [1,] 5 7
## [2,] 6 8# データとして引き渡したベクトルの長さが次元の長さの積と等しくない場合、配列が埋まるまでデータのベクトルが反復される
array(data = 1:3, dim = c(2,2,2))## , , 1
##
## [,1] [,2]
## [1,] 1 3
## [2,] 2 1
##
## , , 2
##
## [,1] [,2]
## [1,] 2 1
## [2,] 3 2ただし、報告者の個人的な経験からすると、array()を用いて分析対象とするクロス集計表を直接作成することはあまりない。むしろ、array()が活躍するのは以下のような場面である。
- 既存の配列の次元の長さを増やす
- 仮想的な配列のデータを作成して、自身が作成した関数が意図した通りの挙動で動くかを確認する
4.2 abind()を用いた配列の作成
クロス集計表がraw dataとして与えられている場合、それらをもとに多元クロス集計表を表す配列を作成する方法は、複数のクロス集計表を行列として読み込み、それらを結合することである。その際に用いる関数がabind()である。以下、この関数の使い方を見ていこう。
はじめに問題になるのは複数のクロス集計表を読み込む方法である。最も愚直な方法は、CSVファイルを読み込むためのコードを調査年ごとに書くことである。
EduWH1980_raw <- read.csv(paste0(CensusDir, "EduWH1980.csv"), header = FALSE)
EduWH1990_raw <- read.csv(paste0(CensusDir, "EduWH1990.csv"), header = FALSE)
EduWH1980 <- as.matrix(EduWH1980_raw)
EduWH1990 <- as.matrix(EduWH1990_raw)
dimnames(EduWH1980) <-
dimnames(EduWH1990) <- list(c("JHS", "HS", "VC/JC", "UNI"),
c("JHS", "HS", "VC/JC", "UNI"))
names(dimnames(EduWH1980)) <-
names(dimnames(EduWH1990)) <- c("EduW", "EduH")読み込むデータファイルが数個程度であればこの方法でも大きな問題は生じないが、歴史の長い反服横断調査やパネル調査ほど読み込むべきデータファイルの数は増えていく。その場合、調査回ごとにデータ読み込みのプログラムを書くのではなく、データ読み込みのための自前の関数を作ってしまうと便利である。例えば、以下ではreadEduWH()という関数を定義している。
readEduWH <- function(path, census_year){
TargetCSV <- paste0(path, "EduWH", as.character(census_year), ".csv") # 例えば1980年のデータファイル名は"EduWH1980.csv"にしておく
EduWH_raw <- read.csv(TargetCSV, header = FALSE)
EduWH_mat <- as.matrix(EduWH_raw)
EduWH <- EduWH_mat
dimnames(EduWH) <- list(c("JHS", "HS", "VC/JC", "UNI"),
c("JHS", "HS", "VC/JC", "UNI"))
return(EduWH)
}この関数は引数としてcensus_yearを持ち、ここに読み込みたい国勢調査の調査年を引き渡すと、その調査年の夫婦の学歴組み合わせのCSVファイルが読み込まれて（ただし、そのためにはCSVのファイル名に調査年を含めておく必要がある）返り値として出力される。例えば、1980年の国勢調査データを読み込むためには、
readEduWH(path = CensusDir, census_year = 1980)## JHS HS VC/JC UNI
## JHS 1825957.2 692306.05 25961.477 43269.13
## HS 744229.0 3106723.41 164422.687 856728.74
## VC/JC 34615.3 207691.82 69230.605 493268.06
## UNI 0.0 25961.48 8653.826 337499.20とすればよい。この自作関数とabind()とを組み合わせることで、以下のように1980年から2010年までの4回の国勢調査データをひとつの配列に格納することができる。
library(abind)
CensusYearList <- seq(1980, 2010, 10) # 国勢調査の調査年が格納されたベクトル
EduWH <- NULL # 空のオブジェクトを作っておく
for(iYear in CensusYearList){
EduWH_temp <- readEduWH(path = CensusDir, census_year = iYear) # iYear年の国勢調査のデータを読み込む
EduWH <- abind(EduWH, EduWH_temp, along = 3) # その結果を事前に作っておいたEduWHに順次結合していく
}
dimnames(EduWH)[[3]] <- CensusYearList # abindで結合した次元にはラベルがついていないのでつけておくと便利
names(dimnames(EduWH)) <- c("EduW", "EduH", "CensusYear") # 各次元の変数名を付与オブジェクトEduWHに4回の国勢調査における夫婦の学歴組み合わせのクロス集計表が保存されているか確認してみよう。
EduWH## , , CensusYear = 1980
##
## EduH
## EduW JHS HS VC/JC UNI
## JHS 1825957.2 692306.05 25961.477 43269.13
## HS 744229.0 3106723.41 164422.687 856728.74
## VC/JC 34615.3 207691.82 69230.605 493268.06
## UNI 0.0 25961.48 8653.826 337499.20
##
## , , CensusYear = 1990
##
## EduH
## EduW JHS HS VC/JC UNI
## JHS 434278.720 244281.78 13571.21 20356.81
## HS 474992.350 2497102.64 189996.94 712488.53
## VC/JC 61070.445 474992.35 162854.52 821058.20
## UNI 6785.605 61070.44 13571.21 597133.24
##
## , , CensusYear = 2000
##
## EduH
## EduW JHS HS VC/JC UNI
## JHS 91394.56 114243.20 11424.32 17136.48
## HS 222774.24 1902149.28 228486.40 605488.96
## VC/JC 45697.28 542655.20 274183.68 925369.92
## UNI 5712.16 74258.08 34272.96 611201.12
##
## , , CensusYear = 2010
##
## EduH
## EduW JHS HS VC/JC UNI
## JHS 48150.360 72225.54 9630.072 14445.11
## HS 134821.008 1136348.50 211861.584 361127.70
## VC/JC 52965.396 587434.39 390017.916 789665.90
## UNI 9630.072 134821.01 86670.648 775220.80このように、妻の学歴（EduW）×夫の学歴（EduH）×調査年（CensusYear）からなる三元クロス表が配列として表現できたことがわかる。
4.3 apply()を用いた配列の処理
apply()は、行列あるいは配列の特定の次元ごとに何らかの処理を並列的に施すための関数である。例えば、夫婦の学歴組み合わせの分布（すなわち全体パーセント）を調査年ごとに算出するためには、
apply(EduWH, MARGIN = 3, FUN = proportions, simplify = FALSE) # "MARGIN =", "FUN = "の部分は省略してもよい## $`1980`
## EduH
## EduW JHS HS VC/JC UNI
## JHS 0.211422846 0.080160321 0.003006012 0.00501002
## HS 0.086172345 0.359719439 0.019038076 0.09919840
## VC/JC 0.004008016 0.024048096 0.008016032 0.05711423
## UNI 0.000000000 0.003006012 0.001002004 0.03907816
##
## $`1990`
## EduH
## EduW JHS HS VC/JC UNI
## JHS 0.064 0.036 0.002 0.003
## HS 0.070 0.368 0.028 0.105
## VC/JC 0.009 0.070 0.024 0.121
## UNI 0.001 0.009 0.002 0.088
##
## $`2000`
## EduH
## EduW JHS HS VC/JC UNI
## JHS 0.016016016 0.02002002 0.002002002 0.003003003
## HS 0.039039039 0.33333333 0.040040040 0.106106106
## VC/JC 0.008008008 0.09509510 0.048048048 0.162162162
## UNI 0.001001001 0.01301301 0.006006006 0.107107107
##
## $`2010`
## EduH
## EduW JHS HS VC/JC UNI
## JHS 0.010 0.015 0.002 0.003
## HS 0.028 0.236 0.044 0.075
## VC/JC 0.011 0.122 0.081 0.164
## UNI 0.002 0.028 0.018 0.161とすればよい。なお、配列の次元にラベルをつけている場合、そのラベルを用いて引数MARGINに引き渡すこともできる。
apply(EduWH, "CensusYear", proportions, simplify = FALSE)## $`1980`
## EduH
## EduW JHS HS VC/JC UNI
## JHS 0.211422846 0.080160321 0.003006012 0.00501002
## HS 0.086172345 0.359719439 0.019038076 0.09919840
## VC/JC 0.004008016 0.024048096 0.008016032 0.05711423
## UNI 0.000000000 0.003006012 0.001002004 0.03907816
##
## $`1990`
## EduH
## EduW JHS HS VC/JC UNI
## JHS 0.064 0.036 0.002 0.003
## HS 0.070 0.368 0.028 0.105
## VC/JC 0.009 0.070 0.024 0.121
## UNI 0.001 0.009 0.002 0.088
##
## $`2000`
## EduH
## EduW JHS HS VC/JC UNI
## JHS 0.016016016 0.02002002 0.002002002 0.003003003
## HS 0.039039039 0.33333333 0.040040040 0.106106106
## VC/JC 0.008008008 0.09509510 0.048048048 0.162162162
## UNI 0.001001001 0.01301301 0.006006006 0.107107107
##
## $`2010`
## EduH
## EduW JHS HS VC/JC UNI
## JHS 0.010 0.015 0.002 0.003
## HS 0.028 0.236 0.044 0.075
## VC/JC 0.011 0.122 0.081 0.164
## UNI 0.002 0.028 0.018 0.161上の出力結果を見ると調査年の前に$がついているが、これは返り値がリストになっていることを示している。$の後に続く文字列がリストの要素の名前になっているのでこれを用いて要素にアクセスできる。
TotProp <- apply(EduWH, 3, proportions, simplify = FALSE)
class(TotProp)## [1] "list"TotProp$`1980`## EduH
## EduW JHS HS VC/JC UNI
## JHS 0.211422846 0.080160321 0.003006012 0.00501002
## HS 0.086172345 0.359719439 0.019038076 0.09919840
## VC/JC 0.004008016 0.024048096 0.008016032 0.05711423
## UNI 0.000000000 0.003006012 0.001002004 0.039078164.4 apply()と自作関数を組み合わせる
最後に、同類婚、すなわち夫婦の学歴が同じカテゴリの組み合わせ、の割合を調査年ごとに算出してみよう。一例として、これは以下の3ステップに分けて考えればよい。
- 任意の調査年の夫婦の学歴組み合わせのクロス表に対して、対角セルを取り出す
- 対角セルの各度数を全体度数で割る
- 1-2を全ての調査年について繰り返す
まず、例として1980年のクロス表を対象にすると、1および2は以下のコードで実行できる。
ObsDiag <- diag(EduWH1980) # diag():行列の対角成分を取得
PropDiag <- ObsDiag / sum(EduWH1980)問題はこうした処理を全ての調査年ごとに繰り返す（ステップ3）にはどうすればよいか、である。今回の例では非常に単純な処理なので、同じコードを調査年ごとに書き起こしても大きな問題はなさそうである。しかしながら、プログラミングの有名な原則「DRY（Don’t Repeat Yourself）原則」にも言われているように、できるだけできるだけコードの重複は避ける習慣はつけておいた方がよい（自戒）。そこで、上記のステップ1-2の処理を行う関数を作成してみよう。
prop_diag <- function(x){
ObsDiag <- diag(x)
PropDiag <- ObsDiag / sum(x)
return(PropDiag)
}これにより、引数として夫婦の学歴組み合わせからなる行列を引き渡すと、対角成分の比率が返り値として出力される関数が定義された。この自作関数とapply()を組み合わせることで、調査年ごとに同類婚の占める割合を算出することができる。
PropDiag_byCensus <- apply(EduWH, 3, prop_diag)
PropDiag_byCensus## CensusYear
## EduW 1980 1990 2000 2010
## JHS 0.211422846 0.064 0.01601602 0.010
## HS 0.359719439 0.368 0.33333333 0.236
## VC/JC 0.008016032 0.024 0.04804805 0.081
## UNI 0.039078156 0.088 0.10710711 0.161この場合、オプションのsimplify = FALSEを外しておくと、返り値が行列となり、その後の処理（出力結果を可視化する、CSVに書き出すなど）も容易になる。出力結果を行列ではなくてリストで得たい場合は、simplify = FALSEをつけておくとよい。