はじめに
ちょっと色々とやっていて R の survival パッケージの tt 関数の挙動が気になるという話になりました。
survival パッケージでは時間変動共変量の解析が実施可能であり、時間変動共変量と時間変動効果（回帰係数）についての vignette には詳細な解説が書かれています。
vignette 中のコードには色々な事例が書かれているのですが、tt 関数を使っているものがあります。
tt 関数は、時間変換関数 を組み込むための機能です。
なんでこんなものが必要になるかというと、モデル中に共変量と時間の交互作用項（例えば、$\beta(t)x=\beta_x x+\beta_{tx}\log(t)x$ のようなもの）が含まれる場合には、部分尤度関数を計算する際にリスク集合毎に分母も含めてイチから全計算をする必要があり、単純な Cox 回帰モデルの場合と計算方法が変わるためです（おまけでも触れていますが、時間交互作用項がある場合でもデータをイベント時間毎に層を持つ形に再構成して、層化統合解析をするのと同じ計算になります）。
モデル中に共変量と時間の交互作用項がないときは、分母の計算にはリスク集合に含まれる部分までの和 $\sum_{j \in R(t_i)} \exp(\mathbf{X}_j^T \boldsymbol{\beta})$ が必要なのでソートして累積和を1セット計算するだけでよいのですが、交互作用項がある場合は $\sum_{j \in R(t_i)} \exp\{\mathbf{X}_j(t_i)^T \boldsymbol{\beta}\}$ になって毎回全計算が必要です（これを自動化するために tt 関数が用意されているということかと思います）。
この違いは、自力で Cox 回帰の部分尤度関数の計算を実装してみるとよく分かるので、やってみると良いでしょう（簡単です）。
この tt 関数は、時間関数形を
coxph(Surv(time, status) ~ ph.ecog + tt(age), data = lung, tt = function(x, t, ...) pspline(x + t/365.25))
こんな感じに定義して使うのですが、関数の定義を書かずに
coxph(Surv(time, status) ~ ph.ecog + tt(age), data = lung)
のようにしても結果が出てきます。
このように定義を書かない場合には一体何が出てきているのか、という話題になり、これを確認した結果をまとめてみようかなと思いました。
tt 関数のデフォルトの挙動
これについてはヘルプファイルには書いていないです（確認した限り、そのような記載は見つけられませんでした）。
正解は vignette をよく読むと書いてありました。
There are other interesting uses for the time-transform capability. One example is O'Brien's logit-rank test procedure [6]. He proposed replacing the covariate at each event time with a logit transform of its ranks. This removes the influence of any outliers in the predictor x. For this case we ignore the event time argument and concentrate on the groupings to create the following tt function.
> function(x, t, riskset, weights){
obrien <- function(x) {
r <- rank(x)
(r-.5)/(.5+length(r)-r)
}
unlist(tapply(x, riskset, obrien))
}This relies on the fact that the input arguments to tt() are ordered by the event number or risk set. This function is used as a default if no tt argument is present in the coxph call, but there are tt terms in the model formula. (Doing so allowed me to depreciate the survobrien function).
ざっくりいうと「tt 関数は時間関数形の指定とは別の面白い使い方ができる、例えば、O'Brien's logit-rank testの実装にも使えて、これは coxph の tt 引数に何も指定しない場合の tt 関数のデフォルトになっている（これにより、私は survobrien 関数を非推奨にできた）」という感じです。
どうやら、tt 引数の指定をしないで実行すると、リスク集合毎に共変量の順位を変換した変数を生成するため、通常よくやられているような時間関数形を指定した解析とは目的が違う結果が出てくるということが分かります（Therneau 先生...、これはちょっとやめてほしいかも...）。
これで計算される回帰係数は共変量の順位に基づいた結果になっていて時間は無視されるので、時間関数形については評価できていません！
tt (time-transform) 関数のデフォルトが、全然 time-transform と関係がない仕様というのは若干思うところがあります。
通常の tt 関数の使い方
これも念のため紹介しておきます。
一番よくある使い方は、比例ハザード性の検定だったり、比較的簡単な構造を持つ時間変動効果のモデリングではないかと思います。
例えば、対数時間のスケールにおける時間変動効果 $\beta(t)x=\beta_x x+\beta_{tx}\log(t)x$ を持つモデルは tt 関数を使って当てはめることが可能です。
また、これと同じモデルにおいて、時間変動する部分の回帰係数のスコア検定を実施できる cox.zph という関数もあります。
fit1 <- coxph(Surv(time, status) ~ karno + tt(karno), data = veteran, tt = function(x, t, ...) x*log(t)) fit2 <- coxph(Surv(time, status) ~ karno, data = veteran) summary(fit1) cox.zph(fit2, transform = "log")
両者の実行結果は、
> summary(fit1) Call: coxph(formula = Surv(time, status) ~ karno + tt(karno), data = veteran, tt = function(x, t, ...) x * log(t)) n= 137, number of events= 128 coef exp(coef) se(coef) z Pr(>|z|) karno -0.083723 0.919686 0.016783 -4.989 6.08e-07 *** tt(karno) 0.013408 1.013498 0.004196 3.195 0.0014 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 exp(coef) exp(-coef) lower .95 upper .95 karno 0.9197 1.0873 0.8899 0.9504 tt(karno) 1.0135 0.9867 1.0052 1.0219 Concordance= 0.709 (se = 0.024 ) Likelihood ratio test= 52.99 on 2 df, p=3e-12 Wald test = 49.81 on 2 df, p=2e-11 Score (logrank) test = 56.97 on 2 df, p=4e-13 > cox.zph(fit2, transform = "log") chisq df p karno 10.5 1 0.0012 GLOBAL 10.5 1 0.0012
となりました。
tt 関数を使った検定は Wald 検定（$P = 0.0014$）で、zph 検定はスコア検定（$P = 0.0012$）のため若干値は異なります。
参考文献
- Therneau T, Crowson C, Atkinson E. Using time dependent covariates and time dependent coefficients in the Cox model. [Link]
- O'Brien P. A nonparametric test for association with censored data. Biometrics 1978;34:243–250. DOI: 10.2307/2530014.
おまけ
一応 tt 関数の中身を見てみて、何を計算しているのかを確認しました。
library(survival) temp <- subset(pbc, id <= 312, select = c(id:sex, stage)) pbc2 <- tmerge(temp, temp, id = id, death = event(time, status)) pbc2 <- tmerge(pbc2, pbcseq, id = id, lbili = tdc(day, log(bili))) cox_model0 <- coxph(Surv(tstart, tstop, death == 2) ~ lbili, pbc2) cox_model1 <- coxph(Surv(tstart, tstop, death == 2) ~ tt(lbili), pbc2) tt <- function(x, riskset) { obrien <- function(x) { r <- rank(x) (r - 0.5)/(0.5 + length(r) - r) } unlist(tapply(x, riskset, obrien)) } Y <- cbind(pbc2$tstart, pbc2$tstop, as.numeric(pbc2$death == 2)) sort.end <- order(-Y[, 2], Y[, 3]) sort.start <- order(-Y[, 1]) newstrat <- c(1L, rep(0, nrow(Y) - 1)) counts <- .Call(survival:::Ccoxcount2, Y, as.integer(sort.start - 1L), as.integer(sort.end - 1L), as.integer(newstrat)) tindex <- counts$index istrat <- rep(1:length(counts$nrisk), counts$nrisk) Y <- Surv(rep(counts$time, counts$nrisk), counts$status) tv <- pbc2$lbili[tindex] ttv <- tt(tv, istrat) pbc3 <- data.frame(time = Y[, 1], status = counts$status, ttv, istrat) cox_modelttv <- coxph(Surv(time, status) ~ ttv + strata(istrat), data = pbc3) cox_model1 cox_modelttv
これが tt 関数をそのまま実行した結果で、
> cox_model1 Call: coxph(formula = Surv(tstart, tstop, death == 2) ~ tt(lbili), data = pbc2) coef exp(coef) se(coef) z p tt(lbili) 0.0073115 1.0073383 0.0005712 12.8 <2e-16 Likelihood ratio test=74.29 on 1 df, p=< 2.2e-16 n= 1807, number of events= 125
こっちは tt 関数が作っているデータを確認するために検算した方です。
> cox_modelttv Call: coxph(formula = Surv(time, status) ~ ttv + strata(istrat), data = pbc3) coef exp(coef) se(coef) z p ttv 0.0073115 1.0073383 0.0005712 12.8 <2e-16 Likelihood ratio test=74.29 on 1 df, p=< 2.2e-16 n= 24422, number of events= 125
データ数以外は一致していることが分かります。
検算用のデータ（pbc3）を見ると、リスク集合毎に層が作られていることが確認できます。
> pbc3 time status ttv istrat 1 4191 0 1.36363636 1 2 4191 0 0.36842105 1 3 4191 0 0.23809524 1 4 4191 0 2.71428571 1 5 4191 0 0.08333333 1 6 4191 0 0.62500000 1 7 4191 0 0.62500000 1 8 4191 0 1.88888889 1 9 4191 0 4.20000000 1 10 4191 0 1.00000000 1 11 4191 0 0.08333333 1 12 4191 0 25.00000000 1 13 4191 1 7.66666667 1 14 4079 0 0.78947368 2 15 4079 0 0.13333333 2 16 4079 0 0.54545455 2 17 4079 0 3.85714286 2 18 4079 0 0.30769231 2 19 4079 0 1.00000000 2 20 4079 0 1.61538462 2 21 4079 0 2.09090909 2 22 4079 0 5.80000000 2 23 4079 0 1.26666667 2 24 4079 0 0.30769231 2 25 4079 0 10.33333333 2 26 4079 0 0.54545455 2 27 4079 0 0.03030303 2 28 4079 0 0.13333333 2 29 4079 0 2.77777778 2 30 4079 1 33.00000000 2 31 3853 0 0.61290323 3
