R vs Python:データ解析を比較

主観的な観点からPythonとRの比較した記事は山ほどあります。それらに私たちの意見を追加する形でこの記事を書きますが、今回はこの2つの言語をより客観的な目線で見ていきたいと思います。PythonとRを比較をしていき、同じ結果を引き出すためにはそれぞれどんなコードが必要なのかを提示していきます。こうすることで、推測ではなく、それぞれの言語の強みと弱みの両者をしっかりと理解できます。Dataquestでは、PythonとRの両方の言語のレッスンを行っていますが、データサイエンスのツールキットの中では両者ともそれぞれに適所があります。

この記事では、NBA選手の2013/2014年シーズンの活躍を分析したデータセットを解析していきます。ファイルはここからダウンロードしてください。解析はまずPythonとRのコードを示してから、その後に2つの異なるアプローチを解説し議論していきます。つまらない説明はこれくらいにして、データ解析の比較を早速スタートさせていきましょう!

csvファイルの読み込み

R

nba <- read.csv("nba_2013.csv")

Python

import pandas
nba = pandas.read_csv("nba_2013.csv")

このコードを実行すると、NBA選手の2013/2014年シーズンのデータの入ったcsvファイルnba_2013.csvがそれぞれの言語で変数nbaに読み込まれます。1つだけ違うのは、Pythonではデータフレームにアクセスできるようにするためpandasライブラリをインポートしなければいけないことです。DataFrameはRとPythonの両者に使える2次元の配列(マトリックス)で、それぞれのカラムは別のデータ型にすることができます。最終段階では、csvファイルは両者の言語によってデータフレームに読み込まれます。

選手の数を割り出す

R

dim(nba)

[1] 481  31

Python

ba.shape

(481, 31)

このコードで選手の数とそれぞれのカラムの数が出力されます。上記によると481個の行、つまり選手の数と、選手のデータを含んだ31個のカラムがあるということになります。

データの最初の行を見る

R

head(nba, 1)

      player pos age bref_team_id
1 Quincy Acy  SF  23          TOT
[output truncated]

Python

nba.head(1)

       player pos  age bref_team_id
0  Quincy Acy  SF   23          TOT
[output truncated]

この2つはほぼ同じですね。2つとも最初の行を出力していて、構文がよく似ています。Pythonはよりオブジェクト指向で、headはデータフレームオブジェクトのメソッドです。Rでは単独のhead関数になります。これは、RとPythonで解析をするときに初めに直面する典型的なテーマです。Pythonはよりオブジェクト指向で、Rはより関数的なのです。

各統計量から平均を割り出す

それぞれの統計量から平均値を割り出してみましょう。ご覧のように、カラムにはfg(フィールドゴール成功数)やast(アシスト数)などの名前がついています。これらは選手たちのシーズンの統計量です。統計量のより詳細なデータはこちらをご覧ください。

R

sapply(nba, mean, na.rm=TRUE)

player NA
pos NA
age 26.5093555093555
bref_team_id NA
[output truncated]

Python

nba.mean()

age             26.509356
g               53.253638
gs              25.571726
[output truncated]

今回はアプローチに大きな違いがいくつか出ました。両者ともデータフレームのカラムに関数を適用させました。Pythonでは、データフレーム上のmeanメソッドがデフォルトで各カラムの平均を割り出してくれます。

Rでは、文字列の値の平均はNA(値なし)という結果になっています。しかし、平均を求める時にはNAの値は無視しなければいけません(na.rm=TRUEmean関数に渡す必要があります)。そうしないと、x3p.のようなカラムの平均は結局、NAになってしまいます。このカラムはスリーポイントシュート率を表しています。スリーポイントを一度も決めてない選手の欄には数字が入っていません。Rでmean関数を試してみると、NAが返されます。なので、平均を取るときはNA値を無視するna.rm=TRUEを記述しましょう。Pythonでのmean()メソッドは、デフォルトでNA値を無視するようになっています。

対散布図を作成する

データセットを調べる一般的な方法は、それぞれのカラムが他のカラムとどのような相互関係にあるかを調べることです。
astfgtrbのカラムを比べてみましょう。

R

library(GGally)
ggpairs(nba[,c("ast", "fg", "trb")])

Python

import seaborn as sns
import matplotlib.pyplot as plt
sns.pairplot(nba[["ast", "fg", "trb"]])
plt.show()

比べるとかなり似た図になっています。これは、Rのデータサイエンスのエコシステムが多くのより小さなパッケージ(GGallyggplot2のヘルパーパッケージで、Rでプロットする際に最も使われるパッケージ)やより視覚的なパッケージを持っていることを一般的に示しています。Pythonでは、matplotlibが一番重要なプロットのパッケージで、seabornはmatplotlibよりも広く使われるレイヤです。Pythonの可視化では、何かやるときの主な方法は通常1つですが、Rでは異なる方法をサポートするたくさんのパッケージがあります(例えば、対になるプロットを作るのに少なくとも6つくらいのパッケージがあります)。

