児童虐待対応の仕事、いわゆるケースワーカー（CW）の仕事は日々とても辛いものがあります。
でもどうして辛いんでしょうか。
仕事の内容の激しさももちろんありますが、その他のとてもシンプルな理由に、仕事の多さ、次々と終わらない担当の受理などがあります。
要するに、この前の仕事が終わっていない(すぐ終わるものでもないし)のに次々くる感じです。
これはもう誰でも嫌になっちゃいますよね。
そんなシンプルなつらさってどうやってモデルで表現するのでしょうか。
今回、色々と小難しいことは抜きにした「CW嫌々モデル」を考えてみました。
架空データとして、「嫌だなって気持ちを点数化したもの（高いほど嫌）」である嫌得点（iyapoint）と、仕事が増えた・殺害予告された・個人名で訴えると怒鳴られた・殴られた（どれも日々のあるあるだね！！）のイベントが発生したかを示す2値データ(ivent)が、日々の時系列データでどのように変化するかを見ていきます。
day iyapoint ivent
1 2019/4/1 7 0
2 2019/4/2 7 0
3 2019/4/3 9 1
4 2019/4/4 8 0
5 2019/4/5 9 1
6 2019/4/6 9 0
7 2019/4/7 8 0
↑こんな感じのデータです。
早速モデルを考えていきます。
◆モデル１：トレンド＋嫌イベントストレスモデル
１．放っておくと嫌々度は下がる；-Trend（負の傾きの直線をイメージしてください）２．定期的に負のイベントがある：CoefIvent*Iventt (トレンドの係数＝一定)×嫌なイベントの有無
１，２よりイベントの影響を含めたトレンド：δt＝Trend＋CoefIvent*Iventt
３．μt＝μt-1＋δt-1+wt, wt~(0,σ2w)：μtはｔ時点の嫌々度の状態
４．μnoiset=μt+Vt：Yt,Vtは各々ｔ時点の過程誤差と観測誤差
５．λt=exp(μt)
６．Yt~Poisson(λt)
嫌なイベントの発生はまれ(だといいなという想定)という前提で、ポアソン分布を使用。でも嫌々は得点化したから、ポアソンのような離散分布でなく連続データの分布にした方がよかったかなぁと今更ながら思った。
少し雑念が入りましたが、上記1～6を元に最初の数理モデルを組み、統計処理をしてみようと思います。
統計ソフト：R-3.5.2
使用パッケージ：Rstan他
連続データを用いた、状態空間モデルによるベイズ推定を行いました。
以下Stanコードを記します。
stancode <- "
data {
int n; // サンプルサイズ
int ivent[n]; // 嫌イベントの有無
int iyapoint[n]; // 嫌得点
}
parameters {
real trend; // 嫌得点減少トレンド
real coef_ivent; // 嫌イベント発生がトレンドに与える影響
real mu_zero; // 状態の初期値
real mu[n]; // 状態の推定値
real mu_noise[n]; // 観測誤差の入った状態の推定値
real<lower=0> s_w; // 過程誤差の分散
real<lower=0> s_v; // 観測誤差の分散
}
transformed parameters{
real delta[n]; // イベントの影響が入ったトレンド
real lambda[n]; // ポアソン分布の期待値
for(i in 1:n){
delta[i] = -trend + coef_ivent * ivent[i];
}
for(i in 1:n){
lambda[i] = exp(mu_noise[i]);
}
}
model {
// 状態の初期値から最初の時点の状態が得られる
mu[1] ~ normal(mu_zero, sqrt(s_w));
// 状態方程式に従い、状態が遷移する
for (i in 2:n){
mu[i] ~ normal(mu[i - 1] + delta[i - 1], sqrt(s_w));
}
// 観測誤差が加わる
for(i in 1:n){
mu_noise[i] ~ normal(mu[i], sqrt(s_v));
}
// ポアソン分布に従って観測値が得られる
for(i in 1:n){
iyapoint[i] ~ poisson(lambda[i]);
}
}
generated quantities{ //新たにサンプリングする変数を作る
vector[n] log_lik;
for(i in 1:n){
//モデル1の対数尤度
log_lik[i] = poisson_lpmf(iyapoint[i] | lambda[i]);//ポアソンの場合、lpdfではなくてlpmfであることに注意。https://www.slideshare.net/simizu706/stan-64926504の42
}
}
"
fit <- stan(model_code = stancode, iter = 8000, chains = 1, data = standata, seed=123)
◆モデル２：CW嫌イベント確率的ボラティリティ変動モデル
確率的ボラティリティモデルとは、金融界でよく使われるモデルです。ボラティリティ（分散的な）が常に一定であるランダムウォークとは異なり、ボラティリティ自体もランダムウォークするというモデルです。株価みたいに、上がる時はじわじわ上昇トレンドで上がるくせに、下がる時は理不尽にガタッとくる値動きをよく表現できているモデルだと思います。
何もなければ（ありえないけど）すこーしずつストレスが和らぐけど、ちょっとしたことで一気に限界突破する児童虐待の仕事と似てると思いませんか？
以下、モデル２のStanコード
stancode2 <- "
data {
int n;
real iyapoint[n];
}
parameters {
real h[n];
real mu;
real<lower=-0.999,upper=0.999> phi;//ボラティリティの持続性
real<lower=0.0001> sigma;
}
model {
h[1] ~ normal(mu,sigma/sqrt(1-phi*phi));
iyapoint[1] ~ normal(0,exp(h[1]/2));
for (i in 2:n){
h[i] ~ normal*1;
}
}
generated quantities{ //新たにサンプリングする変数を作る
vector[n] log_lik;
for(i in 1:n){
//モデル2の対数尤度
log_lik[i] = normal_lpdf(iyapoint[i] | 0,exp(h[i]/2));
}
}
"
fit2 <- stan(model_code = stancode2, iter = 8000, chains = 1, data = standata, seed=1234)
そして、このモデル１とモデル２のどちらがデータに適合したモデルなのでしょうか。
適したモデルこそ、児童虐待CWの辛さをよく表現したモデルってことになります。
モデルの比較はWAICを使用します。
適合度比較では、構造方程式モデリングではAICとかGFIとかRMSEAとかがスタンダードですが、ベイズではWAICがよく使用されています。以下にRで実行する際のコードを示します。
log_lik1 <- extract_log_lik(fit)#waic用コードで計算したfit12の対数尤度の値を取り出す
waic1 <- waic(log_lik1)#変数名を付ける
print(waic1 , digits = 3)#waic1の結果表示
log_lik2 <- extract_log_lik(fit2)
waic2 <- waic(log_lik2)
print(waic2 , digits = 3)
compare(waic1,waic2)#waicの比較.昇順に並べ替える
print(waic1 , digits = 3)#waic1の結果表示
Estimate SE
elpd_waic -198.555 1.014
p_waic 1.131 0.087
waic 397.111 2.028
log_lik2 <- extract_log_lik(fit2)
waic2 <- waic(log_lik2)
print(waic2 , digits = 3) #waic2の結果表示
Estimate SE
elpd_waic -337.736 1.597
p_waic 0.161 0.013
waic 675.471 3.195
compare(waic1,waic2)#waicの比較
elpd_diff se
-139.2 0.9
★結論：モデル１(トレンド＋嫌イベントストレスモデル)のWAICの方が小さいため、予測の観点から良いモデルといえることが分かりました。
とはいえ、モデル２の確率的ボラティリティな感じも、日常業務である「落ち着いてきていけると思ったけど炎上してやっぱ無理」な感じをよく表現できていると思います。
確率的ボラティリティを用いた、より良いモデルに改良していけるのではないかと思います。今回のはシンプルでしたしね。
