はじめに：論文紹介「心理学的研究における重回帰分析の適用に関わる諸問題」

こんにちは S です。
重回帰分析（以降、重回帰=線形重回帰とする）はある入力とそれに対応する出力の関係を、入力値の線形和で表す予測手法であり、予測の仕組みのわかりやすさや、「データから何かを予測、説明する」という汎用的な使い勝手から広く使われています。
線形重回帰モデルを数式で表すと以下のようになります。
$$
y_i = \beta_0 + x_{i1}\beta_1 + \ldots + x_{ip}\beta_p + \epsilon_i, \quad i = 1, \ldots, n.
$$
$y_i$：従属変数
$x_{ij},,(j = 1, \ldots, p)$：独立変数
$\beta_0, \beta_j$：偏回帰係数（推定対象となるパラメータ）
$\epsilon_i$：互いに独立に平均0、分散 $\sigma^2$ の正規分布に従う誤差項（あくまで仮定）
$n$：サンプル数
$p$：特徴量数（次元）

これを見ればわかる通り、線形重回帰モデルは「従属変数は独立変数の線形和に誤差が加わったもの」という非常にわかりやすいモデルで成り立っています。しかし、どんな独立変数をどう組み込むべきか、多重共線性、$x_1 \times x_2$ といった交互作用項を入れるべきかなど、実は考えるべきことが多いとても奥の深い手法であり、故に誤用・誤解が散見されます。これに対し、研究のプロでさえ誤用してしまうという現状に警鐘を鳴らすべく立ち上がった方々がいました、それが以下の論文の著者です。

心理学的研究における重回帰分析の適用に関わる諸問題

2021年5月に発表されたこの論文は、心理学分野の論文における重回帰分析の誤用を名指しでひたすら指摘しまくるという耳痛世直し論文です。「心理学的研究における」とありますが、重回帰分析というあらゆる問題に対して適用される手法を題材にしているため、広い読者層に対して非常に学びのあるものだと思います。

まず、こちらの論文の概要を簡単に紹介してから指摘されている事項を列挙し、それぞれ例を交えながら解説します。最後に、簡単に検証できそうなものを選び、実際に誤用した場合どんな問題が生じるかを確認して行きたいと思います。

用語について

様々な表現がある用語について、ここではどの言い方を使うのか、またそれらの意味を簡単に説明します。基本的に論文に準拠した用語を使います。

独立変数...気温、風量、カテゴリなど、何かを数値で表したもの。多くの場合、ある関心のある量に対する影響や相関を見るために収集される。別名：特徴量、説明変数
従属変数...独立変数によって予測・説明される変数。別名：目的変数、基準変数、被説明変数
交絡変数...従属変数、独立変数の双方に対して相関を持つとされる変数。別名：共変量、統制変数
媒介変数...「独立変数→媒介変数→従属変数」といったように、独立変数と従属変数を媒介する変数。独立変数と従属変数の関係をより詳細に表現するために使われる。別名：中間変数、メディエータ

概要と誤用の指摘

概要

心理学的研究に関する専門誌「心理学研究」に掲載されている重回帰分析を行なった論文に対し、「説明」と「予測」という観点から問題点を指摘する。また、「説明」は因果関係を追求することとし、主に「交絡変数を考慮しているか」「逆方向の因果は存在するか」「媒介変数と独立変数を分離しているか」という3点に対する議論が十分になされているか否かを問題の基準とする。

誤用の指摘（登場した順）

