打ち切りのあるデータ
野外のデータ、あるいは実験データを扱う際、「ある特定の値以上であることはわかるけれども、本当の値はわからない」という状況に出くわすことがあります。このようなデータタイプは「打ち切り」と呼ばれ、統計モデルを組むうえで厄介な問題です。該当する例はとして、以下のようなものが挙げられます。
- ある一定量のエサを与え、一定時間後に無くなった量を調べる（全部消費された場合、もともと与えたエサ量「以上」を消費したかったことは確か）
- 標識再捕獲を行い、ある時点で生きていたかどうかを調べる（最後の再捕獲で生きていた場合、その時点より”長く”生きていたことは確か）
他にも例は挙げればたくさんありますが、こういう性質のデータを扱うのは頭を抱えます。
センサリング（Censoring）とは
センサリングとは、そのまま「打ち切り」を意味する言葉なのですが、このテクニックにより上記の問題に対応することが可能です。センサリングを適用する場合、上限値に達した観測値は一時的に欠測、NAとして扱われます。BUGSソフトウェアでは、（通常の線形モデルの場合）応答変数の欠測は自動的に推定されるべきパラメーター（つまり、ほかの偏回帰係数などの推定には用いられない情報）として扱われ、ランダムな値があてがわれます。しかし、打ち切りの場合は部分的な情報（y >= 60など）を持っていますので、上限値に達した値に対しては、こちらで定義する上限値以上の値が（ランダムに？）割り振られます。そうした値が割り振られたうえで、線形モデルの各種パラメーターを推定することになります。
JAGSによる実装
基本的な実装の手順は、普通のBUGSモデルのコードに1行追加するだけで済みます。こちらのモデルファイルを参考にBUGSコードを書いてみます。
model{
#Likelihood
for(i in 1:N.sample){ # loop for i
isCensored[i] ~ dinterval(y[i],C[i])
y[i] ~ dpois(lambda[i]) # poisson model
log(lambda[i]) <- alpha + beta * X[i] # linear model structure
}
#Prior
alpha ~ dnorm(0, 0.001)
beta ~ dnorm(0, 0.001)
}
dintervalが打ち切り関数になっており、isCensored[i]（こちらでRから与えるデータ）はデータiが打ち切りかどうかを示す変数（0:打ち切り無し、1:打ち切りあり）、C[i]は打ち切った場合の上限値（こちらデータを与える；任意の変数）になります。なお、打ち切りありの場合、yは”NA”として与えることになります。例えば、観測の上限がy=7だった場合、データは下記のようになります。
> head(D) y isCensored C x 1 NA 1 7 12.056546 2 NA 1 7 11.205303 3 4 0 7 9.865988 4 NA 1 7 12.633539 5 NA 1 7 9.116709 6 6 0 7 12.067588
図にするとこんなふうになります。
センサリングモデルと普通のポアソンモデルで推定値を比較してみます。
#true values; alpha=0.3, beta=0.15
#with censoring model
Lower95 Median Upper95 Mean SD Mode MCerr MC%ofSD SSeff AC.200 psrf
alpha 0.135 0.59 1.08 0.59 0.241 NA 0.0167 6.9 208 0.100 1
beta 0.083 0.13 0.18 0.13 0.024 NA 0.0017 7.0 202 0.097 1
#with no censoring model
Lower95 Median Upper95 Mean SD Mode MCerr MC%ofSD SSeff AC.200 psrf
alpha 0.587 0.970 1.41 0.972 0.212 NA 0.010 5 407 -0.0027 1
beta 0.026 0.071 0.11 0.071 0.021 NA 0.001 5 403 -0.0037 1
センサリング関数を使わない場合、y>=7の値はy=7とされてしまうので、推定値がかなり過小評価されています。それに対し、センサリング関数を用いたモデルでは、推定値の偏りが解消され、その95%CIも真の値を含んでいます。
注意点：センサリング関数を用いる場合、打ち切りされたデータについては、「打ち切り上限値以上の値」を初期値として与える必要があります。上の例で挙げると、y=NAとなっているデータについては、y=8などの値を初期値として与えます。そうしないとエラーメッセージが帰ってきます。
SPONSOR LINK