Education

1. 交絡因子による擬似相関について

(私の実行環境は、Windows 10 Pro for Workstations、Rのバージョンは3.4.4です。)

●交絡因子って

交絡因子の定義は資料によりやや異なるみたいなのですが、例えば一般社団法人日本疫学会のHPでは、『交絡因子は、2つの集団のアウトカムを比較する際に、1)アウトカムに影響を与える、2)要因と関連がある、3)要因とアウトカムの中間因子^※でない、という3つの条件を満たす』とされています。

※中間因子(中間変数)は、XとYの因果の流れの中間に位置付けられる変数のことを指します。例えば、X(長距離輸送) → M(馬体重減) → Y(着順下がる)におけるMが中間変数です。

そのような交絡因子Zが存在している場合、要因XとアウトカムYの間に因果関係がない場合でも、XとYの間には関連が生じてしまうので、安易に「XとYの間に因果関係あり」と結論付けてしまわないように注意が必要です。

以下、交絡因子による擬似相関をRを使って簡単に説明しようと思います。

Rで乱数を発生させて、以下の4つの条件(右図A)を満たすような架空の野外調査データを作成します。

化学物質Aの生態影響評価を目的として、全50地点で調査を実施 (めっちゃキツいですね)
化学物質A(その濃度をX₁とします)は、底生動物の個体数(Yとします)には影響を与えない (当然調査者はそのことを知りません)
一方、化学物質B(その濃度をX₂とします)は底生動物の個体数に負の影響を与える
化学物質AとBはともにある業種の工場に由来の化学物質であり、Aの濃度が高いところ(周辺にその業種の工場の数(Zとします)が多い地点)ではBの濃度も高い (つまり正の相関を示します)

このようなシチュエーションにおいて、Zは、1)アウトカム(底生生物の個体数)に対して(化学物質Bを通じて)影響を与え、2)着目する要因(化学物質Aの濃度)と関連があり、3)アウトカムと要因の中間変数ではないので、交絡因子となります。

#以上のような条件をRでの乱数生成で表現すると...

n.sample <- 50　#調査地点数

set.seed(307); seed <- sample(1:100000, 4, replace = F)　#乱数の種

#まずはZ (工場の数)を生成 (元論文では「工業地帯の面積」にしていたので、正規分布になっています...)

set.seed(seed[1]); factory <- rnorm(n.sample, mean = 100, sd = 15)

# X₁, X₂, Y を生成するときの誤差項を生成

set.seed(seed[2]); e1 <- rnorm(n.sample, mean = 0, sd = 15)

set.seed(seed[3]); e2 <- rnorm(n.sample, mean = 0, sd = 15)

set.seed(seed[4]); e3 <- rnorm(n.sample, mean = 0, sd = 15)

#X₁, X₂, and Y を生成

#物質A(X₁)の係数は0(影響なし)です。一方、物質B(X₂)の係数は-2(負の影響あり)としています。

X1 <- 1.0 * factory + e1; X2<- 0.8 * factory + e2

Abundance <- ceiling(0 * X1 + (-2.0) * X2 + 250 + e3)

そして、以上のような架空データセットに用いて、Y を目的変数、X₁ を説明変数とする単回帰モデルを構築すると...

mod <- glm(Abundance ~ X1, family = "gaussian")

summary(mod)

-------------------------------------------------------------------------------------------------------------------------------

Call:

glm(formula = Sp.richness ~ Substance.X1, family = "gaussian")

Deviance Residuals:

Min 1Q Median 3Q Max

-68.239 -23.656 -1.101 21.079 85.541

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 145.6117 27.3958 5.315 2.73e-06 ***

X1 -0.5343 0.2640 -2.024 0.0485 *

(---以下略---)

-------------------------------------------------------------------------------------------------------------------------------

となり、X₁ の係数が負の値、しかもp値が0.05より小さくなります。(p値はサンプルサイズに依存するのでここを強調するのはよくないですが...)

この結果を基に「化学物質Aは底生動物に負の影響を与えている！濃度を下げろ！」という結論を導くのは誤りなわけです。

ではこのような擬似相関に惑わされないようにするにはどうすればいいのかという話ですが、そもそもYに対するX₁の単回帰モデルにおいて、

X₁とYの間に有意な関連が生じてしまうのは、

X₁の値が変化するときには、Zの値が変化している

⇒Zの値が変化すると、X₂の値が変化する

⇒X₂の値が変化すると、Yの値が変化する

ということが裏で起こっているためです。

なので、Zの値を変化させなければ良いのです。

つまり、Zをある値で固定した時のX₁とYの関連の強さを評価すれば、乱数生成時に指定したように、(ほぼ)0になるはずです。回帰分析において、△△をある値で固定した時の〇〇と××の関連の強さというのは、「説明変数を○○と△△とした重回帰モデル(目的変数は××)における〇〇の偏回帰係数」をみることに他ならないので、実際に上記のデータセットを用いてそのような重回帰モデルを構築すると...

ex1 <- glm(Abundance ~ X1 + factory, family = "gaussian")

summary(ex1)

-------------------------------------------------------------------------------------------------------------------------------

Call:

glm(formula = Abundance ~ X1 + factory, family = "gaussian")

Deviance Residuals:

Min 1Q Median 3Q Max

-65.378 -18.940 -3.777 22.382 70.762

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 220.06130 33.18697 6.631 2.97e-08 ***

X1 0.08327 0.30101 0.277 0.78326

factory -1.31891 0.39008 -3.381 0.00146 **

(---以下略---)

-------------------------------------------------------------------------------------------------------------------------------

となり、X₁ の係数はほぼ0、p値も0.78と大きな値を示します。

なお、データ生成メカニズムが図Aのような因果構造である場合は、ZによるYへの影響は全てX₂を介するので、以下のようにX₂の値を固定しても、Zの値を固定したことと同じ意味になります。

ex2 <- glm(Abundance ~ X1 + X2, family = "gaussian")

summary(ex2)

-------------------------------------------------------------------------------------------------------------------------------

Call:

glm(formula = Abundance ~ X1 + X2, family = "gaussian")

Deviance Residuals:

Min 1Q Median 3Q Max

-25.843 -7.876 -1.043 10.322 38.577

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 239.77255 11.48447 20.88 <2e-16 ***

X1 0.06736 0.10357 0.65 0.519

X2 -1.89956 0.10899 -17.43 <2e-16 ***

(---以下略---)

-------------------------------------------------------------------------------------------------------------------------------

X₁ の係数はここでもほぼ0、p値も0.519です。

また、X₂ の係数もおおよそ2に近い値になっています。

●擬似相関が生じるメカニズムを数式を使って解説

(作成中)