二値の時間依存性変数（一度変化したら一定）を考慮した生存時間の生成

ことの始まり

生存時間解析だけでなく、longitudinal dataの解析や因果推論など、いろいろな場面で時間依存性の変数を持つデータセットを用いる機会は増えている。

数理統計苦手マンとしては、サンプルデータを生成できたら事前に色々な評価ができて便利なので、時間依存性の変数を持つデータの生成ができるようになりたい。

また、そのような生成関数を作っておいて、諸条件を変えていけるような便利コードを作りたい。

生存時間の生成

下記の要領で、シンプルな事例からトライした。

指数分布を使う。指数分布だとrexp()があるので、ハザード関数の逆関数を用いて生存時間を得る必要はない。
時間依存性変数X3は二値変数で、どこかで1となったらずっと1。
簡単のため時間依存変数以外の変数は独立した2つの標準正規分布。
マトリックス表記でないので変数増やしたい時は行列形式にする。
365日目で打ち切りとしている。ランダムな打ち切りに変える／追加することもできると思う。
実際のデータだと生存時間はintegerなのでfloor関数を使っている。
- ChatGPTに見せるとこの点は不要ではないかといつも指摘される。
X3の変化点までの時間とそこからイベントが起きるまでの時間を足しているのは、指数分布のmemoryless propertyによるものかもしれない？
X3が変化する時点の候補から365日目を抜かないとstarttime = endtime = 365 となってしまうことがあるので調整している。max(1, hoge)なども同様の理由で、integerにしている関係で生成された時間が0になってstartime = endtime = 0となることを避けるための条件。
X3は時間依存性変数としてCoxモデルに入れたいので計数過程形式のデータにしている。
ベースとなるコードはChatGPTに言って書いてもらった。
- なので自分ならこうは書かないというものもある。逆に参考になる書き方もあるのでいい勉強になる。

library(survival)
library(tidyverse)

# Prep --------------------------------------------------
# Set parameters
set.seed(123)  # For reproducibility
N <- 500        # Number of subjects
T_max <- 365   # Maximum follow-up time in days
true_beta <- c(0.2, -0.3, 0.5)  # True coefficients for the covariates (log hazard)

# Approach using rexp() - binary static TD variable ----------------------------
## data generating function -------
gen_data <-
  function() {
    # Generate static covariates
    X1 <- rnorm(N)  # Covariate 1
    X2 <- rnorm(N)  # Covariate 2
    
    # Initialize dataset
    data <- data.frame(id=integer(0), starttime=double(0), endtime=double(0),
                       X1=double(0), X2=double(0), X3=integer(0), status=integer(0))
    
    # Hazard function and survival time simulation
    lambda_base <- 0.01  # Base hazard rate
    
    for (i in 1:N) {
      # Random time point to switch the time-varying covariate from 0 to 1
      # T_max - 1 is to ensure after the switch there is at least one day follow-up,
      # preventing starttime = endtime
      switch_time <- sample(1:(T_max-1), 1)
      
      # Hazard rate for this subject before the switch
      hazard_rate_pre_switch <- lambda_base * exp(true_beta[1] * X1[i] + true_beta[2] * X2[i])
      
      # Simulate survival time before the switch
      survival_time_pre_switch <- max(1,floor(rexp(1, hazard_rate_pre_switch)))
      
      if (survival_time_pre_switch < switch_time) {
        # Subject experiences the event before the switch
        data <- bind_rows(data, c(id=i, starttime=0, endtime=survival_time_pre_switch, X1=X1[i], X2=X2[i], X3=0, status=1))
      } else {
        # Subject does not experience the event before the switch
        data <- bind_rows(data, c(id=i, starttime=0, endtime=switch_time, X1=X1[i], X2=X2[i], X3=0, status=0))
        
        # Hazard rate for this subject after the switch
        hazard_rate_post_switch <- lambda_base * exp(true_beta[1] * X1[i] + true_beta[2] * X2[i] + true_beta[3])
        
        # Simulate survival time after the switch
        survival_time_post_switch <- max(1,floor(rexp(1, hazard_rate_post_switch)))
        
        # Determine end time and status after the switch
        end_time_post_switch <- min(switch_time + survival_time_post_switch, T_max) 
        status_post_switch <- as.integer(switch_time + survival_time_post_switch <= T_max)
        
        # Add data for the period after the switch
        data <- bind_rows(data, c(id=i, starttime=switch_time, endtime=end_time_post_switch, X1=X1[i], X2=X2[i], X3=1, status=status_post_switch))
      }
      
    }
    
    return(data)
    
  }

