見出し画像

『医学のための因果推論』のための統計学入門：7.デザイン行列とコーディング

Hikaru Ooba

『医学のための因果推論Ⅰ　一般化線形モデル』を勉強した時のまとめです。
結構難しい本でしたので自分なりの解釈を入れながらまとめました。

統計モデルは、データの中からパターンを見つけ出し、要因間の関係を理解するための重要なツールです。
気象予報や医学研究など、日常生活の多くの場面で統計モデルが活用されています。

特に一般化線型モデル（Generalized Linear Model; GLM）は多様なデータタイプや状況に適用できる柔軟なモデルとして広く利用されています。
GLMを正しく構築し、解釈するために必要なものがデザイン行列とコーディングの概念です。


デザイン行列とコーディング

デザイン行列

デザイン行列は、統計モデルにおいてデータ中の説明変数（独立変数、共変量とも呼ばれる）を整理して行列表現にしたものです。
各行は個々の観測データ（対象者や実験単位）を表し、各列は説明変数を表します。

デザイン行列を使うことで、複数の説明変数を含むモデルを行列の形で簡潔に表現できます。
また、線形代数の手法を用いて計算を効率化できます。

例として、ある薬の効果を調べるために100人の被験者に対して実験を行ったとします。
被験者の年齢と薬の投与量が結果に影響を与えると考え、これらを説明変数としてモデルに組み込みたいとします。

被験者ごとに年齢 $${X_{i1}}$$​ と投与量 $${X_{i2}}$$のデータがある場合、デザイン行列 $${\boldsymbol{X}}$$ は次のようになります。

$$
\boldsymbol{X} = \begin{pmatrix} 1 & X_{11} & X_{12} \\ 1 & X_{21} & X_{22} \\ \vdots & \vdots & \vdots \\ 1 & X_{N1} & X_{N2} \end{pmatrix}
$$

ここで、最初の列の「1」は切片項を表し、モデル全体の平均効果を捉えるために使われます。

コーディング

コーディングは、デザイン行列を構築する際に説明変数をどのような数値で表現するかを決めることです。
特にカテゴリカルデータ（分類データ）を数値に変換する際に重要になります。

適切なコーディングを行うことで、モデルの結果を正しく解釈でき、統計解析がスムーズに進みます。

コーディングの具体例

  • 連続データの場合：
    そのまま数値としてデザイン行列に含めます。
    ただし単位やスケールによっては平均を0に、標準偏差を1に調整するなどの前処理を行うことがあります。

  • カテゴリカルデータの場合：
    グループやカテゴリーを数値に変換します。
    例えば性別を「男性=0」「女性=1」とするなどのダミー変数を用いてコーディングします。

一般化線型モデル

一般化線型モデル（GLM）は、従属変数（目的変数）の種類に応じて適切な確率分布とリンク関数を用いることで線形回帰モデルを拡張したものです。

GLMは以下の3つの要素で構成されます。

確率分布：従属変数 $${Y}$$ の分布
（例：正規分布、二項分布、ポアソン分布）

線形予測子：説明変数とパラメータの線形結合

$$
\eta_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_p X_{ip}
$$

リンク関数：線形予測子と従属変数の期待値を結びつける関数

$$
g(\mu_i) = \eta_i
$$

モデルの数式表現

個々の観測値 $${i}$$について、一般化線型モデルは次のように表されます。

$$
{g\left( \mathrm{E}\left( Y_i \mid \boldsymbol{X}_i \right) \right) = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_p X_{ip}}
$$

ここで

$${\mathrm{E}\left( Y_i \mid \boldsymbol{X}_i \right)}$$は従属変数の条件付き期待値です。

$${g(\cdot)}$$ はリンク関数で、従属変数の期待値と線形予測子を結びつけます。

連続データの扱い

単回帰モデル

最も基本的なモデルとして、連続変数を1つだけ用いる単回帰モデルがあります。
例えば、身長 $${X}$$ が体重 $${Y}$$ に与える影響を調べたいとします。

リンク関数として恒等関数 $${g(x) = x}$$ を用いると、モデルは次のようになります。

$$
\mathrm{E}\left( Y_i \mid X_i \right) = \beta_0 + \beta_1 X_i
$$

デザイン行列とパラメータベクトルは以下のようになります。

$$
\boldsymbol{X} = \begin{pmatrix} 1 & X_1 \\ 1 & X_2 \\ \vdots & \vdots \\ 1 & X_N \end{pmatrix}, \quad \boldsymbol{\beta} = \begin{pmatrix} \beta_0 \\ \beta_1 \end{pmatrix}
$$

高次の項や交互作用項の導入

高次の項

データが直線的な関係でない場合、変数の2次の項（またはそれ以上）をモデルに含めることで、より適切なフィッティングが可能になります。

2次の項を含むモデル

