Sageでグラフを再現してみよう：データ解析のための統計モデリング入門第６章

この企画は、雑誌や教科書にでているグラフをSageで再現し、 グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。

前回に続き、データ解析のための統計モデリング入門 （以下、久保本と書きます） の第6章の例題をSageを使って再現してみます。

数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるのが大きな特徴です。 この機会にSageを分析に活用してみてはいかがでしょう。

前準備

最初に必要なライブラリーやパッケージをロードしておきます。

上限のあるカウントデータの回帰分析

6章の最初に二項分布を使った回帰分析の例がでてきます。ポアソン分布では上限のないカウントデータでしたが、 上限のある場合には、二項分布を使うのだと説明がありました。 （自然界ではこのように上限のあるカウントデータが多いので、とても参考になりました）

二項分布用のデータについて

二項分布用のデータは、サポートページにあるdata4a.csvです。 これを読み込んで、データの性質と分布を表示します。

直感的に施肥を施す(f=T)ことでサイズxが大きくなっているのが分かります。 また、サイズと種子数の右上がりの関係が見られます。（種子数の上限は8です。）

二項分布

ここで、Sageのプロット関数を使って二項分布の確率分布がどのような形になっているか 久保本の図6.3と同じ条件でブロットしてみます。

二項分布を $$ p(y | N, q) = \binom{N}{y} q^y (1 -q)^{N-y} $$ 以下の様に_p関数で定義します。

Sageでプロットしてみます。Sageではグラフ（Graphics）にグラフを足すことで、重ね合わせのプロットができます。 とても直感的にグラフを書くことができます。（ggplotも同じです）

ロジスティック関数

ロジスティック関数 $$ q_i = logistic(z_i) = \frac{1}{1 + exp(-z_i)} $$ の関数を_logisticとして定義し、その分布をSageを使ってプロットします。

ロジスティック関数のような曲線を持つデータの場合には、二項分布の確率分布があると予測して回帰分析を行います。 リンク関数は、ロジスティック関数の逆関数であるロジットリンク関数を指定します。

ロジット関数は、以下の様に表されます(式a)。 $$ logit(q_i) = log \frac{q_i}{1 - q_i} $$

パラメータ推定

残念ながら、statmodelsではcbindに相当する二項分布の解析方法が分からず、 Rを使って計算することにしました。

Sageの中からRの関数を使うことができるので、計算を中断することなく進めることができます。

回帰の結果、$\beta_1 = -19.536, \beta_2 = 1.95, \beta_3=2.02$と求まりました。

結果の図化

ggplot2のstat_smoothで結果の曲線を表示しようとしたのですが、上手くできませんでした。

そこで、Sageのプロット機能を使って結果を表示してみました。

結果の解釈

久保本のすごいところは、単に回帰をすることで終わらず、その解釈を説明しているところです。 リンク関数がロジット関数である(式a)であることから $$ log \frac{q_i}{1 - q_i} = \beta_1 + \beta_2 x_i + \beta_3 f_i $$ となり、両辺の指数を取ると、 $$ \frac{q_i}{1 - q_i} = exp(\beta_1 + \beta_2 x_i + \beta_3 f_i) = exp(\beta_1) exp(\beta_2 x_i) exp(\beta_3 f_i) $$ となります。$\frac{q_i}{1 - q_i}$がオッズと呼ばれる量です。

施肥の有無による違いは、exp(2.02)=7.5倍であると推定されました。

割り算を使わない方法はすごい！

久保本の6章で感動したのは、割り算を使わない方法です。 私のようなものは、割り算を使って正規化してしまいますが、 このように割り算をしないで、回帰分析する方法としてオフセット項の使い方を 紹介しています。

オフセット用のデータには、data4b.csvを使います。 このデータは、

• 調査地iの面積を$A_i$
• 調査地iの明るさを$x_i$
• 調査地iにおける植物個体数$y_i$

オフセット項は、offset=np.log(オフセットのカラム)のように指定します。ここでは、オフセット値Aのlogを 取った値が使われていますが、これは以下のような理由からです。

