店舗売上 - 時系列予測（機械学習を使って食料品の売上を予測する）

## はじめに
プログラミングスクールで学んだ、データ分析手法を活用して、
Kaggle内のテーマ「店舗売上 - 時系列予測（機械学習を使って食料品の売上を予測する）」の分析を実装しましたので、この記事にて、共有させていただきます。

## 解決したい社会課題
昨今、企業が保有しているデータを使って、マーケティングや売上予測を実施する事例が増え、今後もあらゆる場所でデータ分析の需要が高まると言われています。

特にこのテーマは、原油価格や地震によって経済が左右されやすい日本と似ていることもあり、また今後の参考データになると考え、選定した。

## 分析するデータ
Kaggle提供のデータ6つ

・holidays_events.csv
・oil.csv
・stores.csv
・test.csv
・train.csv
・transactions.csv

## 実行環境
パソコン：Windows
開発環境：Kaggle Notebook
言語：Python
ライブラリ：Pandas、Numpy、Matplotlib

## 分析の流れ
1.データ読み込み
2.データ確認
3.前処理
4.モデル学習
5.予測・評価
6.データ修正
7.再モデル学習
8.再予測・評価

## 分析の過程
### 1.データ読み込み
実行したコード
```python:python
df_train = pd.read_csv('/kaggle/input/store-sales-time-series-forecasting/train.csv')
df_test = pd.read_csv('/kaggle/input/store-sales-time-series-forecasting/test.csv')

df_oil = pd.read_csv('/kaggle/input/store-sales-time-series-forecasting/oil.csv')
df_sample_submission = pd.read_csv('/kaggle/input/store-sales-time-series-forecasting/sample_submission.csv')
df_holidays_events = pd.read_csv('/kaggle/input/store-sales-time-series-forecasting/holidays_events.csv')
df_stores = pd.read_csv('/kaggle/input/store-sales-time-series-forecasting/stores.csv')
df_transactions = pd.read_csv('/kaggle/input/store-sales-time-series-forecasting/transactions.csv')
````

### 2.データ確認
実行したコード
```python:python
display(df_train.head(), df_test.head(), df_oil.head(), df_sample_submission.head(), df_holidays_events.head(), df_stores.head(), df_transactions.head())

display(len(df_train), len(df_test), len(df_oil), len(df_sample_submission), len(df_holidays_events), len(df_stores), len(df_transactions))
````

![スクリーンショット 2024-03-27 190035.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/0d1ac8f0-4ee4-50d5-64c0-2e6d1586cc9d.png)
![スクリーンショット 2024-03-27 190044.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/ba778679-d451-f01e-b229-00ef1ba26e25.png)
![スクリーンショット 2024-03-27 190054.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/c5795de2-dda8-f01a-5a9a-9ba9b7631dff.png)
![スクリーンショット 2024-03-27 190104.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/0c8e6986-0329-764e-d079-52e4acb86d01.png)
![スクリーンショット 2024-03-27 190121.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/544dc330-6e9b-7cb8-114a-08246c39b470.png)
![スクリーンショット 2024-03-27 190131.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/891d17d6-99a3-348c-037b-72eb963d829f.png)
![スクリーンショット 2024-03-27 190139.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/c6c77cdc-65e1-33e2-fd26-1b83840adbb8.png)

3000888
28512
1218
28512
350
54
83488

→ 「type(休日)、locale、locale_name、description、transferred、city、state、type(店のタイプ？)、cluster」は、trainデータに対して、データ数が少ないため、 trainには結合しないこととする。

#### ・欠損値を除いた原油価格の可視化
```python:python
import matplotlib.pyplot as plt
import matplotlib.style
matplotlib.style.use('ggplot')

is_not_null = df_merged["dcoilwtico"].notnull()
df_filtered = df_merged[is_not_null]

fig, ax = plt.subplots()
ax.plot(df_filtered['date'], df_filtered['dcoilwtico'], color='blue')

fig.set_size_inches(16, 8)  # 16インチx8インチのサイズに
plt.show()
df_filtered.shape
````
![IMG_1006.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/013e5bc0-89bf-d67c-4c38-718b1234a350.png)

→ 100前後と50前後と大きく2つ分けられることが見えるため、前処理では、直前値と直後値で補完することにする。

#### ・df_transactionsの可視化と傾き確認
```python:python
# 取引数(transactions)を日付ごとに合算する
daily_transactions = df_transactions.groupby("date")["transactions"].sum()

