Partial cross mapping eliminates indirect causal influences
교차사상은 비분리 비선형계에서 인과 방향을 읽어내지만 간접 인과를 직접 인과로 오독한다. PCM이 편상관으로 매개 경로의 몫을 덜어내는 방식과, 그 유비가 비선형 환경에서 어디까지 성립하는지를 검토합니다.
Paper: Leng, S., Ma, H., Kurths, J., Lai, Y.-C., Lin, W., Aihara, K., & Chen, L. (2020). Partial cross mapping eliminates indirect causal influences. Nature Communications, 11, 2632. https://doi.org/10.1038/s41467-020-16238-0
Abstract: 시계열에서 인과를 읽는 방법은 크게 세 계보로 나뉜다. Granger 인과, 엔트로피 기반 측도, 그리고 교차사상이다. 앞의 두 계보는 한 변수의 정보를 다른 변수에서 떼어낼 수 있다는 분리성을 전제하는데, 비선형 동역학계에서는 이 전제가 대체로 깨진다. Takens–Mañé 임베딩 정리에 따르면 계 전체의 정보가 단일 변수의 관측만으로 재구성될 수 있기 때문이다. 교차사상 계열은 바로 이 성질을 뒤집어 쓴다. Y의 재구성 다양체에서 X의 상태를 복원할 수 있다면 X의 정보가 Y에 들어가 있는 것이고, 따라서 X가 Y의 원인이라고 판정한다. 이 논문이 지적하는 문제는 그 판정이 직접과 간접을 구별하지 못한다는 것이다. X \to Z \to Y인 사슬에서 X의 정보는 Z를 거쳐 Y에 도달하므로 X와 Y 사이에도 강한 교차사상이 성립하고, 기존 방법은 이를 직접 링크로 읽는다. 저자들은 세 가지 도구, 즉 위상공간 재구성과 mutual cross mapping(MCM), 그리고 편상관을 결합한 partial cross mapping(PCM)을 제안한다. 전체 정보를 담은 사상 \hat{X}^{Y}와 매개변수를 거친 간접 정보만 담은 사상 \hat{X}^{\hat{Z}^{Y}}를 각각 만든 뒤, 후자를 조건으로 놓은 편상관 \varrho_D = |\mathrm{Pcc}(X, \hat{X}^{Y} \mid \hat{X}^{\hat{Z}^{Y}})|을 직접 인과의 지표로 삼는 것이다. 세 종 로지스틱 벤치마크에서 사슬 구조의 X \to Y에 대해 MCM 지표는 0.8681이지만 PCM 지표는 0.1871로 떨어지고, 여덟 종 네트워크와 DREAM4 유전자 조절망, 플랑크톤 먹이망, 홍콩 대기오염 자료에도 적용된다. 다만 이 방법이 덜어내는 것은 관측된 매개변수를 거친 몫뿐이며, 편상관이라는 선형 조작을 비선형 재구성 위에 얹었을 때 무엇이 남는지는 논문이 완전히 답하지 않는다.
핵심 요약
| 항목 | 설명 |
|---|---|
| 연구 질문 | 변수들이 비분리적이고 결합이 약하거나 중간 정도인 비선형 동역학계에서, 시계열만으로 직접 인과 링크를 간접 인과 링크와 구별할 수 있는가? |
| 진단 | 교차사상 계열은 비분리성 문제를 해결하지만 인과 이행성 때문에 간접 링크를 직접 링크로 오검출한다. 조건부 Granger 인과나 편전달 엔트로피는 직접·간접을 가르지만 분리성을 요구하므로 비선형계에 적용되지 않는다. |
| 해법 | 위상공간 재구성 + MCM + 편상관의 결합. 전체 정보 사상 \hat{X}^{Y}에서 매개 경로를 거친 사상 \hat{X}^{\hat{Z}^{Y}}의 몫을 편상관으로 덜어낸 \varrho_D를 직접 인과 지표로 쓴다. |
| 판정 규칙 | 문턱값 1 > T \gg 0에 대해 \varrho_C \ge \varrho_D \ge T면 직접 링크, \varrho_C \ge T \gg \varrho_D면 간접 링크만 존재, T > \varrho_C \ge \varrho_D면 링크 없음. 보조 지표 \gamma = \varrho_D / \varrho_C를 함께 본다. |
| 대표 결과 | 저자 보고 기준 세 종 사슬에서 \varrho_C^{X\to Y} = 0.8681(오검출) 대 \varrho_D^{X\to Y} = 0.1871(정확), 루프에서 0.8052 대 0.4467, DREAM4 다섯 망의 AUROC 0.7907·0.6543·0.7353·0.6639·0.8629, 플랑크톤 먹이망의 직접 링크 복원, 홍콩 대기오염에서 NO₂와 호흡성 부유입자만 심혈관질환의 원인으로 검출. |
| 명시된 범위 밖 | 동기화된 강결합 변수, 자율계가 아닌 계 전반, 큰 잡음, 그리고 조건화 집합에 들어가지 못한 미관측 공통원인. |
목차
- 비분리성이라는 출발점
- 교차사상은 어떻게 방향을 읽는가
- 인과 이행성: 간접 링크는 왜 실재하는가
- PCM: 매개 경로의 몫을 덜어내기
- 판정 규칙과 보조 지표
- 고차 확장과 큰 네트워크에서의 우회
- 세 종 벤치마크
- 여덟 종 네트워크와 문턱값 선택
- 기존 방법과의 비교
- 실데이터 세 종
- 편상관 유비는 비선형 환경에서 정확한가
- 문턱값과 유의성 처리에 대한 검토
- 이 방법이 하지 못하는 것
- 결론
1. 비분리성이라는 출발점
논문은 서론에서 인과 검출 문헌을 네 갈래로 정리한다. Granger 인과와 그 확장, 엔트로피 기반 측도, 동역학적 베이즈 추론, 그리고 mutual cross mapping이다. 이 중 앞의 두 갈래는 예측 프레임워크 위에 서 있다. 어떤 변수의 과거를 예측 모형에서 빼거나 넣었을 때 예측 오차가 달라지는지를 보는 방식이다.
이 방식이 성립하려면 한 변수의 정보를 다른 변수에서 분리할 수 있어야 한다. 논문은 이 조건이 비선형 동역학계에서 일반적으로 깨진다고 보고, 방법 절에서 그 이유를 정리한다. 연속시간 계 \dot{x} = F(x)의 상태가 컴팩트 다양체 위의 흡인자 A 안에서 움직이고 그 상자 세기 차원이 d_A라 하자. Takens–Mañé 임베딩 정리와 그 프랙탈 일반화에 따르면, 양의 지연 \tau와 매끄러운 관측함수 h에 대해 지연좌표 사상
\Gamma_{h,\varphi,\tau}(x) = [\,h(x),\ h(\varphi_\tau(x)),\ h(\varphi_{2\tau}(x)),\ \dots\,]^{\top}은 임베딩 차원이 2 d_A보다 크면 확률 1로 임베딩이 된다. 관측함수를 좌표함수 h(x) = x_i로 잡으면 원 다양체 M_x가 그림자 다양체 M_y로 사상되고, 임베딩 사상이 일대일이므로 두 다양체 위의 동역학은 위상적으로 공액이다.
y(t+\tau) = \psi_\tau(y(t)) = \Gamma \circ \varphi_\tau \circ \Gamma^{-1}(y(t))여기서 논문이 끌어내는 결론이 비분리성이다. 원래 계에서 성분 x_j의 미래는 모든 성분 x_1, \dots, x_n의 이력에 의존한다. 그런데 공액 관계 때문에 같은 미래가 단 하나의 변수 x_i의 이력과 임베딩 사상만으로도 결정된다. 즉 한 변수의 관측만으로 이루어진 예측이 모든 변수를 다 쓴 예측만큼 정확할 수 있다. 논문은 이것이 "Granger 인과의 발상을 무력화한다"고 적고, 어떤 변수의 정보를 다른 변수에서 제거하는 조작 자체가 일반적으로 불가능하다는 뜻으로 비분리성 개념을 도입한다.
교차사상 계열은 이 성질을 약점이 아니라 도구로 쓴다. 정보가 변수들 사이에 스며들어 있다면, 그 스며듦의 방향을 읽어 인과를 판정하면 된다는 것이다.
2. 교차사상은 어떻게 방향을 읽는가
세 변수 X = {x_t}{t=1}^{L}, Y = {y_t}{t=1}^{L}, Z = {z_t}_{t=1}^{L}를 지연좌표 임베딩해 세 개의 그림자 다양체 M_X, M_Y, M_Z를 만든다. 벡터는 다음과 같다.
\mathbf{x}_t = (x_t,\ x_{t-\tau_x},\ \dots,\ x_{t-(E_x-1)\tau_x})\mathbf{y}t, \mathbf{z}t도 같은 형태이며, E\xi는 임베딩 차원, \tau\xi는 시간 지연, r = \max_{\xi} {1 + (E_\xi - 1)\tau_\xi}가 시작 시점이다. 논문은 E를 false nearest neighbor(FNN), \tau를 delayed mutual information(DMI)으로 정하며 더 진전된 기법도 쓸 수 있다고 적는다.
교차사상의 핵심 정의는 이웃 집합의 대응이다. 임의의 \xi, \eta \in {x, y, z}에 대해
\hat{\mathcal{N}}^{\xi}(\eta_t) = \{\eta_{t'} \mid \xi_{t'} \in \mathcal{N}(\xi_t)\}로 둔다. 여기서 \mathcal{N}(\xi_t)는 M_\xi 위에서 \xi_t에 가장 가까운 점들의 집합이며, 개수는 보통 E_\xi + 1을 쓴다. E_\xi차원 공간에서 유계 심플렉스를 만드는 데 필요한 최소 점 개수라는 이유에서다. \xi = \eta이면 이 집합은 그냥 \mathcal{N}(\eta_t)가 되고, \xi \ne \eta이면 \mathcal{N}(\xi_t)에서 건너온 교차사상 이웃이 된다.
방향의 논리는 여기서 나온다. \mathcal{N}(\eta_t)에서 \hat{\mathcal{N}}^{\xi}(\eta_t)로의 의존성은, \eta_t를 만든 변수에서 \xi_t를 만든 변수로 가는 인과적 영향을 특징짓는다. 구체적으로는 M_Y 위에서 \mathbf{y}_t의 이웃들을 찾고, 그 시간 첨자를 그대로 M_X로 옮겨 대응점들의 가중평균을 취해 \hat{x}^{y}_t = \mathbb{E}[\hat{\mathcal{N}}^{y}(x_t)]를 만든다. 가중치는 이웃점과 \mathbf{y}_t 사이의 거리로 정한다. 그리고
\varrho_C = \left| \mathrm{Corr}(x_t,\ \hat{x}^{y}_t) \right|가 경험적 문턱값 T보다 크면 X에서 Y로 인과적 영향이 있다고 판정한다. 방향이 뒤집혀 보이는 이 규칙이 이 계열의 특징이다. Y의 다양체로 X를 복원할 수 있다는 것은 X의 정보가 Y의 궤적에 각인되어 있다는 뜻이고, 정보가 그 방향으로 흘렀다는 뜻이므로 X \to Y로 읽는다.
논문은 이 계열을 mutual cross mapping(MCM)이라 부르고 여섯 편의 선행 연구를 그 안에 묶는다. Sugihara 등의 2012년 Science 논문이 그중 하나이며, 위의 상관계수 지표를 그 논문에서 가져왔다고 명시한다. 이웃 개수 E+1의 근거는 Sugihara와 May의 1990년 비선형 예측 연구다.
절차는 방법 절에 Procedure A로 정리되어 있다. U에서 V로의 인과를 검출한다고 할 때, (1) 두 시계열을 지연좌표 임베딩하고, (2) 각 시점 t마다 \mathbf{v}_t의 이웃 \mathcal{N}(\mathbf{v}_t)를 E_v + 1개 찾고, (3) 같은 시간 첨자를 갖는 M_U 위의 점들을 찾아 거리 기반 가중평균으로 \hat{u}^{v}_t를 얻고, (4) |\mathrm{Corr}(u_t, \hat{u}^{v}_t)| 같은 지표로 \hat{U}^{V}와 U의 일치도를 재는 네 단계다.
여기에 시간 지연 탐색이 붙는다. 인과는 일정한 지연을 두고 나타나므로, 후보 지연 \tau_i들에 대해 평행이동한 Y와 X 사이의 상관을 최대화하는 지연을 찾아 쓴다.
절차에는 들어 있지 않은 것도 하나 짚어 둘 만하다. 이 논문의 판정은 고정된 시계열 길이에서 계산한 상관계수 하나를 문턱값과 비교하는 것이며, 라이브러리 길이를 늘려 가며 지표가 수렴하는지를 확인하는 단계는 Procedure A에 포함되어 있지 않다.
3. 인과 이행성: 간접 링크는 왜 실재하는가

