Data Science by R and Python

統計学を、広く、深く、わかりやすく。

pythonで一般化線形モデルと公式リファレンスの罠 - using statsmodels -

寒くなってきました。最近、pythonでデータの解析をすることにいそしんでおります。
Rでできることをpythonでやりたいなと思っていろいろ調べてみると、まぁなかなかできるようになっていなかったりするわけで、その辺を整備し始めたので、ここに書いていこうと思います。

一般化線形モデル(Generalized Linear Model)

一般化線形モデルは、回帰の問題ではよく利用される方法で、Rではglmという関数で使うことができます。有名なものはロジスティック回帰です。図的には、以下のような感じになります。

f:id:tomoshige_n:20141107184034p:plain

これをpythonで行うときには、多くの場合scikit-learnを使うような気がするのですが、scikit-learnのアウトプットは、統計的に結果を見たい場合に十分とはいえません。できればRのロジスティック回帰の結果みたいなの(下図)、、、を出力したいんです。

f:id:tomoshige_n:20141107185723p:plain

statsmodelsを使ってみよう。

そこで、そんな要望に答えるために、statsmodelsというモジュールが提供されています。どうもこれを使用すれば、Rのglm的なコトができるらしいと聞きつけて、やってみました。以下のコードをのせているので興味のある人は是非やってみて欲しいと思います。最後にRの対応するコードも付けておきました。使用しているデータはstatsmodelsの中のものですが、以下のリンクからもダウンロードできるようにしました。
https://dl.dropboxusercontent.com/u/2133906/sample.csv
※最後に注意事項も書いていますので、そこもできれば、読んでください。

## 一般化線形モデル
from __future__ import print_function
import numpy as np
import pandas as pd
from pandas import DataFrame, Series
from matplotlib import pyplot as plt

#glmを可能にするために必要なものなど
#statmodelsのGLMを使えるようにするモジュール
import statsmodels.api as sm
#GLMの中で用いる統計処理関数のインポート, chi**2の値などの算出に使用している 
from scipy import stats 
#Rのglmを使用可能にするための関数
import statsmodels.formula.api as smf 

# ## OLS, GLM: Gaussian response data

#sampleデータのロード
#形式はpandas.DataFrame, naは削除する必要あり
df = sm.datasets.get_rdataset("Guerry", "HistData").data
df = df[['Lottery', 'Literacy', 'Wealth', 'Region']].dropna()
df.head()

#線形回帰モデルの構築, 最小自乗法を使う場合
#modelの構築
ols_model = 'Lottery ~ Literacy + Wealth + Region'
#モデルの構築
mod = smf.ols(formula=ols_model, data=df)
#fitに回帰した結果が入っているので、これをresに代入する
res = mod.fit()
#結果の参照
res.summary()

"""
                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                Lottery   No. Observations:                   85
Model:                            GLM   Df Residuals:                       78
Model Family:                Gaussian   Df Model:                            6
Link Function:               identity   Scale:                   436.433771252
Method:                          IRLS   Log-Likelihood:                -375.30
Date:                Fri, 07 Nov 2014   Deviance:                       34042.
Time:                        18:59:36   Pearson chi2:                 3.40e+04
No. Iterations:                     4                                         
===============================================================================
                  coef    std err          z      P>|z|      [95.0% Conf. Int.]
-------------------------------------------------------------------------------
Intercept      38.6517      9.456      4.087      0.000        20.118    57.186
Region[T.E]   -15.4278      9.727     -1.586      0.113       -34.493     3.637
Region[T.N]   -10.0170      9.260     -1.082      0.279       -28.167     8.133
Region[T.S]    -4.5483      7.279     -0.625      0.532       -18.815     9.718
Region[T.W]   -10.0913      7.196     -1.402      0.161       -24.195     4.013
Literacy       -0.1858      0.210     -0.886      0.376        -0.597     0.225
Wealth          0.4515      0.103      4.390      0.000         0.250     0.653
===============================================================================
"""

#線形回帰モデルの構築, glmによるニュートン法を用いる場合
#modelの構築
glm_gauss_model = 'Lottery ~ Literacy + Wealth + Region'
#モデルの構築
mod = smf.glm(formula=glm_gauss_model, data=df, family=sm.families.Gaussian())
#fitに回帰した結果が入っているので、これをresに代入する
res = mod.fit()
#結果の参照
res.summary()


# ## GLM: Binomial response data
# 
# ### Load data
# 
#  In this example, we use the Star98 dataset which was taken with permission
#  from Jeff Gill (2000) Generalized linear models: A unified approach. Codebook
#  information can be obtained by typing: 

#一応データセットに関するディスクリプションを参照しておく。
print(sm.datasets.star98.NOTE)

"""
簡単に言えば、目的変数に「数学の成績が国の平均以上」だった人と、
「以下だった人」の数がたぶん学校毎に入っていて、それを学校に関係する
その他の変数で説明しましょうってことだと思います。
"""

#データセットの読み込み
star98 = sm.datasets.star98.load()

#pandas形式にしたい場合は
#star98 = sm.datasets.star98.load_pandas().data

"""
大事:
本来なら上記のロードだけで良いのだが、、、
pythonのstatsmodelsの通常のGLM = sm.GLMは
どういうわけか、intercept項を入れないまま回帰するので
sm.GLMを用いる場合には定数項を入れる必要があります。
"""
#data.exog = sm.add_constant(data.exog, prepend=False)

#datasetの形式を確認
type(star98.data)

# >>> numpy.core.records.recarray
#データセットをpandasのDataFrameにしておく。
star98 = DataFrame(star98.data)

#データセットの最初の5件を確認
star98.head()

#一般化線形モデルの構築
#2項回帰モデル、link:logit

"""
Binomialでモデリングする際には、
通常のRでは目的変数にcbind(成功数, 失敗数)を取るようにする。
これをpythonのstatsmodelsでは、~の前に以下のように記述する
'成功数 + 失敗数 ~ 説明変数'

"""

#モデルを定義する。
formula_binom = 'NABOVE + NBELOW ~  LOWINC+ PERASIAN + PERBLACK + PERHISP + PCTCHRT + PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF'
#成功割合の変数を作成する
star98['SUCCESS'] = star98['NABOVE'] / (star98['NABOVE'] + star98['NBELOW'])
star98['TOTAL'] = star98['NABOVE'] + star98['NBELOW']


