読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

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

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

PersonIDYTreat
138.30
1<200
152.30
140.80
1<200
226.40
2<300
2<300
230.10
2250
228.50
345.21
3<351
357.31
359.91
467.21
448.71
460.81
4711
5104.21
5851
591.71
561.31
572.11

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

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

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

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

PersonIDTreat
10
20
31
41
51


PersonIDYY.cen
138.3NA
152.3NA
140.8NA
226.4NA
230.1NA
225NA
228.5NA
345.2NA
357.3NA
359.9NA
467.2NA
448.7NA
460.8NA
471NA
5104.2NA
585NA
591.7NA
561.3NA
572.1NA
1NA20
1NA20
2NA30
2NA30
3NA35

上が個人ごとの情報(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
model10.3
model20.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.trMCMCサンプルの分布および、処理の影響があるかどうか(処理の影響がプラスである確率)です。これらは以下のように求めることができます。

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

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