交絡変数はなるべく多くモデルに組み込むべきであるが、これに対する議論が足りていない
...極端な話、「アイスの売上高と熱中症患者数に正の相関がある」という結果をえたとき、「気温」というどちらに対しても影響している変数を考慮すべきであるのにもかかわらず、どんな交絡変数がありうるかという議論が弱い、または無いといった指摘です。
逆の因果が考えられるのにもかかわらず、一方向のみのを想定した分析・議論になっている
...ニワトリタマゴ問題に似たような状況を感じます。人口増加率と GDP 成長率、住宅の価格と部屋数（予算としてある程度価格が決まっている場合など）、など色々な状況で直面すると思います。分析のデザイン段階で、理論的に因果方向があっているかどうかは気をつけて議論すべき点です。
間接効果と直接効果を分離するために、媒介変数を独立変数の1つとして扱うべきではない
...これを全て無視することで、風が吹いたら桶屋が儲かります。便宜上、とりあえず「風量と桶屋の売上には正の相関がある」とします。このとき、「風量」→「桶屋の売上」という直接効果のみを評価していることになりますが、本来は「風量」→...→「猫の数」→「ネズミの数」→「桶屋の売上」といった具合で、別の変数を介して効果が波及します（間接効果）。媒介変数を独立変数と並列に扱わず、きちんと媒介変数として使うことで因果関係をより詳細に評価し、「風量」→「桶屋の売上」が有意でないというより現実的な結果が得られることが期待できます。ちなみに、媒介変数によって有意だった直接効果が有意でなくなる場合は「完全媒介モデル」とよばれるみたいです。
独立変数と従属変数の内容が重複してる
...本文では「サイコパシー特性」で「非道徳的意図ないし行動」を回帰する例を挙げて指摘しています。従属変数と独立変数の重複は、いわば同じような対象同士の相関を見ているため、なんら示唆に繋がらないということだと思います。
心理学的研究の一般的な目的としては個人内変動に基づく検討を行う必要があるのに、大体は個人間変動を調べている
...「〇〇が高い人ほど□□が高い傾向がある」といった個人間での比較ではなく、個人内での変動因に基づく検討を行う必要がある、という指摘。例えば、「サイコパシー特性を有する人はお団子が好きな傾向がある」ではなく、「サイコパシー特性を有することは食の好みにどう影響するか」を調査すべき、といったことでしょうか。横断的に相関を見て傾向を述べるのではなく、縦断的に個人の変化を調査すべきだという主張だと捉えています（自信ないので、ご指摘お願いします、、）。
$j$ 番目の偏回帰係数は「$j$ 番目の独立変数を他の独立変数で回帰した結果の残差」と「従属変数を $j$ 番目以外の独立変数で回帰した結果の残差」の関係を表すが、その変数そのものの値と目的変数との関係とみなしている
...要するに偏相関として考えるか、単相関として考えるかです。もし独立変数同士がその名の通り全て互いに独立であれば偏相関と単相関は一致するはずですが、そんな状況は現実的にはほぼないと思います。重回帰であれば重回帰なりに、偏回帰係数は他の変数の影響を取り除いた上での数値であることを踏まえて解釈を行うべきです。
重回帰での有意性から「効果が認められない」と結論づけているが、目的によっては単相関係数を見るべき場合もある
...検討したいことは「他の独立変数を統制したときの従属変数との関係」なのか、「独立変数そのものと従属変数との関係」なのかによって、単回帰の結果、重回帰の結果、または両方の結果を使い分けるべきという指摘です。単回帰の結果が重要となる問題設定ってどんなものですかね？これに関しては正直、理解と経験が追いついていません。
「A は B を予測するか」という題目にもかかわらず、「A は B を促進させることが確認された」と帰結しているが、これは予測と説明が混在した解釈であり、因果関係を特定したことになっている
...その通りです、これはついつい口が滑ってしまった方もいらっしゃるのではないでしょうか、、。本当に因果だと決めつけているかは指摘された論文を読んでいないのでわかりませんが、これは分析手法の誤用というか、言葉が定義であることを意識して慎重に使うべき、といったところでしょうか。気をつけます。
予測が目的なら変数選択してもいいが、説明が目的なら変数選択すべきでない
...交絡変数を入れたり、より目的変数の分散を表現するために変数選択すべきでないというのはなんとなくわかりますが、これは状況と変数選択方法にもよると思います。ちなみに、変数選択には解釈性を向上させるといった動機があると思いますが、解釈性は人間が見て予測の過程を理解できるかの度合いであって、説明性とは別の話です。
検定の際は有意性だけじゃなくて効果量を示すべき
...今後はますますこの傾向が強くなっていくのではないでしょうか。例えば $t$ 検定における検定統計量 $t$ は以下で表されます。
$$t = \frac{\bar{x} - \mu}{\sqrt{V / n}}$$
着目していただきたいのはサンプルサイズ $n$ だけなので他の定義は割愛しますが、式から分かる通りサンプルサイズを大きくすれば $t$ も大きくなり、帰無仮説が棄却されやすくなります。このことからも、どれだけの差（効果量）があったかを報告することが大切です。効果量の詳しい解説は別に譲りますが、ざっくりいうと効果量は帰無仮説で決められた量と検定統計量の間に生じた差の大きさです。ちなみに、効果量は逆算的に設計されるものでもあり、「これくらいの効果量があったとすれば、高い検出力で帰無仮説を棄却したい」といった設定のもと、それを満たすようなサンプルサイズを決めることができます。効果量が大きいのであれば、こんな小さいサンプル数でも帰無仮説を棄却できるような差が生まれるはず、といったイメージです。
測定値の単位の意味が定かでない場合が多いので、独立変数は標準化して標準偏回帰係数で解釈を行うべき
...測定値の単位の意味に関わらず、スケールを揃えた方が他の変数と比較しやすいと思います。ちなみに、線形重回帰の場合、予測精度という観点では標準化はしてもしなくても変わりません。ただし、正則化を使う際はスケールを合わせて回帰係数に対する罰則をフェアにするために標準化は必須です。

