Kwangmin Kim - Censoring의 중요성과 IPW MSM

1 정의

정의: 두 종류의 Censoring

Administrative censoring (행정): 추적 종료 시점의 도달로 사건 미관측. 모든 사람의 행정 종료 시점 동일 (또는 staggered entry 시 가입 일자에 따라 다름). - 외생적 — 실제 사건 시간과 무관 가정 합리적.

Non-administrative censoring (비행정): 추적 손실 (loss to follow-up), 다른 연구로 이동, 경쟁 사건 (사망 외 다른 사건). 무작위 아닐 수 있음. - 내생적 — 처치·결과 위험 인자와 연관 가능.

직관 — 두 censoring 의 결정적 차이: 행정은 언제 추적 끝나는지 정해진 시점. 비행정은 왜 사라졌는지 모르는 시점. 후자는 사라진 사람과 남아 있는 사람의 위험이 다를 수 있어 selection bias 위험. 행정은 이 위험 없음 (시간 의존 변수 외).

정의: 잠재 결과 표기 with Censoring

\(D_k^{a, \bar{c} = \bar{0}}\): 처치 \(A=a\) + 모든 시점 censoring 안 됨 (\(\bar{c} = (c_1, \ldots, c_{k_\text{end}}) = \bar{0}\)) 의 가상 시나리오에서 시점 \(k\) 의 사망 여부.

\(\Pr(D_k^{a, \bar{c}=\bar{0}} = 0 | A = a)\): “모든 사람이 censoring 없이 추적되고 처치 \(a\) 받았을 때 시점 \(k\) 의 생존률”.

이것이 인과 분석의 진짜 추정 대상. 단순 관측 생존률과 다름.

직관 — Superscript \(\bar{c}=\bar{0}\) 의 의미: “censoring 도 인과 처치처럼 다룬다” 는 발상. 우리가 묻는 인과 질문은 “처치 받았고 + censoring 안 됐다면” 의 가상 시나리오. 두 종류의 가상 개입이 결합 된 양.

2 17.3 Censoring 이 만드는 편향

2.1 단순 잔존 비율의 함정

가상의 작은 연구로 시작하자. \(k_\text{end} = 2\), 다음 시나리오:

Hernan 의 작은 사례 (Ch.17.3)

\(\Pr(C_1 = 0) = 1\): 시점 1 까지 아무도 censoring 안 됨.
\(\Pr(D_1 = 0 | C_0 = 0) = 0.9\): 90% 시점 1 까지 생존.
\(\Pr(C_2 = 0 | D_1 = 0, C_1 = 0) = 0.5\): 생존자의 50% 가 시점 2 에서 censoring (무작위).
\(\Pr(D_2 = 0 | C_2 = 0, D_1 = 0, C_1 = 0) = 0.9\): 남은 사람의 90% 시점 2 까지 생존.

잔존 비율: 시점 2 에서 살아 있고 censoring 안 된 사람 비율 = \(1 \times 0.9 \times 0.5 \times 0.9 = 0.405 = 40.5\%\).

진짜 (가상) 생존률: 만약 censoring 없었다면 = \(1 \times 0.9 \times 1 \times 0.9 = 0.81 = 81\%\).

직관 — 두 숫자의 차이가 큰 이유: 잔존 비율 (40.5%) 은 censoring 안 됐고 살아 있는 사람의 비율. 실제 50% 가 censoring 으로 사라져 분모가 절반으로 줄어듦. 이 절반의 사람들이 잘못된 분자에서 빠지면서 분모도 같이 줄어 생존률 추정이 망가짐. 진짜 생존률 81% 와 잔존 40.5% 의 두 배 차이는 censoring 이 만든 인공물.

직관 — 이 사례가 단순한데도 보여주는 것: censoring 이 무작위 (50% 가 무작위로 사라짐) 라도 잔존 비율이 진짜 생존률을 향하지 않음. 무작위 censoring 도 생존률 직접 추정 불가 — 보정 절차 필요.

2.2 Hazard 곱셈 공식의 정확성

같은 사례를 hazard 곱셈으로 다시 보자:

