時系列データにt 検定を行う

MikuHatsune2016-05-30

読んだ Nature (2016) doi:10.1038/nature18294
PD1/PDL1 系の薬剤は高いので、それが良く効く患者集団を選びたいというのが最近のトレンドォ!!!
3' 側のUTR がないことで、PD-L1 の発現があがり、それによって免疫機構から逃れるPD1/PD-L1 axis がクッソすごいことになっているが、ニボルマブのようなPD1/PD-L1 を阻害するような薬剤で腫瘍増大が抑えられるという話。
 
Figure 4b で、正常マウス(Mock) にPBS(タダの水みたいなもの)と、自然免疫応答を促進して抗腫瘍効果を持つ(らしい)Poly(I:C)を投与した場合の、腫瘍体積の時間変化の実験と
sgPd-l1 (single guide Pd-l1, 3'UTR に干渉してaberrant, つまりこの部分の遺伝子領域が機能しない状態を模倣している)マウスにPBS とPoly(I:C)を投与した実験
を行っている。

 
Mock では腫瘍移植後、9日目以降からPBS とPoly(I:C) で差が出始めている。これらの実験時点ではBrunner–Munzel test使っている。
Brunner–Munzel test を使っているなんて、少し、玄人臭がする(実際にここの研究室はドライガチ勢がいるとかいないとか。風のうわさ)。
sgPdl1 では、PD-L1 の発現亢進により、Poly(I:C) による抗腫瘍効果はなくなっている。Poly(I:C) はTLR? INF? を介したシグナルがうんぬんかんぬんで私が完全に弱い免疫領域なので教えてエロい人。
載せていないけどFigure 4c では、Mock のPoly(I:C) はCD8+ 細胞(細胞傷害性T細胞といい、悪いやつをやっつける代表的なT細胞)がよく増えていて、局所における抗腫瘍効果が発揮されているのがわかる一方で、sgPdl1 のPoly(I:C) はCD8+ がコントロールPBS レベルしかおらず、T細胞による抗腫瘍効果がなくなるのもわかる、という感じ。
 
しかし、こういう図を見ると毎回思うのが、時系列でデータをとっているのに、各点でのt検定(今回はBrunner–Munzel test だけど)を使っているのが、なんかもにょる。
t検定をするときは、前後の時系列からは独立であってほしいと思う。しかも、データをとってこのように図示している以上、頭では時系列を思っているはずなので、その仮定に適した検定というものがあってほしい…
sgPd-l1 の実験の方は、たぶん、各点でt検定をせずに、増殖速度で回帰したら多分、有意差が出る。腫瘍体積の増殖速度は、r^3に比例するから、対数を取れば直線回帰できるだろうし、見た目でも、Poly(I:C) でちょっとは抗腫瘍効果がある感じ…まあ、この実験で示したいのはこれじゃないからいいんだけども。
 
Mock で9日目から体積に差が出始める、というのも、モデルは難しいから、ノンパラメトリック的に、各点でt検定(ここはパラメトリックでもノンパラでもよい)にせざるを得ないのは、仕方ないと思う。
 
こんな感じに近い。

 
シミュレーションした。けどいい検定が思いつかなかったので結局t検定している。
f(x)=x^2exp(-x/a)/c-b
の関数を容易して、適当に極大点をとるように係数a, b, cを調整する。
x=0,1で同じ値を取るように適当に係数を調整する。
こんな感じ。2日目で差が出るのは、まあ、見た目から明らかなのだが、これを各点でt検定は、はやり、もにょる。

 
x=1ではどちらの曲線も値は同じで、緑はまだ増加傾向、黄色はx=2で極大になる。
イデアとしては、ふたつの群を一緒くたにしてリサンプリングか、ふたつの折れ線が作る面積がどれほど増えるか、なるだろうか。
そうすれば、ふたつの曲線が異なることと、変化点がどこからか、が検出できるのかも知れない。

x <- seq(0, 10, length=100)

b <- 0.5
g <- exp(-1/2.5)/exp(-1)
f0 <- function(x) ((x)^2)*exp(-x/1) + b
f1 <- function(x) ((x)^2)*exp(-x/2.5)/g + b

y0 <- replicate(n, f0(d) + rnorm(length(d), sd=0.2))
y1 <- replicate(n, f1(d) + rnorm(length(d), sd=0.2))
y01 <- list(y0, y1)

d <- 0:5
n <- 10
cols <- c("yellow", "green")
yl <- c(0, 3.5)
matplot(x, cbind(f0(x), f1(x)), type="l", xlab="day", ylab="Measure", xlim=range(d), ylim=yl, col=cols, lty=1, lwd=3)
for(i in seq(y01)) boxplot(t(y01[[i]]), add=TRUE, boxwex=0.12, at=d+0.1*c(-1, 1)[i], axes=FALSE, col=cols[i])


p <- mapply(function(i) t.test(y01[[1]][i,], y01[[2]][i,])$p.value, seq_along(d))
p <- round(replace(p, p < 0.001, 0.001), 3)

matplot(d, sapply(y01, rowMeans), type="o", pch=15, xlab="day", ylab="Measure", ylim=yl, col=cols, lty=1, lwd=3)
for(i in seq(y01)) boxplot(t(y01[[i]]), add=TRUE, boxwex=0.12, at=d+0.1*c(-1, 1)[i], axes=FALSE, col=cols[i])
text(d, apply(do.call(cbind.data.frame, y01), 1, max), c("***", "*", "NS")[cut(p, c(1, 0.05, 0.001, 0))], pos=3)
legend("topleft", legend=paste(c("NS", "*", "***"), c("Not significant", "p < 0.05", "p < 0.001"), sep="\t"))

 
stanの神がやってた。
二つの時系列データの間に「差」があるか判断するには - StatModeling Memorandum