東京都心部におけるシェアサイクルの貸出状況分析および重回帰モデルの構築（2024/11)

１．分析概要

docomoバイクシェアとOpenStreetのリアルタイムデータを使用し、東京都心部におけるシェアサイクルの利用状況を分析しました。アウトプットとしては、PythonおよびQGISを用いて、貸出数のヒートマップを作成し、潜在需要予測モデルを構築しました。

２．目的と背景

シェアサイクル市場状況

シェアサイクルは1965年にオランダで始まり、日本では2005年に世田谷区で導入されて以降、全国269都市に普及しています。平均貸出距離は2～3kmで、通勤や買い物などの短距離移動に利用され、主要都市では1000以上のポートが設置されています。2022年度には年間4000万件以上の利用があり、毎年約10％の増加が見られます。

課題としてはポート不足や交通事故リスクへの対策が挙げられ、ヘルメット着用の努力義務や罰則強化も進んでいます。都市部の混雑緩和やCO2削減（年間約1万トン）への期待があり、国交省も補助金やインフラ整備を推進中です。

個人的背景

自転車が好きで地元や旅行先でもよくシェアサイクルを利用してきたことと、Pythonを用いたデータ分析の実力を試したいという思いから今回実施しました。

分析の目的

東京都心部のシェアサイクル利用状況を分析し、潜在需要や課題を明らかにすることを目指します。

３．分析対象

本分析では、GBFS（General Bikeshare Feed Specification）フォーマットに準拠したバイクシェアのオープンデータを使用しました。対象となるデータは、国内企業が提供するdocomoバイクシェアとOpenStreet（HELLO CYCLING）のリアルタイムデータです。

データの概要

docomoバイクシェア：東京のみ。1分単位でデータが更新されます。
OpenStreet（HELLO CYCLING）：全国を対象。5分単位での更新。
今回は東京のデータに絞って分析を行いました。

■主な使用データ内容

• ポートID
• ポートの緯度経度
• ポート情報の更新タイムスタンプ
• 各ポートにおけるリアルタイムの貸出可能数と返却可能数
  ※貸出可能数と返却可能数の変化から貸出数と返却数を推定
  データ取得期間: 9月6日（金）から9月12日（木）までの1週間

■サンプル数
総サンプル数：約2000万
1時間あたり1ポートのサンプル数：docomoは60サンプル、OpenStreetは12サンプル
このデータを基に、東京都心部のシェアサイクルの利用状況を分析しました。

４．実行環境

OS: Windows 10 Home 64-bit
CPU: Intel Core i7-8700 @ 3.20GHz
RAM: 16GB
Python バージョン: 3.9.7
IDE: Anaconda Navigator 2.5.2 / Spyder IDE 5.5.1
使用ライブラリ:（2024/11/5時点）
pandas(2.2.2), numpy(1.26.4), matplotlib(3.9.0), scikit-learn(1.4.2), statsmodels(0.14.2), geopandas(1.0.1)
GIS:QGIS(3.34.12)

５．今回の結論

貸出数の時系列グラフとヒートマップから、シェアサイクルの利用は主に通勤時間帯である7時～8時台と17時～18時台にピークが見られました。田町や大手町・神田エリアなど、複数路線が交わり、駅間が数kmあるエリアで、ラストワンマイルの通勤手段として頻繁に利用されていることがわかります。田町駅や品川駅の現場観察では、利用者の8割以上がスーツ姿のビジネスマンでした。
貸出数、駅乗降客数、従業者数データから重回帰モデルを構築しモデル評価を行いました。
今後は、潜在需要の予測と比較から、需要ギャップのあるエリアの分析を進めます。

６．分析内容

6-1.データ取得

今回は自転車シェアリングのオープンデータからデータを取得します。
以下が初期情報の取得コードとステータス情報取得コードです。
初期情報は主に、ポートID、緯度経度などの情報になります。
ステータス情報は主に、ポートID、貸出可能数、返却可能数、ステータス最終更新日時（UNIXタイムスタンプ）などの情報になります。
データはdocomoバイクシェアとOpenStreet（HELLO CYCLING）のデータがありますが、以下コードは分かり易さを重視しdocomoデータに絞って記載しています。
分析としては、OpenStreetの東京エリアデータも使用しています。
なお、データが更新される頻度は、docomoが1分ごと、OpenStreetは5分ごとです。

初期情報取得コード

#%%初期情報取得
#ライブラリインポート
import pandas as pd
import requests

#データ取得関数
def read_gbfs(bike_url):
    #リアルタイムデータを取得するためにHTTPリクエスト
    url = requests.get(bike_url)
    #リクエスト内容をテキスト変換
    text = url.text
    #JSONのリクエスト内容を読み込み
    df_json = pd.read_json(text)
    #JSONリスト構造から、DataFrame化（正規化と平坦化）
    df_data = pd.json_normalize(df_json["data"][0])
    return df_data

#docomo　station_infomation基本情報の取得    
df_docomo_stinfo = read_gbfs("https://api-public.odpt.org/api/v4/gbfs/docomo-cycle-tokyo/station_information.json")
df_docomo_stinfo.head()
df_docomo_stinfo.columns
#(1460, 6)
#Index(['lat', 'lon', 'name', 'capacity', 'region_id', 'station_id'], dtype='object')

#docomo　station_status　レンタル状況における情報の取得    
df_docomo_ststat = read_gbfs("https://api-public.odpt.org/api/v4/gbfs/docomo-cycle-tokyo/station_status.json")
df_docomo_ststat.shape
df_docomo_ststat.columns
#(1480, 7)
#Index(['is_renting', 'station_id', 'is_installed', 'is_returning','last_reported', 'num_bikes_available', 'num_docks_available'],dtype='object')

ステータス情報取得コード（1分ごと）

#%%ステータス情報の定期取得
#ライブラリインポート
import requests
import pandas as pd
import schedule
import time
from datetime import datetime

#スケジュールをクリア
schedule.clear()
#※実行の度にクリアしないと過去実行時のジョブが残り重複するので注意。

#GBFSのstation_statusデータURLを設定
status_url = "https://api-public.odpt.org/api/v4/gbfs/docomo-cycle-tokyo/station_status.json"

#必要なカラムの指定
columns = [
    "station_id", 
    "is_renting", 
    "is_installed", 
    "is_returning", 
    "last_reported", 
    "num_bikes_available", 
    "num_docks_available"
]

#取得する回数の変数と最大値を設定
count = 0
max_count = 30240