원논문 Figure 1. (a) X와 Y 사이에만 방향성 상호작용이 있고 Z는 독립인 경우. (b) 세 변수가 단방향 사슬을 이루어 X에서 Y로 가는 간접 링크가 생기는 경우. (c) 인과 루프에서 이웃한 두 변수가 한 방향으로는 직접, 반대 방향으로는 간접 링크를 갖는 경우. (d) 변수가 많아지면 더 많은 간접 링크가 직접 링크로 오검출된다. Leng et al. (2020), Figure 1을 비평·설명을 위해 원형 그대로 인용했으며 원논문은 CC BY 4.0으로 공개되어 있다.
논문은 세 변수의 세 가지 구조로 문제를 세운다. Figure 1a처럼 두 변수만 한 방향으로 상호작용하고 세 번째가 고립돼 있으면 기존 방법으로 충분하다. Figure 1b처럼 세 변수가 단방향 사슬을 이루면, 이웃하지 않은 X와 Y 사이에 거짓 직접 링크가 검출된다. Figure 1c의 루프에서는 이웃한 모든 쌍이 한 방향으로 직접 링크를, 반대 방향으로 간접 링크를 갖게 되어 기존 방법이 간접을 직접으로 읽는다. 그리고 Figure 1d처럼 관측 변수가 늘어나면 이런 오검출의 가능성이 크게 늘어난다.
여기서 중요한 것은 간접 링크가 통계적 착시가 아니라는 점이다. 논문은 방법 절에서 이를 대수적으로 보인다. 세 종이 사슬로 연결된 로지스틱 모형
x_t = x_{t-1}(\alpha_x - \alpha_x x_{t-1}),\quad
z_t = z_{t-1}(\alpha_z - \alpha_z z_{t-1} - \beta_{zx} x_{t-1}),\quad
y_t = y_{t-1}(\alpha_y - \alpha_y y_{t-1} - \beta_{yz} z_{t-1})을 생각하자. 두 번째 식을 한 시간 단계 밀어 세 번째 식에 대입하면 y_t가 z_{t-2}를 통해 x_{t-2}에 의존하는 형태가 되고, 세 번째 식을 z에 대해 풀어 얻은
z_{t-2} = \frac{1}{\beta_{yz}}\left(\alpha_y - \alpha_y y_{t-2} - y_{t-1}/y_{t-2}\right)를 다시 대입하면 y_t가 y_{t-1}, y_{t-2}, 그리고 x_{t-2}만의 함수로 표현된다. 즉 X에서 Y로 가는 직접 결합항이 방정식에 없는데도, 궤적 수준에서는 X가 Y의 미래를 결정하는 데 실제로 참여한다. 논문은 이 의존이 이행성에서 유도된 간접적인 것이며, 이산시간계에서는 시간 지연 효과를 동반한다고 정리한다.
이 유도는 이 논문이 겨냥하는 표적이 무엇인지도 규정한다. 논의 절에서 저자들은 동역학계의 인과가 "X에 개입하면 Y가 달라지는가"라는 통계학의 전통적 관점과 다르며, 가장 단순하게는 계의 방정식에서 X로부터 Y로 가는 결합항으로 해석될 수 있다고 적는다. 그러니까 PCM이 찾으려는 직접 인과는 궤적에서 읽히는 의존이 아니라 방정식에 적힌 항이다. 모형을 모르는 상태에서 모형 수준의 대상을 겨냥한다는 이 구도가 이 논문의 야심이자 어려움의 근원이다.
4. PCM: 매개 경로의 몫을 덜어내기

