## Settings for RMarkdown http://yihui.name/knitr/options#chunk_options
opts_chunk$set(comment = "", warning = FALSE, message = FALSE, tidy = FALSE,
echo = TRUE, fig.width = 10, fig.height = 8)
options(width = 116, scipen = 10)
## Introductory Statistics with R 2nd ed. (Peter Dalgaard)のデータを使っています。
## 1版は和訳 Rによる医療統計学 (岡田 昌史) もあります。
## 使い方
## このコードをコピーして手元のRのエディタに貼付ける。
## Windows版では左上 [ファイル] → [新しいスクリプト]でエディタ。
## Mac版では真っ白い紙の形のアイコンでエディタ。
## 実行する行、あるいは範囲を選択して
## WindowsであればCtrl＋R、MacであればCommand＋Returnで実行。
## インストール用の関数定義
## rownames(installed.packages()でインストールされているパッケージの一覧
## PKG %in% ... という表現は...の中にPKGが含まれればTRUEを返します。!をつけて評価を逆転。
inst <- function (PKG) {
if (!(PKG %in% rownames(installed.packages()))) {
install.packages(pkgs = PKG, dependencies = TRUE)
}
}
## Introductory Statistics with Rのデータ集(ISwR)をインストール。
inst("ISwR")
## library()はパッケージを読み込むコマンド。ISwRが上記の本のデータ集パッケージです。
library(ISwR)
## EstoniaのTartuという都市の脳血管障害のデータを使用します。
## J. Korv, M. Roose, and A.E. Kaasik (1997). Stroke Registry of Tartu, Estonia, from 1991 through 1993. Cerebrovascular Disorders 7:154–162.
## http://content.karger.com/ProdukteDB/produkte.asp?Aktion=ShowAbstract&ArtikelNr=108182&Ausgabe=232923&ProduktNr=224153
data(stroke)
## データベースの説明は?strokeでみられます。
?stroke
## strokeから必要な変数のみかつcomplete casesを抜き出す n=814
## stroke[ ,c(-2,-3)] というのはstroke[行, 列]という形で抜き出します。
## マイナス指定をすると2,3列目以外という意味です。
## Strokeという名前の新しいデータベースにしました。
Stroke <- stroke[complete.cases(stroke[,c(-2,-3)]),c(-2,-3)]
## 先頭10行をhead()で表示してみます。hanは高血圧既往のようです。
head(Stroke, n = 10)
sex age dgn coma diab minf han dead obsmonths
1 Male 76 INF No No Yes No TRUE 0.16340
2 Male 58 INF No No No No FALSE 59.60784
3 Male 74 INF No No Yes Yes TRUE 4.73856
4 Female 77 ICH No Yes No Yes TRUE 0.06536
5 Female 76 INF No Yes No Yes FALSE 59.28105
6 Male 48 ICH Yes No No Yes TRUE 0.10000
7 Female 81 INF No No No Yes TRUE 34.37908
8 Male 53 INF No No Yes Yes TRUE 10.84967
9 Female 73 ID No No No Yes FALSE 59.21569
10 Female 69 INF No No No Yes TRUE 33.66013
## column名をわかりやすく書き換える。diagnosisはICH, ID (unidentified), INF (INFarction), SAH
## ちょっとわかりにくいですが、names(Strokes)でcolumn名の一覧が出るのでそれを置き換えるという方法。
names(Stroke) <- c("sex","age","diagnosis","coma","diabetes","MI","HTN","dead","time_mo")
## サマリー機能を使ってみましょう。集計やquartileがでます。
summary(Stroke)
sex age diagnosis coma diabetes MI HTN dead time_mo
Female:500 Min. : 1.0 ICH: 79 No :735 No :717 No :718 No :405 Mode :logical Min. : 0.03
Male :314 1st Qu.:61.0 ID :190 Yes: 79 Yes: 97 Yes: 96 Yes:409 FALSE:343 1st Qu.: 0.60
Median :71.0 INF:498 TRUE :471 Median :21.68
Mean :69.7 SAH: 47 NA's :0 Mean :20.93
3rd Qu.:81.0 3rd Qu.:37.12
Max. :96.0 Max. :59.61
## 生存分析用のsurvivalパッケージを読み込みます。元々入っているのでインストール不要。
library(survival)
## survplotというnumber at riskを表示してくれるグラフかきのコマンドのためrmsパッケージインストール。
## 昔はDesignを使っていましたがパッケージがなくなったようです。
inst("rms") ; library(rms)
## 性別で群分けしただけのKM estimator。 outcome 〜 sexの単説明変数モデルを作ります。
## survfit(Surv(時間変数, イベント変数) ~ 層別化因子, data = データ名)の書式です。
## 群分けしない時は ~ 1としてください。
fit.1 <- survfit(Surv(time_mo, dead) ~ sex, data = Stroke)
## プロットします。
## CONFidence intervalなし、Number at RISKあり、時間の刻み6ヶ月です。
survplot(fit.1,
conf = "none", n.risk = TRUE, xlab = "Months", time.inc = 6)
## log rank検定では有意差あり。ただしこれは交絡因子を考慮していない。
## survdiff(Surv(時間変数, イベント変数) ~ 層別化因子, data = データ名)の書式です。
survdiff(Surv(time_mo, dead) ~ sex, data = Stroke)
Call:
survdiff(formula = Surv(time_mo, dead) ~ sex, data = Stroke)
N Observed Expected (O-E)^2/E (O-E)^2/V
sex=Female 500 312 271 6.17 14.7
sex=Male 314 159 200 8.36 14.7
Chisq= 14.7 on 1 degrees of freedom, p= 0.000125
## 年齢の影響を階層化で取り除けるか見てみます。
## age年齢の変数を65歳で分割して因子にしています。FALSE,TRUEになるのでYoung, Oldとラベルを付けます。
## with(Stroke, ...)とすることでStroke$ageと書かないで済みます。
Stroke$old <- with(Stroke, factor(age >= 65, levels = c(FALSE, TRUE), labels=c("Young","Old")))
fit.2 <- survfit(Surv(time_mo, dead) ~ sex + old, data = Stroke)
survplot(fit.2, conf = "none", n.risk = TRUE, xlab = "Months", time.inc = 6)
## グラフを詳細に指定するデモ。
## そのまま表示すると各群のn、イベント数、median survival(50%死んでいないとでない)がでる。
fit.2
Call: survfit.formula(formula = Surv(time_mo, dead) ~ sex + old, data = Stroke)
records n.max n.start events median 0.95LCL 0.95UCL
sex=Female, old=Young 114 114 114 37 NA NA NA
sex=Female, old=Old 386 386 386 275 3.81 2.09 7.25
sex=Male, old=Young 166 166 166 58 NA 51.90 NA
sex=Male, old=Old 148 148 148 101 13.01 5.56 21.44
## summary()を使うと、各群のイベント発生時ごとのsurvival estimationがでる。
summary(fit.2)
Call: survfit.formula(formula = Surv(time_mo, dead) ~ sex + old, data = Stroke)
sex=Female, old=Young
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0.0327 114 1 0.991 0.00873 0.974 1.000
0.0654 113 2 0.974 0.01499 0.945 1.000
0.0980 111 2 0.956 0.01918 0.919 0.994
0.1000 109 5 0.912 0.02649 0.862 0.966
0.1961 104 1 0.904 0.02765 0.851 0.959
0.2288 103 2 0.886 0.02977 0.829 0.946
0.2614 101 3 0.860 0.03253 0.798 0.926
0.2941 98 3 0.833 0.03490 0.768 0.905
0.3268 95 2 0.816 0.03631 0.748 0.890
0.4248 93 1 0.807 0.03696 0.738 0.883
0.4575 92 1 0.798 0.03759 0.728 0.875
0.8824 91 1 0.789 0.03818 0.718 0.868
1.1111 90 1 0.781 0.03875 0.708 0.860
1.3399 89 1 0.772 0.03930 0.699 0.853
1.6340 88 1 0.763 0.03982 0.689 0.845
2.3203 87 1 0.754 0.04032 0.679 0.838
2.8758 86 1 0.746 0.04079 0.670 0.830
3.9216 85 1 0.737 0.04124 0.660 0.822
6.0458 84 1 0.728 0.04167 0.651 0.815
20.3268 83 1 0.719 0.04208 0.641 0.807
21.6667 82 1 0.711 0.04248 0.632 0.799
26.6013 72 1 0.701 0.04302 0.621 0.790
31.8627 55 1 0.688 0.04408 0.607 0.780
33.3007 51 1 0.674 0.04523 0.591 0.769
34.7059 45 1 0.659 0.04664 0.574 0.758
sex=Female, old=Old
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0.0327 386 12 0.969 0.00883 0.952 0.986
0.0654 374 13 0.935 0.01253 0.911 0.960
0.0980 361 6 0.920 0.01383 0.893 0.947
0.1000 355 8 0.899 0.01534 0.869 0.930
0.1307 347 8 0.878 0.01664 0.846 0.911
0.1634 339 5 0.865 0.01738 0.832 0.900
0.1961 334 7 0.847 0.01832 0.812 0.884
0.2288 327 7 0.829 0.01916 0.792 0.867
0.2614 320 11 0.801 0.02034 0.762 0.841
0.2941 309 7 0.782 0.02100 0.742 0.825
0.3268 302 8 0.762 0.02169 0.720 0.805
0.3595 294 3 0.754 0.02192 0.712 0.798
0.3922 291 7 0.736 0.02244 0.693 0.781
0.4248 284 6 0.720 0.02285 0.677 0.766
0.4575 278 4 0.710 0.02310 0.666 0.757
0.4902 274 2 0.705 0.02322 0.661 0.752
0.5229 272 4 0.694 0.02345 0.650 0.742
0.5882 268 3 0.687 0.02361 0.642 0.734
0.6209 265 1 0.684 0.02366 0.639 0.732
0.6536 264 2 0.679 0.02377 0.634 0.727
0.6863 262 3 0.671 0.02392 0.626 0.720
0.7516 259 3 0.663 0.02406 0.618 0.712
0.7843 256 4 0.653 0.02423 0.607 0.702
0.8170 252 4 0.642 0.02439 0.596 0.692
0.8497 248 1 0.640 0.02443 0.594 0.690
0.8824 247 2 0.635 0.02451 0.588 0.685
0.9150 245 2 0.630 0.02458 0.583 0.680
0.9477 243 1 0.627 0.02462 0.581 0.677
0.9804 242 3 0.619 0.02472 0.573 0.670
1.0131 239 1 0.617 0.02475 0.570 0.667
1.0784 238 3 0.609 0.02484 0.562 0.659
1.1111 235 3 0.601 0.02492 0.554 0.652
1.1765 232 3 0.593 0.02500 0.546 0.644
1.2418 229 3 0.585 0.02507 0.538 0.637
1.2745 226 3 0.578 0.02514 0.530 0.629
1.3072 223 2 0.573 0.02518 0.525 0.624
1.3725 221 1 0.570 0.02520 0.523 0.622
1.4379 220 2 0.565 0.02523 0.517 0.616
1.5686 218 2 0.560 0.02527 0.512 0.611
1.6340 216 2 0.554 0.02530 0.507 0.606
1.6667 214 2 0.549 0.02533 0.502 0.601
2.0915 212 1 0.547 0.02534 0.499 0.599
2.2222 211 2 0.541 0.02536 0.494 0.594
2.2549 209 1 0.539 0.02537 0.491 0.591
2.4510 208 1 0.536 0.02538 0.489 0.588
2.4837 207 1 0.534 0.02539 0.486 0.586
2.5817 206 2 0.528 0.02541 0.481 0.581
2.7451 204 1 0.526 0.02542 0.478 0.578
2.7778 203 2 0.521 0.02543 0.473 0.573
2.8105 201 1 0.518 0.02543 0.471 0.570
2.9085 200 1 0.516 0.02544 0.468 0.568
3.0719 199 1 0.513 0.02544 0.465 0.565
3.1046 198 1 0.510 0.02544 0.463 0.563
3.3333 197 2 0.505 0.02545 0.458 0.558
3.6275 195 1 0.503 0.02545 0.455 0.555
3.7582 194 1 0.500 0.02545 0.453 0.552
3.8562 193 1 0.497 0.02545 0.450 0.550
4.0196 192 1 0.495 0.02545 0.447 0.547
4.6078 191 1 0.492 0.02545 0.445 0.545
4.6405 190 2 0.487 0.02544 0.440 0.540
4.6732 188 1 0.484 0.02544 0.437 0.537
4.7059 187 2 0.479 0.02543 0.432 0.532
4.8039 185 1 0.477 0.02542 0.429 0.529
4.9346 184 1 0.474 0.02542 0.427 0.527
5.0654 183 1 0.472 0.02541 0.424 0.524
5.5882 182 1 0.469 0.02540 0.422 0.521
5.6536 181 1 0.466 0.02539 0.419 0.519
6.0131 180 1 0.464 0.02538 0.417 0.516
6.1765 179 1 0.461 0.02537 0.414 0.514
6.4379 178 1 0.459 0.02536 0.411 0.511
6.6013 177 1 0.456 0.02535 0.409 0.508
7.0915 176 1 0.453 0.02534 0.406 0.506
7.2222 175 2 0.448 0.02531 0.401 0.501
7.2549 173 1 0.446 0.02530 0.399 0.498
7.3203 172 1 0.443 0.02528 0.396 0.495
7.5490 171 1 0.440 0.02527 0.394 0.493
8.5294 170 1 0.438 0.02525 0.391 0.490
9.5098 169 1 0.435 0.02523 0.388 0.488
9.9346 168 1 0.433 0.02522 0.386 0.485
10.1961 167 1 0.430 0.02520 0.383 0.482
10.4575 166 1 0.427 0.02518 0.381 0.480
12.8105 165 1 0.425 0.02516 0.378 0.477
13.3660 164 1 0.422 0.02514 0.376 0.475
13.5621 163 1 0.420 0.02512 0.373 0.472
13.6275 162 1 0.417 0.02510 0.371 0.469
13.9869 161 1 0.415 0.02507 0.368 0.467
15.6863 160 1 0.412 0.02505 0.366 0.464
15.7190 159 1 0.409 0.02503 0.363 0.461
15.9150 158 1 0.407 0.02500 0.361 0.459
16.0131 157 1 0.404 0.02498 0.358 0.456
16.2745 156 1 0.402 0.02495 0.356 0.454
16.3725 155 1 0.399 0.02492 0.353 0.451
17.0261 154 1 0.396 0.02490 0.350 0.448
18.2353 153 1 0.394 0.02487 0.348 0.446
18.5621 152 1 0.391 0.02484 0.345 0.443
18.6275 151 1 0.389 0.02481 0.343 0.440
18.9542 150 1 0.386 0.02478 0.340 0.438
20.0327 149 1 0.383 0.02475 0.338 0.435
20.3922 148 1 0.381 0.02472 0.335 0.432
20.5556 147 1 0.378 0.02468 0.333 0.430
20.8497 146 1 0.376 0.02465 0.330 0.427
21.6993 145 1 0.373 0.02462 0.328 0.425
22.4837 144 1 0.370 0.02458 0.325 0.422
23.5948 143 1 0.368 0.02454 0.323 0.419
23.6601 142 1 0.365 0.02451 0.320 0.417
24.2157 139 1 0.363 0.02447 0.318 0.414
24.3137 137 1 0.360 0.02444 0.315 0.411
24.3464 136 1 0.357 0.02440 0.313 0.409
24.4444 135 1 0.355 0.02436 0.310 0.406
24.5098 133 1 0.352 0.02432 0.307 0.403
25.0327 131 1 0.349 0.02429 0.305 0.400
25.1961 129 1 0.347 0.02425 0.302 0.398
26.4379 125 1 0.344 0.02421 0.300 0.395
26.5359 124 1 0.341 0.02418 0.297 0.392
27.5163 120 1 0.338 0.02414 0.294 0.389
27.8758 119 1 0.335 0.02411 0.291 0.386
27.9085 118 1 0.333 0.02407 0.289 0.383
28.1046 116 1 0.330 0.02403 0.286 0.380
28.1373 115 1 0.327 0.02399 0.283 0.377
29.1503 109 1 0.324 0.02396 0.280 0.374
30.4248 105 1 0.321 0.02393 0.277 0.371
30.7516 103 1 0.318 0.02390 0.274 0.368
31.8627 97 1 0.314 0.02387 0.271 0.365
32.8105 93 1 0.311 0.02386 0.268 0.361
32.9739 92 1 0.308 0.02383 0.264 0.358
33.6601 88 1 0.304 0.02382 0.261 0.355
34.3791 85 1 0.301 0.02381 0.257 0.351
35.1634 79 1 0.297 0.02381 0.254 0.347
36.4379 72 1 0.293 0.02383 0.249 0.343
38.4641 66 1 0.288 0.02388 0.245 0.339
38.7908 64 1 0.284 0.02393 0.240 0.335
39.8039 62 1 0.279 0.02397 0.236 0.330
41.3399 57 1 0.274 0.02405 0.231 0.326
42.8431 54 1 0.269 0.02413 0.226 0.321
47.1242 35 1 0.261 0.02464 0.217 0.314
sex=Male, old=Young
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0.0327 166 3 0.982 0.0103 0.962 1.000
0.0654 163 3 0.964 0.0145 0.936 0.993
0.0980 160 1 0.958 0.0156 0.928 0.989
0.1000 159 8 0.910 0.0223 0.867 0.954
0.2614 151 1 0.904 0.0229 0.860 0.950
0.2941 150 1 0.898 0.0235 0.853 0.945
0.3595 149 1 0.892 0.0241 0.845 0.940
0.3922 148 2 0.880 0.0253 0.831 0.930
0.6863 146 1 0.873 0.0258 0.824 0.926
0.8170 145 1 0.867 0.0263 0.817 0.921
1.3725 144 1 0.861 0.0268 0.810 0.916
1.4052 143 1 0.855 0.0273 0.804 0.911
1.6667 142 1 0.849 0.0278 0.797 0.906
1.9281 141 1 0.843 0.0282 0.790 0.901
2.1895 140 1 0.837 0.0286 0.783 0.895
3.3007 139 1 0.831 0.0291 0.776 0.890
3.4641 138 1 0.825 0.0295 0.770 0.885
3.9542 137 2 0.813 0.0302 0.756 0.875
4.3791 135 1 0.807 0.0306 0.749 0.870
5.5229 134 1 0.801 0.0310 0.743 0.864
6.4706 133 1 0.795 0.0313 0.736 0.859
7.2222 132 1 0.789 0.0317 0.729 0.854
7.3203 131 1 0.783 0.0320 0.723 0.848
7.4837 130 1 0.777 0.0323 0.716 0.843
7.7451 129 1 0.771 0.0326 0.710 0.838
9.0196 128 1 0.765 0.0329 0.703 0.832
9.5098 127 1 0.759 0.0332 0.697 0.827
10.0980 126 1 0.753 0.0335 0.690 0.822
10.8497 125 1 0.747 0.0337 0.684 0.816
10.8824 124 1 0.741 0.0340 0.677 0.811
11.2418 123 1 0.735 0.0343 0.671 0.805
15.9150 122 1 0.729 0.0345 0.664 0.800
19.7059 121 1 0.723 0.0347 0.658 0.794
21.9608 120 1 0.717 0.0350 0.652 0.789
22.1569 119 1 0.711 0.0352 0.645 0.783
25.7843 111 1 0.704 0.0354 0.638 0.777
33.3333 89 1 0.697 0.0359 0.630 0.771
34.5098 88 1 0.689 0.0364 0.621 0.764
37.4510 73 1 0.679 0.0371 0.610 0.756
43.6275 51 1 0.666 0.0387 0.594 0.746
49.2810 38 1 0.648 0.0414 0.572 0.735
51.1765 31 1 0.627 0.0451 0.545 0.722
51.2092 30 1 0.607 0.0482 0.519 0.709
51.8954 27 1 0.584 0.0514 0.492 0.694
53.6601 17 1 0.550 0.0587 0.446 0.678
sex=Male, old=Old
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0.0327 148 3 0.980 0.0116 0.957 1.000
0.0654 145 2 0.966 0.0149 0.938 0.996
0.0980 143 4 0.939 0.0196 0.901 0.978
0.1000 139 4 0.912 0.0233 0.868 0.959
0.1307 135 1 0.905 0.0241 0.859 0.954
0.1634 134 2 0.892 0.0255 0.843 0.943
0.1961 132 5 0.858 0.0287 0.804 0.916
0.2288 127 4 0.831 0.0308 0.773 0.894
0.2614 123 2 0.818 0.0317 0.758 0.882
0.2941 121 1 0.811 0.0322 0.750 0.876
0.3268 120 3 0.791 0.0334 0.728 0.859
0.3595 117 3 0.770 0.0346 0.705 0.841
0.3922 114 1 0.764 0.0349 0.698 0.835
0.4575 113 1 0.757 0.0353 0.691 0.829
0.4902 112 1 0.750 0.0356 0.683 0.823
0.5229 111 1 0.743 0.0359 0.676 0.817
0.5556 110 2 0.730 0.0365 0.662 0.805
0.6536 108 2 0.716 0.0371 0.647 0.793
0.7190 106 2 0.703 0.0376 0.633 0.780
0.7516 104 1 0.696 0.0378 0.626 0.774
0.8170 103 1 0.689 0.0380 0.619 0.768
0.8824 102 1 0.682 0.0383 0.611 0.762
1.0458 101 1 0.676 0.0385 0.604 0.755
1.0784 100 1 0.669 0.0387 0.597 0.749
1.3399 99 1 0.662 0.0389 0.590 0.743
1.7320 98 1 0.655 0.0391 0.583 0.737
1.7974 97 1 0.649 0.0392 0.576 0.730
2.1569 96 1 0.642 0.0394 0.569 0.724
2.3856 95 1 0.635 0.0396 0.562 0.718
2.8431 94 1 0.628 0.0397 0.555 0.711
2.9085 93 1 0.622 0.0399 0.548 0.705
3.9216 92 1 0.615 0.0400 0.541 0.698
4.2157 91 1 0.608 0.0401 0.534 0.692
4.2810 90 1 0.601 0.0402 0.527 0.686
4.3464 89 1 0.595 0.0404 0.521 0.679
4.6078 88 1 0.588 0.0405 0.514 0.673
4.7386 87 1 0.581 0.0406 0.507 0.666
5.5556 86 1 0.574 0.0406 0.500 0.660
5.6209 85 1 0.568 0.0407 0.493 0.653
6.6340 84 1 0.561 0.0408 0.486 0.647
7.6144 83 2 0.547 0.0409 0.473 0.634
8.0719 81 1 0.541 0.0410 0.466 0.627
8.4967 80 1 0.534 0.0410 0.459 0.621
9.3464 79 1 0.527 0.0410 0.452 0.614
10.2614 78 1 0.520 0.0411 0.446 0.607
10.2941 77 1 0.514 0.0411 0.439 0.601
11.3725 76 1 0.507 0.0411 0.432 0.594
12.9085 75 1 0.500 0.0411 0.426 0.587
13.1046 74 1 0.493 0.0411 0.419 0.581
13.2353 73 1 0.486 0.0411 0.412 0.574
13.8889 72 1 0.480 0.0411 0.406 0.567
13.9542 71 1 0.473 0.0410 0.399 0.561
14.9673 70 1 0.466 0.0410 0.392 0.554
15.9477 69 1 0.459 0.0410 0.386 0.547
15.9804 68 1 0.453 0.0409 0.379 0.540
17.0261 67 1 0.446 0.0409 0.373 0.534
18.0392 66 1 0.439 0.0408 0.366 0.527
19.6732 65 1 0.432 0.0407 0.360 0.520
19.9020 64 1 0.426 0.0406 0.353 0.513
20.0000 63 1 0.419 0.0406 0.347 0.506
21.4379 62 1 0.412 0.0405 0.340 0.500
22.4510 61 1 0.405 0.0404 0.334 0.493
23.6275 60 1 0.399 0.0402 0.327 0.486
24.1830 59 1 0.392 0.0401 0.321 0.479
24.6405 57 1 0.385 0.0400 0.314 0.472
24.6732 56 1 0.378 0.0399 0.308 0.465
24.9020 54 1 0.371 0.0398 0.301 0.458
29.3137 39 1 0.362 0.0399 0.291 0.449
29.7059 37 1 0.352 0.0400 0.282 0.440
29.7712 36 1 0.342 0.0400 0.272 0.430
33.2026 30 1 0.331 0.0403 0.260 0.420
35.9150 25 1 0.317 0.0408 0.247 0.408
35.9804 24 1 0.304 0.0412 0.233 0.397
39.8693 22 1 0.290 0.0416 0.219 0.384
42.5490 20 1 0.276 0.0419 0.205 0.372
## そのままグラフ化するといまいち。X軸も単位が間違っている。
survplot(fit.2)
## これぐらい指定すればとりあえず見れる。
survplot(fit.2,
xlab = "Month", # X軸が標準でDaysになってしまうのでmonthにするあわせる
conf = "none", # 信頼区間は複数グラフのときは外す
time.inc = 6, # 何ヶ月ごとにX軸に数字をだすか
n.risk = TRUE # number at riskを表示する
)
## さらに詳細な設定で見やすいグラフを作成。
par()$mar # グラフ周囲の空白の幅。デフォルト表示。5.1 4.1 4.1 4.1。
[1] 5.1 4.1 4.1 4.1
par(mar=c(12,4,4,10)) # グラフ周囲のMARginの幅。下、左、上、右の順番で指定。下を広く取る。
par(las=1) # Y軸の数字を横書きにLabel of Axis Styleの略。
survplot(fit.2, # survfit()で作ったオブジェクト
xlab = "Month", # X軸が標準でDaysになってしまうのでmonthにするあわせる
conf = "none", # 信頼区間は複数グラフのときは外す
time.inc = 6, # 何ヶ月ごとにX軸に数字をだすか
n.risk = TRUE, # n at riskを表示する
srt.n.risk = 0, # n at riskの左端の開始位置、デフォルトでよいので書かないでよい
sep.n.risk = 0.056, # n at riskの行の隙間、デフォルトが0.056、狭くしたければ書く。
adj.n.risk = 1, # n at riskの左端の数字の位置。デフォルトでよいので書かないでよい。
y.n.risk = -0.5, # n at riskのY軸位置。-0.5にするとX軸の下側にいく
cex.n.risk = 1.0, # n at riskの文字の大きさ。1だと軸の数字と同じ大きさになる。
label.curves = list(method = "arrow", cex = 1.5) # 矢印でラベル指定、list()で方法と文字サイズを入れます。
)
## titleをつけます。
title("Survival curves from stroke dataset")
## おまけ
## listは一つのオプションに二つ以上の指定をするときに使われることがあります。
## 下記のごとく二つの見出しにそれぞれ値が指定されています。
list(method = "arrow", cex = 0.8)
$method
[1] "arrow"
$cex
[1] 0.8
## おまけ2
## EZRの中にnrisk()という関数があって、
## nrisk(SURVFIT, c(0, 6, 8, 12)) とかやると任意の時間でのn at riskが取り出せます。
## nrisk() from EZR by Dr. Yoshinobu Kanda
## http://www.jichi.ac.jp/saitama-sct/SaitamaHP.files/statmed.html
## Originally in survplot package. Changed to work with single stratum analysis
nrisk <- function (x, times = pretty(x$time)) {
if (!is.null(x$strata)){
ns <- length(x$strata)
idx <- rep.int(1:ns, x$strata)
str.n.risk <- split(x$n.risk, idx)
str.times <- split(x$time, idx)
m <- sapply(times, function(y) {
sapply(1:ns, function(i) {
w <- which(str.times[[i]] >= y)[1]
ifelse(is.na(w), 0, str.n.risk[[i]][w])
})
})
rownames(m) <- names(x$strata)
colnames(m) <- times
} else {
str.n.risk <- x$n.risk
str.times <- x$time
m <- sapply(times, function(y) {
w <- which(str.times >= y)[1]
ifelse(is.na(w), 0, str.n.risk[w])
})
}
m
}
nrisk(fit.2, seq(from = 0, to = 60, by = 3))
0 3 6 9 12 15 18 21 24 27 30 33 36 39 42 45 48 51 54 57 60
sex=Female, old=Young 114 85 84 83 83 83 83 82 81 71 60 51 41 33 24 19 12 8 5 2 0
sex=Female, old=Old 386 199 180 169 165 160 153 145 140 122 108 91 74 62 55 40 31 22 15 8 0
sex=Male, old=Young 166 139 133 128 122 122 121 120 118 106 96 89 81 66 59 46 39 33 16 4 0
sex=Male, old=Old 148 92 84 79 75 69 66 62 59 48 35 31 23 22 20 17 12 8 7 2 0