# 状態推移

• Rmdファイル
```# 状態推移 Status-transition

## 基本 Basics

- 確率的に何かが起きる Consider a stochastic process
- 次に何が起きるかは「今」の影響を受けるが、「前」の影響を受けない "Next event" depends on "Now" but not on "Past".

## 簡単な例

### ただの確率事象 A simple stochastic process

ある確率で、機能性細胞が死滅するとする。その確率はpである。初め、細胞は全部でn個あった。生存細胞数をシミュレーションする。

Assume that cells die stochasitically. The probability of death is p. Simulate for n cells at the begininng.
```{r}
n <- 1000
p <- 0.1
s.or.d <- sample(0:1,n,replace=TRUE,prob=c(1-p,p))
s.or.d[1:100]
sum(s.or.d) # 死細胞数 No. dead cells
sum(s.or.d)/n # 死細胞割合 Fraction of dead cells
```
### 増えたり減ったり Increase and decrease

Cells may die but living cells will divide into two. Both are stochastic, with probabilities p and q, respectively.
At the begininng there are n cells.

いくつの細胞が死に、いくつの細胞が分裂するかは、その時々の生細胞数に依存するので、「今」の影響を受けている。
しかしながら、一段階前に何個の細胞が生きていたかには影響を受けていない。
No. death/division depend on the No. cells "NOW", but not on "One period previous".

```{r}
p <- 0.05
q <- 0.03
weeks <- 1:300
# No. living cell from 1st through 100th week
n.sv.1 <- rep(0,length(weeks))
n.sv.1[1] <- 100
for(i in 2:length(weeks)){
# death
tmp1 <- sample(0:1,n.sv.1[i-1],replace=TRUE,prob=c(1-p,p))
n.death <- sum(tmp1)
# No. living cells for now
tmp2 <- n.sv.1[i-1] - n.death
# Cell division
tmp3 <- sample(0:1,tmp2,replace=TRUE,prob=c(1-q,q))
# No. new born cells
tmp4 <- sum(tmp3)
# No. living cells for now
n.sv.1[i] <- tmp2 + tmp4
}

plot(n.sv.1)
```

## 加速する Acceration

When there are m cells, each cell has to do 1/m job and the death rate is positively correlated with quantity of job. Division rate is independent from job.

```{r}
p <- 0.05 # p depends on number of living cells
q <- 0.03
weeks <- 1:300
# No. living cell from 1st through 100th week
n.sv.2 <- rep(0,length(weeks))
n.sv.2[1] <- 100
for(i in 2:length(weeks)){
# death rate
p. <- p + 1/n.sv.2[i-1]
# p. <= 1
p. <- min(p.,1)
# death
tmp1 <- sample(0:1,n.sv.2[i-1],replace=TRUE,prob=c(1-p.,p.))
n.death <- sum(tmp1)
# No. living cells for now
tmp2 <- n.sv.2[i-1] - n.death
# Cell division
tmp3 <- sample(0:1,tmp2,replace=TRUE,prob=c(1-q,q))
# No. new born cells
tmp4 <- sum(tmp3)
# No. living cells for now
n.sv.2[i] <- tmp2 + tmp4
}

plot(n.sv.2)
```

Plot two simulation results together.

```{r}
matplot(cbind(n.sv.1,n.sv.2),type="b")
```

## 「今」だけとは限らないが、それでもマルコフ連鎖 MCMC may depend on "NOW" and previous

Baby cells just born do not divide but they become dividable one period later.

この場合は、生細胞を「赤ん坊」か「そうでないか」の２通りで記録する必要がある。

You have to distinct baby cells and non-baby cells.

```{r}
p <- 0.05
q <- 0.03
weeks <- 1:300
# No. living cell from 1st through 100th week
# baby and non-baby
n.sv.3.nb <- n.sv.3.bb <- rep(0,length(weeks))
n.sv.3.nb[1] <- 100
for(i in 2:length(weeks)){
# death of non-baby
tmp1.nb <- sample(0:1,n.sv.3.nb[i-1],replace=TRUE,prob=c(1-p,p))
n.death.nb <- sum(tmp1.nb)
# death of baby
tmp1.bb <- sample(0:1,n.sv.3.bb[i-1],replace=TRUE,prob=c(1-p,p))
n.death.bb <- sum(tmp1.bb)
# No. living non-baby cells for now
tmp2 <- n.sv.3.nb[i-1] - n.death.nb
# Cell division
tmp3 <- sample(0:1,tmp2,replace=TRUE,prob=c(1-q,q))
# No. new born cells
tmp4 <- sum(tmp3)
# No. living non-baby cells. Living baby cells are now non-baby
n.sv.3.nb[i] <- tmp2 + n.sv.3.bb[i-1]-n.death.bb
# No. newborns
n.sv.3.bb[i] <- tmp4
}

