RでCOX比例ハザードモデル

seed <- 20250625
set.seed(seed)

n <- 200 # 患者数

# 1. 共変量の生成
# treatment: 0 = 標準治療, 1 = 新薬
treatment <- rbinom(n, 1, 0.5)
# age: 40-70歳の一様分布から生成し、平均を0に標準化
age_raw <- runif(n, 40, 70)
age <- scale(age_raw, center = TRUE, scale = FALSE)

# 2. 真のパラメータ(β)を設定
# 新薬はリスクを約40%に下げる (exp(-0.9) ≈ 0.40)
# 年齢が1歳上がるとリスクが約3%上がる (exp(0.03) ≈ 1.03)
beta_true <- c(treatment = -0.9, age = 0.03)

# 3. 生存時間の生成 (Weibull分布を利用した逆関数法)
# ベースラインハザード h0(t) をWeibull分布のハザードと仮定してデータを生成する
# h(t) = shape * (lambda^shape) * t^(shape-1)
# このモデルは推定時には使わない。データ生成のためだけ。
weibull_shape <- 1.5
weibull_lambda <- 0.01

# 各個人のハザード比
hazard_ratio <- exp(treatment * beta_true["treatment"] + age * beta_true["age"])
# 逆関数法を用いて生存時間を生成
u <- runif(n)
T_event <- (-log(u) / (weibull_lambda * hazard_ratio))^(1 / weibull_shape)

# 4. 打ち切り時間の生成
# 追跡期間は最大60ヶ月とし、一部はそれ以前にランダムに打ち切られる
T_censor <- runif(n, 30, 60)

# 5. 観測データを作成
time <- pmin(T_event, T_censor) # 観測時間
status <- as.numeric(T_event <= T_censor) # 1: イベント発生, 0: 打ち切り

# データフレームにまとめる
sim_data <- data.frame(time, status, treatment, age)

# 生成されたデータの確認
print(head(sim_data))
cat("\nイベント発生率:", mean(status))

# COXモデル実装
fun_coxph <- function(time, status, X) {
  # 行列形式に変換
  X <- as.matrix(X)
  p <- ncol(X)

  # データを時間でソート（リスクセットの計算を効率化するため）
  ord <- order(time)
  time <- time[ord]
  status <- status[ord]
  X <- X[ord, , drop = FALSE]

  # ニュートン・ラフソン法でβを推定
  beta <- rep(0, p) # 初期値
  max_iter <- 20

  for (iter in 1:max_iter) {
    # 現在のbetaでのリスクスコア
    risk_score <- exp(X %*% beta)

    # リスクセットの合計を効率的に計算 (時間でソート済みなので逆順に累積和を取る)
    S0 <- rev(cumsum(rev(risk_score)))

    # スコアベクトル U(beta) と情報行列 I(beta) を計算
    U <- numeric(p)
    I <- matrix(0, p, p)

    # イベントが発生した時点のみでループ
    event_indices <- which(status == 1)
    for (i in event_indices) {
      # リスクセット内の個体を取得
      risk_set_indices <- i:length(time)

      # S1, S2の計算
      S1_i <- colSums(X[risk_set_indices, , drop = FALSE] * risk_score[risk_set_indices])
      E_i <- S1_i / S0[i]

      # スコアベクトルの更新
      U <- U + (X[i, ] - E_i)

      # 情報行列の更新
      S2_i <- t(X[risk_set_indices, , drop = FALSE] * risk_score[risk_set_indices]) %*% X[risk_set_indices, , drop = FALSE]
      I <- I + (S2_i / S0[i] - E_i %*% t(E_i))
    }

    # betaを更新
    delta_beta <- solve(I, U)
    beta <- beta + delta_beta

    # 収束判定
    if (all(abs(delta_beta) < 1e-6)) {
      cat("収束しました (iteration =", iter, ")\n")
      break
    }
  }

  # 結果の計算
  coef <- as.vector(beta)
  names(coef) <- colnames(X)
  se <- sqrt(diag(solve(I)))
  z_val <- coef / se
  p_val <- 2 * pnorm(-abs(z_val))
  hr <- exp(coef)

  return(list(
    coef = coef,
    hazard_ratio = hr,
    se = se,
    z = z_val,
    p_value = p_val,
    var = solve(I)
  ))
}

# 実装を実行
X_vars <- sim_data[, c("treatment", "age")]
fitresult <- fun_coxph(sim_data$time, sim_data$status, X_vars)

# 実装の結果を表示
cat("\n--- 結果 ---\n")
print(data.frame(
  coef = fitresult$coef,
  exp_coef = fitresult$hazard_ratio,
  se_coef = fitresult$se,
  z = fitresult$z,
  p_value = fitresult$p_value
))