Pythonは非常に汎用性の高い言語であり、幅広い分野の様々なタスクに有用です。そのような分野の1つがデータセットの統計分析であり、SPSSと並んでPythonは統計のための最も一般的なツールの1つです。
Pythonのユーザーフレンドリーで直感的な性質は、特にstatsmodelsライブラリの使用により、統計テストの実行や分析手法の実装を容易にしています。
Python の statsmodels ライブラリを紹介します。
statsmodelsライブラリはPythonのモジュールで、統計的なテストやデータ探索のための様々な統計ツールに簡単にアクセスすることができるようになります。このライブラリでは、通常の最小二乗（OLS）回帰、一般化線形モデル、ロジットモデル、主成分分析（PCA）、自己回帰統合移動平均（ARIMA）モデルなど、数多くの統計テストや関数にアクセスすることができるようになっています。
モデルの結果が正確であることを保証するために、常に他の統計パッケージと比較してテストされています。SciPyやPandasと組み合わせれば、データの可視化、統計テストの実行、関係の有意性の確認が簡単にできます。
データセットの選択
Pythonで統計の練習をする前に、データセットを選択する必要があります。今回は、Gapminder Foundationによって編集されたデータセットを使用します。
Gapminder データセットは、世界中の国の人々の一般的な健康と幸福を評価するために使用される多くの変数を追跡します。このデータセットは、非常によく文書化され、標準化され、完全であるため、私たちはこのデータセットを使用することにしました。このデータセットを利用するために、前処理をする必要はあまりない。
回帰分析、ANOVA、その他のテストを実行するためにデータセットを準備する必要があることはいくつかありますが、全体としてデータセットはすぐに使える状態になっています。
Gapminderデータセットの統計解析の出発点は、探索的データ解析です。MatplotlibとSeabornのいくつかのグラフとプロット関数を使って、いくつかの興味深い関係を視覚化し、どのような変数関係を調べたいかのアイデアを得ます。
探索的データ解析と前処理
まず、考えられる関係を可視化することから始めます。Seaborn と Pandas を使って、データセットの変数間の相関の強さを見る回帰をいくつか行い、どの変数関係が研究する価値があるかのアイデアを得ることができます。
この2つのライブラリと、その他に使用するライブラリをインポートします。
import statsmodels.formula.api as smf
import statsmodels.stats.multicomp as multi
import scipy
from scipy.stats import pearsonr
import pandas as pd
from seaborn import regplot
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
前処理はあまり必要ありませんが、いくつか必要なことがあります。まず、データの欠落やnullをチェックし、非数値の項目は数値に変換しておきます。また、変換後のデータフレームのコピーを作成し、それを使って作業することになります。
# Check for missing data
def check_missing_values(df, cols):
for col in cols:
print("Column {} is missing:".format(col))
print((df[col].values == ' ').sum())
print()
# Convert to numeric
def convert_numeric(dataframe, cols):
for col in cols:
dataframe[col] = pd.to_numeric(dataframe[col], errors='coerce')
df = pd.read_csv("gapminder.csv")
print("Null values:")
print(df.isnull().values.any())
cols = ['lifeexpectancy', 'breastcancerper100th', 'suicideper100th']
norm_cols = ['internetuserate', 'employrate']
df2 = df.copy()
check_missing_values(df, cols)
check_missing_values(df, norm_cols)
convert_numeric(df2, cols)
convert_numeric(df2, norm_cols)
以下はその出力です。
Null values:
Column lifeexpectancy is missing:
22
Column breastcancerper100th is missing:
40
Column suicideper100th is missing:
22
Column internetuserate is missing:
21
Column employrate is missing:
35
一握りの欠損値がありますが、数値変換によってそれらを NaN 値に変換し、データセットに対して探索的なデータ分析を実行できるようにします。
具体的には、インターネット利用率と平均寿命の関係や、インターネット利用率と就業率の関係などを分析してみることができます。SeabornとMatplotlibを使って、これらの関係のいくつかを個別にグラフ化してみましょう。
sns.lmplot(x="internetuserate", y="breastcancerper100th", data=df2, fit_reg=False)
plt.title("Internet Use Rate and Breast Cancer Per 100k")
plt.show()
sns.lmplot(x="internetuserate", y="lifeexpectancy", data=df2, fit_reg=False)
plt.title("Internet Use Rate and Life Expectancy")
plt.show()
sns.lmplot(x="internetuserate", y="employrate", data=df2, fit_reg=False)
plt.title("Internet Use Rate and Employment Rate")
plt.show()
グラフの結果は以下の通りです。
さらに調べると、興味深い関係がありそうです。興味深いことに、インターネット利用率と乳がんの間には、かなり強い正の関係があるように見えますが、これは、テクノロジーへのアクセスがより良い国々でより良いテストが行われたことによる単なるアーチファクトである可能性があります。
また、平均寿命とインターネット利用率の間にも、直線的ではないものの、かなり強い関係があるように思われます。
最後に、インターネット利用率と就業率には放物線を描くような非線形の関係があるように見える。
適切な仮説の選択
さらに調べる価値のある関係を選び出したい。いろいろな関係が考えられるので、仮説を立てて、統計的検定で関係を探ってみる。仮説を立て、2つの変数の間で相関検定を行い、相関検定が有意であった場合、その相関がどれくらい強いのか、2つの変数の間の相関が単なる偶然以上のものであると確実に言えるのか、統計的検定を行う必要があります。
使用する統計検定の種類は、説明変数と応答変数（独立変数と従属変数）の性質に依存します。ここでは、3つの異なるタイプの統計的検定を実行する方法を説明します。
- ANOVAs（分散分析
ANOVAs（分散分析） * Chi-Square Tests（カイ二乗検定 - 回帰
私たちは、上記で可視化したものを使って、インターネット使用率と平均寿命の関係を調べることにします。
帰無仮説は、インターネット利用率と平均寿命の間に有意な関係はないというもので、我々の仮説は、この2つの変数の間に関係があるというものです。
これからデータセットに対して様々なタイプの仮説検証を行うことになります。使用する仮説検定の種類は、説明変数と応答変数の性質に依存します。説明変数と応答変数の組み合わせが異なると、異なる統計的検定が必要になります。たとえば、1つの変数がカテゴリーで、1つの変数が量的な性質を持っている場合、分散分析が必要です。
分散分析(ANOVA)
分散分析（ANOVA）は、分散分析によって決定された2つ以上の平均を比較するために採用される統計検定です。一元配置分散分析は、グループ間の差異を分析し、その差異が統計的に有意であるかどうかを判断するために利用される。
一元配置分散分析は、2つ以上の独立したグループの平均を比較しますが、実際には、少なくとも3つの独立したグループがあるときに最もよく使われます。
GapminderデータセットでANOVAを実行するために、いくつかの特徴を変換する必要があります。データセット中のこれらの値は連続ですが、ANOVA分析は、1つの変数がカテゴリーで、1つの変数が量的である状況に適しているからです。
我々は、カテゴリを選択して、問題の変数をビン化し、それをパーセンタイルに分割することによって、データを連続から量的に変換することができます。独立変数はカテゴリ変数に変換され、従属変数は連続のままとなります。Pandasの qcut() 関数を使ってデータフレームをビンに分割することができます。
def bin(dataframe, cols):
# Create new columns that store the binned data
for col in cols:
new_col_name = "{}_bins".format(col)
dataframe[new_col_name] = pd.qcut(dataframe[col], 10, labels=["1=10%", "2=20%", "3=30%", "4=40%", "5=50%", "6=60%", "7=70%", "8=80", "9=90%", "10=100%"])
df3 = df2.copy()
# This creates new columns filled with the binned column data
bin(df3, cols)
bin(df3, norm_cols)
変数が変換され、分析の準備ができたら、statsmodelライブラリを使って、選択された特徴量のANOVAを実行することができます。ANOVAの結果を出力し、2つの変数の関係が統計的に有意であるかどうかを確認します。
anova_df = df3[['lifeexpectancy', 'internetuserate_bins', 'employrate_bins']].dropna()
relate_df = df3[['lifeexpectancy', 'internetuserate_bins']]
anova = smf.ols(formula='lifeexpectancy ~ C(internetuserate_bins)', data=anova_df).fit()
print(anova.summary())
# We may also want to check the mean and standard deviation for the groups
mean = relate_df.groupby("internetuserate_bins").mean()
sd = relate_df.groupby("internetuserate_bins").std()
print(mean)
print(sd)
以下はモデルの出力です。
OLS Regression Results
==============================================================================
Dep. Variable: lifeexpectancy R-squared: 0.689
Model: OLS Adj. R-squared: 0.671
Method: Least Squares F-statistic: 38.65
Date: Mon, 11 May 2020 Prob (F-statistic): 1.71e-35
Time: 17:49:24 Log-Likelihood: -521.54
No. Observations: 167 AIC: 1063.
Df Residuals: 157 BIC: 1094.
Df Model: 9
Covariance Type: nonrobust
======================================================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------------------------------
Intercept 56.6603 1.268 44.700 0.000 54.157 59.164
C(internetuserate_bins)[T.2=20%] 1.6785 1.870 0.898 0.371 -2.015 5.372
C(internetuserate_bins)[T.3=30%] 5.5273 1.901 2.907 0.004 1.772 9.283
C(internetuserate_bins)[T.4=40%] 11.5693 1.842 6.282 0.000 7.932 15.207
C(internetuserate_bins)[T.5=50%] 14.6991 1.870 7.860 0.000 11.005 18.393
C(internetuserate_bins)[T.6=60%] 16.7287 1.870 8.946 0.000 13.035 20.422
C(internetuserate_bins)[T.7=70%] 17.8802 1.975 9.052 0.000 13.978 21.782
C(internetuserate_bins)[T.8=80] 19.8302 1.901 10.430 0.000 16.075 23.586
C(internetuserate_bins)[T.9=90%] 23.0723 1.901 12.135 0.000 19.317 26.828
C(internetuserate_bins)[T.10=100%] 23.3042 1.901 12.257 0.000 19.549 27.060
==============================================================================
Omnibus: 10.625 Durbin-Watson: 1.920
Prob(Omnibus): 0.005 Jarque-Bera (JB): 11.911
Skew: -0.484 Prob(JB): 0.00259
Kurtosis: 3.879 Cond. No. 10.0
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
このモデルは、1.71e-35というとても小さなP値（Prob F-statistic）を与えていることがわかります。これは、通常の有意水準である0.05よりはるかに小さいので、平均余命とインターネット利用率の間に有意な関係があると結論づけられます。
相関のP値は有意であるように見えますが、10種類のカテゴリーがあるので、ポストホック検定を実行して、タイプ1エラーをチェックした後でも平均値の差が有意であることをチェックしたいと思います。Tukey Honestly Significant Difference (Tukey HSD) 検定を利用して、multicompモジュールの助けを借りて、ポストホック検定を実行することができます。
multi_comparison = multi.MultiComparison(anova_df["lifeexpectancy"], anova_df["internetuserate_bins"])
results = multi_comparison.tukeyhsd()
print(results)
以下はその結果です。
Multiple Comparison of Means - Tukey HSD, FWER=0.05
=======================================================
group1 group2 meandiff p-adj lower upper reject
-------------------------------------------------------
10=100% 1=10% -23.3042 0.001 -29.4069 -17.2015 True
10=100% 2=20% -21.6257 0.001 -27.9633 -15.2882 True
10=100% 3=30% -17.7769 0.001 -24.2097 -11.344 True
10=100% 4=40% -11.7349 0.001 -17.9865 -5.4833 True
10=100% 5=50% -8.6051 0.001 -14.9426 -2.2676 True
10=100% 6=60% -6.5755 0.0352 -12.913 -0.238 True
10=100% 7=70% -5.4241 0.2199 -12.0827 1.2346 False
10=100% 8=80 -3.4741 0.7474 -9.9069 2.9588 False
10=100% 9=90% -0.2319 0.9 -6.6647 6.201 False
1=10% 2=20% 1.6785 0.9 -4.3237 7.6807 False
1=10% 3=30% 5.5273 0.1127 -0.5754 11.6301 False
1=10% 4=40% 11.5693 0.001 5.6579 17.4807 True
1=10% 5=50% 14.6991 0.001 8.6969 20.7013 True
1=10% 6=60% 16.7287 0.001 10.7265 22.7309 True
1=10% 7=70% 17.8801 0.001 11.5399 24.2204 True
1=10% 8=80 19.8301 0.001 13.7274 25.9329 True
1=10% 9=90% 23.0723 0.001 16.9696 29.1751 True
2=20% 3=30% 3.8489 0.6171 -2.4887 10.1864 False
2=20% 4=40% 9.8908 0.001 3.7374 16.0443 True
2=20% 5=50% 13.0206 0.001 6.7799 19.2614 True
2=20% 6=60% 15.0502 0.001 8.8095 21.291 True
2=20% 7=70% 16.2017 0.001 9.6351 22.7683 True
2=20% 8=80 18.1517 0.001 11.8141 24.4892 True
2=20% 9=90% 21.3939 0.001 15.0563 27.7314 True
3=30% 4=40% 6.042 0.0678 -0.2096 12.2936 False
3=30% 5=50% 9.1718 0.001 2.8342 15.5093 True
3=30% 6=60% 11.2014 0.001 4.8638 17.5389 True
3=30% 7=70% 12.3528 0.001 5.6942 19.0114 True
3=30% 8=80 14.3028 0.001 7.87 20.7357 True
3=30% 9=90% 17.545 0.001 11.1122 23.9778 True
4=40% 5=50% 3.1298 0.8083 -3.0237 9.2833 False
4=40% 6=60% 5.1594 0.1862 -0.9941 11.3129 False
4=40% 7=70% 6.3108 0.0638 -0.1729 12.7945 False
4=40% 8=80 8.2608 0.0015 2.0092 14.5124 True
4=40% 9=90% 11.503 0.001 5.2514 17.7546 True
5=50% 6=60% 2.0296 0.9 -4.2112 8.2704 False
5=50% 7=70% 3.181 0.8552 -3.3856 9.7476 False
5=50% 8=80 5.131 0.2273 -1.2065 11.4686 False
5=50% 9=90% 8.3732 0.0015 2.0357 14.7108 True
6=60% 7=70% 1.1514 0.9 -5.4152 7.718 False
6=60% 8=80 3.1014 0.8456 -3.2361 9.439 False
6=60% 9=90% 6.3436 0.0496 0.0061 12.6812 True
7=70% 8=80 1.95 0.9 -4.7086 8.6086 False
7=70% 9=90% 5.1922 0.2754 -1.4664 11.8508 False
8=80 9=90% 3.2422 0.8173 -3.1907 9.675 False
-------------------------------------------------------
これで、比較の中でどのグループが統計的に有意な差を持っているのか、よりよく理解できるようになりました。
もし reject 列に False というラベルがあれば、帰無仮説を棄却して、比較される2つのグループの間に有意差があると仮定することが推奨されることがわかります。
独立性のカイ二乗検定
ANOVAは、1つの変数が連続的で、もう1つがカテゴリ的であるような場合に適切です。今、我々は独立性のカイ2乗検定を実行する方法を見ます。
独立性のカイ2乗検定は、説明変数と応答変数の両方がカテゴリカルであるときに利用されます。また、説明変数が量的で、応答変数がカテゴリで、説明変数をカテゴリに分割することによってできる場合にも、カイ2乗検定を使用したいと思うかもしれません。
独立性のカイ2乗検定は、2つのカテゴリ変数の間の関係がどれだけ有意であるかを分析するために使用される統計的検定である。カイ2乗検定を実行すると、1つの変数のすべてのカテゴリが、2番目の変数のカテゴリに対して、その頻度を比較されます。これは、データが度数表として表示され、行が独立変数、列が従属変数を表すことを意味します。
独立変数を（ビニングして）カテゴリ変数に変換したように、ANOVA検定では、カイ2乗検定を実行するために、両方の変数をカテゴリ変数にする必要があります。この問題の仮説は、前の問題の仮説と同じで、「平均寿命とインターネット利用率の間に有意な関係がある」というものです。
ここではシンプルに、インターネット利用率という変数を2つに分けてみますが、もっと多くすることも可能です。そのための関数を書きます。
タイプ1エラー（偽陽性）を防ぐために、ボンフェローニ調整と呼ばれるアプローチでポストホック比較を行うことにします。これを行うには、応答変数の異なる可能性のあるペアについて比較を実行し、次にそれらの調整された有意性をチェックします。
ここでは、すべての異なる可能なペアの比較を実行するのではなく、どのように実行できるかを示します。再符号化スキームを用いていくつかの異なる比較を行い、レコードを新しい特徴列にマッピングします。
その後、観測されたカウントをチェックし、それらの比較のテーブルを作成することができる。
def half_bin(dataframe, cols):
for col in cols:
new_col_name = "{}_bins_2".format(col)
dataframe[new_col_name] = pd.qcut(dataframe[col], 2, labels=["1=50%", "2=100%"])
half_bin(df3, ['internetuserate'])
# Recoding scheme
recode_2 = {"3=30%": "3=30%", "7=70%": "7=70%"}
recode_3 = {"2=20%": "2=20%", "8=80": "8=80"}
recode_4 = {"6=60%": "6=60%", "9=90%": "9=90%"}
recode_5 = {"4=40%": "4=40%", "7=70%": "7=70%"}
# Create the new features
df3['Comp_3v7'] = df3['lifeexpectancy_bins'].map(recode_2)
df3['Comp_2v8'] = df3['lifeexpectancy_bins'].map(recode_3)
df3['Comp_6v9'] = df3['lifeexpectancy_bins'].map(recode_4)
df3['Comp_4v7'] = df3['lifeexpectancy_bins'].map(recode_5)
カイ二乗検定とポストホック比較の実行は、まずクロスタブ比較表を作成することを含みます。クロスタブ比較表は、説明変数の異なる水準での応答変数の出現のパーセンテージを表示します。
これがどのように機能するかを知るために、すべての平均余命ビンの比較の結果をプリントアウトしてみましょう。
# Get table of observed counts
count_table = pd.crosstab(df3['internetuserate_bins_2'], df3['lifeexpectancy_bins'])
print(count_table)
lifeexpectancy_bins 1=10% 2=20% 3=30% 4=40% ... 7=70% 8=80 9=90% 10=100%
internetuserate_bins_2 ...
1=50% 18 19 16 14 ... 4 4 1 0
2=100% 0 0 1 4 ... 15 11 16 19
クロスタブの比較は、ある変数のカテゴリの頻度を、2番目の変数でチェックすることがわかります。上では、我々が作成した2つのビンのうちの1つに入る状況での平均寿命の分布がわかります。
そして、これはカイ2乗検定を実行するものなので、我々は上で作成したさまざまなペアのクロス集計を計算する必要があります。
count_table_3 = pd.crosstab(df3['internetuserate_bins_2'], df3['Comp_3v7'])
count_table_4 = pd.crosstab(df3['internetuserate_bins_2'], df3['Comp_2v8'])
count_table_5 = pd.crosstab(df3['internetuserate_bins_2'], df3['Comp_6v9'])
count_table_6 = pd.crosstab(df3['internetuserate_bins_2'], df3['Comp_4v7'])
カイ二乗検定を実行できるように変数を変換したら、 statsmodelの chi2_contingency 関数を使用して検定を実行します。
カイ二乗検定の結果と同様に、列のパーセンテージも出力したいので、そのための関数を作成します。そして、作成した関数を使用して、作成した4つの比較表に対してカイ二乗検定を行います。
def chi_sq_test(table):
print("Results for:")
print(str(table))
# Get column percentages
col_sum = table.sum(axis=0)
col_percents = table/col_sum
print(col_percents)
chi_square = scipy.stats.chi2_contingency(table)
print("Chi-square value, p-value, expected_counts")
print(chi_square)
print()
print("Initial Chi-square:")
chi_sq_test(count_table)
print(" ")
chi_sq_test(count_table_3)
chi_sq_test(count_table_4)
chi_sq_test(count_table_5)
chi_sq_test(count_table_6)
その結果がこちらです。
Initial Chi-square:
Results for:
lifeexpectancy_bins 1=10% 2=20% 3=30% 4=40% ... 7=70% 8=80 9=90% 10=100%
internetuserate_bins_2 ...
1=50% 18 19 16 14 ... 4 4 1 0
2=100% 0 0 1 4 ... 15 11 16 19
[2 rows x 10 columns]
lifeexpectancy_bins 1=10% 2=20% 3=30% ... 8=80 9=90% 10=100%
internetuserate_bins_2 ...
1=50% 1.0 1.0 0.941176 ... 0.266667 0.058824 0.0
2=100% 0.0 0.0 0.058824 ... 0.733333 0.941176 1.0
[2 rows x 10 columns]
Chi-square value, p-value, expected_counts
(102.04563740451277, 6.064860600653971e-18, 9, array([[9.45251397, 9.97765363, 8.9273743 , 9.45251397, 9.45251397,
9.97765363, 9.97765363, 7.87709497, 8.9273743 , 9.97765363],
[8.54748603, 9.02234637, 8.0726257 , 8.54748603, 8.54748603,
9.02234637, 9.02234637, 7.12290503, 8.0726257 , 9.02234637]]))
-----
Results for:
Comp_3v7 3=30% 7=70%
internetuserate_bins_2
1=50% 16 4
2=100% 1 15
Comp_3v7 3=30% 7=70%
internetuserate_bins_2
1=50% 0.941176 0.210526
2=100% 0.058824 0.789474
Chi-square value, p-value, expected_counts
(16.55247678018576, 4.7322137795376575e-05, 1, array([[ 9.44444444, 10.55555556],
[ 7.55555556, 8.44444444]]))
-----
Results for:
Comp_2v8 2=20% 8=80
internetuserate_bins_2
1=50% 19 4
2=100% 0 11
Comp_2v8 2=20% 8=80
internetuserate_bins_2
1=50% 1.0 0.266667
2=100% 0.0 0.733333
Chi-square value, p-value, expected_counts
(17.382650301643437, 3.0560286589975315e-05, 1, array([[12.85294118, 10.14705882],
[ 6.14705882, 4.85294118]]))
-----
Results for:
Comp_6v9 6=60% 9=90%
internetuserate_bins_2
1=50% 6 1
2=100% 13 16
Comp_6v9 6=60% 9=90%
internetuserate_bins_2
1=50% 0.315789 0.058824
2=100% 0.684211 0.941176
Chi-square value, p-value, expected_counts
(2.319693757720874, 0.12774517376836148, 1, array([[ 3.69444444, 3.30555556],
[15.30555556, 13.69444444]]))
-----
Results for:
Comp_4v7 4=40% 7=70%
internetuserate_bins_2
1=50% 14 4
2=100% 4 15
Comp_4v7 4=40% 7=70%
internetuserate_bins_2
1=50% 0.777778 0.210526
2=100% 0.222222 0.789474
Chi-square value, p-value, expected_counts
(9.743247922437677, 0.0017998260000241526, 1, array([[8.75675676, 9.24324324],
[9.24324324, 9.75675676]]))
-----
全数表の結果だけを見ると、P値は6.064860600653971e-18のように見えます。
しかし、異なるグループが互いにどのように乖離しているかを確認するために、データフレーム内の異なるペアでカイ二乗検定を実行する必要があります。我々は、選択した異なるペアのそれぞれについて、統計的に有意な差があるかどうかを確認します。有意な結果を示すP値は、いくつの比較を行うかによって変わることに注意してください。
6対9の比較ではP値が0.127となり、これは0.05のしきい値を超えているため、そのカテゴリーの差は有意でない可能性があることを示しています。比較の違いを見ることで、なぜ異なるレベルを互いに比較する必要があるのかを理解することができます。
ピアソン相関
我々は、カテゴリ説明変数と量的応答変数があるときに使用するテスト（ANOVA）と、2つのカテゴリ変数があるときに使用するテスト（カイ2乗）をカバーしました。
ここで、量的説明変数と量的応答変数があるときに使用する適切なタイプの検定、Pearson Correlationを見ましょう。
Pearson 相関検定は、2つの提供された変数（両方とも量的変数）の間の関係の強さを分析するために使用されます。ピアソン相関の値、または強さは、+1 から -1 の間になります。
相関が1の場合は、変数間の完全な関連性を示し、相関は正または負のいずれかである。0に近い相関係数は、非常に弱い、ほとんど存在しない相関を示します。Spearman相関やKendall Rank相関など、2つの変数間の相関を測定する他の方法もありますが、Pearson相関はおそらく最もよく使われる相関テストでしょう。
Gapminderデータセットは、その特徴を量的変数で表現しているので、Pearson Correlationを実行する前に、データのカテゴリー変換をする必要はありません。両方の変数が正規分布で、データセットに有意な外れ値があまりないことを仮定していることに注意してください。ピアソン相関を実行するためには、SciPyにアクセスする必要があります。
平均寿命とインターネット利用率、インターネット利用率と就業率の関係をグラフ化し、別の相関グラフがどのように見えるかを確認します。グラフ作成後、SciPyの関数 personr() を使って相関をとり、その結果を確認します。
df_clean = df2.dropna()
df_clean['incomeperperson'] = df_clean['incomeperperson'].replace('', np.nan)
def plt_regression(x, y, data, label_1, label_2):
reg_plot = regplot(x=x, y=y, fit_reg=True, data=data)
plt.xlabel(label_1)
plt.ylabel(label_2)
plt.show()
plt_regression('lifeexpectancy', 'internetuserate', df_clean, 'Life Expectancy', 'Internet Use Rate')
plt_regression('employrate', 'internetuserate', df_clean, 'Employment Rate', 'Internet Use Rate')
print('Assoc. - life expectancy and internet use rate')
print(pearsonr(df_clean['lifeexpectancy'], df_clean['internetuserate']))
print('Assoc. - between employment rate and internet use rate')
print(pearsonr(df_clean['employrate'], df_clean['internetuserate']))
以下はその出力です。
Assoc. - life expectancy and internet use rate
(0.77081050888289, 5.983388253650836e-33)
Assoc. - between employment rate and internet use rate
(-0.1950109538173115, 0.013175901971555317)
最初の値は相関の方向と強さ、2番目はP値です。この数値は、平均寿命とインターネット利用率の間に、偶然ではない、かなり強い相関があることを示唆しています。一方、就業率とインターネット利用率の間には、弱いながらも有意な相関があることがわかる。
なお、カテゴリーデータに対してピアソン相関を行うことも可能ですが、結果は多少異なります。その気になれば、所得水準をグループ化してピアソン相関を実行することも可能です。興味のある関連性に影響を与える可能性のある緩和変数の存在を確認するために使用することができます。
モデレータと統計的相互作用
ここでは、複数の変数間の統計的相互作用（別名：モデレーション）を説明する方法を見てみましょう。
モデレーションとは、第3の（またはそれ以上の）変数が、独立変数と従属変数の間の関連性の強さに影響を与えることです。
第3の変数と独立/従属変数の間のモデレーション/統計的交互作用を検定するさまざまな方法があります。たとえば、ANOVA検定を実施した場合、可能な限りのモデレーションを検定するために、二元配置ANOVA検定を行うことでモデレーションを検定することができます。
しかし、どのような種類の統計検定（ANOVA、カイ2乗、ピアソン相関）を実行しても、モデレーションを検定する信頼できる方法は、第3の変数のサブグループ/レベルごとに説明変数と応答変数の間に関連性があるかどうかを確認することです。
より具体的に言うと、もしANOVA検定を実施しているなら、3番目の変数（あなたが調査している関係性にモデレート効果を持つかもしれないと疑う変数）のすべてのカテゴリについてANOVAを実行すればよいのです。
もし、カイ2乗検定を使用しているなら、モデレーティング変数のカテゴリ内で見つかったすべてのデータポイントを保持している新しいデータフレームで、カイ2乗検定を実行すればよいでしょう。
統計的検定がPearson相関であれば、モデレーティング変数のカテゴリまたはビンを作成して、それらのビンの3つすべてについてPearson相関を実行する必要があります。
それでは、モデレーティング変数についてPearson相関を実行する方法を簡単に見てみましょう。我々は、連続特徴量から人工的なカテゴリ/レベルを作成します。他の2つのテストタイプ（カイ2乗とANOVA）のモデレーションを検定するプロセスは、とても似ていますが、その代わりに既存のカテゴリ変数があります。
我々は、我々のモデレーティング変数として作用する適切な変数を選びたいのです。1人あたりの所得レベルを試して、3つの異なるグループに分けてみましょう。
def income_groups(row):
if row['incomeperperson'] <= 744.23:
return 1
elif row['incomeperperson'] <= 942.32:
return 2
else:
return 3
# Apply function and set the new features in the dataframe
df_clean['income_group'] = df_clean.apply(lambda row: income_groups(row), axis=1)
# Create a few subframes to try test for moderation
subframe_1 = df_clean[(df_clean['income_group'] == 1)]
subframe_2 = df_clean[(df_clean['income_group'] == 2)]
subframe_3 = df_clean[(df_clean['income_group'] == 3)]
print('Assoc. - life expectancy and internet use rate for low income countries')
print(pearsonr(subframe_1['lifeexpectancy'], subframe_1['internetuserate']))
print('Assoc. - life expectancy and internet use rate for medium income countries')
print(pearsonr(subframe_2['lifeexpectancy'], subframe_2['internetuserate']))
print('Assoc. - life expectancy and internet use rate for high income countries')
print(pearsonr(subframe_3['lifeexpectancy'], subframe_3['internetuserate']))
以下はその出力です。
Assoc. - life expectancy and internet use rate for low income countries
(0.38386370068495235, 0.010101223355274047)
Assoc. - life expectancy and internet use rate for medium income countries
(0.9966009508278395, 0.05250454954743393)
Assoc. - life expectancy and internet use rate for high income countries
(0.7019997488251704, 6.526819886007788e-18)
もう一度言いますが、最初の値は相関の方向と強さで、2番目はP値です。
結論
statsmodels` はPythonユーザーがデータを分析し、データセットに対して統計的検定を実行できるようにする非常に便利なライブラリである。ANOVA、カイ二乗検定、ピアソン相関、そしてモデレーションの検定を行うことができます。
これらの検定を実行する方法に慣れると、変数のカテゴリや連続の性質に適応して、従属変数と独立変数の間の有意な関係を検定することができるようになります。