Edited at

PythonでstepAIC

More than 1 year has passed since last update.


概要

えっPythonにstepAICないの.

…と思ってたらStackOverflowにこんな記事が.

回答者が紹介しているリンク (Forward Selection with statsmodels; しかし線形回帰のみに対応し,また指標はAICではなく決定係数) を参考にstepAICを書くことにしました.


step_aic

def step_aic(model, exog, endog, **kwargs):

"""
This select the best exogenous variables with AIC
Both exog and endog values can be either str or list.
(Endog list is for the Binomial family.)

Note: This adopt only "forward" selection

Args:
model: model from statsmodels.formula.api
exog (str or list): exogenous variables
endog (str or list): endogenous variables
kwargs: extra keyword argments for model (e.g., data, family)

Returns:
model: a model that seems to have the smallest AIC
"""

# exog, endogは強制的にリスト形式に変換しておく
exog = np.r_[[exog]].flatten()
endog = np.r_[[endog]].flatten()
remaining = set(exog)
selected = [] # 採用が確定された要因

# 定数項のみのAICを計算
formula_head = ' + '.join(endog) + ' ~ '
formula = formula_head + '1'
aic = model(formula=formula, **kwargs).fit().aic
print('AIC: {}, formula: {}'.format(round(aic, 3), formula))

current_score, best_new_score = np.ones(2) * aic

# 全要因を採択するか,どの要因を追加してもAICが上がらなければ終了
while remaining and current_score == best_new_score:
scores_with_candidates = []
for candidate in remaining:

# 残っている要因を1つずつ追加したときのAICを計算
formula_tail = ' + '.join(selected + [candidate])
formula = formula_head + formula_tail
aic = model(formula=formula, **kwargs).fit().aic
print('AIC: {}, formula: {}'.format(round(aic, 3), formula))

scores_with_candidates.append((aic, candidate))

# 最もAICが小さかった要因をbest_candidateとする
scores_with_candidates.sort()
scores_with_candidates.reverse()
best_new_score, best_candidate = scores_with_candidates.pop()

# 候補要因追加でAICが下がったならば,それを確定要因として追加する
if best_new_score < current_score:
remaining.remove(best_candidate)
selected.append(best_candidate)
current_score = best_new_score

formula = formula_head + ' + '.join(selected)
print('The best formula: {}'.format(formula))
return model(formula, **kwargs).fit()


使い方

説明変数がxとfだった場合はこんな感じ.


(['y']じゃなくて'y'でもいけます)

20170225.png


あとがき

合ってるか…合ってるよな…?

みどりぼん (データ解析のための統計モデリング入門) の二項分布やロジスティック回帰の章とは回答が完全一致しました.

でも間違ってたら教えてください.