Python3ではじめるシステムトレード：分散分析と直交表

統計検定準一級の内容は２級と比べて格段に難しくなります。実験計画法では２級ではない直交表の習得が求められます。また、実験計画法の内容も乱塊、交絡など難しい概念が現れ困惑するのではないでしょうか？また、実務経験を積んだ人も多く、経験が数学的概念の受け入れを邪魔すのではないでしょか？私も困惑しました。そこでまとめてみました。

間違いがあると思うのでコメントをいただけると助かります。
2023/8/25 1270

基本的な理解

実験の結果とその要因の関係を推定することを考えます。実験の結果を応答($y$)といい、応答に影響を及ぼすと考えられる変数の内で、実際に実験で操作をするものを因子といいます。AとかBで表します。因子の条件を水準といいます。これは$A_1,...A_a$で表します。主効果は、一つの因子が応答に持つ効果を指します。これは、他の因子が一定の条件である場合に、ある因子が応答にどれくらいの影響を与えるかを示すものです。

主効果をAとBとし、その水準をそれぞれ1,2として分散分析を行います。応答$y$は正規分布にしたがうと仮定します。

$y_{i,j}\sim N(\mu_{i,j},\sigma^2)$

これを線形の方程式で表すと

$\mu_{i.j}=\mu+\alpha_i+\beta_j+\gamma_{i,j}$

となります。$\mu$は観測値の総平均、$\alpha_i$は主効果Aの水準$i$の効果(additive effect)、$\beta_j$は主効果Bの水準$j$の効果(additive effect)、$\gamma_{i,j}$は交互作用を表し因子Aの水準$i$の影響が因子Bの水準$j$により変わるような効果を示すものです。一般に$(A_iB_j)$の形で表現されます。しかし、交互作用効果がない場合、この項はモデルから省略されます。

応答の式はノイズを含めて次のようになります。

$y_{i,j}=\mu_{i,j}+\varepsilon_{i,j}\text{ with } \varepsilon_{i,j}\sim N(0,\sigma^2)$

$\mu$: すべての観測データの平均値です。
$\alpha_i$: これは、因子Aの水準$i$が応答にどれだけの影響を与えるかを示します。
$\beta_j$: これは、因子Bの水準$j$が応答にどれだけの影響を与えるかを示します。
$\gamma_{i,j}$: これは、因子AとBの線形関係を超えた影響を示すものです。
$\varepsilon_{i,j}$: これは、他の要因やランダムな誤差など、モデルで説明できない変動を捉えるためのものです。

また、データ点が、調査中の科学的問題に関して適切であることが求められ、応答変数の平均が、因子によって加法的（交互作用項がない場合）かつ直線的に影響されることも仮定されています。また、パラメータの識別可能性を確保するために、以下の制約を加えることがあります。

$\sum_i \alpha_i=\sum_j \beta_j=\sum_i \gamma_{i,j}=\sum_j \gamma_{i,j}=0$

ここで、パラメータの識別可能性（identifiability）とは、統計学や機械学習において、あるモデルのパラメータがデータから一意に決定される性質を指します。つまり、同じデータセットに対して、推定されたパラメータが一意であるということです。

水準の効果の和がゼロ、すなわち $\sum_i \alpha_i=0$ とする制約を加えることによって、モデルのパラメータを一意に識別することができます。この制約がない場合、モデルのパラメータは複数の異なる値を取ることがあり、それによって統計的推論が困難になります。

これを例で説明してみましょう。線形モデルが以下のように表されるとします：

応答の平均（水準1） = 基本効果 + $\alpha_1$
応答の平均（水準2） = 基本効果 + $\alpha_2$
応答の平均（水準3） = 基本効果 + $\alpha_3$

この場合、$\alpha_i$ の和がゼロであるとする制約を加えないと、$\alpha_i$ と基本効果が同時に任意の値を取ることができます。例えば、基本効果が1増えて、$\alpha_i$ がすべて1減っても、応答の平均は変わりません。このような状況を回避するために、制約を加えることで、基本効果と水準の効果が一意に識別されるようにします。

統計的な推定の文脈では、識別可能性はパラメータが一意に決定されることを保証する重要な性質です。パラメータが識別不可能な場合、同じデータセットを持っていても、異なるパラメータセットが同じように良くデータを説明することがあります。この場合、どのパラメータセットが「真の」ものであるかを判断することができないため、統計的な推論が不確実になります。

主効果、交互作用の有無は、総平方和、$A$間平方和$S_A$、$B$間平方和$S_B$、誤差平方和$S_E$、交互作用$A \times B$による平方和$S_{AxB}$を用いた分散分析により判断します。

トイモデル

ANOVA（Analysis of Variance、分散分析）では、群間平方和（Between-group Sum of Squares、SSB）と群内平方和（Within-group Sum of Squares、SSW）を計算します。この二つの平方和はデータ内の変動（分散）の元を識別するのに役立ちます。

群間平方和（SSB）

群間平方和は各群（または処理、カテゴリなど）の平均と全体平均との差を示します。これは、各群が全体からどれだけ異なっているかを評価する指標です。

なぜ群間平方和が重要か？

効果の大きさ: 群間平方和が大きいと、各群の平均が全体の平均から大きく異なる可能性が高く、処理効果（またはカテゴリ間効果）が強いことを示す可能性があります。

統計的検定: 群間平方和はF検定に用いられ、群間の平均に有意な違いがあるかどうかを評価します。

データの解釈: 群間平方和が大きい場合、その変数（または処理、カテゴリ）が応答変数に与える影響が大きいと解釈されることが多くなる傾向があります。

以上の理由から、ANOVA分析では群間平方和を計算し、評価します。

ANOVA（Analysis of Variance、分散分析）において、AとBが因子で、それぞれが水準A1, A2, B1, B2を持っている場合、群はこれらの水準のすべての組み合わせになります。