\[\Pr(D_k = 0 | A = a) = \prod_{m=1}^{k} \Pr(D_m = 0 | D_{m-1} = 0, C_m = 0, A = a)\]

각 시점의 조건부 hazard (살아 있고 censoring 안 된 사람 중 사망 비율) 의 곱:

시점 1 hazard: \(1 - 0.9 = 0.1\) → 생존 0.9
시점 2 hazard (살아 있고 censoring 안 된 사람만): \(1 - 0.9 = 0.1\) → 생존 0.9
누적 생존: \(0.9 \times 0.9 = 0.81\) ✓

직관 — 곱셈 공식이 작동하는 이유: 각 시점의 hazard 는 “그 시점에 살아 있고 censoring 안 된 사람” 만 사용. Censoring 으로 사라진 사람은 그 시점부터 분모에서 제외 — 이전 시점의 생존 추정에는 영향 안 줌. 시점별 조건부 추정 → 곱셈 이 censoring 의 인공물 회피.

직관 — Kaplan-Meier 가 이 공식의 비모수 구현: KM 추정량의 매 시점 “살아 있는 사람 / 사망 - 자” 비율이 곧 conditional hazard. 곱셈이 KM 의 본질.

2.3 무작위 vs 비무작위 Censoring

Censoring 가정의 검증 어려움

무작위 censoring 가정 (\(\bar{C} \perp\!\!\!\perp T\)): 사라진 사람과 남은 사람의 사망 위험이 같음.

이 가정 위반 시: - Loss to follow-up 이 처치-실패 환자에 집중 → 잔존 분석에서 처치 효과 과대 추정. - 부작용으로 인한 dropout → 부작용 위험 사람 사라져 처치 안전성 과대 추정.

→ 측정된 covariate \(L\) 로 censoring 모형 적합 + IPW 보정 필요.

직관 — Loss to follow-up 의 위험: 임상시험에서 처치가 효과 없는 환자가 시험 포기하고 다른 처치로 이동 → censoring. 남은 사람 = 처치 효과 있는 사람 → 처치 효과 과대 추정. Censoring 의 비무작위성이 처치 평가의 왜곡 원천.

직관 — Staggered Entry 의 추가 복잡성: 가입 시기가 다르면 행정 censoring 시점도 다름. 후기 가입자는 짧은 추적 → 행정 censoring 빠름. 시기가 결과와 연관 (예: 진단 기술 발전) 이면 staggered entry 자체가 confounding 원천. Baseline 시기 변수 보정 권장.

3 17.4 IPW MSM 으로 인과 추정

3.1 Confounding + Censoring 의 결합 보정

NHEFS 사례에서 두 종류의 보정 동시 필요:

Confounding 보정: 측정된 9 변수 (\(L\)) 로 처치 (\(A\)) 의 외생화.
Censoring 보정: 여기서는 행정 censoring 만 (모든 사람 120 month) — 무작위 가정 합리적.

만약 비행정 censoring (loss to follow-up) 있으면 추가 IPW 가중치 필요.

3.2 잠재 결과 표기

\(D_k^{a, \bar{c}=\bar{0}}\) → 표기 단순화 위해 \(D_k^a\) 로 표기 (이 챕터에서는 \(\bar{c}=\bar{0}\) 가 자명).

추정 대상:

\[\Pr(D_{k+1}^{a=1} = 0) \quad \text{vs.} \quad \Pr(D_{k+1}^{a=0} = 0) \quad \text{for } k = 0, 1, \ldots, k_\text{end}-1\]

직관 — 인과 질문의 명시: “모든 사람이 흡연 중단 했다면 시점 \(k\) 의 생존률” vs “모든 사람이 흡연 계속 했다면 시점 \(k\) 의 생존률”. 단순 비교는 confounding 으로 진짜 차이를 측정 못 함.

3.3 IPW MSM 의 두 단계

Step 1: 안정화 IP 가중치 \(SW^A\)

Ch.12 와 동일 절차:

\[SW_i^A = \frac{f(A_i)}{f(A_i | L_i)}\]