# plot
plt.figure(figsize=(20, 6))
plt.plot(daily_transactions.index, daily_transactions.values, marker="o", linestyle="-")
plt.xlabel("date")
plt.ylabel("total_transactions")
plt.title("total_transactions vs date")
plt.grid()
plt.show()

# 傾き確認
from sklearn.linear_model import LinearRegression
x = daily_transactions.index
y = daily_transactions.values
x = np.array(daily_transactions.index).reshape(-1, 1) # n行1列
y = np.array(daily_transactions.values)
model = LinearRegression()
model.fit(x, y)
print('intercept: ', model.intercept_)
print('coefficient; ', model.coef_[0])

# 相関係数
print(df_merged[['dcoilwtico','sales']].corr())
````

![スクリーンショット 2024-03-27 193537.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/315e0471-0767-c871-0999-ee9ffeee309b.png)

intercept:  -23321.907040687525
coefficient;  7.51381796055891e-14
→日を追うごとに、微増している。

|  | dcoilwtico | sales |
|:-----------|------------:|:------------:|
| dcoilwtico       | 1.000000        | -0.074808         |
| sales     | -0.074808      | 1.000000       |

→原油価格とsalesとの相関はなし。

### 3.前処理

#### ・日付型(datetime型)に変換
```python:python
df_train['date'] = pd.to_datetime(df_train['date'])
df_oil['date'] = pd.to_datetime(df_oil['date'])
df_holidays_events['date'] = pd.to_datetime(df_holidays_events['date'])
df_transactions['date'] = pd.to_datetime(df_transactions['date'])
df_test['date'] = pd.to_datetime(df_test['date'])
````

#### ・年月日別に分割
```python:python
# df_mergedの年月日分け
df_merged['year'] = df_merged['date'].dt.year
df_merged['month'] = df_merged['date'].dt.month
df_merged['day'] = df_merged['date'].dt.day

# df_test_mergedの年月日分け
df_test_merged['year'] = df_test_merged['date'].dt.year
df_test_merged['month'] = df_test_merged['date'].dt.month
df_test_merged['day'] = df_test_merged['date'].dt.day
````

#### ・欠損値(dcoilwtico)の補完
```python:python
# 欠損値を直前の値で補完
df_merged['dcoilwtico'] = df_merged['dcoilwtico'].ffill()
df_merged['dcoilwtico'].isnull().sum() # 1782

# 欠損値を直後の値で補完
df_merged['dcoilwtico'] = df_merged['dcoilwtico'].bfill()
df_merged['dcoilwtico'].isnull().sum() # 0
````

#### ・欠損値補完後の可視化
```python:python
fig, ax = plt.subplots()
total_coilwtico.plot(x='date', y='dcoilwtico', color='blue', ax=ax)

# 背景色を変える日付を指定
target_date1 = '2016-04-16'
# 指定した日付以降の領域に背景色を塗る
ax.axvspan(target_date1, total_coilwtico.index[-1], color='gray', alpha=0.5)

# 背景色を変える日付を指定
target_date2 = '2017-08-16'
# 指定した日付以降の領域に背景色を塗る
ax.axvspan(target_date2, total_coilwtico.index[-1], color='gray', alpha=0.3)

fig.set_size_inches(16, 8)
ax.set_xlim('2012-12-20', '2017-10-01')
ax.set_ylim(0, 130)
plt.show()
````
![IMG_1008.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/bc5e0ddc-4659-6dc9-a5b0-b8afb30366d1.png)


→ 地震の影響で原油価格が大きく変わることはなかったと考えられる。

## 4.モデル学習
##### 実行したコード
```python:python
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error,r2_score

x = df_merged[['store_nbr', 'family', 'onpromotion', 'dcoilwtico', 'year',
       'month', 'day']]
y = df_merged['sales']

x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)

model = LinearRegression()
model.fit(x_train, y_train)
````
## 5.予測・評価
##### 実行したコード
```python:python
# r^2
print("r^2 (train):", r2_score(y_train, y_train_pred))
print("r^2 (test):", r2_score(y_test, y_pred))
````
r^2 (train): 0.19438794330824127
r^2 (test): 0.19377239220436604