# binomialの回帰モデルを構築する
glm_binom = smf.glm(formula=formula_binom, data=star98, family=sm.families.Binomial())
# res_binomにglmの結果を保存する
res_binom = glm_binom.fit()
# summaryを表示
res_binom.summary()

                  Generalized Linear Model Regression Results                   
================================================================================
Dep. Variable:     ['NABOVE', 'NBELOW']   No. Observations:                  303
Model:                              GLM   Df Residuals:                      282
Model Family:                  Binomial   Df Model:                           20
Link Function:                    logit   Scale:                             1.0
Method:                            IRLS   Log-Likelihood:                -2998.6
Date:                  Fri, 07 Nov 2014   Deviance:                       4078.8
Time:                          19:58:26   Pearson chi2:                 4.05e+03
No. Iterations:                       7                                         
============================================================================================
                               coef    std err          z      P>|z|      [95.0% Conf. Int.]
--------------------------------------------------------------------------------------------
Intercept                    2.9589      1.547      1.913      0.056        -0.073     5.990
LOWINC                      -0.0168      0.000    -38.749      0.000        -0.018    -0.016
PERASIAN                     0.0099      0.001     16.505      0.000         0.009     0.011
PERBLACK                    -0.0187      0.001    -25.182      0.000        -0.020    -0.017
PERHISP                     -0.0142      0.000    -32.818      0.000        -0.015    -0.013
PCTCHRT                      0.0049      0.001      3.921      0.000         0.002     0.007
PCTYRRND                    -0.0036      0.000    -15.878      0.000        -0.004    -0.003
PERMINTE                     0.2545      0.030      8.498      0.000         0.196     0.313
AVYRSEXP                     0.2407      0.057      4.212      0.000         0.129     0.353
PERMINTE:AVYRSEXP           -0.0141      0.002     -7.391      0.000        -0.018    -0.010
AVSALK                       0.0804      0.014      5.775      0.000         0.053     0.108
PERMINTE:AVSALK             -0.0040      0.000     -8.450      0.000        -0.005    -0.003
AVYRSEXP:AVSALK             -0.0039      0.001     -4.059      0.000        -0.006    -0.002
PERMINTE:AVYRSEXP:AVSALK     0.0002   2.99e-05      7.428      0.000         0.000     0.000
PERSPENK                    -1.9522      0.317     -6.162      0.000        -2.573    -1.331
PTRATIO                     -0.3341      0.061     -5.453      0.000        -0.454    -0.214
PERSPENK:PTRATIO             0.0917      0.015      6.321      0.000         0.063     0.120
PCTAF                       -0.1690      0.033     -5.169      0.000        -0.233    -0.105
PERSPENK:PCTAF               0.0490      0.007      6.574      0.000         0.034     0.064
PTRATIO:PCTAF                0.0080      0.001      5.362      0.000         0.005     0.011
PERSPENK:PTRATIO:PCTAF      -0.0022      0.000     -6.445      0.000        -0.003    -0.002
============================================================================================

# ### 興味のある統計量を表示する方法
#試行回数
print("Total number of trials:" + str(star98[['NABOVE','NBELOW']].sum().sum(axis=1)))
#パラメータの表示
print('Parameters: ', res_binom.params)
#各パラメータに対するt値
print('T-values: ', res_binom.tvalues)



# ### Plots
# ロジスティック回帰の結果について図で確認する。
# プロットを通して、データの予測精度について分析する
# ・予測値と実測値のプロット
# ・予測値と残差のプロット


#観測値
y = star98['SUCCESS']
#推定値
yhat = Series(res_binom.mu)
#linear_predictor
y_lp = Series(np.log(res_binom.mu/(1-res_binom.mu)))
#残差
residuals = Series(res_binom.resid_pearson)

#dataframe:
dd = DataFrame({
	'OBS':y,
	'FITTED':yhat,
	'LINEAR_PRED':y_lp,
	'RESIDUALS':residuals
	})

# 予測値と観測値の差を確認する
from ggplot import *
p_fitted_vs_obs = ggplot(aes(x='FITTED',y='OBS'),data=dd)+geom_point(colour="red")+geom_abline(intercept=0,slope=1)
# 予測値と残差の関係性を確認する
p_fitted_vs_residuals = ggplot(aes(x='FITTED',y='RESIDUALS'),data=dd)+geom_point(colour="red")+stat_smooth()

f:id:tomoshige_n:20141107202424p:plain

Rのコード

以下のようになります。よって、結果が一致しました!

> star98 = read.csv("sample.csv")
> 
> glm_binom = glm(cbind(NABOVE,NBELOW) ~ LOWINC + PERASIAN + PERBLACK + PERHISP + PCTCHRT + PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF, data=star98, family="binomial")
> summary(glm_binom)

Call:
glm(formula = cbind(NABOVE, NBELOW) ~ LOWINC + PERASIAN + PERBLACK + 
    PERHISP + PCTCHRT + PCTYRRND + PERMINTE * AVYRSEXP * AVSALK + 
    PERSPENK * PTRATIO * PCTAF, family = "binomial", data = star98)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-9.5217  -2.6505  -0.3357   2.1658  15.4200  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               2.959e+00  1.547e+00   1.913   0.0557 .  
