以前、segmentedパッケージを使って変化点を探す方法を調べました。
necostat.hatenablog.jp
segment( )はglmオブジェクト（=glm関数で当てはめたモデル）とlmeオブジェクト（=変量効果を含む線形回帰モデル）には対応しているんですが、lme( )は変量効果を含んだ一般化線形モデルの当てはめが出来ません。馴染みのあるglmer( )をsegmented( )に渡せると楽なんですが...（実は方法があるんでしょうか？）。
Stanで実装すればマルチレベルモデルでの変化点検出もできるかもと思ってやってみました。
階層構造のないデータ
まずはsegmented( )のまとめのときと同じく、downデータで年齢についての変化点を探してみます。
library(segmented)
library(rstan)
data(down)
使用したStanコード
次のようなコードを使いました（//より後ろはコメントなので削除可）。
data{
int I; //total no. of obs.
int N[I]; //no. of trial
int n[I]; //no. of success
real Z[I]; //explanatory variable
}
parameters{
real<lower=min(Z)> psi; //breakpoint of Z
real b01; //intercept1
real b02; //intercept2
real b1; //baseline slope
real b2; //difference-in-slope
}
model{
for(i in 1:I){
if (Z[i] < psi){
target += binomial_logit_lpmf(n[i]|N[i], b01 + b1*Z[i]);
}else{
target += binomial_logit_lpmf(n[i]|N[i], b02 + (b1+b2)*Z[i]);
}
}
}
各ブロックごとに見ていきます。
まずdataブロックでは、モデルに渡す変数として以下のものを定義しています。
- I：総観察数（データフレームの行数）
- N：試行数（ここでは出生総数
births）
- n：成功数（ここでは発生数
cases）
- Z：説明変数（ここでは母体年齢
age）
parametersブロックでは、変化点（psi）と変化点の前後で描かれる2本の直線の切片と傾きを定義しています。変化点後の直線の傾きを直接与えるのではなく、傾きの変化量（b2）として定義しました（この方がb2=0かどうか判断しやすそうだから）。切片はあまり興味がないのでそのまま定義してます。
modelブロックでは、変化点の前後で異なったモデルを定義しています。つまり、変化点（psi）よりも低い年齢のときは、
を当てはめ、変化点（psi）よりも高い年齢のときは、
を当てはめるようにしています。
bayes_segmented.stan という名前で保存しておきます。
サンプリング実行コード
datalist <- list(
I = nrow(down),
N = down$births,
n = down$cases,
Z = down$age
)
model <- stan_model("bayes_segmented.stan")
fit <- sampling(model, data=datalist,
iter=10000, warmup=5000, thin=10,
chain=4, cores=4, seed=1234)
結果
一部のパラメータのトレースプロットです。
stan_trace(fit, pars=c("psi","b01","b02","b1","b2"))
一部のパラメータのヒストグラムです。
stan_hist(fit, pars=c("psi","b01","b02","b1","b2"))
パラメータ推定値のサマリーです。
> round(summary(fit)$summary,digit=2)
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
psi 30.70 0.07 1.65 28.19 29.32 30.61 31.69 34.74 544.12 1.00
b01 -6.83 0.03 0.49 -7.76 -7.16 -6.83 -6.50 -5.88 284.48 1.02
b02 -15.33 0.03 0.69 -16.94 -15.72 -15.26 -14.86 -14.14 654.03 1.00
b1 -0.01 0.00 0.02 -0.05 -0.03 -0.01 0.00 0.03 280.82 1.02
b2 0.27 0.00 0.02 0.22 0.26 0.27 0.29 0.32 297.88 1.02
lp__ -91.09 0.04 1.55 -94.89 -91.87 -90.75 -89.93 -89.12 1336.07 1.00
変化点の平均は30.7歳でした。segmented( )を使ったときの推定値が 31.081なので、まずまず許容範囲でしょうか。
傾き（b1, b2）の推定値もほとんど同じ（-0.01341, 0.27422 ）なので、上手く出来てそうです。
階層構造のあるデータ
次に階層構造のあるデータでやってみます。
set.seed(1234)
n_cluster <- 10
b_RE <- rnorm(n=n_cluster, mean=0, sd=2)
data_RE <- NULL
for(i in 1:n_cluster){
d <- down %>%
mutate(id = i,
theta = plogis(qlogis(cases/births) + b_RE[i]),
cases = rbinom(n=nrow(down), size=births, prob=theta))
data_RE <- bind_rows(data_RE, d)
}
10個のクラスターそれぞれに、正規分布に従うb_REを変量効果として加えています。平均0なので全体の平均にはあまり影響していないはずです。
> b_RE
[1] 1.5104121 -0.4452811 2.6536869 -2.5044597 -1.6837192 -0.3739502 -1.1919180 1.5708385 2.9564148 0.8196293
出来上がったデータを確認してみます。
data_RE %>%
ggplot(aes(x=age, y=cases/births)) +
geom_line(aes(colour=as.factor(id), group=id)) +
theme_bw()
使用したStanコード
階層構造がないときと比べて、以下のデータ、パラメータを追加しています。
- J：クラスターの数
- ID：クラスターID
- RE：各クラスターの変量効果
- sigma：クラスター変量効果の標準偏差
data{
int I; //total no. of obs.
int J; //no. of cluster
int N[I]; //no. of trial
int n[I]; //no. of success
real Z[I]; //explanatory variable
int ID[I]; //ID of cluster
}
parameters{
real<lower=min(Z)> psi; //breakpoint of Z
real b01; //intercept1
real b02; //intercept2
real b1; //baseline slope
real b2; //difference-in-slope
real RE[J]; //effect of cluster
real<lower=0> sigma; //SD of RE
}
model{
RE ~ normal(0,sigma);
for(i in 1:I){
if (Z[i] < psi){
target += binomial_logit_lpmf(n[i]|N[i], b01 + b1*Z[i] + RE[ID[i]]);
}else{
target += binomial_logit_lpmf(n[i]|N[i], b02 + (b1+b2)*Z[i] + RE[ID[i]]);
}
}
}
bayes_segmented_multilevel.stan という名前で保存しておきます。
サンプリング実行コード
渡すデータが増えた他は、先程とほぼ同じです。
datalist <- list(
I = nrow(data_RE),
J = length(unique(data_RE$id)),
N = data_RE$births,
n = data_RE$cases,
Z = data_RE$age,
ID = data_RE$id
)
model <- stan_model("bayes_segmented_multilevel.stan")
fit <- sampling(model, data=datalist,
iter=10000, warmup=5000, thin=10,
chain=4, cores=4, seed=1234)
結果
一部のパラメータのトレースプロットです。まずまずでしょうか。
stan_trace(fit, pars=c("psi","b01","b02","b1","b2","sigma"))
一部のパラメータのヒストグラムです。psiは30.5〜31.5でだいたい一様ですが、元のデータの年齢が1歳刻みだからでしょう。
stan_hist(fit, pars=c("psi","b01","b02","b1","b2","sigma"))
ちょっと見にくいですが、パラメータ推定値のサマリーを全部載せてしまいます。
> round(summary(fit)$summary,digit=2)
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
psi 31.01 0.01 0.31 30.52 30.75 30.99 31.26 31.49 1372.54 1.00
b01 -6.51 0.08 0.71 -7.91 -6.99 -6.51 -6.05 -5.08 81.66 1.02
b02 -15.02 0.08 0.72 -16.45 -15.51 -15.02 -14.56 -13.58 81.35 1.02
b1 -0.01 0.00 0.00 -0.02 -0.01 -0.01 -0.01 -0.01 1380.29 1.00
b2 0.27 0.00 0.00 0.27 0.27 0.27 0.28 0.28 1582.65 1.00
RE[1] 1.20 0.08 0.71 -0.22 0.75 1.21 1.68 2.61 80.47 1.02
RE[2] -0.70 0.08 0.71 -2.14 -1.14 -0.69 -0.22 0.71 79.85 1.02
RE[3] 2.33 0.08 0.71 0.89 1.87 2.33 2.80 3.71 80.18 1.02
RE[4] -2.73 0.08 0.72 -4.22 -3.18 -2.72 -2.24 -1.32 80.37 1.02
RE[5] -1.94 0.08 0.71 -3.41 -2.39 -1.93 -1.46 -0.54 83.70 1.02
RE[6] -0.77 0.08 0.71 -2.19 -1.22 -0.76 -0.29 0.64 79.81 1.02
RE[7] -1.50 0.08 0.71 -2.94 -1.95 -1.50 -1.02 -0.10 80.41 1.02
RE[8] 1.22 0.08 0.71 -0.21 0.77 1.23 1.70 2.62 80.68 1.02
RE[9] 2.63 0.08 0.71 1.20 2.18 2.64 3.10 4.02 79.79 1.02
RE[10] 0.52 0.08 0.71 -0.91 0.07 0.53 1.00 1.91 80.51 1.02
sigma 2.15 0.02 0.62 1.33 1.73 2.03 2.43 3.72 764.72 1.00
lp__ -1807.01 0.08 2.80 -1813.28 -1808.68 -1806.67 -1805.00 -1802.62 1110.61 1.00
変化点・傾き・REの分散などを比べると、まずまず良い出来でしょうか。
おわりに
- データとパラメータとの大小関係で場合分けする方法で上手くいきました。
- 事前に変化点の個数が分かっていないといけないのが弱点です（単調増加に見えるなら1点を見つけたいという状況が多いように思いますが）
- キャットタワーの爪研ぎ部分の麻紐がとれて、芯材が剥き出しになってしまいました。崩壊する前に買い替えを検討中です。
![Statistical Methods In Medical Research, 4Ed [Paperback] [Jan 01, 2017] Armitage P.](https://m.media-amazon.com/images/I/41kQ8DLXpPL._SL500_.jpg "Statistical Methods In Medical Research, 4Ed [Paperback] [Jan 01, 2017] Armitage P.")