## 6.データ修正
残差プロットを確認。

##### 実行したコード
```python:python
# 残差プロット
y_pred = model.predict(x_test)
plt.scatter(y_pred, y_pred - y_test, color='blue')
plt.hlines(y=0, xmin=0, xmax=30000, color='black')
plt.title('Residual plot')
plt.xlabel('Predicted Values')
plt.ylabel('Residuals')
plt.grid()
plt.show()
````
![__results___31_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3743369/69af0947-0600-d6fd-4ad6-14a8a9f54796.png)

→外れ値を省く。
##### 実行したコード
```python:python
residuals = y_pred - y_test

# 残差が-20000以下または20000以上のデータポイントのみ抽出
filtered_residuals = residuals[(residuals <= -20000) | (residuals >= 20000)]

# 抽出したデータポイントのインデックスを取得
filtered_indices = filtered_residuals.index
len(filtered_indices) # 24

# データフレームから指定されたデータポイントを含む行を除外する
df_merged = df_merged[~df_merged.index.isin(filtered_indices)]
````

## 7.再モデル生成
ダミー数値化とランダムフォレストで、再度モデル学習。
##### 実行したコード
```python:python
# ダミー変数化（train）
add_dum = pd.get_dummies(df_merged['family'], dtype='uint8', columns=['family'])
dummy_train = pd.concat([df_merged, add_dum], axis=1)
dummy_train.drop(columns=['family'], inplace=True)
dummy_train.tail()

# ダミー変数化（test）
add_dum = pd.get_dummies(df_test_merged['family'], dtype='uint8', columns=['family'])
dummy_test = pd.concat([df_test_merged,add_dum], axis=1)
dummy_test.drop(columns=['family'], inplace=True)
dummy_test.tail()

# モデル学習
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score

train_features = dummy_train[['store_nbr', 'onpromotion', 'dcoilwtico', 'year',
               'month', 'day', 'AUTOMOTIVE', 'BABY CARE', 'BEAUTY', 'BEVERAGES', 'BOOKS', 'BREAD/BAKERY',
               'CELEBRATION', 'CLEANING', 'DAIRY', 'DELI', 'EGGS', 'FROZEN FOODS', 'GROCERY I',
               'GROCERY II', 'HARDWARE', 'HOME AND KITCHEN I', 'HOME AND KITCHEN II',
               'HOME APPLIANCES', 'HOME CARE', 'LADIESWEAR', 'LAWN AND GARDEN', 'LINGERIE',
               'LIQUOR,WINE,BEER', 'MAGAZINES', 'MEATS', 'PERSONAL CARE', 'PET SUPPLIES', 
               'PLAYERS AND ELECTRONICS', 'POULTRY', 'PREPARED FOODS', 'PRODUCE', 'SCHOOL AND OFFICE SUPPLIES', 'SEAFOOD']].values
train_target = dummy_train['sales'].values

train_data, test_data, train_labels, test_labels = train_test_split(train_features, train_target, test_size=0.2, random_state=0)

model = RandomForestRegressor(n_jobs=-1, random_state=0)
model.fit(train_data, train_labels)
````

## 8.再予測
##### 実行したコード
```python:python
predicted_labels = model.predict(test_data)
y_train_pred = model.predict(train_data)

# r^2
print("r^2 (train):", r2_score(train_labels, y_train_pred)) 
print("r^2 (test):", r2_score(test_labels, predicted_labels))
````
##### 可視化したデータ
r^2 (train): 0.9875103196178642
r^2 (test): 0.91339862331505

## 分析したデータから得られた情報(考察)
・原油価格によって売上は大きく影響しなかった。
・一時、地震発生により、食料品、飲み物、お肉（主に食料品）の需要が高くなった。
・線形回帰モデルから、モデル変更と外れ値の除外により、決定係数がかなり改善された。

## まとめ
第一に、モデル学習の前の「データ確認の重要性」をとても感じた。
具体的な策として、
・標準化をして各特徴量の影響度を調べる
・店舗差、エリア別の傾向を調べる
が挙げられる。

またモデル選定についても、私の引き出しが少なく、ベストなモデル選定ができたかどうかの自信がないため、統計検定準1級レベルの統計学の勉学に努めようと思います。