初探存活分析(3)

如何估計Cox model 的係數

Survial Analysis
R
作者

紙魚

發佈於

2026年1月16日

摘要
探討 Cox model 的係數如何導出，以及觀察它怎麼被刻出來。

緣起

本來這學期的課程加分作業有一個手刻 Cox model 係數，但是我卡在刻不出來QQ。1當然最後也沒交，但是又覺得這個題目有趣加上都寫一半了，不要浪費任何可以交 kpi 的機會，所以只好事後叫 AI 刻刻來觀察，順便看看能不能偷學幾招刻函數的技巧。

還好寫了其他的加分作業有過，存活大成功

Cox model 概念

先來看看 Cox model 的概念，它假設每個個體\(i\)的 hazard function如下：

\[h(t|Z_i) = h_0(t) e^{\beta^\top Z_i}\]

  • \(h(t|Z_i)\): 就是在時間點 \(t\) 且事件尚未發生的時候，發生事件的機率，例如：如果這裡的事件發生是指死亡，那就是特定時間點\(t\)的死亡率。

  • \(h_0(t)\): 單一個體的 base-line hazar function，白話來說就是個體身上其他影響發生事件與否的風險因子。Cox model 的理念是剃除時間及其他變因，以便觀察主要風險因子的影響，所以這個\(h_0(t)\)在估計的過程中會被消除。

  • \(Z_i\): 是個向量，包含\(p\)個影響發生事件與否的主要風險因子，也是研究者收集回來的資料們。

  • \(\beta\): 估計係數\((\beta_0,\beta_1,\dots,p)^T\)，是\(p\times1\)的矩陣，影響主要風險因子的權重(weight)，是研究者最關心的東西，也是本文的目標。

對第\(i\)個事件發生的個體\(r_i\)（即第\(i\)個死亡的個體），給定已經發生的過去事件，這名個體上事件發生的機率是：

\[p(r_i | \mathcal{H}_i) = \frac{e^{\beta^\top Z_{(i)}}}{\sum_{j \in R(t_{(i)})} e^{\beta^\top Z_j}}\]

  • \(R(t_{(i)})\) 是在 \(t_{(i)}\) 時仍在風險集的個體集合（未發生事件或 censor）。白話來說就是還在時間點\(t_{(i)}\)中，資料中的存活個數。

這個東西才是我們真正要拿來估計\(\beta\)時會用到的式子。

Cox model 係數數學上是如何得出的

Cox model 是用 likelihood method 得出係數，先來導一下：

\[L(\beta )= \prod \limits_ {i=1}^n p(r_i|\mathcal{H} _i)= \prod_{i=1}^n \left( \frac{e^{\beta^\top Z_{(i)}}}{\sum_{j \in R(t_{(i)})} e^{\beta^\top Z_j}} \right)^{\delta_{(i)}}\]

取 log

\[ \log L(\beta) = l(\beta)= \sum_{i=1}^{n} \left[ \beta^{T} z_{(i)} - \log \left\{ \sum_{j \in R(t_{(i)})} e^{\beta^{T} z_j} \right\} \right] \delta_{(i)} \]

到這裡都是非常正規的 likelihood method 流程，接下來就是比較不一樣的地方了！

\[ U(\beta) = \frac{\partial l(\beta)}{\partial \beta} = \sum_{i=1}^{n} \left[ z_{(i)} - \frac{ \sum_{j \in R(t_{(i)})} e^{\beta^{T} z_j } z_j }{ \sum_{j \in R(t_{(i)})} e^{\beta^{T} z_j} } \right] \delta_{(i)} \]

這個\(U(\beta)\)是甚麼東東呢？數學上是對 log 化後的 likelihood 偏微，是一種 score function，但它也可寫成：

\[U(\beta) = \sum_{i=1}^n \delta_{(i)} \left[ Z_{(i)} - \bar Z_{(i)} \right]\]

看起來很像什麼東西？當第 i 個的事件發生時，\(\delta_{(i)}\)會標記為 1 ，這個式子便是在衡量，當實際發生事件者的共變數(covariates)與平均的差距。

