16.8 생존분석과 인과추론 기초 · AI for Science 로드맵

이 장은 성격이 다른 두 주제를 묶는다. 공통점은 데이터가 질문에 대한 답을 직접 담고 있지 않다는 것이다. 수명 시험에서 시험이 끝날 때까지 고장 나지 않은 셀은 "수명이 최소 이만큼"이라는 불완전한 정보만 준다. 관찰 데이터에서 첨가제를 쓴 셀의 열화가 더 컸다는 사실은 "첨가제를 쓰면 열화가 커진다"는 뜻이 아닐 수 있다. 두 경우 모두 데이터를 표준 회귀에 그대로 넣으면 답이 나오지만, 그 답은 체계적으로 틀린다.

생존분석(survival analysis)은 관측이 중간에 끊긴 시간 데이터를 버리지 않고 쓰는 방법이고, 인과추론(causal inference)은 "개입하면 무엇이 바뀌는가"라는 질문을 관찰 데이터로 답할 수 있는 조건과 방법을 다룬다. 둘 다 한 장으로 끝낼 수 있는 분야가 아니므로, 이 장은 연구에서 가장 자주 마주치는 함정을 계산으로 확인하고 올바른 기본 도구를 익히는 데 집중한다.

학습 목표

censoring을 포함한 문제를 정의한다
상관과 인과를 구분한다
인과 그래프로 교란을 식별한다

생존분석

왜 필요한가. 배터리 셀 수백 개의 사이클 수명 시험을 생각한다. 시험 장비는 한정되어 있고 연구 기간도 정해져 있어서, 많은 셀이 고장 나기 전에 시험이 끝난다. 이런 셀의 수명은 "관측 시간 이상"이라는 것만 알 수 있다. 이것을 중도절단(censoring, 정확히는 오른쪽 중도절단)이라 한다. 같은 구조가 임상시험(연구 종료 시 생존한 환자), 기계 신뢰성(아직 고장 나지 않은 장비), 고객 이탈(아직 떠나지 않은 고객)에 나타난다.

정의. 수명을 확률변수 $T$ 라 하고, 셀 $i$ 마다 관측 시간 $t_i = \min(T_i, C_i)$ 와 고장 지표 $\delta_i = \mathbb{1}[T_i \le C_i]$ 를 얻는다. $C_i$ 는 중도절단 시점이다. 생존분석의 기본 양은 다음과 같다.

S(t) = P(T > t), \qquad h(t) = \lim_{\Delta t \to 0}\frac{P(t \le T < t + \Delta t \mid T \ge t)}{\Delta t} = \frac{f(t)}{S(t)}.

$S(t)$ 는 생존함수(survival function), $h(t)$ 는 위험함수(hazard function)로 " $t$ 까지 살아남은 것 중 바로 다음 순간 고장 날 비율"이다. $f(t)$ 는 수명의 밀도다. 누적 위험 $H(t) = \int_0^t h(u)\,du$ 에 대해 $S(t) = e^{-H(t)}$ 가 성립한다. 위험함수가 시간에 따라 커지면 노화(마모), 줄어들면 초기 불량이 걸러지는 과정이다. Weibull 분포 $S(t) = \exp(-(t/\lambda)^k)$ 는 형상 모수 $k > 1$ 이면 위험이 증가하고 $k < 1$ 이면 감소해 신뢰성 공학에서 널리 쓰인다.

거의 모든 표준 방법은 독립 중도절단(non-informative censoring)을 가정한다. 중도절단 시점이 수명과 무관해야 한다는 뜻이다. 장비 일정 때문에 시험을 끝낸 것은 이 가정을 만족하지만, "곧 고장 날 것 같아 보이는 셀을 미리 빼낸 것"은 만족하지 않는다.

중도절단을 무시하면. 셀 300개의 수명이 Weibull 분포(형상 3, 척도 1000 cycle)를 따르고, 셀마다 400–1200 cycle 사이의 임의 시점에 시험이 끝난다고 하자. 흔한 두 가지 처리, 즉 중도절단된 셀을 버리는 것과 시험 종료 시점을 수명으로 취급하는 것을 참값과 비교한다.

import numpy as np

