13강 / 20강

moderated t 와 p값 — 함수 하나로 묶기

stats.t.sf, 양측검정, 그리고 def

수축된 분산 s2p 를 얻었으니 t 를 다시 만듭니다. 이것을 moderated t(조절된 t) 라고 불러요.

t = lfc / np.sqrt(s2p * (1/na + 1/nb)) dft = d0 + d # 자유도도 빌려온 만큼 늘어난다 p = 2 * stats.t.sf(np.abs(t), dft)

자유도가 d(5)에서 d0 + d(≈7.65)로 늘어난 점을 보세요. 이웃 유전자에게서 정보를 빌려왔으니, 그만큼 더 많은 관측을 한 셈으로 쳐 주는 겁니다. 이것이 limma 가 적은 반복에서도 강한 이유입니다.

p값 읽는 법

p값은 "이 유전자가 변했을 확률"이 아닙니다. 정반대에 가까워요.

p = 정말 아무 차이도 없다면, 이만큼 큰 t 가 우연히 나올 확률

p 가 작다는 것은 "우연으로 설명하기 힘들다"는 뜻입니다. 그래서 우리는 차이가 없다는 가정(귀무가설) 을 의심하게 됩니다.

고등학생p 값이 나오는 자리

드디어 p 값이 나옵니다. /stat 9강에서 배운 그것이에요.

p 값 = 귀무가설이 참일 때 이만한 t 가 나올 확률

여기서 귀무가설은

"이 유전자는 두 군에서 차이가 없다"

그러니 p 가 작다는 것은 "차이가 없다고 보기엔 이 결과가 너무 드물다" 입니다.

퀴즈가 이걸 묻습니다

"p = 0.001 을 바르게 읽은 것은?"

해석	맞나
"차이가 없을 확률 0.1%"	✘ — 조건을 뒤집었다
"차이가 있을 확률 99.9%"	✘
"차이가 없다면 이런 결과가 나올 확률 0.1%"	✔

/stat 9강의 그 함정 — 조건을 뒤집지 마라

그리고 여기가 특별합니다

보통의 t 검정이 아니라 moderated t 를 씁니다. 12강에서 분산을 빌려 왔으니까요.

	보통 t	moderated t
분산	이 유전자 것만	전체와 섞은 것
자유도	d	d₀ + d (크다)
안정성	불안정	안정

같은 자료인데 더 믿을 만한 p 값이 나온다

고등학생sf 와 2를 곱하는 이유

note 가 양측검정을 설명합니다. 코드에 2 * 가 붙어 있는 이유예요.

p = 2 × stats.t.sf(|t|, df)

sf 가 무엇인가

survival function — "이보다 클 확률" 입니다. 꼬리의 넓이예요.

왜 2를 곱하나

우리가 묻는 것이 "달라졌나" 이지 "늘었나" 나 "줄었나" 가 아니기 때문입니다.

물음	꼬리	계산
"늘었나?"	오른쪽만	sf
"달라졌나?"	양쪽	2 × sf

양쪽 꼬리를 다 더해야 "달라졌다" 를 잰다

비유 — 과녁에서 벗어남

화살이 과녁 중심에서 벗어났는지 봅니다. 왼쪽으로 벗어나도, 오른쪽으로 벗어나도 벗어난 것이죠.

한쪽만 세면 절반을 놓칩니다.

언제 한쪽만 보나

방향을 미리 정했을 때만요. "이 처리는 반드시 늘린다" 를 실험 전에 정한 경우요.

결과를 보고 방향을 정하면 안 된다

/stat 11강의 p 해킹입니다. 단측검정은 p 가 절반이 되니 유혹이 커요.

이 강의는 양측을 씁니다

이 논문은 유전자 전체를 훑습니다. 어느 방향으로 변할지 모르는 상태죠.

그러니 양측이 맞습니다. 정직한 선택이에요.

🎮 t 분포와 p값 면적

t 값과 자유도를 움직여 보세요. 꼬리 면적(= p값)이 얼마나 빨리 줄어드는지, 자유도가 커지면 곡선이 어떻게 날렵해지는지 보입니다.

지금까지의 모든 것을 함수 하나로

10~13강에서 배운 것을 코드랩은 함수 하나에 담았습니다. 이제 이 함수의 모든 줄을 읽을 수 있습니다.

코드랩의 그 함수 그대로입니다. 이제 한 줄도 낯설지 않습니다.

Python 브라우저에서 실행

import numpy as np
from scipy import stats
from scipy.special import digamma, polygamma
from scipy.optimize import brentq