$$
\mathrm{E}\left( Y_i \mid X_i \right) = \beta_0 + \beta_1 X_i + \beta_2 X_i^2
$$

この時、デザイン行列は以下となります。

$$
\boldsymbol{X} = \begin{pmatrix} 1 & X_1 & X_1^2 \\ 1 & X_2 & X_2^2 \\ \vdots & \vdots & \vdots \\ 1 & X_N & X_N^2 \end{pmatrix}
$$

交互作用項

複数の変数間の相互作用をモデルに組み込みたい場合、交互作用項を導入します。

変数 $${X_1}$$​ と $${X_2}$$​ の交互作用を含むモデル

$${\mathrm{E}\left( Y_i \mid X_{i1}, X_{i2} \right) = \beta_0 + \beta_1 X_{i1} +\beta_2 X_{i2} + \beta_3 X_{i1} X_{i2}}$$​

この場合、デザイン行列は以下になります。

$$
\boldsymbol{X} = \begin{pmatrix} 1 & X_{11} & X_{12} & X_{11} X_{12} \\ 1 & X_{21} & X_{22} & X_{21} X_{22} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & X_{N1} & X_{N2} & X_{N1} X_{N2} \end{pmatrix}
$$

分類データの扱いとダミー変数

カテゴリカルデータの数値化

カテゴリカルデータ（例：性別、地域、グループなど）は、そのままでは数値として扱えないため、ダミー変数を用いて数値化します。

ダミー変数

ダミー変数は、カテゴリカルデータの各カテゴリーを0または1などの数値で表現したものです。
例えば、性別を「男性=0」「女性=1」とコーディングします。
複数のカテゴリーがある場合は、それぞれのカテゴリーに対応するダミー変数を作成します。

具体例：2グループの比較

コントロール群と試験治療群の2つのグループで平均の差を比較したいとします。
この時、デザイン行列、パラメータベクトル、モデル式は以下のようになります。

デザイン行列：

$$
\boldsymbol{X} = \begin{pmatrix} 1 & 0 \\ 1 & 0 \\ \vdots & \vdots \\ 1 & 1 \\ 1 & 1 \\ \vdots & \vdots \end{pmatrix}
$$

パラメータベクトル：

$$
\boldsymbol{\beta} = \begin{pmatrix} \beta_0 \\ \beta_1 \end{pmatrix}
$$

モデル式：

$$
\mathrm{E}\left( Y_i \mid \text{group} \right) = \beta_0 + \beta_1 X_{\text{group}, i}
$$

ここで、$${X_{\text{group}, i} = 0}$$ ならコントロール群、$${X_{\text{group}, i} = 1}$$ なら試験治療群を表します。

この場合、$${\beta_0}$$​ はコントロール群の平均、$${\beta_1}$$​ はグループ間の平均差（試験治療群の効果）と解釈することができます。

コーディングのポイント

  • 基準カテゴリーを設定する：一つのカテゴリーを基準（参照）として、他のカテゴリーをそれに対する差として表現します。

  • ダミー変数の数：カテゴリー数が $${k}$$ のときは、ダミー変数の数は $${k - 1}$$ 個とします。
    すべてのカテゴリーに対してダミー変数を作ると、モデルが正しく動作しません（多重共線性）。

モデルの適用と解釈

デザイン行列の重要性

デザイン行列の構築方法（コーディング）は、モデルの結果解釈に直接影響します。
同じデータでもデザイン行列のコーディングが異なれば、パラメータの意味や推定される値が変わります。

適切なコーディングの選択においては、合目的性(研究の目的や仮説に応じて、どのようなコーディングが適切か)と解釈性(パラメータが何を表しているか)を明確に理解できるようにすべきです。

実データを用いた分析

簡単なデータを用いて、デザイン行列と一般化線型モデルの概要を見てみます。

データ

  • サンプルサイズ：100

  • 説明変数：

    • $${X_1}$$​：連続変数（ここでは年齢とします）

    • $${X_{\text{group}}}$$​：カテゴリカル変数（ここでは0：コントロール群、1：試験治療群とします）

  • 従属変数：

    • $${Y}$$：連続変数（ここでは血圧とします）

デザイン行列

カテゴリカル変数（グループ）と連続変数（年齢）を含むデザイン行列は以下のように表されます：

$$
\boldsymbol{X} = \begin{pmatrix} 1 & X_{11} & X_{\text{group}, 1} \\ 1 & X_{12} & X_{\text{group}, 2} \\ \vdots & \vdots & \vdots \\ 1 & X_{1N} & X_{\text{group}, N} \end{pmatrix}
$$

1 ：切片項を表す列

$${X_{1i}}$$​ ：各対象者の年齢を表す連続変数

$${X_{\text{group}, i}}$$​ ：カテゴリカル変数（グループ：0または1）を示すダミー変数