원논문 Figure 2. (a) 세 변수가 단방향 사슬로 상호작용하는 설정에서 MCM은 \mathcal{N}(\mathbf{y}_t)를 M_X의 왼쪽 원형 영역 \hat{\mathcal{N}}^{y}(\mathbf{x}_t)로 사상하며, 이때 추정된 \hat{X}^{Y}가 참 X에 가까워 간접 링크가 직접 링크로 오판된다. (b) PCM은 \mathcal{N}(\mathbf{y}_t)에서 M_Z를 거쳐 M_X의 오른쪽 원형 영역으로 사상된 \hat{X}^{\hat{Z}^{Y}}를 조건으로 놓고 X와 \hat{X}^{Y}의 편상관을 계산한다. 기하적으로 \varrho_C는 전체 공간에서 X와 \hat{X}^{Y}가 이루는 각의 코사인이고, \varrho_D는 그것을 \hat{X}^{\hat{Z}^{Y}}에 직교하는 부분공간으로 사영한 값이다. Leng et al. (2020), Figure 2를 비평·설명을 위해 원형 그대로 인용했으며 원논문은 CC BY 4.0으로 공개되어 있다.
핵심 발상은 한 문장으로 요약된다. 한 시계열과 다른 시계열에서 만든 교차사상 예측의 일치도를 보되, 제3의 변수를 거쳐 전달된 몫을 조건으로 걸어 두는 것이다.
구현은 세 개의 사상을 만드는 일이다. 먼저 시계열 Y를 후보 지연 \tau_i(i = 1, \dots, m)만큼 평행이동해 Y^{\tau_i} = {y_{t+\tau_i}}를 만든다.
첫째, Z와 Y^{\tau_i} 쌍에 MCM을 적용해 \hat{Z}^{Y^{\tau_i}}를 얻고 상관을 최대화하는 지연을 고른다.
i_1 = \arg\max_{1 \le i \le m} \mathrm{Corr}(Z,\ \hat{Z}^{Y^{\tau_i}})이때의 사상을 간단히 \hat{Z}^{Y}로 쓴다. 이것은 Y의 다양체에서 복원한 Z다.
둘째, 이 \hat{Z}^{Y}를 다시 평행이동해 X와 짝지어 같은 절차를 반복한다.
i_2 = \arg\max_{1 \le i \le m} \mathrm{Corr}(X,\ \hat{X}^{\hat{Z}^{Y^{\tau_i}}})이렇게 얻은 \hat{X}^{\hat{Z}^{Y}}가 Z를 거친 간접 정보 흐름을 담는다. Y에서 Z로, 다시 Z에서 X로 두 번 건너간 사상이기 때문이다.
셋째, X와 평행이동한 Y^{\tau_i}에 직접 MCM을 적용해 \hat{X}^{Y}를 얻는다. 이것은 X에서 Y로 전달된 모든 인과 정보를 담는다.
i_3 = \arg\max_{1 \le i \le m} \mathrm{Corr}(X,\ \hat{X}^{Y^{\tau_i}})그리고 직접 인과의 지표를 다음과 같이 정의한다.
\varrho_D = \left| \mathrm{Pcc}\!\left(X,\ \hat{X}^{Y} \,\middle|\, \hat{X}^{\hat{Z}^{Y}}\right) \right|\mathrm{Pcc}(\cdot,\cdot\mid\cdot)는 세 번째 변수의 정보를 제거한 뒤 앞의 두 변수의 연관 정도를 재는 편상관계수이며, 논문은 그 정의를 그대로 적는다.
\mathrm{Pcc}(X, Y \mid Z_1) = \frac{\mathrm{Corr}(X,Y) - \mathrm{Corr}(X,Z_1)\,\mathrm{Corr}(Y,Z_1)}{\sqrt{\left(1 - \mathrm{Corr}(X,Z_1)^2\right)\left(1 - \mathrm{Corr}(Y,Z_1)^2\right)}}조건이 둘 이상이면 재귀적으로 정의된다.
\mathrm{Pcc}(X, Y \mid Z_1, Z_2) = \frac{\mathrm{Pcc}(X,Y\mid Z_1) - \mathrm{Pcc}(X,Z_2\mid Z_1)\,\mathrm{Pcc}(Y,Z_2\mid Z_1)}{\sqrt{\left(1 - \mathrm{Pcc}(X,Z_2\mid Z_1)^2\right)\left(1 - \mathrm{Pcc}(Y,Z_2\mid Z_1)^2\right)}}기하적 해석은 Figure 2b가 보여 준다. \varrho_C = |\mathrm{Corr}(X, \hat{X}^{Y})|는 전체 공간에서 X와 \hat{X}^{Y}가 이루는 각의 코사인이고, \varrho_D는 그 값을 간접 정보 \hat{X}^{\hat{Z}^{Y}}에 직교하는 부분공간으로 사영한 값이다. Figure 2a의 사슬 구조에서는 \hat{X}^{Y}가 담은 정보가 사실상 전부 Z를 거쳐 온 것이므로, 사영이 직각에 가까워지고 \varrho_D가 0 근처로 내려간다.
방법 절의 Procedure B는 이 과정을 여섯 단계로 정리한다. Y를 여러 지연으로 평행이동하고, Z에서 Y^{\tau_i}로 Procedure A를 돌려 \hat{Z}^{Y}를 고정하고, 그 \hat{Z}^{Y}를 다시 평행이동해 X와 짝지어 \hat{X}^{\hat{Z}^{Y}}를 고정하고, X와 Y^{\tau_i}로 \hat{X}^{Y}를 고정한 뒤, 마지막에 편상관을 계산한다.
5. 판정 규칙과 보조 지표
논문은 세 변수가 인과적으로 상호작용할 때 일반적으로 \varrho_C \ge \varrho_D가 성립한다고 적고, 문턱값 1 > T \gg 0을 두어 세 경우로 나눈다.
\varrho_C \ge \varrho_D \ge T:X에서Y로 직접 인과 링크가 있다.\varrho_C \ge T \gg \varrho_D: 직접 링크 없이 간접 링크만 있다.T > \varrho_C \ge \varrho_D: 어떤 인과 링크도 없다.
실제 자료에서는 전달되는 인과 신호가 충분히 강하지 않아 \varrho_C가 T를 겨우 넘고 \varrho_D도 T 근처에 놓이는 일이 생긴다. 이때 판정이 T 값에 민감해지므로, 논문은 두 지표의 근접도를 재는 보조 지표를 도입한다.
\gamma = \varrho_D / \varrho_C\gamma가 1에 가까울수록 직접 링크일 가능성이 높다는 해석이다. 방법 절에서는 고차 지표까지 포함한 일반형
\gamma = \left(\prod_{n=1}^{s} \varrho_{Dn}\right) \big/ \varrho_C^{\,s}를 제시하며, 본문의 정의는 s = 1인 경우에 해당한다. 보충자료에 따르면 여덟 종 모형에서 \gamma \ge 70%를 함께 요구하면 작은 T에서도 검출 정확도가 80% 이상으로 올라간다.
통계적 신뢰성은 다중검정 보정으로 확보한다. 보충자료는 편상관이 0이라는 귀무가설을 0이 아니라는 대립가설에 대해 검정해 p값을 얻고, 네트워크의 엣지 수 n에 대해 p < \alpha/n의 Bonferroni 조정을 적용했다고 밝힌다.
6. 고차 확장과 큰 네트워크에서의 우회
변수가 X, Y, Z_1, \dots, Z_s(s \ge 2)로 늘어나면, 매개 후보 전부를 조건에 넣은 1차 지표를 쓴다.
\varrho_{D1} = \left| \mathrm{Pcc}\!\left(X,\ \hat{X}^{Y} \,\middle|\, \left\{\hat{X}^{\hat{Z}_i Y}\ \middle|\ i = 1, \dots, s\right\}\right) \right|간접 경로가 두 변수를 거치는 경우(X \to Z_1 \to Z_2 \to Y)를 위해 2차 지표가 정의되고, 일반적으로 n개 매개변수의 모든 조합에 대해 n차 지표 \varrho_{Dn}이 정의된다. 다만 논문은 n이 커지면 조합 수가 매우 커진다는 이유로 고차 방법의 계산과 응용을 후속 연구로 미루고, 이 논문에서는 1차 문제만 다룬다고 명시한다.
네트워크 규모가 커져 조건화 집합이 커지면 편상관 계산 자체에 문제가 생긴다는 점도 인정한다. 논문이 제시하는 우회는 두 가지다. 하나는 \varrho_C^{X \to Z_i} + \varrho_C^{Z_i \to Y}(또는 두 값의 최솟값)를 최대화하는 노드 몇 개만 골라 조건화하는 것이다. Z_i를 거친 간접 링크가 존재할 확률이 높은 노드를 우선한다는 뜻이다. 다른 하나는 네트워크가 희소하다는 사전 지식이 있을 때 Z_1, \dots, Z_s를 하나씩 조건화하고 \varrho_D^{X \to Y \mid Z_i}의 최솟값을 최종 결과로 삼는 것이다.
방법 절 말미에는 프레임워크의 변형 가능성도 적혀 있다. 편상관 대신 조건부 의존성을 특징짓는 다른 측도를 쓸 수 있으며, 결정계수 r^2가 한 후보로 언급된다. 또 간접 영향 X \to Z \to Y를 없애려면 두 링크 중 하나만 끊어도 충분하다는 발상도 변형의 근거로 제시된다.
7. 세 종 벤치마크

