こーひーをぶんなぐれ！ - CSVデータからのt検定

本項でもGoogle Colaboratoryを使用する．（もちろんファイル読み込みの部分以外はJupyter Notebookでもよいし，一気にPythonのソースコードをまとめて書いて実行しても問題はない）

東京の異なる2地点，同じ地点での1月と8月の気温

回帰分析の項と同じく，気象庁のオープンデータを使う．

今回は，1977年〜2018年の，江戸川臨海，東京，青梅，八王子，練馬，府中の5箇所の，1月の月平均気温データと，8月の月平均気温データをCSVとしてダウンロードしておく．

これをやるためには2つのファイルに分けるのが一番楽っぽいので，同時点で期間の月のところを変えて，2つCSVファイルをダウンロードする．

ダウンロードした状態は，それぞれ以下から参照すること．

• 東京の複数地点の1月の平均気温データ
• 東京の複数地点の8月の平均気温データ

t検定とは

今回は，東京湾の海風が吹き付ける江戸川区と，内陸の練馬区で，月平均気温が違うかどうかをみる．

小中学校までなら，単純に「1977年から2018年までの平均を出してそれを比較すれば良い」と考えるかもしれないが，回帰分析の項で見た様に，平均気温は年単位で周期的な動きを見せるなど，安定しておらず，単に平均をとって比較するだけでは「それが本当に正しいのか」には大きく疑問が出てしまう．

「2つの気温データ群の間の差を，統計的に（ある程度）保証してくれる」のがt検定である．

「絶対的な保証」ではないが，「ある程度は」保証してくれる，と考えると良い．

2つ群間の検定には，

1. 2群間に対応がある場合(2つの変数の間に何らかの従属関係がある場合)
2. 2群間に対応ない場合の2つ(2つの変数は全く独立)

の2つがある．

1.の「2群間に対応がある場合」は，例えば「食事の前の血糖値と食事後の血糖値との比較」を考えてみると良い．この場合，データとしては「Aさんの食事前と食事後」「Bさんの食事前と食事後」「Cさんの食事前と食事後」……という対になるデータが必要である．

気温で考えてみると「同じ年の（1月の）平均気温で，江戸川と練馬に差があるか」になる．

この場合，月平均気温データが年ごとにペア対応している必要がある．

ペアのデータはこんな感じ「1977年の練馬の1月の気温, 1977年の江戸川の1月の気温」,「1978年の練馬の1月の気温, 1978年の江戸川の1月の気温」,「1979年の練馬の1月の気温, 1979年の江戸川の1月の気温」……このように年ごとに違う地点が並んでいれば良い．（順番はちゃんと並んでいなくても良い）

帰無仮説と対立仮説，p値

この様な統計処理の場合，「否定したい」帰無仮説と，帰無仮説が否定されることで（ある程度）保証される対立仮説の2つで，仮説を定義する．

「観測地点により，年ごとに考えても平均気温に差がある」ことを示したい場合，帰無仮説が「地点での平均気温差がない」であり，その帰無仮説が否定されると（ある程度保証される）のが「観測地点での平均気温の差がある」である．

この「ある程度保証する」の「ある程度」を表すのが，検定の結果出てくるp値という値である．p値が0.05「以下」の時，「5%水準で，有意差がある」（つまり「差がある」）なんて言い方をする．一般的に5%が統計的に「差がある」という水準．これが10%になると「値は色々違うけれど，統計的に立ちかめてみると，実は大した差がない」ぐらいの感覚で捉えておくと良い．

内部計算的には信頼区間とかが出くるが，詳しい説明は統計の教科書を読むこと．

データ前処理と2群間に対応があるt検定のコード例

from google.colab import files
uploaded = files.upload()

これを2回実行して，1月と8月のCSVファイル，tokyo-temp-data-multi-point-jan.csv, tokyo-temp-data-multi-point-aug.csvを，それぞれGoogle Colab上にアップロードしておく．

2つ目のセル．回帰分析の時と同じく，pandasのDataFrameを使って読み込む．

その後余計な行，余計な列を削除する．

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月と8月の配列に，それぞれしておく．Numpyの1次元配列．

値に欠損があるときには，Nanという値が入っている．2群間に関係性のある検定ではデータが対になっている必要があるので，Nanだけを取り除くと対の関係が壊れてしまう．そこで行ごと，つまり対ごと取り除いてしまう．

