フェノタイプの整理メモ８ 確率的に考える３

• 浸透率
  • 発病する遺伝因子を有する個体のうち実際に発病する割合
  • こちらの記事で、「イベント累積型発病」の場合に、発病年齢が人生の後半になり、ヒトの寿命上限までに、十分にイベントが累積しないことがあるとこを示した
  • 実際に浸透率を考える場合には、「イベント累積型発病」が「死亡というイベント」よりも先に起きる個体の割合を考えることとなる
  • 以下のプロットはその様子
    • 左側が「死亡しない場合」、右側が「死亡イベント以前の発病」の場合

# ペネトランス

# 発病前に死亡するような条件のイベント生起確率と累積回数の疾病と
# 死亡とを組み合わせて、年齢別発病人数比率、その累積を見る

NL<-1 # 座位数
TL<-120 # 一生は120年

# イベントがしょっちゅう起きて(均質化して)いる場合には
# Dxのタイミングにはばらつきが小さい

Fs<-0.01*2^(10) # イベントが起きる平均時間間隔(年)

ThresNums<-10
par(ask=TRUE)
for(fi in 1:length(Fs)){
	F<-Fs[fi]
	ThresNum<-ThresNums[1]

	Niter<-1000
	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)



}
plot(ecdf(DxTimes[,1]),xlim=c(0,TL))


# 死亡モデル
NL<-1 # 座位数
TL<-120 # 一生は120年

F<-c(1) # イベントが起きる平均時間間隔(年)

ThresNums<-86
par(ask=TRUE)
for(ti in 1:length(ThresNums)){
	ThresNum<-ThresNums[ti]

	Niter<-1000
	DxTimes2<-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]]
		}
		DxTimes2[ii,]<-DxTime
	}

	hist(DxTimes2[,1],breaks=0:ceiling(max(DxTimes2[,1])),xlim=c(0,TL),
	xlab="Age of Dx(year)",ylab="Fraction of Dx",col=2,freq=FALSE)


}

DxBeforeDeath<-which(DxTimes[,1]<DxTimes2[,1])

DxTimes3<-rep(Inf,length(DxTimes[,1]))
DxTimes3[DxBeforeDeath]<-DxTimes[DxBeforeDeath,1]

par(mfcol=c(1,2))
	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=TRUE,ylim=c(0,20))

	hist(DxTimes3,breaks=0:ceiling(max(DxTimes[,1])),xlim=c(0,TL),
	xlab="Age of Dx(year)",ylab="Fraction of Dx",col=2,freq=TRUE,ylim=c(0,20))
par(mfcol=c(1,1))

par(mfcol=c(1,2))
plot(ecdf(DxTimes[,1]),xlim=c(0,TL),ylim=c(0,1),cex=0.1)

plot(ecdf(DxTimes3),xlim=c(0,TL),ylim=c(0,1),cex=0.1)
par(mfcol=c(1,1))