データの中身はこんな感じ。

> sample_data <- gen_data()
> head(sample_data)
  id starttime endtime         X1          X2 X3 status
1  1         0      14 -1.1998415 -0.08959434  0      1
2  2         0      79  1.8164182 -0.28290468  0      1
3  3         0      71  0.2776409  2.16411666  0      1
4  4         0      62 -1.1545589 -0.90907695  0      0
5  4        62      90 -1.1545589 -0.90907695  1      1
6  5         0      46  0.4880977 -0.56249520  0      1

データ生成関数の検証

バイアスとカバレッジを確認する。
このやり方は{simsurv}のvignetteを参考にした。

sim <-
  function() {

    data <- gen_data()
    
    fit <- coxph(Surv(starttime, endtime, status) ~ X1+X2+X3, data)
    summary(fit)
    
    est_cox <- fit$coefficients
    ses_cox <- sqrt(diag(fit$var))
    cil_cox <- est_cox + qnorm(.025) * ses_cox
    ciu_cox <- est_cox + qnorm(.975) * ses_cox
    
    c(bias_cox = est_cox - true_beta, 
      coverage_cox = ((true_beta > cil_cox) & (true_beta < ciu_cox))
    )
    
  }

# Get the means of replications
sim_res <-rowMeans(replicate(100, sim()))

実行結果はこんな感じ。多分うまくできている。

> sim_res
    bias_cox.X1     bias_cox.X2     bias_cox.X3 coverage_cox.X1 coverage_cox.X2 coverage_cox.X3 
  -0.0015067926   -0.0003750131    0.0058340024    0.9400000000    0.9800000000    0.9500000000

ironman

逆変換を使って前と同じことをやってみる

やること

次は同じことをAustin (2012)の式1から生成する。

この文献では乱数生成の一つである逆変換を用いている。「生存時間を表す確率変数T、その累積密度関数を $S$ としたときに、 $S(T)$ が一様分布に従い（S(t)ではなくS(T)であることがポイント）、 $F(T)=1-S(T)$ も一様分布に従うこと」を利用する。ここで $F(T)$ がイベント発生時間の関数であるので、これからデータを生成する。

式で書くと：

