R で競合リスク・モデル
太郎丸 博
2018年11月06日
1 Introduction
Paul D. Allison 2014, Event History and Survival Analysis Sage の Chapter 5 のデータを使い、その分析結果を再現する。
1.1 Downloading Sample Data
https://statisticalhorizons.com/resources/data-sets の “tarp.dta” をダウンロードして、ワーキング・ディレクトリに置いておく。
1.2 Reading Stata Data
# setwd("C:/Users/taroh/Dropbox/18HI_Course/SampleDataAnalyses") # Working Directory の指定
library(foreign) # Stata のデータを読み込むためのパッケージ
tarp.d0 <- read.dta("tarp.dta") # read.dta() で dta 形式のファイルを読み込める
names(tarp.d0) # 使わない変数がたくさんあるので、使う変数だけのデータを作る## [1] "id" "age" "race" "sex" "mar" "ed" "prvarrst"
## [8] "curcrim" "crime" "yrsent" "depart" "relyr" "relmo" "expgrp"
## [15] "sexint" "paro" "higrade" "age1st" "numarst" "numconv" "emp1"
## [22] "emp2" "emp3" "emp4" "emp5" "emp6" "emp7" "emp8"
## [29] "emp9" "emp10" "emp11" "emp12" "emp13" "emp19" "frstarrmo"
## [36] "frstarrday" "frstarryr" "frstcharge" "arr2mo" "arr2dy" "arr2yr" "arr3mo"
## [43] "arr3dy" "arr3yr" "arr4mo" "arr4dy" "arr4yr" "arr5mo" "arr5dy"
## [50] "arr5yr" "arr6mo" "arr6dy" "arr6yr" "crimper" "crimprop" "crimoth"
## [57] "numper" "numprop" "numoth" "white" "arstdate" "reldate" "arrstday"
## [64] "arrest" "paroled" "married" "edcomb" "fin" "male" "type"
## [71] "charge" "arrstcount"tarp.d1 <- tarp.d0[, c("arrstday", "type", "fin", "age", "white", "male",
"married", "paro", "numprop", "crimprop","numarst", "edcomb")]
summary(tarp.d1)## arrstday type fin age white male
## Min. : 1.0 Min. :0.0000 Min. :0.0000 Min. :16.6 Min. :0.0000 Min. :0.0000
## 1st Qu.:223.0 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:21.5 1st Qu.:0.0000 1st Qu.:1.0000
## Median :365.0 Median :0.0000 Median :1.0000 Median :24.9 Median :0.0000 Median :1.0000
## Mean :294.3 Mean :0.5697 Mean :0.5912 Mean :27.8 Mean :0.4195 Mean :0.9431
## 3rd Qu.:365.0 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:31.0 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :365.0 Max. :2.0000 Max. :1.0000 Max. :71.8 Max. :1.0000 Max. :1.0000
## married paro numprop crimprop numarst
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. : 1.000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 2.000
## Median :0.0000 Median :0.0000 Median :0.0000 Median :1.0000 Median : 4.000
## Mean :0.3069 Mean :0.4431 Mean :0.2758 Mean :0.6685 Mean : 7.441
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.: 8.000
## Max. :1.0000 Max. :1.0000 Max. :5.0000 Max. :1.0000 Max. :98.000
## edcomb
## Min. : 0.000
## 1st Qu.: 8.000
## Median :10.000
## Mean : 9.458
## 3rd Qu.:11.000
## Max. :18.0002 Aalen-Johnson Estimation
通常のカプランマイヤー推定は複数のタイプのイベントを想定したものではないが、survival パッケージの survfit() を使えば、複数のタイプのイベントでも推定できる。これは Aalen-Johnson 推定と呼ばれるもので、Kaplan Meyer 推定や後の Cummulative Incidence 推定を一般化したものと位置づけられている（競合リスクモデルだけでなく、もっと複雑な移動パターンも扱えるらしい）。 Therneau Crowson and Atkinson 2018 によると、以下のように定義できる。今、複数の状態 \(i=1, \ ,2, \, \ldots k\) がある。例えば、Allison (2014) の例では、逮捕されていない状態、財産犯罪で逮捕された状態、財産犯罪以外の罪で逮捕された状態、の3つがある。\(p (t) = (p_1(t), \ p_2 (t), \, \ldots p_k (t) )\) は \(t\)時点において、それぞれの状態にある確率を示すベクトルである。さらに移行行列 \(T(t)\) は \(k\) 行 \(k\) 列の行列でその \(i\) 行 \(j\) 列目の要素は \(\lambda_{it}(t)\) と表記し、\(t\) 時点において状態 \(i\) から \(j\) に移行した人の比率である。最初にどの状態にいるかの確率（あるいは比率でもよいと思うが）を \(p(0)\)とすると、 \[ p(t) = p(0)\prod_{s \leq t}T(s) \qquad (1) \] と推定される。Allison (2014) の例では、釈放された時点では全員逮捕されていないので、\(p(0) = (1, 0, 0)\) である。\(t=1\)のときに釈放された932人のうち財産犯で二人、財産犯以外で一人逮捕されたので、 \[ T(1) = \begin{pmatrix} 929/932 & 2/932 & 1/932 \\ 0 & 1& 0 \\ 0 & 0& 1 \end{pmatrix} = \begin{pmatrix} .997 & .002 & .001 \\ 0 & 1& 0 \\ 0 & 0& 1 \end{pmatrix} \] である。実際には逮捕された後、死刑や終身刑にならない限りいつかは釈放されるが、Allison (2014) では釈放を無視して単純な競合リスクモデルを仮定しているので、2行目と3行目は常に主対角線上が 1 でそれ以外は 0 に必ずなる。
これらを (1) 式に代入すると \[ p (1) = p(0) T(1) = (1,\ 0, \ 0) \begin{pmatrix} .997 & .002 & .001 \\ 0 & 1& 0 \\ 0 & 0& 1 \end{pmatrix} = (.997, .002, .001) \] である。同様にして、\(t=2\) では、逮捕されていない 929人のうち一人が財産犯で逮捕されたので、 \[ T(2) = \begin{pmatrix} 928/929 & 1/929 & 0 \\ 0 & 1& 0 \\ 0 & 0& 1 \end{pmatrix} = \begin{pmatrix} .999 & .001 & 0 \\ 0 & 1& 0 \\ 0 & 0& 1 \end{pmatrix} \] である。さらに同様にして \[ p (2) = p(0) T(1) T (2) = (1,\ 0, \ 0) \begin{pmatrix} .997 & .002 & .001 \\ 0 & 1& 0 \\ 0 & 0& 1 \end{pmatrix} \begin{pmatrix} .999 & .001 & 0 \\ 0 & 1& 0 \\ 0 & 0& 1 \end{pmatrix} = (.997, .002, .001) \begin{pmatrix} .999 & .001 & 0 \\ 0 & 1& 0 \\ 0 & 0& 1 \end{pmatrix} = = (.996, .003, .001) \] である。以下同様にしてフォローアップ期間の最後まで同じ計算を続ければよいだけである。
R で Aalen-Johnson 推定をする場合、イベントに二種類のタイプがあるとすると、 Surv() 関数ではイベントの生起を示すダミー変数ではなく、“センサー”、“タイプ1のイベント発生”、“タイプ2のイベント発生” の3つのカテゴリからなる因子を指定する。このときセンサーが最初のカテゴリでなければならない。
tarp.d1$type <- factor(tarp.d1$type, labels = c("no_arrest", "non_property", "property")) # 何故か 1が non - property
summary(tarp.d1$type) ## no_arrest non_property property
## 598 137 197tarp.d1$type <- factor(tarp.d1$type, levels = c("no_arrest", "property", "non_property")) # 順番をテキストに合わせる
summary(tarp.d1$type)## no_arrest property non_property
## 598 197 137library(survival)
km0 <- survfit(Surv(arrstday, type) ~ 1, data = tarp.d1)
summary(km0, c(1, 50, 100, 150, 200, 250, 300, 350, 365)) # 出力が長くなりすぎるので、特定の時点を指定して出力## Call: survfit(formula = Surv(arrstday, type) ~ 1, data = tarp.d1)
##
## time n.risk n.event P(property) P(non_property) P()
## 1 932 3 0.00215 0.00107 0.997
## 50 892 38 0.03326 0.01073 0.956
## 100 827 64 0.06867 0.04399 0.887
## 150 770 60 0.10730 0.06974 0.823
## 200 718 51 0.14485 0.08691 0.768
## 250 682 35 0.16309 0.10622 0.731
## 300 646 35 0.18348 0.12339 0.693
## 350 611 35 0.20172 0.14270 0.656
## 365 607 13 0.21137 0.14700 0.642plot(km0, col = 2:3, lty = 1:2)
legend("topleft", col = 2:3, lty = 1: 2, legend = c("property", "non_property"))
競合リスク・モデルの場合は、生存率だけでなく、それぞれのイベントの累積発生確率（その時点までにイベントが発生する確率）も計算される。plot() ではデフォルトでは生存関数ではなく、それぞれのイベントの累積発生確率がプロットされる。これはテキストとは順番が前後するが、p. 65 の Figure 5.1 に対応する。ただし、同じ推定法なのかどうかは未確認。
もちろんグループ別に推定することもできる。
km1 <- update(km0, ~ age < 25)
summary(km1, 1 + 50 * (1:7) )## Call: survfit(formula = Surv(arrstday, type) ~ age < 25, data = tarp.d1)
##
## age < 25=FALSE
## time n.risk n.event P(property) P(non_property) P()
## 51 446 19 0.0323 0.00862 0.959
## 101 415 30 0.0625 0.04310 0.894
## 151 394 21 0.0905 0.06034 0.849
## 201 373 21 0.1164 0.07974 0.804
## 251 354 19 0.1358 0.10129 0.763
## 301 340 14 0.1487 0.11853 0.733
## 351 325 15 0.1595 0.14009 0.700
##
## age < 25=TRUE
## time n.risk n.event P(property) P(non_property) P()
## 51 445 23 0.0363 0.0128 0.951
## 101 412 33 0.0748 0.0449 0.880
## 151 373 40 0.1261 0.0791 0.795
## 201 343 29 0.1731 0.0940 0.733
## 251 327 16 0.1902 0.1111 0.699
## 301 306 21 0.2179 0.1282 0.654
## 351 286 20 0.2436 0.1453 0.611plot(km1, col= 2 : 3, lty = c(1, 1, 2, 2))
legend("topleft", col = 2:3, lty = c(1, 1, 2, 2), legend = c("non_property, age < 25", "property, age < 25", "non_property, age => 25", "property, age => 25"))
3 Cox Regression
競合リスク・モデルを Cox 回帰で推定する場合、イベントのタイプごとに Cox 回帰モデルを推定すればよい。このとき他のタイプのイベントはセンサーとみなす。それゆえ、財産犯での逮捕ハザードを従属変数とする場合は、その他の犯罪で逮捕されたケースは、逮捕された時点でセンサーされたとみなす。以下はテキスト Table 5.1 の推定。
cox0 <- vector("list", 3) # 3回コックス回帰分析して表にまとめるので、結果を収納するリストを用意
names(cox0) <- c("All Arests", "Property Arrests", "Non-Property Arrests") # 名前を付ける
library(survival)
cox0[[1]] <- coxph(Surv(arrstday, type != "no_arrest") ~ . , data = tarp.d1)
cox0[[2]] <- coxph(Surv(arrstday, type == "property") ~ . , data = tarp.d1)
cox0[[3]] <- coxph(Surv(arrstday, type == "non_property") ~ . , data = tarp.d1)
library(texreg)## Version: 1.36.23
## Date: 2017-03-03
## Author: Philip Leifeld (University of Glasgow)
##
## Please cite the JSS article in your publications -- see citation("texreg").screenreg(cox0, digits = 3, single.row = TRUE)##
## =============================================================================
## All Arests Property Arrests Non-Property Arrests
## -----------------------------------------------------------------------------
## fin 0.121 (0.114) 0.201 (0.149) 0.007 (0.176)
## age -0.034 (0.009) *** -0.041 (0.012) *** -0.027 (0.012) *
## white -0.241 (0.118) * -0.353 (0.156) * -0.077 (0.182)
## male 0.502 (0.309) 0.111 (0.346) 1.388 (0.715)
## married -0.222 (0.130) -0.332 (0.177) -0.073 (0.193)
## paro -0.211 (0.118) -0.147 (0.154) -0.302 (0.186)
## numprop 0.310 (0.072) *** 0.319 (0.097) ** 0.308 (0.107) **
## crimprop 0.425 (0.137) ** 0.883 (0.206) *** -0.069 (0.192)
## numarst 0.018 (0.005) *** 0.019 (0.006) ** 0.016 (0.007) *
## edcomb -0.068 (0.025) ** -0.053 (0.034) -0.083 (0.038) *
## -----------------------------------------------------------------------------
## AIC 4365.169 2556.328 1805.483
## R^2 0.087 0.082 0.031
## Max. R^2 0.991 0.940 0.857
## Num. events 334 197 137
## Num. obs. 932 932 932
## Missings 0 0 0
## PH test 0.001 0.011 0.376
## =============================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.053.1 Tests
サンプルサイズが十分に大きければ Cox 回帰の係数の z 値は正規分布に近似する。それゆえ、二つの係数が独立に分布するならば、テキスト60頁の式 (5.3) のような計算で、二つの係数が等しいかどうか検定できる。以下は、property arrests と non-property arrests の係数の差をすべて検定する計算。
coef.prop <- summary(cox0[[2]])$coefficients[, c(1, 3)] # Cox回帰の係数とSE だけとってくる
coef.nonprop <- summary(cox0[[3]])$coefficients[, c(1, 3)]
# 係数の差の検定統計量 z
z <- (coef.prop[, 1] - coef.nonprop[, 1]) /
sqrt(coef.prop[, 2] ^ 2 + coef.nonprop[, 2] ^ 2)
p <- 2 * (1 - pnorm(abs(z)) ) # 両側の p 値
kable(cbind(z, p), digits = 3)| z | p | |
|---|---|---|
| fin | 0.840 | 0.401 |
| age | -0.829 | 0.407 |
| white | -1.151 | 0.250 |
| male | -1.608 | 0.108 |
| married | -0.990 | 0.322 |
| paro | 0.642 | 0.521 |
| numprop | 0.072 | 0.943 |
| crimprop | 3.380 | 0.001 |
| numarst | 0.253 | 0.801 |
| edcomb | 0.574 | 0.566 |
crimprop の z 値がテキストと少し違うのだが、これは Property Arrests の標準誤差が異なる (テキストは .296 だが、coxph()の計算結果は .206 である) のが原因である。なぜこのような違いが生じたのかは不明。Allison の転記ミスか?
続いて尤度比検定。
(ll.cox0 <- sapply(cox0, logLik) ) # sapply() は apply のリスト版。logLik() はモデルの対数尤度を返す## All Arests Property Arrests Non-Property Arrests
## -2172.5847 -1268.1640 -892.7417(
dev.cox0 <- -2 * (
ll.cox0[1] - sum(ll.cox0[2 : 3])
)
)## All Arests
## 23.357891 - pchisq(dev.cox0, df= 10)## All Arests
## 0.009499923微妙にテキストと値が異なるが、許容範囲だろう。
4 Cumulative Incidence Functions
イベントのタイプ別に、累積確率を計算するだけならば、すでに Aalen-Johnson Estimation 推定の節で説明した。
Allison (2014) の Table 5.2 のような分析をするためには、cmprsk パッケージの crr() 関数を使う。
crr(ftime = フォローアップ時間（ベクトル）, fstatus = イベントのタイプ（因子）, cov1 = 独立変数行列, failcode = 推定したいイベントの文字例, cencode = センサーの文字列)
で、通常の回帰分析とはだいぶ書き方が違う。最も大きなポイントは data = 引数 が使えない点と、 ~ x1 + x2 といった式を使えない点である。以下が Table 5.2 の2列目と3列目のモデルの推定結果の再現。
covariates0 <- model.matrix(~., data = tarp.d1[, -1:-2]) # model.matrix() はモデル式とデータから、因子をダミー変数に変換し、切片や連続変数を含めた共変量の行列を返す
covariates0 <- covariates0[, -1] # 最初の列は切片なので、取り除く
library(cmprsk)
cif01 <- crr(ftime = tarp.d1$arrstday, # 財産犯での逮捕の推定
fstatus = tarp.d1$type,
cov1 = covariates0,
failcode = "property",
cencode = "no_arrest")
summary(cif01)## Competing Risks Regression
##
## Call:
## crr(ftime = tarp.d1$arrstday, fstatus = tarp.d1$type, cov1 = covariates0,
## failcode = "property", cencode = "no_arrest")
##
## coef exp(coef) se(coef) z p-value
## fin 0.2358 1.266 0.15187 1.553 1.2e-01
## age -0.0404 0.960 0.01366 -2.960 3.1e-03
## white -0.3514 0.704 0.16165 -2.174 3.0e-02
## male 0.0426 1.043 0.37151 0.115 9.1e-01
## married -0.3076 0.735 0.17607 -1.747 8.1e-02
## paro -0.1079 0.898 0.15386 -0.701 4.8e-01
## numprop 0.2755 1.317 0.10215 2.697 7.0e-03
## crimprop 0.8875 2.429 0.20837 4.259 2.1e-05
## numarst 0.0151 1.015 0.00697 2.165 3.0e-02
## edcomb -0.0490 0.952 0.03738 -1.311 1.9e-01
##
## exp(coef) exp(-coef) 2.5% 97.5%
## fin 1.266 0.790 0.940 1.705
## age 0.960 1.041 0.935 0.986
## white 0.704 1.421 0.513 0.966
## male 1.043 0.958 0.504 2.161
## married 0.735 1.360 0.521 1.038
## paro 0.898 1.114 0.664 1.214
## numprop 1.317 0.759 1.078 1.609
## crimprop 2.429 0.412 1.615 3.654
## numarst 1.015 0.985 1.001 1.029
## edcomb 0.952 1.050 0.885 1.025
##
## Num. cases = 932
## Pseudo Log-likelihood = -1288
## Pseudo likelihood ratio test = 73 on 10 df,cif02 <- crr(ftime = tarp.d1$arrstday, # 財産犯での逮捕の推定
fstatus = tarp.d1$type,
cov1 = covariates0,
failcode = "non_property",
cencode = "no_arrest")
summary(cif02)## Competing Risks Regression
##
## Call:
## crr(ftime = tarp.d1$arrstday, fstatus = tarp.d1$type, cov1 = covariates0,
## failcode = "non_property", cencode = "no_arrest")
##
## coef exp(coef) se(coef) z p-value
## fin 0.00366 1.004 0.17735 0.0206 0.980
## age -0.02185 0.978 0.01131 -1.9318 0.053
## white -0.01028 0.990 0.17374 -0.0592 0.950
## male 1.43598 4.204 0.71293 2.0142 0.044
## married -0.03460 0.966 0.19997 -0.1731 0.860
## paro -0.27939 0.756 0.17861 -1.5643 0.120
## numprop 0.25794 1.294 0.10233 2.5207 0.012
## crimprop -0.18365 0.832 0.19203 -0.9564 0.340
## numarst 0.01133 1.011 0.00678 1.6699 0.095
## edcomb -0.07513 0.928 0.03699 -2.0312 0.042
##
## exp(coef) exp(-coef) 2.5% 97.5%
## fin 1.004 0.996 0.709 1.421
## age 0.978 1.022 0.957 1.000
## white 0.990 1.010 0.704 1.391
## male 4.204 0.238 1.039 17.001
## married 0.966 1.035 0.653 1.429
## paro 0.756 1.322 0.533 1.073
## numprop 1.294 0.773 1.059 1.582
## crimprop 0.832 1.202 0.571 1.213
## numarst 1.011 0.989 0.998 1.025
## edcomb 0.928 1.078 0.863 0.997
##
## Num. cases = 932
## Pseudo Log-likelihood = -915
## Pseudo likelihood ratio test = 23 on 10 df,以下は Figure 5.2 の再現。
preddata <- colMeans(covariates0) # 共変量の平均値
plot(
predict(cif02, preddata),
xlab = "Days since Release",
ylab = "Cumulative Incidence"
)