- こちらの続き
- 医学・生物学データの扱いに慣れるための、データシミュレーションセミナーの資料(医学研究科大学院コース ゲノム遺伝学＠京大 2014:こちら)
- htmlファイルはこちら(Login as a guest でログインすると見られます)
- こちらを参考にすればhtml化可能です
次のグラフの形を見てみよう。
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)
```
* 自由な設定で2群の生存解析データを作成しプロットせよ。Simulate survival data set of two groups as you like.
