見出し画像

【Pythonで統計モデル】Rライクにモデルを書く:Patsy

はじめに

株式会社GA technologies(以下GA)の福中です。今、僕は、AISCという部署でChief Data Scientistとして働いています(詳しい自己紹介はAISC WEBサイトのプロフィールをご覧ください)。今回はPythonでRライクにモデルを書いて分析する方法についてご紹介したいと思います。


なぜこの記事を書こうと思ったのか?

僕はもともとデータ解析にはRを使っていたのですが、GAに来てからはPythonを使うようになりました。昔、Rを使っていた理由は、こと統計モデリングに限って言えばPythonの出力が貧弱だったからです。このあたりの経緯は以下の記事をご参照ください。

しかし時代が進み、現在ではPyhonでも問題なく統計モデルが構築でき、統計学者でも満足できる出力が得られるようになりました。とはいえ、若干の不満点もありました。それは統計モデルにかけるまでのデータセットの準備に多少の手間があったことです。例えば、回帰モデルにおいて、切片項を入れる際に固定変数1を追加したり、カテゴリカル変数を扱う際にはダミー化する必要があったり、あるいは交差項を入れたいときには手動での変数の準備が必要だったりすることなどです。

Rではこれらのことはある程度自動でやってくれるので、その感覚があるとどうしてもPythonで統計分析をすることに面倒さを感じてしまいます。出力結果がリッチになった現在でも、統計学界隈でRがよく選択される理由の1つがこのあたりにあるのではないかと思っています。良くも悪くもPythonは「プログラミング言語寄り」なのです。

このようにちょっとしたことではあるのですが、このラスト1マイルのせいでPythonによる気軽な統計モデリングが阻害されているのではないかと思いました。これはとても残念なことです。そこで、ある程度この面倒さを緩和してくれるPatsyというライブラリを紹介する記事を書こうと思いました。

テキスト

本記事では岩波書店から出版されている以下の書籍「データ解析のための統計モデリング入門-一般化線形モデル・階層ベイズモデル・MCMC-」の3章で行われていた分析を例にとって解説していきます。

分析用データは著者である久保拓弥先生のサポートページからダウンロードできますので、興味のある方はぜひ試してみてください。

環境

本記事の分析は、Dockerの「datascience-notebook:x86_64-ubuntu-22.04」というイメージを使って環境を構築し、その中で行いました。環境構築に興味のある方は、以下のDS Lab(データサイエンスの実験的環境)に関する記事をご参照ください。

また、今回使用するライブラリは以下です。あらかじめ読み込んでおきます。

import numpy as np
import pandas as pd
import statsmodels.api as sm
from patsy import dmatrices

一般化線形モデルによる分析

Pythonで一般化線形モデルを利用する場合、statsmodelsが使いやすく、よく利用されます。例えばテキスト第3章のデータ(data3a.csv)を使って、説明変数が2つ(数量型+因子型)のポアソン回帰モデルを構築してみましょう。手順は以下の通りです。

まずはデータのCSVファイルを読み込みます。

# データの読み込み
d = pd.read_csv('data3a.csv')

次にデータフレームから目的変数だけを取り出します。

# 目的変数の作成
outcome = d.y

そして説明変数の作成を行います。次のように因子型の変数がTの場合1、Cの場合0に変換し、ダミー変数として追加します。また、切片項のために定数1の列を追加しておきます。

# 説明変数の作成
# f列の値をTの場合1に、Cの場合0に変換し、fT列として追加
d['fT'] = pd.get_dummies(d.f)['T'].astype(int)

# x列とfT列のみ取り出して分析用の説明変数を作成
predictor = d[["x", "fT"]]

# 切片項のための1の定数列を追加
predictor = sm.add_constant(predictor)

print(predictor)
    const      x  fT
0     1.0   8.31   0
1     1.0   9.44   0
2     1.0   9.50   0
3     1.0   9.07   0
4     1.0  10.16   0
..    ...    ...  ..
95    1.0   9.15   1
96    1.0   8.52   1
97    1.0  10.24   1
98    1.0  10.86   1
99    1.0   9.97   1