LOWINC                   -1.682e-02  4.340e-04 -38.749  < 2e-16 ***
PERASIAN                  9.925e-03  6.014e-04  16.505  < 2e-16 ***
PERBLACK                 -1.872e-02  7.436e-04 -25.182  < 2e-16 ***
PERHISP                  -1.424e-02  4.339e-04 -32.818  < 2e-16 ***
PCTCHRT                   4.917e-03  1.254e-03   3.921 8.81e-05 ***
PCTYRRND                 -3.580e-03  2.255e-04 -15.878  < 2e-16 ***
PERMINTE                  2.545e-01  2.995e-02   8.498  < 2e-16 ***
AVYRSEXP                  2.407e-01  5.714e-02   4.212 2.53e-05 ***
AVSALK                    8.041e-02  1.392e-02   5.775 7.69e-09 ***
PERSPENK                 -1.952e+00  3.168e-01  -6.162 7.19e-10 ***
PTRATIO                  -3.341e-01  6.126e-02  -5.453 4.95e-08 ***
PCTAF                    -1.690e-01  3.270e-02  -5.169 2.36e-07 ***
PERMINTE:AVYRSEXP        -1.408e-02  1.905e-03  -7.391 1.46e-13 ***
PERMINTE:AVSALK          -4.005e-03  4.740e-04  -8.450  < 2e-16 ***
AVYRSEXP:AVSALK          -3.906e-03  9.624e-04  -4.059 4.92e-05 ***
PERSPENK:PTRATIO          9.171e-02  1.451e-02   6.321 2.60e-10 ***
PERSPENK:PCTAF            4.899e-02  7.452e-03   6.574 4.89e-11 ***
PTRATIO:PCTAF             8.041e-03  1.499e-03   5.362 8.22e-08 ***
PERMINTE:AVYRSEXP:AVSALK  2.220e-04  2.989e-05   7.428 1.10e-13 ***
PERSPENK:PTRATIO:PCTAF   -2.249e-03  3.490e-04  -6.445 1.15e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 34345.4  on 302  degrees of freedom
Residual deviance:  4078.8  on 282  degrees of freedom
AIC: 6039.2

Number of Fisher Scoring iterations: 4

公式リファレンスの罠。

githubにあるロジスティック回帰のコードが入っていますが、これはおそらく間違いです。

statsmodels/glm_formula.py at master · statsmodels/statsmodels · GitHub

どういうわけか、以下のコマンドでロジスティック回帰が実行できますということになっていますが...

## Generalized Linear Models (Formula)

# This notebook illustrates how you can use R-style formulas to fit Generalized Linear Models.
# 
# To begin, we load the ``Star98`` dataset and we construct a formula and pre-process the data:

from __future__ import print_function
import statsmodels.api as sm
import statsmodels.formula.api as smf
star98 = sm.datasets.star98.load_pandas().data
formula = 'SUCCESS ~ LOWINC + PERASIAN + PERBLACK + PERHISP + PCTCHRT +            PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF'
dta = star98[['NABOVE', 'NBELOW', 'LOWINC', 'PERASIAN', 'PERBLACK', 'PERHISP',
              'PCTCHRT', 'PCTYRRND', 'PERMINTE', 'AVYRSEXP', 'AVSALK',
              'PERSPENK', 'PTRATIO', 'PCTAF']]
endog = dta['NABOVE'] / (dta['NABOVE'] + dta.pop('NBELOW'))
del dta['NABOVE']
dta['SUCCESS'] = endog


# Then, we fit the GLM model:

mod1 = smf.glm(formula=formula, data=dta, family=sm.families.Binomial()).fit()
mod1.summary()

問題があるのは、formulaの部分です。この形式では「目的変数」にSUCCESS(成功割合)のみを取ればよいということになっていますが、ロジスティック回帰の結果というのは「どれだけの回数試行を行ったか」に依存するため、成功の割合のみを回帰すれば良いというのは間違いです。実際、この結果を確認すると、以下のようになります。

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                SUCCESS   No. Observations:                  303
Model:                            GLM   Df Residuals:                      282
Model Family:                Binomial   Df Model:                           20
Link Function:                  logit   Scale:                             1.0
Method:                          IRLS   Log-Likelihood:                -189.70
Date:                Fri, 07 Nov 2014   Deviance:                       380.66
Time:                        20:11:04   Pearson chi2:                     8.48
No. Iterations:                     7                                         
============================================================================================
                               coef    std err          z      P>|z|      [95.0% Conf. Int.]
--------------------------------------------------------------------------------------------
Intercept                    0.4037     25.036      0.016      0.987       -48.665    49.472
LOWINC                      -0.0204      0.010     -1.982      0.048        -0.041    -0.000
PERASIAN                     0.0159      0.017      0.910      0.363        -0.018     0.050
PERBLACK                    -0.0198      0.020     -1.004      0.316        -0.058     0.019
PERHISP                     -0.0096      0.010     -0.951      0.341        -0.029     0.010
PCTCHRT                     -0.0022      0.022     -0.103      0.918        -0.045     0.040
PCTYRRND                    -0.0022      0.006     -0.348      0.728        -0.014     0.010
PERMINTE                     0.1068      0.787      0.136      0.892        -1.436     1.650
AVYRSEXP                    -0.0411      1.176     -0.035      0.972        -2.346     2.264
PERMINTE:AVYRSEXP           -0.0031      0.054     -0.057      0.954        -0.108     0.102
AVSALK                       0.0131      0.295      0.044      0.965        -0.566     0.592
PERMINTE:AVSALK             -0.0019      0.013     -0.145      0.885        -0.028     0.024
AVYRSEXP:AVSALK              0.0008      0.020      0.038      0.970        -0.039     0.041
PERMINTE:AVYRSEXP:AVSALK  5.978e-05      0.001      0.068      0.946        -0.002     0.002
PERSPENK                    -0.3097      4.233     -0.073      0.942        -8.606     7.987
PTRATIO                      0.0096      0.919      0.010      0.992        -1.792     1.811
PERSPENK:PTRATIO             0.0066      0.206      0.032      0.974        -0.397     0.410
PCTAF                       -0.0143      0.474     -0.030      0.976        -0.944     0.916
PERSPENK:PCTAF               0.0105      0.098      0.107      0.915        -0.182     0.203
PTRATIO:PCTAF               -0.0001      0.022     -0.005      0.996        -0.044     0.044
PERSPENK:PTRATIO:PCTAF      -0.0002      0.005     -0.051      0.959        -0.010     0.009
============================================================================================


回帰係数が異なっている分には実際、問題はなくて、問題があるのは「Yの推定値」に差があるかどうかです。それを確認するために、結果を確認します。まず、正しい結果はこのようになっています。

In [280]: yhat.head()
Out[280]: 
0    0.583312
1    0.751447
2    0.500583
3    0.685345
4    0.322510
dtype: float64

一方で間違っている結果は、このようになります。

In [279]: Series(mod1.mu).head()
Out[279]: 
0    0.578461
1    0.769802
2    0.449346
3    0.683891
4    0.247756
dtype: float64

