メインコンテンツへスキップ
見出し画像

Pythonで学ぶデータ分析２：推測統計学×Pythonによる体重分析

    KIYO

    １．概要

    　夏に向けて体重を減らしたい思う人は多いと思いますが、体重は測定のタイミングによる条件で変化します。
    　今回は推測統計学における検定を手法を用いて「各月の体重に変化があるか（太った／痩せた）」を分析します。

    １－１．サンプルデータ

    　私の2~5月までの体重データを使用しました。参考用してPickleファイルを入れましたので、自分でも処理したい方はPandasから読み込みできます。

    ２．理論

    　検定で実施することを簡単に記載しました。詳細は下記記事の通りです。

    ２－１．検定の手順

    　一般的な検定のフローは下記の通りです。

    1. 仮説の設定：帰無仮説H0H_{0}と対立仮説H1H_{1}を立てる

    2. 検定法の選定：分析目的に合った検定方法（検定手法や片側・両側など）を分析者が選定する

    3. 有意水準α・棄却域の決定：有意水準αの決定(一般的にはα=0.05、厳しいときはα=0.01を選定）

    4. 検定統計量と有意点を算出：標本から検定統計量、有意水準に対する有意点を算出

    5. 仮説の判定：検定統計量が棄却域内か確認、検定統計量に対するp値、有意水準を比較

    　今回は各月で体重変化の有無を確認するため下記の通りです。

    1. 仮説の設定

      • 帰無仮説（H0H_{0}）：平均体重に差はない（体重変化無し）

      • 対立仮説（H1H_{1}）：平均体重に差がある（体重変化あり）

    2. 検定法の選定

      • 各月の体重のばらつき（分散）は異なるためウェルチのt検定

      • 体重変化の検定のため両側検定

    3. 有意水準α・棄却域の決定

      • 有意水準α=0.05を選定

    4. 検定統計量と有意点を算出

      • 検定統計量（t値）とp値を計算

    5. 仮説の判定

      • p値が有意水準0.05未満であれば、帰無仮説を棄却して対立仮説を採択

    ２－２．ウェルチのt検定

    　今回の検定法は各月の体重のばらつきが異なることを考慮して、ウェルチのt検定を使用します。計算式は下記の通りです。
    　なおt値の数式変換に関して、通常のウェルチのt検定では母平均（母集団の平均）μ1\mu_1 と μ2\mu_2 の差に関心があり(μ1−μ2=0)(\mu_1 - \mu_2 = 0) と仮定するため変換できます。

    t=(Xˉ1−μ1)−(Xˉ2−μ2)s12n1+s22n2=Xˉ1−Xˉ2s12n1+s22n2t = \frac{(\bar{X}_1 - \mu_1) - (\bar{X}_2 - \mu_2)}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}}=\frac{\bar{X}_1 - \bar{X}_2}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}}

    自由度ν=(s12n1+s22n2)2(s12n1)2n1−1+(s22n2)2n2−1自由度\nu = \frac{(\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2})^2}{\frac{(\frac{s_1^2}{n_1})^2}{n_1 - 1} + \frac{(\frac{s_2^2}{n_2})^2}{n_2 - 1}}

    X1ˉ\bar{X_1},X2ˉ\bar{X_2},：各群のデータの平均値、 s12s_1^2, s22s_2^2：各群の分散、n1n_1, n2n_2：各群のデータ数、ν\nu：ウェルチの自由度

    ３．データの前処理／EDA

    　まずはデータの前処理をして可視化によるEDA（探索的データ分析）を実行します。生データのため前処理の含めて処理していきます。

    ３－１．前処理

    　まずはデータの前処理を実施します。手順は下記の通りです。

    1. 欠損値可視化：可視化用関数（plot_missing_values_ratio）の作成

    2. データ読み込み（下記コードはExcelからの読み込みのため、自分で実行したい方はPickleファイルをDLして"pd.read_pickle('df_weights.pkl')"）

    3. データを月ごとに抽出：月ごとの体重変化を確認するため

    4. 欠損値は全て削除：前日の食事で大きく変わるため補間はしない

    　なおデータの欠損値が多いのは体重計で測るのを忘れていた日であり、体重のばらつきは前日の食事に依存します（特に金・土で外食をする日は1kg近い増加あり）。

    [IN]
