フェノタイプの整理メモ６　確率的に考える１ - ryamadaの遺伝学・遺伝統計学メモ

死亡というフェノタイプ
- 0/1のフェノタイプ
- すべての人が必ず1に転ずるフェノタイプ
- 0歳時を例外とすれば、加齢とともに率が上昇し、ピークを迎えて減少するという一峰性の確率分布となる
- 生命表を例に
  - 年齢別死亡者数をプロットしている
  - ヒトが生まれたら、何歳で亡くなる確率がいくつか、というプロットがほしいのだが、細かいところは無視して、おおまかにこんな感じ、ということで

# だんだんに事件が蓄積するモデルによる、死亡年齢分布
# 平成１８年簡易生命表（男）からhttp://www.mhlw.go.jp/toukei/saikin/hw/life/life06/hyo-m.html
death<-c(280,41,30,21,16,14,13,12,10,9,9,9,10,13,17,22,28,36,44,50,54,57,59,60,61,62,63,65,68,70,71,73,76,80,85,90,96,104,113,124,135,146,159,176,195,216,237,257,279,302,329,359,392,430,471,516,566,614,661,712,771,836,903,968,1023 ,1080 ,1153 ,1242 ,1344 ,1457 ,1585 ,1731 ,1895 ,2068 ,2244 ,2426 ,2608 ,2783 ,2950 ,3110 ,3267 ,3405 ,3523 ,3616 ,3692 ,3750 ,3768 ,3735 ,3629 ,3463 ,3244 ,2980 ,2683 ,2366 ,2040 ,1719 ,1414 ,1133 ,884,669,492,350,240,159,101,138)
plot(1:length(death),death/sum(death),type="l",xlab="Age(year)",ylab="Fraction of death")

死亡というフェノタイプをシンプルにモデル化する
- 一定の確率で生起する事件が特定回数、繰り返されると、「死亡」というフェノタイプに転じるというモデル

２つのプロットの大きな違いは、実データが左右非対称の程度が強いのに対して、シミュレーションは左右対称度が高いことであるが、傾向はつかめている、と考えることにする

par(mfcol=c(1,2))
plot(1:length(death),death/sum(death),type="l",xlab="Age(year)",ylab="Fraction of death")

NL<-1 # 座位数
TL<-120 # 一生は120年
Fs<-1 # イベントが起きる平均時間間隔(年)(平均して毎年、１回、事件が起きる)
# 86回蓄積すると「死亡フェノタイプ」に転ずるとする
# 86は生命表でのピーク年齢を取った
ThresNums<-c(86)
par(ask=TRUE)
for(fi in 1:length(Fs)){
	F<-Fs[fi]
for(ti in 1:length(ThresNums)){
	ThresNum<-ThresNums[ti]
# Niter人でシミュレーション
	Niter<-100000
	DxTimes<-matrix(0,Niter,NL)
	for(ii in 1:Niter){
		NR<-ThresNum
		Rmat<-matrix(0,NL,NR)
		for(i in 1:NL){
			Rmat[i,]<-rexp(NR,1/F[i])
		}
		CumEventMat<-t(apply(Rmat,1,cumsum))
		DxTime<-rep(0,NL)
		for(i in 1:NL){
			DxTime[i]<-CumEventMat[i,ThresNum[i]]
		}
		DxTimes[ii,]<-DxTime
	}
hist(DxTimes[,1],breaks=0:ceiling(max(DxTimes[,1])),xlim=c(0,TL),
	xlab="Age of Dx(year)",ylab="Fraction of Dx",col=2,freq=FALSE)
}
}
par(mfcol=c(1,1))