NHEFS 9 변수 + quadratic logistic 으로 \(f(A | L)\) 추정. NHEFS 사례에서 \(SW^A\) 평균 1.00, 범위 0.33~4.21 — 합리적 분포.

이 가중치를 person-time 데이터에 적용 → 가상 모집단에서 \(A \perp\!\!\!\perp L\) 회복.

Step 2: 가중 logistic hazard 모형 (MSM)

\[\mathrm{logit}\, \Pr(D_{k+1}^a = 0 | D_k^a = 0) = \beta_{0,k} + \beta_1 a + \beta_2 a k + \beta_3 a k^2\]

이 모형은 counterfactual hazard 의 모형. IPW 가중 logistic 회귀로 \(\beta\) 추정.

\(\beta_{0,k}\): 시간변동 절편 (예: \(\beta_0 + \beta_4 k + \beta_5 k^2\)).
\(\beta_1, \beta_2, \beta_3\): 처치 효과의 시간 의존성.

→ \(\widehat{\beta}\) 로 시점별 hazard 예측 → 곱셈으로 survival 도출.

직관 — MSM 의 의미: 일반 logistic 은 관측된 hazard 모형. IPW 가중 logistic 은 가상 모집단의 hazard 모형 = counterfactual hazard. 가중치가 \(L\)-편향 제거 → 모수가 인과 효과로 해석 가능.

직관 — 시간 의존성 모수: \(\beta_2 a k\) 는 처치 효과의 선형 시간 변화, \(\beta_3 a k^2\) 는 비선형 시간 변화. 단순 모형 (\(\beta_2 = \beta_3 = 0\)) 은 처치 효과가 모든 시점 동일 가정 — Cox PH 와 비슷. 시간 의존 product term 이 PH 가정 완화.

3.4 NHEFS IPW MSM 결과 (Hernan, Program 17.3)

NHEFS IPW MSM Survival 추정

가중 적합 후 시점별 survival 도출:

측정	값
120 month \(\widehat{\Pr}(D^{a=1}=0)\)	80.7%
120 month \(\widehat{\Pr}(D^{a=0}=0)\)	80.5%
차이	+0.2%
95% CI (bootstrap 500)	\((-4.1\%, +3.7\%)\)
추적 기간 최대 차이	\(-1.4\%\) (CI \(-3.4\%, +0.7\%\))

직관 — 보정 전후의 정반대 결론: - 보정 전: quitters 76% < non-quitters 82% → 흡연 중단이 사망 증가시키는 듯. - 보정 후: 차이 거의 0 → 흡연 중단의 직접 사망 영향 없음.

보정으로 흡수된 6% 차이는 나이 confounding — 늙은 사람이 더 많이 끊고 더 많이 사망. 보정 후 진짜 처치 효과는 거의 0.

직관 — 0% 효과의 의미: 흡연 중단이 사망에 영향 없다는 결론은 의학적으로 의외. 그러나 추적 기간이 10 년 — 흡연으로 인한 만성 질환 (COPD, 폐암) 의 효과 발현에는 짧을 수 있음. 더 긴 추적이면 다른 결과 가능. 추적 기간이 인과 결론에 본질적 영향.

직관 — 그러나 흡연 중단이 체중 증가시킴 (Ch.12): 같은 NHEFS 데이터에서 체중 증가는 분명히 보정 후 3.5kg 추정. 사망 영향은 약하지만 체중 영향은 강함. 결과 변수에 따라 처치 효과의 메커니즘이 다름 — 인과 분석은 결과 specific.

3.5 IPW MSM 의 가정 (전체 정리)

IPW MSM Survival 의 식별 가정

가정	영향
조건부 교환가능성 (\(Y^a \perp\!\!\!\perp A \| L\))	측정된 \(L\) 이 모든 처치-결과 confounder 포함
양의 확률	모든 \(L\) 에서 처치-비처치 모두 존재
일관성	\(A=a\) 의 사람의 \(D_k\) 가 \(D_k^a\) 와 일치
무작위 censoring (또는 \(L\) 로 보정)	행정 censoring 은 외생적 가정
Hazard 모형 specification	logistic 함수 형태 옳음
처치 모형 specification	propensity logistic 옳음