共變數也叫做共變量，直觀就是餵進去模型的資料，不過常常會看到定義相同，但是術語不同的狀況，例如：在計量經濟學裡會稱呼為獨立變數(independent variable)；機器學習會叫做 inputs 。這是因為統計是各種領域交會的領域，所以有些用語不一定會統一，如果要再算上中文翻譯，可能就會更加多樣了…。

順便一提，這個用語最早是由統計學的奠基者—費雪(R.A Fisher)於 1940 年代發明。

要求得\(\beta\)，就要求 \(U(\beta)\) 的 MLE，也就是利用一皆微分

\[U(\beta) \overset{\text{set}}{=} 0\]

並且估出 Fisher information matrix，為了簡化式子，定義以下：

令 \(Z_j\) 為 p 維的共變數。

  • \(\Gamma_0(\beta) = \sum_{j \in R(t_i)} e^{\beta^\top Z_j},\)

  • \(\Gamma_1(\beta) = \sum_{j \in R(t_i)} e^{\beta^\top Z_j} Z_j,\)

  • \(\Gamma_2(\beta) = \sum_{j \in R(t_i)} e^{\beta^\top Z_j} Z_j^{\otimes 2},\)

    • 其中\(Z_j^{\otimes 2} = Z_j Z_j^\top\)

代入 \(U(\beta)\)：

\[U(\beta) = \sum_{i=1}^n \left[z_{(i)}-\frac{\Gamma_1(\beta)}{\Gamma_0(\beta)}\right] \delta_{(i)}\]

因此 Fisher information matrix 為：

\[ \mathcal{I}(\beta) = - \frac{\partial^2 \ell(\beta)}{\partial \beta^2} = \sum_{i=1}^n \left[ \frac{\Gamma_2(\beta)}{\Gamma_0(\beta)} - \frac{\Gamma_1(\beta)\Gamma_1(\beta)^\top}{\Gamma_0(\beta)^2} \right] \delta(i)\]

不過很遺憾，努力導出來的這些東西是沒有辦法直接算出來的QQ，因此，我們需要牛頓法來迭代，求出近似值。

牛頓法的運用

回顧微積分牛頓法的概念，設定初始值\(x_0\)，在可微函數\(f(x)\)上求以下：

\[0=(x_{n+1}-x_n)f'(x_n)+f(x_n) ~~~~~~~ n=0,1,\dots\]

也就是

\[x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}\]

以上是微的變數只有一個的版本，要微很多個變數是用泰勒展開式表示。

牛頓法也相當於泰特展開式的只取到一階微分的版本，運用在這裡的話就是把\(x_0\)換成\(\beta\)，理想的正確答案為\(\hat{\beta}\)：

\[\begin{aligned}U(\hat{\beta})-U(\beta)& \approx \frac{\partial U(\beta)}{\partial \beta}(\hat{\beta}-\beta)\\ & = -\mathcal{I}(\beta)(\hat{\beta}-\beta) \end{aligned}\]

因為\(\hat{\beta}\)是理想的正確答案(a.k.a MLE)，所以\(U(\hat{\beta})=0\)，然後再經過一點小小的線性代數運算，求得\(\hat{\beta}\)：

\[\hat{\beta}\approx \beta +U(\beta)\mathcal{I}^{-1}(\beta)\]

這就是牛頓法的運用。

實作

模擬係數

先模擬一份假資料：

  • 存活分配來自 \(\exp(\lambda=0.1)\)

  • 含有來自 Bernoulli \(Ber(0.5)\)的 \(X_1\)(類別變數) 及來自常態分配\(N(0,1)\)的 \(X_2\)(連續變數)

  • 樣本為 1000，censor的機率來自\(U(10,20)\)

  • 分別給予 \(X_1\) 跟 \(X_2\)初始的\(\beta_1=0.7\)與\(\beta_1=0.5\)

set.seed(20260116)

