ベイズ法・ポワソンモデルを用いたクロス集計表の分析

クロス表データのベイズ分析のために使用できるいくつかの可能性があるモデルがある。

クロス表の各セルでの頻度に対するいわゆる対数線形モデル（ポアソン・モデル）。各セルのデータの比率の推定値を得るための二項モデルである。カイ二乗検定よりもさらに詳細にデータのパターンを探索する可能性がもたらすだろう。

データをデータフレーム形式のクロス表に並べる。つまり、クロステーブルの各セルの頻度が変数となり、行名と列名も別々の変数で与えられる。

タイタニックデータの呼び出し

生存・性別を使用する。

library(titanic)
data(titanic_train)
dat <- subset(titanic_train, select = c(Survived, Sex))

リコード

数値で入っている生存・死亡を文字にリコード。

library(memisc)
dat$Survived <- memisc::recode(dat$Survived, "Survived" <- 1, "Death"<- 0)
dat$Survived <- as.factor(dat$Survived)
dat$Sex <- as.factor(dat$Sex)

データタイプを確かめる

str(dat)

'data.frame': 891 obs. of 2 variables: $ Survived: Factor w/ 2 levels "Survived","Death": 2 2 2 2 2 2 2 2 2 2 ... $ Sex : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...

スタック形式に変換

dat.stack0 <- table(dat$Sex, dat$Survived)
dat.stack <-as.data.frame(dat.stack0)
colnames(dat.stack) <- c("Sex","Survived","freq")
print(dat.stack)

スタック形式のクロス表。

     Sex Survived freq
1 female Survived  233
2   male Survived  109
3 female    Death   81
4   male    Death  468

対数線形モデル log-linear model

library(arm)
nsim <- 5000

mod <- glm(freq ~ Sex + Survived + Sex:Survived, 
           data=dat.stack, family=poisson)
bsim <- sim(mod, n.sim=nsim)
round(t(apply(bsim@coef, 2, quantile, 
              prob=c(0.025, 0.5, 0.975))),2)

出力

                       2.5%   50% 97.5%
(Intercept)            5.32  5.45  5.58
Sexmale               -0.99 -0.76 -0.54
SurvivedDeath         -1.31 -1.06 -0.81
Sexmale:SurvivedDeath  2.19  2.51  2.84

相互作用パラメータは、相関の強さを測定するものである。パラメータ値2.51の意味を定量的に理解するためには、すべてのパラメータ値の解釈を見る必要がある。

切片5.45は、セル「female」と「Survived」（女性で生存した数）の対数に相当する。

$\exp(β_0)=233$

第2パラメータの指数は，セル「femaleとSurvived」と「maleとSurvived」乗法的な差に相当する。

$exp(\beta_0)exp(\beta_1)= exp(5.45)exp(-0.76) = 109$

となる。第3のパラメータは、「femaleとSurvived」と「femaleとDeath」のセル間の乗法的な差、すなわち「femaleとDeath」のセル内の頻度である。

$exp(\beta_0)exp(\beta_2) = exp(1.79)exp(-1.06) = 81$

第3パラメータは、女性で生存した者と女性で死亡した者の対数値の差である。交互作用パラメータは、この差を、男性と女性の間で差分化したものである。これは直感的に理解するのが難しい。ここで、もう一回、定式化を試みる。交互作用パラメータは、クロス表の行と列の差のカウントの対数の差を測定している。モデルパラメータから「maleとDeath」のセル内の頻度を取り出すと、明らかになるかもしれない。

$exp(\beta_0)exp(\beta_1)exp(\beta_2)exp(\beta_3) = exp(5.45)exp(-0.76)exp(-1.06)exp(2.51) = 468$

参照したのはこちらの解説。

https://tobiasroth.github.io/BDAEcology/addbasics.html