児童虐待相談件数の二次分析
Poisson-Gammaモデルの階層ベイズ


PoissonGammaBar.png
author: Mr.Unadon (見習い飯炊き兵)
動作環境:Mac OS Sierra 10.12.1; R version3.3.1; rstan 2.10.1



はじめに


児童虐待について考える第一歩。

私にできることと言えば、データをさわって、事実を掴んで、さらに見えない現場を想像すること。

まずはそこからはじめていきたいと思います。

悲しいものは、できれば見たくないものです。

しかし、データを通じて知識を得て、

そこから現場の一瞬一瞬へと想像力を働かせれば、

そこに、あるいはその先に、あたたかい世界が見えてくるかもしれません。

今回は、Poisson-Gammaモデルを用いた階層ベイズで、

  1. 都道府県別の児童虐待相談件数の平均
  2. 都道府県別の人口に占める児童虐待相談件数の比率

この二つを別々のモデルで推定していきます。


使用データ


二つのデータを使用します。 1つは、児童虐待相談件数の総数データ。 平成26年度福祉行政報告例にある、都道府県別の児童虐待相談件数の総数データです。 児童相談所、福祉事務所、保健センター、市町村の福祉事務所や保健センター、保育所、児童福祉施設、 指定医療機関、保健所、医療機関、幼稚園、学校、教育委員会等に寄せられた案件の総数です。

もう1つは、都道府県別人口のデータ。 政府統計で入手可能な「第4表 都道府県,男女別人口及び人口性比-総人口,日本人人口(平成26年10月1日現在)」です。

データを読み込んでざっと眺めてみましょう。 reportが都道府県別の児童虐待相談件数の総数で、 populationが県別の人口です。

#############################################################
# Poisson-Gamma Model with maltreatment report data
#############################################################
#Rのワーキングスペースを一度クリアにする
invisible({rm(list=ls());gc();gc()})

#パッケージの読み込み
library(tidyverse) #データ成形
library(formattable) #表での出力パッケージ
library(rstan) #ベイズ推定パッケージ
library(reshape2) #少し古いデータ成形パッケージ

#Stanのチェーンを複数コアで並列化
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())


#データの読み込み
Dat<-read.csv("AbusePopulation.csv",fileEncoding = "shift-jis")%>%
    dplyr::mutate(population=人口*1000,ID=1:47)%>%
    dplyr::rename(prefecture=都道府県)%>%
    dplyr::rename(report=虐待相談件数)%>%
    dplyr::select(prefecture,
                  population,
                  report,
                  ID)

#
#
#  ローデータの確認
#
#

#読み込んだデータの確認
formattable(Dat)

今回使用するデータ「都道府県別人口と虐待相談件数」

データの可視化と特徴の把握


続いて、上記データをグラフにして把握していきます。

2種類出します。県別相談件数の棒グラフと、相談件数のヒストグラムです。

#虐待相談件数県別でグラフ化してデータを把握する
g1<-ggplot(Dat)+theme_set(theme_classic(base_family = "HiraKakuPro-W3"))+
    geom_bar(aes(x=reorder(prefecture,ID), y=report),stat = "identity", position = "dodge") + 
    theme(axis.title = element_text(family = "serif", 
    size = 16), axis.text = element_text(family = "serif", 
    size = 14), axis.text.x = element_text(angle = 90, vjust = 0.5), 
    plot.title = element_text(family = "serif", 
    size = 18,hjust = 0.5)) +labs(title = "都道府県別虐待相談件数", 
    y = "虐待相談件数",x="都道府県")+
    scale_y_continuous(expand = c(0,0))
#グラフを表示
g1

こちらが一つ目。

PoissonGammaBar.png

Figure1 都道府県別の虐待相談件数


東京都と大阪府が飛び抜けています。

これは一体、何人のどんな支援者によって、対応されているのでしょうか。


#相談件数の分布を確認する
g2<-ggplot(Dat)+theme_set(theme_classic(base_family = "HiraKakuPro-W3"))+
    geom_histogram(aes(x=report),colour="black",fill="darkgray",bins=20)+
    scale_y_continuous(expand=c(0,0)) +
    theme(axis.title = element_text(family = "serif", 
    size = 16), axis.text = element_text(family = "serif", 
    size = 14), axis.text.x = element_text(family = "serif"), 
    plot.title = element_text(family = "serif", 
    size = 18, hjust = 0.5)) +labs(title = "全国の虐待相談件数度数分布", 
    y = "観測度数")
    