先ほども説明した通り、公式リファレンスの罠ですね。下の結果は正しいロジスティック回帰の結果ではありません。この原因を作っているのは、先ほども言ったように成功割合のみでロジスティック回帰を実行したせいです。なぜなら、本来の成功回数Y(i)は、試行数N(i)、成功確率Pi(i)の2項分布に従います。しかしながら、Y(i)/N(i)の分布は当たり前ですが2項分布ではありません。確率変数変換をすれば自明です。なぜなら、P(i)の値域は0-1ですからね。結果的に、ロジスティック回帰は結果が違うことになってしまうというわけです。
注意)実際、成功割合のみを与えるようなロジスティック回帰が実行される場合には「全てのケースでの試行回数」が同じであるという仮定のもとで実行されていますので、結果は異なって当然です。

なので、pythonで行う際にはその点はご注意ください。

最後に

今回、pythonでglmを行う方法を考えましたが、まぁ実際Rでしようよっていうのが正直なところです。残念なことに笑

メモ@Python - sys.pathの問題「Extra/lib/python」という罠

夏休み以来、一度も投稿していませんでした。完全に、ブログ書くのも忘れて、研究にいそしんでおります...最近、Rからpythonに色々と移行しつつありまして、環境設定で、原因不明のエラーに苦しめられております。そろそろ、またいろんなlibraryを開拓したのでちょこちょこそういうのも載せていこうと思います。

python sys.path問題

ここ数日、pythonでnumpy, pandasなどをアップデートしたのに、アップデートしたものが、どういうわけだかimportすると昔のバージョンが使用される問題がありまして、かなり困っておりました。どうも、sys.pathに含まれている、以下のpathが悪さをしていることが分かりました。

sys.pathの中身

/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip
/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7
/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin
/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac
/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages
/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python
/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk
/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old
/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload
/Library/Python/2.7/site-packages

悪さしてるpath

'/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python'

なんだろこのpath...って思って、色々調べてみたのですが、こんな記事が出てきました...
Which Python · MacPython/wiki Wiki · GitHub

タイトルをみれば分かる通り...

「Why you should not use system Python for OSX

--OS XPythonシステムを使うべきではない理由

ということです。

で、その中にこんなことが書かれていました。
OSX 10.9のシステムでは、numpy, scipy, matplotlibや、ウェブ系で便利なパッケージが入っています。そして、そのpathとして、"/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python"というものを採用しているようです。

その結果、新しいパッケージをインストールすると、
"/Library/Python/2.7/site-packages"
に本来ならインストールされて、pythonのpath内にこれは含まれているので、参照されるわけですが、このpathが"/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python"のpathよりも「下位」のものとしてリストに含まれています。なので、まずは/Extras/lib/pythonを先に参照してから、/site-packagesを読むため、2つのサイトに「numpy」が含まれていれば、先に読まれる方が優先されます。

そして、pip install numpyでインストールすると「site-packages」に入ってしまうため、いくらpipで最新版をインストールしても"/Extra"に入っているnumpy1.6.1が参照されます。などということで、標準搭載のpythonは非常に使いにくいということでした。原文は以下。

The interesting entry is "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python" because this entry is specific to system Python. Call this path "Extras/lib". On a standard OSX 10.9 machine this directory contains some Python packages including the numerical packages numpy, scipy and matplotlib and other useful web packages including twisted.

If you, the user, install new packages in system Python, the packages will get written to /Library/Python/2.7/site-packages. As you can see from the path listing above, your newly installed package will therefore be lower in the Python module path than the packages in Extras/lib. That means that if you try to install a new version of one of the packages in Extras/lib into system Python, system Python will ignore the version you installed and continue to use the version in Extras/lib.

Apple presumably chose this setup because they use Python for system tasks, and they want to be able to depend on working versions of the packages installed in Extras/lib. This is another reason to prefer another Python for your daily Python work.

まとめ

ということだそうなので、標準搭載のpythonを捨てて、macの皆様は是非homebrewでpythonをぶちこむか、以下のリンクなどを使って、構築した方がいいのかもしれません。
pyenvでPython開発環境を構築 | youria blog

ちなみに、僕はhomebrewでpythonを入れるのをまだしていなくて、標準搭載のpythonの'Extra/lib'のフォルダーを削除してやりました笑。←ダメだと思いつつ笑

ということで、これぐらいにします。

相関の概念をさらにもう一歩深めるために。-ケンドールのタウを考えてみる-

昨日のエントリーで相関係数の有意性を確かめるためにどうすればいいのかということを書きました。

相関係数の有意性を確かめる方法について -相関係数について1歩踏み込む- - Data Science by R and Python

ですが、最後に書いたようにこの指標の有効性が検証できるのは「2つの変数が正規分布」に従うときに限られるということも最後に書きました。変数で正規分布に従わない変数はたくさんあります。身近な変数では「体重」は正規分布には従わないとされています。裾の重い分布、下の図のようになります。


f:id:tomoshige_n:20140927181754g:plain
URLBody Weight Distribution of Adult CT Scan Patients

この図、リンクを見ればわかる通りですが、右に引きづられています。この結果、分布は正規分布しないことがわかります。このようなデータ、変数同士の相関はどのようにして確認すればよいかという点を説明してみようと思います。

注意)以降の議論で使う「相関」というのは、これまでのような「直線的関係性」を意味していません。相関=関連性を意味するものとします。

