差分を取る

• 昨日の記事(こちら)で、介入には期間があることがあることを書き、観察は時点であることを書いた
• 期間を持つ介入は連続する期間について、それぞれの値の介入を行うものとする
  • 
• 他方、観察はのように、個の時点で観察値があるとする
• 介入と観察は、同じ時空間(同じ個人に流れる時間)にあるので、その間に関係があるかを量に変換することができる
• 介入は時間のすべてに連続的に値があり、観察は離散的に値がある
• 観察の折れ線グラフを考え、その傾きを考えると、観察量にも時間軸上に連続的な値を与えることができる
• ２つの尺度がペアをなしたので、ここに関連があるかを統計量に変換することにしよう
• 暫定版

k<-4
h<-5


K<-sort(runif(k+1))
H<-sort(runif(h))

V<-sort(runif(h))
U<-rep(1:0,k)[1:k]

MakeTandemPeriod<-function(t){
	tmpt<-sort(t)
	cbind(tmpt[1:(length(t)-1)],tmpt[2:length(t)])
}


MakeDiff<-function(v,t){
	d<-diff(v)
	tdiff<-diff(t)
	Period<-MakeTandemPeriod(t)
	list(d=d/tdiff,Period=Period)
}

Ks<-MakeTandemPeriod(K)
hout<-MakeDiff(V,H)
Hs<-hout$Period
Vs<-hout$d
Us<-U

Us
Ks
Vs
Hs

Therapy<-list(v=Us,t=Ks)
Phenotype<-list(v=Vs,t=Hs)

Commons<-function(t,p){
	rangeT<-range(t$t)
	rangeP<-range(p$t)
	st<-max(rangeT[1],rangeP[1])
	end<-min(rangeT[2],rangeP[2])
	period<-sort(c(t$t[,1],t$t[length(t$t[,1]),2],p$t[,1],p$t[length(p$t[,1]),2]))
	period<-period[which(period>=st & period<=end)]
	period<-MakeTandemPeriod(period)
	t.value<-p.value<-rep(0,length(period[,1]))
	comp.t<-outer(period[,1],t$t[,1],FUN=">=")
	comp.p<-outer(period[,1],p$t[,1],FUN=">=")
	max.t<-apply(comp.t,1,sum)
	max.p<-apply(comp.p,1,sum)
	list(t.value=t$v[max.t],p.value=p$v[max.p],period=period)
}

IntegralValue<-function(cmm){
	sum(cmm$t.value*cmm$p.value*(cmm$period[,2]-cmm$period[,1]))
}

k<-6
K<-sort(runif(k+1))
U<-rep(1:0,k)[1:k]
Ks<-MakeTandemPeriod(K)

Nt<-100
tmpt<-seq(from=0,to=1,length=Nt)
kv<-1
tmpV<-rep(0,length(tmpt))
for(i in 1:length(tmpt)){
	for(j in 1:length(Ks[,1])){
		if((Ks[j,1]-tmpt[i])*(Ks[j,2]-tmpt[i])<0){
			#if(U[j]==1){
				tmpV[i]<-kv*U[j]
			#}else{
				#tmpV[i]<-kv*(-1)
			#}
			
		}
	}
}



h<-Nt/2
selectedt<-sort(sample(1:length(tmpt),h))
H<-tmpt[selectedt]
V<-cumsum(tmpV[selectedt])
Ks<-MakeTandemPeriod(K)
hout<-MakeDiff(V,H)
Hs<-hout$Period
Vs<-hout$d
Us<-U

Us
Ks
Vs
Hs

Therapy<-list(v=Us,t=Ks)
Phenotype<-list(v=Vs,t=Hs)

cmm.out<-Commons(Therapy,Phenotype)
Stat.Ori<-IntegralValue(cmm.out)

Niter<-1000

Stat.Perm<-rep(0,Niter)

for(i in 1:Niter){
	tmpP<-Phenotype
	tmpP$v<-Phenotype$v[sample(1:length(Phenotype$v))]
	tmp.cmm.out<-Commons(Therapy,tmpP)
	Stat.Perm[i]<-IntegralValue(tmp.cmm.out)
}

par(mfcol=c(1,2))
plot(c(Stat.Ori,sort(Stat.Perm)))
abline(h=Stat.Ori,col=2)

xlim=range(c(Ks,Hs))
plot(c(t(Phenotype$t)),c(t(cbind(Phenotype$v,Phenotype$v))),type="b",col=2,xlim=xlim)
par(new=TRUE)
plot(c(t(Therapy$t)),c(t(cbind(Therapy$v,Therapy$v))),type="b",xlim=xlim)