rng = np.random.default_rng(6)
n, shape, scale = 300, 3.0, 1000.0                             # Weibull 수명 [cycle]
life = scale * rng.weibull(shape, n)                           # 참 수명 (실험에서는 일부만 보인다)
stop = rng.uniform(400, 1200, n)                               # 셀마다 실험이 끝난 시점 [cycle]
time = np.minimum(life, stop)                                  # 관측 시간
event = life <= stop                                           # True = 고장 관측, False = 중도절단
true_median = scale * np.log(2) ** (1 / shape)
print(f"고장 관측 {event.sum()}개, 중도절단 {(~event).sum()}개")
print(f"참 중앙 수명                    {true_median:6.0f} cycle")
print(f"중도절단을 버린 중앙값          {np.median(time[event]):6.0f} cycle")
print(f"중도절단을 고장으로 친 중앙값   {np.median(time):6.0f} cycle")

고장 관측 135개, 중도절단 165개
참 중앙 수명                       885 cycle
중도절단을 버린 중앙값             698 cycle
중도절단을 고장으로 친 중앙값      687 cycle

두 처리 모두 중앙 수명을 약 200 cycle(20% 이상) 과소추정한다. 중도절단된 셀은 대개 오래 버틴 셀이다. 이들을 버리면 일찍 고장 난 셀만 남고, 시험 종료 시점을 수명으로 치면 오래 버틴 셀의 수명을 깎는다. 어느 쪽이든 편향의 방향은 "수명이 짧다"이고, 중도절단 비율이 높을수록 커진다.

Kaplan–Meier 추정량. 중도절단된 셀도 "그 시점까지 살아 있었다"는 정보를 준다. Kaplan과 Meier(1958)의 곱-극한(product-limit) 추정량은 이 정보를 위험 집합(risk set)으로 쓴다. 고장이 관측된 시각 $t_{(1)} < t_{(2)} < \cdots$ 마다

\hat{S}(t) = \prod_{j:\,t_{(j)} \le t}\left(1 - \frac{d_j}{n_j}\right)

이다. $d_j$ 는 시각 $t_{(j)}$ 에 고장 난 수, $n_j$ 는 그 직전까지 고장도 중도절단도 되지 않은 수다. 중도절단된 셀은 중도절단 시점까지 분모 $n_j$ 에 들어가 기여하고 그 뒤로 빠진다. 각 단계의 $1 - d_j/n_j$ 는 "그 시각까지 살아남은 것이 그 시각을 넘길 조건부 확률"의 추정이고, 곱하면 생존함수가 된다.

def kaplan_meier(time, event):
    t_j = np.unique(time[event])                               # 고장이 관측된 시각들, (J,)
    d_j = np.array([np.sum((time == t) & event) for t in t_j]) # 그 시각의 고장 수
    n_j = np.array([np.sum(time >= t) for t in t_j])           # 그 직전까지 살아 있는(위험 집합) 수
    return t_j, np.cumprod(1 - d_j / n_j)                      # S(t_j)

def S_at(t_j, S, t):                                           # 계단 함수 S(t)
    k = np.searchsorted(t_j, t, side="right")
    return 1.0 if k == 0 else S[k - 1]

t_j, S = kaplan_meier(time, event)
print("   t    참 S(t)   KM   중도절단 버림   중도절단=고장")
for t in (400, 600, 800, 1000):
    drop = np.mean(time[event] > t)                            # 고장 관측만으로 센 생존 비율
    as_fail = np.mean(time > t)
    print(f"{t:5d}   {np.exp(-(t / scale) ** shape):.3f}   {S_at(t_j, S, t):.3f}       {drop:.3f}          {as_fail:.3f}")
print(f"KM 중앙 수명 {t_j[np.argmax(S <= 0.5)]:.0f} cycle (참값 {true_median:.0f})")

   t    참 S(t)   KM   중도절단 버림   중도절단=고장
  400   0.938   0.943       0.874          0.943
  600   0.806   0.840       0.674          0.643
  800   0.599   0.593       0.319          0.313
 1000   0.368   0.308       0.067          0.087
KM 중앙 수명 868 cycle (참값 885)

Kaplan–Meier 곡선은 모든 시각에서 참 생존함수를 표본 변동 범위 안에서 따라가고, 중앙 수명도 868 cycle로 참값 885에 가깝다. 중도절단을 버리거나 고장으로 친 추정은 800 cycle에서 생존 비율을 절반 정도로, 1000 cycle에서는 4분의 1 이하로 과소추정한다. 다만 Kaplan–Meier도 관측 기간 밖은 말해 주지 않는다. 이 예에서 가장 늦은 시험 종료가 1200 cycle이므로 그 이후의 $S(t)$ 는 추정할 수 없다.

Cox 비례위험 모델. 온도나 충전 속도 같은 공변량 $\mathbf{x}$ 가 수명에 미치는 영향을 보려면 회귀 모델이 필요하다. Cox(1972)의 비례위험 모델은