目立つものを選びましたが、本文中にはまだまだ指摘があります。
次に、「ちゃんと交絡変数を組み込まないとどうなるか」という検証を、簡単な設定で行たいと思います。

実際に、間違ってみた。

分析の設定

「$y$ を説明するために変数 $x_1, x_2, x_3$ を収集して回帰分析したが、実は $c_1, c_2$ という交絡変数が存在した」

従属変数 $y_i$ は、ある交絡変数 $c_1, c_2$ と独立変数 $x_3$ の線形和に誤差が加わった値とします（真のモデル）。$c_1, c_2$ はそれぞれ $x_1, x_2$ と交絡しています。設定を数式で書いたものが以下になります。

$y_i = \frac{1}{3}c_{i1} + \frac{1}{3}c_{i2} + \frac{1}{3}x_{i3} + \epsilon_{i1}, \quad \epsilon_{i1} \sim\mathcal{N}(0, 0.5^2),$
$c_{i1}, c_{i2} \sim \mathcal{N}(0, 1^2),$
$x_{i1} = c_{i1} + \epsilon_{i2}, \quad \epsilon_{i2}\sim \mathcal{N}(0, 1^2),$
$x_{i2} = c_{i2} + \epsilon_{i3}, \quad \epsilon_{i3} \sim \mathcal{Po}(1),$
$x_{i3} \sim \mathcal{U}(-\sqrt{3}, \sqrt{3}),$
$i=1,\ldots, 100$

$x_1 \sim x_3$ はそれぞれ正規分布、ポアソン分布、一様分布に従い、分散はそれぞれ1にしています。変数間の関係を図で表すと以下のような感じです。

例えば $y$ がアイスの売上高、 $x_1, x_2, x_3$ がそれぞれ熱中症患者数、スポーツドリンクの売上高、アイスの広告費、$c_1, c_2$ が気温、人口密度、といった感じです（例えです、上記の数字とは整合しません）。

この設定のもと、交絡変数を入れない場合、また入れた場合はどんなことが起こるのかを見ていきます。なお、検証の主旨とは関係しないため、変数の標準化は行っていません。

R を使った検証

まずライブラリの読み込みと、1. 交絡変数、2. 独立変数、3. 従属変数、のそれぞれを乱数を使って生成します。

library(tidyverse)
set.seed(0)
# 交絡変数
C = tibble(c1 = rnorm(100, mean = 0, sd = 1), 
               c2 = rnorm(100, mean = 0, sd = 1))
# 独立変数
X = tibble(x1 = C$c1 + rnorm(100, mean = 0, sd = 1), 
               x2 = C$c2 + rpois(100, lambda = 1), 
               x3 = runif(100, min = -sqrt(3), max = sqrt(3)))
# 従属変数：3つの変数を均等に重み付けした値に誤差を加える
y = C$c1 / 3 + C$c2 / 3 + X$x3 / 3 + rnorm(100, mean = 0, sd = 0.5)

data = tibble(y, X, C)

出来上がった変数を確認するために、相関係数や分布、散布図と単回帰による直線を見てみます。

library(GGally)
fig_pairs = ggpairs(data = data, 
                    lower = list(continuous = wrap("smooth", size = 1.5)),
                    upper=list(continuous=wrap("cor",size = 7)), 
                    diag = list(continuous = wrap('barDiag', bins = 30))) + 
  theme(axis.text = element_text(size = 15), strip.text = element_text(size = 20)) 
ggsave(file = 'pairs.png', plot = fig_pairs, dpi = 400, width = 10, height = 10)

上三角成分：相関係数と、相関係数の信頼区間をもとにした有意性
対角成分：ヒストグラム（相対頻度）
下三角成分：散布図と回帰直線+信頼区間

係数や散布図を見ると、どの独立変数も従属変数 $y$ と正の相関があります。交絡変数である $c_1, c_2$ は $y$ および $x_1, x_2$ とそれぞれ相関していることが確認できます（そういう作り方したので当然ですが）。
この時点で安直に判断すれば「多重共線性っぽさもなく $x_1 \sim x_3$ 全部が $y$ と関係しそう」となります。