#まず，江戸川と練馬だけにする
df_jan = df_jan.drop(columns=['東京', '青梅', '八王子', '府中'])
df_aug = df_aug.drop(columns=['東京', '青梅', '八王子', '府中'])
print(df_jan)
print(df_aug)
#両方確認してみると，2009年8月にNanが入っている，
# つまり値が欠損している．

# 欠損していると2群間に対応がある検定はかけられないので，1月と8月の両方からその行を削除する
# （2群間に対応がない検定は，行番号で対応づける必要があるどころか，母数が違っていてもかけられるので，片方だけ抜いても良い）
df_jan = df_jan.drop(34)
df_aug = df_aug.drop(34)

#それぞれ列を取り出して，データの次元を確認する
edogawa_jan = df_jan.loc[:, ['江戸川臨海']]
edogawa_aug = df_aug.loc[:, ['江戸川臨海']]
print('Edogawa Jan. shape:', edogawa_jan.shape)
# 1次元にする
edogawa_jan = edogawa_jan.values.flatten()
edogawa_aug = edogawa_aug.values.flatten()
# 次元を確認
print('Edogawa Jan. shape:', edogawa_jan.shape)
# 中身を確認
print(edogawa_jan)
# printした結果，中身が文字列なので，実数に変更
edogawa_jan = edogawa_jan.astype('float32')
edogawa_aug = edogawa_aug.astype('float32')
# 中身を確認する
print(edogawa_jan)
print(edogawa_aug)
print(edogawa_jan.shape)
print(edogawa_aug.shape)

# 練馬のデータも同じ様に処理する
nerima_jan = df_jan.loc[:, ['練馬']]
nerima_aug = df_aug.loc[:, ['練馬']]
nerima_jan = nerima_jan.values.flatten()
nerima_aug = nerima_aug.values.flatten()
nerima_jan = nerima_jan.astype('float32')
nerima_aug = nerima_aug.astype('float32')
print(nerima_jan)
print(nerima_jan.shape)
print(nerima_aug)
print(nerima_aug.shape)

4つ目のセル．データが揃ったので，とりあえず平均を出してみる．

# データが揃ったので，とりあえず平均を出してみる
print('1977年から1988年までの1月の江戸川の平均気温', np.average(edogawa_jan))
print('1977年から1988年までの8月の江戸川の平均気温', np.average(edogawa_aug))
print('1977年から1988年までの1月の練馬の平均気温', np.average(nerima_jan))
print('1977年から1988年までの8月の練馬の平均気温', np.average(nerima_aug))

4つ目のセルの出力は以下の様になる

1977年から1988年までの1月の江戸川の平均気温 5.597561
1977年から1988年までの8月の江戸川の平均気温 26.360973
1977年から1988年までの1月の練馬の平均気温 4.453658
1977年から1988年までの8月の練馬の平均気温 27.117075

1977年〜2018年の間の月平均気温は，1月は江戸川の方が暖かく，8月は練馬の方が暑い，ということになる．

ここからt検定を行なっていく．5つ目のセル．「2群間に関係性のあるt検定」を行う．

同じ月の江戸川と練馬の平均気温のペアが必要であるが，並びが揃っている2つ1次元配列，すなわち2つの配列の同じ場所に比較したい同じ年の平均気温が入っている形で与えてやる．

それを検定にかける．

# 帰無仮説は「二つの集団の間に差がない」 = 「同じ年の江戸川と練馬の1月（8月）の月平均気温の間には，差がない」
# 対立仮説は「二つの集団の間に差がある」 = 「同じ年の江戸川と練馬の1月（8月）の月平均気温の間には，差がある」
# であり， 結果のp値が0.05を下回っていれば，「差がある」と見るべきである．

# この場合，
# 1. 同じ年で違う地点でのデータが，二つの地点の配列の並びの中で，同じ場所にある
# 2. 「何年」なのかは順不同になる．つまり古い順とか新しい順に並んでいる必要はない
# というデータの形にして食べさせれば良い

import scipy.stats

# これだけで検定を実行できる
result_jan = scipy.stats.ttest_rel(edogawa_jan, nerima_jan)
result_aug = scipy.stats.ttest_rel(edogawa_aug, nerima_aug)

#結果を出力する
print('同じ年の江戸川と練馬の1月の月平均気温の検定', result_jan)
print('同じ年の江戸川と練馬の8月の月平均気温の検定', result_aug)

このセルの出力結果は以下の様になる．

同じ年の江戸川と練馬の1月の月平均気温の検定 Ttest_relResult(statistic=27.58108051178886, pvalue=1.203404241196893e-27)
同じ年の江戸川と練馬の8月の月平均気温の検定 Ttest_relResult(statistic=-11.70687214652232, pvalue=1.698531176507447e-14)