ケンドールのタウ( \tau

ケンドールの \tauという指標(ケンドールの順位相関係数)を取り上げます。この指標は変数X, Yの単調性の指標です。単調性というのは「直線的」な増加や減少を指していません。「Xが増加すれば、Yは増加する」「Yが減少すれば、Xも減少する」などの傾向を指します。それが「直線的な増加」や「直線的な減少」ではなくてもよいのです。指数的、対数的な関係であっても、単調でありさえすれば良いという便利な指標です。

ただし、どのような関係で変動するのかはわかりませんし、さらに「増加、減少」などが周期的に起こるような場合にも検出ができません。この点に付いては注意する必要があります。

しかしながら、この指標の優れているところは、相関係数とは異なり、2つの変数X,Yに特別な分布(正規分布など)を仮定しないという点です。このような手法を「ノンパラメトリック」と呼びます。

ケンドールの \tauの詳しい説明

まず、ケンドールの \tauの定義を書くと、 (X_1, Y_1), (X_2, Y_2)を同一の2変量分布に従う独立な組であるとします。このとき、ケンドールのタウは次のように定義されます。

 \tau = P[sgn{(X_1 - X_2)(Y_1 - Y_2)} = 1] - P[sgn{(X_1 - X_2)(Y_1 - Y_2)} = -1]

さて、sgn{x}というのはxが正であれば1を返し、負であれば-1となるような関数です。この点を踏まえると、組の間に sgn{(X_1 - X_2)(Y_1 - Y_2)} = 1が成立すれば、増加傾向があり、逆に-1が成立するときには、減少の傾向があるということがわかります。これが定義です。

しかしながら、母集団のことは私たちにはわからないので、次のようなケンドールのタウの不偏推定量を用いてこの値をサンプルから推定するということを行います。なので、利用する際にはこちらを使います。

  K = \frac{(符号が一致する数) - (符号が一致しない数)}{n(n-1)/2}


この値は実際、期待値を取ると上で定義した \tauになることが示されます(簡単なので、興味のある人はやってみてください)。

このケンドールのタウという指標は、相関係数とにたような性質を持っています。例えば、変数XとYが独立であるとき、ケンドールの \tau = 0が成立することが示されます。もちろん、逆は成立しません。相関が0 -> XとYは独立は成立しません。

ただし、ケンドールのタウが1となるような場合は相関係数の場合と状況が異なります。先ほども説明した通り、ケンドールのタウは「単調性」に関する指標なので、XとY=X**2のタウの値は1になります。一方で、相関係数は1にはなりません。

また、データのノイズが大きいと、ケンドールのタウは検出力が大幅に低下します。この点は注意が必要になります。

実際に実行する方法

これを実際に使うためのコマンドが以下になります。データについては、ニューヨークの気象データがRに初期装備されていたので、これを使います。まずは、データを図示するところから始めます。

head(airquality)
    Ozone Solar.R Wind Temp Month Day
1      41     190  7.4   67     5   1
2      36     118  8.0   72     5   2
3      12     149 12.6   74     5   3
4      18     313 11.5   62     5   4
5      NA      NA 14.3   56     5   5

#必要な変数は最初の4つ
df = airquality[,1:4]

#散布図とヒストグラムを書く
pairs3 <- function (data) {
  oldpar <- par(no.readonly = TRUE) # 現在のグラフィックスパラメータ退避
  on.exit(par(oldpar)) # (関数がエラー中断しても)パラメータ復帰
  panel.hist <- function(x, ...) # ヒストグラムを描くための関数を定義
  { usr <- par("usr") # 現在のユーザー領域座標情報を得る
    on.exit(par(usr)) # 関数終了時に usr パラメータ復帰
    par(usr = c(usr[1:2], 0, 1.5) ) # ユーザー領域座標情報を変更
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)
  }
  pairs(data, panel=panel.smooth,
        cex = 1.0, pch = 3, bg="light blue",
        diag.panel=panel.hist, cex.labels = 2, font.labels=2)
}

#散布図 + HISTOGRAM
pairs3(df)

f:id:tomoshige_n:20140927195254p:plain


この図からわかるように、Ozoneと、Solarは正規分布に近くなさそうです。この2つを前回のエントリーの関数でqqplot, histogram重ね合わせた視覚的確認、さらにコルモゴロフ・スミルノフ検定を行ってみます。すると、以下の図と結果が得られます(結果はOzoneのみ確認します)。

f:id:tomoshige_n:20140927195457p:plain

> ks.test(na.omit(df[,1]),pnorm,mean=mean(na.omit(df[,1])),sd=sd(na.omit(df[,1])))

	One-sample Kolmogorov-Smirnov test

data:  na.omit(df[, 1])
D = 0.148, p-value = 0.01243
alternative hypothesis: two-sided

 警告メッセージ: 
