分析② ロジスティック回帰分析を導入

 

分析① 距離帯別係数を用いた分析で得られたグラフ

前々回の記事「考え方と回帰モデルの基本形」では、「後発路線(=既存の鉄道網の中に、後から開通した路線)jの単独駅であるi駅」の乗車人員を推計する式の基本形として、以下の(1)のような式を立てました。

 
 乗車人員Pi = 定数C ×  住勤鉄道需要Ji ×  路線別係数Lj ×  関数f(d) × 誤差項εi …(1)
 
前回の記事「分析①:「距離帯別係数」を用いる」では、(1)の式の 関数f(d) の部分に「距離帯別係数」を導入して回帰分析を行い、以下のようなグラフを得ました。
 
図 分析①:関数f(d)に距離帯別係数を用いた際のグラフ
 

そして、 関数f のグラフの概形が、以下のようなものであることを確認しました。

  • 「他の路線の最寄駅までの距離」が遠いほど、数値が大きくなる。
  • 但し、2500m程度以上離れると傾きがほぼ0となり、それ以上には数値が大きくならない。
 
このうち1点目は、「他路線の駅から離れた方が、駅同士の需要を食い合いが少ない分、乗車人員が多くなりやすい」という傾向があることを、実証したものだと考えられます。
一方2点目は、「他路線の駅から離れれば離れるほど、無制限に乗車人員が多くなり続ける」というわけではなく、駅の需要には上限がある(どこかの地点で飽和する)、ということを示すものだと考えられます。
 

出力に上限がある関数として、ロジスティック関数を導入

dの値を大きくしつづけたときに、 関数f の出力値 関数f(d) がどこかで飽和するということであれば、関数fには、出力値に上限があるような回帰モデルを適用することが妥当でしょう。
 
そこで導入するのが、ロジスティック関数です。
ロジスティック関数は、以下のような特徴を持つ、S字型の曲線を描く関数です。
  • 入力値が小さい時には、出力値が0に限りなく近い。
  • 入力値を大きくしてゆくと、どこかで出力値が跳ね上がる。
  • しかし、出力値には上限がある。
  • 上限に近づくと、いくら入力値を増やしても、出力値は増えなくなる。
 
図 ロジスティック関数のグラフの例
 

ロジスティック関数の式

ロジスティック関数は、次のような式で表されます(eは自然対数の底)。
 