h(t \mid \mathbf{x}) = h_0(t)\exp(\boldsymbol{\beta}^\top\mathbf{x})

로 쓴다. 기저 위험 $h_0(t)$ 의 모양은 정하지 않고, 공변량이 위험을 일정한 배율 $e^{\beta_k}$ (위험비, hazard ratio)로 곱한다고만 가정한다. 계수는 고장이 일어난 시각마다 "위험 집합 중 하필 그 셀이 고장 날 조건부 확률"을 곱한 부분 likelihood

L(\boldsymbol{\beta}) = \prod_{i:\,\delta_i = 1}\frac{\exp(\boldsymbol{\beta}^\top\mathbf{x}_i)}{\sum_{j:\,t_j \ge t_i}\exp(\boldsymbol{\beta}^\top\mathbf{x}_j)}

를 최대화해 추정한다. $h_0(t)$ 가 분자와 분모에서 약분되므로 추정할 필요가 없다. 핵심 가정은 위험비가 시간에 따라 일정하다는 비례위험성이다. 두 집단의 Kaplan–Meier 곡선이 교차하거나, 시간에 따라 계수가 달라지는 잔차 패턴이 보이면 가정을 의심한다.

모수적 대안으로 가속수명 모델(AFT, accelerated failure time)은 $\log T = \mathbf{b}^\top\mathbf{x} + \varepsilon$ 처럼 공변량이 수명의 시간 척도를 늘리거나 줄인다고 본다. 계수의 해석이 "수명이 $e^{b}$ 배"로 직관적이다. Weibull 분포는 비례위험 모델과 AFT 모델 양쪽에 속하는 유일한 경우로, 형상 모수가 $k$ 이면 AFT 계수 $b$ 와 위험비 계수 $\beta$ 사이에 $\beta = -k\,b$ 가 성립한다.

수명 예측과의 연결

왜 필요한가. 배터리·베어링·터빈 같은 장비의 잔여 수명(RUL, remaining useful life) 예측은 공학 ML에서 가장 흔한 과제 중 하나다. 가장 단순한 정식화는 "고장까지 남은 사이클 수"를 목표로 하는 회귀다. 그런데 회귀의 목표값을 만들려면 고장 시점을 알아야 하므로, 학습 데이터에는 고장까지 관측된 장비만 들어가게 된다. 이것은 앞에서 본 "중도절단을 버리는" 처리와 같다.

정식화	목표	중도절단된 장비	결과
회귀 (고장 장비만)	고장까지 남은 시간	버림	일찍 고장 나는 장비 쪽으로 편향, 수명과 공변량의 관계도 약화
회귀 (중도절단 = 고장)	관측 시간	수명으로 취급	오래 가는 장비의 수명을 깎음
생존분석	생존함수·위험	"최소 이만큼"으로 likelihood에 기여	독립 중도절단 가정 아래 일관된 추정
고정 horizon 분류	" $H$ 사이클 안에 고장 나는가"	$H$ 전에 끊기면 label 불명	단순하지만 $H$ 전 중도절단 장비 처리가 다시 문제

공변량의 효과를 추정할 때 차이가 분명해진다. 온도가 수명에 미치는 효과를 세 방법으로 추정한다. 참 모델은 온도 1 °C 상승마다 log 수명이 0.03 줄어드는(10 °C 상승 시 수명 0.74배) Weibull AFT 모델이고, 셀마다 300–1500 cycle 사이의 임의 시점에 시험이 끝난다. 생존분석 방법은 고장 셀에는 밀도 $\log f(t_i)$ 를, 중도절단 셀에는 생존확률 $\log S(t_i)$ 를 likelihood에 넣어 최대화한다.

from scipy.optimize import minimize

temp_c = rng.uniform(15, 45, n)                                # 운전 온도 [°C]
b_true = -0.03                                                 # 1 °C당 log 수명 변화
life2 = np.exp(np.log(1500) + b_true * (temp_c - 25)) * rng.weibull(shape, n)
stop2 = rng.uniform(300, 1500, n)
t2, e2 = np.minimum(life2, stop2), life2 <= stop2              # 관측 시간, 고장 여부

def neg_loglik(th):                                            # Weibull AFT: log λ = a + b (온도 - 25)
    a, b, log_k = th
    k, lam = np.exp(log_k), np.exp(a + b * (temp_c - 25))
    z = (t2 / lam) ** k
    log_f = np.log(k) - np.log(lam) + (k - 1) * np.log(t2 / lam) - z    # 고장: 밀도
    return -np.sum(np.where(e2, log_f, -z))                    # 중도절단: 생존확률 log S = -z

