シミュレーションして生存曲線をプロットする

MikuHatsune2011-09-16

遺伝学コースの、ただRだけをひたすらやる集中講義
やりたい検定、したい解析を見つける。
そのためにデータセットをシミュレーションして作る。
ということで、生存曲線を描こう。
 
とりあえず、初心者が使える関数などを知らずに生存解析とプロットまでにたどり着くまでの
とてもちいさな、とてもおおきな、とてもたいせつな、あいとゆうきのおとぎばなし
 
気になるブログ。 

データの確認。
最終的には、曲線同士を比較したい。
g1g2の2群比較とする。
それぞれの群は行列にデータが収まっている。
行はサンプルID。
列1は、サンプルの状態が、生存なら0、死亡なら1、が入っている。
列2は、時刻が入っている。今回は連続量とする。
サンプル状態が0なら、少なくともその時刻には生存していた、1なら、死亡時刻、とする。
 
シミュレーションデータの作成。
生きている人について:中央値を定めて正規分布に従う期間。観察期間。
死んでいる人について:ある時刻Tに、ほぼ均等に死亡時刻が決まる。一様乱数。
 
という訳で中身はできた。あとは、生存解析をしてくれる関数を探そう。
グーグル先生に聞くと、Kaplan Meier法とかいうのがあるらしいので、Rでもsurvival packageの中にSurvという関数があるらしい。
そこで、いろいろ検索すると、先にやっている人に行き着く。あ、これ使える、と思ってやろうとすると、関数がほしがる形にデータを整えなければならないようで、この人のやり方を見ると
別に群を2つの行列に分ける必要はなく、行列内にグループ分けを表す列3を作っておけばいいようだ。
という訳で、cbindおよびrbind関数を使って細工する。
あとは、Surv関数に引数をいれればよさげだ。
 
帰無仮説としては、生存期間に差はない、として‥‥ここらあたりを聞き逃した。
 
と思って関数を動かすと、p値を返してくれない。
これまた調べると、logrank検定ということをすればいいらしい。
これをするのはsurvdiffらしい。
だがしかし、p.valueが取り出せない事態。
ここでもグーグル先生の協力を得て、p.valueが取り出せる関数を発見した。ここをコピペして使用する。
 
あとは、このシミュレーションを複数回行い、p値が一様に分布しているか確かめよう。

# 人数の設定
n1<- 20
n2<- 20
# 生存中央値、分散の設定。月の気分。
t_mean<- 6
t_sd<- 0.8
t_death<- 6
# 繰り返し回数
Niter<- 1000
p_values<- rep(0,Niter)

# 生死の振り分け確率
p1<- c(0.5,0.5)
p2<- c(0.5,0.5)

# 生死の振り分け
s1<- sample(0:1,size=n1,replace=TRUE)
s2<- sample(0:1,size=n2,replace=TRUE)
# 2群の作成
g1<- matrix(0,nr=n1,nc=2)
g2<- matrix(0,nr=n2,nc=2)
g1[,1]<- s1
g2[,1]<- s2
# 生きている人を抜き出して‥‥
g1_0<- sum(g1[,1]==0) #length(which(g1[,1]==0))
g2_0<- sum(g2[,1]==0) #length(which(g2[,1]==0))
# 時刻を正規分布から作成して‥‥
g1_0_t<- rnorm(g1_0,t_mean,t_sd)
g2_0_t<- rnorm(g2_0,t_mean,t_sd)
# 入れる。
g1[which(g1[,1]==0),2]<- g1_0_t
g2[which(g2[,1]==0),2]<- g2_0_t
# 死んでいるほうもやる。
g1_1<- sum(g1[,1]==1) #length(which(g1[,1]==1))
g2_1<- sum(g2[,1]==1) #length(which(g2[,1]==1))
# 時刻を一様分布から作成して‥‥
g1_1_t<- runif(g1_1,0,t_death)
g2_1_t<- runif(g2_1,0,t_death)

g1[which(g1[,1]==1),2]<- g1_1_t
g2[which(g2[,1]==1),2]<- g2_1_t

# 時刻がそれぞれの分布に従っているかヒストグラム。
par(mfcol=c(1,2))
hist(g1[which(g1[,1]==0),2],xlab="time",main="Alive people observing time")
hist(g1[which(g1[,1]==1),2],xlab="time",main="Dead people dead time")

# boxplotなんてしてみる。
par(mfcol=c(1,1))
boxplot(g1[which(g1[,1]==1),2],g1[which(g1[,1]==0),2])

# 生存解析用にあったデータ型にする。
g12<- rbind(cbind(g1,rep(1,n1)),cbind(g2,rep(2,n2)))

library(survival)
result<- survfit(Surv(g12[,2],g12[,1])~g12[,3],type="kaplan-meier")

# この関数のmethodには3種類あるらしいけど何か違うのかな。
par(mfcol=c(1,3))
surv_type<- c("kaplan-meier", "fleming-harrington", "fh2")
for(j in 1:length(surv_type)){
	plot(survfit(Surv(g12[,2],g12[,1])~g12[,3],g12,type=surv_type[j]),main=surv_type[j])
}

# p.valueを返してくれない関数。
survdiff(Surv(g12[,2],g12[,1])~g12[,3])

Res<- logrank(g12[,3],g12[,1],g12[,2])