#データ取得とCSV保存関数
def get_status():
    #設定したcountを使用
    global count
    #最大値になるまで取得
    if count < max_count:
        #HTTPリクエストのレスポンス取得
        response = requests.get(status_url)
        #レスポンスが正常（200）の場合
        if response.status_code == 200:
            #JSONリクエスト内容を読み込み
            data = response.json()
            #DataFrame化と必要なカラムの指定
            df = pd.json_normalize(data["data"]["stations"])
            df = df[columns]
            #現在時刻の追加
            df["timestamp"] = datetime.now().strftime("%Y-%m-%d %H:%M:%S")
            #CSVに追記
            df.to_csv("data/docomo_station_status_data_1109.csv", mode="a", header=False, index=False)
            #何時に何回目か実行中に表示
            print(f"Data saved at {datetime.now()}(count {count + 1})")
            #取得完了する度に+1カウント
            count += 1
        else:
            #データ取得失敗を通知
            print("Failed to get data")
    else:
        #取得終了通知
        print("Max count. Stopping.")
        #スケジュールを停止し、取得を終了
        return schedule.CancelJob
        
#1分毎にデータを取得するスケジュール設定
schedule.every(1).minutes.do(get_status)

#スケジュールジョブ（スクリプト）の実行
if __name__ == "__main__":
    #CSVファイルにヘッダーを追加（初回のみ）。その他の取得方法はget_status関数と同様。
    response = requests.get(status_url)
    data = response.json()
    df = pd.json_normalize(data["data"]["stations"])
    df = df[columns]
    df["timestamp"] =datetime.now().strftime("%Y-%m-%d %H:%M:%S")
    #header=Trueでヘッダ追加
    df.to_csv("data/docomo_station_status_data_1109.csv", mode="w", header=True, index=False)
    #countが最大値になるまでループでスケジュール実行
    while count < max_count:
        #スケジュールジョブの実行
        schedule.run_pending()
        #CPU負担軽減のために1秒待機
        time.sleep(1)

6-2. 取得データの前処理とグラフ化

取得したデータ（docomo）から貸出数と返却数の推移を見るために行う前処理は以下になります。

• 廃止ポートの削除と休止ポートの欠損値補完
• 取得タイミングの誤差による重複行の削除
• 欠損時間の確認と欠損値補完
• 貸出数・返却数の計算とグラフ化

廃止ポートの削除と休止ポートの欠損値補完

取得したCSVを読み込み、計算のために日時変換やデータ型の変換をした後、
データ取得期間中に停止または廃止したポートについて、サンプル数から有効なポートのみ使用するため削除と欠損値補完をします。
廃止および停止ポートは欠損値の数や欠損時間からポートIDに着目し、docomoバイクシェアなどの公式情報から手で確認しました。
以下が削除および欠損値補完のコードです。

廃止ポートの削除と休止ポートの欠損値補完コード

#%%廃止ポートの削除と休止ポートの欠損値補完
#ライブラリインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
from datetime import datetime, timedelta


#docomoデータのCSV読み込み
docomo_stat = pd.read_csv(r"C:\\Users\\user\\Desktop\\sharecycle\\data\\docomo_station_status.csv")


##タイムスタンプ変換
#Unixタイムスタンプ（秒単位）を日時形式（UTC）に変換
docomo_stat["last_reported_date"] = pd.to_datetime(docomo_stat["last_reported"], unit="s", utc=True)

#日時形式UTC→JSTに変換
docomo_stat["last_reported_date"] = docomo_stat["last_reported_date"].dt.tz_convert("Asia/Tokyo")

#不要なタイムゾーン情報を削除
docomo_stat["last_reported_date"] = docomo_stat["last_reported_date"].dt.tz_localize(None)

#データ型object→datetime64[ns]変換
docomo_stat["timestamp"] = pd.to_datetime(docomo_stat["timestamp"], format='%Y-%m-%d %H:%M:%S')

#グループ化した計算をしやすくするため、秒を0に統一
docomo_stat["last_reported_date"] = pd.to_datetime(docomo_stat["last_reported_date"]).dt.floor("min")


##休止・廃止ポートのデータ補完・削除
#廃止ポート（3344）の削除
docomo_stat = docomo_stat[docomo_stat["station_id"] != 3344]

#停止ポートの欠損値補完
#補完するIDの開始時刻と終了時刻を定義
stop_port_info = [
    {"station_id": 2431, "start_time": "2024-09-10 19:00:00"},
    {"station_id": 10245, "start_time": "2024-09-12 11:00:00"}]
end_time = "2024-09-12 23:59:00"

#停止ポートの欠損値補完関数
def stop_port_fill(station_id, start_time, end_time):
    #開始時刻
    start_dt = datetime.strptime(start_time, "%Y-%m-%d %H:%M:%S")
    #終了時刻
    end_dt = datetime.strptime(end_time, "%Y-%m-%d %H:%M:%S")
    #作成したデータ用リスト
    filled_row = []
    #時刻に開始時刻を設定
    current_time = start_dt
    #時刻が終了時刻になるまでループ
    while current_time <= end_dt:
        #時刻をUNIXタイムスタンプ（秒）に変換
        last_reported = int(current_time.timestamp())
        #補完データ内容
        row = {
            #ポートID
            "station_id": station_id,
            #貸出可否（停止中なのでFalse固定）
            "is_renting": False,
            #ポート設置有無（物理的に存在するためTrue固定）
            "is_installed": True,
            #返却可否（停止中なのでFalse固定）
            "is_returning": False,
            #最後にステータスが更新された時刻（UNIXタイムスタンプ）
            "last_reported": last_reported,
            #貸出可能数（0固定）
            "num_bikes_available": 0,
            #返却可能数（0固定）
            "num_docks_available": 0,
            #データ取得時のタイムスタンプ（年月日時分秒）
            "timestamp": current_time.strftime("%Y-%m-%d %H:%M:%S"),
            #最後にステータスが更新された時刻（年月日時分秒）
            "last_reported_date": current_time.strftime("%Y-%m-%d %H:%M:%S")
        }
        #補完データ行を追加
        filled_row.append(row)
        #時刻を1分進める
        current_time += timedelta(minutes=1)
    #DataFrame化して返す
    return pd.DataFrame(filled_row)

#停止ポートの欠損値補完
#補完するデータがある分ループ
for info in stop_port_info:
    #停止ポート欠損値補完関数の実行
    new_data = stop_port_fill(info["station_id"], info["start_time"], end_time)
    #既存DFと結合
    docomo_stat = pd.concat([docomo_stat, new_data], ignore_index=True)