b_drop = np.polyfit(temp_c[e2] - 25, np.log(t2[e2]), 1)[0]    # 고장만으로 log 수명 회귀
b_fail = np.polyfit(temp_c - 25, np.log(t2), 1)[0]            # 중도절단도 고장으로 취급
b_mle = minimize(neg_loglik, x0=[np.log(t2.mean()), 0.0, 0.0], method="Nelder-Mead",
                 options={"xatol": 1e-6, "fatol": 1e-6, "maxiter": 4000}).x[1]
print(f"중도절단 {np.mean(~e2):.0%}.  10 °C 상승 시 수명 배율 exp(10 b)")
for name, b in [("참값", b_true), ("고장만 회귀", b_drop), ("중도절단=고장 회귀", b_fail), ("중도절단 포함 Weibull MLE", b_mle)]:
    print(f"  {name:22s} b = {b:+.3f}   배율 {np.exp(10 * b):.2f}")

중도절단 68%.  10 °C 상승 시 수명 배율 exp(10 b)
  참값                     b = -0.030   배율 0.74
  고장만 회귀                 b = -0.002   배율 0.98
  중도절단=고장 회귀             b = -0.010   배율 0.90
  중도절단 포함 Weibull MLE    b = -0.028   배율 0.75

고장 셀만으로 회귀하면 온도 효과가 거의 사라진다(10 °C당 수명 0.98배). 고온 셀은 빨리 고장 나 대부분 데이터에 남지만, 저온 셀은 오래 버티는 것들이 중도절단으로 빠지고 일찍 고장 난 것만 남는다. 그래서 저온 셀의 "관측된 수명"이 고온 셀과 비슷해 보인다. 중도절단을 고장으로 취급해도 효과가 크게 약해진다(0.90배). 중도절단을 likelihood에 넣은 Weibull 최대우도는 0.75배로 참값 0.74배를 되찾는다. 중도절단을 무시한 회귀는 수명의 수준뿐 아니라 공변량의 효과까지 약화시킨다. "온도는 수명에 별 영향이 없다"는 잘못된 결론으로 이어질 수 있다.

연구에서의 의미. RUL 모델을 만들 때는 데이터에 중도절단된 장비가 얼마나 있는지, 왜 중도절단되었는지(시험 종료, 다른 이유의 제거, 데이터 수집 중단)부터 기록한다. 그리고 생존분석 모델(Cox, AFT, 또는 생존분석용 트리 앙상블과 신경망)을 baseline에 포함한다. 중도절단 셀을 버린 회귀 모델이 시험 데이터에서 좋은 점수를 받더라도, 그 시험 데이터 역시 고장 장비만 모은 것이라면 같은 편향을 공유해 문제가 드러나지 않는다.

인과추론 기초

왜 필요한가. 지금까지의 모델은 모두 예측을 했다. $\mathbf{x}$ 를 보면 $y$ 가 얼마일지 맞히는 일이다. 그러나 연구 질문의 상당수는 개입에 관한 것이다. "첨가제를 넣으면 열화가 줄어드는가", "온도를 10 °C 낮추면 수명이 얼마나 늘어나는가". 예측을 잘하는 모델이 개입의 효과를 맞힌다는 보장은 없다.

상관과 인과. 관찰 데이터에서 $T$ 와 $Y$ 가 함께 변하는 이유는 여러 가지다. $T$ 가 $Y$ 의 원인일 수도, $Y$ 가 $T$ 의 원인일 수도, 둘 다 제3의 변수 $Z$ 의 결과일 수도 있다. 마지막 경우의 $Z$ 를 교란 변수(confounder)라 한다. 교란이 있으면 상관의 크기는 물론 부호도 인과 효과와 다를 수 있다.

인과 DAG. 변수 사이의 인과 관계를 화살표로 그린 방향성 비순환 그래프(DAG, directed acyclic graph)로 가정을 명시한다. 이 장의 예는 다음과 같다.

        Z (운전 온도)
       ↙          ↘
T (첨가제) ──────→ Y (용량 감소)

더운 현장에서 첨가제를 더 많이 쓰고( $Z \to T$ ), 더운 현장의 셀은 더 빨리 열화한다( $Z \to Y$ ). 우리가 알고 싶은 것은 $T \to Y$ 화살표의 크기다. $T$ 와 $Y$ 사이의 경로 중 $T \leftarrow Z \to Y$ 처럼 $T$ 로 들어오는 화살표로 시작하는 경로를 뒷문 경로(backdoor path)라 한다. 뒷문 경로가 열려 있으면 $T$ 와 $Y$ 의 상관에 인과가 아닌 성분이 섞인다.