In ks.test(na.omit(df[, 1]), pnorm, mean = mean(na.omit(df[, 1])),  :
   コルモゴロフ・スミノフ検定において、タイは現れるべきではありません 

ここで、警告メッセージが出されます。が、これは詳しくやり始めると大変なので無視しておきます。ks.testは値の中に「タイ(同じ数値)」があっては本当はいけないんです。

ただ、結果からもわかるように、5%有意正規分布であるという帰無仮説を棄却することができ、かつ視覚的確認からも、qqnormの裾の方がかなり外れていることがわかります。よって、正規分布と仮定するのは適切ではないという結論を得ることができました。

変数が正規分布に従わないときの検定

このような場合、例えばOzoneとTempの関係性の数値を確認したいと思って「相関係数」を計算しても、この同時分布は前回も書きましたが「2変量正規分布」ではありませんので、検定をおこなうことはできません。なので、前回の手法で有意かどうかを確認することはできないということになります。

そこで、ケンドールのタウの登場です。ケンドールのタウを実際に計算すると以下のような値が得られます。

> cor(df$Ozone,df$Temp,use="pairwise",method="kendall")
[1] 0.5862988

そして、この統計量が有意であるかどうかを調べるためには、前回同様検定を行う必要があります。このとき「検定」のためには以下の統計量を使います(この統計量は、漸近的に標準正規分布に従います。またKは上で求めたケンドールのタウです)。

 Z_{k} = \frac{K}{\sqrt{2(2n+1)/9n(n-1)}} \sim Normal(0,1)

この値を、実際に計算して検定するための関数が次の関数です。

> cor.test(df$Ozone,df$Temp,use="pairwise",method="kendall")

	Kendall's rank correlation tau

data:  df$Ozone and df$Temp
z = 9.1599, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.5862988 

以上の結果より、「z = 9.1599, p-value < 2.2e-16」であるから、帰無仮説「ケンドールのタウは0」という仮説が成立する可能性は極めて低く棄却することができ、対立仮説「ケンドールのタウは0ではない」を採用することができます。よって、オゾン(Ozone)レベルと、気温(Temp)の間には「関連性」があることが認められました。

まとめ

さて、今回はケンドールのタウを説明してきました。おそらく、これまでの「相関係数」の考え方と随分離れており、取っ付きにくい方もいるかもしれません。ただ、正規分布しないデータに対して、正規分布だと決めつけて分析を行うのは、やはり危険ですし、少し難しくても別の手法を引き出してくる方が良いと思います。特に、ケンドールのタウはノンパラメトリックですから、ある意味どんな分布に従うデータに対しても使うことができます。

ただ、直感的に理解しづらい部分もありましたので、次回は「相関係数」とよく似ているけど、「ノンパラメトリック」な考え方である「スピアマンの順位相関係数」について書こうと思います。読んでくださってありがとうございました。

(告知)イベントを開催します。

統計の初心者の方(少しやったことある人も対象)向けに、データフェストというイベントを開催します。2日間でRの操作+実際の解析を行ってみようという試みです。興味のある方は、以下のリンクから詳細をご確認ください。(サイトが90年代後半感があるのは使用です笑)
http://datafest.jp/

p.s.

今回の、ブログごちゃごちゃしてる...
ブログで複雑な話するのはやっぱり無理があるのかな...

相関係数の有意性を確かめる方法について -相関係数について1歩踏み込む-

相関係数とは?

相関係数については、8月のエントリーを参照してもらうと良いかと思います。

「相関係数」ってなんですか? -意味と利点と欠点をわかりやすく- - Data Science by R and Python

簡単に説明しておくと、相関というのは「2つの変数の間に存在する、直線的関係」を指しています。相関係数というのは、その直線度合いを0~1の数字で表したものと言えます。ただし、値には大小関係しかありません(相関係数0.8は、相関係数0.4のときよりも2倍相関が強いという意味ではないので注意してください)。

さて、データの分析・解析に関わらず、いくつかの変数があるデータでは、相関係数を計算して、変数の間に関係があるかをチェックするのが一般的です。前のエントリーにもあるようにこれは「直線的な関係」しか調べられないですし、色々と制約もあり、いちいち確認するのは面倒ですが、それでも使い勝手がいいのでみんな使っています。

さて、今日のモチベーションは相関係数を説明することではありません。「相関係数」を求めたときに、果たして求めた相関係数が「有意」であるか、「有意」ではないかに興味があるとします。つまり、求められた「相関係数」は統計的に意味がある数字であるのか、それとも「偶然」や「見かけ上」相関があるようにみえているだけで、本当は相関がない可能性もあるのか。この点を調べていきたいとします。つまり、「求めた相関係数が0である可能性があるか、ないか」を確認する方法についてです。

特に、相関係数を計算したあとに、その相関係数が「意味があるのか、ないのか」をチェックせずに進めてしまったりすると、問題になるかもしれませんし、「0.31の相関がありました!」のような、統計的に意味があることが保証されない報告をしてしまうことになるかもしれません。

なので、今回は求めた相関係数が有意かどうかを確かめる方法について考えていきます。


相関係数が統計的に有意(意味があるか)を検定する方法の前準備

まず、相関係数(ピアソンの積率相関係数)が有意かどうかを調べるためには、相関を計算する変数の同時分布が「2変量正規分布」に従うのかを調べておく必要があります(この仮定が成立しない場合は相関係数が有意かどうかを統計的に確かめることは難しいです)。

そのためには、次の事実を利用します。「2つの確率変数が、それぞれ正規分布(平均・分散は異なっても良い)に従うとき、その同時分布は正規分布に従う」。つまり、相関を確認したい2つの変数それぞれが正規分布に従うことを確認すればよいということになります。

直感的な方法(図で診断する!)

まず1つ目の簡単な方法は、「ヒストグラム」です。ヒストグラムの外形を見たときに、それが正規分布の外形と一致しているのかを確認します。そのためには、得られたデータから平均と分散を計算して、それらをパラメータにする理論的な正規分布ヒストグラムの上から書くのが良いでしょう。この方法はとても直感的で、おそらく確認の方法として粗くはありますが、強力です。

それから、もう1つはqqnormを用いる方法です。qqplotは、正規分布の累積分布関数を用います。細かい説明は、英語のwikipediaを確認してください。(Q–Q plot - Wikipedia, the free encyclopedia)。簡単に説明をしておくと、累積分布の25%点と75%点を直線で結び、その直線上にデータが乗るかを確認することで、正規分布に従っているかを確かめる方法です。

f:id:tomoshige_n:20140926220232p:plain

ここまで紹介してきた2つの方法を、手軽にできるようにするために、以下にRで作成した関数を添付しておきます。

x_norm = rnorm(1000,10,5)
x_unif = runif(1000)

range = function(val){
	if((class(val) == 'numeric')| (class(val)=='vector')){
		diff = (max(val) - min(val))/10
		return(seq(min(val)-diff,max(val)+diff,(max(val)-min(val))/1000))
	}else{
		print("only one variable available")
	}
}

norm_diag = function(x){
	par(mfcol=c(1,2))
	diff = (max(x) - min(x))/10
	xrange = c(min(x)-diff,max(x)+diff)
	hist(x,breaks=floor(length(x)/20),freq=F,xlim=xrange)
	par(new=T)
	points(range(x),dnorm(range(x),mean=mean(x), sd = sd(x)),type="l",xlim=xrange)
	segments(x0=mean(x),y0=0,x1=mean(x),y1=dnorm(mean(x),mean(x),sd(x)),lty=2,col="red",lwd=1.5)	
	qqnorm(x,cex=0.5)
	qqline(x,col="red",lty=2)
}

#変数が正規分布の場合
norm_diag(x_norm)

#結果はOKだと視覚的に確認できるはず

#変数が一様分布の場合
norm_diag(x_unif)

#結果はNOだと視覚的に確認できるはず

理詰めしたいときの検定

ですが、このままだと視覚的確認なので不安だという人のために、もう1つ統計的仮説検定を用いて正規分布に従うかどうかを検定する方法を紹介します。この検定は「コルモゴロフ・スミルノフ検定」と言います。仮説の設定をまずは確認しましょう。

ここで、注意が必要なのは「帰無仮説が棄却されれば、正規分布を仮定するのは適切ではない」という結論になることです。なので、正規分布を仮定しても、大きな問題にならないケースであることを示すには、P値が大きいことが必要です。

この検定の中身への言及はさけますが、簡単に言えば、経験分布関数と、理論的な分布関数の差の最大値を検定統計量として用いています。詳しくは、医学統計学の事典(丹後俊郎, 小西貞則 2010)を参考にしてください。この他、wikipediaなどにも記述があります。この検定をRで行うためには、以下のコマンドを実行すればOKです。簡単ですね。

#平均と分散をきちんと指定する必要があります
#正規分布から発生させたデータの場合
> ks.test(x_norm,"pnorm",mean=mean(x_norm),sd=sd(x_norm))

#	One-sample Kolmogorov-Smirnov test

#data:  x_norm
#D = 0.02, <b><span style="color: #ff2600">p-value = 0.8185</span></b>
#alternative hypothesis: two-sided

#一様分布から発生させたデータの場合
> ks.test(x_unif,"pnorm",mean=mean(x_unif),sd=sd(x_unif))

	One-sample Kolmogorov-Smirnov test

data:  x_unif
D = 0.0599, p-value = 0.001522
alternative hypothesis: two-sided


結果を見れば、(上)正規分布から発生させた乱数では帰無仮説は棄却できませんが、(下)一様分布から発生したものは帰無仮説が棄却されてしまいます。よって、(上)は正規分布を仮定しても問題にはならなそうだということが推察され、(下)は正規分布を仮定するのは適切ではないという結論を得ます。

ちなみに、この他にも、検定手法はいろいろ提案されています。医学統計学の事典には、アンダーソン・ダーリング検定などが紹介されています。検定には、それぞれ特徴がありますので、他の検定手法では有意ではないなどの結果が出ることもあります。

※注意)ここで、注意ですが、「帰無仮説が棄却できなかった = 正規分布に従っている」は間違いです。帰無仮説が棄却できないというのは、帰無仮説が有意水準5%で棄却できないだけで、「変数が正規分布に従っている」と安易に結論づけてはいけません。