#データ型をobject→datetime64[ns]変換
docomo_stat["timestamp"] = pd.to_datetime(docomo_stat["timestamp"], format='%Y-%m-%d %H:%M:%S')
docomo_stat["last_reported_date"] = pd.to_datetime(docomo_stat["last_reported_date"], format='%Y-%m-%d %H:%M:%S')

#DFのソート
docomo_stat = docomo_stat.sort_values(by=["timestamp", "station_id"], ascending=[True, True]).reset_index(drop=True)

これで廃止ポートの削除と停止ポートの欠損値補完が完了します。

取得タイミングの誤差による重複行の削除

docomoデータの更新タイミングおよび取得時CPU負荷軽減のための１秒待機のために、取得タイミングに秒単位の誤差が生まれます。
そのため僅かに欠損する時間と、重複したデータ取得時間が発生します。
先に重複行について１行のみにします。
まずは存在する超副業の一覧を確認した後、重複行の削除、最後に１行のみになっているか確認します。

重複行の削除コード

#%%取得タイミングの誤差による重複行の削除
#重複行の確認（複数行存在を確認）
docomo_duptest_pre = docomo_stat[docomo_stat.duplicated(subset=['station_id', 'is_renting', 'is_installed', 'is_returning', 'last_reported', 'num_bikes_available', 'num_docks_available'], keep=False)]

#重複行の削除（最初の行を残し他の重複行を削除）
docomo_stat_dup = docomo_stat.drop_duplicates(subset=['station_id', 'is_renting', 'is_installed', 'is_returning', 'last_reported', 'num_bikes_available', 'num_docks_available'], keep="first")

#重複行の再確認（重複行が1行のみになっていることを確認）
docomo_duptest_post = docomo_stat_dup[docomo_stat_dup.duplicated(subset=['station_id', 'is_renting', 'is_installed', 'is_returning', 'last_reported', 'num_bikes_available', 'num_docks_available'], keep=False)]
docomo_stat_dup[(docomo_stat_dup["station_id"] == 10001) & (docomo_stat_dup["last_reported"] == 1725551583)].head()

これで重複した行の削除が完了します。

欠損時間の確認と欠損値補完

続いて、欠損時間の確認および欠損値補完をします。
以下、欠損時間の確認および欠損値補完のコードです。

欠損時間の確認と欠損値補完コード

#%%欠損時間の確認
#欠損時間確認関数
def miss_times_get(df):
    #結果リストの設定
    result = []
    #ポートID分ループ
    for station_id in df["station_id"].unique():
        #特定のポートIDのデータを取得
        df_st = df[df["station_id"] == station_id].copy()
        #ステータス最終時刻の分単位列を作成
        df_st["last_reported_min"] = df_st["last_reported_date"].dt.floor("min")
        #上限/下限となるlast_reported_dateを取得
        min_time = df_st["last_reported_min"].min()
        max_time = df_st["last_reported_min"].max()
        #１分ごとの時間を生成
        all_times = pd.date_range(min_time, max_time, freq="min")        
        #last_reported_dateをsetに変換
        exist_times = set(df_st["last_reported_min"])        
        #存在しない時間を抽出
        miss_times = [time for time in all_times if time not in exist_times]
        #欠損時間の数
        count = len(miss_times)
        #特定のポートIDにおける欠損時間を結果に追加
        output = f"station_id is {station_id}. miss_time_count{count}"
        result.append(output)
    print(f"【END】 station_id:{station_id} is miss_times{miss_times}")
    return miss_times, result

#欠損時間の確認
miss_times_pre, result_pre = miss_times_get(docomo_stat_dup)


#%%欠損時間における欠損値補完
#欠損値補完関数
def fill_miss_time(df, miss_times):
    #追加する行のリスト変数
    filled_rows = []
    #各欠損時間ごとにループ
    for time in miss_times:
        #ポートIDごとにループ
        for station_id in df["station_id"].unique():
            #該当のポートIDにおけるデータを取得
            station_data = df[df["station_id"] == station_id]
            #ステータス最終取得時間ごとにソート
            station_data_sort = station_data.sort_values(by="last_reported_date")
            #端のデータか確認用に前後のデータを取得
            prev_data = station_data_sort[station_data_sort["last_reported_date"] < time]
            next_data = station_data_sort[station_data_sort["last_reported_date"] > time]
            if not prev_data.empty and not next_data.empty:
                #欠損を埋める時間の前後データを取得
                prev_data = prev_data.iloc[-1]
                next_data = next_data.iloc[0]
                #前後のデータの平均値を計算（貸出可能数、返却可能数）
                avg_bikes = (prev_data["num_bikes_available"] + next_data["num_bikes_available"]) / 2
                avg_docks = (prev_data["num_docks_available"] + next_data["num_docks_available"]) / 2
                #UNIXタイムスタンプ変換
                last_reported = int(time.timestamp())
                #欠損行の作成
                new_row = {
                    "station_id": station_id,
                    "is_renting": True,
                    "is_installed": True,
                    "is_returning": True,
                    "last_reported": last_reported,
                    "num_bikes_available": avg_bikes,
                    "num_docks_available": avg_docks,
                    "timestamp": time,
                    "last_reported_date": time
                    }
                #欠損行をリストに追加
                filled_rows.append(new_row)
            else:
                print(f"prev/next_data is empty.time is {time},station_id is {station_id}")
    #欠損補完リストをDFに変換
    filled_df = pd.DataFrame(filled_rows)
    #欠損補完DFを元のDFに結合しそーと
    df = pd.concat([df, filled_df]).sort_values(by=["station_id", "last_reported_date"]).reset_index(drop=True)
    return df

#欠損行の補完
docomo_stat_fill = fill_miss_time(docomo_stat_dup, miss_times_pre)

#欠損値補完後の欠損時間チェック
miss_times_post, result_post = miss_times_get(docomo_stat_fill)

#ステータス最終取得時間でソート
docomo_stat_fill = docomo_stat_fill.sort_values(by=["station_id", "last_reported_date"]).reset_index(drop=True)

#CSV出力
docomo_stat_fill.to_csv(r"C:\\Users\\user\\Desktop\\sharecycle\\data\\docomo_station_status_filled.csv", index=False)

これで取得タイミングの誤差による欠損値の補完が完了します。 ここまででグラフ化するために取得したデータの前処理は終了です。

貸出数・返却数の計算とグラフ化

