クロス集計からモンテカルロシュミレーションまで
UPDATE: 2023-10-08 17:37:15
はじめに
ここでは、基本的なクロス集計の見方から多項分布を用いたクロス集計のモンテカルロシュミレーションまでをまとめている。また、Rでクロス集計を行うためのパッケージについてもまとめておく。
クロス集計のためのパッケージ
クロス集計の見方を説明する前に、クロス集計を行うためのパッケージをまとめておく。まずはサンプルデータの準備しておく。大阪と東京の人を対象に、たこ焼きの好き嫌いを聞いたアンケートみたいなイメージ。大阪のほうがたこ焼きが好きな割合を多くして、東京の方は意図的に少なくしておく。
library(tidyverse)
set.seed(1989)
n_tokyo <- 200
n_osaka <- 100
x_tokyo <- rbinom(n = n_tokyo, size = 1, prob = 0.4)
x_osaka <- rbinom(n = n_osaka, size = 1, prob = 0.7)
like <- c(x_tokyo, x_osaka)
like[like == 1] <- 'Yes'
like[like == 0] <- 'No'
df <- tibble(
area = rep(c('tokyo', 'osaka'), times = c(n_tokyo, n_osaka)),
is_like = like
)
df## # A tibble: 300 × 2
## area is_like
## <chr> <chr>
## 1 tokyo Yes
## 2 tokyo No
## 3 tokyo Yes
## 4 tokyo Yes
## 5 tokyo No
## 6 tokyo No
## 7 tokyo No
## 8 tokyo No
## 9 tokyo No
## 10 tokyo No
## # ℹ 290 more rowsクロス集計といえば、まずはtable()関数。ベクトルを渡せばクロス集計してくれる。
table(df$area, df$is_like)##
## No Yes
## osaka 29 71
## tokyo 123 77prop.table()関数を使えば割合で表示してくれる。横％、縦％、全体％かは`marginで指定できる。
prop.table(table(df$area, df$is_like), margin = 1) # 横％ ##
## No Yes
## osaka 0.290 0.710
## tokyo 0.615 0.385他にはxtabs()関数がある。引数をモデル式(formula)で渡す。
# prop.tableの引数に下記を渡すことで、割合は計算可能
xtabs(~ area + is_like, df)## is_like
## area No Yes
## osaka 29 71
## tokyo 123 77他にはftable()関数がある。引数をモデル式(formula)で渡す。多層のクロス集計をするときは便利。
# prop.tableの引数に下記を渡すことで、割合は計算可能
ftable(df, row.vars = 1, col.vars = 2)## is_like No Yes
## area
## osaka 29 71
## tokyo 123 77クロス集計用の関数を使わなくても集計するだけなので、パイプをつなげることでクロス集計はできる。
df %>%
group_by(area, is_like) %>%
# tally()でもよい
count() %>%
pivot_wider(names_from = is_like, values_from = n)## # A tibble: 2 × 3
## # Groups: area [2]
## area No Yes
## <chr> <int> <int>
## 1 osaka 29 71
## 2 tokyo 123 77janitorパッケージを利用すればデータフレームをベースにクロス集計が作成できる。詳細は下記のドキュメントを参照願います。
library(janitor)
df %>%
tabyl(area, is_like) %>%
adorn_totals(c('row', 'col'))## area No Yes Total
## osaka 29 71 100
## tokyo 123 77 200
## Total 152 148 300adorn_*()関数を繋げていくことで、出力形式を柔軟に変更できる。
df %>%
tabyl(area, is_like) %>%
adorn_totals(c('row', 'col')) %>%
adorn_percentages('row') %>%
adorn_pct_formatting(digits = 1)## area No Yes Total
## osaka 29.0% 71.0% 100.0%
## tokyo 61.5% 38.5% 100.0%
## Total 50.7% 49.3% 100.0%実数とパーセントを併記することも可能。
ct <- df %>%
tabyl(area, is_like) %>%
adorn_totals(where = c('row', 'col'))%>%
adorn_percentages('row') %>%
adorn_pct_formatting(digits = 1) %>%
# 'count (%)'表記
adorn_ns(position = 'front') %>%
adorn_title(
row_name = 'Area',
col_name = 'IsLike'
)
ct## IsLike
## Area No Yes Total
## osaka 29 (29.0%) 71 (71.0%) 100 (100.0%)
## tokyo 123 (61.5%) 77 (38.5%) 200 (100.0%)
## Total 152 (50.7%) 148 (49.3%) 300 (100.0%)データフレームがベースなので、そのままggplotに渡せば可視化も簡単にできる。
df %>%
tabyl(area, is_like) %>%
adorn_totals(where = 'row')%>%
adorn_percentages('row') %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns(position = 'front') %>%
pivot_longer(cols = -area, names_to = 'is_like', values_to = 'n') %>%
separate_wider_delim(n, ' (', names = c('n','percent')) %>%
mutate(
n = parse_number(n),
percent = parse_number(percent)/100
) %>%
ggplot(., aes(x = area, y = n, fill = is_like)) +
geom_bar(position = 'fill', stat = 'identity') +
geom_text(
aes(label = paste0(scales::percent(percent), ' (', n, ')')),
position = position_fill(vjust = 0.5),
size = 5) +
coord_flip() +
scale_fill_brewer(palette = 'Set1') +
theme_bw() +
scale_y_continuous(labels = scales::percent) +
labs(title = '横%ベースのクロス集計', x = NULL, y = NULL) +
theme(text = element_text(size = 15, family = "Fira Sans"))
クロス集計表の見方
言わずもがなではあるが、クロス集計表(contingency table)は2つの変数の関係を見ている表である。contingencyはそもそも「偶然性、不測」などを意味する言葉であり、クロス集計表の割合をもって2つの変数間の偶然性を示している表ともいえる。一方で、変数間に全く関係がない場合、統計的独立と呼ばれる。
右端の列を行周辺合計(row marginal total)
、下端の行を列周辺合計(column marginal
total)と呼ぶ。右下端の角にあたる部分を総計(grand
total)と呼ぶ。また、クロス表の側面(ここではArea)は表側と呼び、頭を(ここではIsLike)は表頭と呼ぶ。
ct## IsLike
## Area No Yes Total
## osaka 29 (29.0%) 71 (71.0%) 100 (100.0%)
## tokyo 123 (61.5%) 77 (38.5%) 200 (100.0%)
## Total 152 (50.7%) 148 (49.3%) 300 (100.0%)探索的にクロス集計表を作成することもあるが、基本的には仮説(因果的な関係)を想定してクロス集計表を作ることが多い。クロス集計から因果関係を導くのは簡単ではない点は注意が必要。今回のケースでは、出身地によってたこ焼きの好き嫌いが変わるだろう、というような関係性を仮定している。もちろん逆でも良く、たこ焼きが好きな人は地域によって偏りがある、でもよい。
どちらが正しいというのはあるわけではないが、時間的な関係性を考慮するほうが個人的にはじっくりくる。つまり、今回のケースであれば、人が生まれた時に出身地は決まるので、その後にたこ焼きの好き嫌いが形成されるような時間的な順序を想定している。たこ焼きの好き嫌いが出身地を決めるとは考えにくい一方で、たこ焼きが好きな集団の特徴を見たいのであれば、列方向にみることで現象の結果としての集団が持つ性質をみるのもありだと思われる。
出身地によってたこ焼きの好き嫌いが変わるだろう、と仮定するのであれば、割合に関しても行％をみればよい。出身地が独立変数(X)、たこ焼きの好き嫌いが目的変数(Y)を想定しているためである。好き嫌い率は、大阪、東京でそれぞれ何パーセントが好きなのかということなので、行パーセントをみるのが適切だと思われる。
このケースでは「大阪は71%、東京は38%より、大阪の人は東京の人よりもたこ焼きが好き」と考えられる(因果関係ではない)。また、行のTotalには東京、大阪の情報がない状態でのたこ焼きの好き嫌いに関する回答が記録されている。東京、大阪の情報がない状態であれば、49%がたこ焼きが好きではあるが、大阪に限定すると、71%がたこ焼き好きであるため、東京と大阪の出生地の関係性はありそうだとわかる。当たり前ではあるが、たこ焼きの好き嫌いを出身地がすべて説明しているわけではない。
対数線形モデル
クロス集計表はいくつかの問題がある。それは、変数内のカテゴリが多くなったり、多層のクロス集計表となると解釈が難しくなる。この問題を対処する方法として、(多重)コレスポンデンス分析などがあるが、ここでは対数線形モデルをまとめておく。
対数線形モデルは、1つの変数を目的変数として扱うモデルではなく、すべての変数を同じように扱うことで、変数間の関係を明らかにすることを目的にしている。ここではあまりありがたくないが、説明のために2×2表の先ほどのデータを利用する。
クロス集計表では2つの変数(X,Y)が「独立」の場合、Xがどの状態であってもYの割合は同じであり、Yがどの状態であってもXの割合は同じになる。独立が成り立つ場合、各セル\(n_{ij}\)の期待値\(\mu_{ij}\)は、
\[ \mu_{ij} = N p^{X}_{i} p^{Y}_{j} \]
となり、\(n_{ij}\)は個数データであるため、対数をとって下記のモデルを考えることができる。ポアソン回帰分析と同じく、このモデルは定数項、\(i\)に依存する項、\(j\)に依存する項に分けることができる。
\[ log \mu_{ij} = logN + logp^{X}_{i} + logp^{Y}_{j} \]
対数線形モデルでは\(log \mu_{ij}\)に対して、パラメタ\(\lambda\)を使って、\(log \mu_{ij}\)を定数項、変数\(X\)のカテゴリに対応するパラメタ\(\lambda^{X}_{i}\)、変数\(Y\)のカテゴリに対応するパラメタ\(\lambda^{Y}_{j}\)の和として表現する。\(\lambda^{X}_{i}\)や\(\lambda^{Y}_{j}\)は主因子項と呼ばれる。
\[ log \mu_{ij} = \lambda + \lambda^{X}_{i} + \lambda^{Y}_{j} \]
そして、対数線形モデルでは下記の制約のもとで、パラメタを推定する。
\[ \sum^{I}_{i=1} \lambda^{X}_{i} = 0, \sum^{J}_{j=1} \lambda^{Y}_{j} = 0 \]
Rではvcdパッケージのloglin()関数で対数線形モデルを実行できる。引数にはクロス集計表やモデル構造を渡す必要がある。margin = list(c(1), c(2))は独立で交互作用がないモデルを想定している。
library(vcd)
crosstable <- xtabs(~ area + is_like, df)
# MASS::loglm(~ area + is_like, crosstable)でフィットしてcoef()でも同様に係数が計算できる
fit <- loglin(crosstable, margin = list(c(1), c(2)), param = TRUE, fit = TRUE)## 2 iterations: deviation 2.842171e-14areaはtokyoが多く、is_likeはnoが多いことがわかる。各変数のカテゴリが、数値の大小で判断できるので、表が巨大で複雑になると便利なのがわかる。
fit$para## $`(Intercept)`
## [1] 4.258508
##
## $area
## osaka tokyo
## -0.3465736 0.3465736
##
## $is_like
## No Yes
## 0.01333412 -0.01333412期待度数、尤度比統計量、χ二乗値を取り出すことも出来る。
list(
ExpectedValue = fit$fit,
LikelihoodRatioTestStatistics = fit$lrt,
# chisq.test(xtabs(~ area + is_like, df), correct = F)と同じ
ChiSquare = fit$pearson,
DegreeOfFreedom = fit$df
)## $ExpectedValue
## is_like
## area No Yes
## osaka 50.66667 49.33333
## tokyo 101.33333 98.66667
##
## $LikelihoodRatioTestStatistics
## [1] 28.82108
##
## $ChiSquare
## [1] 28.17167
##
## $DegreeOfFreedom
## [1] 1変数間の交互作用を想定することできる。むしろクロス集計表では交互作用を想定するほうが自然だと思われる。margin = list(c(1, 2))は交互作用があるモデルを想定している。交互作用項の係数を見ると、osaka-Yes、tokyo-Noが多いことがわかる。つまり、大阪の方ははたこ焼きが好きで、東京の人はたこ焼きが好きではないという関係がわかる。
fit2 <- loglin(crosstable, margin = list(c(1, 2)), param = TRUE, fit = TRUE)## 2 iterations: deviation 0fit2$para$area.is_like## is_like
## area No Yes
## osaka -0.3409407 0.3409407
## tokyo 0.3409407 -0.3409407モンテカルロシュミレーション
最後に。クロス集計表の結果を用いて、モンテカルロシュミレーションを行ってみる。多項分布はこんな感じの分布。
\[ \begin{eqnarray*} f(x_{1}, x_{2}, …, x_{k}) &=& \displaystyle \frac{n!}{x_{1}! x_{2}! … x_{k}!} p_{1}^{x_{1}} p_{2}^{x_{2}} … p_{k}^{x_{k}} ~~ (x_{i} \geq 0, ~~ x_{1} + … + x_{k} = n) \end{eqnarray*} \\ p_{i}>0 ~~ (i = 1, 2, …, k), ~~ p_{1} + p_{2} + … + p_{k} = 1 \]
天下り的にクロス集計表を用意しておく。何らかの変数間の関係が下記のクロス集計として得られたとする。
# クロス集計表の値
ct <- matrix(
c(11, 25,
35, 31),
nrow = 2,
byrow = TRUE
)
ct## [,1] [,2]
## [1,] 11 25
## [2,] 35 31χ二乗検定を行ったところ、5%で有意な結果が得られた。ただ、この結果はたまたまかもしれないので、モンテカルロシュミレーションを行ってみる。
chisq.test(ct)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: ct
## X-squared = 3.8878, df = 1, p-value = 0.04864まずはクロス集計表から、多項分布のパラメタを計算する。パラメタを計算するといっても、各セル度数を総度数で割ることで各セルの確率を算出する。
# 多項分布のパラメータを推定
N <- sum(ct)
prob <- ct / N
list(
N = N,
Prob = prob,
Params = c(prob[1,], prob[2,])
)## $N
## [1] 102
##
## $Prob
## [,1] [,2]
## [1,] 0.1078431 0.2450980
## [2,] 0.3431373 0.3039216
##
## $Params
## [1] 0.1078431 0.2450980 0.3431373 0.3039216あとはこのパラメタをもとに多項分布から乱数を生成する。シュミレーション回数は1万回とする。rmultinom()関数は列で1回の乱数生成の結果を表す。
# サンプリング
n_samples <- 10000
x <- rmultinom(n = n_samples,
size = N,
prob = c(prob[1,], prob[2,])
)
# 10回分の結果
x[,1:10]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 12 11 10 10 12 9 13 7 16 9
## [2,] 28 21 25 21 27 28 21 28 25 18
## [3,] 26 34 42 33 35 30 32 37 39 45
## [4,] 36 36 25 38 28 35 36 30 22 30生成されたシュミレーション結果のデータを用いて、1万回、χ二乗検定を行い、有意かどうか判定する。
res <- vector(mode = 'logical', length = n_samples)
for (i in seq_len(n_samples)) {
m <- matrix(
c(x[1,i], x[2,i],
x[3,i], x[4,i]),
nrow = 2,
byrow = TRUE
)
res[i] <- chisq.test(m, correct = FALSE)$p.value
}
hist(res,
main = "モンテカルロシミュレーション結果",
xlab = "p value",
breaks = 100)
シミュレーション結果のp値の平均を計算する。今回の結果であれば、クロス集計表を得た直後のχ二乗検定ではたまたま有意となったが、シュミレーションして繰り返した結果を見ると、たまたまそうなった可能性が高い。
cat("モンテカルロシミュレーションによるp値の平均:", mean(res), "\n")## モンテカルロシミュレーションによるp値の平均: 0.1110431