# 模擬參數
n <- 1000
lambda <- 0.1       # 基礎事件率
beta1 <- 0.7        # X1 的 log-hazard 影響
beta2 <- 0.5        # X2 的 log-hazard 影響

# 模擬自變量 X1、X2
X1 <- rbinom(n, 1, 0.5)   # 50%機率
X2 <- rnorm(n, mean = 0, sd = 1)   # 標準常態分布

# 存活時間，考慮 X1 與 X2 的影響
true_surv_time <- rexp(n, rate = lambda * exp(beta1*X1 + beta2*X2))

# censor 時間
censor_time <- runif(n, min = 10, max = 20)

# 觀察時間與事件狀態
obs_time <- pmin(true_surv_time, censor_time)
status <- ifelse(true_surv_time <= censor_time, 1, 0)

# 整理成 data frame
data <- data.frame(
  id = 1:n,
  X1 = round(X1, 2),
  X2 = round(X2, 2),
  time = round(obs_time, 2),
  status = status,
  true_surv = round(true_surv_time, 2)
)

# 查看資料結構
str(data)
'data.frame':   1000 obs. of  6 variables:
 $ id       : int  1 2 3 4 5 6 7 8 9 10 ...
 $ X1       : num  1 0 1 1 1 1 1 1 0 1 ...
 $ X2       : num  -0.09 -0.5 1.03 0.25 -0.31 -0.53 -0.73 1.18 0.91 0.67 ...
 $ time     : num  0.01 6.49 1.05 1.8 0.83 0.39 8.74 2.33 0.42 3.62 ...
 $ status   : num  1 1 1 1 1 1 1 1 1 1 ...
 $ true_surv: num  0.01 6.49 1.05 1.8 0.83 0.39 8.74 2.33 0.42 3.62 ...

基本資料

[1] "資料收集部分："
[1] "censor 數: 183"
[1] "event 數 : 817"
[1] "----------------"
[1] "  total  : 1000"
[1] "變數部分："
[1] "X1=0 數: 533"
[1] "X1=1 數: 467"

X2是連續變數，且已知放入的平均值跟標準差，故省略。

我們可以用 R 裡面的 package survival估計正確答案。

library(survival)
cox.ph <- coxph(Surv(time,status) ~ as.factor(X1)+X2, data=data)
s2 <- summary(cox.ph)

s2$coefficients
                    coef exp(coef)   se(coef)        z     Pr(>|z|)
as.factor(X1)1 0.9017075  2.463806 0.07247418 12.44178 1.550194e-35
X2             0.4892881  1.631155 0.03712873 13.17815 1.172262e-39

裡面的0.9017075跟0.4892881分別為X1與X2正確係數答案。

這是因為存活資料多了 censor 的狀態，資料不齊全的情況注定無法給出完全正確的答案，在樣本不夠多但 censor 的比例高的情況下就會影響真正估計的係數；反之樣本夠多或censor 的比例夠低，估計出來的係數就會非常接近正確答案。

為了加快網頁的渲染，本例設1000，有興趣的人可以把程式碼複製起來把樣本數調大，就會很接近設定的

AI 給的手刻程式

AI 總共給了 2 個手刻程式解，第一個是比較長，但是比較直觀的版本：

牛頓法版本

先將資料按時間先後排序以及轉矩陣跟設變數等，以便進行後面的刻算法。

data <- data[order(data$time), ]
X <- as.matrix(data[, c("X1", "X2")]) # 大矩陣，裡面的其中一行就是 Z_j
time <- data$time
status <- data$status
n <- nrow(X)
p <- ncol(X)

\(U(\beta)\)的函數版

score_cox <- function(beta, X, time, status) {
  eta <- as.vector(X %*% beta) # 注意 beta在這裡是個向量，所以(beta^TZ)^T=Xbeta
  U <- rep(0, length(beta))　# 初始係數矩陣(U(beta))
  
# sum from i=1 to n  且 delta =1 式子才能成立
  for (i in which(status == 1)) {
    risk <- which(time >= time[i]) # risk set 數量
    w <- exp(eta[risk]) 
    # 算Z bar
    X_bar <- colSums(X[risk, , drop = FALSE] * w) / sum(w)
    # U 是 sum 起來的所以要加上
    U <- U + X[i, ] - X_bar
  }
  
  U
}

