中断時系列分析における多重共線性の考察

Interrupted Time Series Analysis(中断時系列分析)でよく利用される segmented regression のパラメータについて、なぜ従来のモデル（Equation 1）が多重共線性の問題を引き起こすのか、そしてその問題をどのように解決するのかを数学的な観点から解説します。

以下のアンサー記事のようなものになっています。
https://qiita.com/Jwellyfish/items/efd3a61fb7b7eec51545

本記事作成にあたって以下の論文を参考にしています。(この論文には多重共線性ではなく、モデルのバイアスについて記載されています)
https://pubmed.ncbi.nlm.nih.gov/33097937/

モデルの定式化と問題の概要

まず、一般的なモデル設定として、次の Equation 1 を考えます。

$$Y_t = \beta_0 + \beta_1 T + \beta_2 X_t + \beta_3 (T \times X_t)$$

ここで、

$T$は経過時間を表します。
$X_t$は介入前は 0、介入後は 1 となるダミー変数です（介入開始時刻を$T_i$ とする）。
$T \times X_t$ は$T$と$X_t$の交互作用項です。

多重共線性とランク落ちの数学的意味

交互作用項の性質

介入前は$X_t = 0$ のため、交互作用項は

$$T \times X_t = 0$$

となりますが、介入後 $X_t = 1$ では

$$T \times X_t = T \times 1 = T$$

となります。
すなわち、介入後のデータにおいては、変数$T$ と交互作用項は常に同じ値となります。
より形式的には、介入後の観測に対して $Z = T \times X_t$ とおくと、

$$Z = T$$

となります。このため、モデル内でこれら二つの変数は同じ情報を持つことになり、数学的には $T \times X_t$は $T$ の 1 倍（定数項 0 の線形結合）で表現できるため、完全な線形依存が生じます。

ただし介入後 $X_t = 1$から、

$$Z = T$$

となるので介入前 $X_t = 0$では線形依存は発生しません。

rank落ちによる問題

回帰分析では、デザイン行列が $X$ として表されますが、もしある説明変数が他の変数の完全な線形結合で表される場合、デザイン行列はフルランクではなくなります。
具体的には、Equation 1 の場合、説明変数として

定数項 $1$
時間変数 $T$
介入ダミー $X_t$
交互作用項 $T \times X_t$

が含まれますが、介入後では $T \times X_t$ = $T$ となるため、$T$ と $T \times X_t$ の両方が同じ情報を持ちます。

その結果、行列のランクが本来の列数よりも低くなり、逆行列を求められなくなります。
このrank落ちは、回帰係数が一意に推定できない（識別できない）状態を招き、推定の信頼性に深刻な影響を与えます。

再パラメータ化による解決策

上記の問題を解決するため、交互作用項を介入開始時点 $T_i$ で中心化する再パラメータ化を行います。具体的には、次の恒等式を利用します。

$$T \times X_t = (T-T_i)X_t + T_iX_t$$

これを Equation 1 に代入すると、

$$Y_t = \beta_0 + \beta_1 T + \beta_2 X_t + \beta_3 [(T-T_i)X_t + T_iX_t ]$$

となり、整理すると、

$$Y_t = \beta_0 + \beta_1 T + (\beta_2 + \beta_3 T_i) X_t + \beta_3 (T-T_i)X_t$$

となります。ここで新たに、

$$\beta'_2 = \beta_2 + \beta_3 T_i,\quad \beta'_3 = \beta_3$$

と定義すると、モデルは以下の Equation 2 として書き換えられます。

$$Y_t = \beta_0 + \beta_1 T + \beta'_2 X_t + \beta'_3 (T-T_i)X_t$$

この形式では、介入直後（$T = T_i$）に

$$(T-T_i)X_t = 0$$

となるため、即時効果が正しく $\beta'_2$ により評価されます。さらに、中心化された交互作用項 $(T-T_i)X_t$ は $T$ とは独立した変動を持つため、行列がフルランクに近い状態となります。これにより、推定が一意に行えるようになり、ランク落ちの問題が解消されます。

シミュレーション

library(car)

# シミュレーションパラメータ
T <- 0:150             # 0から150までの時間
T_i <- 75             # 介入開始時刻
X <- ifelse(T >= T_i, 1, 0)  # 介入ダミー：T < T_i は 0、T >= T_i は 1

# Equation 1 の交互作用項: T * X
interaction1 <- T * X

# Equation 2 の交互作用項: (T - T_i) * X
interaction2 <- (T - T_i) * X

# Equation 1 のデザイン行列 (Intercept は lm() が自動で追加)
df_eq1 <- data.frame(
  T = T,
  X = X,
  Interaction = interaction1
)

# Equation 2 のデザイン行列
df_eq2 <- data.frame(
  T = T,
  X = X,
  Interaction = interaction2
)

# ダミーの応答変数を生成（例：正規乱数）
set.seed(123)
df_eq1$Y <- rnorm(nrow(df_eq1))
df_eq2$Y <- rnorm(nrow(df_eq2))

# 線形モデルのフィッティング
model_eq1 <- lm(Y ~ T + X + Interaction, data = df_eq1)
model_eq2 <- lm(Y ~ T + X + Interaction, data = df_eq2)

 # VIF の計算
vif_eq1 <- vif(model_eq1)
vif_eq2 <- vif(model_eq2)

cat("Equation 1 の VIF:\n")
print(vif_eq1)

cat("\nEquation 2 の VIF:\n")
print(vif_eq2)     

以上からEquation2ではVIFが改善されているのが分かると思います。

まとめ

以上、Interrupted Time-Series Analysisにおける多重共線性について考察しました。
すごく簡単な式展開でわかる話でしたが、知らないと間違えたモデリングをしてしまいそうですね。
ジャーナルに乗ってるからと言ってそれが正しい保証はない典型的な例でした。
自戒。