相関係数の検定

さて、ここまでの結果、相関を見たい変数が正規分布に従うと仮定しても問題ないとわかったとしましょう。すると、始めて相関係数が有意かどうかを検定することができます。

まずは、仮説の設定から始めます。
「帰無仮説:相関係数が0である」
「対立仮説:相関係数は0ではない」

ここで、「相関係数が0である」というのは、確率変数には正規分布を仮定していますから、独立性の検定と同等です。これを検定するためには、2つの方針があります。1つ目は尤度比検定です。もう1つが、帰無仮説の元での標本相関係数の分布を求めるというやり方です。今回は後者を考えます。

さて、帰無仮説の元での標本相関係数の分布は、t分布に従います。この証明は、非常に難解です。興味のある人は以下の資料をご覧下さい。おそらく、ほとんどの人がやらないでもいいんじゃないでしょうか...
Pearson product-moment correlation coefficient - Wikipedia, the free encyclopedia
Correlation Coefficient--Bivariate Normal Distribution -- from Wolfram MathWorld


こんな難しいことを知らなくても、Rを使えば簡単に検定を行うことができます。2つの正規分布に従う確率変数の間に相関があるかないかは、次のコマンドで実行できます。

#無相関
> cor.test(x_norm, rnorm(1000))

#Pearson's product-moment correlation
#
#data:  x_norm and rnorm(1000)
#t = 0.097, df = 998, p-value = 0.9227
#alternative hypothesis: true correlation is not equal to 0
<b><span style="color: #ff2600">#95 percent confidence interval:
# -0.05893338  0.06505162</span></b>
#sample estimates:
#        cor 
#0.003070921 

#相関あり
> cor.test(x_norm, x_norm + rnorm(1000,0,8))

#Pearson's product-moment correlation
#
#data:  x_norm and x_norm + rnorm(1000, 0, 8)
#t = 20.5626, df = 998, p-value < 2.2e-16
#alternative hypothesis: true correlation is not equal to 0
<span style="color: #ff2600"><b>#95 percent confidence interval:
# 0.5004476 0.5876367</b></span>
#sample estimates:
#      cor 
#0.5455164 

これを確認して、有意であれば95%信頼区間が0をまたがない出力を得ます。前者は0を跨いでしまっているので、有意ではないですね。一方で、後者は0をまたいでないので、有意に相関があることがわかります。

まとめ

今回は、相関係数が有意であるかどうかの検定について書きました。相関係数が有意であるかどうかを統計的に確認するのは、様々な仮定が必要です。そして、満たさないデータもたくさんあります。ですが、満たすデータもあるはずですし、そのような場合には、上記のような検定は効果を発揮すると思います。

でも、満たさない場合も扱いたい!そんな人もいらっしゃるかもしれません。次回は、正規分布に従わない場合に変数同士の関係性(相関ではなく)を確認する方法の1つである「スピアマンの順位相関」について、それから見かけ上の相関(疑似相関)を取り除いて純粋に2つの変数の関係性を確認するための考え方である「偏相関」について書きたいと思います。

この記事の続き

tomoshige-n.hatenablog.com


(告知)イベントを開催します。

統計の初心者の方(少しやったことある人も対象)向けに、データフェストというイベントを開催します。2日間でRの操作+実際の解析を行ってみようという試みです。興味のある方は、以下のリンクから詳細をご確認ください。(サイトが90年代後半感があるのは使用です笑)
http://datafest.jp/


参考文献

PypeRでpythonからRを使う -- データ解析をpythonで!

PythonからRを使うための方法.

pythonからRを使うためには、pythonにpypeRというパッケージ?をインストールする必要があります。こちらです:PypeR

$pip install PypeR

で導入して使ってみてください。

実際にPypeRを使ってみます。

PypeRを使ってみた結果は以下から参照できます。ipython notebookが便利なので、これからはここにプログラムを書くのをやめて、こんな感じで公開するのもいいかな。。。と思ったりしています。
いや、本当に便利なんですnotebookが。そういえば、Rpubsとか、どうなんでしょう?使いやすいんでしょうか?PypeRを使ってみた結果を解説付きでこちらで公開しました。

http://nbviewer.ipython.org/github/tomoshige/python_blog/blob/master/pypeR/pypeR_sample.ipynb

f:id:tomoshige_n:20140909024509p:plain

まとめ

うーん、一応pythonでRは使えるけど、やっぱりまだまだ不便な部分もあるなぁと。。。慣れれば楽になるんですかね?特にいちいちRに渡して、結果を戻したりするのが一手間かかかるんですよね。やはり、pythonはデータ取得・データ加工に特化させて、Rで解析・解釈するのがいいんでしょうかね。でも、全部pythonで完結できるようにした方が後々いいんじゃないかなとか思ったりもしますね。

補足