g2

さらに、二つ目。

PoissonGammaHistogram.png

Figure2 虐待相談件数のヒストグラム

以下、混乱がおきないよう修正致します(2017/05/01 7:56 現在、上記を一部修正致します。)。

平均に比して、分散が大きそうです。 これを”めったに起こらない事象”と捉えるかどうかは難しいところですが、 平均より分散のほうが圧倒的に大きいデータです。 一応確認します。

2017/05/01 7:56 上記を一部修正致します。

同一の平均パラメータにもとづいて全国の相談件数が生成されていると考えるならば、上記のようなヒストグラムによる確認は重要かもしれません。しかし、「県別で異なる平均パラメータをもつポアソン分布から、各県の相談件数が生成されている」と想定して階層ベイズを組む以上、上記ヒストグラムがいかにポアソン分布に似ていようが、各都道府県の相談件数にポアソン分布を設定する必然性はありません。理屈がリンクしていないからダメというわけではなく、「別にポアソン分布でもそうでなくともよい」ということです。

今回は、したがってヒストグラムの形がうんぬんだからではなく、「各県の虐待相談件数がそれぞれポアソン分布に従うと仮定する」という見方をベースに、階層ベイズのモデリングに進んでいきます。今あるデータでポアソン分布より使えそうなものがあれば、是非教えてください。比率の方は、二項ベータ分布が使えそうかな。色々やってみる必要がありそうです。

@berobero11先生、教えていただいて、本当にありがとうございます。


階層ベイズモデリング1【虐待相談件数の推定】


ここから、虐待相談件数の推定を行っていきます。 まずはデータ等の宣言をします。

#
#
#  rstanによるパラメータ推定の実行
#
#

#解析用データ・セットの作成
#Poisson-Gamma Model用
d=Dat$report #相談件数を変数dに
r=length(Dat$ID)#都道府県の数

datastan=list(d=d,r=r) #ベイズ推定パッケージrstanに渡すデータのリスト

#推定するパラメータ
parameters=c("lambda", #相談件数数推定平均 i県分
             "alpha", #ガンマ分布ハイパーパラメータ
             "beta", #ガンマ分布ハイパーパラメータ
             "log_lik"
)

#.stanファイルに書いた階層ベイズモデルコードを読み込み
#PoissonGammaModel1<-stan_model("PoissonGammaModel1.stan")
#読みこんだ階層ベイズモデルコードC++言語にコンパイル
#saveRDS(PoissonGammaModel1,"PoissonGammaModel1.rds")

#上のコードで一回コンパイルすれば、次からはこのコードで読み込む。
PoissonGammaModel1<-readRDS("PoissonGammaModel1.rds")

#rstanで推定
fit<-sampling(PoissonGammaModel1,data=datastan,pars=parameters,iter=4000,warmup=2000,thin=2)


使用するモデル(PoissonGammaModel1)はとても単純です。 都道府県別(添字i)の相談件数(di)がポアソン分布に従うと仮定し、 ポアソン平均のλiにガンマ分布を仮定します。

ガンマ分布のパラメータαとβは、全国で共通。それぞれ指数分布とガンマ分布を事前分布に設定します。

以下が、.stanファイルのモデルコードです。

data{
    int r;
    int d[r];
}
parameters{
    real <lower=0> lambda[r];
    real <lower=0> alpha;
    real <lower=0> beta;
}
model{
    for(i in 1:r){
        d[i] ~ poisson(lambda[i]);
        lambda[i] ~ gamma(alpha,beta);
    }
    alpha ~ exponential(1);
    beta ~ gamma(0.1,1);
}
generated quantities{
    real log_lik[r];
    for (i in 1:r){
        log_lik[i] = poisson_lpmf(d[i]|lambda[i]);
    }
    
}

なお、事前分布に一様分布を仮定して感度分析を行いましたが、 以下と結果は同様でした。

早速、可視化します。


#
#
# 階層ベイズ推定結果の整理
#
#

#MCMC サンプルの取り出しと整理
lambdas<-rstan::extract(fit)$lambda #予測相談報告数
results<-matrix(NA,47,4)#空の入れ物