就是直白的把\(U(\beta)=\sum_{i=1}^n \delta_{(i)} \left[ Z_{(i)} - \bar Z_{(i)} \right]\)做出來。

再來是\(I(\beta)\)的函數版

info_cox <- function(beta, X, time, status) {
  eta <- as.vector(X %*% beta)
  I <- matrix(0, ncol(X), ncol(X))
  
# sum form i = 1 to n 且 delta =1 式子才能成立
  for (i in which(status == 1)) {
    risk <- which(time >= time[i])
    w <- exp(eta[risk])
    
    X_bar <- colSums(X[risk, , drop = FALSE] * w) / sum(w)
    # 相當於算 bar Z _{(i)}
    X_centered <- sweep(X[risk, , drop = FALSE], 2, X_bar)
    I <- I + t(X_centered) %*% (X_centered * w) / sum(w)
  }
  
  I
}
語法字典

  • colSums()：對矩陣中每一欄進行加總

  • sweep()：用一個向量掃過矩陣的某個方向，做逐元素運算，塞入的input 為：sweep(x, MARGIN, STATS, FUN = "-")

參數 意思
x 矩陣
MARGIN 1 = row，2 = column
STATS 要套用的向量
FUN 運算（預設是減法）

再來刻牛頓法：

newton_cox <- function(X, time, status,
                       beta_init = c(0, 0),
                       tol = 1e-6,
                       maxit = 50) {
  
  beta <- beta_init
  
  for (iter in 1:maxit) {
    U <- score_cox(beta, X, time, status)
    I <- info_cox(beta, X, time, status)
    
    step <- solve(I, U) #求UI^{-1}
    beta_new <- beta + step
    
    # 如果係數幾乎沒有改變就停止迭代
    if (max(abs(beta_new - beta)) < tol) {
      cat("Converged at iteration", iter, "\n")
      return(list(beta = beta_new, I = I))
    }
    
    beta <- beta_new
  }
  
  stop("Did not converge")
}

刻完了！來算算看

fit_newton <- newton_cox(X, time, status)
Converged at iteration 4 
fit_newton$beta
       X1        X2 
0.9009551 0.4887990

對照正確答案

                    coef exp(coef)   se(coef)        z     Pr(>|z|)
as.factor(X1)1 0.9017075  2.463806 0.07247418 12.44178 1.550194e-35
X2             0.4892881  1.631155 0.03712873 13.17815 1.172262e-39

非常接近，不錯不錯。

AI 給的第二種答案

這個答案的寫法比較精簡

cox_ploglik <- function(beta, time, status, X) {
  eta <- as.vector(X %*% beta)
  loglik <- 0
  
  for (i in which(status == 1)) {
    risk_set <- which(time >= time[i])
    loglik <- loglik +
      eta[i] -
      log(sum(exp(eta[risk_set])))
  }
  
  loglik
}

就這樣，因為它使用時用了一個現成的optim()來求解，就不用手刻牛頓法。

X <- as.matrix(data[, c("X1", "X2")])

neg_pl <- function(beta) {
  -cox_ploglik(beta, data$time, data$status, X)
}

fit <- optim(par = c(0, 0),
  fn = neg_pl,
  method = "BFGS",
  hessian = TRUE)

fit$par
[1] 0.9009556 0.4887988

也是相當接近。

結論

刻函數還是好難TT，不過 AI 裡面用到的語法跟順序結構滿值得一學的，只是最難的應該還是初始矩陣跟變數要怎麼設、跟數學式的差異，還有迴圈的應用時機吧。

參考資料

余日彰. (2025). 存活分析. 未公開課堂講義.
林建甫. (2008). 存活分析 (初版). 雙葉書局.
無符合的項目

腳註

  1. 老師甚至在課堂上表演了一次，結果回家還是刻不出來，我要哭了…。↩︎