statsmodels

statsmodels

statsmodelsは、統計モデルの推定や統計検定、データの可視化などを行うためのPythonライブラリです。

以下に、statsmodelsの便利な使い方をいくつか紹介します。

1. 線形回帰モデル

statsmodelsを使用して線形回帰モデルを構築し、データの関係性を解析できます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import statsmodels.api as sm
import pandas as pd

# データの準備
df = pd.DataFrame({
    'X': [1, 2, 3, 4, 5],
    'Y': [1, 2, 1.3, 3.75, 2.25]
})

X = df['X']
Y = df['Y']

# 定数項の追加
X = sm.add_constant(X)

# 回帰モデルの構築とフィッティング
model = sm.OLS(Y, X).fit()

# モデルの概要
print(model.summary())

[実行結果]

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      Y   R-squared:                       0.393
Model:                            OLS   Adj. R-squared:                  0.191
Method:                 Least Squares   F-statistic:                     1.942
Date:                Sun, 07 Jul 2024   Prob (F-statistic):              0.258
Time:                        02:26:29   Log-Likelihood:                -5.6369
No. Observations:                   5   AIC:                             15.27
Df Residuals:                       3   BIC:                             14.49
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.7850      1.012      0.776      0.494      -2.434       4.004
X              0.4250      0.305      1.393      0.258      -0.546       1.396
==============================================================================
Omnibus:                          nan   Durbin-Watson:                   3.369
Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.554
Skew:                           0.641   Prob(JB):                        0.758
Kurtosis:                       1.993   Cond. No.                         8.37
==============================================================================

2. ロジスティック回帰

カテゴリカルデータに対してロジスティック回帰を適用する例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import statsmodels.api as sm
import numpy as np

# データの準備
np.random.seed(42)
df = pd.DataFrame({
    'X1': np.random.rand(100),
    'X2': np.random.rand(100),
    'Y': np.random.randint(0, 2, 100)
})

X = df[['X1', 'X2']]
Y = df['Y']

# 定数項の追加
X = sm.add_constant(X)

# ロジスティック回帰モデルの構築とフィッティング
model = sm.Logit(Y, X).fit()

# モデルの概要
print(model.summary())

[実行結果]

Optimization terminated successfully.
         Current function value: 0.677050
         Iterations 4
                           Logit Regression Results                           
==============================================================================
Dep. Variable:                      Y   No. Observations:                  100
Model:                          Logit   Df Residuals:                       97
Method:                           MLE   Df Model:                            2
Date:                Sun, 07 Jul 2024   Pseudo R-squ.:                 0.01611
Time:                        02:27:12   Log-Likelihood:                -67.705
converged:                       True   LL-Null:                       -68.814
Covariance Type:            nonrobust   LLR p-value:                    0.3299
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.4792      0.524      0.915      0.360      -0.547       1.506
X1             0.4080      0.688      0.593      0.553      -0.939       1.756
X2            -0.9362      0.701     -1.336      0.181      -2.309       0.437

3. 時系列解析

statsmodelsを使用してARIMAモデルなどの時系列解析を行うことができます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import statsmodels.api as sm
import pandas as pd

# 時系列データの準備
data = sm.datasets.co2.load_pandas().data
data = data.fillna(data.bfill())

# ARIMAモデルのフィッティング
model = sm.tsa.ARIMA(data['co2'], order=(1, 1, 1))
result = model.fit()

# モデルの概要
print(result.summary())

# 予測
forecast = result.forecast(steps=10)
print(forecast)

[実行結果]

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                    co2   No. Observations:                 2284
Model:                 ARIMA(1, 1, 1)   Log Likelihood               -1536.855
Date:                Sun, 07 Jul 2024   AIC                           3079.710
Time:                        02:27:49   BIC                           3096.910
Sample:                    03-29-1958   HQIC                          3085.984
                         - 12-29-2001                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.8959      0.022     40.360      0.000       0.852       0.939
ma.L1         -0.7587      0.034    -22.509      0.000      -0.825      -0.693
sigma2         0.2250      0.006     39.341      0.000       0.214       0.236
===================================================================================
Ljung-Box (L1) (Q):                  55.87   Jarque-Bera (JB):                74.67
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               1.13   Skew:                            -0.14
Prob(H) (two-sided):                  0.10   Kurtosis:                         3.84
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
2002-01-05    371.652800
2002-01-12    371.789687
2002-01-19    371.912318
2002-01-26    372.022178
2002-02-02    372.120597
2002-02-09    372.208766
2002-02-16    372.287753
2002-02-23    372.358514
2002-03-02    372.421906
2002-03-09    372.478695
Freq: W-SAT, Name: predicted_mean, dtype: float64

4. ANOVA（分散分析）

複数グループ間の平均値の差異を検定するためにANOVAを使用します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import statsmodels.api as sm
from statsmodels.formula.api import ols

# データの準備
df = sm.datasets.get_rdataset("iris").data

# 列名の変更
df.columns = [col.replace('.', '_') for col in df.columns]
print(df.columns)  # 変更後の列名の確認

# モデルの構築
model = ols('Sepal_Length ~ C(Species)', data=df).fit()

# ANOVAテーブル
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)

[実行結果]

Index(['Sepal_Length', 'Sepal_Width', 'Petal_Length', 'Petal_Width',
       'Species'],
      dtype='object')
               sum_sq     df           F        PR(>F)
C(Species)  63.212133    2.0  119.264502  1.669669e-31
Residual    38.956200  147.0         NaN           NaN

5. 残差プロット

モデルの適合度を視覚的に評価するために、残差プロットを作成できます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import statsmodels.api as sm
from statsmodels.formula.api import ols
import matplotlib.pyplot as plt

# データの準備
df = sm.datasets.get_rdataset("iris").data
df.columns = [col.replace('.', '_') for col in df.columns]

# モデルの構築（Sepal_LengthとPetal_Lengthの線形回帰）
model = ols('Sepal_Length ~ Petal_Length', data=df).fit()

# 残差のプロット
fig, ax = plt.subplots()
sm.graphics.plot_regress_exog(model, 'Petal_Length', fig=fig)
plt.show()

[実行結果]

6. 様々な統計検定

t検定やカイ二乗検定など、様々な統計検定を行うことができます。

t検定

1
2
3
4
5
6
7
8
9
10
from scipy import stats
import statsmodels.api as sm

# データの準備
df = sm.datasets.get_rdataset("iris").data
df.columns = [col.replace('.', '_') for col in df.columns]

# 2つの独立した標本のt検定
t_stat, p_value = stats.ttest_ind(df['Sepal_Length'], df['Petal_Length'])
print(f'T-statistic: {t_stat}, P-value: {p_value}')

[実行結果]

T-statistic: 13.098353108960858, P-value: 2.857104069581941e-31

カイ二乗検定

1
2
3
4
5
6
7
8
9
10
11
12
import pandas as pd
from scipy.stats import chi2_contingency
import statsmodels.api as sm

# データの準備
df = sm.datasets.get_rdataset("iris").data
df.columns = [col.replace('.', '_') for col in df.columns]

# カイ二乗検定
table = pd.crosstab(df['Species'], df['Sepal_Length'] > df['Sepal_Length'].median())
chi2, p, dof, ex = chi2_contingency(table)
print(f'Chi2: {chi2}, P-value: {p}')

[実行結果]

Chi2: 78.64285714285712, P-value: 8.37376079899524e-18

これらの例は、statsmodelsの基本的な使用方法を示していますが、このライブラリは他にも多くの機能を提供しており、統計解析やモデリングの幅広いニーズに応えることができます。