f(x) = C / (1 + e-A(x-B) …(2)
  • 上図のうち、黒の二重線のグラフは、A=1, B=0, C=1 とした時のグラフですが、この時、f(x)は0から1の間をとります。
  • Aが大きくなると、xを増やした際のf(x)の跳ね上がりの傾きが、急となります
  • Bが大きくなると、S字を描くグラフが右方向へずれてゆきます(S字曲線の中心点(変曲点)の座標が、 (x , y) = (B , 0.5C) となる)。
  • Cが大きくなると、下限値は0のまま、上限値がCとなり、グラフがタテ方向に引き伸ばされます
 

乗車人員の回帰式に、ロジスティック関数を組み込んでみる

(2)の式のうち、C(=ロジスティック関数の上限値)の役割は、(1)の式では、定数Cがこれを果たします。よって、関数f(d)に入れるロジスティック関数には、改めてCを入れる必要はありません。上限値を1とした、下記の式を用いることとなります。
 
 関数f(d) 1 / (1+ e-A(d-B))
 
したがって、(1)の回帰式は、次のような回帰式となります。
 
 乗車人員Pi = 定数C × 住勤鉄道需要Ji × 路線別係数Lj × 1 / (1+ e-A(d-B)) × 誤差εi …(3)
 

Excelのアドイン「ソルバー」を使って最適化

(3)のようなロジスティック回帰を組み入れた分析について、当初、統計解析用のフリーのプログラム言語「R言語」の使用を考えました。
しかし、実際に用いる関数やコマンドは、R初心者である私にとっては検索サイト頼りです。
 
まず、ロジスティック回帰にも使える関数として、「glm関数」を見つけたのですが、被説明変数に「成功する=1、失敗する=0」というようなゼロイチのデータを用いる際の関数であるようで、今回の分析には使えそうにありません。
他にも使えそうな関数を試したのですが、どうにもプログラムが回りません。
 
そこで、さらに検索を進めると、Excelのアドイン「ソルバー」が使えるらしいことを知りました。
ソルバーは、「様々な条件や制約を与えた時に、その制約の範囲内で、出力値を最大化・最小化するような入力値を探索する」という最適化を行なってくれるツールのようです。
 
これを使って、下記を満たすような 定数ABC 路線別係数Lj の組み合わせを探索しました。
  • 最適化の内容:(3)の回帰式のdに各駅のデータを入力したときに、乗車人員の理論値Pと実測値との誤差εが最小化される
  • 制約条件1:定数C>0
  • 制約条件2:それぞれの 路線別係数Lj >0
 誤差εの最小化には、(統計学的・数理的に正しい自信がないのですが)「乗車人員の対数の残差平方和∑εi2を最小化する」という、最小二乗法を用いました。
 

最適化を行う上での留意点

なお、この種の最適化では、ABC 路線別係数Lj に初期値を与えると、初期値をスタートとして計算を開始し、ABC 路線別係数Lj の値を少しずつ動かしながら、最小化させたい値(ここでは残差平方和∑εi2)が最小化されるポイントを探索する、という計算手順がとられるようです。
 
但し、与える初期値が異なると、異なる探索結果が出てきてしまう場合もあるようです。ある値を最大化させる分析を登山にたとえるならば、「最高峰に登りたいのに、登山口を間違うと、別の山の山頂に連れていかれる」というようなイメージでしょうか。
 
そこで今回は、3通りほどの初期値パターンを与え、3パターンでほぼ同じ値に収束することを確認した上で、最終結果としました。

 

  分析② ロジスティック回帰分析の結果

 

A,B,Cの値と回帰式

ソルバーを回した結果、(1)の回帰式の係数は、A = 0.01478B = 1425.2C = 1.3914 となりました。

 

従って、得られた回帰式は、以下のようになります。

 

 乗車人員Pi 1.3914 × 住勤鉄道需要Ji × 路線別係数Lj × 1 /(1 + e-0.001478(d -1425.2)) × 誤差εi

 

距離に関する関数 f(d) のグラフ

距離に関する 関数f(d) = 1 /(1 + e-0.001478(d -1425.2)) のグラフを描くと、以下のようになります。

 

図 分析②:関数f(d)にロジスティック関数を用いた際のグラフ

 

最大値1に対して、他路線の最寄駅からの距離が500m地点では0.203、1000m地点では0.348であるのに対し、1425.2m(=B)地点で0.5を超え、3000m地点では0.9を上回る0.911となっています。

 

分析①と②との比較 ―路線別係数Lj

ここでは、今回の分析②で得られた19路線それぞれの路線別係数Ljを表にまとめ、分析①で得られた路線別係数と比較します。

路線別係数(表のうち黄色のグラフの欄)は、便宜的にJR埼京線の係数を1.0とおいて算出したものに過ぎません。そこで、正確な背比べのため、分析対象179駅の係数の相乗平均値を、分析①と分析②とで1.0に揃えてから比較します(表のうち白地の欄)。

 

表 分析②と①の路線別係数の比較

 

結果を見ると、分析①と分析②とで、ほぼ同様の結果が出ていることが分かります。

 

分析①と②との比較 ―距離に関する関数 f(d) のグラフ

ここでは、

分析①による回帰式:

 乗車人員Pi = 0.221 × 住勤鉄道需要Ji × 路線別係数Lj × 距離帯別係数Zk × 誤差εi

 

分析②による回帰式:

 乗車人員Pi = 1.3914 × 住勤鉄道需要Ji × 路線別係数Lj × 1/(1 + e-0.001478(d - 1425.2)) × 誤差εi

 

の両者の関数fのグラフ同士を比較してみます。

 

分析①における関数fに当たるのは 距離帯別係数Zk ですが、これは「他路線の最寄駅からの距離が0~500m以内」の係数を基準値=1.0とした時に、他の距離帯の係数がどのような比率となるかを示したものにすぎません。そのため、分析②における関数fのグラフ(出力値が0~1の間をとるロジスティック関数)とは、グラフの高さが異なることとなります。

 

グラフの高さを合わせるため、分析①による回帰式の定数(=0.221)が分析②による回帰式の定数(=1.3914)と等しくなるよう、また 路線別係数Lj の179駅の相乗平均値(=0.470)が 路線別係数Lj の相乗平均値(=0.489)と等しくなるように式変形すると、分析①による回帰式は、以下のようになります。

 

 乗車人員Pi 1.3914 × 住勤鉄道需要Ji × 分析②の路線別係数Lj × ((0.221/1.3914)×(0.470/0.489)× 距離帯別係数Zk × 誤差εi

 

 乗車人員Pi 1.3914 × 住勤鉄道需要Ji × 分析②の路線別係数Lj × 0.1525× 距離帯別係数Zk × 誤差εi

 

そこで、分析①の 0.1525× 距離帯別係数Zk と、分析②のロジスティック回帰分析 1/(1 + e-0.001478(d - 1425.2)とを比較してみることとします。

 

その結果は、以下の通りです。

 

 

分析②のグラフは、分析①のグラフと比べて、2000~2500m区間でやや高く、2500~3000m区間でやや低くなっていますが、おおむね両者は同じような軌跡を描いていることが分かります。

分析①において、2500mを境に落差が大きくなっていた箇所を、分析②では、滑らかな曲線でつないでいる、という言い方もできるでしょう。

 

分析①と②との比較 ―残差平方和

分析によって得た残差平方和は、分析①が36.7059であるのに対し、分析②が37.9807となりました。

 

グラフの形に特段の前提をおかずに(制限をかけずに)分析を行った①の方が残差平方和を小さく抑えられていますが、ロジスティック関数を当てはめた②も、当てはまりの良さでは遜色がないといえそうです。

 

完成した回帰式と、モデルの当てはまり状況

分析①の時と同様に、分析②についても、得られた関数fのグラフの当てはまりのよさを確認するために、

  • 横軸…「他の路線の最寄駅までの距離」
  • 縦軸…「定数・路線別係数で補正した各駅の鉄道需要・乗車人員比 Ri /(C × Lj)  = 関数f(d) × 誤差εi

を置き、これを 関数fのグラフと重ね合わせてみます。

 

図 分析②による関数fと補正済み鉄道需要・乗車人員比の重ね合わせ

 

分析①の散布図と同様に、各駅の分布のバラツキは抑えられ、その分布のほぼ中央を通る形で、関数fのグラフが描かれています。