具体的には、以下のような4つの群が形成されます。

A1とB1の組み合わせ
A1とB2の組み合わせ
A2とB1の組み合わせ
A2とB2の組み合わせ
これらの4つの群に分かれたデータに対してANOVAを行うことで、因子Aと因子Bそれぞれ、そして因子Aと因子Bの交互作用が応答変数（実験結果、アウトカムなど）に与える影響を解析することができるのです。

総平方和は

$S_T = \sum_i\sum_j(y_{i,j}-\bar{y})^2$
で表されます。これは

$= \sum_i\sum_j[(\mu_{A_i}-\bar{y})+(\mu_{B_j}-\bar{y})+(y_{i,j}-\mu_{A_i,B_j})+(\mu_{A_i,B_j}-\mu_{A_i}-\mu_{B_j}+\bar{y})]^2$
$=\sum_{i}\sum_j(\mu_{A_i}-\mu)^2=\sum_{i}\sum_j(\mu_{B_j}-\mu)^2+\sum_i\sum_j(y_{i,j}-\mu_{A_i,B_j})+\sum_i\sum_j(\mu_{A_iB_j}-\mu_{A_i}-\mu_{B_j}+\bar{y})^2$

$=S_A + S_B + S_E + S_{A \times B}$

となります。

これは

$\sum_i\sum_j[(\mu_{A_i}-\bar{y})(\mu_{B_j}-\bar{y})=\sum_i(\mu_{A_i}-\bar{y})\sum_j(\mu_{B_j}-\bar{y})=0$
(因子Aの各レベルの平均値$\mu A_i$と全体の平均値$\bar{y}$との差、および因子
B の各レベルの平均値$\mu A_i$と全体の平均値 $\bar{y}$との差をそれぞれ取り、これらの総和がゼロになる。)

$\sum_i\sum_j[(\mu_{A_i}-\bar{y})(y_{i,j}-\mu_{A_i}\mu_{B_j}+\bar{y})=\sum_i(\mu_{A_i}-\bar{y})\sum_j(\mu_{A_i,B_j}-\mu_{A_i}-\mu_{B_j}+\bar{y})=0$
(因子Aの各レベルの平均値$\mu_{A_i}$と全体の平均値$\bar{y}$との差に、交互作用による差分を掛け合わせた総和がゼロになる。)

$\sum_i\sum_j[(\mu_{B_j}-\bar{y})(y_{i,j}-\mu_{A_i}\mu_{B_j}+\bar{y})=\sum_j(\mu_{B_j}-\bar{y})\sum_j(\mu_{A_i,B_j}-\mu_{A_i}-\mu_{B_j}+\bar{y})=0$
(因子Bの各レベルの平均値$\mu_{B_i}$と全体の平均値$\bar{y}$との差に、交互作用による差分を掛け合わせた総和がゼロになる。)

より成り立ちます。

A1とA2を切削機械、B1,B2を切削工法とします。単位時間当たりの切削量データを

	B1	B2
A1	339	132
A2	138	173

とします。

$\mu=1/4\sum_i\sum_j y_{i,j}=(339+132+138+173)/4=195.5$
$\mu_{A_1}=1/2\sum_{1,j} y_{i,j}=(339+132)/2=235.5$
$\mu_{A_2}=1/2\sum_{2,j} y_{i,j}=(138+173)/2=155.5$
$\mu_{B_1}=1/2\sum_{i,1} y_{i,j}=(339+138)/2=238.5$
$\mu_{B_2}=1/2\sum_{i,2} y_{i,j}=(132+173)/2=152.5$

今回は、繰り返しの部分がないので平均は観測値そのものになります。

$\mu_{A_1,B_1}=1/1\sum_{1,1} y_{1,1}=339=339$
$\mu_{A_1,B_2}=1/1\sum_{1,2} y_{1,2}=132=132$
$\mu_{A_2,B_1}=1/1\sum_{2,1} y_{2,1}=138=138$
$\mu_{A_2,B_2}=1/1\sum_{2,2} y_{2,2}=173=173$

各種平方和は総平方和の分解に相当します。

$S_T=\sum_{i=1}^2\sum_{j=1}^2(y_{i,j}-\bar{y})^2$
$ =(339-195.5)^2+(132-195.5)^2+(138-195.5)^2+(173-195.5)^2=28,437$
$S_A=\sum_{i}\sum_j(\mu_{A_i}-\mu)^2=(235.5-195.5)^2+(155.5-195.5)^2=3,200$
$S_B=\sum_{i}\sum_j(\mu_{B_j}-\mu)^2=(238.5-195.5)^2+(152.5-195.5)^2=3,698$
$S_E=\sum_{i}\sum_j(y_{i,j}-\mu_{A_i,B_j})^2$
$ =(339-339)^2+(132-132)^2+(138-138)^2+(173-173)^2=0$

この観測値においては$S_T=S_A+S_B+S_E$になっていないので、交互作用があります。
それは

$S_{A\times B}=\sum_i\sum_j(\mu_{A_iB_j}-\mu_{A_i}-\mu_{B_j}+\bar{y})^2$

の存在を示しています。これに観測値を入れてみると

$S_{A \times B}=\sum_i\sum_j(\mu_{A_i,B_j}-\mu_{A_i}-\mu_{B_j}+\bar{y})^2$
$=(339-235.5-238.5+195.5)^2+$
$(132-155.5-238.5+195.5)^2+$
$(138-235.5-152.5+195.5)^2+$
$(173-155.5-152.5+195.5)^2$
$=14,641$

分散分析（ANOVA）において「群間」と「群内」は、データのバラつき（分散）を理解するための重要な概念です。

群間（Between-groups）

群間バラつきは、各群の平均値が全体の平均からどれくらい異なるかを示します。言い換えれば、各群（ここではA1B1, A2B1, A1B2, A2B2）の平均値と、全体の平均値との差を用いて計算されます。

全体の平均（全体の合計値/全体のサンプル数）:

(339+138+132+173)/4=196

(339+138+132+173)/4=196

各群の平均と全体の平均との差を計算

これを各群について計算し、それらを合計または平均することで群間のバラつきを計算します。

群内（Within-groups）

群内バラつきは、各群内でのデータがその群の平均値からどれくらいバラついているかを示します。ただし、この特定のケースでは、各群には単一のデータ点しかない（例：A1B1の群には「339」という一つのデータ点しかない）ため、群内バラつきはゼロです。

まとめ

この特定のケースにおいては、群間バラつきはありますが、群内バラつきはありません（各群には1つのデータ点しかないため）。このような状況では、一般的なANOVAのアプローチが適用できない場合もあります。通常は各群に複数のデータ点が必要です。

Pythonによる計算

AとBの因子を1,2のダミー変数で表します。

y	A	B
339	1	1
132	1	2
138	2	1
173	2	2

import pandas as pd
import itertools
import numpy as np
# 因子の水準を定義
levels = [1,2]

# すべての水準の組み合わせを生成
dumy=[(1,1),(1,2),(2,1),(2,2)]

df = pd.DataFrame(dumy, columns=['A', 'B'])
print("ダミー変数 \n",df)
#交互作用の追加
df['AB'] = df['A'] * df['B']
df= df.replace({4: 1})
design_matrix=df.values
print("デザイン行列\n",design_matrix)

response=np.array([339,132,138,173])
# 水準
levels = [1, 2]

# 主効果と交互作用効果の組み合わせ
effects = ['A', 'B', 'AB']

# 群間平方和と群内平方和の初期化
ss_between = {}
ss_within = {}

# 各効果ごとに群間平方和と群内平方和を計算
print('効果　群間平方和　水準平均　全体平均')
for effect in effects:
    # グループごとの応答変数の抽出
    tmp=0
    for lvl in levels:
        group_responses = response[design_matrix[:,effects.index(effect)]==lvl]
        # グループの水準ごとの平均
        group_mean = np.mean(group_responses)
        # 群間平方和
        tmp0=len(levels) * (group_mean - np.mean(response)) ** 2
        tmp+=tmp0
        print("  ",effect,lvl,"    ",tmp0,group_mean,np.mean(response))
    ss_between[effect] = tmp

    # 群間平方和の合計を計算
T = sum(ss_between.values())

# 結果の表示
# 全体の平方和を計算
S_T = np.square(response - response.mean()).sum()

# 結果の表示
for effect in effects:
    print(effect + ":")
    print("SS_between =", ss_between[effect])

# 群間平方和の合計と全体の平方和が一致することを確認
print("T (SS_between total) =", T)
print("S_T (Total SS) =", S_T)
print("T == S_T?", T == S_T)

全体平方和と群間平方和は異なるので、その理由としては以下のようなものが考えられます。

交互作用の存在: 因子Aと因子Bの間に交互作用が存在する場合、群間平方和の合計と全体の平方和が一致しないことがあります。この場合、交互作用項をモデルに加えることで、この差異を説明することができます。

誤差の影響: 実験の実施中に生じるランダムな誤差や外的な要因が影響している場合、群間平方和の合計と全体の平方和が一致しないことがあります。

モデルの不適切さ: 使用しているモデルがデータを適切に説明できていない場合、群間平方和の合計と全体の平方和が一致しないことがあります。

水準の不適切さ: 実験の水準が不適切であったり、水準の選択が限定的であったりすると、群間平方和の合計と全体の平方和が一致しないことがあります。

交互作用があるかどうかを判断するためには、交互作用項を含むモデルを使ってANOVAを実施するか、または交互作用のプロットを作成して視覚的に評価することが重要です。これによって、交互作用が実際に存在するのか、それとも他の要因が影響しているのかを判断することができます。

交絡

交絡とは、ある変数が他の一つ以上の変数と関連しているため、特定の変数間の真の関係が偽の関係に見えることを指します。たとえば、アイスクリームの売上と溺死の件数の関係を考えます。これらの変数は高い相関を示すことがありますが、これがアイスクリームの売上が溺死を引き起こすという意味ではありません。実際には、気温が高くなるとアイスクリームの売上が増え、同時にプールやビーチでの遊泳が増えるため、溺死の件数も増えます。この場合、気温は交絡因子となります。

この例のように、交絡は特定の変数間の関係を誤解させる可能性があるため、統計的解析や研究の設計段階で交絡の影響を考慮することが重要です。交絡の影響を取り除くための一般的な方法としては、ランダム化、層別化などがあります。

交絡(wikiの簡易翻訳 by DeepL)

統計学において、交絡因子（交絡変数、交絡因子、外来決定因子、潜伏変数とも）とは、従属変数と独立変数の両方に影響を与え、偽の関連性を引き起こす変数のことです。交絡は因果関係の概念であり、そのため相関関係や関連性で説明することはできません。交絡因子の存在は、相関が因果を意味しない理由を定量的に説明する重要なものです。交絡因子は内部妥当性を脅かすものです。

メディエーターが因果関係の連鎖の要因となるのに対し（上）、交絡因子は因果関係を誤って示唆する偽因子です（下）。

定義

交絡は、データ生成モデルの観点から定義されます。Xをある独立変数とし、Yをある従属変数とします。XのYに対する効果を推定するために、統計学者はXとYの両方に影響を与える無関係な変数の効果を抑制しなければなりません。ZがXとYの両方に因果的に影響するとき、XとYは他の変数Zによって交絡されると言いいます。

仮想的な介入X = xの下での事象Y = yの確率を$P(y \mid {\text{do}}(x))$ とします。
XとYが交絡しないのは、以下が成立する場合のみです：

$$P(y \mid {\text{do}}(x))=P(y \mid x)　(1)$$

は、X = x と Y = y のすべての値について、X = x を見たときの条件付き確率です。直感的には、この等式は、XとYの間に見られる関連が、xをランダムにした対照実験で測定される関連と同じであれば、XとYは交絡されないことを示します。

この文脈でいうと、式（１）が成立するということは、ランダム化された対照実験で測定されたXの効果$P(y∣do(x))$が、観測データで測定されたXとYの関連
$P(y∣x)$と一致するということです。

つまり、この等式が成立する場合、交絡因子がない（または効果的に制御されている）と判断できるため、観測データだけでXとYの間の因果関係を正確に推定できる、ということになります。

これは理想的な状況であり、実際のデータ解析でこのような状況が成立するかどうかは、様々な要因に依存します。しかし、この概念は因果推論の理論的基盤の一つであり、交絡が問題になる状況を理解する上で非常に有用です。

ある条件下で$P(y∣do(x))$（介入xが行われた後のyの確率）と$P(y∣x)（xが観測された場合のyの確率）が等しくなる場合、そのモデルが交絡（confounding）の影響を受けていない、すなわち因果関係をきちんと推定できる、ということを示しています。

この式が成り立つかどうかを検証するには、データ生成モデル（つまり、変数がどのように生成され、相互にどのように影響し合うのかを説明する数学的モデル）が完全に指定され、その上で関連する確率分布が知られている必要があります。

具体的には、変数間の関係を示す方程式（例えば、線形回帰モデルの方程式など）と、各変数が従う確率分布（例えば、正規分布、二項分布など）が分かっている場合、理論的にはこの式が成り立つかどうかを計算やシミュレーションで検証できます。

このような検証は、多くの場合、因果グラフや構造方程式モデリング、ベイジアンネットワークなどを用いて行われます。これらの手法は、変数間の因果関係を明示的にモデリングすることで、交絡因子の影響を評価または除去する能力を持っています。

簡単に言えば、この文は「全ての関連情報（方程式と確率分布）が揃っていれば、交絡の有無を数学的に検証できる」と言っているのです。

しかし、等式$P(y\mid {\text{do}}(x))=P(y \mid x)$を検証するには、グラフ構造だけで十分であることがわかっています。

直交表による実験計画

まずL８直交表を作ってみよう。そしてAとBとCの主効果からなる応答を作って、分散分析を行ってみよう。

import statsmodels.api as sm
from statsmodels.formula.api import ols

np.set_printoptions(precision=4, suppress=True)
# 因子の水準を定義
levels = [1,2]

# すべての水準の組み合わせを生成
combinations = list(itertools.product(levels, repeat=3))

# 直交表の作成
df = pd.DataFrame(combinations, columns=['A', 'B', 'C'])
df['AB'] = df['A'] * df['B']
df['AC'] = df['A'] * df['C']
df['BC'] = df['B'] * df['C']
df['ABC'] = df['A'] * df['B'] * df['C']
df = df.replace({4: 1, 8: 2})
y=[np.sum(a[:3]) for a in df.values]#割付
data = {"factorA": df.A,
        "factorB": df.B,
        "factorC": df.C,
        "factorAB": df.AB,
        "factorAC": df.AC,
        "factorBC": df.BC,
        "factorABC": df.ABC,        
        "y": y}
df2 = pd.DataFrame(data)
print(df2)
model = ols('y ~ C(factorA) + C(factorB) + C(factorC) + C(factorAB)+ C(factorAC)+ C(factorBC) ', data=df2).fit()
anova_table = sm.stats.anova_lm(model)
print(anova_table)

有意であるのが、主効果A,B、Cであるのが分かります。

つぎに割付を考えてみましょう。もし因子CをfactorABで割り付けてしまったらどうなるでしょうか。
やってみましょう。

y=[np.sum(a[0]+a[1]+a[3]) for a in df.values]#割付
data = {"factorA": df.A,
        "factorB": df.B,
        "factorC": df.C,
        "factorAB": df.AB,
        "factorAC": df.AC,
        "factorBC": df.BC,
        "factorABC": df.ABC,        
        "y": y}
df2 = pd.DataFrame(data)
print(df2)
model = ols('y ~ C(factorA) + C(factorB) + C(factorC) + C(factorAB)+ C(factorAC)+ C(factorBC) ', data=df2).fit()
anova_table = sm.stats.anova_lm(model)
print(anova_table)

主効果Cが有意ではなく、相互作用ABが有意になりました。つぎにすべての効果を応答に加えてみましょう。

y=[np.sum(a) for a in df.values]
data = {"factorA": df.A,
        "factorB": df.B,
        "factorC": df.C,
        "factorAB": df.AB,
        "factorAC": df.AC,
        "factorBC": df.BC,
        "factorABC": df.ABC,        
        "y": y}
df2 = pd.DataFrame(data)
print(df2)
model = ols('y ~ C(factorA) + C(factorB) + C(factorC) + C(factorAB)+ C(factorAC)+ C(factorBC) ', data=df2).fit()
anova_table = sm.stats.anova_lm(model)
print(anova_table)

どの効果も有意でなくなってしまいました。これは全ての効果から生成した応答なので、へんですが、残差の平方和が他の平方和と同じ大きさになってしまっているのが原因です。つまり、結果の解釈には注意が必要だということです。

統計学実践ワークブック20章問20.5 p.177

問題は直交表とは関係ありませんが、直交表を作って、平方和を求めてみます。

タイヤの製造について考える。A1とA2のゴムの種類とB1,B2のタイヤのトレッドパターンがあり、これを(A1,B1),(A1,B2),(A2,B1),(A2,B2)の組み合わせで５組のタイヤを製造する。それを５台の四輪の車両V1,V2,V3,V4,V5に１本ずつランダムに決めた箇所に装着し、一定期間走行させた。 その耐久性は以下のとおりである。

V1	A1B1	132	V2	A1B1	100	V3	A1B1	95
	A1B2	116		A1B2	86		A1B2	84
	A2B1	112		A2B1	87		A2B1	78
	A2B2	117		A2B2	86		A2B2	84
V4	A1B1	114	V5	A1B1	110
	A1B2	103		A1B2	97
	A2B1	100		A2B1	94
	A2B2	102		A2B2	97

SA=320, SB=125,SV=2862.2
SAxB=320, SAxV=5.0, SBxV=13.00
SAxBXV=11.0, ST=3656.2

# 因子の水準を定義
levels_V = [1, 2,3,4,5]
levels_A = [1, 2]
levels_B = [1, 2]

# すべての水準の組み合わせを生成
combinations = list(itertools.product(levels_V, levels_A, levels_B))

# 直交表の作成
df4 = pd.DataFrame(combinations, columns=['V', 'A', 'B'])
df4['response']=[132,116,112,117,100,86,87,86,95,84,78,84,114,103,100,102,110\
,97,94,97]
df4

つぎに回帰分析を行いますが、ABVについては指定していません。これを指定しまうと残差がゼロになりうまく動きません。この場合にはABVは残差になります。

model = ols('response ~ C(A) + C(B) + C(V) \
        + C(A)*C(B) +C(A)*C(V) + C(B)*C(V) ', data=df4).fit()

anova_table = sm.stats.anova_lm(model,typ=2) 
print(anova_table)

他モデルとの比較

他のモデルと比較することで分散分析の適用範囲がより明確になることがあります。

線形回帰モデル

線形回帰モデルは、応答変数と説明変数の間の関係を線形的にモデル化する手法です。通常、最小二乗法を使用してモデルパラメータを推定し、回帰係数や切片などのパラメータを求めます。線形回帰モデルは、連続型の応答変数を予測するために広く使用されています。

一方、分散分析（ANOVA）は、グループ間の平均値の差を評価する統計的手法です。分散分析では、要因とその水準に基づいてデータをグループ分けし、それぞれのグループ間の差異を検定します。分散分析は、カテゴリカルな要因に関連する応答変数を分析するために広く使用されています。

線形回帰モデルと分散分析の共通点は、両方の手法がデータの変動を説明し、要因や説明変数の効果を評価することにあります。ただし、線形回帰モデルは連続型の応答変数を予測するために使用され、説明変数の効果を数値化することが一般的です。一方、分散分析はカテゴリカルな要因に関連する応答変数を分析し、要因間の差異を評価することが主な目的です。

固定効果モデル

固定効果モデルは、カテゴリカルな要因（または連続的な予測子）の水準間の差異を推定するための統計モデルです。データを複数のグループに分け、各グループにおける平均値の差を評価します。具体的には、各水準の効果をパラメータとしてモデル化し、最尤推定などの手法を用いてパラメータを推定します。固定効果モデルは、カテゴリカルな要因に基づく効果を評価するために使用され、分散分析などの手法と結びつけて使用されることがあります。

一方、分散分析は、カテゴリカルな要因に基づく効果を評価するための統計手法です。特に、一つの応答変数を複数のグループに分けて、グループ間の差異を検定します。分散分析では、グループ間の平均値の差異を推定するために、群間平方和と群内平方和を計算します。そして、これらの平方和を用いてF検定などの統計検定を行い、要因の効果の有意性を評価します。

要点としては、固定効果モデルは効果をパラメータとしてモデル化し、パラメータの推定に焦点を当てています。一方、分散分析は群間平方和と群内平方和を用いて効果の有意性を評価するための手法です。固定効果モデルは個々の効果の推定や統計的な検定に重点を置いており、分散分析はグループ間の差異を検定するための手法としてよく知られています。

固定効果モデルと分散分析はしばしば混同されることがありますが、厳密には異なる手法です。固定効果モデルはより一般的で柔軟なモデリング手法であり、分散分析は特定のデザインや問題に特化した手法と考えることができます。

ランダム効果モデル

ランダム効果モデルは、データの背後にある階層的な構造や個体差を考慮してモデル化する統計モデルです。データが複数のグループに分けられている場合、各グループにおける個体間の差異をランダム効果として扱います。具体的には、ランダム効果を確率分布に従うパラメータとしてモデル化し、最尤推定やベイズ推定などの手法を用いてパラメータを推定します。ランダム効果モデルは、個体差やグループ間の異質性を考慮するため、データのばらつきをより正確にモデル化することができます。

要点としては、ランダム効果モデルは個体差やグループ間の異質性を考慮してデータをモデル化し、ランダム効果の推定や効果の統計的な検定に焦点を当てています。一方、分散分析は群間平方和と群内平方和を用いて効果の有意性を評価するための手法です。

ランダム効果モデルと分散分析はしばしば混同されることがありますが、厳密には異なる手法です。ランダム効果モデルは階層モデリングや個体差のモデリングにより適しており、分散分析は特定のデザインや問題に特化した手法と考えることができます。

混合効果モデル

混合効果モデルは、個体間の異質性や階層構造を考慮した統計モデルです。データセットが階層構造を持つ場合や、観測値が相互に依存する場合に適用されます。混合効果モデルは、固定効果とランダム効果の両方を含んでおり、個体間の異質性やクラスター間の差異を推定することができます。また、混合効果モデルは、個体間の相関構造をモデル化することができるため、観測値の相互作用や重複を考慮して分析することができます。

要約すると、分散分析は主にグループ間の差異を評価するために使用され、個々の観測値に対する効果の推定は行いません。一方、混合効果モデルは、個体間の異質性や階層構造を考慮して、個体間の差異や相関構造を推定することができます。

多重共線性、交絡因子、交互作用

交絡因子と交互作用

交絡因子（Confounders）

交絡因子は研究で主要な関心を持つ独立変数と従属変数の両方に影響を与える可能性がある第三の変数です。これがあると、独立変数と従属変数の関係は偽りのものかもしれません。つまり、実際には第三の変数が影響を与えている可能性があるのです。交絡因子が存在する場合、それを制御しないと、誤った因果関係の結論に至る可能性があります。

交互作用（Interaction）

交互作用は、二つ以上の独立変数が従属変数に与える影響が、単純な合成作用を超えている場合に発生します。言い換えれば、一つの独立変数の影響が、もう一つの独立変数のレベルによって変わる場合、交互作用があると言います。

多重共線性との関係

多重共線性（Multicollinearity）は統計学と特に回帰分析において使用される用語で、この概念は交絡因子（Confounders）や交互作用（Interaction）とは異なります。

多重共線性とは

多重共線性とは、回帰モデルにおいて、独立変数（説明変数）間に高度な相関関係が存在する状態を指します。多重共線性がある場合、それぞれの独立変数が従属変数（目的変数）に与える影響を正確に評価することが困難になります。さらに、係数の推定が不安定になる場合もあります。

交絡因子との違い

交絡因子は、独立変数と従属変数の両方に影響を与える第三の変数です。多重共線性と交絡は異なる問題であり、交絡因子がいくつかの独立変数に影響を与える場合でも、それが自動的に多重共線性を引き起こすわけではありません。

交互作用との違い

交互作用は、二つ以上の独立変数の影響が非線形に組み合わさることを指します。多重共線性は独立変数間の線形関係に関する問題であり、交互作用の存在が多重共線性を引き起こすわけではありません。

総括

多重共線性は主に独立変数間の相関に起因する問題であり、交絡因子と交互作用はそれぞれ異なる統計的な問題を解決するための概念です。しかし、これらは全てモデリングの際に考慮すべき重要な側面であり、特に複雑な研究デザインやデータセットでは、これらの問題が同時に発生する可能性もあります。

線形回帰と交絡因子

線形回帰において交絡因子が期待される場合、その交絡因子を適切に制御することが重要です。以下は一般的なアプローチです。

交絡因子をモデルに組み込む

交絡因子を回帰モデルに追加することで、その影響を制御します。これにより、目的変数に対する他の説明変数の「純粋な」効果を推定できるようになります。

$Y=f(X_i)+e=\beta_0+\beta_1 X_1 + \beta_2 X_2 + e $

$X_2$が$X_1$と$Y$に影響を与える交絡因子であると仮定します。

マッチングまたはプロペンシティスコア

交絡因子が多い場合や、交絡因子が連続変数である場合には、マッチングやプロペンシティスコアの手法を用いることで、交絡因子を制御することができます。

層別解析

交絡因子ごとに層を作り、それぞれで独立した回帰分析を行う方法もあります。しかし、この方法は解釈が複雑になりやすく、またデータの利用効率が低下する可能性があります。

交互作用項を加える

交絡因子が主効果だけでなく交互作用も持っている可能性がある場合、交互作用項もモデルに加えると良いでしょう。

$Y=f(X_i)+e=\beta_0+\beta_1 X_1 + \beta_2 X_2 + \beta_3(X_1+X_2)+e$

ここで、$\beta_3(X_1+X_2)$が相互作用項です。

階層型モデル

相互作用が多数ある場合や、階層型のデータ構造がある場合は、階層型モデルを考慮することもあります。

線形混合モデルまたは一般化線形混合モデル

ランダム効果項が必要な場合（例えば、多階層データや繰り返し測定データなど）、線形混合モデル（LMM）や一般化線形混合モデル（GLMM）も有用です。

モデル選択

AIC（赤池の情報量規準）、BIC（ベイズ情報量規準）などのモデル選択基準を用いて、相互作用項を含むモデルと含まないモデルとを比較することもあります。

モデル診断と検証

相互作用項を加えた後も、モデル診断（残差分析、多重共線性のチェックなど）と検証（交差検証、ブートストラップなど）は必須です。

相互作用がある場合、その効果を適切に制御することが非常に重要です。過度に複雑なモデルを避け、必要な相互作用項だけを含めるように注意が必要です。モデルの解釈可能性と予測精度のバランスを考慮しながら、最も適したモデルを選択する必要があります。

フルモデル

フルモデル（Full Model）は、統計モデリングや回帰分析において、すべての予測変数（説明変数）を含めた完全なモデルを指します。つまり、調査や実験のすべての要因や変数を考慮したモデルです。

フルモデルの意義

• 因子や変数の重要性の特定: フルモデルにはすべての予測変数が含まれているため、それぞれの因子や変数の結果への寄与を評価できます。重要な要因や変数を特定することで、現象やプロセスの理解を深めることができます。

• モデルの解釈性: フルモデルを使用すると、異なる予測変数の間の相互作用や非線形関係など、より複雑なモデル構造を考慮できます。これにより、実世界の現象や関係をより正確にモデル化することができます。

ただし、フルモデルは過剰適合のリスクがあるため、データセットのサイズやモデルの目的に応じて適切なバランスを見つける必要があります。過剰適合を避けるためには、交差検証や情報基準などの手法を使用してモデルの適切な複雑さを決定することが重要です。

交絡因子とフルモデル

フルモデルを用いたシミュレーションは交絡効果（confounding effects）の理解に有用です。特に、実世界のデータでは観測できないような条件や設定を仮想的に作成して、その影響を分析することができます。以下は、シミュレーションを交絡効果の理解にどう生かすかのいくつかのアイデアです：

効果の分解: フルモデルを用いてシミュレーションを行うことで、主効果、交互作用、交絡効果を個別に評価することができます。これにより、どの要因が結果にどれだけ影響を与えるのか、明確にすることができます。

仮説検証: シミュレーションを通じて、特定の交絡因子がどれだけ影響を与えるのかの仮説を検証することができます。これは、後の実世界での実験設計やデータ分析に有用な知見を提供します。

実験設計の改善: シミュレーションで得られた知見を用いて、実験設計をより効率的、または有効なものに改善することができます。

リスク評価: 交絡効果が強いと判明した場合、それがもたらす可能性のあるリスクを事前に評価することができます。

制御戦略の評価: 異なる種類の統計的制御戦略（例：マッチング、傾向スコア、インバリアントセットなど）をシミュレーション内で試し、どれが最も効果的に交絡を制御できるのかを評価できます。

教育とコミュニケーション: シミュレーションの結果は、交絡効果をより広く理解するための教材や説明資料として用いることができます。

パラメータの感度分析: フルモデルに基づくシミュレーションは、モデルのパラメータが結果に与える影響を理解するのにも役立ちます。これは、モデルの信頼性を評価する一つの方法です。

政策/治療効果の予測: 交絡効果を正確に捉えた上でシミュレーションを行えば、異なる政策や治療がもたらす効果をより正確に予測することが可能です。

シミュレーションは理論と実践のギャップを埋める強力なツールですが、その設計と解釈には慎重さが求められます。特に、シミュレーションの結果が現実世界の複雑性を適切に反映しているかどうかを常に評価する必要があります。

交互作用の理解

フルモデルは交互作用の意味を理解するために役立ちます。交互作用は、予測変数同士が相互に影響し合うことを意味します。つまり、予測変数の組み合わせによって、結果に与える影響が変化する場合に交互作用が存在します。

フルモデルでは、すべての予測変数やその組み合わせを含めて分析するため、交互作用の効果を明確に評価することができます。予測変数同士の相互作用が結果にどのような影響を与えるのかを分析することで、現象やプロセスの理解が深まります。

残差と交互作用

フルモデルを用いると、残差と交互作用の区別が一部難しくなる場合があります。これは、フルモデルではすべての予測変数やその組み合わせが含まれており、予測変数同士の相互作用も考慮されるためです。

フルモデルにおける残差は、モデルが観測されたデータに対してどれだけ適合していないかを表す指標です。残差は、観測された応答変数とモデルによって予測された応答変数との差を表します。ただし、フルモデルでは予測変数同士の相互作用も考慮されているため、残差には交互作用の影響も含まれています。

したがって、フルモデルを用いる場合、交互作用の影響と残差の区別を明確にすることは難しくなるかもしれません。ただし、統計的な手法やグラフィカルな手法を用いることで、残差と交互作用の影響をより細かく解釈することができます。また、統計モデルの仮定やモデルの評価指標なども考慮しながら解釈を行うことが重要です。

2因子、2水準、応答=[339, 132, 138, 173]の例

この例では主効果をAとB、交互作用をABとして平方和を計算しました。しかし、総平方和から主効果の平方和を差し引いたものが、残差になるか交互作用によるものかは、現象そのものがなんであるかがよくわからなければ、判断できません。

繰り返しのない実験で、応答を均衡値と考えるとそれらは平均値です。そこで、フルモデルを頭に描くとそのモデルの本質が理解しやすくなります。また、実際の実験では繰り返し観測を行いますが、その数が少なければ、どのような実験に直交表が役に立つのだろうかという点も理解しやすくなります。

statsmodelsの利用

statsmodelsではANOVA専用のライブラリがあります。線形OLSModelを用いたANOVA分析のためのanova_lmと、反復測定ANOVA用のAnovaRMを含む分散分析モデルです。

このコードはols(最小二乗法:数式形式)を使用して線形回帰モデルを推定します。

例

levels=[1,2]
combinations=list(itertools.product(levels,repeat=2))

df = pd.DataFrame(dumy, columns=['A', 'B'])
print("ダミー変数 \n",df)
#交互作用の追加
df['AB'] = df['A'] * df['B']
df= df.replace({4: 1})
design_matrix=df.values
print("デザイン行列\n",design_matrix)

response=np.array([339,132,138,173])

# 主効果と交互作用効果の組み合わせ
effects = ['A', 'B', 'AB']
data = {"A": df.A,
        "B": df.B,
        "AB": df.AB,
        "response": response}
df2 = pd.DataFrame(data)
print(df2)

model = ols('response ~ C(A) + C(B) + C(AB) ', data=df2).fit()

anova_table = sm.stats.anova_lm(model)
print(anova_table)

total_ss = np.square(response - response.mean()).sum()
print("Total SS =", total_ss)
print("ANOVA table sum of squares =", anova_table.sum_sq.sum())
print("Total SS == ANOVA table sum of squares?", total_ss ==anova_table.sum_sq.sum())

data = {"A": df.A,
        "B": df.B,
        "response": response}
df2 = pd.DataFrame(data)
print(df2)
model = sm.OLS(df2.response,df2.iloc[:,:-1]).fit()
model = ols('response ~ C(A) + C(B) ', data=df2).fit()

anova_table = sm.stats.anova_lm(model)
print(anova_table)
np.var(response)*4,anova_table.sum_sq.sum()

olsのRスタイルの公式

C(A)、C(B)、およびC(AB)は独立変数としてモデルに含められています。C(A)、C(B)、C(AB)はカテゴリ変数を示しており、これらの因子がカテゴリ変数であるため、C()関数はカテゴリ変数をダミー変数（指示変数）に変換します。ダミー変数化により、カテゴリ変数の効果を線形モデルで評価できるようになります。

C(A)は因子Aのダミー変数化を表します。
C(B)は因子Bのダミー変数化を表します。
C(AB)は交互作用項AB（AとBの組み合わせ）のダミー変数化を表します。

0,1の水準

直交表で水準を0と1にするメリットとデメリットは次のようになります。

メリット：

簡潔な表現：水準を0と1に限定することで、直交表のセルには0か1の値が入るため、表が簡潔で見やすくなります。
解釈の容易さ：0と1の水準は対比関係を表現しており、因子の効果の解釈が直感的になります。
計算の容易さ：0と1は二進数で表現されるため、計算や統計的な処理が簡単になります。
デメリット：

水準の制約：水準を0と1に限定することで、実際の要因の水準範囲を表現できない場合があります。そのため、実験の条件や研究対象に応じて、より適切な水準を選ぶ必要があります。
情報の欠落：水準が0と1のみの場合、因子の効果の詳細な情報や非線形の関係を捉えることができない場合があります。
応用の制約：一部の実験や研究において、水準を0と1に制限することが適切でない場合があります。その場合は、より多様な水準を持つ直交表や他の実験計画法を使用する必要があります。

水準を0と1にするかどうかは、実験や研究の目的や制約、解析方法などによって異なります。適切な水準設定を選ぶことが重要です。

OLSの使用

次のコードは、水準が0と1の2水準因子A、B、および交互作用項ABに対する線形回帰モデルを構築するためのものです。responseとして、指定された値が入力され、sm.OLSを使用して、因子A、B、ABと応答変数responseとの関係を線形回帰モデルで推定します。このコードはstatsmodelsライブラリを使用しており、OLS(最小二乗法)を使用して線形回帰モデルを推定します。この際には、モデルの要約には、推定された係数、標準誤差、t値、p値、決定係数（R-squared）などの統計的な情報が含まれます。

levels = [0,1]

# すべての水準の組み合わせを生成
combinations = list(itertools.product(levels, repeat=2))

# 直交表の作成
df = pd.DataFrame(combinations, columns=['A', 'B'])
response = np.array([339, 132, 138, 173])
data = {"A": df.A,
        "B": df.B,
        "response": response}
df2 = pd.DataFrame(data)
print(df2)
model = sm.OLS(df2.response,df2.iloc[:,:-1]).fit()
model = ols('response ~ C(A) + C(B) ', data=df2).fit()


anova_table = sm.stats.anova_lm(model)
print(anova_table)
np.var(response)*4,anova_table.sum_sq.sum()

levels = [0,1]

# すべての水準の組み合わせを生成
combinations = list(itertools.product(levels, repeat=2))

# 直交表の作成
df = pd.DataFrame(combinations, columns=['A', 'B'])
df['AB'] = df['A'] * df['B']
response = np.array([339, 132, 138, 173])
#response = np.array([0, 1, 1, 3])
data = {"A": df.A,
        "B": df.B,
        "AB": df.AB,
        "response": response}
df2 = pd.DataFrame(data)
print(df2)
#model = sm.OLS(df2.response,df2.iloc[:,:-1]).fit()# 定数項のあるなしは適宜判断
model = sm.OLS(df2.response, sm.add_constant(df2.iloc[:, :-2])).fit()
model.summary()

直交表の行の和とresponse（モデルの作用)

要因の水準の和を取ることで、異なる要因や水準が応答変数（response）に与える効果を調査できます：

効果の定量化：要因の水準を合計することで、それぞれの要因が応答変数にどれだけの効果を持つかを定量的に評価できます。これにより、要因の重要性を比較し、最も影響の大きい要因を特定できます。

相互作用の検出：複数の要因の相互作用がある場合、単一の要因だけでなく、その相互作用も考慮する必要があります。要因の水準の和を取ることで、相互作用が応答変数にどれだけの効果を持つかを評価できます。

効果の分解：要因の水準の和を取ることで、応答変数の変動を異なる要因や相互作用による効果に分解できます。これにより、どの要因が応答変数に最も影響を与えるかを把握し、モデルの解釈を容易にします。

要因の水準の和を取ることで、実験デザインや研究の解釈を支援し、要因間の関係や応答変数への効果をより深く理解することができます。これは、効果的な実験計画や解析を行うための重要な手法です。

response = np.array([0, 1, 1, 3])
data = {"A": df.A,
        "B": df.B,
        "AB": df.AB,
        "response": response}
df2 = pd.DataFrame(data)
print(df2)
model = sm.OLS(df2.response,df2.iloc[:,:-1]).fit()
model.summary()

水準が[1,2]の場合と[0,1]の場合で、responseに要因の和を使用した場合のOLSモデルの結果は異なります。具体的には、水準が[1,2]の場合、因子AとBの値は1または2になります。水準が[0,1]の場合、因子AとBの値は0または1になります。

要因の水準が変わることで、デザイン行列（要因の値を含む行列）が異なり、結果的にOLSモデルの推定係数が変わります。しかし、どちらの水準設定でも、交互作用の有無や主効果の方向性についての結論は変わらない可能性があります。

繰り返し実験

import pandas as pd
import itertools

# 因子の水準を定義
levels = [1, 2]
repeats = 9

# すべての水準の組み合わせを生成
combinations = list(itertools.product(levels, repeat=2))
# 直交表の作成
df = pd.DataFrame(combinations*repeats, columns=['A', 'B'])
# 結果の表示
df['AB'] = df['A'] * df['B']
df= df.replace({4: 1})
design_matrix=df.values
response=np.array([339,132,138,173,190,59,197,157,423,381,368,235,419\
                   ,178,142,192,201,77,126,168,384,283,383,230,289,202,206,166,120,25,204,194,312,247,345,171])
data = {"A": df.A,
        "B": df.B,
        "AB": df.AB,
        "response": response}
df2 = pd.DataFrame(data)
print(df2)
model = sm.OLS(df2.response,df2.iloc[:,:-1]).fit()
model = ols('response ~ C(A) + C(B) +C(AB) ', data=df2).fit()

anova_table = sm.stats.anova_lm(model)
print(anova_table)