for(i in 1: 47){
    results[i,1]<-mean(lambdas[,i])#予測平均
    #確率密度最大の値(Maximum a Posteriori)
    results[i,2]<-density(lambdas[,i])$x[which(density(lambdas[,i])$y == max(density(lambdas[,i])$y))]
    results[i,3]<-quantile(lambdas[,i],0,025)#95%信用区間下限
    results[i,4]<-quantile(lambdas[,i],0.975)#95%信用区間上限
}
df.res<-data.frame(results)#データフレーム化
df.res$都道府県<-Dat$prefecture
names(df.res)<-c("Mean","MAP","Lower95CI","Upper95CI","都道府県")#列に名前をつける
df.res<-df.res%>%select(都道府県, Mean, MAP,Lower95CI,Upper95CI)

#csvファイルで結果を保存
write.csv(df.res,"PoissonGammaModelResults.csv",fileEncoding = "shift-jis")

#可視化
gg1<-ggplot(df.res)+theme_bw()+
    geom_hline(yintercept = mean(d),linetype=2)+
    geom_hline(yintercept = 0,linetype=1)+
    
        geom_errorbar(aes(x=reorder(df.res$都道府県,as.numeric(Dat$ID)),ymax=Upper95CI,ymin=Lower95CI))+
    geom_point(aes(x=reorder(df.res$都道府県,as.numeric(Dat$ID)),y=Mean))+
    theme(axis.title = element_text(family = "serif", 
    size = 16), axis.text = element_text(family = "serif", 
    size = 14, vjust = 0.5), axis.text.x = element_text(angle = 90), 
    plot.title = element_text(family = "serif", 
    size = 18,hjust = 0.5)) +labs(title = "都道府県別虐待相談件数の推定平均値", 
    y = "虐待相談件数の推定平均値",x="都道府県")+
    scale_y_continuous()+
    theme(axis.text = element_text(hjust = 1), 
    axis.text.x = element_text(vjust = 0.5))
gg1

Poisson-Gammaモデルを用いた、 都道府県別、児童虐待報告件数の階層ベイズ推定結果です。

PoissonGammaResult.png

Figure3 虐待相談件数の推定値
黒点はEAP推定量、エラーバーは95%信用区間を示す。0を実線、全国平均は点線。

どこが多い、どこが少ない、 というのも一つの重要な論点ですが、

ゼロじゃないことにも、目を向けてみたく思います。


階層ベイズモデリング2【人口に占める虐待相談率の推定】


上記のデータは、何より人口に依存したものと考えられます。 東京も大阪も、人口が多いので相談件数も多くなっているはずです。

ここからは、人口に占める虐待相談率の推定を行っていきます。 なお、相談件数は同一個人のものを1と数えているとは到底思えないので、あくまで件数という数値ではあります。

データの整理からStanの推定まで一気に行きます。 なお、便宜上、このモデルをPoisson-Gamma-Thetaモデルと呼ぶことにします。


#
#
# 人口に占める虐待相談件数の比率推定
#  Poisson-Gammaモデルをベースに
#
#

#解析用データ・セットの作成
#Poisson-GammaTheta Model用
N=Dat$population #各都道府県人口をNに
d=Dat$report #相談総数を変数dに
r=length(Dat$ID)#都道府県の数

datastan=list(d=d,r=r,N=N) #rstanに渡すデータリスト
parameters=c("theta", #相談率 i県分
             "alpha", #ベータ分布ハイパーパラメータ
             "beta", #ベータ分布ハイパーパラメータ
             "log_lik")

#.stanファイルの読み込み
#PoissonGammaThetaModel1<-stan_model("PoissonGammaThetaModel1.stan")
#読み込んだ.stanファイル(階層ベイズモデルコード)をC++言語に変換して解析を高速化
#saveRDS(PoissonGammaThetaModel1,"PoissonGammaThetaModel1.rds")
#一度C++変換したら、次からは以下のコードで読み込むだけ
PoissonGammaThetaModel1<-readRDS("PoissonGammaThetaModel1.rds")

#rstanによる階層ベイズ推定の実行
fit2<-sampling(PoissonGammaThetaModel1,data=datastan,pars=parameters,iter = 4000,warmup = 2000,thin=2)

