本項でもGoogle Colaboratoryを使用する．（もちろんファイル読み込みの部分以外はJupyter Notebookでもよいし，一気にPythonのソースコードをまとめて書いて実行しても問題はない）
異なる2地点以上の比較
前回のt検定は，2つの群の間で，統計的に本当に違いがあるかをいうためのものだった．，
今回は，1977年〜2018年の，青梅，八王子，府中の3箇所，つまり3群で違いがあるかを見る．
前回と全く同じ気象庁のオープンデータを使う．
分散分析(一元配置分散分析)とは
今回の場合は「3つ以上の地点の気温データ群の間の差を，統計的に（ある程度）保証してくれる」のが分散分析である．要するに3群以上の検定．
分散分析にもt検定と同じく，対応がある場合と対応がない場合があるので，前回を復習しておくこと．今回は対応がない場合のみを取り扱う．
Bartlett検定
t検定では，2群の分散（散らばり方）が等しい時とそうでない時で，適用する方法が違っていた．（Studentのt検定と，Welchのt検定）3群以上の場合に等分散性をみるのがBartlett検定である．
F検定と同じく，帰無仮説は「与えられた3群（以上）の分散は等質である」，対立仮説は「与えられた3群（以上）の分散は等質であるとは言えない」である．したがって結果のp値の見方も同じく，0.05以上であれば帰無仮説が棄却されず「分散が糖質である」となる．
分散分析のコード例
とはいっても，データの整え方などはt検定の時と同じである．
1つ目のセル．2回実行してCSVファイルを2つ読み込んでおく．t検定の時と全く同じ
from google.colab import files
uploaded = files.upload()2つ目のセル．これもt検定の時と前回と全く同じ．
import numpy as np
import pandas as pd
from pandas import Series, DataFrame
# 1月の月平均気温を読み込む
df_jan = pd.read_csv('tokyo-temp-data-multi-point-jan.csv', encoding='SHIFT-JIS', header=1)
# 8月の月平均気温を読み込む
df_aug = pd.read_csv('tokyo-temp-data-multi-point-aug.csv', encoding='SHIFT-JIS', header=1)
# headerの行数1で「ダウンロードした日時」を削除
# 確認のために，一度全体をprint
print(df_jan)
print(df_aug)
# printした内容から，余計なものを削除していく
# 品質情報とかの行を削除
df_jan = df_jan.drop(1)
# 同じ処理を8月に対しても適用
df_aug = df_aug.drop(1)
# 確認
print(df_jan)
# 「年月」「平均気温」．．．．．となっている行を削除,
df_jan = df_jan.drop(0)
df_aug = df_aug.drop(0)
#再度確認
print(df_jan)
# 品質情報とかが入っていた列を削除する．
# 列を削除するときには，columns=[削除する列名の配列]とする
df_jan = df_jan.drop(columns=['江戸川臨海.1', '江戸川臨海.2', '東京.1', '東京.2', '青梅.1', '青梅.2', '八王子.1', '八王子.2', '練馬.1', '練馬.2', '府中.1', '府中.2'])
df_aug = df_aug.drop(columns=['江戸川臨海.1', '江戸川臨海.2', '東京.1', '東京.2', '青梅.1', '青梅.2', '八王子.1', '八王子.2', '練馬.1', '練馬.2', '府中.1', '府中.2'])
#次元を確認
print(df_jan.values.shape)
#最終確認
df_jan.head()3つ目のセル．今度は青梅，八王子，府中にだけにして，それぞれ1次元配列で取り出す．
最後は関数チェーンを使って短く書いている．3回やらなければならないので，1回目はちゃんと確認しながら進め，2回目3回目は同じ手順なので，省略して書いている．
# 多摩の内陸部だけを比較したい
# 青梅，八王子，府中だけにする
df_jan_tama = df_jan.drop(columns=['東京', '江戸川臨海', '練馬'])
df_aug_tama = df_aug.drop(columns=['東京', '江戸川臨海', '練馬'])
print(df_jan_tama)
print(df_aug_tama)
#それぞれ列を取り出して，データの次元を確認する
ome_jan = df_jan_tama.loc[:, ['青梅']]
ome_aug = df_aug_tama.loc[:, ['青梅']]
print('Ome. shape:', ome_jan.shape)
# 1次元にする
ome_jan = ome_jan.values.flatten()
ome_aug = ome_aug.values.flatten()
# 次元を確認
print('Ome Jan. shape:', ome_jan.shape)
# 中身を確認
print(ome_jan)
# printした結果，中身が文字列なので，実数に変更
ome_jan = ome_jan.astype('float32')
ome_aug = ome_aug.astype('float32')
# 中身を確認する
print(ome_jan)
print(ome_aug)
print(ome_jan.shape)
print(ome_aug.shape)
# 今回の分散分析は項目間に関係性がないt検定の3群以上版なので，サンプル数が違っていても良い
# なので欠損地nanだけ取り除いても大丈夫
# 今回対象にしている1次元配列からnanだけ取り除くコードは以下の通り．（nanが含まれていないなら何も起こらない)
ome_jan = ome_jan[~np.isnan(aume_jan)]
#もしnanが含まれていたら，shapeが違う
print(ome_jan.shape)
# 8月についても行う
ome_aug = ome_aug[~np.isnan(aume_aug)]
# 配列名 = 配列名[~np.isnan(配列名)]
# 八王子，府中のデータも同じ様に処理する
hachioji_jan = df_jan_tama.loc[:, ['八王子']]
hachioji_aug = df_aug_tama.loc[:, ['八王子']]
hachioji_jan = hachioji_jan.values.flatten()
hachioji_aug = hachioji_aug.values.flatten()
hachioji_jan = hachioji_jan.astype('float32')
hachioji_aug = hachioji_aug.astype('float32')
hachioji_jan = hachioji_jan[~np.isnan(hachioji_jan)]
hachioji_aug = hachioji_aug[~np.isnan(hachioji_aug)]
print(hachioji_jan)
print(hachioji_jan.shape)
print(hachioji_aug)
print(hachioji_aug.shape)
# 「関数チェーン」でつなげて書くと短くなる
fuchu_jan = df_jan_tama.loc[:, ['府中']].values.flatten().astype('float32')
fuchu_jan = fuchu_jan[~np.isnan(fuchu_jan)]
fuchu_aug = df_aug_tama.loc[:, ['府中']].values.flatten().astype('float32')
fuchu_aug = fuchu_aug[~np.isnan(fuchu_aug)]
print(fuchu_jan)
print(fuchu_jan.shape)
print(fuchu_aug)
print(fuchu_aug.shape)4つ目のセル．これも前回と同じく，単純に全ての年での平均を出している．
# データが揃ったので，とりあえず平均を出してみる
print('1977年から1988年までの1月の青梅の平均気温', np.average(ome_jan))
print('1977年から1988年までの8月の青梅の平均気温', np.average(ome_aug))
print('1977年から1988年までの1月の八王子の平均気温', np.average(hachioji_jan))
print('1977年から1988年までの8月の八王子の平均気温', np.average(hachioji_aug))
print('1977年から1988年までの1月の府中の平均気温', np.average(fuchu_jan))
print('1977年から1988年までの8月の府中の平均気温', np.average(fuchu_aug))4つ目のセルの結果は以下の通り．
1977年から1988年までの1月の青梅の平均気温 2.804762
1977年から1988年までの8月の青梅の平均気温 25.483334
1977年から1988年までの1月の八王子の平均気温 3.0928574
1977年から1988年までの8月の八王子の平均気温 26.07857
1977年から1988年までの1月の府中の平均気温 4.102381
1977年から1988年までの8月の府中の平均気温 26.4571425つ目のセル．Bartlett検定を行う．
帰無仮説「3つの地点，青梅，八王子，府中の1977〜2018年の1月（8月）の気温の分散は等しい」 対立仮説「分散が等しくない」
import scipy.stats
# 2群の時のF検定に相当する等分散性の検定の3群以上版が，Bartlett検定
result_bart_jan = scipy.stats.bartlett(ome_jan, hachioji_jan, fuchu_jan)
print('1月のBartlett検定の結果', result_bart_jan)
result_bart_aug = scipy.stats.bartlett(ome_aug, hachioji_aug, fuchu_aug)
print('8月のBartlett検定の結果', result_bart_aug)5つ目のセルの結果は以下のようになる．
1月のBartlett検定の結果 BartlettResult(statistic=0.009294277545815795, pvalue=0.9953636424693709)
8月のBartlett検定の結果 BartlettResult(statistic=0.1273940963804959, pvalue=0.9382892128921961)1月の場合も8月の場合も，p値が十分に大きいので，青梅，八王子，府中の3地点に関しては帰無仮説が棄却されず，分散が等しいと言うことができる．
6つ目のセル．分散分析を行う．
帰無仮説「3つの地点，青梅，八王子，府中の1977〜2018年の1月（8月）の気温に差がない」
対立仮説「どこか1つ（以上）の地点で差がある」
である．
from scipy import stats
# 3群以上の1元配置分散分析
result_jan = stats.f_oneway(ome_jan, hachioji_jan, fuchu_jan)
print(result_jan)
result_aug = stats.f_oneway(ome_aug, hachioji_aug, fuchu_aug)
print(result_aug)6つ目のセルの結果は以下のようになる．
F_onewayResult(statistic=15.758491657818265, pvalue=8.07470716601788e-07)
F_onewayResult(statistic=5.972340163640083, pvalue=0.0033464517636488693)p値が0.05より十分に低いので，帰無仮説が棄却され，有意差があるという結果が出たが，対立仮説は「どこか1つ以上の地点で差がある」であることに注意．
statsmodelsのライブラリを使う場合，ちょっと複雑．
statsmodelsの分散分析は「値とラベルがペアになったものがたくさんある」というデータ形式で与えなければならない．ここでいう「ラベル」とは，この値がどの群のものなのか，を表す．それをDataFrameにしたものを与える．
気温の配列は全部繋げた上で，観測地点ラベルの配列も作り，それを2次元配列にまとめる．その後次元を変更し 「気温と観測地点ラベルが対応したペアが，126個存在する」という感じに変更する．
このような，値とラベルが順不同で積み重なっているデータ形式を「スタック形式」という．それに対して先に作った表のデータは「アンスタック形式」と言う．表にするとわかりやすい．
アンスタック形式のデータのイメージ
| 江戸川 | 青梅 | 八王子 |
|---|---|---|
| 3.1 | 0.7 | 0.9 |
| 5.4 | 2.6 | 2.9 |
スタック形式のデータのイメージ
| 気温 | 地点ラベル |
|---|---|
| 3.1 | 江戸川 |
| 0.7 | 青梅 |
| 0.9 | 八王子 |
| 5.4 | 江戸川 |
| 2.6 | 青梅 |
| 2.9 | 八王子 |
まずスタック形式の配列を作り，それをもとにDataFrameを作成，，気温データを文字列から実数にし，分散分析に食わせる，という流れになる．
ols, anova_lmという名前が出てきてわかったと思うが，分散分析は原理的には線形回帰モデルを使っている．
statsmodelsを使った分散分析のコード，7つ目のセルは以下の通り．
import statsmodels.api as sm
import statsmodels.formula.api as smf
# 気温の配列を全部つなげる
all_point_jan = np.concatenate([ome_jan, hachioji_jan, fuchu_jan])
# 観測地点ラベルの配列を作り，全部連結する
ome_label = np.repeat('Ome', 42)
hachioji_label = np.repeat('Hachioji', 42)
fuchu_label = np.repeat('Fuchu', 42)
all_point_label = np.concatenate([ome_label, hachioji_label, fuchu_label])
# 気温の配列と観測地点ラベルの配列を，重ねて2次元配列の表にする．
all_point_table = np.vstack([all_point_jan, all_point_label])
# 配列の次元を入れ替える
# 「気温と観測地点ラベルが対応したペアが，126個存在する」という感じに
all_point_table = all_point_table.transpose()
# 「気温と観測地点ラベルが対応したペアが126個」の2次元配列からDataFrameを作る
# 列の名前はそれぞれ，TempとPointにする
df_statsmodels_anova = pd.DataFrame(data=all_point_table, columns=['Temp', 'Point'])
# Temp列を文字列から実数に置き換える
df_statsmodels_anova['Temp'] = df_statsmodels_anova['Temp'].astype(float)
# DataFrameの確認
print(df_statsmodels_anova.head())
print(df_statsmodels_anova.values)
print(df_statsmodels_anova.values.shape)
# formula(式)を与えて線形回帰を呼び出す
ols_lm = smf.ols(formula='Temp ~ Point', data=df_statsmodels_anova).fit()
# 線形回帰の結果をもとにType2の分散分析をかける
anova_result = sm.stats.anova_lm(ols_lm, typ=2)
print(anova_result)7つ目のセルの結果は以下の通り．
sum_sq df F PR(>F)
Point 39.003333 2.0 15.758491 8.074710e-07
Residual 152.216667 123.0 NaN NaN
PR(>F)がp値なので，scipyを使った時と同じ結果が出ている．
ちなみに分散分析のType2とは，説明変数が2つ以上の重回帰分析のモデルで順番を入れ替えても成り立つタイプの手法を言う．説明変数が2つ以上時（分散分析では二元配置分散分析と呼ぶ）の時は，Type2にする必要がある．
多重検定
この「有意差が出た」時点で平均だけを見て「全部の地点で差がある！」といってしまいがちだが，そうではなく「どこか1つ以上の地点で差がある」だけである．
分散分析自体は，どの群とどの群に差があるのかを出してくれないため，それを見つけるために多重検定というのを行う．
全ての地点を「クロス」で検定するみたいな言葉のイメージで，一番簡単なのが「分布の調整による多重比較法であるTukey-Kramerの多重検定」と呼ばれる手法である．
「Tukey-Kramerの多重検定」は，「各群の母分散も等しい(分散の等質性)」時につかえる．各群のサイズが等しくなくてもよい（この場合だと全ての地点が42個で同じサイズだが）．
「各群の母分散も等しい」とは，簡単に言えば「全ての地点の気温データのばらつきの形が，だいたい同じ」と思っておくと良い（これは統計的にも確かめることができる．母分散の定義は「母集団の特徴を表す母数」参考: 母分散と不偏分散)
ちなみに名前の似ている「Tukeyの方法」という多重検定方法もあるが，これは各群のサイズが同じである必要がある．(pythonのstatmodelsの中に入っているのはTukey-Kramerらしいので，気にせず使える)
多重比較のコード例
Tukey-Kramerの多重比較はstatsmodelsライブラリの中に入っていて，比較的簡単に実行できる．7つ目のセルは以下の通り．
statsmodelを使った分散分析の時と同じく，全てのデータを1次元の配列にし，「観測地点ラベル」の配列を別に作ってペアで与えてやるという形．分散分析の時と同じくデータの配列の各点と地点ラベルは対応していなければならないが，2つの配列で与えてやればよく，分散分析のようにDataFrame形式になっている必要はない
from statsmodels.stats.multicomp import pairwise_tukeyhsd
# まず全ての地点の気温の配列を1つの1次元配列に連結する
# 連結方法は以下のどちらでも良い
all_point_jan = np.concatenate([ome_jan, hachioji_jan, fuchu_jan])
#print(all_point_jan)
all_point_jan = np.hstack([ome_jan, hachioji_jan, fuchu_jan])
#print(all_point_jan)
# 観測地点ラベルの1次元配列を作る
# 平均気温の数だけ，Omeという文字列をリピートした配列を作る
point_label_ome = np.repeat('Ome', len(ome_jan))
#print(point_label_ome)
# 同様に八王子と府中も
point_label_hachioji = np.repeat('Hachioji', len(hachioji_jan))
point_label_fuchu = np.repeat('Fuchu', len(fuchu_jan))
# ラベルの配列を連結する
all_point_label = np.concatenate([point_label_ome, point_label_hachioji, point_label_fuchu])
tukey_result_jan = pairwise_tukeyhsd(all_point_jan, all_point_label)
print(tukey_result_jan)8つ目のセルの結果は以下のようになる．
Multiple Comparison of Means - Tukey HSD, FWER=0.05
=======================================================
group1 group2 meandiff p-adj lower upper reject
-------------------------------------------------------
Fuchu Hachioji -1.0095 0.001 -1.5855 -0.4336 True
Fuchu Ome -1.2976 0.001 -1.8736 -0.7217 True
Hachioji Ome -0.2881 0.465 -0.8641 0.2879 Falseこの結果を見ると，青梅-八王子間のみが，調整p値(出力表中p-adj)が0.05を上回っており，「reject」がFalseになっている．
すなわち，青梅-八王子間は「同じ」という帰無仮説が棄却されなかったことになり，逆に府中-八王子間，府中-青梅間では棄却されていることから「1977年から2018年までの1月平均気温では，この3地点のみでは府中だけが有意に気温が高い」ということができる．
分散分析と多重比較
「多重比較の方でも有意差を出してくれるなら，多重比較だけやればいいのではないか？」という疑問が出てくる．
Tukeyの方法は，分散分析とは内部計算が違う（分散分析はF統計量というものを用いているが，Tukeyの方法では用いていない）ために，「分散分析では有意差が出なかったが，Tukeyの多重比較では有意差が出た」なんてこともありうるし，その逆，つまり「分散分析で有意差が出たが，Tukeyの多重比較では有意差が出なかった」ということもありうる．特に群が増えてくるとTukeyの多重比較では「帰無仮説の数が増えて有意になりにくくなる」．
つまり，分散分析と多重比較は完全に別のものであり，従属関係にはない．「分散分析で有意差がある→多重比較を行う」ではなく，両方を独立させて目的に応じて適用を考える必要がある．
t検定を3回やってはダメな理由
多重比較を使わずに，最初からt検定を3回やればいいのでは？と思うかもしれないが，これはやってはいけない．
青梅，八王子，府中で，2群間の検定をそれぞれ「青梅-八王子間」「八王子-府中間」「府中-青梅間」でそれぞれ「有意水準5%で行う」とすると，「この組み合わせの少なくとも1つは有意差がある」となる確率の計算式は以下のようになる．
1 - (1 - 0.05) * (1 - 0.05) * (1 - 0.05)つまり「0.05を基準に検定したつもりが実際には0.14で検定している」ことになってしまう．0.14の基準は0.05よりもはるかに甘いので「非常に甘く見ている」ということである．
これは当然「比較する群が増えれば増えるほど甘くなる」ということであるから，3群以上の場合は必ず分散分析と多重比較を行う．
二元配置分散分析
今回行なった一元配置分散分析は，前述の通りアンスタック形式の表で言えば以下のようになる．
| 青梅 | 八王子 | 府中 |
|---|---|---|
| 22.8 | 23.7 | 23.8 |
| 26.6 | 27.4 | 27.6 |
| 24.9 | 25.9 | 26.1 |
| 21.7 | 22.4 | 22.5 |
これに対し二元配置分散分析というものもあり，2つの要素の違いを見る．例えば観測地点の名称ではなく，観測時点でのその地点の天気で見ると，以下のようなデータを取り扱うことになる．（平均気温ではなく，その観測時の気温）
| 観測時の天気\観測地点 | 青梅 | 八王子 | 府中 |
|---|---|---|---|
| 晴れ | 26.6 | 27.4 | 27.6 |
| 晴れ | 26.6 | 27.4 | 27.6 |
| 雨 | 24.9 | 25.9 | 26.1 |
| 雨 | 21.7 | 22.4 | 22.5 |
これを分析する時には，以下の3点に注目することになる．
- 地点によって気温に違いはあるのか
- 天気によって気温に違いはあるのか
- 地点と天気の2つの要素による相乗効果はあるのか
この3つ目の「相乗効果」を「交互作用」と呼び，交互作用を分析するのがわざわざ天気を入れた理由であり，これを分析することが二元配置分散分析となる．
二元配置分散分析は前述したようにstatsmodels.stats.anova_lmを線形回帰したモデルに対して適用することで出すことができる．(復習しておくこと)
