「統計学関連なんでもあり」の過去ログ--- 047

No.21900　KM図のat risk表示方法　　【赤岳】　2016/02/01(Mon) 20:46

いつもご教示ありがとうございます。
EXZでKM図を描いて，その図の下にAt riskを表示させました。
下記がそのプログラムです。
しかし，At riskは，0, 2, 4, 6, 8yearと，2年ごとにしか表示されません。
このプログラムを使って，0,1,2・・・8yearと，1年ごとに表示させたいのですが，どこをどのように修正すれば1年ごとにAt riskが表示されるようになるのでしょうか。
どなたか，ご教示お願いします。

> #####他の因子で調整した生存曲線の表示#####
> coxmodel <- coxph(Surv(day, event==1)~ age + skill + strata(TTT), data=Dataset,
+ method="breslow")

以下削除

No.21901　Re: KM図のat risk表示方法　　【青木繁伸】　2016/02/01(Mon) 21:12

スクリプト（プログラム）だけをコピー・ペーストしただけであなたの分析を再現できるスクリプトを提示してください。

> library(survival)
> #####他の因子で調整した生存曲線の表示#####
> coxmodel <- coxph(Surv(day, event==1)~ age + skill + strata(TTT), data=Dataset,
+ method="breslow")
terms.formula(formula, special, data = data) でエラー:
オブジェクト 'Dataset' がありません

以下削除

No.21902　Re: KM図のat risk表示方法　　【赤岳】　2016/02/01(Mon) 22:49

青木先生，お世話になっております。早速のご対応いただきありがとうございます。
EZRを使って図を描きました。そのときのスクリプトは，以下のとおりです。

> #####他の因子で調整した生存曲線の表示#####
> coxmodel <- coxph(Surv(day, event==1)~ age + skill + strata(TTT), data=Dataset,
+ method="breslow")
> cox <- survfit(coxmodel)
> windows(width=7, height=7); par(lwd=1, las=1, family="sans", cex=1, mgp=c(2.5,1,0))
> group.levels <- names(cox$strata[cox$strata>0])
> mar <- par("mar")
> mar[1] <- mar[1] + length(cox$strata) + 0.5
> mar[2] <- mar[2] + 2
> par(mar=mar)
> opar <- par(mar = mar)
> on.exit(par(opar))
> plot(cox, col=1:32, lty=1, lwd=1, bty="l", mark.time=TRUE, xlab="year", ylab="Probability",
+ xscale=365.25)
> xticks <- axTicks(1)
> n.atrisk <- nrisk(cox, xticks * 365.25)
> for (i in 1:length(cox$strata)){axis(1, at = xticks, labels = n.atrisk[i,], line=3+i, tick =
+ FALSE)}
> for (i in 1:length(cox$strata)){mtext(group.levels[i], at=-(xticks[2]-xticks[1])/2, side=1,
+ line=4+i, cex=1)}
> title(xlab = "Number at risk", line = 3.5, adj = 0)
> legend("topright", group.levels, title="TTT", col=1:32, lty=1, lwd=1, box.lty=0)
> title("Survival curve")
> summary(cox)

No.21903　Re: KM図のat risk表示方法　　【青木繁伸】　2016/02/02(Tue) 06:35

Dataset というデータフレームがないんです

No.21904　Re: KM図のat risk表示方法　　【赤岳】　2016/02/02(Tue) 10:07

青木先生，確かにEZRのスクリプトでは隠れたコードがあったようです。
以下に，EZRを経ずに書けるスクリプトにしました。
なお，ダミーデータも添付します。データは，"Dataset”というネーミングにしています。
ご指導よろしくお願いします。
＊＊＊＊
library(survival)
coxmodel <- coxph(Surv(year, event==1)~ age + skill + strata(TTT), data=Dataset, method="breslow")
cox <- survfit(coxmodel)
windows(width=7, height=7); par(lwd=1, las=1, family="sans", cex=1, mgp=c(2.5,1,0))
group.levels <- names(cox$strata[cox$strata>0])
mar <- par("mar")
mar[1] <- mar[1] + length(cox$strata) + 0.5
mar[2] <- mar[2] + 2
par(mar=mar)
opar <- par(mar = mar)
on.exit(par(opar))
plot(cox, col=1:32, lty=1, lwd=1, bty="l", mark.time=TRUE, xlim=c(0, 10), xlab="year", ylab="Survival Rate")
xticks <- axTicks(1)
n.atrisk <- nrisk(cox, xticks)
for (i in 1:length(cox$strata)){axis(1, at = xticks, labels = n.atrisk[i,], line=3+i, tick = FALSE)}
for (i in 1:length(cox$strata)){mtext(group.levels[i], at=-(xticks[2]-xticks[1])/2, side=1, line=4+i, cex=1)}
title(xlab = "Number at risk", line = 3.5, adj = 0)
title("Survival curve")
legend(locator(1), group.levels, title="TPC", col=1:32, lty=1, lwd=1, box.lty=0)
summary(cox)
remove(cox)
remove(coxmodel)

Dataset
subject TTT age event skill year
No1 S 51 1 a 0.09
No2 S 50 1 b 0.10
No5 S 47 1 b 0.73
No6 S 50 1 b 0.99
No7 S 71 1 c 1.15
No8 S 42 1 c 3.06
No9 S 47 1 a 4.42
No10 S 42 0 b 6.98
No11 S 32 0 c 7.46
No12 S 74 0 a 8.50
No13 S 56 0 a 8.59
No14 S 48 0 a 9.22
No15 S 50 0 a 9.53
No16 M 56 1 c 0.10
No17 M 57 0 b 1.19
No18 M 48 0 b 1.22
No19 M 42 0 b 1.24
No25 M 74 1 b 1.82
No26 M 39 0 c 1.92
No27 M 72 0 a 2.20
No28 M 45 0 a 2.22
No29 M 54 0 a 2.22
No30 M 41 0 a 2.60
No31 M 57 0 c 2.62
No32 M 41 0 b 2.93
No33 M 47 0 b 2.99
No34 M 68 0 b 3.09
No35 M 56 0 b 3.38
No36 M 41 0 b 3.52
No37 M 45 1 c 3.73
No38 M 47 0 c 4.14
No39 M 74 0 a 4.73
No40 M 80 0 b 4.87
No41 M 44 0 c 4.94
No42 M 50 0 a 4.98
No43 M 50 0 a 5.25
No44 M 66 0 a 5.45
No45 M 54 0 a 5.96
No46 M 48 0 c 6.46
No47 M 66 0 b 6.57
No48 M 50 0 b 6.97
No49 M 69 0 b 7.58

No.21905　Re: KM図のat risk表示方法　　【青木繁伸】　2016/02/02(Tue) 11:26

> n.atrisk <- nrisk(cox, xticks)
エラー: 関数 "nrisk" を見つけることができませんでした