さて、そろそろ、本格的にscikit-learn使い始めようとおもっているのと、最近クローリングとスクレイピングの勉強をしてるのですが、なかなか進みません(泣。。。おなかすいた...

pythonでggplotを使おう! -pythonによるデータ分析(2)-

ggplot2をpythonで使う。

Rを使っている人にとっては、ggplot2はおなじみな感じがしますが、pythonでもこれ使えるといいなということで、pythonで使える方法を探しておりましたら、"ggplot"はpythonでもあるんですね。

ggplot 0.6.5 : Python Package Index

ということで、これを積極的に使っていくためにセットアップします。あんまり日本語の記事もないですし。ちなみに、Rではggplot2でしたが、pythonではggplotです。

依存パッケージなど

ggplotは以下のパッケージに依存しています、pythonでのデータ解析で、使うモノばかりなのでインストールしている人も多いかと思いますが、ない人はインストールをしてください。

  • matplotlib
  • pandas
  • numpy
  • scipy
  • statsmodels
  • patsy

インストール

インストールは、いつも通りpipを使ってできます。

$ pip install ggplot

基本操作

基本操作はRとほとんど変わらないので、気にすることはないと思います。これで、pythonでggplotができるようになりました。

pythonでデータ分析(1)

f:id:tomoshige_n:20140907233155j:plain

R -> python

さて、pythonでRと同程度以上のデータ解析をできるようになろうと決意して、はや1週間が経ちました。まだ、1週間かよ!という突っ込みはやめてください笑。ここ数日は、全力でpythonによるデータ分析入門を読み込んでいます。というか、写経して、基本的なデータ分析の操作を覚えています。
O'Reilly Japan - Pythonによるデータ分析入門

本の内容に沿って、やってみたこと。

まず、Movielensの映画評価データを使います。Movielensは映画評価システムの1つで、ミネソタ大学のコンピュータ理工学研究科のグループレンズが運営しています。映画評価データは1990年代後半から、2000年代前半に書けて、Movielensユーザーを対象に収集されてきました。各データは、その映画への評価値、映画の属性(ジャンル・公開年度)、評価者の属性(年齢・郵便番号・性別・職業)から構成されています。

このデータは、データ数は100万件あります。4000本の映画に対する、6000人の評価をまとめたものになっています(つまり、それぞれが全ての映画にレビューを書くわけではないということです)。データは3つの表(評価・ユーザー情報・映画情報)から構成されています。

These files contain 1,000,209 anonymous ratings of approximately 3,900 movies
made by 6,040 MovieLens users who joined MovieLens in 2000.


MovieLens | GroupLens

それぞれのデータを確認すると、まずは「ユーザー」情報は。

In [52]: users
Out[52]: 
      user_id gender  age  occupation    zip
0           1      F    1          10  48067
1           2      M   56          16  70072
2           3      M   25          15  55117
3           4      M   45           7  02460
4           5      M   25          20  55455
5           6      F   50           9  55117
6           7      M   35           1  06810
7           8      M   25          12  11413
8           9      M   25          17  61614
9          10      F   35           1  95370

となっていて、評価結果はこのようになっています。これらのデータはuseridでひもづいています。

In [53]: ratings[:10]
Out[53]: 
   user_id  movie_id  rating  timestamp
0        1      1193       5  978300760
1        1       661       3  978302109
2        1       914       3  978301968
3        1      3408       4  978300275
4        1      2355       5  978824291
5        1      1197       3  978302268
6        1      1287       5  978302039
7        1      2804       5  978300719
8        1       594       4  978302268
9        1       919       4  978301368

そして、映画情報はこんな感じ。これは評価データとmovie_idでひもづいているデータです。

In [54]: movies[:10]
Out[54]: 
   movie_id                               title                        genres
0         1                    Toy Story (1995)   Animation|Children's|Comedy
1         2                      Jumanji (1995)  Adventure|Children's|Fantasy
2         3             Grumpier Old Men (1995)                Comedy|Romance
3         4            Waiting to Exhale (1995)                  Comedy|Drama
4         5  Father of the Bride Part II (1995)                        Comedy
5         6                         Heat (1995)         Action|Crime|Thriller
6         7                      Sabrina (1995)                Comedy|Romance
7         8                 Tom and Huck (1995)          Adventure|Children's
8         9                 Sudden Death (1995)                        Action
9        10                    GoldenEye (1995)     Action|Adventure|Thriller

このデータをまずは、pandasにあるmerge関数で統合する必要があります。これはRでもmerge関数があって、どこを引数にするのかを指定すれば、マージしてくれます。このあと、性別毎の平均などを導出して、性別で評価が変わる映画は何かなどを出していきます。コードはこんな感じ。

#data merge
data = pd.merge(pd.merge(ratings,users),movies)

#1行目を参照する
data.ix[0]

#男女/タイトルで結果を出力する
mean_ratings = data.pivot_table('rating',rows='title',cols='gender',aggfunc='mean')

#分析対象を250件以上のレビューのある映画に限る
#titleでグループ分けして、sizeを調べる
rating_by_title = data.groupby('title').size()

#250件以上に絞り込む
active_title = rating_by_title.index[rating_by_title>=250]

#mean_ratingsの中でactive titleのみを抜き出す
mean_ratings = mean_ratings.ix[active_title]

#女性の中で評価が高かった映画の抽出
top_female_ratings = mean_ratings.sort_index(by="F",ascending=False)



#評価の分かれた映画の抽出
mean_ratings['diff'] = mean_ratings['M'] - mean_ratings['F']

#sort by diff
sort_by_diff = mean_ratings.sort_index(by='diff')

1週間経ってみて。

Rで関数を書くよりもpythonで書いてる方が、まずなんというか「楽しい」ということがわかりました。もちろん、不便だなと感じることもあります。記述量が多いです。ただ、Rで処理をすると遅いことも、pythonならやっぱり体感値ですが、速いなぁと感じます。
あとは、プログラム書ける周りの人と話してても、Rは取っ付きにくいし、学習コストが大きいなんて言われたりもします(Rになれると、その感覚はわからないんですけど、そういう感覚もあるようですから)。なので、僕がpythonでデータ処理ができるようになれば、そういう人達にも伝えるお手伝いができるようになるかなとも思うので・・・引き続き、がんばります。