이 모든 가정 위에서 추정량이 일치.

직관 — 6 가지 가정의 누적: 단일 시점 인과 분석의 4 가정 + censoring 가정 + 두 모형 각각. 각 가정의 약한 위반이 결과에 누적적 영향. Sensitivity analysis 가 결과 해석의 필수.

4 비행정 Censoring 의 보정

Loss to Follow-up 보정 절차

비행정 censoring 이 있으면 추가 IPW 가중치:

\[SW_k^C = \frac{\Pr(C_{k+1} = 0 | C_k = 0, A)}{\Pr(C_{k+1} = 0 | C_k = 0, A, L)}\]

분모: \(L\) 보정한 censoring 미발생 확률.
분자: \(L\) 보정 안 한 marginal censoring 확률 (안정화).

총 가중치: \(SW^A \times \prod_{m=1}^{k} SW_m^C\).

각 시점 \(k\) 의 logistic 모형으로 \(\Pr(C_{k+1}=0 | \cdot)\) 추정.

직관 — 시점별 censoring 가중치: 매 시점 그 시점에서 censoring 안 될 조건부 확률 추정. 시점 1 은 1 의 가중치, 시점 2 는 1×2 의 누적, … 시점 \(k\) 는 \(k\) 개 곱. 누적 곱이라 시점이 길수록 가중치 분포 폭발 위험.

직관 — Stable IP 가중치 vs IP 가중치: 안정화 (\(SW^C\)) 는 분자에 \(L\)-비의존 함수 두어 가중치 변동 줄임. 시간변동 처치에서는 안정화가 사실상 필수 — 누적 곱이 너무 폭발하므로.

5 응용 분야

임상 코호트 사망 분석: NHEFS 같은 large cohort
종양학 long-term outcome: 무진행 생존, 전체 생존
HIV 코호트의 ART 효과: 사망까지 시간
고객 retention: subscription 이탈까지 시간
약물 안전성: 부작용 발생까지 시간 + dropout 보정

6 코드 — IPW MSM Survival

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

# Person-time 데이터 + 9 변수 confounder (NHEFS 가정)
nhefs = pd.read_csv("nhefs.csv").reset_index(drop=True)
nhefs["death"] = nhefs["death"].fillna(0).astype(int)
nhefs["survtime"] = np.where(
    nhefs["death"] == 1,
    (nhefs["yrdth"] - 1983) * 12 + nhefs["modth"],
    120
)

# === Step 1: SW^A 추정 ===
ps_formula = (
    "qsmk ~ sex + race + C(education) + age + I(age**2) "
    "+ smokeintensity + I(smokeintensity**2) + smokeyrs + I(smokeyrs**2) "
    "+ C(exercise) + C(active) + wt71 + I(wt71**2)"
)
ps_model = smf.logit(ps_formula, data=nhefs).fit(disp=False)
nhefs["ps"] = ps_model.predict()

p_a = nhefs["qsmk"].mean()
nhefs["sw_a"] = np.where(
    nhefs["qsmk"] == 1,
    p_a / nhefs["ps"],
    (1 - p_a) / (1 - nhefs["ps"])
)
print(f"SW^A: mean={nhefs['sw_a'].mean():.3f}, "
      f"min={nhefs['sw_a'].min():.3f}, max={nhefs['sw_a'].max():.3f}")

# === Step 2: Person-time 변환 ===
def to_person_time(df):
    rows = []
    for _, r in df.iterrows():
        n_periods = int(r["survtime"])
        for k in range(n_periods):
            row_dict = {
                "id": r["seqn"],
                "k": k,
                "qsmk": r["qsmk"],
                "sw_a": r["sw_a"],
                "D_k1": 1 if (k == n_periods - 1 and r["death"] == 1) else 0
            }
            rows.append(row_dict)
    return pd.DataFrame(rows)

person_time = to_person_time(nhefs)

# === Step 3: 가중 logistic hazard 모형 ===
hazard_formula = (
    "D_k1 ~ qsmk + qsmk:k + qsmk:I(k**2) + k + I(k**2)"
)
hazard_model = smf.glm(
    hazard_formula,
    data=person_time,
    family=sm.families.Binomial(),
    freq_weights=person_time["sw_a"]
).fit()
print(hazard_model.summary())

