statsmodels 統計

主成分分析

 このノートブックでは、私たちは、世界銀行からデータを入手して、192か国の出生率の時系列データを分析するために、主成分分析を使います。主なゴールは、国と国の間で異なる時間経過の出生率のトレンドがどれくらいか理解することです。これは、データが時系列のために、少々、不規則なPCAの解説です。機能的PCAのような方法は、この設定のために開発されました。しかし、出生率のデータはとてもスムースなので、このケースで標準的なPCAを使うことに、実際に不便な点はありません。

%matplotlib inline

import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.multivariate.pca import PCA

plt.rc("figure", figsize=(16, 8))
plt.rc("font", size=14)

 データは、世界現行のサイトにから取得することができます。しかし、ここで、私たちは少しデータを修正したバージョンを使います。

data = sm.datasets.fertility.load_pandas().data
data.head()
Country NameCountry CodeIndicator NameIndicator Code196019611962196319641965...2004200520062007200820092010201120122013
0ArubaABWFertility rate, total (births per woman)SP.DYN.TFRT.IN4.8204.6554.4714.2714.0593.842...1.7861.7691.7541.7391.7261.7131.7011.690NaNNaN
1AndorraANDFertility rate, total (births per woman)SP.DYN.TFRT.INNaNNaNNaNNaNNaNNaN...NaNNaN1.2401.1801.2501.1901.220NaNNaNNaN
2AfghanistanAFGFertility rate, total (births per woman)SP.DYN.TFRT.IN7.6717.6717.6717.6717.6717.671...7.1366.9306.7026.4566.1965.9285.6595.395NaNNaN
3AngolaAGOFertility rate, total (births per woman)SP.DYN.TFRT.IN7.3167.3547.3857.4107.4257.430...6.7046.6576.5986.5236.4346.3316.2186.099NaNNaN
4AlbaniaALBFertility rate, total (births per woman)SP.DYN.TFRT.IN6.1866.0765.9565.8335.7115.594...2.0041.9191.8491.7961.7611.7441.7411.748NaNNaN

5 rows × 58 columns

 


ここで、データフレームを構築します。それは、出生率の数値データだけを含んでおり、国名にインデックスを設定します。

columns = list(map(str, range(1960, 2012)))
data.set_index("Country Name", inplace=True)
dta = data[columns]
dta = dta.dropna()
dta.head()

1960196119621963196419651966196719681969...2002200320042005200620072008200920102011
Country Name
Aruba4.8204.6554.4714.2714.0593.8423.6253.4173.2263.054...1.8251.8051.7861.7691.7541.7391.7261.7131.7011.690
Afghanistan7.6717.6717.6717.6717.6717.6717.6717.6717.6717.671...7.4847.3217.1366.9306.7026.4566.1965.9285.6595.395
Angola7.3167.3547.3857.4107.4257.4307.4227.4037.3757.339...6.7786.7436.7046.6576.5986.5236.4346.3316.2186.099
Albania6.1866.0765.9565.8335.7115.5945.4835.3765.2685.160...2.1952.0972.0041.9191.8491.7961.7611.7441.7411.748
United Arab Emirates6.9286.9106.8936.8776.8616.8416.8166.7836.7386.679...2.4282.3292.2362.1492.0712.0041.9481.9031.8681.841

5 rows × 52 columns

 対角の行列を分析するためにPCAを使うふたつの方法があります。私たちは、縦列をオブジェクト、横列を変数、逆もまた同様に、整えることができます。ここで私たちは、出生率の尺度を変数として整え、国の尺度をオブジェクトとして使います。このようにゴールは、年間の出生率の値を、異なる国の時間を超えた多くの種類を取り込んだ少数の出生率"プロファイル"または"基礎機能"に減少させることです。 

 平均のトレンドは、PCAから除外されますが、それを観測することは役に立ちます。出生率は、このデータでカバーする時間経過に応じてしっかりと減少していくことを示しています。平均値は、分析ユニットとして、人口の大きさを無視した国を使って計算していることに注意してください。これもまた、以下で導出される主成分分析の真実です。もっと洗練された分析は、例えば1980年の人口によって、国を重みづけするかもしれません。