次に、1. 交絡変数を入れなかった場合（間違い例）、2. 全ての変数を使った場合（正解例）、3. 交絡変数のみを使った場合、4. 正解モデルの場合、の4通りについて、偏回帰係数の有意性や自由度調整済み決定係数を見てみます。自由度調整済み決定係数は、説明変数を増やせば増やすほど大きくなってしまう決定係数を、説明変数の数をペナルティとして調整したものです。なお、結果は今回関係するところに着目できるように一部省略しています。

1. 交絡変数を入れなかった場合（間違いの例）

summary(lm(y ~ x1 + x2 + x3, data = data))

# 実行結果(一部省略)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.23651    0.07984  -2.962 0.003848 ** 
x1           0.20288    0.04648   4.365  3.2e-05 ***
x2           0.18320    0.04713   3.887 0.000187 ***
x3           0.26316    0.06615   3.978 0.000135 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.3837,	Adjusted R-squared:  0.3645

真に関係するのは $x_3$ のみですが、一応見かけ上は相関する $x_1, x_2$ の変回帰係数もしっかり0.1%水準で有意となっています。ただし、自由度調整済み決定係数（'Adjusted R-squared'）は0.3645と、分野にもよりますがあまり高いとはいえない値となっています。
このモデルを信じるとすれば $x_1, x_2$ という誤った変数に因果関係を言及してしまう可能性があります。

2. 全ての変数を使った場合（交絡に対するより適切な対策）

しっかり $c_1, c_2$ の存在に気づき、モデルに入れて変数を統制しようと改善させたケースです。

summary(lm(y ~ ., data = data))

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.006543   0.080848   0.081  0.93567    
x1          -0.091558   0.059145  -1.548  0.12498    
x2          -0.029679   0.052564  -0.565  0.57367    
x3           0.333805   0.060440   5.523 2.95e-07 ***
c1           0.333411   0.078991   4.221 5.61e-05 ***
c2           0.244881   0.072350   3.385  0.00104 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.395,	Adjusted R-squared:  0.3629

統制を行い、見事 $x_1, x_2$ のみ有意でないという結果となり、最終的な解釈としては正しいものとなっています。間違いモデルと比較すると、決定係数はより高く、自由度調整済み決定係数はより低くなっています。サンプル数100程度だと2つの変数の増加は大きくモデル評価指標に影響することが分かると同時に、使用する変数の数が違うのであれば自由度調整済みで見ることの大切さも伺えます。また、間違いモデルと係数の符号を比較すると、$x_1, x_2$ の符号が逆転していることがわかります。これも統制を行った際によくみられる現象です。

3. 交絡変数のみを使った場合（おまけ1）

論文で指摘されていることとはずれますが、比較として交絡変数のみの場合も見てみます。

summary(lm(y ~ c1 + c2, data = data))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.07855    0.06409  -1.226 0.223351    
c1           0.17831    0.06634   2.688 0.008464 ** 
c2           0.26797    0.06787   3.948 0.000149 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.1749,	Adjusted R-squared:  0.1579

これは少し意外な結果でしたが、ちゃんと有意にはなっているものの、 $c_1$ に関しては5%水準、自由度調整済み決定係数は0.1579とかなり低いです。ただ、決定係数の低さは $x_3$ を落としている点、そもそも2変数しか使っていない点からも納得です。

4. 正解モデルの場合（おまけ2）

これは本来気付けるはずもないので、あくまで比較のためですが。

summary(lm(y ~ x3 + c1 + c2, data = data))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.02983    0.05664  -0.527 0.599641    
x3           0.33058    0.05916   5.588 2.15e-07 ***
c1           0.25538    0.05955   4.289 4.28e-05 ***
c2           0.22399    0.05979   3.746 0.000306 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.3774,	Adjusted R-squared:  0.3579

（自由度調整済み決定係数が間違いモデルよりも低い、、、頑張れ真のモデル、、、）
これも意外だったのが、2のモデルよりも真の偏回帰係数から遠いという点。正規乱数のせいといえばそれまでですが、ここは理解できていない部分です。

まとめ

重回帰分析の誤用・誤解に対する注意喚起となる論文を紹介しました。この内容を踏まえ、あえて重回帰分析を誤った方法で使うことで起こりうる問題を確認しました。
「先人たちの行いを鵜呑みにするな」といった内容の記述が文中にありますが、鵜呑みにしないようにするためには間違いや違和感に気づくための理論的な基礎知識が重要になるのではないかと思います。