U = exp(−H_{0}(T) \cdot e^{\beta' X}) \sim U[0,1]

T = H_{0}^{−1} (\cfrac {−log(U )} {e^{\beta' X}})

となる。

仮定している生存時間によってハザードの逆関数の形は変わってくるが、指数分布の場合：

T = \cfrac {−log(U)} {\lambda \cdot e^{\beta' X} }

となる。時間依存性共変量が1になる前はこの式からのサンプリングになる。
なお、ある確率変数にその分布関数を適用すると一様分布に変換されることをProbability integral transformと言い、帰無仮説が成り立つときのP値が一様分布であることや、回帰モデルの診断プロットであるQQプロットなどにも関連する。

データ生成

上のコードとは違ってfloorではなくceilingを使ったり、ちょっと変わっているところがあるが気にしない。

[2025-05-05 追記] 以下のコードの間違いを修正し、正しく動くようになった。

修正点1
時間依存性バイナリ変数が1になった以降の生存時間を計算する部分で、 $h_{post} * t_0$ が余分に足されてしまっていたのを削除した。

# Event occurs after t0
        survival_time_post_t0 <- (-log(u) - hazard_pre_t0*t0 + hazard_post_t0*t0) / hazard_post_t0

<after>

# Event occurs after t0
        raw_time_post_t0 <-  (-log(u) - hazard_pre_t0*t0 ) / hazard_post_t0

修正点2
survival_time_post_t0は $t_0$ 以降の生存時間なので、最終的に $T_{max}$ と比較するのは $t_0$ を足し込んだ時間になるので、survival_time_post_t0 + t0に変更した。

survival_time_post_t0 <- min(survival_time_post_t0, T_max)

<after>

        # survival_time_post_t0 is survival time after t0 so t0 is added to get a new overall survival T (Fixed)
        survival_time_post_t0 <- ceiling(min(raw_time_post_t0 + t0, T_max))

補足
chatGPTがどんどん賢くなっているのでコードを投げて色々と確認をした結果の修正。修正作業をしながら自分で整理したこととして、この共変量のスイッチ前後で異なるハザードを用いて生存時間を生成する方法は、二段階サンプリングとして行うのが妥当だということ。なぜなら、最初の生存時間の生成u <- runif(1)のタイミングでは時間依存性共変量X3=0の場合の生存解析なので、この生存時間が最終的な生存時間にはならないから。なので、X3のスイッチ（０→1）が起こる時間 $t_0$ よりも生存時間が長いかどうかをまず一段階目のサンプリングで決めて、その後もう一度一様分布からuをサンプリングしてそれとX3の回帰係数を考慮した $h_{post}$ とでスイッチ後の生存時間を生成し、 $t_0$ と足し合わせるという方法になる。

これを1段階目のサンプリングのみで行うこともできて、それが今回の実装方法。具体的には、

T = \frac{-log(u)}{h_{pre}} \Rarr h_{pre} = \frac{-log(u)}{T}

なので、 $T$ までの累積ハザードは $-log(u)$ ということになる。これを利用すると、この累積ハザードから $t_0$ までの累積ハザードを引いた部分が $t_0$ 以降の累積ハザードということになるので、その累積ハザードを $h_{post}$ で除してあげれば $T_1$ = $T - t_0$ ということになる。そして、この $T_1 + t_0$ とすることで生存時間を更新して最終的な生存時間の生成ができる。

コードの全貌

# Set parameters
set.seed(123)  # For reproducibility
N <- 500       # Number of subjects
T_max <- 365  # Maximum time in days
lambda <- 0.01  # Base hazard rate
true_beta <- c(0.2, -0.3, 0.5)  # Coefficients for covariates

## data generating function -------
gen_data <-
  function() {
    
    # Generate static covariates
    X1 <- rnorm(N)  # Covariate 1 (standard normal)
    X2 <- rnorm(N)  # Covariate 2 (standard normal)
    
    # Initialize dataset
    data <- data.frame(id=integer(0), starttime=double(0), endtime=double(0),
                       X1=double(0), X2=double(0), X3=integer(0), status=integer(0))
    
    # Simulate survival times
    for (i in 1:N) {
      t0 <- sample(1:(T_max-1), 1)  # Random time point for covariate change
      
      # Calculate hazard before and after t0
      hazard_pre_t0 <- lambda * exp(true_beta[1] * X1[i] + true_beta[2] * X2[i])
      hazard_post_t0 <- lambda * exp(true_beta[1] * X1[i] + true_beta[2] * X2[i] + true_beta[3])
      
      # Simulate survival time
      u <- runif(1)
      if (-log(u) < hazard_pre_t0 * t0) {
        # Event occurs before t0
        raw_time_pre_t0 <- -log(u) / hazard_pre_t0
        status_pre_t0 <- as.integer(raw_time_pre_t0 <= T_max)
        survival_time_pre_t0 <- ceiling(min(raw_time_pre_t0, T_max))
        data <- bind_rows(data, c(id=i, starttime=0, endtime=survival_time_pre_t0, X1=X1[i], X2=X2[i], X3=0, status=status_pre_t0))

      } else {
        # Event occurs after t0 (Fixed)
        raw_time_post_t0 <- (-log(u) - hazard_pre_t0*t0) / hazard_post_t0 
        # survival_time_post_t0 is survival time after t0 so t0 is added to get a new overall survival T (Fixed)
        status_post_t0 <- as.integer(raw_time_post_t0 + t0 <= T_max)
        survival_time_post_t0 <- ceiling(min(raw_time_post_t0 + t0, T_max))
        data <- bind_rows(data, c(id=i, starttime=0, endtime=t0, X1=X1[i], X2=X2[i], X3=0, status=0))
        data <- bind_rows(data, c(id=i, starttime=t0, endtime=survival_time_post_t0, X1=X1[i], X2=X2[i], X3=1, status=status_post_t0))
      }
    }
    
    return(data)
  }

生成コードの検証

sim2 <-
  function() {
    
    data <- gen_data2()
    
    fit <- coxph(Surv(starttime, endtime, status) ~ X1+X2+X3, data)
    summary(fit)
    
    est_cox <- fit$coefficients
    ses_cox <- sqrt(diag(fit$var))
    cil_cox <- est_cox + qnorm(.025) * ses_cox
    ciu_cox <- est_cox + qnorm(.975) * ses_cox
    
    c(bias_cox = est_cox - true_beta, 
      coverage_cox = ((true_beta > cil_cox) & (true_beta < ciu_cox))
    )
    
  }

# Get the means of replications
sim_res2 <- rowMeans(replicate(100, sim2()))

実行結果。~~時間依存変数だけ少しバイアスが大きいような気がしないでもない・・・。~~

<after>

> sim_res2
    bias_cox.X1     bias_cox.X2     bias_cox.X3 coverage_cox.X1 coverage_cox.X2 
   -0.009181562     0.001042621    -0.004604479     0.930000000     0.950000000 
coverage_cox.X3 
    0.940000000

ironman

Austin (2012) 再現の続き

論文リンク
3.2 Continuous time-varying covariateの式(4)を再現。~~しようとしてうまくいっていない。~~
これは、時間依存変数の値が時間tに比例するという仮定を置いた場合に導出される式を用いているとのこと。
(4)式自体は、 $-log(u)=H(t)$ という上のスクラップでも書いた関係を利用して、以下のハザードの式変形を愚直にやって $t=hoge$ の形にすれば得られる。

H(t) = \int_{0}^{t} h_{0}\,e^{\beta_{3} k s}\,ds = \frac{h_{0}}{\beta_{3} k}\Bigl(e^{\beta_{3} k t}-1\Bigr)

データ生成関数の作成

[2025-05-05 追記]
o4-mini-high氏との対話を繰り返すこと数回・・・データ生成の時のX3の作り方を計数過程の時間幅の中点にするというアイデアを提案されその修正を実行したらうまくいったぽい・・・。

<修正後のコード>

# Set parameters -------
N <- 500       # Number of subjects
T_max <- 365  # Maximum time in days
lambda <- 0.01  # Base hazard rate
true_beta <- c(0.2, -0.3, 0.5)  # Coefficients for covariates

## data generating function -------
gen_data3 <- function() {
  X1 <- rnorm(N); X2 <- rnorm(N)
  out <- vector("list", N)

  for (i in seq_len(N)) {
    k  <- runif(1, 0.01, 1) 
    u  <- runif(1)
    h0 <- lambda * exp(true_beta[1]*X1[i] + true_beta[2]*X2[i])

    insidelog <- 1 + (true_beta[3]*k * -log(u)) / h0

    # insidelogが1以下だと不定になるので
    # insidelogはuが[0,1]を動くので-log(u)>=0だが、
    # true_beta[3]を負にしたりする場合、1以下になる可能性もある
    if (insidelog < 1) {
      T_cont   <- T_max; status_i <- 0
    } else {
      T_tmp    <- (1/(true_beta[3]*k)) * log(insidelog)
      status_i <- as.integer(T_tmp <= T_max)
      T_cont   <- min(T_tmp, T_max)
    }

    # 日次区間
    n_days <- ceiling(T_cont)
    startv <- 0:(n_days-1)
    endv   <- 1:n_days

    # 中点で評価 これをstartvにするとカバレッジが下がるしバイアスも出る。一番肝になるところ。
    midv   <- (startv + endv) / 2
    X3v    <- k * midv

    status_vec <- rep(0, n_days)
    status_vec[n_days] <- status_i

    df_i <- data.frame(
      id        = i,
      starttime = startv,
      endtime   = endv,
      X1        = X1[i],
      X2        = X2[i],
      X3        = X3v,
      status    = status_vec
    )

    out[[i]] <- df_i
  }

  data_cp <- do.call(rbind, out)
  rownames(data_cp) <- NULL
  return(data_cp)

  }

検証のための関数の作成

sim3 <-
  function() {
    
    data <- gen_data3()
    
    fit <- coxph(Surv(starttime, endtime, status) ~ X1+X2+X3, data)
    summary(fit)
    
    est_cox <- fit$coefficients
    ses_cox <- sqrt(diag(fit$var))
    cil_cox <- est_cox + qnorm(.025) * ses_cox
    ciu_cox <- est_cox + qnorm(.975) * ses_cox
    
    c(bias_cox = est_cox - true_beta, 
      coverage_cox = ((true_beta > cil_cox) & (true_beta < ciu_cox))
    )
    
  }

検証

rep=100だと数字のばらつきが大きかったのでrep=500にした。いい感じ！！
~~肝心の時間依存変数X3のカバレッジが95%前後にならず、修正方法募集中。~~

> # Get the means of replications
> system.time(
+   sim_res3 <-
+     rowMeans(replicate(500, sim3()))
+ )
   user  system elapsed 
 75.778   0.254  76.227 
> 
> sim_res3
    bias_cox.X1     bias_cox.X2     bias_cox.X3 coverage_cox.X1 coverage_cox.X2 
   -0.003397825     0.005024217    -0.006335695     0.964000000     0.960000000 
coverage_cox.X3 
    0.948000000