貸出可能数から貸出数と返却数を推定し、時系列でグラフ化して利用傾向を見ます。
推定は各１分ごとの貸出可能数について前後の時間の差を取り、１減っていたら１台貸出された、１増えていたら1台返却されたと考えます。
現実的には駐輪台数のバランスをとるために、トラックなどで一斉に運ばれる運用がありますが今回のデータからの特定は困難のため考慮しないこととします。
以下、貸出数と返却数の計算とグラフ化のコードです。

貸出数と返却数の計算とグラフ化のコード

#%%貸出数と返却数のグラフ化
##貸出数と返却数計算
#差分を計算
docomo_stat_fill["bike_diff"] = docomo_stat_fill.groupby("station_id")["num_bikes_available"].diff()

#貸出数と返却数の計算と列作成
docomo_stat_fill["rental_count"] = docomo_stat_fill["bike_diff"].apply(lambda x: -x if x < 0 else 0)
docomo_stat_fill["return_count"] = docomo_stat_fill["bike_diff"].apply(lambda x: x if x > 0 else 0)

#last_reported_date（分）ごとに、貸出数と返却数の合計値を計算
docomo_rental_bikes_min = docomo_stat_fill.groupby("last_reported_date")["rental_count"].sum()
docomo_rental_bikes_min.head()
docomo_rental_bikes_min.index

docomo_return_bikes_min = docomo_stat_fill.groupby("last_reported_date")["return_count"].sum()
docomo_return_bikes_min[docomo_return_bikes_min < 0]
docomo_return_bikes_min.head()


##貸出数と返却数のグラフプロット
#docomo　貸出数
plt.figure(figsize=(10, 5), dpi=500)
plt.plot(docomo_rental_bikes_min.index, docomo_rental_bikes_min.values, marker="o", color="b")
plt.title("docomo東京エリア　時系列別（分）合計貸出数")
plt.xlabel("日時（分）")
plt.ylabel("合計貸出数（台）")
plt.tight_layout()
plt.show()

#docomo　返却数
plt.figure(figsize=(10, 5), dpi=500)
plt.plot(docomo_return_bikes_min.index, docomo_return_bikes_min.values, marker="o", color="r")
plt.title("docomo東京エリア　時系列別（分）合計返却数")
plt.xlabel("日時（分）")
plt.ylabel("合計返却数（台）")
plt.tight_layout()
plt.show()

結果のグラフがこちらになります。

■図1:docomo東京エリア　時系列別（分）合計貸出数　2024/9/6（金）～2024/9/12（木）

■図2:docomo東京エリア　時系列別（分）合計返却数　2024/9/6（金）～2024/9/12（木）

x軸とy軸の範囲は適宜設定してください。
平日5日間に特に1日のなかで、大きなピークが2つ見られます。
土日2日間は日中に貸出返却の傾向が固まって見えます。
貸出数について、2024/9/6（金）に絞ってみます。

■図3:docomo東京エリア　時系列別（分）合計貸出数　2024/9/6（金）

概ね7時～8時と18時～19時ごろに貸出数のピークが見られます。
出勤時と退勤時の傾向が典型的に見て取れます。
概ねのピークの時間帯および平日と休日の傾向が見られたので、続いてGIS（QGIS）を用いて東京エリアにおける利用率が高いエリアを探ります。

6-3.東京都心エリアのヒートマップ

Pythonを使い貸出数と返却数のCSVデータと東京都のマップデータ（小ゾーンのshpデータ）を結合し、QGISを使いヒートマップを作成する。
東京都のマップデータは、東京都市圏交通計画協議会より無償で使用できる。
ゾーンとは東京都23区より細かく都市計画などに沿って区分された考え方で、東京都市圏交通計画協議会より、以下の区分で考えられています。

■図4:ゾーン（出展：東京都市圏交通計画協議会）

今回は都心部におけるサンプル数が十分とれるかという考えのもと、一番細かい小ゾーンでも都心部であれば数百以上とれることが見込まれたため小ゾーンを採用しました。

以下がCSVデータとマップデータを結合し、小ゾーン単位の貸出数と返却数データを作成するコードになります。
（個人都合により再度、貸出可能数から貸出数と返却数計算をしています。）

CSVデータとマップデータの結合コード

#%%小ゾーン単位の貸出数・返却数データ作成（docomo）

#ライブラリインポート
import pandas as pd
import geopandas as gpd
import requests
from shapely.geometry import Point

#CSV読み込み
dcm_status = pd.read_csv(r"C:\\Users\\user\\Desktop\\sharecycle\\data\\docomo_station_status_filled.csv")

#station_idをint→str変換
dcm_status["station_id"] = dcm_status["station_id"].astype(str)

#必要なカラムのみ取得
dcm_status = dcm_status[["station_id", "last_reported_date", "num_bikes_available", "num_docks_available"]]

#last_reported_dateの名前修正
dcm_status = dcm_status.rename(columns={"last_reported_date": "date"})

#dateをobject→datetime型に変換
dcm_status["date"] = pd.to_datetime(dcm_status["date"])

#データ取得関数
def read_gbfs(bike_url):
    #APIエンドポイントからデータ取得
    url = requests.get(bike_url)
    #中身（body）を取得（JSON形式
    text = url.text
    #JSONデータdf化（まだ先頭に入れ子式になってる）
    df_json = pd.read_json(text)
    #normalizeでフラットなdf化に変換
    df_data = pd.json_normalize(df_json["data"][0])
    return df_data

#station_infomation の取得
dcm_info = read_gbfs("https://api-public.odpt.org/api/v4/gbfs/docomo-cycle-tokyo/station_information.json")

#必要情報に絞る
dcm_info = dcm_info[["lat", "lon", "name", "station_id"]]

#statusに結合目的で、station_idの先頭0を削除のために一度int変換して戻す。
dcm_info["station_id"] = dcm_info["station_id"].astype(int)
dcm_info["station_id"] = dcm_info["station_id"].astype(str)

#infoとstatusの結合
dcm_gbfs = pd.merge(dcm_info, dcm_status, on="station_id", how="inner")

##貸出数と返却数
#差分を計算
dcm_gbfs["bike_diff"] = dcm_gbfs.groupby("station_id")["num_bikes_available"].diff()

#貸出数（貸出可能数が減っていたら、つまり差分がマイナスだったら、その絶対値が貸し出された数）
dcm_gbfs["rental_count"] = dcm_gbfs["bike_diff"].apply(lambda x: -x if x < 0 else 0)
#返却数（貸出可能数が増えていたら、つまり差分がプラスだったら、その値が返却された数）
dcm_gbfs["return_count"] = dcm_gbfs["bike_diff"].apply(lambda x: x if x > 0 else 0)

