書記が数学やるだけ#336 二元配置分散分析の概要
分散分析についても，説明変数を増やすことができる。
問題
本問は『一般線形モデルによる生物科学のための現代統計学』で公開されているデータセットに基づいている。
説明
二元配置分散分析では，2つの因子からなるデータを分析する。
重回帰分析と同様，第3の変数を導入することにより実際の効果が変わることがある。
詳しくは後日に回すが，因子同士が直交していると仮定することで，直交表を用いて実験回数を減らすことができる。
解答
今回もPythonを用いる。
#ライブラリの読み込み
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy as sp
import statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf#データの読み込み
data = pd.read_csv('分散分析-2.csv')
data
data.dtypesWATER int64
INITHT float64
FINALHT float64
dtype: object
ここで，WATERをカテゴリカル変数にするために，文字列に変換する。
data['WATER'] = data['WATER'].astype(str)
data.dtypesWATER object
INITHT float64
FINALHT float64
dtype: object
まず，WATERのみを説明変数としてモデルを作ってみる。箱ひげ図で可視化してみると，2と3では差がないようにも見える。
#WATERによるFINALHTの違いを可視化
sns.boxplot(x=data["WATER"], y=data["FINALHT"])
sns.swarmplot(x="WATER", y="FINALHT",data=data,dodge=True)
分散分析表より，有意水準5%で有意でない，つまりWATERはFINALHTに影響を与えるとは言えない。
#モデルFINALHT ~ WATER
anova_mod1 = smf.ols('FINALHT ~ WATER', data=data).fit()
#分散分析表
anova_table1 = sm.stats.anova_lm(anova_mod1)
anova_table1
ところが，ここで第3の変数としてINITHTを説明変数に加えるとどうなるか。
参考までに，INITHTのみを説明変数としたモデルを評価してみると，INITHTはFINALHTに影響を与えると言えそうである。
#モデル化FINALHT ~ INITHT
anova_mod2 = smf.ols('FINALHT ~ INITHT', data=data).fit()
#分散分析表
anova_table2 = sm.stats.anova_lm(anova_mod2)
anova_table2
では，WATERとINITHTを説明変数としたモデルを作る。すると，WATERとINITHTの両方がFINALHTに影響を与えると言えそうである。
これは，見た目上は散水頻度は最終樹高に影響を与えないが，初期樹高を差し引いて考えると散水頻度は最終樹高に影響を与えることを示している，考えてみれば当然のことではある。
#モデル化FINALHT ~ WATER+INITHT
anova_mod3 = smf.ols('FINALHT ~ WATER+INITHT', data=data).fit()
#分散分析表
anova_table3 = sm.stats.anova_lm(anova_mod3)
anova_table3
ここで，さらに見通しを良くするために「FINALHT-INITHT=delta（樹高の変化）」として，WATERを説明変数にdeltaに影響を与えるか考えてみる。
可視化してみると，WATERごとに違いがありそうである。
#WATERによるdeltaの違いを可視化
data['delta'] = data['FINALHT']-data['INITHT']
sns.boxplot(x=data["WATER"], y=data["delta"])
sns.swarmplot(x="WATER", y="delta",data=data,dodge=True)
実際に，WATERとはdeltaに影響を与えると言えそうである。
#モデルdelta ~ WATER
anova_mod4 = smf.ols('delta ~ WATER', data=data).fit()
#分散分析表
anova_table4 = sm.stats.anova_lm(anova_mod4)
anova_table4
本記事のもくじはこちら：