backdoor 기준. 변수 집합 $\mathcal{Z}$ 가 (1) $T$ 에서 $Y$ 로 가는 모든 뒷문 경로를 막고, (2) $T$ 의 후손을 포함하지 않으면, $\mathcal{Z}$ 로 보정해 인과 효과를 얻을 수 있다(Pearl, 1995).

P(Y \mid do(T = t)) = \sum_{z} P(Y \mid T = t, \mathcal{Z} = z)\,P(\mathcal{Z} = z).

$do(T = t)$ 는 $T$ 를 관찰한 것이 아니라 외부에서 $t$ 로 정했다는 뜻이다(다음 절). 우변은 관찰 데이터로 계산할 수 있는 양이다. $\mathcal{Z}$ 의 층마다 처치군과 대조군을 비교한 뒤, 층의 분포로 평균한다.

보정 방법은 여러 가지다. 회귀 보정은 $Y$ 를 $T$ 와 $\mathcal{Z}$ 에 회귀해 $T$ 의 계수를 읽는다(모델이 맞아야 한다). 성향점수(propensity score) $e(z) = P(T = 1 \mid \mathcal{Z} = z)$ 는 $\mathcal{Z}$ 대신 이 스칼라 하나로 보정해도 충분하다는 성질이 있다(Rosenbaum, Rubin, 1983). 처치군을 $1/e$ , 대조군을 $1/(1-e)$ 로 가중하면(IPW, inverse probability weighting) 두 군의 $\mathcal{Z}$ 분포가 전체 모집단과 같아진다. 셀 5000개의 관찰 데이터에서 세 방법과, 같은 현장에서 첨가제를 무작위로 배정한 경우를 비교한다. 참 효과는 용량 감소 $-1.00$ %다.

m = 5000
Z = rng.normal(30, 7, m)                                       # 운전 온도 [°C] (교란 변수)
T = rng.uniform(size=m) < 1 / (1 + np.exp(-0.3 * (Z - 30)))    # 더운 현장일수록 첨가제를 더 쓴다
Y = 0.2 * Z - 1.0 * T + rng.normal(0, 0.5, m)                  # 용량 감소 [%]: 첨가제의 참 효과 -1.0

def logistic_fit(X, t, iters=25):                              # Newton 방법, X: (m, p)
    w = np.zeros(X.shape[1])
    for _ in range(iters):
        p = 1 / (1 + np.exp(-X @ w))
        w += np.linalg.solve(X.T @ (X * (p * (1 - p))[:, None]), X.T @ (t - p))
    return 1 / (1 + np.exp(-X @ w))

naive = Y[T].mean() - Y[~T].mean()
reg = np.linalg.lstsq(np.c_[np.ones(m), T, Z], Y, rcond=None)[0][1]   # Y ~ 1 + T + Z의 T 계수
e = logistic_fit(np.c_[np.ones(m), Z], T.astype(float))         # 성향점수 P(T=1 | Z)
w1, w0 = T / e, ~T / (1 - e)                                  # 역확률 가중치
ipw = np.sum(w1 * Y) / np.sum(w1) - np.sum(w0 * Y) / np.sum(w0)
T_rct = rng.uniform(size=m) < 0.5                              # 같은 현장에서 무작위로 배정했다면
Y_rct = 0.2 * Z - 1.0 * T_rct + rng.normal(0, 0.5, m)
print(f"첨가제 사용 셀의 평균 온도 {Z[T].mean():.1f} °C, 미사용 {Z[~T].mean():.1f} °C")
for name, v in [("단순 비교", naive), ("회귀 보정 (Z 포함)", reg), ("성향점수 가중(IPW)", ipw),
                ("무작위 배정 실험의 단순 비교", Y_rct[T_rct].mean() - Y_rct[~T_rct].mean())]:
    print(f"  {name:26s} 효과 추정 {v:+.2f} %   (참값 -1.00)")

첨가제 사용 셀의 평균 온도 34.3 °C, 미사용 25.6 °C
  단순 비교                      효과 추정 +0.72 %   (참값 -1.00)
  회귀 보정 (Z 포함)               효과 추정 -1.02 %   (참값 -1.00)
  성향점수 가중(IPW)               효과 추정 -0.89 %   (참값 -1.00)
  무작위 배정 실험의 단순 비교           효과 추정 -0.96 %   (참값 -1.00)

