statsmodelsで回帰直線を作る

Python

こんにちは!ウメハラ(plumfield56)です。

この記事ではPythonのstatsmodelsライブラリの使い方を説明していきます。
states(統計)models(モデル)のことで統計データを出したり統計検定するのに活用できます。
公式サイトはこちらです。

Introduction — statsmodels

データの用意

今回はKaggleにある住宅情報を利用して、家賃と部屋の広さとの回帰分析を最小二乗法で書いていきたいと思います。
下記のURLからデータをダウンロードしてください。

kc_house_data
Kaggle is the world’s largest data science community with powerful tools and resources to help you achieve your data science goals.

ダウンロードは下記の箇所からできます。

ダウンロードできるCSVファイルはこちらになるのでここからダウンロードでも問題ないです。

データ内容の確認

カラムの確認

まずはデータを取り込んで、どのようなカラムがあるのか確認するところから始めていきましょう。

import pandas as pd
data = pd.read_csv('kc_house_data.csv')
data.head()

カラム一覧を見てみると下記のようになっています。

英語カラム名日本語
idID
date販売日
price 価格($)
bedroomsベッドルームの数
bathroomsバスルームの数
sqft_livingリビングルームの広さ
sqft_lot駐車場の広さ
floors階数
waterfrontウォーターフロントか (0,1)
viewビューの種類(0,1,2,3,4)
condition状態(1,2,3,4,5)
gradeグレード(1~13)
sqft_above地上階の広さ
sqft_basement地下の広さ
yr_built建築年
yr_renovatedリノベーションの年
zipcode 郵便番号
lat緯度
long経度
sqft_living15 近所15件の平均のリビングの広さ
sqft_lot15近所15件の駐車場の広さ

この中で今回はpriceとsqft_aboveを使用していきます。

欠損値の確認

欠損地は下記のどちらかのコードで確認できます。

data.info()
data.isnull().sum()

データを確認するとこれから使おうと思っているsqft_aboveに欠損地があるので、何かしら処理をしていきます。
処理の方法としては主に下記の2つがあります。

  • 欠損値の削除
  • 特定の値で埋める(0や平均値、中央値、1つ上の値)

今回は削除して進めていきたいと思います。
dropna()で空白行を削除することができます。

data = data.dropna()

実行後にちゃんと削除して欠損値がなくなっているか再度確認してください。

散布図の表示

次に大まかにどのようにプロットされているのか散布図で確認していきます。
matplotlibで書いていきますが、個人的にはseabornの表示のほうがきれいで好きなので表示はseabornになるようにしていきます。
どうするのかというと下記のようにseabornでsetをしてあげます。

import matplotlib.pylab as plt
import seaborn as sns
sns.set()

この後に下記のように記載すると散布図が作成できます。

y = data['price']
x1 = data['sqft_above']

plt.scatter(x1, y)
plt.xlabel('sqft_above')
plt.ylabel('price')
plt.show()

すると下記の散布図が描かれます。
seabornをセットした場合としない場合をどちらも載せておきます。

この散布図をみるとある程度相関はありそうなのでここから本題のstatsmodelsを使って統計データを出して回帰線を引いていきたいと思います。

statsmodelsで統計データを表示

statsmodelsで統計データを出す手順としては下記になります。

  1. 目的変数、説明変数の設定
  2. 分析手法の設定
  3. 表示させる

目的変数、説明変数の設定

statsmodelsで統計データを出すためにはxとyになる値を渡すことで出すことができます。
その際に説明変数となるデータはadd_constantで設定する必要があります。
今回は単回帰分析ですが、複数ある場合は複数データ渡すことも可能です。
ちなみにエラーがでたときにendog,exogという言葉がでてくるので一覧で日本語とセットで英語も確認ください。

日本語1日本語2statsmodels
での英語表記
元の英語英語の意味
目的変数従属変数endogendogenous内因性の
説明変数独立変数exogexogenous外因性の

下記公式の$y$が目的変数、$x_1$が説明変数になります。

$$ y= β_0+ β_1 x_1  + ε $$

分析手法の設定

分析手法には下記のような種類があります。

日本語英語
最小二乗法Ordinary Least SquaresOLS
加重最小二乗法Weighted Least SquaresWLS
一般化最小二乗法Generalized Least SquaresGLS
再帰的最小二乗法Recursive Least SquaresRLS

初見だとわけわからない名前だなと思っていましたが、並べてみると素直な英語でわかりやすい英語でした。
直訳でそのまま意味がわかりそうですね。

今回は統計でよく使われる最小二乗法を使用するのでOLS(Ordinary Least Squares)で目的変数、説明変数を渡してあげます。

統計データを出すコードは下記となります。

y = data['price']
x = sm.add_constant(data['sqft_above'])
results = sm.OLS(y, x).fit()
results.summary()

すると下記の表が出力されます。

この時に下記のエラーが出た方は欠損値処理がうまくいっていないので、上記の欠損値処理の箇所のコードを再度実行してください。

exogの説明変数の中にインフィニティ(∞)かNanの欠損値があるというエラーです。
今回のデータでエラーが出た場合は欠損値のエラーがでているということになります。

回帰表の理解

今回の回帰分析で見るべき項目は下記となります。

P値は0.05より小さい場合、その値は有意であるということです。
今回のケースだと切片も傾きも有意があるということになります。

回帰式の$\hat{y_i}= b_0+ b_1 x_i$の切片$b_0$が59950で$b_1$の傾きが268.4726ということです。
なので今回の回帰式は$\hat{y_i}= 59950+ 268.4726x_i$だとわかったのでその直線式をグラフに追加していきたいと思います。

回帰直線の追加

上記でわかった傾きと切片を使ってyhatという直線$\hat{y_i}= b_0+ b_1 x_i$の関数を作成して表示させていきます。

plt.scatter(x1, y)
yhat = 268.4726*x1 + 59950
flg = plt.plot(x1, yhat, lw=4, c='orange')
plt.xlabel('SAT', fontsize = 20)
plt.ylabel('GPA', fontsize = 20)
plt.show()

実行すると下記のように表が出てきます。

これで最小二乗法での単回帰直線を引くことができました。
上記で重要だと記載した記載したP値と解説と最小二乗法とはどういう直線なのかということは別の記事で解説していきたいと思います。

コメント

タイトルとURLをコピーしました