回帰分析II：重回帰分析の方法、科学的な将来予測：ITエンジニアのためのデータサイエンティスト養成講座（7）（1/3 ページ）
単純な回帰分析によるデータ分析ができたところで、今後は少し複雑な重回帰分析にチャレンジしてみましょう。モデル評価の方法や将来予測に役立つ考え方を交えて紹介します。
はじめに
前回の記事「回帰分析I：回帰分析って何？ から、最小二乗法、モデル評価、妥当性検討の実際まで」は回帰分析について、解法で使われる最小二乗法やモデルの評価、妥当性の検証など基本的なことをご紹介しました。今回は応用編ということで、説明変数が2つ以上のより複雑な重回帰分析や回帰分析を“意思決定やアクションにつなげる”ために、モデルを予測に活用する方法について詳しく紹介していきます。
再び、電力消費量と気温との相関関係
まずは第5回の後半でご紹介した電力消費量と気温との相関関係についてあらためて見ていきたいと思います。データはいつも利用している電力の使用状況データリンクと気象庁が提供している「過去の気象データ・ダウンロード」リンクを使います。過去の気象データは、東京地方の最高気温と最低気温を2013年1月1日から2013年11月15日までの日別値であらかじめCSV形式で“data.csv”というファイル名でダウンロードして利用します。
以下がダウンロードしたCSVデータです。
ダウンロードした時刻：2013/11/17 09:36:40 ,東京,東京 年月日,最高気温(℃),最低気温(℃) 2013/1/1,9.6,3.0 2013/1/2,14.4,4.2 2013/1/3,9.9,1.8 2013/1/4,6.9,1.0 2013/1/5,4.8,0.2 （途中省略） 2013/11/11,16.8,8.3 2013/11/12,12.5,7.9 2013/11/13,14.1,6.4 2013/11/14,14.5,6.2 2013/11/15,12.7,10.1
データの取り込みと変換
まずは電力資料量と気温のデータを取り込んで、分析できる状態に変換します。電力使用量の方は1時間ごとのデータになっていますのでダウンサンプリングして1日ごとのデータに変換します。（In [4]〜In [13]）また、それぞれのデータの期間を5月から9月までの平日のデータに絞り込んでおきます。このあたりの詳細は第5回を参考にしてください。
In [1]: # ライブラリのインポート
In [2]: import pandas as pd
In [3]: import numpy as np
In [4]: # 電力の使用状況データの取り込み
In [5]: df_power = pd.read_csv('http://www.tepco.co.jp/forecast/html/images/juyo-2013.csv', skiprows=3, names=['date', 'time', 'actual'], encoding='Shift_JIS')
In [6]: # 時系列インデックスの生成と付与
In [7]: idx_power = pd.to_datetime(df_power['date']+' '+df_power['time'])
In [8]: df_power.index = idx_power
In [9]: # 必要ない列の削除
In [10]: del df_power['date']
In [11]: del df_power['time']
In [12]: # ダウンサンプリング（最大値）
In [13]: df_power_daily = df_power.resample('D', how='max', kind='period')
In [14]: # 東京都の最高気温と最低気温データの取り込み
In [15]: df_temp = pd.read_csv('data.csv', skiprows=4, names=['date', 'high', 'low'], encoding='Shift_JIS')
In [16]: # 時系列インデックスの生成と付与
In [17]: idx_temp = pd.to_datetime(df_temp['date'])
In [18]: df_temp.index = idx_temp
In [19]: # 必要ない列の削除
In [20]: del df_temp['date']
In [21]: # 5月から9月までの平日のデータに絞り込み
In [22]: dfx = df_temp.ix['2013-05-01':'2013-09-30'].resample('B', how='max', kind='period')
In [23]: dfy = df_power_daily.ix['2013-05-01':'2013-09-30'].resample('B', how='max', kind='period')
電力消費量と最高気温の相関関係を単回帰分析する
単回帰分析ではy＝ ax ＋ bというモデルを仮定して係数と定数を求めてモデルを完成させます。
今回も「1日の電力使用量の最大値は最高気温と相関関係がある」という仮説を立て、モデル（1）を仮定します。
［電力消費量］＝［最高気温］×［係数a］＋［定数b］ ――モデル（1）
今回も前回と同様に最小二乗法（OLS：Ordinary Least Squares）を使って［係数a］と［定数b］を求めます。前回と異なり、今回はPandasの関数ではなくStatsModels（リンク）で提供されているAPIを使って、モデル（1）の係数と定数を求めます。
まずは、必要なライブラリをインポートし（In [24]〜In [28]）、説明変数Xと目的変数yを定義します（In [29]〜In [35]）。前回は"intercept=True"というオプションで［定数b］（切片）が0ではない値を取ると定義していましたが、今回はStatsModelsで提供されているadd_constant関数を用いてXの中に定数を持つことを指定します（In [32]）。
変数の準備ができたらモデルを定義して最小二乗法を実行します（In [36]〜In [41]）
In [24]: # ライブラリのインポート
In [25]: %pylab
Welcome to pylab, a matplotlib-based Python environment [backend: TkAgg].
For more information, type 'help(pylab)'.
In [26]: from statsmodels.sandbox.regression.predstd import wls_prediction_std
In [27]: import statsmodels.api as sm
In [28]: import matplotlib.pyplot as plt
In [29]: # 最高気温を説明変数Xとして定義
In [30]: X = dfx['high']
In [31]: # 定数項をXに加える
In [32]: X = sm.add_constant(X, prepend=False)
In [33]: # 電力消費量を目的変数yとして定義
In [34]: y_true = dfy['actual']
In [35]: y = y_true
In [36]: # モデルを定義
In [37]: model = sm.OLS(y,X)
In [38]: # 回帰分析を実行
In [39]: results = model.fit()
In [40]: # 結果のサマリを確認
In [41]: results.summary()
Out[41]:
<class 'statsmodels.iolib.summary.Summary'>
"""
OLS Regression Results
==============================================================================
Dep. Variable: actual R-squared: 0.777
Model: OLS Adj. R-squared: 0.775
Method: Least Squares F-statistic: 372.3
Date: Sun, 17 Nov 2013 Prob (F-statistic): 1.25e-36
Time: 15:14:59 Log-Likelihood: -753.85
No. Observations: 109 AIC: 1512.
Df Residuals: 107 BIC: 1517.
Df Model: 1
==============================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
high 110.4039 5.721 19.296 0.000 99.062 121.746
const 765.1262 167.654 4.564 0.000 432.772 1097.480
==============================================================================
Omnibus: 3.225 Durbin-Watson: 1.354
Prob(Omnibus): 0.199 Jarque-Bera (JB): 2.547
Skew: -0.243 Prob(JB): 0.280
Kurtosis: 2.430 Cond. No. 209.
==============================================================================
"""
結果を確認すると、モデル（1）の係数a（high）と定数b（const）はそれぞれ、110.4039と765.1262となることが分かります。よって、モデル（1）は以下の数式で表すことができます。
［電力消費量］＝［最高気温］×110.4039＋765.1262
それでは、前回紹介した3つの観点からモデルの妥当性を評価してみましょう。
（1）決定係数（R-squared: 0.777）
決定係数は77％以上であり、よく当てはまっているといえます。
（2）F検定（F-stat: 372.3, p-value: 1.25e-36）
F値が十分に大きく、有意確率（p-value）も十分に小さいので、モデル全体として妥当であると判断できます。
（3）t検定（係数 high:t-stat=19.296、p-value=0.000、定数 const:t-stat=4.564、p-value=0.000）
係数high、定数constの有意確率双方とも十分に小さい値であり係数、定数とも妥当であると判断できます。
これらの結果から「1日の電力使用量の最大値は最高気温と相関関係がある」という仮説に基づくモデル（1）は、かなり良好なモデルであるといえます。
Copyright © ITmedia, Inc. All Rights Reserved.