上一篇文章中，藉由比較存活分析與時間序列之間的差異，了解存活分析關心的課題。這次，終於要進入估計的部分了。
一般的教科書或是存活分析的課程，大多會從有母數方法(Parametric methods)開始教起，也就是從介紹存活分析常見的機率分布開始，推導它們的存活函數或危險函數。因為它們的模型可解釋性強，而且不需要大量樣本，只要能滿足它們的統計假設，實務上的表現會很優異。
不過本次文章略過它，想要從不需統計假設的無母數方法(Nonparametric method)開始，雖然教授的順序比有母數還要後面，但其實有些無母數方法由來已久。本次介紹的生命量表(life table)最早可以追溯到17世紀，使用的技術也不複雜，會數數就可以了。另一個方法 Kaplan-Meier 就比較晚，20 世紀上半才由 Kaplan 提出。
生成資料
我們模擬的資料來自指數分配 \(\text{EXP} (0.1)\) 的資料。指數分配簡單介紹如下：
如果事件時間 \(T \sim \text{Exponential}(\lambda)\)，其中\(\lambda\)為每單位時間發生該事件的次數，那麼：
- 機率密度函數： \(f(t) = \lambda e^{-\lambda t}\)
- 存活函數： \(S(t) = P(T > t) = e^{-\lambda t}\)
- hazard function： \(h(t) = \lambda\)（
這次我們鎖定 存活函數，將\(\lambda\)設為0.1/年，樣本設定為 100 生成資料：
set.seed(1005)
# 模擬參數
n <- 100 # 模擬個體數量
lambda <- 0.1 # 平均存活時間 = 1/lambda = 10
# 模擬真實存活時間 (from exponential)
true_surv_time <- rexp(n, rate = lambda)
# 模擬隨機截尾時間（假設最大觀察時間是 12）
censor_time <- runif(n, min = 5, max = 12)
# 觀察到的時間和事件狀態
obs_time <- pmin(true_surv_time, censor_time)
status <- ifelse(true_surv_time <= censor_time, 1, 0)
# 整理成資料表
data <- data.frame(
id = 1:n,
time = round(obs_time, 2),
status = status,
true_surv = true_surv_time
)
head(data)
id time status true_surv
1 1 9.88 0 16.0519022
2 2 1.57 1 1.5700031
3 3 10.50 0 14.5478572
4 4 0.93 1 0.9324826
5 5 1.20 1 1.1950362
6 6 2.60 1 2.6003750
來看看資料的長相：
Statistic Obs_Time True_Surv
1 Min. 0.2400 0.2364795
2 1st Qu. 2.2725 2.2699674
3 Median 5.1300 5.8076985
4 Mean 5.4698 10.6284341
5 3rd Qu. 8.0175 14.2577409
6 Max. 11.9300 69.9463918
Life table
Life table 的先決條件是要為時間切出一樣等長的單位(\(t_1, t_2, , \dots ,t_j\))，用這些做出時間間隔(區間)。再來是在每個個時間區間裡的位於時間段落裡的人數、死亡人數、存活人數，然後就可以推算出存活函數，具體來說，Life table 關心的統計圖是上次提到，原始資料按實際時間排列的圖：
也因此符號的定義會跟上次提供的有點不同，這裡列出 :
\(t_i\): 為時間點，區間則為 \(I_i = [t_{i-1},t_i)\)
\(z_i\): 第 \(i\) 個區間長度
\(n_i\): 第 \(i\) 個區間內可追蹤的樣本數，等於下面的\(d_i\)跟\(l_i\)的總和
\(d_i\): 第 \(i\) 個區間內死亡的樣本數
\(l_i\): 第 \(i\) 個區間內 censer (活著離開研究)的樣本數
這裡有時會再細分，例如用 \(l_i\)表示中途離開研究，不知後續生死的樣本數，\(w_i\)表示第 \(i\) 個區間內確定活著離開研究的樣本數。
- \(w_i=n_i-l_i*0.5\): 第 \(i\) 個區間具有死亡風險的的平均個體數，概念是：因為不確定將活著離開的人未來是否會在符合研究條件下去世，因此假設 censer 時間服從均勻分配，讓活著離開的人加權值為0.5扣除。化簡後得到\(d_i+l_i*0.5\)
其中，\(w_i\)的計算顯然是估計出來的，這也是會導致後續估計存活函數時會有偏差的原因，繼續看下去就知道了。
然後我們就可以求出：
\(d_{t_j}\): 在時間點 \(t_j\) 時的死亡人數，即\(P(T \in I_i|T>I_{i-1})\)
\(n_{t_j}\): 在 \(t_j\) 時的還活著的人數，即\(P\)
\(\hat q_j = \frac{d_j}{n_j}\): 至特定時間點\(t_j\)的死亡率。
\(\hat p_j= 1- \hat q_j\): 至特定時間點\(t_j\)的存活率 (\(\to\)這是我們的目標之一
為什麼 \(\hat p_j\) 是我們的目標 ? 這是因為存活函數\(S(t)\)可以拆解成:
\(\begin{array}{lll}
S(t) & = & P(T>t_i) \\
& = & P(T>t_1)P(T>t_2|T>t_1)\dots P(T>t_j|T>t_{j-1})
& = \prod _{i=1}^j p_i
\end{array}\)
如此一來，我們會需要每個時間點的 \(p_j\) ，才能湊出 \(S(t)\)，而估計出來的 \(\hat p_j\) 因為是至某個時間點產生的東西，所以隨之估計出來的 \(\hat S(t)\) 會是一個階梯形狀的函數。
# 分區間（每年）
breaks <- seq(0, 13, by = 1)
data$interval <- cut(data$time, breaks = breaks, right = FALSE, labels = FALSE)
# 初始化生命表資料框
life_table <- data.frame(
interval = paste0(breaks[-length(breaks)], "-", breaks[-1]),
n = integer(length(breaks) - 1),
deaths = integer(length(breaks) - 1),
censored = integer(length(breaks) - 1),
w = numeric(length(breaks) - 1),
q = numeric(length(breaks) - 1),
p = numeric(length(breaks) - 1),
lx = numeric(length(breaks) - 1)
)
# 初始存活人數（lx）設為 1，代表比例生命表
lx <- 1
censored <- n
for (i in 1:(length(breaks)-1)) {
in_interval <- data$interval == i
deaths <- sum(data$status[in_interval] == 1, na.rm = TRUE)
if(sum(data$status[in_interval] == 0, na.rm = TRUE) != 0){
censored <- sum(data$status[in_interval] == 0, na.rm = TRUE)
}
w <- deaths + 0.5 * censored
# q <- ifelse(w > 0, deaths / w, 0) # 死亡率(單點)
q <- deaths/w
p <- 1 - q # 存活率(單點)
# if(round(p,1) == 0){
# lx <- 1
# }else {
# lx <- lx*p
# }
lx <- lx * p
life_table$n[i] <- sum(in_interval, na.rm = TRUE)
life_table$deaths[i] <- deaths
life_table$censored[i] <- censored
life_table$w[i] <- round(w, 2)
life_table$q[i] <- round(q, 4)
life_table$p[i] <- round(p, 4)
life_table$lx[i] <- round(lx, 4)
# lx <- lx * p # 更新下一期的 lx
}
interval n deaths censored w q p lx
1 0-1 8 8 10 13.0 0.6154 0.3846 0.3846
2 1-2 14 14 10 19.0 0.7368 0.2632 0.1012
3 2-3 9 9 10 14.0 0.6429 0.3571 0.0361
4 3-4 8 8 10 13.0 0.6154 0.3846 0.0139
5 4-5 8 8 10 13.0 0.6154 0.3846 0.0053
6 5-6 12 5 7 8.5 0.5882 0.4118 0.0022
7 6-7 5 2 3 3.5 0.5714 0.4286 0.0009
8 7-8 9 3 6 6.0 0.5000 0.5000 0.0005
9 8-9 6 2 4 4.0 0.5000 0.5000 0.0002
10 9-10 7 0 7 3.5 0.0000 1.0000 0.0002
11 10-11 7 0 7 3.5 0.0000 1.0000 0.0002
12 11-12 7 0 7 3.5 0.0000 1.0000 0.0002
13 12-13 0 0 7 3.5 0.0000 1.0000 0.0002
欄位解說
interval |
區間（年份） |
觀察時間的分段，例如 0–1 年、1–2 年…等。 |
n |
區間樣本數 |
落在該時間區間的個體數量（包含死亡與審查）。 |
deaths |
死亡人數 |
在該區間內發生事件（如死亡）的個體數。 |
censored |
審查人數 |
在該區間內未發生事件，但中途退出觀察的個體數（例如失聯、試驗結束）。 |
w |
平均風險人數 |
用於估計死亡率的有效樣本數，通常計算為 deaths + 0.5 × censored。 |
q |
區間死亡率 |
在該區間內死亡的機率，公式為 q = deaths / w。 |
p |
區間存活率 |
存活機率，p = 1 - q，表示該區間內個體存活下來的機率。 |
lx |
存活比例 |
到該時間區間結束時，仍然存活的比例（累積存活率），初始為 1（或 100%）。 |
圖為：
time <- 0:(length(life_table$lx)-1)
plot(time, life_table$lx, type = "s", col = "blue", lwd = 2, xlab="觀察時間",ylab="存活函數")
# 加入點
points(time, life_table$lx, pch = 16, col = "blue")
# # 加入真實模擬存活曲線（指數分布）
t_seq <- seq(0, max(time), by = 0.1)
true_surv <- exp(-0.1 * t_seq)
lines(t_seq, true_surv, col = "red", lwd = 2, lty = 2)
# 加上圖例
legend("topright",
legend = c("生命表 (手刻)", "真實指數分布"),
col = c("blue", "red"),
lwd = 2,
lty = c(1, 2),
bty = "n"
)
可以看到距離真實值偏離了不少，實務上會試著調整\(w_i\)的計算方式，不過下一個介紹的 Kaplan-Meier 可以很好的改進這點。
Kaplan-Meier
Kaplan-Meier 跟 Life table 最大的差異在於，Life table 的時間區段是事先定義，Kaplan-Meier 的時間區段會發生在\([\text{發生事件的時間點},\text{發生事件的下一個時間點})\)，這裡的事件指的即是病患死亡或離開研究。也就是說，每次的時間間隔不一定會一樣，產生的 \(\hat S(t)\) 會是大小不一的階梯狀函數。
Kaplan-Meier 所關心的資料畫成圖就像這樣，已經專注在存活長度本身上
實務上， Life table 常用於政府或研究機構定期發布的調查資料，而 Kaplan-Meier 則用在可以具體知道事件發生時間點的臨床研究資料上。
為了與前面的 Life table 法區隔，我們改用\((v_1, \dots ,v_j)\)標示Kaplan-Meier的時間間隔。
\(T_i\) : 按資料存活時間/截尾時間排序後，第\(i\)位對象的存活時間/截尾時間，是個隨機變數。
\(c\) :研究開始到研究結束的時間
\(\delta_i\) : 研究結束時是否過世/死亡
\(U_i = \left\{
\begin{array}{r}
T_i \quad \text{if }T_i \leq C_r \\
C_r \quad \text{if }T_i > C_r
\end{array} \right .\)
\(d_{v_j}\): 在時間點 \(v_j\) 時的死亡人數，即\(P(T \in [v_{j-1},v_j)|T>v_{j-1})\)
\(n_{v_j}\): 在 \(v_j\) 時的還活著的人數，即\(P\)
\(\hat q_j = \frac{d_{v_j}}{n_{v_j}}\): 至特定時間點\(v_j\)的死亡率。
\(\hat p_j= 1- \hat q_j\): 至特定時間點\(v_j\)的存活率 (\(\to\)這是我們的目標之一
而我們的目標 \(\hat S(t) = \prod _{i=1}^j \hat p_i\) ， 原因前面已提過，就不再贅述。
# sort data
data <- data[order(data$time), ]
# Step 1: 只取有事件的時間
event_times <- data$time[data$status == 1]
unique_event_times <- sort(unique(event_times))
# Step 2: 建立表格
km_table <- data.frame(
time = unique_event_times,
n_risk = NA,
n_event = NA,
survival_prob = NA,
cumulative_survival = NA
)
n <- nrow(data)
cum_surv <- 1
for (i in 1:nrow(km_table)) {
t <- km_table$time[i]
at_risk <- sum(data$time >= t)
events <- sum(data$time == t & data$status == 1)
p <- 1 - events / at_risk # 計算存活機率
cum_surv <- cum_surv * p # 計算累積存活機率
km_table$n_risk[i] <- at_risk
km_table$n_event[i] <- events
km_table$survival_prob[i] <- round(p, 3)
km_table$cumulative_survival[i] <- round(cum_surv, 3)
}
# 加入真實值
# 求出"某個時間點 t 的「還活著的機率」"
lambda_t <- exp(-lambda * km_table$time)
km_table$true_survival <- round(lambda_t, 3)
time n_risk n_event survival_prob cumulative_survival true_survival
1 0.24 100 1 0.990 0.990 0.976
2 0.31 99 1 0.990 0.980 0.969
3 0.35 98 1 0.990 0.970 0.966
4 0.37 97 1 0.990 0.960 0.964
5 0.52 96 1 0.990 0.950 0.949
6 0.81 95 1 0.989 0.940 0.922
7 0.93 94 2 0.979 0.920 0.911
8 1.06 92 1 0.989 0.910 0.899
9 1.08 91 1 0.989 0.900 0.898
10 1.15 90 1 0.989 0.890 0.891
11 1.20 89 1 0.989 0.880 0.887
12 1.25 88 1 0.989 0.870 0.882
13 1.27 87 2 0.977 0.850 0.881
14 1.57 85 2 0.976 0.830 0.855
15 1.62 83 1 0.988 0.820 0.850
16 1.79 82 1 0.988 0.810 0.836
17 1.89 81 1 0.988 0.800 0.828
18 1.90 80 1 0.988 0.790 0.827
19 1.93 79 1 0.987 0.780 0.824
20 2.03 78 1 0.987 0.770 0.816
21 2.06 77 1 0.987 0.760 0.814
22 2.19 76 1 0.987 0.750 0.803
23 2.30 75 1 0.987 0.740 0.795
24 2.37 74 1 0.986 0.730 0.789
25 2.47 73 1 0.986 0.720 0.781
26 2.60 72 1 0.986 0.710 0.771
27 2.71 71 1 0.986 0.700 0.763
28 2.94 70 1 0.986 0.690 0.745
29 3.10 69 1 0.986 0.680 0.733
30 3.13 68 1 0.985 0.670 0.731
31 3.29 67 1 0.985 0.660 0.720
32 3.41 66 1 0.985 0.650 0.711
33 3.42 65 2 0.969 0.630 0.710
34 3.52 63 1 0.984 0.620 0.703
35 3.88 62 1 0.984 0.610 0.678
36 4.25 61 1 0.984 0.600 0.654
37 4.28 60 1 0.983 0.590 0.652
38 4.49 59 1 0.983 0.580 0.638
39 4.72 58 1 0.983 0.570 0.624
40 4.75 57 1 0.982 0.560 0.622
41 4.91 56 1 0.982 0.550 0.612
42 4.92 55 1 0.982 0.540 0.611
43 4.99 54 1 0.981 0.530 0.607
44 5.01 53 1 0.981 0.520 0.606
45 5.19 49 1 0.980 0.509 0.595
46 5.80 44 1 0.977 0.498 0.560
47 5.81 43 1 0.977 0.486 0.559
48 5.85 42 1 0.976 0.475 0.557
49 6.38 40 1 0.975 0.463 0.528
50 6.88 37 1 0.973 0.450 0.503
51 7.11 36 1 0.972 0.438 0.491
52 7.30 32 1 0.969 0.424 0.482
53 7.68 29 1 0.966 0.409 0.464
54 8.31 23 1 0.957 0.392 0.436
55 8.99 22 1 0.955 0.374 0.407
欄位解釋
time |
發生事件的時間點（如死亡、復發） |
n_risk |
在該時間點仍處於風險中的人數（尚未發生事件或被截尾） |
n_event |
在該時間點發生事件的人數 |
survival_prob |
單一時間點的存活機率（即 (1 - )） |
cumulative_survival |
累積存活率，表示從起始時間到該時間點的整體存活機率 |
true_survival |
真實存活率（按照模擬資料指定分布 \(\text{EXP}(0.1)\) 反推） |
表格解讀
Kaplan-Meier 曲線是階梯狀的，每次事件發生時下降，未發生事件的時間點則維持不變。這份表格展示了：
- 初始時刻（time = 0.24）：有 100 人在風險中，發生 1 件事件，存活機率為 0.990，累積存活率也為 0.990。
- 隨時間推移：每次事件發生都會使累積存活率下降。例如：
- 在 time = 1.27，有 2 件事件，累積存活率從 0.870 降至 0.850。
- 到 time = 5.01，累積存活率已降至 0.520，表示只有 52% 的人存活至該時間點。
true_survival 欄位：用來與 Kaplan-Meier 所估計的 cumulative_survival 做比較。
畫圖觀察
# 畫 Kaplan-Meier 估計的存活機率（step line）
plot(km_table$time, km_table$cumulative_survival, type = "s",
col = "blue", lwd = 2, ylim = c(0, 1),
xlab = "Time", ylab = "Survival Probability",
main = "Kaplan-Meier vs True Survival Function")
# 加上真實存活函數的線（用對應時間點計算的真實值）
lines(km_table$time, km_table$true_survival, col = "red", lwd = 2, lty = 2)
# 加圖例
legend("topright", legend = c("KM Estimate", "True S(t)"),
col = c("blue", "red"), lwd = 2, lty = c(1, 2))
畫圖觀察後可以發現，Kaplan-Meier 在樣本數為100下，比起 life table ， Kaplan-Meier fiting 得更好。