단순 비교는 첨가제가 용량 감소를 0.72% 늘린다고 말한다. 부호까지 틀렸다. 첨가제를 쓴 셀의 평균 운전 온도가 8.7 °C 높아서, 온도에 의한 열화 차이가 첨가제의 효과를 덮었기 때문이다. 온도로 보정한 회귀(−1.02)와 성향점수 가중(−0.89)은 참값 근처로 돌아온다. IPW가 조금 더 벗어난 것은 아주 덥거나 추운 현장에서 성향점수가 0이나 1에 가까워 일부 셀의 가중치가 매우 커지고, 그만큼 분산이 크기 때문이다. 무작위 배정 실험에서는 단순 비교만으로 −0.96을 얻는다.

이 결과는 두 가지 가정 위에 있다. 첫째, 교란 변수를 모두 관측했다(측정되지 않은 교란이 없다). 운전 온도 외에 현장의 충전 습관이 첨가제 사용과 열화에 함께 영향을 준다면, 온도만 보정해서는 편향이 남는다. 이 가정은 데이터로 검증할 수 없고, 도메인 지식으로 DAG를 정당화해야 한다. 둘째, 양의 확률(positivity)이다. 모든 $\mathcal{Z}$ 값에서 처치와 대조가 모두 일어날 수 있어야 한다. 45 °C 현장에서 첨가제를 쓰지 않은 셀이 하나도 없다면, 그 조건에서의 효과는 데이터가 아니라 모델의 외삽으로만 추정된다.

개입과 반사실

do 연산의 직관. $P(Y \mid T = 1)$ 은 "첨가제를 쓴 것으로 관찰된 셀들"의 결과 분포이고, $P(Y \mid do(T = 1))$ 은 "모든 셀에 첨가제를 넣었다면"의 결과 분포다. 앞의 것은 첨가제를 쓰게 된 이유(더운 현장)를 함께 담고, 뒤의 것은 $T$ 로 들어오는 화살표를 모두 끊은 그래프에서의 분포다. 인과 효과는 $E[Y \mid do(T=1)] - E[Y \mid do(T=0)]$ 로 정의한다.

반사실. 셀 하나에 대해 "첨가제를 넣었을 때의 결과" $Y(1)$ 과 "넣지 않았을 때의 결과" $Y(0)$ 를 생각할 수 있다. 이를 잠재 결과(potential outcomes)라 한다. 실제로는 둘 중 하나만 관측되고, 나머지는 반사실(counterfactual)로 영원히 볼 수 없다. 개별 효과 $Y(1) - Y(0)$ 는 직접 관측할 수 없으므로, 인과추론은 보통 집단 평균 $E[Y(1) - Y(0)]$ 를 추정한다. 이 관점과 DAG 관점은 서로 대응하며, 실무 교과서로는 Hernán, Robins(2020)가 두 관점을 함께 다룬다.

무작위 실험. 처치를 무작위로 배정하면 $T$ 는 어떤 변수와도 무관해지므로 $T$ 로 들어오는 화살표가 모두 사라진다. 관측했든 하지 않았든 모든 교란이 평균적으로 균형을 이루므로, 단순 비교가 인과 효과의 불편 추정이 된다. 앞 셀의 무작위 배정 결과가 이것이다. 실험을 설계할 수 있다면 무작위 배정이 가장 강한 근거다.

관찰 데이터의 한계: 보정하면 안 되는 변수. "변수를 많이 넣을수록 보정이 잘 된다"는 직관은 틀리다. backdoor 기준의 두 번째 조건이 말하듯 처치의 결과를 보정하면 안 된다. 두 가지 경우가 대표적이다.

매개변수(mediator): $T \to M \to Y$ . 첨가제가 SEI 저항 증가( $M$ )를 줄여 열화를 줄인다면, $M$ 을 보정하는 순간 첨가제가 효과를 내는 경로를 막아 "효과 없음"을 얻는다.
충돌변수(collider): $T \to C \leftarrow Y$ . 품질 검사를 통과한 셀만 데이터에 남는데 통과 여부가 첨가제 사용과 열화에 모두 달려 있다면, 통과한 셀만 분석하는 것 자체가 $C$ 로 보정한 것이고 $T$ 와 $Y$ 사이에 가짜 연관을 만든다.

매개변수:  T ──→ M ──→ Y          충돌변수:  T ──→ C ←── Y

M = 2.0 - 0.8 * T + rng.normal(0, 0.3, m)                      # 저항 증가 [mΩ]: 첨가제 → M → Y
Y2 = 0.2 * Z + 1.25 * M + rng.normal(0, 0.5, m)                # 첨가제의 총효과 = 1.25 × (-0.8) = -1.0
S = Y2 < np.quantile(Y2, 0.5) + 1.5 * T                        # 선별: 첨가제 셀은 느슨한 기준으로 통과