원논문 Figure 3. (a) 세 가지 상호작용 구조와 각각의 결합 계수. (b) MCM이 검출한 X에서 Y로의 인과. 두 번째와 세 번째 구조에서 거짓 직접 인과가 나온다. (c) PCM이 검출한 직접 인과. (b)의 거짓 검출이 제거된다. 5000 길이 시계열에서 무작위로 뽑은 길이 1000의 시행 100회를 평균했고, 표본율은 1 Hz여서 시점 수가 계의 시간 단위와 정확히 대응한다. 재구성 파라미터는 E = 4, \tau = 1이다. Leng et al. (2020), Figure 3을 비평·설명을 위해 원형 그대로 인용했으며 원논문은 CC BY 4.0으로 공개되어 있다.
벤치마크는 세 종이 상호작용하는 로지스틱 사상이다.
x_t = x_{t-1}(\alpha_x - \alpha_x x_{t-1} - \beta_{xy} y_{t-1}) + \epsilon_{x,t}y_t = y_{t-1}(\alpha_y - \alpha_y y_{t-1} - \beta_{yx} x_{t-1} - \beta_{yz} z_{t-1}) + \epsilon_{y,t}z_t = z_{t-1}(\alpha_z - \alpha_z z_{t-1} - \beta_{zx} x_{t-1}) + \epsilon_{z,t}파라미터는 \alpha_x = 3.6, \alpha_y = 3.72, \alpha_z = 3.68이고 잡음은 평균 0, 표준편차 0.005의 백색잡음이다. 결합 계수를 바꿔 세 가지 구조를 만들며, Figure 3a에 표시된 값은 모두 0.35다. 첫 구조는 X \to Y 직접 링크만 있는 경우, 둘째는 X \to Z \to Y 사슬, 셋째는 X \to Z \to Y \to X 루프다.
문턱값 T = 0.5에서의 결과는 다음과 같다. 직접 링크만 있는 경우 MCM은 0.9847, PCM은 0.9817로 둘 다 정확하다. 사슬 구조에서는 MCM이 0.8681로 거짓 직접 링크를 보고하는 반면 PCM은 0.1871로 내려간다. 루프 구조에서는 MCM이 0.8052, PCM이 0.4467이다.
보충자료의 표에는 여섯 방향 전부의 값이 실려 있어 조금 더 볼 것이 있다. 사슬 구조에서 MCM이 틀리는 방향은 X \to Y만이 아니다. 참 사슬이 X \to Z \to Y인데 MCM은 역방향인 Z \to X에 0.6565, Y \to Z에 0.5307을 부여해 둘 다 문턱을 넘긴다. PCM은 같은 방향들에 0.3479와 0.3662를 주어 모두 걸러낸다. 즉 이 표에서 PCM이 얻는 이득은 이행성에서 생긴 간접 링크를 지우는 데 그치지 않고, 역방향 거짓 양성까지 함께 낮추는 데서 나온다. 논문의 서술은 간접 링크 제거에 초점이 맞춰져 있지만 표가 보여 주는 개선의 범위는 그보다 넓다.
루프 구조의 0.4467은 다르게 읽을 필요가 있다. 이 값은 T = 0.5 아래에 있으므로 정답이지만 여유가 0.05에 불과하다. 사슬에서의 0.1871과 비교하면 판정이 문턱값 선택에 훨씬 민감한 셈이고, 실제로 루프처럼 양방향 결합이 있는 구조가 이 방법에서 가장 어려운 경우임을 시사한다.
문턱값 강건성도 보충자료에 있다. 세 구조의 모든 링크에 대한 18개 지표로 계산한 검출 정확도를 T에 대해 그리면, PCM은 넓은 구간에서 100%에 가깝고 MCM은 0.85 < T < 0.95라는 좁은 구간에서만 100%에 도달한다. 논문은 실제로 0.85보다 큰 문턱값을 쓰는 것은 합리적이지 않다고 덧붙인다.
세 종의 다른 모든 상호작용 구조, 즉 링크가 0개인 경우부터 6개인 경우까지를 망라한 계통적 분석도 보충자료에 있다. 여기에는 공통원인에서 나오는 fan-out 허위 인과와 사슬에서 나오는 간접 허위 인과, 그리고 fan-in 구조가 포함되며, 논문은 이들이 복잡 네트워크의 기본 구성 단위이므로 이 모드들에서 정확하면 대규모 네트워크 재구성이 가능해진다고 논증한다.
8. 여덟 종 네트워크와 문턱값 선택
여덟 종이 결합된 로지스틱 모형에서 참 직접 링크는 일곱 개다(1 \to 3, 2 \to 3, 2 \to 4, 3 \to 5, 3 \to 6, 6 \to 7, 6 \to 8). 고유 성장률은 3.9, 3.5, 3.62, 3.75, 3.65, 3.72, 3.57, 3.68이고 결합 계수는 모두 0.35, 잡음은 표준편차 0.005다.
저자 보고에 따르면 T = 0.3에서도 PCM은 거의 모든 직접 링크를 복원하고, T = 0.5와 T = 0.7에서는 재구성 오차가 사실상 0이다. 같은 T 값에서 MCM은 직접·간접·거짓 링크가 뒤섞인 조밀한 네트워크를 내놓는다. 일곱 개 직접 링크의 p값은 다중검정 보정을 거쳐 모두 0에 가깝다고 보고된다.
문턱값을 어떻게 고를 것인가에 대해 논문은 두 가지 자동화 방안을 제시한다. 하나는 k-평균 군집화로 검출 결과를 직접 인과가 있는 군과 없는 군으로 나누는 것이다. 여덟 종 모형에서 이 분류기가 내놓은 문턱값은 MCM 결과에 대해 0.4234, PCM 결과에 대해 0.5684다. MCM 쪽에서는 참 직접 링크가 상대적으로 높은 값을 갖긴 해도 두 군이 뚜렷이 갈리지 않아 이 문턱값이 많은 위양성을 만든다. PCM 쪽에서는 두 군이 분리되고 각 군이 내부에 몰려 있어서 100% 정확도를 주는 문턱값의 범위가 꽤 넓다. 다른 하나는 t-검정으로 두 군의 차이를 확인하는 것인데, 이 네트워크의 모든 가능한 방향에서 p < 6.2397 \times 10^{-25}가 보고된다.
여기서 논문이 문턱값에 대해 내놓는 논증은 방향이 하나 더 있다. T를 크게 잡는 것이 항상 정확한 것은 아니라는 주장이다. 비선형 상호작용과 상관·편상관 계산, 그리고 내부·외부 섭동이 겹치면 지표의 절댓값이 어느 정도 줄어들 수밖에 없으므로, 실제 응용에서 큰 T를 고르는 것은 비현실적이라는 것이다. 이 서술은 뒤에서 볼 실데이터 결과와 함께 읽어야 한다.
강건성 시험도 보충자료에 있다. 잡음 규모를 0.002로 고정하고 시계열 길이를 100에서 2500까지 바꾸면 검출 정확도는 길이에 따라 올라가되 전 구간에서 높은 수준을 유지한다. 세로축 범위는 0.85에서 1이다. 길이를 1000으로 고정하고 잡음 규모를 0에서 0.016까지 올리면 정확도가 약간 내려간다. 다만 잡음 규모가 0.01을 넘으면 계 자체가 발산하므로 그 이상은 시험 대상이 아니다. 이 분석의 문턱값은 0.5이고, 음영은 100회 시뮬레이션의 표준편차다.
9. 기존 방법과의 비교
보충자료의 범위 표가 논문의 자리매김을 압축한다. 분리 가능한 계에서는 쌍별 분석에 GC, cGC, pGC, TE, pTE가, 네트워크 분석에 cGC, pGC, pTE가 쓰인다. 비분리 계에서는 쌍별 분석에 MCM과 PCM이 있고, 네트워크 분석에는 PCM만 있다. 논문이 채우려는 칸이 어디인지가 이 표에 그대로 드러난다.
정량 비교는 사슬 모드(\beta_{yz} = 0.1, \beta_{zx} = 0.1)에서 시계열 길이 100·200·500·1000·3000과 잡음 규모 0·0.005·0.01·0.015·0.02를 교차해 수행됐다. 길이 200, 잡음 0.005 조건에서는 PCM만 사슬 구조를 정확히 복원하고 나머지 다섯 방법은 모두 실패한다.
방법별 진단도 구체적이다. GC는 검출 결과가 상당히 불안정한데, 비선형계에 비분리 변수가 있기 때문이다. cGC와 pGC는 시계열이 충분히 길 때만 잘 작동하지만, 논문은 그 결과도 견고한 근거 위에 있지 않다고 본다. 회귀 기반 방법은 첫 단계에서 가정한 회귀 모형을 관측 시계열에 맞추는데, 이 적합이 참 동역학에서 크게 벗어나기 때문이라는 것이다. TE와 pTE는 관련 변수들의 확률분포를 자료에서 재구성해야 하므로 성능이 시계열 길이에 크게 의존하고, 이론적으로는 비분리 변수를 허용하지 않는다. 다만 논문은 강한 확률성을 가진 모형이 만든 자료에 대해서는 GC 계열이 여전히 가장 유용한 기법이라고 인정한다.
동역학적 베이즈 추론(DBI)과의 비교는 우열이 아니라 분업으로 정리된다. DBI는 명시적 방정식을 요구하지 않되 함수 공간의 기저를 세심하게 골라야 하고, 자율계와 비자율계 모두에 적용된다. PCM은 임베딩 정리에만 기대는 모형 자유 방법이라 자율계와 특정 형태의 비자율계에 적용된다. DBI는 결합 함수와 동역학 기제 자체를 추론할 수 있고 PCM은 인과 관계의 검출에 집중한다. 논문은 이 차이를 DBI가 찾는 연결이 실효 연결성(effective connectivity)에, PCM이 찾는 관계가 기능 연결성(functional connectivity)에 가깝다고 표현한다. 그리고 매우 복잡한 네트워크에서 PCM으로 기본 구조를 먼저 찾아 DBI의 초기 회귀 구조를 단순화하는 식의 상보적 사용을 제안한다.
보충자료는 PCM 자신의 약점도 같은 절에서 적는다. 분리성 조건에 매이지 않고 시계열이 충분히 길지 않아도 효과적이지만, 큰 잡음에는 민감하며 위음성 검출이 나오는 것은 불가피하다는 것이다.
10. 실데이터 세 종

