不等号を含むデータの解析

以下のようなデータがあるとします。

PersonID	Y	Treat
1	38.3	0
1	<20	0
1	52.3	0
1	40.8	0
1	<20	0
2	26.4	0
2	<30	0
2	<30	0
2	30.1	0
2	25	0
2	28.5	0
3	45.2	1
3	<35	1
3	57.3	1
3	59.9	1
4	67.2	1
4	48.7	1
4	60.8	1
4	71	1
5	104.2	1
5	85	1
5	91.7	1
5	61.3	1
5	72.1	1

1列目はある人のID、2列目はその人から測定された値（ここでは便宜的にある血中タンパク質の量とします）、3列目は人ごとの処理したか否かの値（処理が1）です。一人当たり何回か測定している場合があるので同じPersonIDが複数行に登場しています。（後述の通りダメなんですが）とりあえず不等号のデータは不等号を取り除いた値で置き換えて図を描きました。色の違いは処理の有無です。

知りたいことは処理の結果、血中のタンパク質の量が上がったどうかです。ただし、測定システムに検出限界がありまして、ある程度低い値は不等号付きで表されてしまっています。こういう状況、他にもありそうですよね。さて、どのように検定しますか？

不等号付きのデータを削る？それは低い値を無視しているのであきらかにダメです。単に不等号を取るとかして、不等号がついているデータを置き換える？それもダメです。考えられる1つの弊害として値がそろってしまって有意差が出やすくなります。困りますよね。僕が知る範囲の基本的な統計学のテキストには対策はどこにも見つからなかったです。こんな状況も自然に扱えます、そうBUGSならね！

背後のメカニズムとしては、PersonIDごとに何か真の値を持っていて、それに測定誤差などのノイズが混ざって測定値が出ると考えます。その測定値がある足きり以下になると不等号がついてしまうという流れです。そこでまずデータを分割します。

PersonID	Treat
1	0
2	0
3	1
4	1
5	1

PersonID	Y	Y.cen
1	38.3	NA
1	52.3	NA
1	40.8	NA
2	26.4	NA
2	30.1	NA
2	25	NA
2	28.5	NA
3	45.2	NA
3	57.3	NA
3	59.9	NA
4	67.2	NA
4	48.7	NA
4	60.8	NA
4	71	NA
5	104.2	NA
5	85	NA
5	91.7	NA
5	61.3	NA
5	72.1	NA
1	NA	20
1	NA	20
2	NA	30
2	NA	30
3	NA	35

上が個人ごとの情報（person.txt）、下が測定値（data.txt）です。data.txtで測定値と足きり値（Y.cen）を導入しました。そして、BUGSコードがラクになるように足きりの値は後ろの方に回しました。

BUGSコードは以下のようになりました。まずは個人差を考慮しないバージョン（model1.bugs）です。

model {
   for (n in 1:N.person) {
      mu[n] <- alpha + beta.tr * Treat[n]
   }
   for (i in 1:N.data.notNA) {
      Y[i] ~ dnorm(mu[PersonID.ix[i]], tau)
   }
   for (i in (N.data.notNA+1):N.data) {
      Y[i] ~ dnorm(mu[PersonID.ix[i]], tau) I( ,Y.cen[i])
   }
   alpha ~ dnorm(0, 1.0E-4)
   beta.tr ~ dnorm(0, 1.0E-4)
   tau <- 1/(s*s)
   s ~ dunif(0, 1.0E+4)

   for (i in 1:N.data) {
      y.pr[i] ~ dnorm(mu[PersonID.ix[i]], tau)
      err.pr[i] <- pow(y.pr[i] - Y[i], 2)
   }
   MSPE <- mean(err.pr[])
}

処理の影響はbeta.trという形で入れています。ポイントは8-10行目です。不等号データのところは切断正規分布というので表します。それをI(lower, upper)で表します。収束が悪くなったりするので多用は禁物ですが、今回のような足きり以下の値を表すにはピッタリです。

16-20行目はモデルの予測の程度を知るためのコードです。muから擬似的に予測値（y.pr; Yのpredict）を発生させ、実測値との二乗誤差を求めて、データあたりの平均値（MSPE; mean squared predictive error）を求めています。

次に個人差を考慮したバージョン（model2.bugs）です。

model {
   for (n in 1:N.person) {
      mu[n] <- alpha + beta.tr * Treat[n] + r[n]
      r[n] ~ dnorm(0, tau[1])
   }
   for (i in 1:N.data.notNA) {
      Y[i] ~ dnorm(mu[PersonID.ix[i]], tau[2])
   }
   for (i in (N.data.notNA+1):N.data) {
      Y[i] ~ dnorm(mu[PersonID.ix[i]], tau[2]) I( ,Y.cen[i])
   }
   alpha ~ dnorm(0, 1.0E-4)
   beta.tr ~ dnorm(0, 1.0E-4)
   for (a in 1:N.tau) {
      tau[a] <- 1/(s[a]*s[a])
      s[a] ~ dunif(0, 1.0E+4)
   }

   for (i in 1:N.data) {
      y.pr[i] ~ dnorm(mu[PersonID.ix[i]], tau[2])
      err.pr[i] <- pow(y.pr[i] - Y[i], 2)
   }
   MSPE <- mean(err.pr[])
}

3-4行目に個人差の項rが入っただけです。

MSPEの結果は以下になります。

	MSPE
model1	0.3
model2	0.2

上記のmodel2.bugsを実行するRコードは以下になります。

source("R2WBwrapper.R")

d <- read.delim("input/data.txt", as.is=TRUE, sep="\t")
person <- read.delim("input/person.txt", as.is=TRUE, sep="\t")
d.notNA <- subset(d, !is.na(Y))
N.tau <- 2

clear.data.param()
set.data("N.data", nrow(d))
set.data("N.data.notNA", nrow(d.notNA))
set.data("N.person", nrow(person))
set.data("N.tau", N.tau)
set.data("Y", d$Y/50)
set.data("Y.cen", d$Y.cen/50)
set.data("PersonID.ix", d$PersonID)
set.data("Treat", person$Treat)

set.param("alpha", 0)
set.param("beta.tr", 0)
set.param("s", rep(1, N.tau))
set.param("MSPE", NA)

post.bugs <- call.bugs(
   # debug = FALSE,
   file = "model/model2.bugs",
   n.iter = 12000, n.burnin = 2000, n.thin = 10
)
post.list <- to.list(post.bugs)
post.mcmc <- to.mcmc(post.bugs)

13-14行目は数値のオーダーを1ぐらいにしておくと収束の速度や精度が良くなるので50で割っています。MCMCサンプルからあとで50を掛けて戻せばOKです。21行目でset.param("MSPE", NA)としておくと、「MSPE」変数のMCMCサンプルも入手できます。MSPE=0.2はおおよそ1測定値あたりsqrt(0.2)*50で±20ぐらいずれているということになります。今回は測定誤差も個人差もかなり大きいのでこれぐらいも納得です。

さて知りたいのはbeta.trのMCMCサンプルの分布および、処理の影響があるかどうか（処理の影響がプラスである確率）です。これらは以下のように求めることができます。

pdf("dist.pdf")
plot(post.list)
dev.off()
beta.tr <- post.mcmc[ ,"beta.tr"]
prob <- sum(beta.tr > 0)/nrow(post.mcmc)

結果は、「確率0.9223で処理の影響あり」になりました。