def coef_T(*cols, y=Y2, rows=slice(None)):                     # y ~ 1 + T + cols의 T 계수
    X_ = np.c_[np.ones(m), T, *cols][rows]
    return np.linalg.lstsq(X_, y[rows], rcond=None)[0][1]

for name, v in [("Z만 보정 (올바름)", coef_T(Z)), ("Z와 매개변수 M 보정", coef_T(Z, M)),
                ("Z 보정, 선별을 통과한 셀만", coef_T(Z, rows=S))]:
    print(f"  {name:26s} 효과 추정 {v:+.2f} %   (총효과 참값 -1.00)")

  Z만 보정 (올바름)                효과 추정 -0.99 %   (총효과 참값 -1.00)
  Z와 매개변수 M 보정               효과 추정 -0.01 %   (총효과 참값 -1.00)
  Z 보정, 선별을 통과한 셀만           효과 추정 -0.54 %   (총효과 참값 -1.00)

온도만 보정하면 총효과 −1을 되찾는다(−0.99). 매개변수 $M$ 까지 보정하면 −0.01로 효과가 사라진다. 이 예에서 첨가제의 효과는 전부 저항 증가를 통해 나타나므로, $M$ 을 고정하면 첨가제가 할 일이 없다. 연구 질문이 "저항 증가 외의 경로로 효과가 있는가"(직접 효과)라면 이것이 맞는 계산이지만, "첨가제를 넣으면 열화가 줄어드는가"(총효과)에는 틀린 답이다. 선별을 통과한 셀만 분석하면 −0.54로 효과가 절반 가까이 줄어든다. 첨가제 셀은 느슨한 기준으로 통과하므로, 통과한 첨가제 셀에는 열화가 큰 셀도 섞이고 통과한 미사용 셀은 열화가 작은 셀만 남는다.

특징 중요도는 인과가 아니다. 16.4장의 permutation importance나 트리의 impurity 중요도는 "이 특징을 섞으면 예측이 얼마나 나빠지는가"를 잰다. 교란 변수의 대리 지표, 결과의 후손(예: 열화 뒤에 측정한 저항), 충돌변수로 선별된 데이터의 특징은 모두 예측에 매우 중요할 수 있지만, 그 변수를 바꾼다고 결과가 바뀌지 않는다. 특징 중요도는 "모델이 무엇을 쓰는가"에 대한 진단이고, "무엇을 바꾸면 결과가 바뀌는가"에 대한 답은 DAG에 근거한 보정이나 실험에서 나온다.

실습

실험이 중단된 셀을 포함한 수명 데이터를 단순 회귀와 생존분석으로 각각 모델링해 결과를 비교한다. 코드는 labs/classical-ml/survival-and-causal/censored_lifetimes.py다. 셀 120개의 사이클 수명이 온도와 C-rate에 따라 달라지는 Weibull AFT 모델을 따르고(온도 계수 −0.03, 형상 3), 셀마다 200–1200 cycle 사이의 임의 시점에 시험이 끝난다. 같은 실험을 200번 반복해 네 방법의 온도 계수 추정값 분포를 비교한다. Cox 모델은 위험비 계수 $\beta$ 를 추정하므로, Weibull 관계 $b = -\beta/k$ 로 AFT 계수에 맞춰 환산했다( $k$ 는 AFT 추정값). 집필 시점의 실행 결과는 다음과 같다.

셀 120개 × 반복 200번, 중도절단 비율 평균 58%
온도 계수 b (참값 -0.030, 10 °C 상승 시 수명 배율 0.741)
  방법                    b 평균    b 표준편차   평균 배율   |b - 참값| ≤ 2sd인 반복 비율
  고장만 회귀                -0.0149    0.0076      0.862      0.53
  중도절단=고장 회귀            -0.0121    0.0046      0.886      0.03
  Weibull AFT MLE       -0.0298    0.0061      0.743      0.95
  Cox (b로 환산)           -0.0296    0.0062      0.744      0.96

하나의 실험(첫 반복과 같은 시드)에서 45 °C, 1C 셀의 중앙 수명 예측 [cycle]
  고장만 회귀                  417   (log 수명 평균의 역변환)
  중도절단=고장 회귀              442   (log 수명 평균의 역변환)
  Weibull AFT MLE         592
  참값                      541

결과를 읽는 법은 다음과 같다.