plot(n.sv.3.nb+n.sv.3.bb)
matplot(cbind(n.sv.3.nb,n.sv.3.bb,n.sv.3.nb+n.sv.3.bb),type="l")
```

## 行列を使う Use matrices

ある疾患には、４つの状態があるという。急性軽症・急性重症・慢性・癌化の４つである。この４状態に、健康・死亡の２状態を加えて、全部で６状態を考える。

Assume a disease has 4 statuses, acureLight,acuteSevere,chronic,cancer. Along these 4, two more statuses, healthy and dead, should be considered.

Call 4+6 statuses with 1,2,...,6 as 1:healthy, 2:acuteLight, 3:acuteSevere, 4:chronic, 5:cancer, 6:dead.

Healthy individual can become acuteLight/acuteSevere/Dead(dead from some reason away from this disease). Let p12,p13,p16 denote their probability.

Healthy individual can not get "chronic/cancer" directry, meaning p14=p15=0. (S)he may stay healthy with probability p11.

このとき、p11 + p12 + p13 + p14 + p15 + p16=1であって、p1iはすべて0以上である。これを、状態１からの推移確率ベクトルと呼ぶ。

p = (p11,...,p16) is called the transition vector from status 1 and p11 + ... + p16 = 1 and p1i >= 0.

Each status has its own transition vector.

6つの状態の状態推移ベクトルを合わせると、6x6行列が作れる。これを状態推移行列という。

With these 6 transition vectors, you can make a 6x6 matrix, called "transition matrix".

Using this transition matrix, history of probabilities in 6 statuses can be simulated.

Transition probability is per-year one.

まず行列を作る。Make a matrix.
```{r}
# From 1:healthy
p1 <- c(0.9,0.08,0.0001,0,0,0.0199)
# From 2:acuteLight
p2 <- c(0.95,0,0,0.04,0,0.01) # 治るか慢性化するか、死亡するか
# From 3:acuteSevere
p3 <- c(0.8,0,0,0.01,0,0.19) # 死亡率0.19
# From 4:chronic
p4 <- c(0.01,0,0,0.95,0.03,0.01) # 癌化する
# From 5:cancer
p5 <- c(0,0,0,0,0.9,0.1)
p6 <- c(0,0,0,0,0,1)

# 行列化 matrix
P <- cbind(p1,p2,p3,p4,p5,p6)
P
# 全列が足して１になるかを確認 check sum of columns being 1
apply(P,2,sum)
```

Initially you are healthy(1), then your status vector is (1,0,0,0,0,0).

```{r}
v1 <- c(1,0,0,0,0,0)
# １年後 One-year later
v2 <- P %*% v1
v2
# ２年後、two-year later
v3 <- P %*% v2
v3

# x年後、x-year later
x <- 1:50
v.all <- matrix(0,length(x),length(v1))
v.all[1,] <- v1
for(i in 2:length(x)){
v.all[i,] <- P %*% v.all[i-1,]
}
matplot(v.all,type="l")
```

## 治療介入してみよう Evaluate therapeutic interventions

Assume a treatment A that increases p41 from 0.01 to 0.1.

```{r}
PA <- P
PA[1,4] <- 0.1
PA[4,4] <- PA[4,4]-(0.1-0.01)
PA
```

```{r}
v.all.A <- matrix(0,length(x),length(v1))
v.all.A[1,] <- v1
for(i in 2:length(x)){
v.all.A[i,] <- PA %*% v.all.A[i-1,]
}
# 治療法Aの効果を死亡状態の割合で比較する
matplot(x,cbind(v.all[,6],v.all.A[,6]),type="l")
```

Assume a treatment B that cures cancer, p51 = 0.1.
```{r}
PB <- P
PB[1,5] <- 0.1
PB[5,5] <- PB[5,5]-0.1
PB
```

```{r}
v.all.B <- matrix(0,length(x),length(v1))
v.all.B[1,] <- v1
for(i in 2:length(x)){
v.all.B[i,] <- PB %*% v.all.B[i-1,]
}
# 治療法Aの効果を死亡状態の割合で比較する
matplot(x,cbind(v.all[,6],v.all.A[,6],v.all.B[,6]),type="l")
```
```
• カメ表紙関数
```make.kame.hyoushi <- function(ttl="",n.kame=30,kame.col="blue",sub=""){

t <- seq(from=0,to=1,length=1000)*2*pi
kora <- cbind(cos(t),sin(t))
s <- 1/4
teashi <- kora*s
thetas <- c(0,1,2,4,5)/6*2*pi
teashi.ctr <- cbind(cos(thetas),sin(thetas))
teashi.coords <- matrix(0,nrow=0,ncol=2)
for(i in 1:length(thetas)){
teashi.coords <- rbind(teashi.coords,cbind(teashi.ctr[i,1]+teashi[,1],teashi.ctr[i,2]+teashi[,2]))
}
teashi.out <- teashi.coords[which(apply(teashi.coords^2,1,sum)>1),]
shippo.st <- c(cos(pi),sin(pi))
shippo.s <- 0.1
shippo.l <- 0.4
shippo.x <- seq(from=0,to=-shippo.l,length=100)
shippo.y <- shippo.s * sin(shippo.x*2*pi/shippo.l)
shippo.coords <- cbind(shippo.st[1] + shippo.x,shippo.st[2]+shippo.y)
kame <- rbind(kora,teashi.out,shippo.coords)
#plot(rbind(kora,teashi.out,shippo.coords))

b <- 0.9
kame2 <- kame
kame2[,2] <- kame2[,2]*b
#plot(kame2,col=grey(1),axes=FALSE)
#points(kame2,cex=0.1,pch=20)

X <- matrix(0,nrow=0,ncol=2)
for(i in 1:n.kame){
theta <- runif(1)*2*pi
R <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),2,2)
new.kame <- t(R %*% t(kame2))
mv <- runif(2)*n.kame/2
new.kame <- cbind(new.kame[,1]+mv[1],new.kame[,2]+mv[2])
X <- rbind(X,new.kame)
}
jpeg(filename = paste(ttl,".jpeg"), width = 5000, height = 5000,
quality = 100, pointsize = 180,bg = "white", res = NA,
restoreConsole = TRUE)
plot(X,col=grey(1),axes=FALSE,xlab="",ylab="",main=paste(ttl,"\n",sub))
points(X,cex=0.5,pch=20,xlim=range(X),ylim=range(X),col=kame.col)
dev.off()
}

# カメ頭数
n.kame <- 15
# カメ色
kame.col <- "orange"
# タイトル
ttl <- "状態推移シミュレーション"
subttl <- "〜臨床的例を用いて〜"
make.kame.hyoushi(ttl,n.kame,kame.col,subttl)
```