#小ゾーンshpファイル読み込み
szone_gdf = gpd.read_file(r"02.GIS/01.shp/H30_gis/H30_szone.shp")

#座標系をEPSG4326に変換
szone_gdf = szone_gdf.to_crs(epsg=4326)

##緯度経度と小ゾーンの結合
#緯度経度をPointオブジェクトに変換し、geometry列を作成
dcm_gbfs["geometry"] = dcm_gbfs.apply(lambda x: Point((x["lon"], x["lat"])), axis=1)

#GeoDataFrameに変換(座標系CRS情報付与、どこにあるか表す幾何情報、地理上の計算の空間操作ができる)
dcm_gdf = gpd.GeoDataFrame(dcm_gbfs, geometry="geometry", crs="EPSG:4326")

#空間結合（緯度経度が小ゾーンのポリゴンに属しているか判定）
dcm_gdf_sjoin = gpd.sjoin(dcm_gdf, szone_gdf, how="left", predicate="within")

#hour単位を作成
dcm_gdf_sjoin["date_hour"] = dcm_gdf_sjoin["date"].dt.floor("H")

#小ゾーンごとに1時間ごとの貸出数と返却数の合計を計算
dcm_szone_hour = dcm_gdf_sjoin.groupby(["szone", "date_hour"]).agg({
    "rental_count": "sum",
    "return_count": "sum"}).reset_index()

#ピボットテーブルを作成して、各ゾーンごとの時間帯別貸出数、返却数を列に展開
dcm_szone_hour_rental = dcm_szone_hour.pivot(index="szone", columns="date_hour", values="rental_count").fillna(0)
dcm_szone_hour_return = dcm_szone_hour.pivot(index="szone", columns="date_hour", values="return_count").fillna(0)

#結果をCSV保存
dcm_szone_hour_rental.to_csv(r"data/zone_csv/docomo_szone_rental_count.csv")
dcm_szone_hour_return.to_csv(r"data/zone_csv/docomo_szone_return_count.csv")

これでCSVデータとマップデータを結合したCSVファイルが完成です。 もう1ステップとして、各小ゾーンごとをマップ上で比較する上で基準を統一するために、各値を小ゾーンの面積で割り正規化します。 正規化するコードは以下になります。

面積による正規化コード

#%%小ゾーン単位の貸出数と返却数（正規化）

#szone shpファイル読み込み
s_gdf = gpd.read_file(r"02.GIS/01.shp/H30_gis/H30_szone.shp")
s_gdf.head()

#面積の計算
s_gdf["面積（km²）"] = s_gdf.geometry.area
s_gdf["面積（km²）"] = s_gdf["面積（km²）"] / 1000000

##貸出数の正規化（面積で割る）
#計算用にDFをコピー
szone_hour_rental_area = dcm_szone_hour_rental.copy()

#面積列を貸出数DFに結合
szone_hour_rental_area["面積（km²）"] = szone_hour_rental_area.index.map(s_gdf.set_index("szone")["面積（km²）"])

#時間帯列の取得
times_columns = szone_hour_rental_area.columns.drop("面積（km²）")

#各時間帯の貸出数を面積で割る
szone_hour_rental_area[times_columns] = szone_hour_rental_area[times_columns].div(szone_hour_rental_area["面積（km²）"], axis=0)

#四捨五入
szone_hour_rental_area[times_columns] = szone_hour_rental_area[times_columns].round()

#面積列の削除
szone_hour_rental_area = szone_hour_rental_area.drop(columns=["面積（km²）"])

#貸出数のCSV出力
szone_hour_rental_area.to_csv(r"02.GIS/02.csv/szone_hour_rental_area.csv")    


##返却数を面積で割る【正規化】
#計算用にDFをコピー
szone_hour_return_area = dcm_szone_hour_return.copy()

#面積列を返却数DFに結合
szone_hour_return_area["面積（km²）"] = szone_hour_return_area.index.map(s_gdf.set_index("szone")["面積（km²）"])

#時間帯列のカラム一覧取得
times_columns = szone_hour_return_area.columns.drop("面積（km²）")

#各時間帯の返却数を面積で割る
szone_hour_return_area[times_columns] = szone_hour_return_area[times_columns].div(szone_hour_return_area["面積（km²）"], axis=0)

#四捨五入
szone_hour_return_area[times_columns] = szone_hour_return_area[times_columns].round()

#面積列の削除
szone_hour_return_area = szone_hour_return_area.drop(columns=["面積（km²）"])

#CSV出力
szone_hour_return_area.to_csv(r"02.GIS/02.csv/szone_hour_return_area.csv")

これでマッピングするためのCSVデータが完成です。 これをGISを用いてマッピングします。 GISとは地理情報システムと言われ、今回の様な地図上でデータを分析等を行う際に使うシステムです。 その中でも今回は無償で使用できるQGISを使います。 QGISで使用するデータは、以下になります

• 結合した貸出数・返却数CSVデータ
• CSVファイルのデータ型を示すCSVTデータ
• マップデータ（小ゾーンshpデータとOpenStreet日本地図データ）

CSVTファイルは、以下の様にテキストベースでCSVデータの各列ごとにデータ型を指定したファイルになるので、テキストエディタを用いて作成し、CSVT形式で保存し作成します。

CSVTファイル内容

"string","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer","integer"

ここからはQGISを使いマッピングしていきます。
ヒートマップ作成までの流れとしては以下になります。

• マップデータの取り込み
• CSVデータの取り込みとテーブル結合
• 各日時におけるヒートマップの作成
• 画像編集と保存

今回QGISの詳細な使い方は割愛しますが、国土交通省にてQGISマニュアルも公開されているのでご参照ください。
今回はマニュアルに書かれていない使い方も込みでヒートマップまで作成しています。
簡単に各流れの操作について説明します。

マップデータの取り込み

初めに国土交通省のマニュアルの通り、ベースとなる日本地図をダウンロードとQGISに取り込みマッピングします。
今回はOpenStreetの地図を使用しましたが、標準地図など好みの地図を使用して問題ありません。
その後、小ゾーンのshpデータをQGIS上にドラッグアンドドロップし、マッピングします。
ここまでで日本地図の上に、関東圏が小ゾーン単位で区分されたマップが表示されます。
■図5:小ゾーンマッピング例

CSVデータの取り込みとテーブル結合