# === Step 4: Counterfactual survival 도출 ===
def counterfactual_survival(model, A_value, k_max=120):
    grid = pd.DataFrame({"qsmk": [A_value]*k_max, "k": range(k_max)})
    h = model.predict(grid).values
    s = np.cumprod(1 - h)
    return s

s1 = counterfactual_survival(hazard_model, 1)
s0 = counterfactual_survival(hazard_model, 0)
print(f"\nIPW MSM 120-month survival:")
print(f"  S(120 | A=1, c=0) = {s1[-1]:.3f}")
print(f"  S(120 | A=0, c=0) = {s0[-1]:.3f}")
print(f"  Difference = {s1[-1] - s0[-1]:.3f}")

# === Step 5: Bootstrap 95% CI ===
import multiprocessing as mp
from functools import partial

def bootstrap_iter(seed, df, formula_ps, formula_h, n_individuals):
    rng = np.random.default_rng(seed)
    boot_ids = rng.choice(df["seqn"].unique(), n_individuals, replace=True)
    boot = df[df["seqn"].isin(boot_ids)].copy()
    # 다시 PS 적합 + 가중 hazard 적합 + survival 도출 (생략)
    return None   # 실제 구현 필요

# 시점별 survival difference + 95% CI 시각화
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(range(120), s1, label="A=1 (Quit)")
ax.plot(range(120), s0, label="A=0 (Continue)")
ax.set_xlabel("Months"); ax.set_ylabel("IPW MSM Survival")
ax.legend()
plt.show()

7 Bootstrap 신뢰구간

500 개 bootstrap 표본 각각에 동일 절차 (PS 추정 → 가중치 → hazard 모형 → survival 곱셈) 적용. 500 개 추정값의 percentile 또는 normal 근사로 95% CI.

직관 — Bootstrap 의 자연스러움: IPW MSM 은 두 모형 (PS + hazard) 의 결합 → 분석적 분산 추정 어려움. Bootstrap 은 가정 없이 두 모형 모두 다시 적합 → 분포 직접 추정. 계산 비용 크지만 분산 추정의 표준.

직관 — Bootstrap 의 한계: 표본이 작거나 가중치 분포가 매우 비대칭이면 Bootstrap 결과도 불안정. 또한 Bootstrap CI 가 일부 영역 (예: hazard 0 근처) 에서 부정확. 표본 크고 가중치 안정적인 NHEFS 같은 데이터에서는 Bootstrap 신뢰성 높음.

8 g-method 도구의 비교

도구	처치 모형	결과 모형	가정 의존
IPW MSM (17.4)	적합	logistic hazard	처치 모형 옳음
g-formula (17.5)	적합 안 함	결과 모형 옳음	결과 모형 옳음
g-estimation (17.6)	적합	semiparametric SNMM	처치 모형 옳음

세 도구 모두 NHEFS 에서 비슷한 결론 (효과 약함) 줄 것이라 예상. 다음 글의 결과 비교.

직관 — 도구의 선택 결정: 결과 모형 더 자신 → g-formula. 처치 모형 더 자신 → IPW MSM. 두 모형 결합 → doubly robust. NHEFS 의 9 변수 + quadratic 모형은 두 모형 모두 합리적이라 어느 도구든 비슷한 결과 기대.

9 한 줄 요약

Censoring 은 단순 잔존 비율이 진짜 생존률을 향하지 못하게 하는 본질적 도전. Hazard 곱셈 공식 이 해결책 — 시점별 조건부 hazard 의 곱이 누적 생존. 행정 censoring 은 외생적, 비행정은 IPW 보정 필요. IPW MSM 은 가중 logistic hazard 모형으로 counterfactual survival 추정. NHEFS 흡연 중단의 사망 효과는 보정 후 거의 0 — confounding (나이) 이 단순 비교의 주요 원천. Bootstrap 이 분산 추정의 표준.

10 관련 주제

선행 지식

후속 주제

다른 카테고리 연결