前期 医科学のためのデータシミュレーション

• こちらの続き
• 医学・生物学データの扱いに慣れるための、データシミュレーションセミナーの資料(医学研究科大学院コース　ゲノム遺伝学＠京大 2014:こちら)
• htmlファイルはこちら(Login as a guest でログインすると見られます)
• こちらを参考にすればhtml化可能です
• epub化も可能ですが、その手間が惜しい向きにはアマゾンkindleへ

Rを用いた医科学のためのデータシミュレーションの基礎2014(前半）

• 過去数日のRmdファイルのただのタンデム連結が以下。

# Rを用いた医科学のためのデータシミュレーションの基礎 Basics of Data Simulation for Medical Sciene with R

## 基本のき ABC

### 起動と終了。Start and quit

q()

### コピー・ペースト Copy and paste

この文書のコマンドの部分は"Ctrl" + "c"でクリップボードにコピーできて、それを自分のRの上で"Ctrl" + "v"でペーストできる。
You can copy the commands on this page with "ctrl" + "c" and paste them onto your R window with "ctrl" + v.

ウィンドウズは"ctrl"+...で、Macはアップル+...

Windows: "ctrl" + ... but Mac: "Apple" + ...

テキストエディタ上にペーストすれば保存できる。
You can paste the commands on your text editor and you can save the commands as a file.

次のコマンドをコピー・ペーストして実行してみよ。

```{r}
boxplot(count ~ spray, data = InsectSprays, col = "lightgray")
```

## データの扱い Handling of Data

### ベクトル Vectors

```{r}
v1 <- c(3,2.5,6)
# Length of vector = number of elements
length(v1)
v1.id <- 1:3
# Same
v1.id.1 <- 1:length(v1)
```

データを作ったらプロットする。Plot after you generate data.
```{r}
plot(v1)
plot(v1.id.1,v1)
```
プロットに少し情報を加える。Some information is added.
プロットを見比べて、どの指定が何をするかを確認する。Compare plots and understand which words determined which changes in the 2nd plot.
```{r}
plot(v1.id.1,v1,xlab="ID",ylab="value",main="test plot")
```

### 乱数列 Random number sequence
```{r}
# Length
n <- 100
v2 <- runif(n) # 0-1の一様乱数 random values in uniform distribution from 0 to 1
v2
```
プロットして理解する。Grab what they are by plotting them.
```{r}
plot(v2)
# ソートするとわかりやすい Sorted sequence might be more informative 
sorted.v2 <- sort(v2)
# 逆順にもソートできる Sorting in descreasing order is also useful.
inv.sorted.v2 <- sort(v2,decreasing=TRUE)
plot(sorted.v2)
plot(inv.sorted.v2)
```

### Help記事と検索 Help article and serching information

関数を知っていてその使い方を知りたいとき。When you know name of functions and you want to know how to use them.
```{r}
help(length)
help(runif)
```

キーワードを知っていてそれに関連する記事を検索したいとき。
When you know keywords and search related articles.
用語の一部を知っているとき。When you can spell a part of terms.

```{r}
?? "length"
?? "uniform"
apropos("chi")
```
### 正規乱数
正規乱数はもう一つの基本。Random numbers in normal distribution are important, too.
```{r}
v3 <- rnorm(n)
plot(v3)
sorted.v3 <- sort(v3)
inv.sorted.v3 <- sort(v3,decreasing =TRUE)
plot(sorted.v3)
plot(inv.sorted.v3)
# 分布のプロット Distribution plot
hist(v3)
plot(density(v3))
boxplot(v3)
# 記述統計する Descriptive statistics.
summary(v3)
```

### データフレーム Data frames : 複数のベクトルの束 Bundle of multiple vectors

2つのベクトルを併せる。Combine two vectors.
```{r}
n <- 20
v4 <- runif(n)
v5 <- sample(c("male","female"), n, replace = TRUE)
df1 <- data.frame(value=v4, sex=v5)
df1
```
プロットしてみる。 Plot the data.frame.
```{r}
plot(df1)
```

追加する。Add vecors.
```{r}
v6 <- rnorm(n)
v7 <- sample(1:4, n, replace = TRUE)
df2 <- data.frame(df1, x = v6, y = v7)
df2
plot(df2)
boxplot(df2)
summary(df2)
```

## 行列 Matrix

すべての項目が数値を持つなら行列の方が便利かもしれない。線形代数処理など。プロットも違うタイプのものがある。When all values are numeric, matrix rather than data frame might be useful. For example, application of linear algebra to it. Different kind of plots are also available.

```{r}
# sex "male","female" are changed into 0,1
v5.num <- sample(0:1,n,replace=TRUE)
m1 <- cbind(v4,v5.num,v6,v7)
m1
plot(as.data.frame(m1))
image(m1)
matplot(m1,type="l")
persp(m1)
```

行列には行ごとの操作と列ごとの操作をすることができる。You can apply a function to each rows or each columns of matrix.
```{r}
m2 <- matrix(c(1,2,3,4,5,6),2,3)
m2
apply(m2,1,sum)
apply(m2,2,max)
```

## リスト List

長さが揃っていないものとまとめる。
When you want to put multiple things that are different in their length, use list.

```{r}
x <- 1:3
y <- 1:4
L <- list(x,y)
L
L[[1]]
L[[2]]
L$x
L$y
L1 <- list(x=x,y=y)
L1[[1]]
L1[[2]]
L1$x
L1$y
```
## テキストファイルの読み書き Read/Write text files

いくつかのオプションがありますが、一つを覚えて使いまわします。
There are multple ways to do but start using one way and in case you really need expand your skills on this.

ファイルのハンドリングをするためには、コンピュータ上でのRの作業場所を確認したり変更したりする必要があります。

When you handle text files on your computer, you should check where you use R on your computer and be able to move to where you want to go.

```{r}
# working directory : wd
apropos("wd")
getwd()
help(setwd)
```

RのGUIを使っているならツールバーからディレクトリ変更をすることができるはず。
If you are using R's RUI, you can change your working directory though its toolbars menu.

さて、ファイルのハンドリング。
Let's start handling your textfiles.

エクセル等のスプレッドシートアプリケーションを使って、数行、数列の数値ばかりのシートを作り、それをタブ区切りテキストファイルで保存時ます。たとえば"mydata.txt"という名前で。

Open excel or similar spread-sheet application and fill a rectangle area with several rows x several columns.
Save the sheet as a tab-delimitated text file, for example with name of "mydata.txt".

```{r}
mydata <- read.table("mydata.txt")
mydata
```

すでに作った行列m1を"m1.txt"という名前で保存してみます。
Let's save the matrix, m1, as a text file "m1.txt".

```{r}
write.table(m1,"m1.txt")
```

## 練習問題 Exercises
* 以下の問題のコマンドをテキストファイルで保存して提出せよ。Save your commands for the following exercises and report the file.
 * Q x軸とy軸のラベルと図のタイトルとを適当な語にして散布図を描け。Draw a coplot with its x-axis and y-axis labels as well as title being what you like.
 * 乱数は分布に従って発生させる。一様分布乱数、正規分布乱数、ポアソン分布乱数、指数分布乱数、ベータ分布乱数のヘルプ記事を表示し、そのExamplesを実行せよ。When you generate random numbers, you need distributions. Open help articles on uniform, normal, poisson, exponential and beta distributions and run the example codes.
 * sample()関数のオプションreplace=TRUEとreplace=FALSEの違いをヘルプ記事で確認し、次のコードが実行することを、言葉で説明せよ。エラーメッセージの異同についても確認すること。そしてその異同の理由についても説明すること。Check help article of sample() for the option relace = TRUE or FALSE and describe what the following codes do. Pay attention to the error messages and check the presence of difference and describe the difference/lack of difference.
 ```{r}
 sample(1:10)
 sample(1:10,8)
 sample(1:10,8,replace=TRUE)
 sample(1:10,8,replace=FALSE)
 sample(1:10,12)
 sample(1:10,12,replace=TRUE)
 sample(1:10,12,replace=FALSE)
 ```
 * 3列のデータフレームを適当に作成せよ。Make 3 column-data frame as you like.
 * 7行5列の整数を要素とする行列を作り、各列の平均要素とするベクトルを作れ。また各行の最小値を要素とするベクトルを作れ。Generate a 7-row 5-column-matrix (7x5 matrix). Make a vector whose elements are mean of each columns of the matrix. Make a vector whose elements are minimum value of each rows of the matrix.
 
 # 関数を作る Make functions

## 関数と入出力 Functions and their inputs and outputs

関数は仕事をします。
Functions do jobs.

関数は
* 名前を持っていて
* 入力を受け取り
* 出力を返します

Functions
* have their name
* take inputs
* returns outputs.

```{r}
# "my1stFx" is "name of function" you make
# x is input(s)
# y is output(s)
my1stFx <- function(x){
  y <- 2*x
  return(y)
}
# 使ってみます。Usage
x <- 1:3
x
my1stFx(x)
# 出力は格納できます You can put the output into an object
my1stOut <- my1stFx(x)
my1stOut
```

入出力は複数にもできる。
Inputs and outputs can be multiple.

```{r}
my2ndFx <- function(x,a){
  y <- x+1
  b <- a + rnorm(length(a))
  return(list(xAdd1=y,aAddz=b))
}
x <- 1:10
a <- x*2+3
my2ndFx(x,a)
my2ndOut <- my2ndFx(x,a)
my2ndOut
my2ndOut$xAdd1
my2ndOut$aAddz
plot(x,a)
plot(my2ndOut$xAdd1,my2ndOut$aAddz,col=2)
```

## if条件分岐とループ If conditions and Looping

### if
```{r}
# 多すぎるデータの場合にだけ一部を取り出す
# When number of records is too many, sample the elements.
my3rdFx <- function(x,n){
  if(length(x)>n){
    y <- sample(x,n)
  }else{
    y <- x
  }
  return(y)
}

x <- rnorm(50)
n <- 100
my3rdOut <- my3rdFx(x,n)
length(my3rdOut)

x.2 <- rnorm(110)
my3rdOut.2 <- my3rdFx(x.2,n)
length(my3rdOut.2)
```

### ループ Looping

A,T,G,Cの塩基配列をn.iter個作る。
長さは20塩基前後とし、A,T,G,Cは等確率でランダムに選ぶとする。

Let's make n.iter base-sequences with A,T,G,and C.
Their length should be around 20 and bases should be chosen with eaqual probability at random.
```{r}
n.iter <- 5
mean.len <- 20
Seqs <- list()
ATGC <- c("A","T","G","C")
pr <- rep(1/4,4)
for(i in 1:n.iter){
  len <- rpois(1,mean.len)
  Seqs[[i]] <- sample(ATGC,len,replace=TRUE,prob=pr)
}
Seqs
```

ランダムな塩基配列を数を増やして作成し、どんな長さ分布かを見てみることとする。
Let's generate much more sequences and check the distribution of length of generated sequences.

```{r}
n.iter <- 1000
mean.len <- 20
Seqs <- list()
ATGC <- c("A","T","G","C")
pr <- rep(1/4,4)
for(i in 1:n.iter){
  len <- rpois(1,mean.len)
  Seqs[[i]] <- sample(ATGC,len,replace=TRUE,prob=pr)
}

len.check <- lapply(Seqs,length)
hist(unlist(len.check))
len.check.s <- sapply(Seqs,length)
hist(len.check.s)
mean(unlist(len.check))
mean(len.check.s)
```

### if とループの組み合わせ Combine if and looping

上記でやったように適当な配列を作りつつ、読み取りエラーを想定に加える。
うまく読めなかった塩基を"N"として加え、"N"の入った配列は捨て、その上で配列を指定の本数だけ作成する。

As we did above, we generate random sequences but this time, bases might not be read well. "N" represents the unread bases and sequences that has one or more "N" should be discarded. We want to have n.iter sequences without "N".

まずは"N"で排除しない場合。
"N" is allowed at first.
```{r}
n.iter <- 10
mean.len <- 20
Seqs <- list()
ATGCN <- c("A","T","G","C","N")
pr <- c(rep(1/4,4),0.05)
pr <- pr/sum(pr)
pr
for(i in 1:n.iter){
  len <- rpois(1,mean.len)
  Seqs[[i]] <- sample(ATGCN,len,replace=TRUE,prob=pr)
}
Seqs
```

"N"があったら登録しない。
When "N" is included, a vacant sequence is registered.

```{r}
n.iter <- 10
mean.len <- 20
Seqs <- list()
ATGCN <- c("A","T","G","C","N")
pr <- c(rep(1/4,4),0.05)
pr <- pr/sum(pr)
pr
for(i in 1:n.iter){
  len <- rpois(1,mean.len)
  tmp.seq <- sample(ATGCN,len,replace=TRUE,prob=pr)
  checker <- which(tmp.seq=="N")
  if(length(checker)>0){
    Seqs[[i]] <- ""
  }else{
    Seqs[[i]] <- tmp.seq
  }
}
Seqs
```

"N"が含まれているときには登録自体を飛ばすことにする。
When "N" is included, no registration.

```{r}
n.iter <- 10
mean.len <- 20
Seqs <- list()
ATGCN <- c("A","T","G","C","N")
pr <- c(rep(1/4,4),0.05)
pr <- pr/sum(pr)
pr
# カウンタを導入する
cnt <- 1
for(i in 1:n.iter){
  len <- rpois(1,mean.len)
  tmp.seq <- sample(ATGCN,len,replace=TRUE,prob=pr)
  checker <- which(tmp.seq=="N")
  if(length(checker)>0){
    #Seqs[[i]] <- ""
  }else{
    Seqs[[cnt]] <- tmp.seq
    cnt <- cnt + 1
  }
}
Seqs
```

指定の数の配列を作る。
Let's generate given-number of sequences.

```{r}
n.iter <- 10
mean.len <- 20
Seqs <- list()
ATGCN <- c("A","T","G","C","N")
pr <- c(rep(1/4,4),0.05)
pr <- pr/sum(pr)
pr
# カウンタを導入する
cnt <- 1
# ループを続けるかどうかの判断のため
# Continue/discontinue looping?
continueLoop <- TRUE
while(continueLoop){
  len <- rpois(1,mean.len)
  tmp.seq <- sample(ATGCN,len,replace=TRUE,prob=pr)
  checker <- which(tmp.seq=="N")
  if(length(checker)>0){
    #Seqs[[i]] <- ""
  }else{
    Seqs[[cnt]] <- tmp.seq
    cnt <- cnt + 1
  }
  # 登録本数が十分になったらループをやめる
  # When adequate sequences are registered, stop looping
  if(cnt==n.iter+1){
    continueLoop <- FALSE
  }
}
Seqs
```

# Rを用いた医科学のためのデータシミュレーションの基礎2 Basics of Data Simulation for Medical Sciene with R 2

## 確率密度分布と累積分布、Probability Density Distribution and Cumulative Density Distribution

値が分布をなしているとき、大きく分けて2つのとらえ方がある。
その一つが密度分布、もう一つが累積分布である。

When values are in a disrtibution, there are two ways to grab them. One is a density distribution and the other is a cumulative distribution.

```{r}
# 正規分布に従う乱数を発生させる
# Generate random numbers in a normal distribution.
n <- 100
x <- rnorm(n)
boxplot(x)
```
```{r}
hist(x)
```
```{r}
plot(density(x))
```
```{r}
plot(sort(x))
```

ソートしてプロットすると累積に関する情報が得られる。
いくつかのやり方がある。
Plot a sorted value set, then you get an idea on their cumulative distribution.
There are several ways for it. 

```{r}
n <- 10
x <- rnorm(n)
s.x <- sort(x)
plot(s.x)
```

x-y 反転。
Exchange x and y axes.
```{r}
s <- seq(from=0,to=1,length=length(x)+1)
s <- s[-1]
plot(s.x,s)
```

階段状に。
Like steps
```{r}
plot(s.x,s,type="s")
```

打点の高さを少しいじる。
Modify the y-axis values a bit.

便利な関数もある。
You may find a useful function.

```{r}
plot(s.x,ppoints(n,a=0.5),type="b")
```
```{r}
plot.stepfun(s.x)
points(s.x,s,col=2)
s.2 <- seq(from=0,to=1,length=length(x)+1)
s.2 <- s.2[-length(s.2)]
points(s.x,s.2,col=3)
points(s.x,ppoints(n,a=0.5),col=4,type="b")
```

## 色々な分布 Various Distributions

### 二項分布 Binomial distribution

治療成功率がpであるとき、n人を治療するとn人のうちk=0,1,...,n人で治療が成功する確率は以下のように計算できる。

Assume a sequence of n independent treatments of patients and the success rate of treatmen is p.
The number of successes, k, in a n treatments ranges from 0 to n and the probability of the number being k is given as below.


$$
\begin{equation}
\begin{pmatrix}
n\\k
\end{pmatrix} p^k (1-p)^{n-k} = \frac{n!}{k!(n-k)!}p^k (1-p)^{n-k}
\end{equation}
$$

```{r}
n <- 5
# 整数列を作る Integer sequence generation
ks <- 0:n
# 二項分布の確率 Probability of binomial distribution
p <- 0.7
pr.binom <-dbinom(ks,n,p)
# 全事象の確率の和は1. Sum of probability of all events should be 1.
sum(pr.binom)
plot(ks,pr.binom,type="h")
```

今、n人の治療をして、k人で成功したとき、成功確率pの推定をしたいとする。
Assume you observed k successes in n treatments and you want to estimate the success probability p.

pの値はわからないので、取りうる値をたくさん仮定してしらみつぶしにn人中k回の成功がおきる確率を計算してみる。

Because we have no idea on the value of p, let's calculate probability to observe k successes in n experiments with many values for p.

pは0から1の範囲にあるから、そこに均等に値を取ることにする。１次元格子を作るとも言う。

p ranges from 0 to 1. Therefore we generate multiple p values with fixed intervals. This is generation of a lttice in 1 dimensional space.

```{r}
n.pt <- 100
ps <- seq(from=0, to=1, length=n.pt)
n <- 5
k<-2
probs <- dbinom(k,n,ps)
plot(ps,probs,type="l")
```

pの値を均等割りにする代わりに、たくさんの乱数を使うこともできる。
Instead of generating p values with fixed intervals, we can use random numbers.

一様分布からの乱数を用いる。
Random numbers in uniform distribution are used.

```{r}
n.pt <- 100
ps.2 <- runif(n.pt)
n <- 5
k<-2
probs.2 <- dbinom(k,n,ps.2)
plot(ps.2,probs.2,type="l",col=2)
```

作成した乱数列をソートすればきちんと描ける。
Sort random numbers first, then you will get the appropriate curve.

```{r}
ps.2.sorted <- sort(ps.2)
probs.2.sorted <- dbinom(k,n,ps.2.sorted)
plot(ps.2.sorted,probs.2.sorted,type="l",col=2)
```

実は、この分布は、ベータ分布と同じ形をしている。
Actually this distribution has the same shape with beta distribution.

```{r}
# kとn-kに1を加えているのがポイント
# Add 1 to k and n-k
beta.prob <- dbeta(ps,k+1,n-k+1)
plot(ps,beta.prob,type="l")
```

Exercise.

* n=50人を治療し、k=42人の治療が成功したという。この治療法の成功率pがどのくらいなのかを、dbeta()関数を使う方法とそうでない方法とでグラフとして示せ。Assume you treated n=50 patients and k=42 patients were successfully treated. Draw a probability graph for success rate p with two different ways; using dbeta() function and not using the function.

## 二項分布で分割表 Contingency tables with binomial distribution

2つの治療法A,Bがあって、それぞれの成功率はpa,pbであるとする。全部でn人を2群にランダムに割り付け治療するとする。

Assume 2 treatments, A and B and their success rates are pa and pb, respectively. Now we treat n patients with random assignment. 

```{r}
n <- 200
pa <- 0.6
pb <- 0.65
na <- rbinom(1,n,0.5)
nb <- n-na
a.success <- rbinom(1,na,pa)
a.failure <- na-a.success
b.success <- rbinom(1,nb,pb)
b.failure <- nb-b.success
fisher.test(matrix(c(a.success,a.failure,b.success,b.failure),2,2))
```

これを何度も繰り返してみる。

```{r}
n <- 200
pa <- 0.6
pb <- 0.65
n.iter <- 500
pvals <- nas <- nbs <- rep(0,n.iter)
for(i in 1:n.iter){
  nas[i] <- rbinom(1,n,0.5)
  nbs[i] <- n-nas[i]
  a.success <- rbinom(1,nas[i],pa)
  a.failure <- nas[i]-a.success
  b.success <- rbinom(1,nbs[i],pb)
  b.failure <- nbs[i]-b.success
  tmp.out <- fisher.test(matrix(c(a.success,a.failure,b.success,b.failure),2,2))
  pvals[i] <- tmp.out$p.value
}
hist(pvals,main="p.values")
```
```{r}
plot.stepfun(sort(pvals),main="p.values cumulative")
```
```{r}
hist(nas,main="nas")
```

2つの治療法A,Bがあって、それぞれの成功率はpa,pbであるとする。全部でn人を2群にランダムに割り付けながら治療するとする。人数を増やしながら逐一検定するとどうなるかをシミュレーションしてみる。

Assume 2 treatments, A and B and their success rates are pa and pb, respectively. Now we treat n patients with random assignment. Along with the assignment, let's test the difference between A and B.

```{r}
n <- 500
pa <- 0.6
pb <- 0.7
pvals <- rep(0,n)
assignment <- sample(0:1,n,replace=TRUE)
tbl <- matrix(0,2,2)
for(i in 1:n){
  if(assignment[i] == 0){#A
    tmp.success <- sample(0:1,1,prob=c(1-pa,pa))
    if(tmp.success==0){# failure
      tbl[1,1] <- tbl[1,1]+1
    }else{# success
      tbl[1,2] <- tbl[1,2]+1
    }
  }else{#B
    tmp.success <- sample(0:1,1,prob=c(1-pb,pb))
    if(tmp.success==0){# failure
      tbl[2,1] <- tbl[2,1]+1
    }else{# success
      tbl[2,2] <- tbl[2,2]+1
    }
  }
  tmp.test.out <- fisher.test(tbl)
  pvals[i] <- tmp.test.out$p.value
}
plot(pvals,type="l")
```

この推移を1度の治験だとすれば、複数回の治験をすればそれごとのp値の推移がシミュレーションできる。

This broken line is on one clinical trial. Multiple clinical trials will produce multiple lines.

```{r}
n.iter <- 10
n <- 500
pa <- 0.6
pb <- 0.7
#pvals <- rep(0,n)
pvals <- matrix(0,n.iter,n)
for(j in 1:n.iter){
  assignment <- sample(0:1,n,replace=TRUE)
  tbl <- matrix(0,2,2)
  for(i in 1:n){
    if(assignment[i] == 0){#A
      tmp.success <- sample(0:1,1,prob=c(1-pa,pa))
      if(tmp.success==0){# failure
        tbl[1,1] <- tbl[1,1]+1
     }else{# success
       tbl[1,2] <- tbl[1,2]+1
      }
    }else{#B
     tmp.success <- sample(0:1,1,prob=c(1-pb,pb))
     if(tmp.success==0){# failure
        tbl[2,1] <- tbl[2,1]+1
      }else{# success
        tbl[2,2] <- tbl[2,2]+1
     }
    }
    tmp.test.out <- fisher.test(tbl)
    pvals[j,i] <- tmp.test.out$p.value
  }
  
}
matplot(t(pvals),type="l")
```

## ポアソン分布 Poisson distribution

ある確率でランダムに起きる事象がある。
一定期間内に、または一定の区間内にその事象が起きる回数についての分布。
その期間・区間中に何回起きたかは0,1,2,...$\infty$回になるけれど、その確率の分布。

During a fixed period, or in a segment with fixed length, a kind of events occurs at random. The number of events that really occur in a period/segment ranges from 0 to $\infty$. Poisson distribution tells the probability that the events occurs k times in a period/segment. 

今、1人の女性が1生に産む子どもの数の平均がk人だとする。
子どもの数が0,1,2,...人であって、これがポアソン分布に従うと仮定する。

Assume the mean number of offsprings that one female has per her life is k. The number of offsprings can be 0,1,2,... and assume the probability is in Poisson distribution.

```{r}
k <- 1.6
ks <- 0:20
pr.kids <- dpois(ks,k)
sum(pr.kids)
# 平均子ども数 mean of number of children
sum(pr.kids * ks)
k
plot(ks,pr.kids,type="h")
```

この調査を子どものいる女性だけで調査するとどうなるだろうか。

Assume you survey the number of children by asking women with any children.

```{r}
k <- 1.6
ks <- 0:20
pr.kids.2 <- dpois(ks,k)
# 子どものいない女性は除かれる Women without any child is out of count
pr.kids.2[1] <- 0
sum(pr.kids.2)
# 平均子ども数 mean of number of children
# カウントされた女性の相対頻度で計算する
# The relative fraction of women with children should be used for the calculation
pr.kids.2.rel <- pr.kids.2/sum(pr.kids.2)
sum(pr.kids.2.rel * ks)
k
```


## 次世代シークエンシングとポアソン分布 Next-generation sequencing (NGS) and Poisson distribution

次世代シークエンシングにおいて、ゲノムの各箇所が平均nデプスで読まれるとき、デプス別の分布はポアソン分布に近似できる。

NGS returns various loci with various depth. The distibution of number of loci per depths is close to Poisson distribution.

この仮定の下で計算すると、どのくらいの割合で「全く読まれない」のだろうか？

How often do we miss loci completely?

```{r}
mean.depth <- 10
ds <- 0:50
pr.pois <- dpois(ds,mean.depth)
sum(pr.pois)
pr.pois[1]
```

今、NGSを用いてSNVの箇所の検出をしているとする。

Assume we are to detect SNVs with NGS.

SNVでは理想的には、アレルごとのDNA分子数が等量である。
Number of DNA molecules per alleles are identical, theoretically.

アレルごとの平均デプスはローカスことの平均デプスの半分である。

Mean depth per alleles is the half of mean depth per loci.

2アレルM,mごとのデプス別確率を求める。

Probability per depths for 2 alleles, M and m.

処理内容を理解しやすくするため、mean.depthを小さくしてやってみる。

For easier understanding the procedures, let's make mean.depth smaller.

```{r}
mean.depth <- 6
mean.depth.allele <- mean.depth/2
ds <- 0:20
pr.pois.M <- dpois(ds,mean.depth.allele)
pr.pois.m <- pr.pois.M
pr.M.m <- outer(pr.pois.M,pr.pois.m,"*")
image(ds,ds,pr.M.m)
```

2アレルの本数の合計がn.all 以上の場合にデータとして採用し、マイナーアレルのリード数がn.minor以上の場合にSNVであるとみなすとすると、どのくらいの確率でSNVであると認められるだろうか。

When you accept your data only when sum of two alleles is equal or more than n.all and when the number of reads of minor allele is equal or more than n.minor, you judge they are SNV. Then you can calculate the probability to identify the loci SNV.

```{r}
n.both <- outer(ds,ds,"+")
n.both
n.all <- 5
n.minor <- 2
ok.n.all <- which(n.both>=n.all,arr.ind=TRUE)

# ok.minor
pr.pois.M.ok <- pr.pois.m.ok <- pr.pois.M
pr.pois.M.ok[1:n.minor] <- pr.pois.m.ok[1:n.minor] <- 0
pr.M.m.ok <- outer(pr.pois.M.ok,pr.pois.m.ok,"*")
image(ds,ds,pr.M.m.ok)
```

n.minorの基準を満たさない図の左下隅が除かれた。

Left-bottom cells were removed due to n.minor restriction.

```{r}
# n.all okをかぶせる
pr.M.m.ok.all.ok <- matrix(0,length(ds),length(ds))
pr.M.m.ok.all.ok[ok.n.all] <- pr.M.m.ok[ok.n.all]
image(ds,ds,pr.M.m.ok.all.ok)
```

n.all基準により右下がりの対角線より上側のみが採用された。
The upper side of the declining line was accepted due to n.all restriction.

```{r}
# n.allもn.minorも基準を満たした割合
# The fraction with both n.all and n.minor adequate.
sum(pr.M.m.ok.all.ok)
```

## Exercise
* 平均デプスと、アレル合算デプス基準値、アレル別デプス基準値を適当に与え、SNV検出感度を算出せよ。 Set mean depth value and n.all cutoff and n.minor cutoff as you like and calculate the sensitivity to detect SNV.

* 数えられる病変数はポアソン分布で表すことができるかもしれない。たとえば、身体のあちこちに現れる皮膚病変、ランダムに起きているように思われる神経変性部位、相当するうある関節のうちの罹患関節数。何かしら身近な例を取ってシミュレーションしてみよ。 Number of countable lesions can take Poisson distributions, e.g., skin lesions scattered out on the body surface, neural degenerative lesions appearing at random, affected joints among substantial number of them. Simulate a medical/biological phenomenon that may show Poisson distribution.


## 正規分布 Normal distribution

血清コレステロール値が群Aでは平均ma,分散vaの、群Bでは、それぞれmb,vbの正規分布に従っているという。

Assume serum cholesterol concentration in group A is in a normal distribution with mean being ma and variance being va and in group B, mean is mb and variance is vb.

```{r}
na <- 150
nb <- 120
ma <- 200
va <- 100
mb <- 180
vb <- 80
chol.a <- rnorm(na,ma,sqrt(va))
chol.b <- rnorm(nb,mb,sqrt(vb))
h <- hist(c(chol.a,chol.b))
hist(chol.a,density=20,col=2,add=TRUE,breaks=h$breaks)
hist(chol.b,density=17,col=3,add=TRUE,breaks=h$breaks)
```

t検定をしてみる。
Let's t-test these data.

```{r}
t.test(chol.a,chol.b)
```

もし、群Aと群Bとが血清コレステロール値に関して純粋な分類であれば、値のまとまりが良い(分散が小さい,vx)が、さまざまな関与因子があってその影響が加わっているとすると、次のようにばらつきを加えることになる。

```{r}
vx <- 100
chol.a.1 <- chol.a + rnorm(na,0,sqrt(vx))
chol.b.1 <- chol.b + rnorm(nb,0,sqrt(vx))
t.test(chol.a.1,chol.b.1)
h <- hist(c(chol.a.1,chol.b.1))
hist(chol.a.1,density=20,col=2,add=TRUE,breaks=h$breaks)
hist(chol.b.1,density=17,col=3,add=TRUE,breaks=h$breaks)
```

ここでは、群Aに関して、まず分散vaの正規分布で値を生成し、そこにさらなるばらつきの項を分散vxの正規乱数で加えた。

We generated random values in normal distribution with variance va then added normal random values of variance vx.

正規分布では、分散を足すことができるので、以下のようにもできる。

You can add variances for normal distribution.

2つの方法で生成した分布を比べてみる。
両者はほとんど同じである。

Let's compare two distributions.
Both two are almost identical.

```{r}
chol.a.2 <- rnorm(na,ma,sqrt(va+vx))
summary(chol.a.1)
summary(chol.a.2)
h <- hist(c(chol.a.1,chol.a.2))
hist(chol.a.1,density=20,col=2,add=TRUE,breaks=h$breaks)
hist(chol.a.2,density=17,col=3,add=TRUE,breaks=h$breaks)
```

2つの分布が近いことはソートして散布図を描くことで示すこともできる。

Coplot of sorted values can show the similarity of two distributions.

```{r}
plot(sort(chol.a.1),sort(chol.a.2))
abline(0,1,col=2)
```

では、2群の平均は固定し、2群の分散をだんだん大きくして、そのt検定の結果を比べてみる。

Now fix two groups' mean and change their variance and see their effect on the result of t-test.

```{r}
na <- 100
nb <- 100
ma <- 200
mb <- 180
vabs <- 2^(seq(from=0,to=10,length=100))
ps <- rep(0,length(vabs))
for(i in 1:length(vabs)){
  chol.a <- rnorm(na,ma,sqrt(vabs[i]))
  chol.b <- rnorm(nb,mb,sqrt(vabs[i]))
  tmp.out <- t.test(chol.a,chol.b)
  ps[i] <- tmp.out$p.value
}
plot(log(vabs),log(ps))
```

上の例では、2つの群を決める何かがあり、それが群の平均値を決めるというモデルであった。
なぜ正規分布を取るかというと、群(とその平均値)を決める以外の要因が正規分布を作り出しているというモデルであった。

We assumed two groups that have mean value per groups and other factors added normal variation from the means in the above.

今度は、複数の要因がそれぞれ正規分布様の寄与をしているときに、そのうちの1要因について効果のある介入をするというようなシミュレーションをしてみる。

Our next example is as follows. Multiple factors are related to a phenotype and each of them contributes in a normal distribution-like fashion. Toward the phenomenon, an intervention specific to one of the factors is given.

介入は、第1因子をk倍(k<1)に小さくする働きをもち、その効果にはばらつきがあるものとする。

Assume the intervention make the 1st factor x k (k<1) with random variation.


```{r}
n.factor <- 4
ms <- runif(n.factor)
vs <- runif(n.factor)^2

n <- 100
x <- matrix(0,n,n.factor)
for(i in 1:n.factor){
  x[,i] <- rnorm(n,ms[i],sqrt(vs[i]))
}
X <- apply(x,1,sum)
hist(X)
```
```{r}
h <- hist(c(x))
for(i in 1:n.factor){
  hist(x[,i],density = 10 +i*3,add=TRUE, breaks=h$breaks,col=i+1)
}
```
```{r}
# k is interventional effect
k <- 0.4
# v.tr is variance of treatment effect
v.tr <- 1
x1.post <- x[,1] * k + rnorm(n,0,sqrt(v.tr))
x.post <- x
x.post[,1] <- x1.post
X.post <- apply(x.post,1,sum)
h <- hist(c(X,X.post),density =20)
hist(X,density=17,add=TRUE,col=2,breaks=h$breaks)
hist(X.post,density=17,add=TRUE,col=3,breaks=h$breaks)
```

治療効果が表れて、介入後の分布が左にシフトしている。

Intervention is successful and the distribution shifted smaller.

治療前と治療後との分布の違いを検定することはできる。

You can test the difference between the pre- and post-interventional distributions.

```{r}
t.test(X,X.post)
```

しかしながら、取られたデータは以下のように、治療前後の対応関係を持っている。

However the pre- and post-data are paired.

```{r}
xlim <- range(c(X,X.post))
plot(X,X.post,xlim=xlim,ylim=xlim)
```

従って、その対応を考慮して効果を判定する方がよい。

```{r}
t.test(X,X.post,paired=TRUE)
```

p値がぐっと小さくなることがわかる。

## Exercise
* この例で、次の4変数について、個別にその値を振り、その他の変数の値は固定し、値を振った変数がpaired検定に及ぼす影響についてシミュレーション的に評価せよ
 * 第1因子の寄与ms[1]
 * 第1因子のばらつきvs[1]
 * 介入の効果k
 * 介入の効果のばらつきv.tr

# Rを用いた医科学のためのデータシミュレーションの基礎3 Basics of Data Simulation for Medical Sciene with R 3

## ガンマ分布と生存期間 Gamma distribution and survival

次のグラフの形を見てみよう。

Look at the graph below.

```{r}
k <- 80
t <- 1
a <- seq(from=0,to=120,by=1)
cum.d <- pgamma(a,k,t)
plot(a,cum.d,type="l",xlab="age",ylab="positive fraction")
```

横軸の60歳まではほとんど何も起こらず、それ以降、急激に多くの割合が該当するようになり、100歳を超えたあたりでほぼ全員が該当する、というカーブである。

The horizontal axis is age.
Upto 60 years, almost zero meaning nothing happened, then suddenly substantial fraction of people is counted "positive" and over age 100, almost all of them turn to be "positive".

加齢によって発病する疾患のモデルや死亡のモデルとして使うことができるカーブであるが、これは、ガンマ分布の累積分布である。

This curve can be read as a positive fraction of a late-onset disease or death. This curve is a cumulative distribution of gamma distribution.

これは累積分布である。対応する確率密度分布も描ける。

This is the cumulative distribution. You can draw its density distributions.

```{r}
d <- dgamma(a,k,t)
plot(a,d,type="l",xlab="age",ylab="probability")
```

80歳をピークにした発症分布である。

The onset is peaked at the age 80.

分布の関数があるので、乱数も発生できるから、それを使ってシミュレーションをする。

Besides the functions used above, there is a function to generate random numbers in gamma distribution. Now we are to simulate phenomena following gamma distribution.

ある因子を持つか持たないかで、平均発症時期がずれると仮定しよう。

Assume posession of a factor affects on the mean age of onset.

```{r}
n1 <- 100
n2 <- 100
m1 <- 80
m2 <- 70
t <- 1
x1 <- rgamma(n1,m1,t)
x2 <- rgamma(n2,m2,t)
h <- hist(c(x1,x2))
hist(x1,density=20,col=2,add=TRUE)
hist(x2,density=17,col=3,add=TRUE)
```

ガンマ分布は生存解析データシミュレーションにも使える。

Gamma distribution is useful also for data simulation of survival analyses.

今、ある疾患にかかると、平均余命mか月で亡くなるとする。誰もが同じ確率で亡くなるとすると、次のようにシミュレーションできる。

Assume patients who develop a disease are to die at random with the mean months to live being m. This can be simulated as below.

```{r}
m <- 8
k <- 1
t <- k/m
months <- seq(from=0,to=24,length=100)
cum.d <- pgamma(months,k,t)
plot(months,cum.d,type="l")
```

生存に着目すれば、

Survival can be plotted,

```{r}
plot(months,1-cum.d,type="l",xlab="months",ylab="fraction of survivals")
```

今、治療法A,Bがあり、それぞれ平均余命をma,mbか月にするも、カーブの形は変えないとします。
それぞれでna,nb人が治療されたとすれば、それぞれの余命は以下のようにシミュレーションできます。

Assume there are two ways to treat, A and B and their mean months to live are ma and mb, respectively and the shape of curves don't change.
When na and nb patients are treated with A and B, this can be simulated as below.

```{r}
na <- 100
nb <- 80
ma <- 8
mb <- 12
k <- 1
ta <- k/ma
tb <- k/mb
x.a <- rgamma(na,k,ta)
x.b <- rgamma(nb,k,tb)
mean(x.a)
mean(x.b)
library(survival)
censor.a <- rep(1,na)
censor.b <- rep(1,nb)

KMcurve.a <- survfit(Surv(x.a, censor.a)~1)
KMcurve.b <- survfit(Surv(x.b, censor.b)~1)
plot(KMcurve.a, conf.int=TRUE, mark.time=TRUE)
```

```{r}
KMcurve.ab <- survfit(Surv(c(x.a,x.b),c(censor.a,censor.b))~c(rep(1,na),rep(2,nb)))
plot(KMcurve.ab, conf.int=FALSE, mark.time=TRUE,col=1:2)
```

この例では、治療は平均余命にのみ影響を与え、カーブの形には影響を与えないモデルであった。

いずれの治療も初期には死亡を防げるが、先延ばしにした分、加速して死亡が観察されるようなこともあるだろう。
それをシミュレーションしてみる。

In this simulation, A and B only affected the mean months to live but no effect on the shape of curve.
In some cases, deaths in early phases are avoided but the accumulation of death events may pile up later.
Next example is for it.

それをするために、ma,mbは変えずにガンマ分布のパラメタkを変える。

For this purpose, the parameter k for gamma distribution is modified.

```{r}
na <- 100
nb <- 80
ma <- 8
mb <- 12
ka <- 7
kb <- 2
ta <- ka/ma
tb <- kb/mb
x.a <- rgamma(na,ka,ta)
x.b <- rgamma(nb,kb,tb)
mean(x.a)
mean(x.b)
library(survival)
censor.a <- rep(1,na)
censor.b <- rep(1,nb)

KMcurve.a <- survfit(Surv(x.a, censor.a)~1)
KMcurve.b <- survfit(Surv(x.b, censor.b)~1)
plot(KMcurve.a, conf.int=TRUE, mark.time=TRUE)
```
```{r}
plot(KMcurve.b, conf.int=TRUE, mark.time=TRUE,col=2)
```
```{r}
KMcurve.ab <- survfit(Surv(c(x.a,x.b),c(censor.a,censor.b))~c(rep(1,na),rep(2,nb)))
plot(KMcurve.ab, conf.int=FALSE, mark.time=TRUE,col=1:2)
```

## Exercises
* 自由な設定で2群の生存解析データを作成しプロットせよ。Simulate survival data set of two groups as you like.