def plot_missing_values_ratio(df_list, titlenames, threshold=1, verbose=False):
    if isinstance(df_list, pd.DataFrame):
        df_list = [df_list]
        
    fig, axes = plt.subplots(nrows=1, ncols=len(df_list), figsize=(8 * len(df_list), 8))
    plt.rcParams['font.size'] = 16
    
    if len(df_list) == 1:
        axes = [axes]
        
    for i, (df, titlename) in enumerate(zip(df_list, titlenames)):
        missing_ratios = df.isna().sum() / len(df)
        colors = ['red' if ratio == 1 else 
                  'pink' if ratio >= threshold else 
                  'orange' for ratio in missing_ratios.values]
        
        #欠損値を割合による色分けで棒グラフを描画
        axes[i].bar(missing_ratios.index,
                    missing_ratios.values,
                    color=colors)
        
        #閾値の線を描画
        axes[i].axhline(y=threshold, color='red', linestyle='--', lw=0.5)
        
        axes[i].set_xlabel('Feature', fontsize=12)
        axes[i].set_ylabel('Missing values ratio', fontsize=12)
        axes[i].set_yticks(np.arange(0, 1.1, 0.1))
        axes[i].set_title(f'Missing values ({titlename})', fontsize=12)
        axes[i].tick_params(axis='x', labelrotation=90)
        axes[i].legend(handles=[mpatches.Patch(color='orange'),
                                mpatches.Patch(color='pink'),
                                mpatches.Patch(color='red')], 
                       labels=[f'部分欠損 < {threshold}' ,
                               f'部分欠損 >= {threshold}',
                               '完全欠損'],
                       fontsize=12)
        
        #グリッド線を描画
        if verbose:
                axes[i].grid(axis='y')

    plt.tight_layout()
    plt.show()
    [IN]
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib
from datetime import datetime
import os 
import glob 
from scipy import stats
import matplotlib.patches as mpatches


#ファイルの読み込み
path_excel = 'log_weights_KIYO.xlsx'
df = pd.read_excel(path_excel)
df.to_pickle('df_weights.pkl')
plot_missing_values_ratio(df, ['体重'])

df['Month'] = df['Date'].dt.strftime('%m')
df_Feb, df_Mar, df_Apr, df_May = df[df['Month'] == '02'].copy(), df[df['Month'] == '03'].copy(), df[df['Month'] == '04'].copy(), df[df['Month'] == '05'].copy()
df_Feb.dropna(subset=['Weight'], inplace=True)
df_Mar.dropna(subset=['Weight'], inplace=True)
df_Apr.dropna(subset=['Weight'], inplace=True)  
df_May.dropna(subset=['Weight'], inplace=True)
print(f'前処理後のデータ数　２月：{len(df_Feb)}　３月：{len(df_Mar)}　４月：{len(df_Apr)}　５月：{len(df_May)}')

[OUT]
前処理後のデータ数　２月：12　３月：23　４月：17　５月：5
    画像

    ３－２．統計量計算

    　各種統計量を計算します。また検定の計算に必要なため（前処理後の）データ数も計算しておきます。

    [IN]
#統計量計算
mean_Feb, std_Feb, var_Feb = df_Feb['Weight'].mean(), df_Feb['Weight'].std(), df_Feb['Weight'].var()
mean_Mar, std_Mar, var_Mar = df_Mar['Weight'].mean(), df_Mar['Weight'].std(), df_Mar['Weight'].var()
mean_Apr, std_Apr, var_Apr = df_Apr['Weight'].mean(), df_Apr['Weight'].std(), df_Apr['Weight'].var()
mean_May, std_May, var_May = df_May['Weight'].mean(), df_May['Weight'].std(), df_May['Weight'].var()
count_Feb, count_Mar, count_Apr, count_May = len(df_Feb), len(df_Mar), len(df_Apr), len(df_May)
print(f'平均体重　２月：{mean_Feb:.1f}kg　３月：{mean_Mar:.1f}kg　４月：{mean_Apr:.1f}kg　５月：{mean_May:.1f}kg')
print(f'標準偏差 ２月：{std_Feb:.2f}　３月：{std_Mar:.2f}　４月：{std_Apr:.2f}　５月：{std_May:.2f}')
print(f'分散 ２月：{var_Feb:.2f}　３月：{var_Mar:.2f}　４月：{var_Apr:.2f}　５月：{var_May:.2f}')
print(f'データ数　２月：{count_Feb}個　３月：{count_Mar}個　４月：{count_Apr}個　５月：{count_May}個')

    [OUT]