Pythonで作成した貸出数CSVデータをQGIS上にドラッグアンドドロップすることでQGISに取り込まれます。
その後、QGIS上で小ゾーンデータのテーブル結合機能によりCSVデータを結合させます。
これにより、小ゾーン単位の各時間帯における貸出数/返却数のヒートマップを表示できる様になります。
尚今回は、東京都23区データ（太枠）、路線データ取り込んでいます。

各日時におけるヒートマップの作成

作成したヒートマップにおける時間帯を指定し、値の範囲をピーク時に適宜合わせ、1週間分全ての時間帯のヒートマップを作成しました。

画像編集と保存

作成したヒートマップごとに画像を作成し時系列で比較するため、QGIS上で値の凡例や方角等含め編集し、全ての時間帯におけるヒートマップの画像を作成し保存しました。

結果のヒートマップが以下になります。
PCから閲覧頂いている場合、各画像をクリックしてブラウザ上でタブを行き来すると変化が見えます。

■図7:2024/9/10（月）8時台　貸出数ヒートマップ

■図7:2024/9/10(月)18時台　貸出数ヒートマップ

■考察
各平日でエリアの傾向を見た結果、貸出数の推移グラフと同様に同じ傾向が確認されました。
特に赤く貸出しされているエリアに着目して見ると、都心に近い且つ、複数の路線が並んでいるまたは交差しているエリアが多いことが見て取れます。
路線・駅間に徒歩では距離がある場所への移動に利用されている可能性が考えられます。
貸出数の推移グラフからも出退勤ピークに多く利用されていたことから、朝は駅からオフィスまたは家から駅、夕方はその逆のパターンとしてラストワンマイルに活用されていることが考えられます。
今回更に実際の現場がどうなっているか、9/27（金）8時台後半から9時台前半に田町駅と品川駅で現場観察をしました。

6-3.田町駅・品川駅ポート現場観察

■図8:田町駅前ポートのdocomoバイクシェアアプリ（左）、8時50分頃現場写真（右）

田町駅および品川駅を合計約40分ほど観察したところ、目視だけで貸出・返却含め46人が利用し、貸出が27人、返却が19人でした。
都心は特に台数が多く、駅近辺のポートにおいては1台も存在しない状況は今回の時間においては見られませんでした。
利用者の想定年齢および服装から8割以上がスーツを着た人でビジネスマンであることが確認できました。
残り2割も私服の2,30代の人たちのため、ビジネスマンの可能性があり得えます。
一部利用があるのではと想定していた子持ちのママさんについては、自前の子ども乗せ電動自転車で移動している姿が複数台見られ、シェアサイクルの利用は見られませんでした。

現場観察の結論として、ビジネスマンのラストワンマイル利用が習慣化されていることが明確になりました。

次節では貸出数データを使って潜在的に需要があるエリアを探るべくモデル構築を行います。

6-4.潜在需要予測モデルの構築

どの小ゾーンエリアにおける貸出数に潜在的な需要があるか分析します。
目的変数は 貸出数（各小ゾーンにおける単位面積当たり） とします。
説明変数は各小ゾーンにおける、

• 駅乗降客数
• 従業者数
• 曜日（平日）
• 時間帯

の4つとします。
説明変数は今後精度向上のために追加検討しますが、まずはシンプルにモデル構築をします。
まずは各説明変数におけるデータを作成します。

【前提】対象の小ゾーン

今回対象とする小ゾーンは9/10（月）18時においてヒートマップ化した際に、都心から離れた白いエリア、つまり単位面積当たりの貸出数が10未満のエリアについては対象外とし、都心部に絞ったエリアでモデル構築します。
つまりは、図7における赤いエリアが対象となります。

【説明変数】従業者数

働いている人を説明変数に入れるべく、各小ゾーンにおける従業者数を計算します。
統計局のe-Statより、経済センサスの2014年データである、東京都における産業（大分類）別・従業者規模別全事業所数及び男女別従業者数データを今回は使用しました。
本データは小ゾーンではなく、小地域と呼ばれる小ゾーンよりも更に細かいエリア単位の従業者数データになります。
よって、一度QGIS上にデータをマッピングし、小ゾーン単位の従業者数を同じくQGIS上で計算しCSV化しました。
小地域のマップデータ（shp形式）は従業者数データと同じく統計局e-Statの東京都の小地域における境界データを使用しました。
方法の詳細は割愛しますが、簡単な流れとしては、QGIS上で小地域の重心をとり、重心が存在する小ゾーンを紐づけて、小ゾーン単位で従業者数の合計を取りCSV化しました。

【説明変数】駅乗降客数

各小ゾーンにおける駅乗降客数のデータを使うために、東京における駅別乗降客数データを使います。
国土数値情報の駅別乗降客数データより1日の駅別乗降客数を今回は使用しました。
このデータをPythonにて小ゾーンと空間結合させ、小ゾーンごとの駅乗降客数を計算しました。
詳しいプログラムは後述の前処理に含めます。
尚、今回対象とした小ゾーンにおける駅乗降客数データの有無を確認したところ、欠損値があったため、欠損値が存在する小ゾーンについてQGISで駅の有無を目視で確認し、駅が存在する小ゾーンについては手動で乗降客数の値を補完しました。

【モデル構築】前処理とデータ結合

各説明変数の前処理とデータ結合をします。
主な処理としては、

• 駅別乗降客数から小ゾーンごとの駅乗降客数の計算
• 駅乗降客数の欠損値確認と補完
• 駅乗降客数の1時間単位振り分け（貸出数の分布を使用）
• 従業者数の欠損値確認（欠損値無し）
• データの結合
• ダミー変数化とモデル構築・評価

になります。
以下がコードです。

前処理とデータ結合

#%%潜在需要予測モデル 前処理とデータ結合
#ライブラリインポート
import pandas as pd
import geopandas as gpd
import statsmodels.api as sm
from sklearn.model_selection import train_test_split

##【駅別乗降客数から小ゾーンごとの駅乗降客数の計算】
#駅別乗降客数CSV読み込み
stations = gpd.read_file(r"C:\\Users\\user\\Desktop\\sharecycle\\data\\令和3年全国駅別乗降客数\\S12-22_NumberOfPassengers.shp")

#小ゾーン読み込み
szone_gdf = gpd.read_file(r"02.GIS/01.shp/H30_gis/H30_szone.shp")

#CRS一致を確認し、異なる場合は変換
stations = stations.to_crs(szone_gdf.crs)

#空間結合
stations_szone = gpd.sjoin(stations, szone_gdf, how="left", predicate="within")

