statsmodels 主成分分析¶

主要概念： 主成分分析、世界銀行數據、生育率

在本筆記本中，我們使用主成分分析 (PCA) 來分析 192 個國家生育率的時間序列，數據來自世界銀行。主要目標是了解各國生育率隨時間變化的趨勢有何不同。這是一個稍微不典型的 PCA 應用，因為數據是時間序列。雖然已經開發了用於此設定的方法（例如函數型 PCA），但由於生育率數據非常平滑，因此在這種情況下使用標準 PCA 並沒有真正的缺點。

[1]:

%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)

數據可以從世界銀行網站取得，但這裡我們使用稍微清理過的數據版本

[2]:

data = sm.datasets.fertility.load_pandas().data
data.head()

[2]:

在這裡，我們建立一個 DataFrame，其中只包含數值的生育率數據，並將索引設定為國家名稱。我們也刪除所有有任何遺失數據的國家。

[3]:

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

[3]:

使用 PCA 分析矩形矩陣有兩種方法：我們可以將列視為「對象」，將欄視為「變數」，反之亦然。在這裡，我們將把生育率衡量標準視為衡量國家「對象」的「變數」。因此，目標是將每年的生育率值縮減為少量的生育率「輪廓」或「基底函數」，這些輪廓或基底函數能捕捉不同國家隨著時間推移的大部分變化。

PCA 中會移除平均趨勢，但還是值得看一下。它顯示在該數據集中涵蓋的時間段內，生育率穩步下降。請注意，平均值是使用國家作為分析單位計算的，忽略了人口規模。以下進行的 PC 分析也是如此。更複雜的分析可能會對國家進行加權，例如以 1980 年的人口為權重。

[4]:

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

[4]:

(0.0, 51.0)

接下來，我們執行 PCA

[5]:

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

根據特徵值，我們看到第一個 PC 佔據主導地位，第二和第三個 PC 可能捕獲了少量的有意義變化。

[6]:

fig = pca_model.plot_scree(log_scale=False)

接下來，我們將繪製 PC 因子。主導因子是單調遞增的。在第一個因子中得分為正的國家，其增長速度會比上面顯示的平均值更快（或減少速度更慢）。在第一個因子中得分為負的國家，其減少速度會比平均值更快。第二個因子呈 U 形，在 1985 年左右達到正峰值。在第二個因子中得分為正值較大的國家，在數據範圍的開始和結束時的生育率會低於平均水平，但在範圍中間的生育率會高於平均水平。

[7]:

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)

/tmp/ipykernel_3281/427128218.py:3: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  ax.set_xticklabels(dta.columns.values[::10])

為了更好地了解情況，我們將繪製具有相似 PC 分數的國家組的生育率軌跡。以下便利函數會產生這樣的圖表。

[8]:

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

首先，我們繪製在 PC 1 上得分最高的五個國家。這些國家的生育率增長率高於全球平均值（正在下降）。

[9]:

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)

[10]:

labels = dta.index[idx[-5:]]
make_plot(labels)

以下是在因子 2 上得分最高的五個國家。這些國家的生育率在 1980 年左右達到頂峰，比世界其他大部分地區晚，然後生育率迅速下降。

[11]:

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

最後，我們有在 PC 2 上得分最負的國家。這些國家在 1960 年代和 1970 年代的生育率下降速度比全球平均值快得多，然後趨於平穩。

[12]:

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

我們也可以查看前兩個主成分分數的散點圖。我們看到國家之間的變化相當連續，但 PC 2 得分最高的兩個國家可能與其他點有些分離。這兩個國家，阿曼和葉門，在 1980 年左右的生育率出現了明顯的尖峰，這是獨特的。沒有其他國家有這樣的尖峰。相反，在 PC 1 上得分較高（生育率持續上升）的國家是變化的連續體的一部分。

[13]:

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

[13]:

array(['Oman', 'Yemen, Rep.'], dtype=object)