平均体重　２月：76.1kg　３月：77.5kg　４月：77.4kg　５月：77.0kg
標準偏差 ２月：0.67　３月：0.39　４月：0.66　５月：0.26
分散 ２月：0.45　３月：0.16　４月：0.43　５月：0.07
データ数　２月：12個　３月：23個　４月：17個　５月：5個

    ３－３．可視化

    　欠損値処理後のデータを下記２種の方法で可視化しました。グラフだけ見ると２月が一番痩せていて、３～５月はばらつきの範囲内で変化なしにも見えます。
    　これを統計的に判断していきたいと思います。

    • 時系列変化：日付ごとの体重の変化

    • 正規分布：μ\mu=平均体重, σ2\sigma^2=分散の正規分布

    正規分布f(x)=12πσ2exp⁡(−(x−μ)22σ2)正規分布f(x) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left(-\frac{(x - \mu)^2} {2\sigma^2} \right)

    [IN]
fig, ax = plt.subplots(figsize=(10, 5), facecolor='w')
sns.scatterplot(x='Date', y='Weight', data=df, hue='Month', style='Month', ax=ax)
plt.xticks(rotation=90, ticks=df['Date'][::7])#90°回転+7日ごとに目盛りを表示
plt.grid()
plt.show()

[OUT]
    画像
    [IN]
X = np.linspace(73, 80, 100) #体重の範囲を指定

dist_Feb = stats.norm.pdf(X, loc=mean_Feb, scale=std_Feb)
dist_Mar = stats.norm.pdf(X, loc=mean_Mar, scale=std_Mar)
dist_Apr = stats.norm.pdf(X, loc=mean_Apr, scale=std_Apr)
dist_May = stats.norm.pdf(X, loc=mean_May, scale=std_May)

#可視化
fig, ax = plt.subplots(figsize=(10, 5), facecolor='w')
ax.plot(X, dist_Feb, label='2月')
ax.plot(X, dist_Mar, label='3月')
ax.plot(X, dist_Apr, label='4月')
ax.plot(X, dist_May, label='5月')
ax.set(xlabel='体重', ylabel='確率密度', title='２ー５月の体重の正規分布', xticks=np.arange(73, 81, 1))
ax.grid(), ax.legend()
plt.show()

[OUT]
    画像

    ４．検定

    ４－１．t値、自由度、p値の計算

    　t値、自由度の計算に関しては前述の通り下記式を使用します。

    t=Xˉ1−Xˉ2s12n1+s22n2t =\frac{\bar{X}_1 - \bar{X}_2}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}}

    自由度ν=(s12n1+s22n2)2(s12n1)2n1−1+(s22n2)2n2−1自由度\nu = \frac{(\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2})^2}{\frac{(\frac{s_1^2}{n_1})^2}{n_1 - 1} + \frac{(\frac{s_2^2}{n_2})^2}{n_2 - 1}}

    X1ˉ\bar{X_1},X2ˉ\bar{X_2},：各群のデータの平均値、 s12s_1^2, s22s_2^2：各群の分散、n1n_1, n2n_2：各群のデータ数、ν\nu：ウェルチの自由度

    　得られたt値からp値を計算しますがこちらは"stats.t.sf()"を使用しました。なおp値を求めるためのt分布に対して自由度を渡すのは、t分布が自由度の値によって形状が変化するためです。
    　なお下記ではt値、p値を自作関数で作成しておりますが、同様のことを"stats.ttest_ind"でも実行可能です（次節で検算のため使用）。

    [IN]
def welch_t_test(mean1, var1, n1, mean2, var2, n2):
    t_value = (mean1 - mean2) / np.sqrt((var1 / n1) + (var2 / n2))
    d_free = ((var1 / n1 + var2 / n2) ** 2) / (((var1 / n1) ** 2) / (n1 - 1) + ((var2 / n2) ** 2) / (n2 - 1))  # ウェルチの自由度
    return t_value, d_free

def print_welch_t_test_results(t_value, d_free, p_value, months):
    print(f"t値({months}): {t_value:.3f}")
    print(f"ウェルチの自由度({months}): {int(d_free)}")
    print(f"p値({months}): {p_value:.5f}")

    if p_value < 0.05:
        print(f"有意水準 0.05 で帰無仮説を棄却：{months}のデータ平均に有意な差あり\n")
    else:
        print(f"有意水準 0.05 で帰無仮説を棄却できず：{months}のデータ平均に有意な差なし\n")