#必要カラムの選択
stations_szone = stations_szone[["S12_001", "S12_001c", "S12_003", "S12_047", "S12_049", "szone", "CityName"]]

#szoneごとに乗降客数を足し合わせ
stations_szone_total = stations_szone.groupby("szone")["S12_049"].sum().reset_index()

#型変換　szone[float→int]
stations_szone_total["szone"] = stations_szone_total["szone"].astype("int64")

#乗降客数があるszoneを調べるためにCSV化
stations_szone_total = stations_szone_total.reset_index(drop=True)
stations_szone_total.to_csv(r"02.GIS/02.csv/passengers_szone_check.csv", index=False)

#単位面積当たりの貸出数データ　CSV読み込み
szone_hour_rental = pd.read_csv(r"02.GIS/02.csv/szone_hour_rental.csv")

#対象の小ゾーンを指定（貸出数が10未満の行（szone）を排除）
szone_hour_rental_fil = szone_hour_rental[szone_hour_rental["2024-09-10 18:00:00"] >= 10]

#不要なszoneの削除と不足szoneの追加
stations_passengers = stations_szone_total.merge(szone_hour_rental_fil[["szone"]], on="szone", how="right")



##【駅乗降客数の欠損値確認と補完】
#名称変更
stations_passengers = stations_passengers.rename(columns={"szone": "小ゾーン", "S12_049": "駅乗降客数"})

#欠損値確認
stations_passengers["駅乗降客数"].isnull().sum()

#欠損しているszoneの確認（追加szone）
stations_passengers[stations_passengers["駅乗降客数"].isna()]["小ゾーン"]

#欠損値補完
#各エリアについてQGISのマップ上で視覚的に駅の有無を確認し、存在するエリアは手動で値を代入。
#乗客数の全量はエリア内に駅が存在する状態、半分はエリアにまたがっている状態を指す。
#小ゾーン3420、3421の南砂町の半分を代入
stations_passengers.loc[stations_passengers["小ゾーン"].isin([3420, 3421]), "駅乗降客数"] = 49802 / 2

#小ゾーン7202に三鷹の全量を代入
stations_passengers.loc[stations_passengers["小ゾーン"] == 7202, "駅乗降客数"] = 147296

#小ゾーン8221に稲城の半分を代入
stations_passengers.loc[stations_passengers["小ゾーン"] == 8221, "駅乗降客数"] = 16630 / 2

#小ゾーン8222に稲城半分を追加
stations_passengers.loc[stations_passengers["小ゾーン"] == 8222, "駅乗降客数"] += 16630 / 2

#その他駅の無い小ゾーンの駅乗降客数に0を代入
stations_passengers.loc[stations_passengers["駅乗降客数"].isna(), "駅乗降客数"] = 0



##【駅乗降客数の1時間単位振り分け（貸出数の分布を使用）】
##平日貸出データの取得
#日付列の作成
weekday_columns = pd.date_range(start="2024-09-06", end="2024-09-13", freq="H")
#平日のみにフィルタ
weekday_columns = weekday_columns[weekday_columns.weekday < 5]
#日付列を文字列に変換
weekday_columns = weekday_columns.strftime("%Y-%m-%d %H:%M:%S").tolist()
weekday_columns.remove("2024-09-13 00:00:00")
#平日データのみ抽出
szone_week_rental = szone_hour_rental_fil[["szone"] + weekday_columns]

##貸出数データから日毎の分布を計算
#年月日時の値から日付と時間を抽出
szone_week_rental_melted = szone_week_rental.melt(id_vars="szone", var_name="日時", value_name="貸出数")
szone_week_rental_melted["日時"] = pd.to_datetime(szone_week_rental_melted["日時"])
szone_week_rental_melted["日付"] = szone_week_rental_melted["日時"].dt.date
szone_week_rental_melted["時間"] = szone_week_rental_melted["日時"].dt.strftime("%H:%M:%S")

#各日付ごとに時間帯ごとに貸出数をグループ化して合計
daily_time_grouped = szone_week_rental_melted.groupby(["日付", "時間"])["貸出数"].sum().reset_index()

#各日付ごとに貸出数の合計を計算
daily_totals = daily_time_grouped.groupby("日付")["貸出数"].sum().reset_index()
daily_totals.columns = ["日付", "総貸出数"]

#各日付の時間ごとの貸出数割合を計算
daily_time_grouped = daily_time_grouped.merge(daily_totals, on="日付")
daily_time_grouped["貸出数割合"] = daily_time_grouped["貸出数"] / daily_time_grouped["総貸出数"]

#駅乗降客数の貸出数分布を用いた割振り計算
results = []
daily_time_grouped["日付"] = pd.to_datetime(daily_time_grouped["日付"])
for _, row in stations_passengers.iterrows():
    zone = row["小ゾーン"]
    daily_passengers = row["駅乗降客数"]

    for date in ["2024-09-06", "2024-09-09", "2024-09-10", "2024-09-11", "2024-09-12"]:
        # その日の日付に対応する割合を抽出
        day_ratios = daily_time_grouped[daily_time_grouped["日付"] == pd.to_datetime(date)]

        # day_ratios に対してループ処理を行う
        for _, ratio_row in day_ratios.iterrows():
            time = ratio_row["時間"]
            ratio = ratio_row["貸出数割合"]
            hourly_passengers = daily_passengers * ratio

            results.append({
                "小ゾーン": zone,
                "年月日": date,
                "時間帯": time,
                "駅乗降客数": hourly_passengers
            })

# 割り振り結果をデータフレームに変換
hourly_passengers_df = pd.DataFrame(results)

# 結果を表示
hourly_passengers_df[hourly_passengers_df["時間帯"] == "00:00:00"].head(20)

# 小ゾーンの型を変更
hourly_passengers_df["小ゾーン"] = hourly_passengers_df["小ゾーン"].astype("int64")



###【従業者数の欠損値確認（欠損値無し）】
#小ゾーンごとの従業者数CSV読み込み
sarea_employee = pd.read_csv(r"C:\\Users\\user\\Desktop\\sharecycle\\02.GIS\02.csv\szone_employees.csv")

#szone差分確認
set(szone_hour_rental_fil["szone"]) - set(sarea_employee["szone"])
len(set(sarea_employee["szone"]) - set(szone_hour_rental_fil["szone"]))

##不要なszoneの削除
sarea_employee = sarea_employee.merge(szone_hour_rental_fil[["szone"]], on="szone", how="right")
len(set(sarea_employee["szone"]) - set(szone_hour_rental_fil["szone"]))