원논문 Figure 4. (a) GeneNetWeaver에서 얻은 20개 상호작용 유전자로 이루어진 다섯 유전자 조절 네트워크 중 하나. 빨간 화살표는 활성화, 파란 화살표는 억제다. (b) PCM 검출 성능의 ROC 곡선과 AUROC 값(E = 2, \tau = 1). (c) 세 플랑크톤 종의 먹이사슬. 빨간 화살표는 피식자에서 포식자로 가는 상호작용이다. (d) 직접 인과 링크의 검출을 뜻하는 PCM 지표들(빨간 상자)과, 윤충류에서 피코시아노박테리아로 가는 상대적으로 약한 직접 링크(노란 상자). E = 4, \tau = 1. (e) 대기오염 물질과 심혈관질환 사이에서 검출된 상호작용(E = 7, \tau = 1). (f) (e)의 결과로 재구성한 인과 네트워크. 모든 검출 결과는 다중검정 보정으로 검증됐다. Leng et al. (2020), Figure 4를 비평·설명을 위해 원형 그대로 인용했으며 원논문은 CC BY 4.0으로 공개되어 있다.
10.1 DREAM4 유전자 조절 네트워크
DREAM4 in silico 네트워크 챌린지의 다섯 네트워크를 쓴다. 각 네트워크는 100개 유전자를 갖고, GeneNetWeaver로 상호작용하는 20개 유전자를 무작위로 골랐다. 각 유전자마다 21개 시점 발현 시계열의 실현 10벌이 있고, 위상공간 재구성을 위해 모든 실현을 하나의 시계열로 이어 붙였다. 재구성 파라미터는 E = 2, \tau = 1이다.
PCM이 검출한 직접 링크를 다섯 네트워크의 기지 엣지와 비교해 그린 ROC 곡선의 AUROC는 각각 0.7907, 0.6543, 0.7353, 0.6639, 0.8629다. 본문은 다섯 값의 평균이 0.75 정도에 이른다고 적는데, 표시된 값들의 산술평균은 0.7414다. 논문은 이 정도면 작은 자료에서도 직접 링크 검출 정확도가 높은 것이며 이 과제에서 MCM을 능가한다고 평가한다.
10.2 플랑크톤 먹이망
발트해에서 분리한 플랑크톤 군집의 8년 메소코즘 실험 자료에서 피코시아노박테리아, 윤충류, 요각류 세 종을 쓴다. 원 자료는 1990년 7월 12일부터 1997년 10월 20일까지 불균등 간격으로 측정됐고, 등간격으로 보간해 794개 시점을 얻었다. 한 시점의 시간 단위는 3.35일이다.
PCM은 여섯 개 가능한 방향 전부에 지표를 매기고, 값이 10^{-1} 이상인 링크만 남긴다. 그 결과 두 개의 직접 인과 링크가 남았고 원 네트워크의 정답과 일치한다. 특히 피코시아노박테리아에서 요각류로 가는 간접 링크는 배제된다. 보충자료의 p값을 보면 이 판정의 구조가 분명하다. 요각류에서 윤충류로 0.7495, 윤충류에서 요각류로 1.1302 \times 10^{-17}, 요각류에서 피코시아노박테리아로 0.0054, 피코시아노박테리아에서 요각류로 0.8259, 윤충류에서 피코시아노박테리아로 0.0011, 피코시아노박테리아에서 윤충류로 1.2119 \times 10^{-5}다. 배제된 간접 링크, 즉 피코시아노박테리아에서 요각류로의 방향은 p값이 0.8259로 여섯 방향 중 가장 크다.
세 번째 링크가 흥미로운 대목이다. 윤충류에서 피코시아노박테리아로 가는 약한 직접 링크가 검출되는데, 이는 포식자에서 피식자로 가는 방향이다. 논문은 포식자가 사냥을 통해 피식자 개체군에 유의한 영향을 줄 수 있으며, 피식자가 아주 풍부하지 않을 때 특히 그렇다고 설명하고, 실제 포식 위계가 직접 인과 링크와 반드시 일치하지는 않는다고 정리한다.
10.3 홍콩 대기오염과 심혈관질환
홍콩의 대기 측정소에서 얻은 이산화질소, 이산화황, 호흡성 부유입자, 오존의 일별 농도와 주요 병원의 일별 심혈관질환 입원 수를 쓴다. 본문은 1994년부터 1997년까지의 자료라 적고, 보충자료는 1995년 초 병상이 갑자기 늘어난 효과를 피하려고 1995년 3월부터 1997년 11월까지 1000일을 골랐다고 밝힌다. 재구성 파라미터는 E = 7, \tau = 1이고 유의성 기준은 p < 10^{-9}다.
검출 결과는 열한 개 방향이다. NO₂와 호흡성 부유입자에서 심혈관질환으로, NO₂와 SO₂ 사이 양방향, NO₂와 호흡성 부유입자 사이 양방향, SO₂와 호흡성 부유입자 사이 양방향, 오존과 호흡성 부유입자 사이 양방향, 그리고 오존에서 SO₂로 가는 단방향이다. 즉 오염물질 중에서는 NO₂와 호흡성 부유입자만이 질환의 주요 원인으로 검출되고 SO₂와 오존은 원인으로 식별되지 않는데, 논문은 이것이 선행 연구와 일치한다고 적는다. 인정된 오염물질들 사이의 관계는 모두 양방향이고, 오존과 SO₂ 사이에서만 단방향이 나온다.
여기서 논문은 자기 결과의 해석에 단서를 단다. 검출된 인과 관계들이 직접일 수도 간접일 수도 있다는 것이다. 온도, 습도, 풍속 같은 다른 요인의 자료가 완전히 확보되지 않았고, 이들이 일부 오염물질의 공통원인일 수 있기 때문이다. 논문은 이를 보충자료의 fan-out 상호작용 모드로 명시적으로 연결한다.
10.4 보충 실데이터
보충자료에는 두 가지 예가 더 있다. 하나는 실험용 쥐 시교차상핵 배양세포의 마이크로어레이 자료로 재구성한 일주기 리듬 유전자 조절 네트워크다. 16개 시점을 136개 시점으로 보간해 썼고, 17개 핵심 일주기 유전자와 20개 키나아제 유전자로 이루어진 37개 유전자 네트워크에서 PCM의 AUC는 0.6506, GC 계열과 TE 계열은 0.44에서 0.54 사이다. 핵심 일주기 유전자와 PFK 계열로 이루어진 18개 유전자 네트워크에서는 PCM이 0.7259다. 논문은 단일하고 짧은 시계열 하나로 유전자 조절망을 추론하는 것이 어려운 과제이며, 기존 방법이 합성 자료에서 0.7 안팎, 실제 실험 자료에서 0.5 안팎에 머문다는 선행 보고를 인용한다.
다른 하나는 10명 피험자의 19채널 EEG 기록이다. 표본율 256 Hz, 기록 길이 약 1분이며, 4초(1024점) 블록으로 나눠 피험자당 10개 블록을 무작위로 골랐다. 재구성 파라미터는 E = 6, \tau = 6이다. PCM 결과는 신호가 주로 인접 영역을 통해 전달됨을 보여 주고, MCM은 거의 모든 방향으로 연결된 조밀한 그림을 내놓아 특정 경로를 읽을 수 없다는 것이 논문의 정리다.
11. 편상관 유비는 비선형 환경에서 정확한가
이 방법의 이름과 그림은 하나의 유비 위에 서 있다. 선형 회귀에서 제3변수의 몫을 덜어내듯, 재구성 상태공간에서 매개변수를 거친 몫을 덜어낸다는 것이다. 이 유비가 어디까지 정확한지가 이 논문을 읽는 데 가장 중요한 질문이다. 논문 자신의 서술을 근거로 네 지점을 짚을 수 있다.
첫째, 조작 자체는 선형이다. \varrho_D를 만드는 편상관 공식은 상관계수들의 유리식이며, 세 시계열 사이의 선형 연관만을 다룬다. 덜어내지는 것은 \hat{X}^{\hat{Z}^{Y}}가 X와 \hat{X}^{Y} 각각에 대해 갖는 선형 성분이다. 매개 경로의 기여가 비선형이라면 그 잔여는 남는다. 논문이 이 지점을 모르는 것은 아니다. 방법 절 말미에서 편상관을 조건부 의존성을 특징짓는 다른 측도로 바꿀 수 있다고 적고 결정계수를 후보로 든다. 즉 편상관은 이 프레임워크의 필연이 아니라 선택이며, 그 선택의 대가가 무엇인지는 이 논문에서 측정되지 않는다.
둘째, \varrho_C \ge \varrho_D는 정리로 제시되지 않는다. 논문은 "세 변수가 인과적으로 상호작용할 때 일반적으로 \varrho_C \ge \varrho_D를 갖는다"고 쓰고, 기하적 설명에는 "heuristically"와 "intuitively"라는 표현을 붙인다. Figure 2b의 사영 그림은 이 부등식이 성립할 때의 그림이고, 부등식 자체의 유도는 없다. 편상관은 조건화 뒤에 원래보다 커질 수 있는 양이므로, 이 순서는 일반적으로 보장되는 성질이 아니라 이 논문이 다루는 구조들에서 관찰된 규칙성에 가깝다. 판정 규칙의 세 경우가 모두 이 순서를 전제로 쓰여 있다는 점을 함께 보면, 순서가 깨지는 구조에서 규칙이 어떻게 작동하는지는 열려 있다.
셋째, 조건화 변수 자체가 추정물이다. \hat{X}^{\hat{Z}^{Y}}는 M_Y에서 M_Z로, 다시 M_Z에서 M_X로 두 번 건너간 결과다. 두 단계 중 어느 쪽이라도 약하면 이 대리물은 간접 정보를 온전히 담지 못하고, 그만큼 덜어내기가 불완전해져 위양성 쪽으로 기운다. 반대 방향의 문제도 있다. 이 논문이 출발점으로 삼은 비분리성 논변에 따르면, 결정론적 비선형계에서는 계 전체의 정보가 어느 한 변수에도 실려 있다. 그렇다면 M_Z에는 Z를 거치지 않은 X의 정보까지 일반적으로 담겨 있고, 그 위에서 만든 \hat{X}^{\hat{Z}^{Y}}를 조건으로 걸면 직접 성분까지 함께 덜어낼 수 있다. 논문은 이 실패를 보충자료에서 직접 보고한다. 상호작용이 매우 조밀한 비교적 큰 네트워크에서는 과다 조건화 때문에 위음성이 나올 수 있으며, 이를 피할 최적화 기법이 현재 연구 과제라는 것이다. 37개 유전자 네트워크의 AUC 0.6506은 그 상황에서 나온 값이다. 즉 이 방법이 잘 작동하는 영역은 간접 경로가 뚜렷하되 네트워크가 지나치게 조밀하지 않은 중간 지대이며, 그 경계를 사전에 판정할 기준은 제시되지 않았다.
넷째, 지연 탐색이 세 번 반복되는 최대화다. Procedure B의 세 단계는 각각 후보 지연 m개에 대해 상관을 최대화하는 지연을 고른다. 그렇게 고른 최댓값들이 그대로 편상관 계산에 들어가고, 그 편상관의 p값도 지연이 미리 정해져 있었던 것처럼 계산된다. 선택 과정에서 생기는 낙관적 편향은 보정되지 않는다. 게다가 세 사상은 서로 다른 최적 지연 \tau_{i_1}, \tau_{i_2}, \tau_{i_3}에서 정렬된 것이어서, 최종 편상관은 각기 다른 기준으로 맞춰진 세 계열을 한자리에 놓고 계산된다. 논문은 이 전략을 일관성을 위해 모든 MCM 결과에도 동일하게 적용했다고 밝히므로 비교 자체는 공정하지만, 지표의 절대 수준을 해석할 때는 감안해야 한다.
12. 문턱값과 유의성 처리에 대한 검토
이 방법의 판정은 결국 \varrho_D를 T와 비교하는 일이다. 그런데 자료마다 그 비교의 형태가 달라진다.
합성 벤치마크에서는 T = 0.5가 쓰인다. 플랑크톤 자료에서는 지표 값이 10^{-1} 이상인 링크를 남긴다. 대기오염 자료에서는 지표 문턱이 아니라 p값 컷오프 10^{-9}가 판정을 한다. 세 자료에서 판정 규칙의 형태가 세 번 바뀌는 셈이다.
지표의 절대 규모가 자료에 따라 크게 다르다는 사실이 이 변화의 배경이다. Figure 4d의 색 눈금은 0에서 0.3까지이고 Figure 4e는 0에서 0.5까지다. 벤치마크에서 참 직접 링크가 0.95 부근에 놓였던 것과 비교하면, 실데이터에서는 지표 전체가 아래로 눌린다. 논문이 여덟 종 모형 분석에서 큰 T를 고르는 것이 비현실적이라고 논증한 것도 같은 사정을 가리킨다. 다만 그 논증이 실제 응용에서 취하는 형태는 T를 다섯 배 낮추거나(0.5에서 0.1로) 아예 다른 기준으로 갈아타는 것이며, 자료마다 어느 쪽을 택할지에 대한 규칙은 제시되지 않는다.
유의성 검정이 하는 일도 좁게 보아야 한다. p값은 편상관이 0이라는 귀무가설에 대한 것이다. 길이 1000의 시계열을 100회 시행 평균하는 설정에서 이 귀무가설은 거의 언제나 기각된다. 보충자료의 표가 이를 잘 보여 준다. Figure 1a처럼 Z가 고립되고 X \to Y만 있는 구조에서, 존재하지 않는 링크인 Y \to X의 PCM 지표는 0.3322인데 p값은 6.9751 \times 10^{-26}이다. 사슬 구조에서 정확히 제거된 X \to Y의 지표 0.1871에도 p값 6.3838 \times 10^{-9}가 붙는다. 두 경우 모두 "편상관이 0이 아니다"는 참이지만 "직접 링크가 있다"는 거짓이다. 판정을 실제로 수행하는 것은 문턱값 T이고 p값은 그 옆에서 다른 질문에 답하고 있다.
컷오프 값 자체도 자료마다 다르다. 플랑크톤에서는 0.0016, 대기오염에서는 10^{-9}가 쓰인다. 보충자료에 적힌 다중검정 보정 규칙은 네트워크의 엣지 수 n에 대한 p < \alpha/n인데, 다섯 노드 계에서 이 규칙으로 10^{-9}가 나오지는 않는다. 보충자료의 서술도 이 값을 "설정했다"고 적는다. 대기오염 자료에서 이 선택은 실질적인 차이를 만든다. NO₂와 오존 사이의 두 방향은 p값이 2.8376 \times 10^{-9}와 6.9308 \times 10^{-9}로 컷오프 바로 바깥에 있어서 최종 네트워크에서 빠진다. 한 자릿수만 느슨했어도 두 링크가 추가됐을 것이다.
문턱값 선택을 자동화하려는 시도가 보충자료에 있다는 점은 함께 적어 둘 만하다. k-평균 군집화로 두 군을 나누는 방법과 t-검정으로 군 간 차이를 확인하는 방법이 제시되고, 여덟 종 모형에서 k-평균이 준 0.5684가 100% 정확도를 낸다. 다만 이 자동화가 적용된 것은 정답을 아는 합성 모형이고, 세 실데이터의 문턱값은 자동화된 절차가 아니라 개별 판단으로 정해졌다.
13. 이 방법이 하지 못하는 것
논문이 명시적으로 밝힌 범위 밖의 항목들부터 정리하면 다음과 같다.
미관측 공통원인. PCM이 덜어내는 것은 조건화 집합에 들어간 변수를 거친 몫이다. 조건화 집합은 관측된 변수들로만 구성되므로, 보이지 않는 공통원인이 만든 상관은 그대로 남는다. 논문은 이 한계를 대기오염 분석에서 스스로 적용한다. 온도·습도·풍속 자료가 없어서 검출된 관계가 직접인지 간접인지 확정할 수 없다는 것이다. 보충자료의 fan-out 모드 그림이 이 상황의 도식이며, PCM이 fan-out 허위 인과를 가려낼 수 있는 것은 공통원인이 관측되어 조건에 들어갈 때다.
동기화. 논문은 강하게 결합해 동기화된 변수들이 이 프레임워크의 범위 밖이라고 명시한다. 이 경우 계 전체가 원인 계의 부분다양체로 붕괴하고 결과 변수가 원인 계 위의 관측함수가 되어, 계산상 항상 양방향 인과가 검출되기 때문이다.
비자율계. PCM은 Takens–Mañé 정리에 기대므로 이론적으로 자율계에만 적용된다. 논문은 예외를 구체적으로 적는다. 전환점을 찾을 수 있고 전환 사이 구간이 충분히 긴 스위칭 계에서는 구간별 인과를 수치적으로 검출할 수 있으며, 일반화된 임베딩 정리에 기대어 일부 강제 계나 약한·중간 정도 잡음이 있는 계에도 적용된다. 시간에 따라 변하는 결합 함수를 가진 진동자처럼 중요한 비자율계에는 동역학적 베이즈 추론이 실용적 해법을 준다고 인정한다.
큰 잡음과 짧은 자료. 보충자료는 PCM이 큰 잡음에 민감하며 위음성이 불가피하다고 적는다. 다만 짧은 자료에 대해서는 상대적으로 강한 편이어서, 길이 200에서 다른 다섯 방법이 모두 실패하는 조건에서도 사슬 구조를 복원한다.
조밀한 대형 네트워크. 앞 절에서 본 과다 조건화 문제다. 논문은 이를 현재 연구 과제로 남긴다.
2차 이상의 간접 경로. 매개변수를 둘 이상 거치는 간접 경로를 위한 고차 지표가 정의돼 있지만, 조합 폭발을 이유로 계산과 응용은 후속 연구로 미뤄졌고 이 논문의 모든 결과는 1차 방법으로 얻어졌다.
결합의 부호와 크기. PCM은 링크의 존재를 판정하는 방법이며, 결합 함수나 그 부호를 추론하지 않는다. 논문 자신이 DBI와 비교하며 PCM이 찾는 것은 기능 연결성에 가깝다고 위치 짓는다.
여기에 논문의 설정에서 따라 나오는 항목을 두 가지 덧붙일 수 있다.
하나는 DREAM4에서의 자료 구성이다. 21개 시점의 실현 10벌을 하나의 시계열로 이어 붙여 위상공간을 재구성했는데, 이어붙인 경계를 걸치는 지연 벡터는 하나의 궤적 위의 상태가 아니다. E = 2, \tau = 1이므로 경계마다 그런 벡터가 하나씩 생기는 정도라 영향의 크기는 작겠지만, 이 자료에서의 AUROC 0.65–0.86이 벤치마크의 거의 완벽한 정확도와 크게 벌어진 원인을 논문이 분석하지는 않는다.
다른 하나는 재현성이다. 코드는 공개돼 있고, 사용한 공개 자료의 출처도 보충자료에 정리돼 있다. 다만 이 연구에서 생성·분석된 자료는 합리적 요청 시 교신저자로부터 제공된다는 문구로 처리돼 있다. 본문에 인용되는 대표 수치들은 100회 시행의 평균이며, 그림의 음영으로 표시된 표준편차 외에는 분산 정보가 함께 제시되지 않는다.
14. 결론
이 논문의 기여는 두 문장으로 정리된다. 교차사상 계열이 비분리 비선형계에서 인과 방향을 읽을 수 있게 해 준 대신 직접과 간접을 구별하지 못했는데, 그 구별을 조건화로 회복하는 절차를 만들었다는 것이다. 절차 자체는 단순하다. 전체 정보를 담은 사상과 매개 경로만 담은 사상을 각각 만들고, 후자를 조건으로 놓은 편상관을 지표로 쓴다.
이 단순함에는 장점과 대가가 함께 있다. 장점은 기존 교차사상 구현 위에 얹기 쉽고, 추가로 요구하는 것이 하나의 관측 변수뿐이라는 점이다. 대가는 두 겹이다. 하나는 조건화가 관측된 매개변수에만 미친다는 것이고, 다른 하나는 비선형 상태공간 위에서 선형 편상관을 쓴다는 것이다. 논문은 두 대가를 모두 인지하고 있다. 전자는 대기오염 분석의 단서로, 후자는 편상관을 다른 조건부 의존 측도로 바꿀 수 있다는 서술로 나타난다.
수치가 보여 주는 것도 이 구도와 맞는다. 결합 구조가 뚜렷하고 잡음이 작은 합성계에서는 개선 폭이 크다. 사슬 구조의 X \to Y에서 0.8681이 0.1871로 내려가고, 여덟 종 네트워크는 T = 0.5와 0.7에서 사실상 오차 없이 복원된다. 반면 실데이터로 갈수록 지표의 절대 수준이 내려가고 판정 규칙이 자료마다 달라진다. DREAM4에서 AUROC는 0.65에서 0.86 사이이고, 37개 유전자의 실제 조절망에서는 0.65로 떨어지며, 이 실패를 논문 스스로 과다 조건화로 진단한다.
그래서 이 방법이 검증한 명제는 좁게 읽는 편이 정확하다. 관측 변수 집합 안에서, 결합이 약하거나 중간 정도이고 잡음이 작으며 네트워크가 지나치게 조밀하지 않을 때, 매개 경로를 거친 몫을 덜어내면 직접 링크가 간접 링크와 갈린다는 것이다. 이 조건들 중 어느 하나가 무너질 때 지표가 어떻게 행동하는지, 그리고 그 무너짐을 자료만 보고 알아챌 수 있는지는 논문이 답하지 않은 부분으로 남는다. 문턱값 선택의 자동화와 고차 방법의 계산, 과다 조건화 회피는 저자들 스스로 후속 과제로 지목한 항목들이다.
References
Baba, K., Shibata, R., & Sibuya, M. (2004). Partial correlation and conditional correlation as measures of conditional independence. Australian & New Zealand Journal of Statistics, 46(4), 657–664. https://doi.org/10.1111/j.1467-842X.2004.00360.x
Barnett, L., & Seth, A. K. (2014). The MVGC multivariate Granger causality toolbox: A new approach to Granger-causal inference. Journal of Neuroscience Methods, 223, 50–68. https://doi.org/10.1016/j.jneumeth.2013.10.018
Benincà, E., Huisman, J., Heerkloss, R., Jöhnk, K. D., Branco, P., Van Nes, E. H., Scheffer, M., & Ellner, S. P. (2008). Chaos in a long-term experiment with a plankton community. Nature, 451, 822–825. https://doi.org/10.1038/nature06512
Benincà, E., Jöhnk, K. D., Heerkloss, R., & Huisman, J. (2009). Coupled predator–prey oscillations in a chaotic food web. Ecology Letters, 12(12), 1367–1378. https://doi.org/10.1111/j.1461-0248.2009.01391.x
Clemson, P. T., & Stefanovska, A. (2014). Discerning non-autonomous dynamics. Physics Reports, 542(4), 297–368. https://doi.org/10.1016/j.physrep.2014.04.001
Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica, 37(3), 424–438. https://doi.org/10.2307/1912791
Leng, S., Ma, H., Kurths, J., Lai, Y.-C., Lin, W., Aihara, K., & Chen, L. (2020). Partial cross mapping eliminates indirect causal influences. Nature Communications, 11, 2632. https://doi.org/10.1038/s41467-020-16238-0
Ma, H., Aihara, K., & Chen, L. (2014). Detecting causality from nonlinear dynamics with short-term time series. Scientific Reports, 4, 7464. https://doi.org/10.1038/srep07464
Ma, H., Leng, S., Tao, C., Ying, X., Kurths, J., Lai, Y.-C., & Lin, W. (2017). Detection of time delays and directional interactions based on time series from complex dynamical systems. Physical Review E, 96(1), 012221. https://doi.org/10.1103/PhysRevE.96.012221
Mañé, R. (1981). On the dimension of the compact invariant sets of certain non-linear maps. In D. A. Rand & L.-S. Young (Eds.), Dynamical systems and turbulence, Warwick 1980 (Lecture Notes in Mathematics, Vol. 898, pp. 230–242). Springer. https://doi.org/10.1007/BFb0091916
Marbach, D., Prill, R. J., Schaffter, T., Mattiussi, C., Floreano, D., & Stolovitzky, G. (2010). Revealing strengths and weaknesses of methods for gene network inference. Proceedings of the National Academy of Sciences, 107(14), 6286–6291. https://doi.org/10.1073/pnas.0913357107
Schaffter, T., Marbach, D., & Floreano, D. (2011). GeneNetWeaver: In silico benchmark generation and performance profiling of network inference methods. Bioinformatics, 27(16), 2263–2270. https://doi.org/10.1093/bioinformatics/btr373
Schreiber, T. (2000). Measuring information transfer. Physical Review Letters, 85(2), 461–464. https://doi.org/10.1103/PhysRevLett.85.461
Stankovski, T., Pereira, T., McClintock, P. V. E., & Stefanovska, A. (2017). Coupling functions: Universal insights into dynamical interaction mechanisms. Reviews of Modern Physics, 89(4), 045001. https://doi.org/10.1103/RevModPhys.89.045001
Stark, J. (1999). Delay embeddings for forced systems. I. Deterministic forcing. Journal of Nonlinear Science, 9(3), 255–332. https://doi.org/10.1007/s003329900072
Sugihara, G., & May, R. M. (1990). Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series. Nature, 344, 734–741. https://doi.org/10.1038/344734a0
Sugihara, G., May, R., Ye, H., Hsieh, C.-H., Deyle, E., Fogarty, M., & Munch, S. (2012). Detecting causality in complex ecosystems. Science, 338(6106), 496–500. https://doi.org/10.1126/science.1227079
Takens, F. (1981). Detecting strange attractors in turbulence. In D. A. Rand & L.-S. Young (Eds.), Dynamical systems and turbulence, Warwick 1980 (Lecture Notes in Mathematics, Vol. 898, pp. 366–381). Springer. https://doi.org/10.1007/BFb0091924