これで準備が整ったので、後はリンク関数を設定して、statsmodelsのGLMメソッドを実行するだけです。

# GLMを適用
# 対数リンク関数を指定してポアソン回帰を実行

link = sm.families.links.Log()
family = sm.families.Poisson(link=link)
model = sm.GLM(outcome, predictor, family=family)
res = model.fit()
print(res.summary())

結果は以下の通り、テキストと解が一致していることがわかると思います。

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       97
Model Family:                 Poisson   Df Model:                            2
Link Function:                    Log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -235.29
Date:                Fri, 12 Apr 2024   Deviance:                       84.808
Time:                        05:56:51   Pearson chi2:                     83.8
No. Iterations:                     4   Pseudo R-squ. (CS):            0.04590
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.2631      0.370      3.417      0.001       0.539       1.988
x              0.0801      0.037      2.162      0.031       0.007       0.153
fT            -0.0320      0.074     -0.430      0.667      -0.178       0.114
==============================================================================

このように、一般化線形モデル自体は簡単に適用でき、また結果の出力も十分な情報が得られているのですが、Rに比べるとデータの準備が大変、、、というか面倒なのです。

Patsyについて

そこで、この面倒な手順を簡略化してくれるライブラリがPatsyです。このPatsyは、実際のところstatsmodelsのformula.api(statsmodels.formula.api)の背後で動いているライブラリなので、単体で使用することはあまりないかもしれません。しかし、今回は説明のため、あえてこのライブラリを取り上げ、実行してみたいと思います。ただし、網羅的な解説はしませんので、詳しいことを知りたい場合は以下のマニュアルをご覧ください。

具体例

何はともあれ、仕様の解説より前に、まずは具体例を見ていきましょう。まずはデータの読み込みを行います。

# データの読み込み
d = pd.read_csv('data3a.csv')

次にPatsyのdmatrices()を使って、Rのformulaオブジェクトのようにモデル式を書いて設定します。

# Patsyを使ってRライクにモデル式を書く
outcome, predictor = dmatrices('y ~ x + C(f)', data=d, return_type='dataframe')

そしてstatsmodelsのGLMを適用します。

# GLMを適用
# 対数リンク関数を指定してポアソン回帰を実行

link = sm.families.links.Log()
family = sm.families.Poisson(link=link)
model = sm.GLM(outcome, predictor, family=family)
res = model.fit()
print(res.summary())

結果は以下のようになり、先ほど(テキストと同様)の出力が得られていることがわかります。

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       97
Model Family:                 Poisson   Df Model:                            2
Link Function:                    Log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -235.29
Date:                Mon, 15 Apr 2024   Deviance:                       84.808
Time:                        00:23:21   Pearson chi2:                     83.8
No. Iterations:                     4   Pseudo R-squ. (CS):            0.04590
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      1.2631      0.370      3.417      0.001       0.539       1.988
C(f)[T.T]     -0.0320      0.074     -0.430      0.667      -0.178       0.114
x              0.0801      0.037      2.162      0.031       0.007       0.153
==============================================================================

このようにデータを準備する際に、切片項のための変数を追加したり、カテゴリカルデータをダミー変数に変換したりすることなく、

dmatrices('y ~ x + C(f)', data=d, return_type='dataframe')

の1行でGLMの分析まで持って行けるのです!なんて便利なのでしょう!また、第1引数の「'y ~ x + C(f)'」を見てもらえばわかるように、Rのformulaライクにモデルを定義することが可能になっています。カテゴリカル変数の場合は「C(列名)」とするだけで指定可能です。

Patsyでできる代表的な例

それではPatsyでどのようなことができるのか見ていきましょう。様々なことができるのですが、ここでは代表的な例を示します。

交差項を入れる
交差項を入れたい場合、コードを以下のように書きます。