stanのコードを以下に記載します。 先ほどと違うのは、各都道府県の児童虐待相談件数の推定平均λiの代わりに、 人口(N i)×推定相談比率(theta i)を置いている点です。

人口はデータとして所与しますので、thetaが相談率としての機能を持って推定されます。

data{
    int r;
    int d[r];
    int N[r];
}
parameters{
    real <lower=0,upper=1> theta[r];
    real <lower=0> alpha;
    real <lower=0> beta;
}
model{
    for(i in 1:r){
        d[i] ~ poisson(N[i]*theta[i]);
        theta[i] ~ beta(alpha,beta);
    }
    alpha ~ uniform(0,1000);
    beta ~ uniform(0,1000);
}
generated quantities{
    real log_lik[r];
    for (i in 1:r){
        log_lik[i] = poisson_lpmf(d[i]|N[i]*theta[i]);
    }
    
}

さて、結果の確認に移ります。

なお、感度分析として、 ベータ分布の超パラメータαとβにガンマ分布(gamma(0.01,0.01))を置いてみましたが、 結果も結論も変化しませんでした。

次のコードで、結果を可視化します。

#
#
# rstanによる推定結果整理と可視化(PoissongammaThetaModel)
#
#

#MCMC サンプルの取り出しと整理
thetas<-rstan::extract(fit2)$theta #都道府県別人口に占める相談件数の比率
results<-matrix(NA,47,4)#空の入れ物

for(i in 1: 47){
    results[i,1]<-mean(thetas[,i])#予測平均
    #確率密度最大の値(Maximum a Posteriori)
    results[i,2]<-density(thetas[,i])$x[which(density(thetas[,i])$y == max(density(thetas[,i])$y))]
    results[i,3]<-quantile(thetas[,i],0.025)#95%信用区間下限
    results[i,4]<-quantile(thetas[,i],0.975)#95%信用区間上限
}
df.res<-data.frame(results)#データフレーム化
df.res$都道府県<-Dat$prefecture
names(df.res)<-c("Mean","MAP","Lower95CI","Upper95CI","都道府県")#列に名前をつける
df.res<-df.res%>%select(都道府県, Mean, MAP,Lower95CI,Upper95CI)

#csvファイルで結果を保存
write.csv(df.res,"PoissonGammaThetaModelResults.csv",fileEncoding = "shift-jis")

gg2<-ggplot(df.res)+theme_bw()+
    geom_hline(yintercept = c(mean(d)/mean(N)),linetype=2)+
    geom_hline(yintercept = 0,linetype=1)+
    
        geom_errorbar(aes(x=reorder(df.res$都道府県,as.numeric(Dat$ID)),ymax=Upper95CI,ymin=Lower95CI))+
    geom_point(aes(x=reorder(df.res$都道府県,as.numeric(Dat$ID)),y=Mean))+
    theme(axis.title = element_text(family = "serif", 
    size = 16), axis.text = element_text(family = "serif", 
    size = 14, vjust = 0), axis.text.x = element_text(angle = 90), 
    plot.title = element_text(family = "serif", 
    size = 18,hjust = 0.5)) +labs(title = "都道府県別人口に占める虐待相談率の推定値", 
    y = "児童虐待相談率の推定値",x="都道府県")+
    scale_y_continuous()+
    theme(axis.text = element_text(hjust = 1), 
    axis.text.x = element_text(vjust = 0.5))
gg2

Poisson-Gamma-Thetaモデルを用いた、 都道府県別の人口に占める児童虐待報告件数率の推定結果です。

PoissonGammaThetaResult.png

Figure4 人口に占める虐待相談率の推定値
黒点はEAP推定量、エラーバーは95%信用区間を示す。0を実線、全国平均は点線。

最後に


Poisson-Gammaモデルは、たった1~2点のデータ(e.g. 虐待相談件数、人口)が47都道府県分あるだけで、 感度分析にも耐える推定結果を出してくれます。 このモデルは色々な場面に応用できます。相対リスクの推定にも使用できる。 アイデア次第で二次分析の幅がぐっと広がります。

データを触ってみて、 もっと知りたいと思いました。

あったかいご飯と、あったかい布団と、あったかい笑顔のある環境で、 生きていってほしいな。

わたしも、自分の身の回りから、取り組んで行きます。


Written on May 1, 2017