2013-11-20

エントリータイミングのパターン

生存解析

こちらで生存解析にあたって、初期にエントリーバーストがあって、その後エントリーが定常化することについて書かれている
そういえば生存解析について、以前にも書いた…(こちら)
この記事にアイディアを得て、どんなモデルが作れるかを考えてみる
できるだけ、要素を分解したい
- (1)インターベンション割り付け
- (2)イベント発生までの期間
- (3)ドロップアウトまでの期間
- (4)エントリーの時刻
に分けてそれぞれに確率密度分布を使ったりするのがよいだろうか
- (1)割り付けはsample()
- (2),(3)は指数分布
- (4)はどうして初期バーストがおきて、その後、定常化するかをモデル化し、スタディ開始からの時刻ごとのエントリー発生の相対比を出せるようにして、それに基づいてMetropolis-Hasting sampling(『直接サンプリングするのが難しい確率分布から統計標本の配列を生成』)でよいだろうか
- イベント発生メディアンの確認

median(ev.time[itv.v[[1]]])
median(ev.time[itv.v[[2]]])

- ドロップ時刻(全員、いつかはドロップするものとすれば場合分けが不要)

- イベントもドロップも全員に起こさせる。どちらが先かを解析には使えばよい。２群を横軸で分け、イベントとドロップの時刻を縦軸で表す。線はドロップまで引いてある。線より上に離れて打たれた点は、ドロップごとにおきるイベントで、これは観測されない

- エントリーはスタディ開始前には、候補者がたまる。いざ始めると、すぐにエントリーを開始できる候補者とそうでない候補者が出る。施設の準備状況などにもよるだろう。その開始前のたまった人数を与え、また、その後に新規に候補になる人数(新規診断人数)を時間の関数で与える。施設の準備状況は加速度的に増える。いくつかのモデルが考えられるが、指数関数と時間のr乗とを使って、初めは少なく、ついで高頻度エントリー期を作り、定常化するようにパラメタを振った
- その累積エントリー数のグラフ

- 各時刻のエントリーはこれの微分

- これに応じてサンプリングしたときのエントリー時刻の分布は

- エントリー時刻を考慮して、イベント・ドロップの打点をやりなおすと

# いくつのパラメタが必要か
# 標本作成数(複数回シミュレーションするにしろ一気につくる…)
n <- 1000
# インターベンションの割り付け
# インターベンションタイプ
itv.type <- 1:2
itv <- sample(itv.type,n,replace=TRUE)
itv <- sort(itv)
# インターベンションタイプ別ベクトル
itv.v <- list()
for(i in 1:length(itv.type)){
	itv.v[[i]] <- which(itv==itv.type[i])
}
# インターベンション開始からイベントまでの時間
# インターベンションタイプ別のポアソン分布のパラメタは
# 指数分布の累積分布関数が1-e^(-kx)なので、1-e^(-k*MST)=0.5によって
# kを定め、それに応じた指数乱数を発生させ、その累積をとる
# インターベンション別のMST
mst <- c(17,10)
mst.k <- log(2)/mst
ev.time <- rep(0,n)
for(i in 1:length(itv.type)){
	tmp <- rexp(length(itv.v[[i]]),mst.k[i])
	ev.time[itv.v[[i]]] <- tmp
}


median(ev.time[itv.v[[1]]])
median(ev.time[itv.v[[2]]])

# インターベンション開始からイベントまでの時間
# １か月あたり0.05減ということは(0.95)^monthsが残っている人数
# (0.95)^months = 0.5となるのは
k <- log(0.5)/log(0.95)
# イベントまでの時刻での考え方と同じで、ドロップアウト時刻は
drop.time <- rexp(n,log(2)/k)
hist(drop.time)

plot(1:n,ev.time,col=itv,pch=20)
points(1:n,drop.time,col=itv+2,pch=20)
segments(1:n,rep(0,n),1:n,drop.time)

# さてエントリーのタイミングを変えてみよう
# 既診断者・既把握対象者がn.old人おり、
# 新規診断者・新規把握対象者がp /monthで増えるとする
# 既・新併せての対象者が実際にエントリーされるには指数関数的に
# 拡大する周知範囲と言うものが影響するものとする
# 周知範囲は、1-exp(-u * months^r)で決まるものとすれば
# ru month exp(-u months^r) * ( n+p month) + (1-exp(-u months^r)) * p がmonthsまでのエントリー数
# ある時刻monthsでのエントリー数はこれの微分
u <- 0.001
n.old <- 50
p <- 3
r <- 3
t <- seq(from=0,to=50,length=51)
cumsum.entry <- (1-exp(-u*t^r))*(n.old+p*t)
plot(t,cumsum.entry,type="l")

entry <- r*u * t *exp(-u *t^r) * ( n+p *t) + (1-exp(-u *t^r)) * p
plot(t,entry,type="l")

plot(t,entry,type="l",ylim=c(-5,30))
abline(h=0,col=2)
abline(h=p,col=3)

# このようなスタディ開始からの時刻別のエントリー確率で上記n人のエントリー時刻を
# ランダムに与えてみる
# 時刻別のエントリー数には最大値がある
# 解析的に求まるかもしれないが、
upper <- max(entry)*1.1
# としても、upperはその最大値より大きいだろう
# このupperを使ってMetropolis-Hastingサンプリングをする
# エントリー時刻の最大値(スタディ継続期間で決まる)end.T
end.T <- 50
entry.time <- rep(0,n)
for(i in 1:n){
	loop <- TRUE
	while(loop){
		tmp <- runif(1)*end.T
		tmp.prob <- r*u * tmp *exp(-u *tmp^r) * ( n+p *tmp) + (1-exp(-u *tmp^r)) * p
		if(runif(1) < tmp.prob/upper){
			entry.time[i] <- tmp
			loop <- FALSE
		}
	}
}

hist(entry.time)

# エントリータイムを考慮する
plot(1:n,ev.time+entry.time,col=itv,pch=20)
points(1:n,drop.time+entry.time,col=itv+2,pch=20)
segments(1:n,entry.time,1:n,drop.time+entry.time)