Following the above steps, we generate survival data with exponential distribution.
event.times <- rexp(100, 0.9)
event.times
## [1] 0.031468670 1.699722626 0.012604250 0.181851015 0.325474574 0.600348067
## [7] 0.253768509 0.106565014 0.239200077 4.722760075 0.045875484 0.307674001
## [13] 0.420839034 0.286222322 0.084884797 0.043504038 0.078568163 1.456698330
## [19] 4.382265779 1.383426225 1.634379621 1.224719648 1.102386634 1.590425714
## [25] 0.332396754 1.440037873 0.141629408 0.097233212 3.595864890 0.540339829
## [31] 0.627334975 0.505893539 0.296015002 0.274083299 2.366269630 2.113700263
## [37] 1.241166406 0.116164153 0.215933847 2.883401751 2.758792475 0.007852311
## [43] 3.732606317 0.118259397 1.112956394 1.981260307 0.219232355 1.906241483
## [49] 0.565370676 1.077138985 0.484502960 0.903881538 2.263368155 0.853008879
## [55] 0.825220779 1.633502468 1.021393551 4.096878624 0.129284012 0.507789280
## [61] 0.065935818 5.602531021 0.183929479 0.886052206 0.604041278 2.703122190
## [67] 0.044561327 0.581029905 0.162785934 1.398531164 2.622305360 2.177511017
## [73] 0.053606869 0.302974220 0.187704037 1.078748603 0.710940666 0.271044257
## [79] 1.606070770 0.751399968 0.480695701 3.406351614 0.198755058 0.237405763
## [85] 0.121899657 1.054376699 1.822300313 2.440518850 0.198321935 3.507952251
## [91] 0.516786499 1.255904736 0.051519817 0.808046473 0.039023547 0.119380507
## [97] 0.182533470 0.124188870 0.227555450 0.926034765
2.Suppose the constant cause-specific hazards \(\alpha_{01} = 0.3\) and \(\alpha_{02} = 0.6\). Then the probability that event type 1 will happen will be \(\frac{0.3}{0.3+0.6}=\frac{1}{3}\) and probability of type 2 event is \(\frac{2}{3}\). We generate type 1 events and type 2 event by binomial distribution.
f.event<- rbinom(100, size = 1, prob = 1/3)
f.event<- ifelse(f.event== 0, 2, 1)
f.event
## [1] 2 1 2 2 2 2 1 1 1 1 2 2 2 1 1 2 1 2 2 2 2 2 2 2 2 1 1 2 2 2 1 1 2 2 2 2 1
## [38] 2 2 2 2 2 2 1 1 1 2 2 2 1 2 2 1 1 1 2 2 2 2 1 2 1 1 1 1 2 1 2 2 2 1 1 2 2
## [75] 2 1 2 2 1 2 2 2 1 2 1 1 2 2 2 2 2 2 2 2 1 2 1 2 2 2
cens.times <- runif(100,0,5)
cens.times
## [1] 4.40325563 3.57977152 1.60662329 1.86146735 0.03073766 4.20676978
## [7] 4.55709257 0.05759936 0.39521799 2.96777923 1.84505252 4.82081719
## [13] 2.62180673 3.68355170 3.24768931 4.21777225 0.32273913 0.30977610
## [19] 0.93340404 1.85665686 4.53449560 4.12867487 2.47397466 3.72439910
## [25] 2.19494584 0.72841787 3.35466714 2.76698896 4.96013970 1.17249203
## [31] 2.69579535 1.28077460 2.04829145 4.95434282 3.73735641 0.99198969
## [37] 4.40487373 1.55842751 1.45689337 0.98854275 3.04720414 3.52224209
## [43] 1.66816887 2.76469070 4.92168653 1.24653166 0.43360143 3.88143215
## [49] 2.41625411 2.99369831 2.68472381 0.19853639 1.50070680 2.26652029
## [55] 3.15149814 0.32219304 1.38156118 4.29247548 3.29898014 4.92913267
## [61] 4.17064182 1.95180268 3.55351071 0.15533860 3.93090686 1.62999291
## [67] 0.11450412 3.80276020 1.17126984 3.26734149 4.53380904 2.19252096
## [73] 2.55235383 3.20062564 0.47773395 0.47499943 0.86354950 3.84234087
## [79] 4.09637422 3.41190787 1.29693795 4.30458435 4.16669638 4.07686430
## [85] 4.22044816 1.43214285 3.16798584 0.62323298 0.21968377 1.61619385
## [91] 0.57218792 1.81918569 0.36828312 2.65724331 1.17148941 2.91036744
## [97] 0.23353197 1.72187895 1.10347305 0.23977229
Combine the event time and censored time, and the minimum time is the observed time.
obs.times <- pmin(event.times, cens.times)
obs.times
## [1] 0.031468670 1.699722626 0.012604250 0.181851015 0.030737664 0.600348067
## [7] 0.253768509 0.057599361 0.239200077 2.967779232 0.045875484 0.307674001
## [13] 0.420839034 0.286222322 0.084884797 0.043504038 0.078568163 0.309776096
## [19] 0.933404035 1.383426225 1.634379621 1.224719648 1.102386634 1.590425714
## [25] 0.332396754 0.728417869 0.141629408 0.097233212 3.595864890 0.540339829
## [31] 0.627334975 0.505893539 0.296015002 0.274083299 2.366269630 0.991989690
## [37] 1.241166406 0.116164153 0.215933847 0.988542751 2.758792475 0.007852311
## [43] 1.668168872 0.118259397 1.112956394 1.246531657 0.219232355 1.906241483
## [49] 0.565370676 1.077138985 0.484502960 0.198536394 1.500706801 0.853008879
## [55] 0.825220779 0.322193044 1.021393551 4.096878624 0.129284012 0.507789280
## [61] 0.065935818 1.951802684 0.183929479 0.155338600 0.604041278 1.629992907
## [67] 0.044561327 0.581029905 0.162785934 1.398531164 2.622305360 2.177511017
## [73] 0.053606869 0.302974220 0.187704037 0.474999429 0.710940666 0.271044257
## [79] 1.606070770 0.751399968 0.480695701 3.406351614 0.198755058 0.237405763
## [85] 0.121899657 1.054376699 1.822300313 0.623232982 0.198321935 1.616193851
## [91] 0.516786499 1.255904736 0.051519817 0.808046473 0.039023547 0.119380507
## [97] 0.182533470 0.124188870 0.227555450 0.239772290
Assign event/censor indicators
obs.event <- c(event.times <= cens.times) * f.event
obs.event
## [1] 2 1 2 2 0 2 1 0 1 0 2 2 2 1 1 2 1 0 0 2 2 2 2 2 2 0 1 2 2 2 1 1 2 2 2 0 1
## [38] 2 2 0 2 2 0 1 1 0 2 2 2 1 2 0 0 1 1 0 2 2 2 1 2 0 1 0 1 0 1 2 2 2 1 1 2 2
## [75] 2 0 2 2 1 2 2 2 1 2 1 1 2 0 2 0 2 2 2 2 1 2 1 2 2 0
The final data set will be
final_data<-data.frame(obs.times,obs.event)
final_data
## obs.times obs.event
## 1 0.031468670 2
## 2 1.699722626 1
## 3 0.012604250 2
## 4 0.181851015 2
## 5 0.030737664 0
## 6 0.600348067 2
## 7 0.253768509 1
## 8 0.057599361 0
## 9 0.239200077 1
## 10 2.967779232 0
## 11 0.045875484 2
## 12 0.307674001 2
## 13 0.420839034 2
## 14 0.286222322 1
## 15 0.084884797 1
## 16 0.043504038 2
## 17 0.078568163 1
## 18 0.309776096 0
## 19 0.933404035 0
## 20 1.383426225 2
## 21 1.634379621 2
## 22 1.224719648 2
## 23 1.102386634 2
## 24 1.590425714 2
## 25 0.332396754 2
## 26 0.728417869 0
## 27 0.141629408 1
## 28 0.097233212 2
## 29 3.595864890 2
## 30 0.540339829 2
## 31 0.627334975 1
## 32 0.505893539 1
## 33 0.296015002 2
## 34 0.274083299 2
## 35 2.366269630 2
## 36 0.991989690 0
## 37 1.241166406 1
## 38 0.116164153 2
## 39 0.215933847 2
## 40 0.988542751 0
## 41 2.758792475 2
## 42 0.007852311 2
## 43 1.668168872 0
## 44 0.118259397 1
## 45 1.112956394 1
## 46 1.246531657 0
## 47 0.219232355 2
## 48 1.906241483 2
## 49 0.565370676 2
## 50 1.077138985 1
## 51 0.484502960 2
## 52 0.198536394 0
## 53 1.500706801 0
## 54 0.853008879 1
## 55 0.825220779 1
## 56 0.322193044 0
## 57 1.021393551 2
## 58 4.096878624 2
## 59 0.129284012 2
## 60 0.507789280 1
## 61 0.065935818 2
## 62 1.951802684 0
## 63 0.183929479 1
## 64 0.155338600 0
## 65 0.604041278 1
## 66 1.629992907 0
## 67 0.044561327 1
## 68 0.581029905 2
## 69 0.162785934 2
## 70 1.398531164 2
## 71 2.622305360 1
## 72 2.177511017 1
## 73 0.053606869 2
## 74 0.302974220 2
## 75 0.187704037 2
## 76 0.474999429 0
## 77 0.710940666 2
## 78 0.271044257 2
## 79 1.606070770 1
## 80 0.751399968 2
## 81 0.480695701 2
## 82 3.406351614 2
## 83 0.198755058 1
## 84 0.237405763 2
## 85 0.121899657 1
## 86 1.054376699 1
## 87 1.822300313 2
## 88 0.623232982 0
## 89 0.198321935 2
## 90 1.616193851 0
## 91 0.516786499 2
## 92 1.255904736 2
## 93 0.051519817 2
## 94 0.808046473 2
## 95 0.039023547 1
## 96 0.119380507 2
## 97 0.182533470 1
## 98 0.124188870 2
## 99 0.227555450 2
## 100 0.239772290 0