選手のクラスタを作成する

このようなデータを探すのにいい方法は、クラスタプロットを生成することです。これで、どの選手が一番類似しているかが分かります。

R

library(cluster)
set.seed(1)
isGoodCol <- function(col){
   sum(is.na(col)) == 0 && is.numeric(col) 
}
goodCols <- sapply(nba, isGoodCol)
clusters <- kmeans(nba[,goodCols], centers=5)
labels <- clusters$cluster

Python

from sklearn.cluster import KMeans
kmeans_model = KMeans(n_clusters=5, random_state=1)
good_columns = nba._get_numeric_data().dropna(axis=1)
kmeans_model.fit(good_columns)
labels = kmeans_model.labels_

正しくクラスタリングをさせるには、数字ではないカラム、もしくは値のない(NANanなど)カラムは排除します。Rでは、それぞれのカラムに関数を適用し、値がない、数字ではないものがあれば排除するようにします。そして、k-means法を実行するのにclusterパッケージを使い、データ上で5個のクラスタを割り出します。
set.seedを使ってランダムなシード値を設定し、結果を求められるようにします。

Pythonでは、メインのPython機械学習パッケージのscikit-learnを使って、モデルをクラスタリングするk-means法に適合させ、クラスタのラベルを取得します。Rで使うデータを準備するためには、非常に類似したメソッドを実行します。数字ではないカラムと値のないカラムを排除するためには、get_numeric_datadropnaのメソッドを使います。

クラスタごとに選手をプロットする

では、選手をクラスタごとにプロットしてパターンを見つけましょう。これをやるには、まずPCA(主成分分析)をして2次元データを作成し、プロットします。そして、クラスタの関係に従ってそれぞれの印に影をつけましょう。

R

nba2d <- prcomp(nba[,goodCols], center=TRUE)
twoColumns <- nba2d$x[,1:2]
clusplot(twoColumns, labels)


訳:これらの2つの成分はばらつきを100%説明している

Python

from sklearn.decomposition import PCA
pca_2 = PCA(2)
plot_columns = pca_2.fit_transform(good_columns)
plt.scatter(x=plot_columns[:,0], y=plot_columns[:,1], c=labels)
plt.show()

今回のデータから散布図を作成して、クラスタに従ってデータのアイコンに影をつけるか、もしくは変更します。Rでは、クラスタのライブラリの一部であるclusplot関数が使われています。Rのビルトインのpccomp関数を用いて、PCAを実行します。

Pythonではscikit-learnライブラリのPCAクラスを使っていて、matplotlibを用いて図を作成しています。

訓練データセットとテストデータセットに分割する

教師ありの機械学習を行いたければ、訓練データセットとテストデータセットに分割して過剰適合しないようにするのがいいでしょう。

R

trainRowCount <- floor(0.8 * nrow(nba))
set.seed(1)
trainIndex <- sample(1:nrow(nba), trainRowCount)
train <- nba[trainIndex,]
test <- nba[-trainIndex,]

Python

train = nba.sample(frac=0.8, random_state=1)
test = nba.loc[~nba.index.isin(train.index)]

Rには、floorsampleset.seedなど、データ解析にフォーカスしたビルトイン関数がより多くあるのがお分かりでしょうか。一方で、Pythonではパッケージを用いています(math.floorrandom.samplerandom.seed)。Pythonでは、最新バージョンのpandasにsampleメソッドがあります。これは、ソースのデータフレームからある割合の数の行をサンプルとして抽出して返します。これでコードはより簡単になります。Rでは、サンプリングを単純にするパッケージがありますが、ビルトインのsample関数を使ったほうが簡単です。どちらにしろ、ランダムなシード値を設定して、結果を再現可能にします。

単回帰

次は、選手ごとのフィールドゴールからアシストの数を予測してみます。

R

fit <- lm(ast ~ fg, data=train)
predictions <- predict(fit, test)

Python

from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(train[["fg"]], train["ast"])
predictions = lr.predict(test[["fg"]])

Scikit-learnは、予測を適合させて生成することのできる単回帰モデルを持っています。Rではビルトインのlm関数とpredict関数を用います。predictは、渡された適合モデルの種類によって違う振る舞いをします。つまり、適合したいさまざまなモデルに使えます。

モデルの統計量の概要を計算する

R

summary(fit)