ax = dta.mean().plot(grid=False)
ax.set_xlabel("Year", size=17)
ax.set_ylabel("Fertility rate", size=17)
ax.set_xlim(0, 51)

次に私たちは、PCAを実行します。

pca_model = PCA(dta.T, standardize=False, demean=True)

 固有値を基礎にして、私たちは、2番目と3番目の主成分を取り込んだ、おそらく少量の有意義な(意味のある)変異とともに、最初の強い影響のある主成分を見ます。

fig = pca_model.plot_scree(log_scale=False)

 次に、私たちは、主成分因子を図示します。強い影響のある因子は、一本調子で増加します。最初の因子が正のスコアの国々は、上に示す平均と比較してより速く増加(またはより遅く減少)します。最初の因子が負のスコアの国々は、平均よりもっと速く減少します。2番目の因子は、1985年の周辺でせいのピークを持つU字型をしています。2番目の因子が大きな正の値を持つ国々は、最初から終端までのデータの範囲の平均出生率より低くなります。しかし、中間の範囲の平均出生率より高くなります。

fig, ax = plt.subplots(figsize=(8, 4))
lines = ax.plot(pca_model.factors.iloc[:, :3], lw=4, alpha=0.6)
ax.set_xticklabels(dta.columns.values[::10])
ax.set_xlim(0, 51)
ax.set_xlabel("Year", size=17)
fig.subplots_adjust(0.1, 0.1, 0.85, 0.9)
legend = fig.legend(lines, ["PC 1", "PC 2", "PC 3"], loc="center right")
legend.draw_frame(False)

 何が起きているかよりよく理解するために、同様の主成分スコアのくのセットのために出生率の足跡を、私たちは図示 します。下の便利な関数は、そうした図を提供します。

idx = pca_model.loadings.iloc[:, 0].argsort()

 最初に、主成分1の最も大きいスコアを持つ五つの国を図示します。これらの国は、全体の平均より高い出生率で増加します。

def make_plot(labels):
    fig, ax = plt.subplots(figsize=(9, 5))
    ax = dta.loc[labels].T.plot(legend=False, grid=False, ax=ax)
    dta.mean().plot(ax=ax, grid=False, label="Mean")
    ax.set_xlim(0, 51)
    fig.subplots_adjust(0.1, 0.1, 0.75, 0.9)
    ax.set_xlabel("Year", size=17)
    ax.set_ylabel("Fertility", size=17)
    legend = ax.legend(
        *ax.get_legend_handles_labels(), loc="center left", bbox_to_anchor=(1, 0.5)
    )
    legend.draw_frame(False)
labels = dta.index[idx[-5:]]
make_plot(labels)

 ここは、因子2の最も大きいスコアを持つ五つの国です。残りの多くの世界よりも遅く、1980年の周辺で出生率がピークに達する国です。その後急速に出生率が減少します。

idx = pca_model.loadings.iloc[:, 1].argsort()
make_plot(dta.index[idx[-5:]])

 最後に、因子2のスコアが最も大きな負を持つ国です。これらは、出生率が1960年代と1970年代の間で全体の平均よりもっと速く減少します。その後、平坦になります。

make_plot(dta.index[idx[:5]])

 私たちは、最初の二つに主成分のスコアのスキャッタープロットを見ることができます。私たちは、国の間の多様さが、結構連続的で、おそらく、他の点から何かが分離される主成分2で最も高いスコアを持つふたつの国をのぞきます。これらの国、オマーンとイエメン、は、1980年の周りで出生率に鋭い頂点を持っており、特別です。他の国はそうした頂点は持っていません。対照的に、因子1が高いスコアの国は、(持続して出生率が増加しています。)は、変化が連続している部分です。

fig, ax = plt.subplots()
pca_model.loadings.plot.scatter(x="comp_00", y="comp_01", ax=ax)
ax.set_xlabel("PC 1", size=17)
ax.set_ylabel("PC 2", size=17)
dta.index[pca_model.loadings.iloc[:, 1] > 0.2].values

-statsmodels, 統計
-