パラメータベクトル

それぞれの説明変数に対応するパラメータは次のようになります：

$$
\boldsymbol{\beta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{pmatrix}
$$

$${\beta_0}$$：切片項、すなわちコントロール群の年齢0歳のときの平均値（基準点）。

$${\beta_1β}$$ ：年齢が1歳増加したときの平均的な変化量（年齢の効果）。

$${\beta_2}$$​：グループが変わったときの平均的な変化量（グループ効果）。

モデル式

モデル式は以下のように表現できます。

$$
\mathrm{E}(Y_i \mid X_{1i}, X_{\text{group}, i}) = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{\text{group}, i}
$$

この式は、年齢とグループの両方を考慮した従属変数 $${Y_i}$$ の条件付き期待値を示します。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# データの生成
np.random.seed(0)
N = 100
X1 = np.random.normal(50, 10, N)  # 年齢
group = np.random.binomial(1, 0.5, N)  # グループ（0 or 1）

# 真のパラメータ
beta_0 = 120  # 切片
beta_1 = 0.5  # 年齢の効果
beta_2 = -5   # 試験治療群の効果

# 従属変数の生成
Y = beta_0 + beta_1 * X1 + beta_2 * group + np.random.normal(0, 5, N)

# デザイン行列の作成
X = pd.DataFrame({
    'const': 1,
    'X1': X1,
    'group': group
})

# モデルの適合
model = sm.OLS(Y, X)
results = model.fit()

# 結果の表示
print(results.summary())
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# データ生成
np.random.seed(0)
N = 100
X1 = np.random.normal(50, 10, N)  # 年齢
group = np.random.binomial(1, 0.5, N)  # グループ（0 or 1）

# 真のパラメータ
beta_0 = 120  # 切片
beta_1 = 0.5  # 年齢の効果
beta_2 = -5   # 試験治療群の効果

# 従属変数の生成
Y = beta_0 + beta_1 * X1 + beta_2 * group + np.random.normal(0, 5, N)

# デザイン行列の作成
X = pd.DataFrame({
    'const': 1,
    'X1': X1,
    'group': group
})

# モデルの適合
model = sm.OLS(Y, X)
results = model.fit()

# 回帰直線の描画
plt.figure(figsize=(8, 6))
plt.scatter(X1[group == 0], Y[group == 0], label='Control Group', alpha=0.7, color='blue')
plt.scatter(X1[group == 1], Y[group == 1], label='Treatment Group', alpha=0.7, color='orange')

# コントロール群の回帰直線
X_pred_control = pd.DataFrame({
    'const': 1,
    'X1': np.linspace(X1.min(), X1.max(), 100),
    'group': 0  # コントロール群
})
Y_pred_control = results.predict(X_pred_control)
plt.plot(X_pred_control['X1'], Y_pred_control, color='blue', label='Regression Line (Control Group)')

# 試験治療群の回帰直線
X_pred_treatment = pd.DataFrame({
    'const': 1,
    'X1': np.linspace(X1.min(), X1.max(), 100),
    'group': 1  # 試験治療群
})
Y_pred_treatment = results.predict(X_pred_treatment)
plt.plot(X_pred_treatment['X1'], Y_pred_treatment, color='orange', label='Regression Line (Treatment Group)')

plt.xlabel('Age')
plt.ylabel('Blood Pressure')
plt.title('Linear Regression with Group Effect')
plt.legend()

# グラフの保存と表示
plt.savefig('linear_regression_group.png', dpi=300, bbox_inches='tight')
plt.close()
from IPython.display import Image
Image('linear_regression_group.png')

この場合、以下のような解釈ができます。

  • $${\beta_0}$$​：コントロール群の、年齢が0歳のときの予測血圧(切片項)。119.10 ± 2.50mmHg (95%信頼区間：114.14 - 124.07mmHg)

  • $${\beta_1}$$​：年齢が1歳増加すると、血圧が平均で何mmHg増加するか。0.51 ± 0.05mmHg (95%信頼区間：0.41 - 0.60mmHg)

  • $${\beta_2}$$​：試験治療群はコントロール群に比べて平均で何mmHg血圧が低いか。-5.61 ± 0.96mmHg (95%信頼区間：-7.51 - -3.71mmHg)

限界

多重共線性

デザイン行列の中で列同士が線形に関連していると、パラメータの推定が不安定になります。
これを多重共線性といいます。

対策：不要な変数を除外するか、コーディングを見直す必要があります。

コーディングの影響

コーディングが適切でないと、パラメータの解釈が難しくなります。

他の人が同じデータを分析したときに、再現性が確保されていないと異なる結果になる可能性があります。

参考文献

『医学のための因果推論Ⅰ　一般化線形モデル』

いいなと思ったら応援しよう！