# 交差項を入れる
outcome, predictor = dmatrices('y ~ x + C(f) + x:C(f)', 
                               data=d, return_type='dataframe')

モデル式に「x:C(f)」を追加するだけですね。この辺はRと同じです。ちなみにRと同様に「*」を使って、以下のように書くことも可能です。この場合、xの主効果、C(f)の主効果、xとC(f)の交互作用のすべてをモデルで表現することになります。

# 交差項を入れる
outcome, predictor = dmatrices('y ~ x * C(f)', 
                               data=d, return_type='dataframe')

またRの場合、全変数を指定するときに「.(ピリオド)」が使用できますが、Patsyでは使えないのでご注意ください

切片項を省略
モデルの中に切片項を表現したくない場合は、以下のようにモデル式に「-1」を入れます。

# 切片項を省略
outcome, predictor = dmatrices('y ~ x + C(f) - 1', 
                               data=d, return_type='dataframe')

変数変換(対数変換)
特定の変数を変数変換してモデル表現することが可能です。例えば対数変換したい変数がある場合、「np.log(変数名)」のように記述します。

# 変数変換(対数変換)
outcome, predictor = dmatrices('y ~ np.log(x) + C(f)', 
                               data=d, return_type='dataframe')

変数変換(標準化・中心化)
特定の変数を標準化してモデル表現することが可能です。「standardize(変数名)」のように記述します。

# 変数変換(標準化)
outcome, predictor = dmatrices('y ~ standardize(x) + C(f)', 
                               data=d, return_type='dataframe')

ちなみに平均のみ調整する中心化を行うことも可能です。

# 変数変換(中心化)
outcome, predictor = dmatrices('y ~ center(x) + C(f)', 
                               data=d, return_type='dataframe')

変数変換(独自関数)
独自に定義した関数を使ってモデル式に組み込むことも可能です。

# 変数変換(独自間数)
def linear(x):
  return 2 * x + 1

outcome, predictor = dmatrices('y ~ linear(x) + C(f)', 
                               data=d, return_type='dataframe')

dmatricesの仕様

最後にdmatricesの仕様をまとめておきます。なお、ここで解説するのはPastyのバージョン0.5.1に基づくものになります。また、dmatricesの仕様はdmatrixと同等になりますが、返り値がoutocomeとpredictorの2つになっている点が異なります。詳しくはマニュアルを直接ご参照ください。

dmatrices(formula_like, 
          data={}, 
          eval_env=0, 
          NA_action='drop', 
          return_type='matrix')

formula_like
モデル式を構築するために使用するオブジェクトです。基本的にRのformulaオブジェクトと同等の働きをしますが、完全に一致しているわけではないことに注意が必要です。基本的(代表的)な使い方は上記の説明をご参照ください。より詳しく知りたい方はマニュアルをご参照ください。

data
formula_like で参照されている変数を検索するために使用されるデータを指定します。

eval_env
デフォルトではeval_env=0に設定されており、この場合、formula_likeの中で参照されている変数をdataの中から探索することを意味しています。もし本関数を外部のライブラリから呼び出したい場合はeval_env=1に設定する必要があります。

NA_action
欠損値を含む行をどのように処理するかを指定します。それらをドロップしたり、エラーを発生させたり、カスタマイズのためにNAActionオブジェクトを渡したりすることができます。

return_type
返り値を"matrix"あるいは"dataframe" のどちらで取得したいかを指定します。

おわりに

RにはRの、PythonにはPythonのメリットデメリットがあるので、一概にどちらを選ぶべきかはその時の状況によると思います。ただ、お互いにそれぞれの良いところを吸収するのは、ユーザーにとってはありがたいことなので、RユーザーからPythonユーザーになった僕にはPatsyは大変助かりました。同じようにPythonで統計モデリングを行うことに苦手意識をもっている方もいると思いますので、本記事がそのような方の参考になればうれしく思います。


この記事が気に入ったらサポートをしてみませんか?