pvalueがp値を指し，これは0.05よりもはるかに小さい．(e-27というのは，10の-27乗を掛ける，という意味である）

ということで，1月と8月双方において，年ごとに見ていって，「江戸川と練馬の気温には有意な差がある」と言うことができる．

2群間に対応がない場合のt検定とコード例

2群間に対応がない場合は，本当にただの値の集合2つの差を見る時に行われる．

先ほどは「年ごと」に対応させて検定したが，対応がない場合「（1977年〜2018年まで「全ての年をひっくるめて」）東京湾の海風が吹き付ける湾岸の江戸川は，風が吹き抜けない練馬に比べ）8月の平均気温が低い」という仮説を検証することになる．

（実際にt検定を行うときは，この「対応がない場合」がほとんどである）

先に出した平均気温を見た時に，確かに江戸川の方が練馬よりも8月の平均気温は低いが，それを統計的に保証するための検定をする．

この場合の，帰無仮説は「（1977年から2018年にかけての）江戸川と練馬の8月の月平均気温の間には，差がない」であり，対立仮説は「（1977年から2018年にかけての）江戸川と練馬の8月の月平均気温には差がある」である．

6番目のセルの内容は以下である．（1月もついでに行なっている）

# 次に「2群間に対応がない場合」を行う

# データが「等分散」というのは，母集団の数値の分散，つまり数値のばらつき方が同じか違うか，ということであり，
# 一般的には「等分散ではない」つまり二つの集団の数値のばらつき具合が異なっていても良い方法を使う

# 分散が同じ場合は「Studentのt検定」と呼び，
# 分散が異なる場合は「Welchのt検定」と呼ぶ

# 帰無仮説は「二つの集団の間に差がない」 = 「江戸川と練馬の1月の月平均気温の間には，差がない」
# 対立仮説は「二つの集団の間に差がある」 = 「江戸川と練馬の1月の月平均気温の間には，差がある」
# であり， 結果のp値が0.05を下回っていれば，「差がある」と見るべきである．

# これだけで実行できる
import scipy.stats
result_jan = scipy.stats.ttest_ind(edogawa_jan, nerima_jan)
result_aug = scipy.stats.ttest_ind(edogawa_aug, nerima_aug)

# 結果を出力する
print('1月の江戸川と練馬の月平均気温の検定', result_jan)
print('8月の江戸川と練馬の月平均気温の検定', result_aug)

この6番目のセルの実行結果は以下の様になる．

1月の江戸川と練馬の月平均気温の検定 Ttest_indResult(statistic=4.6782503978085215, pvalue=1.1599005671463611e-05)
8月の江戸川と練馬の月平均気温の検定 Ttest_indResult(statistic=-2.477117763781257, pvalue=0.015352474172098927)

結果の「statistic」はt統計量，pvalueがp値である．

8月のt検定のp値は0.015であり，0.05を十分に下回っているので，「差がない」帰無仮説は否定される，すなわち「（1977年〜2018年にかけての）8月の月平均気温は，江戸川の方が練馬よりも有意に低い」と言うことができる．

ただし，風が通っているからかどうかはこのデータからは当然わからない． 今回取得した他の「風が通る場所」と「風が通らない場所」を比較していくと，その可能性を検証できるかもしれない．試してみること．

ちなみに，scipyだけではなく，statsmodelsにも当然t検定はあり，全く同じように実行できる．結果の情報量は，statsmodelsの方が1つ多いので，statsmodelsを使う方がいいかもしれない．

statsmodelsの中の「対応がないt検定」を使ったコードは以下の通り．

import statsmodels.stats.weightstats
result_jan = statsmodels.stats.weightstats.ttest_ind(edogawa_jan, nerima_jan)
result_aug = statsmodels.stats.weightstats.ttest_ind(edogawa_aug, nerima_aug)
print('1月の江戸川と練馬の月平均気温の検定', result_jan)
print('8月の江戸川と練馬の月平均気温の検定', result_aug)

statsmodelsの中のttest_indを使った上記のコードの結果は以下のようになる．

1月の江戸川と練馬の月平均気温の検定 (4.678248968183461, 1.1599069103452961e-05, 80.0)
8月の江戸川と練馬の月平均気温の検定 (-2.477104608556971, 0.015352998772965724, 80.0)

結果は，1つ目がt統計量，2つ目がp値，3つ目が自由度が帰ってくる．scipy.stats.ttest_ind関数を使った時と全く同じ結果が出ているので，同じ計算をしていることがわかる．