단순 회귀는 온도 효과를 절반 이하로 추정한다. 200번 반복의 평균이 고장만 쓰면 −0.015, 중도절단을 고장으로 치면 −0.012로 참값 −0.030의 절반 이하다. 셀을 더 모아도 사라지지 않는 체계적 편향이다. "중도절단=고장" 회귀는 표준편차가 작아 정밀해 보이지만, 참값이 추정값의 ±2sd 안에 드는 반복이 3%뿐이다. 정밀하지만 틀린 추정이다.
두 생존분석 방법은 편향 없이 참값을 되찾는다. Weibull AFT와 Cox 모두 평균 −0.030 근처이고, 참값이 ±2sd 안에 드는 비율이 약 95%다. Cox 모델은 기저 위험의 모양을 가정하지 않고도 AFT와 거의 같은 정밀도를 냈다. 이 예에서는 Weibull 가정이 맞으므로 AFT가 유리할 이유가 있었는데도 그렇다.
수명 예측값도 달라진다. 한 실험에서 45 °C, 1C 셀의 중앙 수명을 단순 회귀는 417–442 cycle로, Weibull AFT는 592 cycle로 예측했다(참값 541). 생존분석 쪽의 오차는 셀 120개짜리 실험 하나의 표본 변동이고, 단순 회귀의 오차는 방향이 정해진 편향이다. 보증 기간이나 교체 주기를 단순 회귀로 정하면 체계적으로 보수적인(짧은) 값을 얻는다.

이어서 두 연습문제를 푼다. 첫째는 Kaplan–Meier 추정량과 중앙 수명을 구현하고, 둘째는 로지스틱 회귀로 성향점수를 적합해 처치 효과를 역확률 가중으로 추정한다.

실패 사례

중단된 실험을 버리거나 수명이 끝난 것으로 취급하는 것. 이 장의 예에서 두 처리 모두 중앙 수명을 20% 이상 과소추정했고, 온도의 효과를 절반 이하로 약화시켰다. 중도절단 여부와 사유를 기록하고, Kaplan–Meier·Cox·AFT 같은 생존분석 방법을 쓴다.
중도절단이 수명과 독립인지 확인하지 않는 것. 성능이 떨어져 보이는 셀을 일찍 빼냈다면 독립 중도절단 가정이 깨져 생존분석도 편향된다. 시험 종료 규칙을 사전에 정하고 기록한다.
특징 중요도를 인과로 해석하는 것. 중요도는 예측에 대한 기여이지 개입의 효과가 아니다. 교란 변수의 대리 지표나 결과의 후손도 중요도가 높게 나온다. 개입에 대한 주장은 DAG와 보정, 또는 실험으로 뒷받침한다.
관찰 데이터의 단순 비교로 처치 효과를 주장하는 것. 이 장의 예에서 단순 비교는 효과의 부호까지 틀렸다. 처치 배정에 영향을 준 변수를 DAG로 정리하고 보정한다.
"변수를 많이 넣을수록 좋다"며 처치의 결과를 보정하는 것. 매개변수를 보정하면 효과가 사라지고(−1.00 → −0.01), 충돌변수로 선별하면 효과가 왜곡된다(−0.54). 보정 집합은 backdoor 기준으로 고른다.

통과 기준

자신의 연구 문제에서 중도절단이 생기는 지점(시험 종료, 탈락, 데이터 수집 중단)과 그 사유를 식별하고, 독립 중도절단 가정이 성립하는지 논증한 뒤 생존분석 방법으로 수명이나 위험을 추정한다.
자신의 연구 질문에 대한 인과 DAG를 그려 교란 변수·매개변수·충돌변수를 구분하고, backdoor 기준을 만족하는 보정 집합과 남는 가정(측정되지 않은 교란, 양의 확률)을 제시한다.
연습문제 두 개(kaplan-meier, ipw-effect)의 테스트를 통과한다.

참고 문헌

Kaplan, Meier (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282), 457–481. — Kaplan–Meier 곱-극한 추정량
Cox (1972). Regression models and life-tables. Journal of the Royal Statistical Society, Series B, 34(2), 187–220. — 비례위험 모델과 부분 likelihood
Rosenbaum, Rubin (1983). The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1), 41–55. — 성향점수
Pearl (1995). Causal diagrams for empirical research. Biometrika, 82(4), 669–710. — 인과 그래프, backdoor 기준, do 연산
Hernán, Robins (2020). Causal Inference: What If. Chapman & Hall/CRC. — 잠재 결과와 인과 그래프를 함께 다루는 공개 교과서