##欠損値確認
#列名称変更
sarea_employee = sarea_employee.rename(columns={"szone": "小ゾーン", "JUGYOSHA": "従業者数"})
sarea_employee["従業者数"].isnull().sum()
#欠損値なし。



##【データの結合】
#従業者数DFと駅乗降客数DFの結合
hourly_passengers_df = hourly_passengers_df.merge(sarea_employee, on="小ゾーン", how="left")

#列名変更
hourly_passengers_df = hourly_passengers_df.rename(columns={"時間帯": "日時"})

##貸出数データ形式を変更
#sznoe列名変更
szone_hour_rental = szone_hour_rental.rename(columns={"szone": "小ゾーン"})

# データを「szone」「日時」「貸出数」の形式に変換
szone_hour_rental_melted = szone_hour_rental.melt(id_vars="小ゾーン", var_name="日時", value_name="貸出数")

#列名変更
szone_hour_rental_melted = szone_hour_rental_melted.rename(columns={"時間帯": "日時"})

##日時結合
# 「日時」列を datetime 型に変換
szone_hour_rental_melted["日時"] = pd.to_datetime(szone_hour_rental_melted["日時"])

# 「日時」列を datetime 型に変換してから「年月日」と「時間帯」に分割
szone_hour_rental_melted["日時"] = pd.to_datetime(szone_hour_rental_melted["日時"])
szone_hour_rental_melted["年月日"] = szone_hour_rental_melted["日時"].dt.date
szone_hour_rental_melted["日時"] = szone_hour_rental_melted["日時"].dt.strftime("%H:%M:%S")

hourly_passengers_df["年月日"] = pd.to_datetime(hourly_passengers_df["年月日"]).dt.date
szone_hour_rental_melted["年月日"] = pd.to_datetime(szone_hour_rental_melted["年月日"]).dt.date

hourly_passengers_df["日時"] = hourly_passengers_df["日時"].astype(str)
szone_hour_rental_melted["日時"] = szone_hour_rental_melted["日時"].astype(str)

#貸出数DFと従業者数DFを結合
merged_df = pd.merge(
    hourly_passengers_df,
    szone_hour_rental_melted,
    left_on=["小ゾーン", "年月日", "日時"],
    right_on=["小ゾーン", "年月日", "日時"],
    how="left"
)

# 「年月日」列を datetime 型に変換
merged_df["年月日"] = pd.to_datetime(merged_df["年月日"])

# 曜日を抽出して新しい「曜日」列を追加
merged_df["曜日"] = merged_df["年月日"].dt.day_name()

#年月日列とszone列を削除
merged_df = merged_df.drop(columns=["年月日"])

# 列の順番を並べ替える
merged_df = merged_df[["小ゾーン", "従業者数", "曜日", "日時", "駅乗降客数", "貸出数"]]

これで前処理が完了し必要なデータが揃いました。 続けてモデル構築をします。

【モデル構築】モデル構築とモデル評価

今回曜日（平日）と時間帯の2つをダミー変数化しました。
モデルは重回帰モデルを使用し、statsmodelからモデル評価を行いました。
以下がコードとstatsmodelの結果になります。

重回帰モデルの構築

#%%【ダミー変数化とモデル構築・評価】
#ダミー変数化　曜日、日時
dum_label = ["曜日", "日時"]
train_add_dum = pd.get_dummies(data=merged_df, columns=dum_label, dtype=int)

#説明変数と目的変数の取得
train_features = train_add_dum.drop(columns=["小ゾーン", "曜日_Friday", "日時_23:00:00", "貸出数"])
train_target = train_add_dum["貸出数"]

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

#分割
train_features_tr, train_features_va, train_target_tr, train_target_va = train_test_split(train_features, train_target, test_size=0.2, random_state=0)

#モデル構築と学習
model = sm.OLS(train_target_tr, train_features_tr).fit()

#モデル要約
print(model.summary())

■図9:貸出数重回帰モデル　statmodel結果

■考察

【モデル評価】決定係数（R値）

R値はモデルの精度を評価する1つの指標であり、1に近いほど精度が高いといえます。
今回0.242という値のため、1には遠く精度は高くない結果となりました。
今後説明変数の追加などで精度をあげる余地があります。

【モデル評価】回帰係数(coef)のプラス/マイナス方向

従業者数と乗降客数の回帰係数が0以上でプラスに向いています。
これは目的変数の貸出数に対して、プラスの影響を与えることを表しています。
つまり従業者数および乗降客数が増えれば、貸出数は増える傾向であると考えると、真っ当な考えになります。
曜日についてはマイナスに向いていますが、想定としては平日における曜日間の影響は少ないことから、説明変数として適切ではない可能性が考えられます。
時間帯について、深夜帯がマイナスに向き、日中がプラスに向いていることから、利用傾向に沿っていると考えられます。

【モデル評価】P値

P値は統計モデルにおいて0.05以下の説明変数が統計的に有意（意味がある）な変数であることを示します。
今回曜日以外の変数の殆どは統計的に有意であり、曜日は有意ではないということが言えます。

7.結論（再掲）と総括

貸出数の時系列グラフとヒートマップから、シェアサイクルの利用は主に通勤時間帯である7時～8時台と17時～18時台にピークが見られました。
田町や大手町・神田エリアなど、複数路線が交わり、駅間が数kmあるエリアで、ラストワンマイルの通勤手段として頻繁に利用されていることがわかります。
田町駅や品川駅の現場観察では、利用者の8割以上がスーツ姿のビジネスマンでした。
貸出数、駅乗降客数、従業者数データから重回帰モデルを構築しモデル評価を行いました。
今後は、潜在需要の予測と比較から、需要ギャップのあるエリアの分析を進めます。

今回分析した結果から、都心部において通勤で多く利用されているものの、都心部以外ではまだポートが少ないことや利用可能性のある人の母数が減ることから広く多く利用されている状況ではないことがわかりました。
電車・バス・タクシー等の交通は既に多くのルールがつくられていますが、シェアサイクルはまだまだルールが他交通に比較して整ってはいません。
ルール形成とポート密度向上、更には事業者のビジネスモデルが拡大できるモデルとして成立させることが、大きな交通課題への解決につながる道であることが考えられます。
今回は交通課題へ向かうまでの現状として都心部におけるビジネスマンの利用状況を明らかにしました。
引き続き予測モデルを構築しつつ、発展していくシェアサイクル市場について分析していきたいと思います。