def moderated_ttest(M, ca, cb):
    a, b = M[:, ca], M[:, cb]                         # 두 그룹 자르기
    na, nb = a.shape[1], b.shape[1]
    lfc = b.mean(1) - a.mean(1)                       # 10강 · 신호
    d = (na - 1) + (nb - 1)                           # 10강 · 자유도
    s2 = np.maximum((a.var(1, ddof=1)*(na-1) + b.var(1, ddof=1)*(nb-1)) / d, 1e-8)
    e = np.log(s2) - digamma(d/2) + np.log(d/2)       # 12강 · 사전분포 적합
    tri = np.var(e, ddof=1) - polygamma(1, d/2)
    if tri <= 0:                                      # 유전자 간 차이가 없다면
        d0, s0 = np.inf, np.exp(e.mean())             #   모두 같은 분산을 쓴다
    else:
        d0 = brentq(lambda x: polygamma(1, x/2) - tri, 1e-3, 1e6)
        s0 = np.exp(e.mean() + digamma(d0/2) - np.log(d0/2))
    if np.isinf(d0):
        s2p, dft = np.full_like(s2, s0), np.full_like(s2, 1e6)
    else:
        s2p, dft = (d0*s0 + d*s2) / (d0 + d), np.full_like(s2, d0 + d)
    t = lfc / np.sqrt(s2p * (1/na + 1/nb))            # 13강 · moderated t
    return lfc, t, 2 * stats.t.sf(np.abs(t), dft)     # 13강 · 양측 p값

print("moderated_ttest 준비 완료")

고등학생함수 하나로 묶기

코드가 지금까지 배운 걸 함수 하나로 묶습니다.

12강의 수축 + 13강의 p 값 = 함수 하나

함수로 묶는 이유

/py 7강에서 배운 그것입니다.

이유	여기서는
덜 쓴다	20줄을 한 줄로 부른다
덜 틀린다	고칠 때 한 군데만
읽기 쉽다	`moderated_t(a, b)` 는 뜻이 분명

캡션이 말한다 — "이제 한 줄도 낯설지 않습니다"

정말 그런지 확인해 보세요

함수 안의 줄들을 하나씩 짚어 보면 —

줄	어느 강
평균의 차이	10강
분산 계산 (ddof=1)	10강
합동분산	10강
수축 (s₀ · d₀)	12강
brentq	12강
t 계산	13강
2 × sf	13강

네 강에서 배운 것이 한 함수에 모였다

이게 이 강의의 구조입니다

한 강에 한 조각씩 배우고, 모아서 도구를 만들고, 그 도구로 본론을 합니다.

14강부터가 본론이에요 — 결과를 해석하는 부분입니다.

고등학생if tri <= 0 은 왜 있나

note 가 가장자리 처리를 설명합니다. /py 7강의 그 이야기예요.

수학적으로 불가능한 값이 나올 때를 막는다

무슨 일이 생기나

tri 라는 값이 0 이하가 되면 그다음 계산(brentq)이 실패합니다.

막지 않으면	막으면
프로그램이 멈춘다	적절한 기본값으로 넘어간다
13,095개 중 한 번만 나도 전체가 죽는다	계속 돈다

2만 번 도는 계산에서 한 번의 실패가 전체를 멈춘다

언제 그런 값이 나오나

자료가 극단적일 때입니다 — 분산이 거의 0 이거나, 값이 다 같거나요.

13,095개나 되니 그런 유전자가 있을 수 있습니다. 9강에서 걸렀지만 완벽하지 않아요.

비유 — 자판기의 예외 처리

자판기가 정상 동전만 처리하면? 구겨진 지폐 하나에 멈춰 버립니다.

잘 만든 자판기는 "이 지폐는 못 받습니다" 하고 계속 돌아갑니다.

전체를 멈추는 대신 그 하나만 다르게 처리한다

이게 실무 코드와 예제 코드의 차이입니다

교과서 예제는 깨끗한 자료를 가정합니다. 실무 코드는 이상한 것이 반드시 온다고 가정해요.

"이런 일은 안 일어날 것" 은 2만 번 돌면 일어난다

수축 덕분에 '우연히 조용한' 유전자가 더 이상 1등이 아닙니다.

Python 브라우저에서 실행

import numpy as np
import bio

L = bio.logcpm()
cols = list(range(7))
lfc, t, p = bio.moderated_ttest(L.values, cols[:4], cols[4:])

print(f"{'유전자':10s}{'log2FC':>9s}{'t':>9s}{'p':>12s}")
for g in ["Ddx54", "Cd47", "Cd38", "Myc", "Rassf3"]:
    i = L.index.get_loc(g)
    print(f"{g:10s}{lfc[i]:+9.2f}{t[i]:+9.2f}{p[i]:12.2e}")

print()
print(f"p < 0.05 인 유전자: {int((p < 0.05).sum()):,} / {len(p):,}")
print()
print("11강의 Rassf3 를 보세요. t 가 −66.3 → −8.56 으로 제자리를 찾았습니다.")

🤔 확인 문제

p = 0.001 을 바르게 읽은 것은?

🤔 확인 문제

2 * stats.t.sf(np.abs(t), dft) 에서 2 를 곱하는 이유는?

🤔 확인 문제

moderated t 에서 자유도가 d0 + d 로 늘어나는 것이 뜻하는 바는?