# 2月-3月の検定結果
tvalue_Feb_Mar, d_free_Feb_Mar = welch_t_test(mean_Feb, var_Feb, count_Feb, mean_Mar, var_Mar, count_Mar)
pvalue_Feb_Mar = stats.t.sf(np.abs(tvalue_Feb_Mar), d_free_Feb_Mar) * 2 # p値
print_welch_t_test_results(tvalue_Feb_Mar, d_free_Feb_Mar, pvalue_Feb_Mar, "2月-3月")
# 3月-4月の検定結果
tvalue_Mar_Apr, d_free_Mar_Apr = welch_t_test(mean_Mar, var_Mar, count_Mar, mean_Apr, var_Apr, count_Apr)
pvalue_Mar_Apr = stats.t.sf(np.abs(tvalue_Mar_Apr), d_free_Mar_Apr) * 2 # p値
print_welch_t_test_results(tvalue_Mar_Apr, d_free_Mar_Apr, pvalue_Mar_Apr, "3月-4月")
# 4月-5月の検定結果
tvalue_Apr_May, d_free_Apr_May = welch_t_test(mean_Apr, var_Apr, count_Apr, mean_May, var_May, count_May)
pvalue_Apr_May = stats.t.sf(np.abs(tvalue_Apr_May), d_free_Apr_May) * 2 # p値
print_welch_t_test_results(tvalue_Apr_May, d_free_Apr_May, pvalue_Apr_May, "4月-5月")
    [OUT]

前処理後のデータ数　２月：12　３月：23　４月：17　５月：5
平均体重　２月：76.1kg　３月：77.5kg　４月：77.4kg　５月：77.0kg
標準偏差 ２月：0.67　３月：0.39　４月：0.66　５月：0.26
分散 ２月：0.45　３月：0.16　４月：0.43　５月：0.07
データ数　２月：12個　３月：23個　４月：17個　５月：5個
t値(2月-3月): -6.281
ウェルチの自由度(2月-3月): 15
p値(2月-3月): 0.00001
有意水準 0.05 で帰無仮説を棄却：2月-3月のデータ平均に有意な差あり

t値(3月-4月): 0.160
ウェルチの自由度(3月-4月): 24
p値(3月-4月): 0.87412
有意水準 0.05 で帰無仮説を棄却できず：3月-4月のデータ平均に有意な差なし

t値(4月-5月): 2.138
ウェルチの自由度(4月-5月): 17
p値(4月-5月): 0.04708
有意水準 0.05 で帰無仮説を棄却：4月-5月のデータ平均に有意な差あり

    　結果として下記の通りとなりました。

    • 2-3月はp値<0.05のため、体重が増えたと思われる

    • 3-4月はp値>0.05のため、体重に変化がなかったと思われる

    • 2-3月はp値<0.05のため、体重が減ったと思われる

    ４－２．可視化

    　自由度を渡した時のt分布における棄却域αの領域とt値の関係を可視化しました。計算したt値が棄却域に入っていなければ帰無仮説は棄却できないため、差があるとは言えないという結論になります。

    　なお本節ではt値、p値計算のために"scipy.stats.ttest_ind()"を使用しており、自作関数と同じ結果であることを確認しました。

    [API]
scipy.stats.ttest_ind(a, b, axis=0, equal_var=True, nan_policy='propagate', 
                      permutations=None, random_state=None, 
                      alternative='two-sided', trim=0)
    [IN]
def visualize_t_test(t_value, d_free, title):
    X_t = np.linspace(-4, 4, 1000)
    t_dist = stats.t.pdf(X_t, d_free)

    fig, ax = plt.subplots(figsize=(10, 5), facecolor='w')
    ax.plot(X_t, t_dist, label='t分布')

    alpha = 0.05 # 有意水準
    critical_value = stats.t.ppf(1 - alpha / 2, d_free) # 両側検定の場合は、有意水準を2で割る
    reject_region_left = np.where(X_t < -critical_value) # 棄却域
    reject_region_right = np.where(X_t > critical_value) # 棄却域
    ax.fill_between(X_t[reject_region_left], t_dist[reject_region_left], color='red', alpha=0.3, label='棄却域（α）')
    ax.fill_between(X_t[reject_region_right], t_dist[reject_region_right], color='red', alpha=0.3)

    if abs(t_value) > 4:
        max_xtick, min_xtick = int(np.ceil(abs(t_value))), -int(np.ceil(abs(t_value)))
    else:
        max_xtick, min_xtick = 4, -4
    
    ax.set(xlabel='t値', ylabel='確率密度', title=title, xticks=np.arange(min_xtick, max_xtick + 1, 1))
    ax.axvline(x=t_value, color='green', linestyle='--', label=f"t値: {t_value:.3f}")
    ax.legend(), ax.grid()
    plt.show()