Call:
lm(formula = ast ~ fg, data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-228.26  -35.38  -11.45   11.99  559.61
[output truncated]

Python

import statsmodels.formula.api as sm
model = sm.ols(formula='ast ~ fga', data=train)
fitted = model.fit()
fitted.summary()

OLS Regression Results                            
============================
Dep. Variable:                    ast   
R-squared:                       0.568
Model:                            OLS   
Adj. R-squared:                  0.567
[output truncated]

R Square(決定係数)のように、適合に関して要約統計量を出したいのであれば、RよりPythonのほうがやることが少し多くなります。Rではビルトインのsummary関数を使うことで、モデル上の情報を引き出せます。Pythonでは、Pythonで使われる多くの統計メソッドを実行可能にしてくれるstatsmodelsパッケージを使う必要があります。どちらもほぼ同じ結果が得られますが、大抵の場合Pythonで統計分析をするほうが少しだけ複雑になり、Rに存在する統計メソッドのいくつかはPythonには存在しません。

ランダムフォレストモデルを適合させる

線形回帰は変数が1つの場合に効果的ですが、ここで私たちはデータに非線形があるか探ってみます。ということで、ランダムフォレストモデルを適合させます。

R

library(randomForest)
predictorColumns <- c("age", "mp", "fg", "trb", "stl", "blk")
rf <- randomForest(train[predictorColumns], train$ast, ntree=100)
predictions <- predict(rf, test[predictorColumns])

Python

from sklearn.ensemble import RandomForestRegressor
predictor_columns = ["age", "mp", "fg", "trb", "stl", "blk"]
rf = RandomForestRegressor(n_estimators=100, min_samples_leaf=3)
rf.fit(train[predictor_columns], train["ast"])
predictions = rf.predict(test[predictor_columns])

2つの大きな違いは、アルゴリズムを使うためにRではランダムフォレストのライブラリを使う必要があるということです。一方で、Pythonではscikit-learnにビルトインされています。scikit-learnはPythonのさまざまな種類の機械学習アルゴリズムと一緒に使える統合インタフェースを持っていて、通常はPythonの各アルゴリズムの主要な実装が1つあるだけです。Rでは、個々のアルゴリズムを含んだ小さなパッケージが多くあり、アクセスするのに一貫性を持たないこともあります。これによってアルゴリズムに多様性が出てきますが(多くはいくつかの実装があり、研究所から出たばかりです)、使い勝手は良くありません。

誤差を計算する

では、2つのモデルが適合したところで、次は誤差を計算しましょう。MSE(平均平方誤差)を使います。

R

mean((test["ast"] - predictions)^2)

4573.86778567462

Python

from sklearn.metrics import mean_squared_error
mean_squared_error(test["ast"], predictions)

4166.9202475632374

Pythonでは、scikit-learnライブラリには、使用できるさまざま誤差のメトリックスがあります。RではMSEを計算するより、小さないくつかのライブラリを持つ傾向があります。しかし、どちらの言語にしても手動で行うほうが断然簡単です。誤差には小さな相違がありますが、これは明らかにパラメータチューニングが原因です。でもこれは大きな問題ではありません。

Webページをダウンロードする

2013/2014シーズンのNBA選手のデータが手元にありますが、次はこのデータを補う追加のデータを取得してみましょう。時間を節約するために、ここにあるNBA決勝戦のボックススコアを1つ使うことにしましょう。

R

library(RCurl)
url <- "http://www.basketball-reference.com/boxscores/201506140GSW.html"
data <- readLines(url)

Python

import requests
url = "http://www.basketball-reference.com/boxscores/201506140GSW.html"
data = requests.get(url).content

Pythonでは、全てのリクエストの種類に対して一貫性のあるAPIを使用することでRequestsパッケージがWebページをダウンロードするのを簡単にしてくれます。Rでも同じようにRCurlがリクエストの作成を容易にしてくれます。どちらもWebページを文字列のデータ型にダウンロードします。注釈:このステップは、次のRのステップに必要がありませんが、比較のために提示しました。

選手のボックススコアを抽出する

Webページがダウンロードできたので、次に選手のスコアを抽出するためにこれをパースする必要があります。

R

library(rvest)
page <- read_html(url)
table <- html_nodes(page, ".stats_table")[3]
rows <- html_nodes(table, "tr")
cells <- html_nodes(rows, "td a")
teams <- html_text(cells)

extractRow <- function(rows, i){
    if(i == 1){
        return
    }
    row <- rows[i]
    tag <- "td"
    if(i == 2){
        tag <- "th"
    }
    items <- html_nodes(row, tag)
    html_text(items)
}

scrapeData <- function(team){
    teamData <- html_nodes(page, paste("#",team,"_basic", sep=""))
    rows <- html_nodes(teamData, "tr")
    lapply(seq_along(rows), extractRow, rows=rows) 
}

data <- lapply(teams, scrapeData)

Python

from bs4 import BeautifulSoup
import re
soup = BeautifulSoup(data, 'html.parser')
box_scores = []
for tag in soup.find_all(id=re.compile("[A-Z]{3,}_basic")):
    rows = []
    for i, row in enumerate(tag.find_all("tr")):
        if i == 0:
            continue
        elif i == 1:
            tag = "th"
        else:
            tag = "td"
        row_data = [item.get_text() for item in row.find_all(tag)]
        rows.append(row_data)
    box_scores.append(rows)

このコードは、2つのリストを含んだリストを作成します。1つ目のリストは、CLEのボックススコア、2つ目のリストは、GSWのボックススコアとなります。どちらのリストにもヘッダー、選手名、そして彼らのゲームでの統計量が含まれています。ここでは、訓練データに変換することはしませんが、簡単にnbaデータフレームに落とし込める形式に変換することができます。

RのコードはPythonのコードに比べて、より複雑になっています。これは、項目を選択する際の正規表現の使用において、便利な方法がないからです。よって、HTMLからチーム名を取得するために追加のパースが必要となります。またRでは、ベクトルに従って関数を適用するのに有利なforループを使っていません。代わりにlapplyを使っていますが、その行がヘッダーなのかどうかに応じて各行を取り扱う必要があるので、必要な項目のインデックスを渡し、rowsリスト全体を関数に入れ込みました。

また、最新かつ広く使われているRのWebスクレイピングパッケージ、rvestを使っています。このコードで、必要なデータを抽出することが可能です。

Pythonでは、最も良く使われているWebスクレイピングパッケージ、BeautifulSoupを使いました。これにより、容易にタグを通してループしたり、リストのリストを構成したりすることが可能になります。

結論

これまでRとPythonを比較して、それぞれのデータセットをどう解析するかを見てきました。ここでは深く掘り下げなかった多くのタスクがまだあります。例えば、解析結果を主張したり、他の人と共有したり、テストや本番準備を行ったり、また、より視覚化したりするなどです。これらについては追って対応していきたいと思います。そうすることで、より確実な結論を打ち出すことができるでしょう。ひとまず、ここで話した内容いついて結論を記述しておきます。

Rはより関数的で、Pythonはよりオブジェクト指向である

lmpredict、その他の関数を見てきましたが、Rではほとんどの処理で関数を使用しているのに対し、PythonではLinearRegressionクラスやデータフレームのsampleメソッドを使用しています。

Rはより多くのデータ分析がビルトインされており、Pythonはパッケージに依存している

統計データの概要を見てみると、Rではsummaryというビルトイン関数を使用していますが、Pythonでは、statsmodelsというパッケージをインポートしています。Rの場合、データフレームはビルトインされていますが、Pythonの場合は、pandasパッケージからインポートしなければなりません。

Pythonには、データ解析のタスクにおいて”メイン”パッケージがありますが、Rには、小さなパッケージのより大規模なエコシステムがあります。

Pythonでは、線形回帰やランダムフォレスト、scikit-learnパッケージを使ったより多くの対応が可能です。一貫性のあるAPIが提供されており、よく管理されています。Rには、多様なパッケージがより多くあるのですが、Pythonに比べるとフラグメンテーションが起こりやすく、一貫性も劣っています(線形回帰はビルトイン、lmrandomForestは別のパッケージである、など)。

Rは一般的に、より統計をサポートしている

Rは、統計を目的とした言語であり、実際にそれを体現しています。Pythonのstatsmodelsや他のパッケージでは、統計メソッド向けの十分な適用を提供していますが、Rのエコシステムは、これらよりもはるかに大規模です。

統計以外のタスクを行う場合は、Pythonの方が単純である

PythonのBeautifulSoupといった十分に管理されたライブラリやリクエスト、Webスクレイピングは、Rよりもはるかに単純です。今回、この記事で述べていないようなタスクにも同じことが言えます。例えば、データベースへの保存やWebサーバーのデプロイ、また複雑なワークフローの実行などが挙げられます。

両言語のデータ解析のワークフローには多くの類似性がある

R、Pythonのどちらも相互的に刺激しあっている点があります(pandasのデータフレームはRのデータフレームの影響を受けていますし、rvestパッケージはBeatufiulSoupの影響を受けています)。そして、それぞれのエコシステムもより強力になってきています。シンタックスやアプローチ方法など、多くの類似した共通タスクが両者の言語にあるというのは、注目するに値します。

最後に

Dataquestでは、元々、Pythonのレッスンだけを行っていましたが、最近はRも行っています。両者の言語は、お互いに補い合う関係だと思っています。Pythonの方が、多くの領域で強力だと感じてはいますが、Rも有効な言語です。データ調査や解析、または単独のデータ解析ツールなど、Pythonを補う言語としても使用できるかもしれません。今回の検証で証明されたように、両者の言語には、シンタックスやアプローチ方法といった共通点が多くあることが分かりました。これについては、どちらの言語においても疑う余地のない事実です。