人口密度を以下の様に表し、 $$ \frac{平均個体数\lambda_i}{A_i} = 人口密度 $$ 平均個体数と明るさに指数関数的な関係があると仮定すると、 $$ \lambda_i = A_i \times 人口密度 = A_i exp(\beta_1 + \beta_2 x_i) $$ ここで、面積を掛けている部分をlogを使うことで $$ \lambda_i = exp(\beta_1 + \beta_2 x_i + log A_i) $$ と変形することができます。つまり、$log A_i$だけずれた値（オフセット項）、久保本では「げた」をはかせたと考えると よいとありました。

推定値

推定された値をプロットすると以下のようになります。 データよりもはっきりと明るさによる影響がきれいに分かります。（そのようにモデルを作ったので当たり前なのですが）

オフセット項を使わないと

オフセット項を使わないと面積による影響が回帰分析に含まれないため、予測された結果は、ぼんやりとしたものになります。

正規分布と尤度

Sageにも確率分布を扱う関数RealDistirubtionがVer.5から導入されました。（まだ一部の関数のサポートですが） これを使って、正規分布の確率密度関数をプロットしていました。

確率密度関数から確率を得る

確率密度関数から指定された範囲に含まれる確率は、その範囲の面積になります。（その範囲で積分した値） cum_distribution_functionは、マイナス無限大から指定された値までの確立を計算するので、 この関数の差を計算することで、指定範囲の確率（面積）を求めることができます。

正規分布の尤度

久保本に従って正規分布の尤度をと最尤尤度を求めてみます。 $y_iがy_i - 0.5\Delta y \leq y \leq y_i + 0.5\Delta y$である確率は、 $$ \begin{eqnarray} L(\mu, \sigma) & = & \prod_i p(y_i | \mu, \sigma) \Delta y \\ & = & \prod_i \frac{1}{\sqrt{2 \pi \sigma^2}} exp\left\{ \frac{(y_i - \mu)^2}{2\sigma^2} \right\} \Delta y \\ & = & \prod_i (2 \pi \sigma^2)^{-\frac{1}{2}} exp\left\{ \frac{(y_i - \mu)^2}{2\sigma^2} \right\} \Delta y \end{eqnarray} $$ 両辺をlogを取って $$ \begin{eqnarray} log L(\mu, \sigma) & = & -\frac{1}{2} \sum_i log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - \mu)^2 + \sum_i log(\Delta y) \\ & = & -\frac{1}{2} N log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - \mu)^2 + N log(\Delta y) \end{eqnarray} $$ $\Delta y$の項を無視すると、 $$ log L(\mu, \sigma) = -\frac{1}{2} N log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - \mu)^2 $$ となります。これを$\mu$で偏微分し、最小となるのは、$y = \mu$と求まります。（最小自乗法の解と同じ）

ガンマ分布

どんなときにガンマ分布を使うのかGoogleでみると、正の値で、寿命や待ち時間など経過時間に無関係一定値の場合に、 使われるとありました。

ガンマ分布をSageでプロットしてみます。

ガンマ分布を使った回帰分析

例題では、ある個体の花の重量を$y_i$が平均$\mu_i$のガンマ分布にしたがっていると仮定します。 $$ \mu_i = A x_i^b $$ で、A=exp(a)とすると、 $$ \mu_i = exp(a)x_i^b = exp(a + log x_i^b) = exp(a + b log x_i) $$ と表され、リンク関数は、expの逆関数のlogとなり、 $$ log \mu_i = a + b log x_i $$ となります。

ガンマ分布の例題データは、RData形式なので、Rで解析します。family=Gamma(link="log")でリンク関数logのガンマ分布を指定します。

glmを使うと複雑な分布の回帰が簡単に計算できます。しかし、どのモデルが良いかは、 プロット結果を見ながら判断するのがよいと実感しました。人間の目で分からない違いは、尤度やAICを使うのが現実できなのでは？