- パッケージとデータの準備
- in_np( )を使ってモデルに当てはめずに生存曲線を描く
- in_sp( )を使ってセミパラメトリックモデルを当てはめる
- in_par( )を使ってパラメトリックモデルに当てはめる
- おわりに
- 参考資料
生存時間データにおいて正確な発生時点がわからないことを打ち切り（censoring）といい、次のような種類がある。
- Left censored：ある観察時点より前（=時間軸の左側）で発生したことは分かる
- Right censored：ある観察時点より後（=時間軸の右側）で発生したことは分かる
- Interval censored：ある観察区間内で発生したことは分かる
- Uncensored：正確な発生時点が分かる
区間打ち切りデータには2つのタイプがある。
1つ目はcurrent status data（あるいはcase I interval censored data）と呼ばれ、 例えばネズミに腫瘍が発生しているかを解剖して調べたデータなどがこれに相当する。 この場合は対象集団は、調査時点で既に腫瘍が発生している例（= left censoring）か、 調査時点では腫瘍が発生していなかった例（= right censoring）のいずれかに分類される。
もう1つはmixed case censoring（あるいはcase II interval censored data）と呼ばれるもので、 ある状態が発生しているかどうか定期的に観察を行うデータが相当する。 例えば定期検診で糖尿病と診断された場合、前回の検診から今回の検診の間のどこかで発症したことになるが、正確な時点は分からない。
アウトカムが「死亡」のように、発生の正確な時点*1が分かるものはright censoringを使った解析が行われるが、「糖尿病発症」や「腫瘍再発」などのように発生した区間しか分からない場合はinterval censoringとして扱った方が適切である。
パッケージとデータの準備
非発生例をright censoringとして扱う生存時間解析では、Surv(time, status)のように時間とイベント発生有無を組み合わせて生存時間データとして扱った。
Interval censoringを扱う場合は下のように区間の両端を指定する。Left censoring, right censoring, uncensoredも同じフレームワークで扱うことが可能。
ここではicenRegパッケージ内のmiceDataデータとIR_diabetesデータを使う。
miceData：肺腫瘍があるかどうかマウスを解剖して調べたcurrent status datal：観察期間の左端u：観察期間の右端grp：マウスの群分け（ce= conventional environment,ge= germ-free environment）
IR_diabetes：糖尿病性腎症の発症時期を調べたデータleft：観察期間の左端right：観察期間の右端gender：性別
library(icenReg) library(tidyverse)
> data(miceData) > head(miceData,3) l u grp 1 0 381 ce 2 0 477 ce 3 0 485 ce > tail(miceData,3) l u grp 142 913 Inf ge 143 942 Inf ge 144 986 Inf ge
後者は元々は区間打ち切りデータだけだが、上下位10%ずつをleft/right censoredに改変したデータdも作成しておく。
- Left censored dataにするとき：leftを
0にする - Right censored dataにするとき：rightを
InfまたはNAにする（"は不要）。 - Uncensored dataにするとき：leftとrightを同じ値にする
data(IR_diabetes) d <- IR_diabetes %>% mutate(left = replace(left, left < quantile(left,prob=0.1), 0)) %>% mutate(right = replace(right, right > quantile(right,prob=0.9), Inf))
> head(d,5) left right gender 1 24 Inf male 2 22 22 female 3 37 Inf male 4 20 20 male 5 0 16 male
症例5はleft censoring、症例1,3はright censoringになった。
in_np( )を使ってモデルに当てはめずに生存曲線を描く
ノンパラメトリック最尤推定法（non-parametric maximum likelihood estimator, NPMLE）を使ってモデルを仮定しないで打ち切り区間データの生存曲線を描く。
in_np( )を使う（interval censoring, non-parametric）。変数の指定の仕方は次のとおり。
- 区間打ち切りデータは
cbind(left, right)もしくはSurv(left, right, type="interval2")として渡す - 群分けしないときは
~0とする
ggplot系のggsurvplot()には対応していないみたい。
fit <- ic_np(cbind(l,u)~grp, data = miceData) plot(fit)
NPMLEで推定される生存曲線は1つに決まらない。表示される2つの曲線の間に入る曲線は全て尤度が等しくなる。
in_sp( )を使ってセミパラメトリックモデルを当てはめる
Coxモデルのようなセミパラメトリックモデルを打ち切り区間データに当てはめる。
in_sp( )を使う（interval censoring, semi-parametric）で次の引数を指定する。
model：比例ハザードモデルを使う場合は="ph"、比例オッズモデルを使う場合は="po"を指定bs_samples：信頼区間をブートストラップで推定するときのサンプリング回数
fit <- ic_sp(Surv(left, right, type="interval2") ~ gender, data = d, model = "ph", bs_samples = 500) summary(fit)
> summary(fit) Model: Cox PH Dependency structure assumed: Independence Baseline: semi-parametric Call: ic_sp(formula = Surv(left, right, type = "interval2") ~ gender, data = IR_diabetes, model = "ph", bs_samples = 500) Estimate Exp(Est) Std.Error z-value p gendermale -0.1402 0.8692 0.08802 -1.593 0.1111 final llk = -1964.96 Iterations = 38 Bootstrap Samples = 500
ハザード比はExp(Est)に表示されるが、信頼区間はEstimate±1.96Std.Errorで計算する必要あり。
複数コアで並列計算したいとき
doParallelパッケージのmakeCluster( )でコア数を指定した後、ic_sp( )で引数useMCores=TRUEを指定する。
library(doParallel) myCluster <- makeCluster(2) registerDoParallel(myCluster) fit <- ic_sp(Surv(left,right,type="interval2") ~ gender, data = d, model = "ph", useMCores=TRUE, bs_samples = 500)
in_par( )を使ってパラメトリックモデルに当てはめる
in_par( )を使う（interval censoring, parametric）。変数の指定の仕方は次のとおり。
model：使用するモデルを指定"ph"：比例ハザードモデル"po"：比例オッズモデル"aft"：加速死亡モデル（accelerated failure time model, AFT model）
dist：ベースラインに使用する分布関数を指定"exponential"：指数分布"gamma"：ガンマ分布"weibull"：ワイブル分布（デフォルト）"lnorm"：対数正規分布"loglogistic"：対数ロジスティック分布"generalgamma"：一般化ガンマ分布
fit <- ic_par(Surv(left, right, type="interval2") ~ gender, data = d, model = "ph", dist = "weibull")
グラフを描きたいときはnewdataを指定する必要があるみたい。
newdata <- data.frame(gender=c("male", "female")) rownames(newdata) <- c('males', 'females') plot(fit, newdata, lgdLocation = "topright")
おわりに
* 明確な発生時点が分からないtime-to-eventデータではinterval censoringを考慮してみる。
参考資料
- Anderson-Bergman C (2017). “icenReg: Regression Models for Interval Censored Data in R.” Journal of Statistical Software, 81(12), 1–23. doi: 10.18637/jss.v081.i12.
*1:正確な時点が分からなくても無視できるほど短い区間の区間打ち切りはuncensoredと扱ってさしつかえない。死亡時刻が分からなくても死亡日が分かっていればinterval censoringとして扱う必要はないということ。
