■

• 疾病を次のようにモデルにする
  • 病気の活動が潜んでいる
  • その活動性を測定することができる
  • 活動性の測定値によって、治療を一定期間投入する
  • 完全寛解にはいたらず、治療をやめると、再燃する
  • 疾患活動性がある程度高いと、不可逆変化も残る
• こんな場合に経時的に、疾患活動性・不可逆変化・治療の有無を記録しているとする

# 変化を作る
# 個人
# 状態
# テスト
# 治療

Tmax<-100
Tlen<-1000
T<-seq(from=0,to=Tmax,length=Tlen)

# 治療無しの場合の経過
V<-(runif(Tlen))
Vnatural<-cumsum(V)
# VがVtより大きいと治療が一定期間、入る
# ただし、そのVtはぶれている
Vt<-50-runif(Tlen)*20
Vperiod<-100
# 治療効果
r<-0.95
Vintervention<-rep(0,Tlen)
Vintervention[1]<-V[1]
Intervention<-rep(0,Tlen)
for(i in 2:Tlen){
	if(max(Vintervention[max(1,i-Vperiod-1):i-1]-Vt[max(1,i-Vperiod):i])>0){
		Vintervention[i]<-Vintervention[i-1]*r+V[i]
		Intervention[i]<-1
	}else{
		Vintervention[i]<-Vintervention[i-1]+V[i]
	}
}
# 不可逆変化閾値
Ct<-30
tmp<-Vintervention
tmp[tmp<Ct]<-0
tmp[tmp>=Ct]<-runif(1)
Cumul<-cumsum(tmp)

CumulSt<-Cumul/max(Cumul)
VinterventionSt<-Vintervention/max(Vintervention)
Tst<-T/max(T)
X<-cbind(CumulSt,VinterventionSt,Intervention,Tst)

matplot(X,type="l")

• 複数の人のデータだと…

# 変化を作る
# 個人
# 状態
# テスト
# 治療
library(MCMCpack)
Tmax<-100
Tlen<-2000
T<-seq(from=0,to=Tmax,length=Tlen)

Nsample<-5
VtParam<-sample(10:20,Nsample,replace=TRUE)

Disease<-function(ID=1,T,Tmax,Tlen,VtParam){
	# 治療無しの場合の経過
	V<-runif(Tlen)
	#V<-(rnorm(Tlen))^2
	#V<-c(rdirichlet(1,rep(1,Tlen)))*100
	Vnatural<-cumsum(V)
	# VがVtより大きいと治療が一定期間、入る
	# ただし、そのVtはぶれている
	Vt<-150-(rnorm(Tlen))*50
	Vperiod<-100
	# 治療効果
	r<-0.95
	Vintervention<-rep(0,Tlen)
	Vintervention[1]<-V[1]
	Intervention<-rep(0,Tlen)
	for(i in 2:Tlen){
		if(max(Vintervention[max(1,i-Vperiod-1):(i-1)]-Vt[max(1,i-Vperiod-1):(i-1)])>0){
			Vintervention[i]<-Vintervention[i-1]*r+V[i]
			Intervention[i]<-1
		}else{
			Vintervention[i]<-Vintervention[i-1]+V[i]
		}
	}
	# 不可逆変化閾値
	Ct<-30
	tmp<-Vintervention
	tmp[tmp<Ct]<-0
	tmp[tmp>=Ct]<-runif(1)
	Cumul<-cumsum(tmp)

	CumulSt<-Cumul/max(Cumul)
	VinterventionSt<-Vintervention/max(Vintervention)
	Tst<-T/max(T)
	X<-cbind(rep(ID,Tlen),CumulSt,VinterventionSt,Intervention,Tst)
	X
}

Out<-list()
for(i in 1:Nsample){
	Out[[i]]<-Disease(ID=i,T,Tmax,Tlen,VtParam[i])
}

Out2<-NULL
for(i in 1:Nsample){
	tmp<-Out[[i]][sample(1:length(Out[[i]][,1]),length(Out[[i]][,1])/100),]
	print(tmp)
	Out2<-rbind(Out2,tmp)
}

matplot(Out2,type="l",ylim=c(0,1))
pairs(Out2)