# 2月-3月の検定
t_Feb_Mar, pvalue_Feb_Mar = stats.ttest_ind(df_Feb['Weight'], df_Mar['Weight'], equal_var=False)
visualize_t_test(t_Feb_Mar, d_free_Feb_Mar, '2月-3月のt検定')
# 3月-4月の検定
t_Mar_Apr, pvalue_Mar_Apr = stats.ttest_ind(df_Mar['Weight'], df_Apr['Weight'], equal_var=False)
visualize_t_test(t_Mar_Apr, d_free_Mar_Apr, '3月-4月のt検定')
# 4月-5月の検定
t_Apr_May, pvalue_Apr_May = stats.ttest_ind(df_Apr['Weight'], df_May['Weight'], equal_var=False)
visualize_t_test(t_Apr_May, d_free_Apr_May, '4月-5月のt検定')

#計算結果の再確認
print(f'再確認：t値（2月-3月）：{t_Feb_Mar:.3f}, p値（2月-3月）：{pvalue_Feb_Mar:.5f}')
print(f'再確認：t値（3月-4月）：{t_Mar_Apr:.3f}, p値（3月-4月）：{pvalue_Mar_Apr:.5f}')
print(f'再確認：t値（4月-5月）：{t_Apr_May:.3f}, p値（4月-5月）：{pvalue_Apr_May:.5f}')
    [OUT]
再確認：t値（2月-3月）：-6.281, p値（2月-3月）：0.00001
再確認：t値（3月-4月）：0.160, p値（3月-4月）：0.87412
再確認：t値（4月-5月）：2.138, p値（4月-5月）：0.04708
    画像
    画像
    画像

    ５．考察

    　今回の結果において、２－３月は直感通り体重が増加していましたが、４－５月は直感とは外れており体重減少が確認されています。再度グラフを並べてみました。

    画像
    画像
    画像

    　これをt値の観点で確認します。「体重変化がある＝帰無仮説が棄却される=t値が大きい」となるため、”差がある”ためには下記が重要です。

    • 平均の差が大きい：Xˉ1−Xˉ2\bar{X}_1 - \bar{X}_2

    • 分散が小さい：s12,s22s_1^2, s_2^2

    • サンプルサイズが大きい：n1,n2n_1, n_2

    t=Xˉ1−Xˉ2s12n1+s22n2t =\frac{\bar{X}_1 - \bar{X}_2}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}}

     　原因として、５月のデータは（GW中で摂生していたため）分散が小さかったところが影響していると思います。ただしp値を見ても棄却域ギリギリであり、データ数も少ないためデータを追加すると（多分体重もばらつくため）有意差が出ない可能性もあります。
    　継続して痩せるのであればここで安心せずに引き続きダイエットに頑張る必要があるとは思います。

    【コラム：p-hacking】
    　p-hackingとは意図的にp値を下げる(t値を上げる)ことで有意差を作り出す不正行為です。例えば新薬に対して性能差があることを説明したいために意図的にp値を調整するような事例があります。具体例は下記の通りです。

    1. 意図的なデータ抽出：特定のデータを選択的に除外または含めることで平均値の差を調整する

    2. 意図的なデータサンプリング：性能差を出すために別の影響がある部分を隠してデータのサンプリングを実施する（研究者が調査対象者の中から特定の属性（年齢、性別、健康状態など）を持つ人々を選択的に調査など）

    3. 解析方法の変更：検定方法やモデルの選択を繰り返し変更することで、望ましい結果が得られるまで解析を続ける。

    4. 多重検定の問題：「統計的仮説検定における多重検定の問題とその対処法」参照

    5. サンプルサイズの調整：サンプルサイズが大きいと標準誤差が小さくなるためt値が増加します。つまり有意差を出すためにサンプルサイズを過剰に増やすことはp-hackingにつながります。

    参考資料

    あとがき

    　夏までに目指せ７５ｋｇ！　

     
     

    KIYO

     
     
    普段は製造業で企画/開発/設計しております。記事はプログラミング・機械学習、IoT関係の記事をメインで作成し、なるべく１つの記事で知りたいことを網羅していきます。内容は学術的より実装・アウトプット（ほしくなるもの）を重視して作成しています。 面白そうな仕事があればやりたいです！

    あなたへのおすすめ