tidymodelsによるtidyな機械学習（その1：データ分割と前処理から学習と性能評価まで）
目次
※この記事をベースにした2019年12月7日に行われたJapan.R 2019での発表資料は以下になります。 tidymodelsによるtidyな機械学習 - Speaker Deck
はじめに
本記事ではtidymodelsを用いたtidyな機械学習フローを紹介したいと思います。
tidyverseはデータハンドリングと可視化のためのメタパッケージでしたが、tidymodelsはtydyverseにフィットするやり方で統計モデリング/機械学習をするためのメタパッケージになります。
tidymodels配下のパッケージは量が多く使い所が限られているパッケージも多いため、一度に全ては紹介できません。 ですので、今回は典型的な
- 訓練データとテストデータの分割
- 特徴量エンジニアリング
- モデルの学習
- モデルの精度評価
という機械学習フローを想定し、各パッケージの機能と連携をコンパクトにまとめたいと思います。
各パッケージの用途は以下になります。
- rsample：訓練データ/テストデータの分割
- recipe：特徴量エンジニアリング
- parsnip：異なるパッケージのアルゴリズムを同じシンタックスで扱えるようにするラッパー
- yardstic：モデルの精度評価
tidyな機械学習フロー
tidyverseとtidymodelsを読み込みます。
library(tidyverse) library(tidymodels)
訓練データとテストデータの分割
まず第一に大元のデータを訓練データとテストデータに分割しましょう。
データはdiamondsを使って、ダイヤの価格を予測することにします。
set.seed(42) df_split = initial_split(diamonds, p = 0.8) df_train = training(df_split) df_test = testing(df_split)
データ分割の役割はrsampleが担っています。
initial_split()で分割するデータと、訓練データが占める割合を指定します。
training()で訓練データを、testing()でテストデータを抽出できます。
特徴量エンジニアリング
次に特徴量エンジニアリングのためのレシピを作成します。作成したレシピを実際にデータに適用することで、処理済みのデータが作成されます。
どのような特徴量エンジニアリングが必要かは機械学習モデルに依存しますが、今回はRandom Forestを用いるための最低限の前処理を行います。
rec = recipe(price ~ carat + color + cut + clarity + depth, data = df_train) %>% step_log(price) %>% step_ordinalscore(all_nominal()) > rec Data Recipe Inputs: role #variables outcome 1 predictor 9 Operations: Log transformation on price Scoring for all_nominal()
特徴量エンジニアリングのための関数はrecipeにまとめられています。
まずはrecipe()で用いる特徴量とデータを指定し、そこに処理をパイプで重ねていきます。
関数は全てstep_*()で統一されていて、今回用いたstep_log()は指定した変数の対数をとる関数で、step_ordinalscore()は順序のあるカテゴリカル変数を順序（1, 2, 3..）に変換します。
（tree系のモデルではなく線形モデルを用いる場合はstep_dummy()でダミー変数にできます。 *1
all_nominal()は全てのカテゴリカル変数を指定しています。他にもたとえばall_numeric()は全ての連続変数を指定できますし、starts_with()などのtidyverseではお馴染みの関数で指定することもできます。詳細は?selectionsをご確認下さい。
レシピができたので実際にデータに適用しましょう。
rec_preped = rec %>% prep(df_train) df_train_baked = rec_preped %>% juice() > df_train_baked # A tibble: 43,153 x 10 carat cut color clarity depth table x y z price <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 0.23 5 2 2 61.5 55 3.95 3.98 2.43 5.79 2 0.21 4 2 3 59.8 61 3.89 3.84 2.31 5.79 3 0.23 2 2 5 56.9 65 4.05 4.07 2.31 5.79 4 0.290 4 6 4 62.4 58 4.2 4.23 2.63 5.81 5 0.31 2 7 2 63.3 58 4.34 4.35 2.75 5.81 6 0.24 3 7 6 62.8 57 3.94 3.96 2.48 5.82 7 0.24 3 6 7 62.3 57 3.95 3.98 2.47 5.82 8 0.26 3 5 3 61.9 55 4.07 4.11 2.53 5.82 9 0.22 1 2 4 65.1 61 3.87 3.78 2.49 5.82 10 0.23 3 5 5 59.4 61 4 4.05 2.39 5.82 # ... with 43,143 more rows
prep()に訓練データを指定し、訓練データの変換を行います。
さらにjuice()で変換されたデータを抽出することができます。
priceに対数がとられてカテゴリカル変数が（cut, color, clarity）が順序になっているのがわかるかと思います。
次に同じレシピを用いてテストデータも変換します。
df_test_beked = rec_preped %>% bake(df_test) > df_test_baked # A tibble: 10,787 x 10 carat cut color clarity depth table price x y z <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 0.22 4 3 3 60.4 61 5.83 3.88 3.84 2.33 2 0.3 3 7 3 62.7 59 5.86 4.21 4.27 2.66 3 0.23 3 2 4 63.8 55 5.86 3.85 3.92 2.48 4 0.23 3 1 4 60.5 61 5.88 3.96 3.97 2.4 5 0.23 3 2 5 59.5 58 6.00 4.01 4.06 2.4 6 0.23 2 3 5 58.2 59 6.00 4.06 4.08 2.37 7 0.33 5 6 2 61.8 55 6.00 4.49 4.51 2.78 8 0.33 5 7 3 61.1 56 6.00 4.49 4.55 2.76 9 0.3 3 6 3 62.6 57 6.00 4.25 4.28 2.67 10 0.3 4 1 3 62.6 59 6.31 4.23 4.27 2.66 # ... with 10,777 more rows
レシピを新しいデータに適用する場合はbake()を使います*2。
モデルの学習
recipeで訓練データ/テストデータをモデルに投入できる形に変換できたので、次はparsnipで学習モデルを定義しましょう。
Rで機械学習を行う際、学習アルゴリズムごとにパッケージが分かれていて、それ故シンタックスが統一されておらず、各々のシンタックスを覚えるのが大変でした（たとえば入力はformulaなのかmatrixなのかなど）。
これを統一的なシンタックスで扱えるようにするメタパッケージとして、たとえばcaretがありましたが、最近よりtidyverseに馴染むパッケージとしてparsnipが開発されました*3。 ちなみに開発者はcaretと同じくMax Kuhnです。
rf = rand_forest(mode = "regression", trees = 50, min_n = 10, mtry = 3) %>% set_engine("ranger", num.threads = parallel::detectCores(), seed = 42) > rf Random Forest Model Specification (regression) Main Arguments: mtry = 3 trees = 50 min_n = 10 Engine-Specific Arguments: num.threads = parallel::detectCores() seed = 42 Computational engine: ranger
rand_forest()は学習アルゴリズムとしてRandom Forestを使用します。今回は回帰問題なのでmode=regressionとしています。
rand_forest()で指定できるハイパーパラメータは3つあり、それぞれ
trees：作成する決定木の数。min_n：ノードに最低限含まれるサンプル数。mtry：各ノードで分割に用いる特徴量の数。
となります。今回は適当に指定しました。
さらにset_engine()でRandom Forestに用いるパッケージを指定でき、ここではranger を指定しています。
set_engine()内でranger固有のパラメータを指定でき、コア数とシードを指定しました。
モデルが定義できたので、学習を行います。これはfit()にモデル式と入力データを指定するだけで完結します。
fitted = rf %>% fit(price ~ ., data = df_train_baked) > fitted parsnip model object Ranger result Call: ranger::ranger(formula = formula, data = data, mtry = ~3, num.trees = ~50, min.node.size = ~10, num.threads = ~parallel::detectCores(), seed = ~42, verbose = FALSE) Type: Regression Number of trees: 50 Sample size: 43153 Number of independent variables: 9 Mtry: 3 Target node size: 10 Variable importance mode: none Splitrule: variance OOB prediction error (MSE): 0.009112244 R squared (OOB): 0.9911622
モデルの精度評価
学習が済んだので、最後に予測精度を評価しましょう。
テストデータへの予測はpredict() で行います。
df_result = df_test_baked %>% select(price) %>% bind_cols(predict(fitted, df_test_baked)) > df_result # A tibble: 10,787 x 2 price .pred <dbl> <dbl> 1 5.83 5.95 2 5.86 5.98 3 5.86 5.97 4 5.88 5.97 5 6.00 6.00 6 6.00 5.99 7 6.00 6.14 8 6.00 6.14 9 6.00 6.03 10 6.31 6.33 # ... with 10,777 more rows
予測結果に対してyardsticの評価指標関数を適用することでモデルの性能が評価できます。
yardsticには様々な評価指標関数が含まれており、個別で指定することもできますが、matrics()は回帰モデルの結果に対してRMSE、MAE、R2をまとめて計算してくれます。
> df_result %>% + metrics(price, .pred) # A tibble: 3 x 3 .metric .estimator .estimate <chr> <chr> <dbl> 1 rmse standard 0.0910 2 rsq standard 0.992 3 mae standard 0.0671
まとめ
本記事ではtidymodelsによる
- rsampleで訓練データ/テストデータの分割し、
- recipeで特徴量エンジニアリングを行い、
- parsnipでモデルを定義・学習し、
- yardsticで性能を評価
する、tidyな機械学習フローを紹介しました。
各パッケージにはまだまだ多くの機能があり、今回紹介できたのは最低限の機能のみです。
次回はtidymodelsで交差検証を行い、ハイパーパラメータをチューニングする方法をご紹介できればと思っています。
参考文献
- tidymodels
- rsample
- recipe
- parsnip
- yardstic
*1: recipeの他のstep関数やselectorは以下をご確認下さい。Function reference • recipes
*2:今回tree系のモデルなので標準化はやっていませんが、たとえば標準化を行う場合は訓練データの平均と分散を用いてテストデータの標準化が行われます。
*3:parsnipが対応しているパッケージは以下を御覧ください：List of Models • parsnip
Ridge回帰の使い所を考える
はじめに
本記事では回帰係数の推定方法としてのRidge回帰（L2正則化）の使い所について考えたいと思います。
結論から言うと、サンプルサイズが不十分な状況下でRidge回帰を用いてより真の値に近い回帰係数を得る確率を高めるという使い方ができそうです。
OLS推定量とRidge推定量
簡単のために、以下の線形回帰モデルを考えます。
\begin{align}
Y = X\beta + \epsilon
\end{align}
ここで、は
の行列で、
は
の列ベクトルとします。
はサンプルサイズ、
は特徴量の数になります。
このとき、OLS推定量は二乗誤差を最小にするようなで与えられます。
\begin{align}
\widehat{\beta}_{\mathrm{OLS}} &= \arg\min_b \left(Y - X b\right)^\top\left(Y - X b\right) \\
&= \left(X^\top X \right)^{-1}X^\top Y
\end{align}
標準的な仮定のもとで、OLS推定量は推定量 は真の値
の一致推定量となることが知られています。*1
これに対して、 のとれる値に制限をかけることで過学習を抑えようとするのがRidge回帰になります。
具体的にはOLSの目的関数に正則化項を付け加えたものを最小化します。
は正則化の強さを表すハイパーパラメータになります。
\begin{align}
\widehat{\beta}_{\mathrm{Ridge}} &= \arg\min_b \left(Y - X b\right)^\top\left(Y - X b\right) + \lambda||b||^2 \\
&= \left(X^\top X + \lambda I_K \right)^{-1}X^\top Y
\end{align}
ここでは
の単位行列です。
と
を比較すれば明らかですが、両者の違いは
のみです。
なので、のときOLS推定量とRidge推定量は完全に一致します。反対に
を大きくするとRidge推定量は0方向に引っ張られます。
つまり、Ridge回帰はパラメータが極端な値を取らないよう安定させる（過学習を抑える）というリターンを得るために、（多少の）バイアスがかかるという犠牲を払うモデルであると言えます。
これは機械学習の文脈でバイアスとバリアンスのトレードオフとしてよく知られています。
シミュレーション
それでは実際にOLSとRidge回帰の挙動の違いについてシミュレーションで確認してみましょう。
まずは必要なパッケージを読み込みます。*2
import pandas as pd import numpy as np from sklearn.linear_model import LinearRegression, Ridge, RidgeCV from sklearn.preprocessing import StandardScaler from sklearn.datasets import make_spd_matrix #分散共分散行列を作る from plotnine import * # 可視化用。
ちょっと極端な方が結果がクリアに見えると思うので、サンプルサイズは100しかないんだけど80個の特徴量について回帰係数を推定したい、というシチュエーションを想定することにします。
下記のコードではランダムにデータを生成してOLSとRidge回帰で回帰係数を推定する、ということを繰り返しています。
np.random.seed(42) I = 1000 N = 100 K = 80 beta = np.random.uniform(-1, 1, K) betas = np.zeros((I, 2, K)) for i in range(I): X = np.random.multivariate_normal(np.zeros(K), make_spd_matrix(K), N) X = StandardScaler().fit_transform(X) U = np.random.normal(0, 1, N) Y = X @ beta + U if i == 0: #一回目のループでRidgeのハイパーパラメータを決める ridge_alpha = RidgeCV(alphas=np.random.uniform(0.1, 20, 100), cv=5).fit(X, Y).alpha_ print(f'ridge: {ridge_alpha}') beta_ols = LinearRegression().fit(X, Y).coef_ beta_ridge = Ridge(alpha=ridge_alpha).fit(X, Y).coef_ betas[i, :, :] = np.concatenate([beta_ols[:, np.newaxis], beta_ridge[:, np.newaxis]], axis=1).T if i % 200 == 0: print(i) # 可視化のためにシミュレーションの結果をデータフレームとして持ち直す df = None for i in range(I): tmp = (pd.DataFrame(data = betas[i, :, :].T, columns=['OLS', 'Ridge']) .assign(Truth = beta, k = lambda df: df.index, i = i)) if df is None: df = tmp else: df = df.append(tmp) df = df.melt(id_vars=['i', 'k', 'Truth'], value_vars=['OLS', 'Ridge'], value_name='coefficient', var_name='method')
シミュレーションの結果を可視化してみましょう。
g = (ggplot(aes('factor(k)', 'coefficient', fill='method', color='method'), df[df.k.between(1, 10)]) + geom_hline(yintercept=0, size=1, color='gray') + geom_violin(alpha=0.5) + geom_point(aes('factor(k)', 'Truth'), size=4, color='black', fill='black') + scale_fill_brewer(name='Method', type='qual',palette='Dark2') + scale_color_brewer(name = 'Method', type='qual',palette='Dark2') + lims(y=(-3, 3)) + labs(x='beta_hat(k)', y='Coefficient') + theme_bw() + theme(figure_size=(15, 7))) print(g)
図は初めの10個の回帰係数を取り出して、シミュレーション1000回分の回帰係数の分布を確認しています。
縦軸が回帰係数の値、横軸が回帰係数の番号です。
真ん中の黒い点が真の回帰係数の値で、緑がOLSで求めた回帰係数の分布、オレンジがRidge回帰での分布になります。
図からは各手法で求めた回帰係数について、以下の性質が見て取れます。
- OLSの場合はバイアスがなく真の値を中心に分布しているが、その分散は大きく、見当はずれの値を取ることもある
- Ridge回帰の場合はゼロ方向にバイアスがかかっている。その一方で回帰係数は非常に安定しており、真の値に近い部分に密度高く分布しており、極端な値を取ることもない。
Ridge回帰の使いどころ
数式とシミュレーションの結果から、サンプルサイズが不十分な状況下でRidge回帰を用いてより真の値に近い回帰係数を得る確率を高めるという使い方が考えられそうです。*3
特に、特徴量のインパクトを見える際により保守的な推定をしたい場合には有効と言えそうです。Ridge回帰は極端な回帰係数をとらず、またRidge回帰のバイアスは0方向にかかるので、Ridge回帰ではより控え目なインパクトが推定されると言えるかと思います。
*1:https://www.ssc.wisc.edu/~bhansen/econometrics/
*2:ちなみにplotnineはRの代表的な可視化パッケージであるggplot2と同じ記法が使えるpython用のパッケージです。若干制限もありますが、ggplot2に慣れたRユーザーにとっては使いやすいかと思います。
*3:サンプルサイズが十分な状況ではOLSで一致推定を行うのが妥当かと思います。
purrrとbroomで複数の回帰モデルを効率的に管理する
私は探索的にデータを見てく段階では、可視化に加えて複数の回帰モデルを作成して比較をする、ということをよくやっています。
モデルの数が少ない場合は個別にモデルを作成してsummary()で見ていく事もできますが、モデルの数が増えるにつれてそのやり方では管理が難しくなってきます。
そこで、本記事では、purrrのmap()とbroomのtidy(), glance()を用いて複数の回帰モデルを効率的に扱う方法を紹介したいと思います。
まずはライブラリを読み込みます。tidyverseはおなじみのデータハンドリングと可視化のためのパッケージ群です。tidymodelsはモデリングをtidyなやり方で統一的に扱えるようにするパッケージ群になります。今回はbroomのみ用いますが、後日他のパッケージの紹介記事も書ければと思っています。
library(tidyverse) library(tidymodels) # 可視化用 theme_scatter = theme_minimal() + theme(panel.grid.minor = element_blank(), axis.text = element_text(color = "black")) theme_minimal2 = theme_scatter + theme(panel.grid.major.x = element_blank())
データセットはこれまたおなじみのdiamondを使います。
詳細は?diamondsを見て頂ければと思いますが、ダイヤの重さ(carat)や透明度(clarity)とその値段(price)などが入ったデータセットになります。
今回はこのデータセットを用いて、ダイヤの属性がダイヤの値段にどんな影響を与えるのかを探索することにします。
可視化によるデータ分析
まずは非説明変数である価格の分布を見てみましょう。
df = diamonds #面倒なので名前をdfに df %>% ggplot(aes(price)) + geom_histogram(fill = "#56B4E9", color = "white") + theme_minimal2
ものすごく歪んだ分布をしています。 歪んだ分布への簡便な対応策として、今回は対数を取ることにします。
df %>% ggplot(aes(log(price))) + geom_histogram(fill = "#56B4E9", color = "white") + theme_minimal2
きれいな正規分布とまではいきませんが先程よりは中心によった分布になりました。価格は対数を取って分析することにします。
次にダイヤの属性と値段の関係を見ていきましょう。 たとえばダイヤが重ければ重いほど（大きければ大きいほど）値段は高くなりそうです。
df %>% sample_frac(0.1) %>% #データが多いので減らす ggplot(aes(carat, log(price))) + geom_point(color = "#0072B2", alpha = 0.5) + theme_scatter
実際に可視化してみると、どうやらこの仮説は正しそうです。
ただ、carat とlog(price)の関係は非線形に見えます。
caratがlog(price)に与えるインパクトはcaratが大きくなるにつれて逓減するという関係を反映するために、caratにも対数をとってみましょう。
df %>% sample_frac(0.1) %>% #データが多いので減らす ggplot(aes(log(carat), log(price))) + geom_point(color = "#0072B2", alpha = 0.5) + theme_scatter
両変数に対数をとることで線形の関係が構築できました！ これなら線形回帰でうまくモデリングできそうです。
同様に透明度と価格の関係も見てみましょう。 直感的には透明度が高ければ高いほど価格が高くなる傾向がありそうに思えます。
df %>% ggplot(aes(clarity, log(price))) + geom_boxplot(fill = "#56B4E9") + theme_minimal2
直感に反する結果となりました。上のグラフは右に行くほど透明度のランクが高いのですが、透明度が高いほどダイヤが安くなっています。
これは一体どういうことでしょうか？ 透明度とカラット数の関係をみることでこの謎が解けます。
df %>% ggplot(aes(clarity, log(carat))) + geom_boxplot(fill = "#56B4E9") + theme_minimal2
透明度が高い場合にはカラット数が小さくなる傾向が見て取れます。
僕はダイヤモンドに詳しくないのですが、大きくて透明なダイヤを作るのは難しいということのようです。*1
つまり、ここでは「ダイヤのカラット数と価格には正の相関がある一方で、カラット数と透明度には負の相関があり、結果として透明度と価格に負の相関があるように見えてしまう」という典型的な交絡の問題が起きています。 この場合は、重回帰分析でカラット数の影響を取り除くことで透明度と価格の正しい関係を見ることができます。 今回のように、単純に一つ一つの説明変数と被説明変数の関係を見ているだけでは間違った結論を下していまう可能性があるので注意が必要です。
モデリングによるデータ分析
ここまでの可視化で得られた仮説を回帰モデルで確かめてみましょう。
回帰モデルを作る前に、カテゴリカル変数がordere=TRUEになっているとlm()の挙動が面倒なので、factor型の変数は全てordered=FALSEにしておきます。
df_input = df %>% mutate_if(is.ordered, factor, ordered = FALSE)
mutate_if()は複数の列に同じ処理を行う際に便利な関数です。
第一引数で条件を指定して（is.ordered）、条件に当てはまった列にのみ第二引数の関数を適用します（factor）。関数のオプションは後ろにくっつけて指定すればOKです（ordered = FALSE）。
複数のモデルを当てはめる場合
今回は
- 価格を透明度のみで説明するモデル
- 透明度とカラット数で説明するモデル
- 上と同じだが、カラット数に対数をとったモデル
の3つのモデルを作成してみます。
formulas = c(log(price) ~ clarity, log(price) ~ clarity + carat, log(price) ~ clarity + log(carat)) %>% enframe("model_no", "formula") formulas # A tibble: 3 x 2 model_no formula <int> <list> 1 1 <S3: formula> 2 2 <S3: formula> 3 3 <S3: formula>
enframe()はベクトルをデータフレームにしてくれる関数です。データフレームとして持っていたほうが後々管理がしやすいと思うので、変換しておきました。
上のデータフレームにmap()を適用することで複数の回帰モデルを一気に作成できます。
df_result = formulas %>% mutate(model = map(formula, lm, data = df_input), #(1) tidied = map(model, tidy), #(2) glanced = map(model, glance)) #(3) df_result # A tibble: 3 x 5 model_no formula model tidied glanced <int> <list> <list> <list> <list> 1 1 <S3: formula> <S3: lm> <tibble [8 × 5]> <tibble [1 × 11]> 2 2 <S3: formula> <S3: lm> <tibble [9 × 5]> <tibble [1 × 11]> 3 3 <S3: formula> <S3: lm> <tibble [9 × 5]> <tibble [1 × 11]>
(1)ではmap()を用いて「formula列の一つ一つの値を引数としてlm()を実行する」ということをやっています。lm()はデータも指定する必要がありますが、map()もmutate_if()同様後ろにくっつけて指定できます。
さらに、(2), (3)では推定されたモデルに対して、それぞれtidy()とglance()を適用しています。
分析結果をデータフレームとして持つことで、どのモデルがどの結果に対応するかを間違うことなく効率的に管理することができます。
tidy()は回帰モデルの係数をtidyなデータフレームとして持ってきてくれる関数です。データフレームの中のデータフレームを取り出すためにunnest()を使います。
df_coef = df_result %>% select(model_no, tidied) %>% unnest() %>% mutate_if(is.double, round, digits=2) df_coef # A tibble: 26 x 6 model_no term estimate std.error statistic p.value <int> <chr> <dbl> <dbl> <dbl> <dbl> 1 1 (Intercept) 8.03 0.04 221. 0 2 1 claritySI2 0.14 0.04 3.69 0 3 1 claritySI1 -0.18 0.04 -4.81 0 4 1 clarityVS2 -0.26 0.04 -7.08 0 5 1 clarityVS1 -0.3 0.04 -7.99 0 6 1 clarityVVS2 -0.5 0.04 -12.8 0 7 1 clarityVVS1 -0.7 0.04 -17.7 0 8 1 clarityIF -0.62 0.04 -14.4 0 9 2 (Intercept) 5.36 0.01 372. 0 10 2 claritySI2 0.570 0.01 40.1 0 # ... with 16 more rows
3つのモデルを比較してみましょう。
モデルを横に並べるためにspread()を使います。
そのままだと回帰係数がアルファベット順になってしまうので、並び順を維持するためにfct_inorder()を使っています（出てきた順に並ぶ）。
df_coef %>% mutate(term = fct_inorder(term)) %>% select(model_no, term, estimate) %>% spread(model_no, estimate) # A tibble: 10 x 4 term `1` `2` `3` <fct> <dbl> <dbl> <dbl> 1 (Intercept) 8.03 5.36 7.77 2 claritySI2 0.14 0.570 0.48 3 claritySI1 -0.18 0.72 0.62 4 clarityVS2 -0.26 0.82 0.78 5 clarityVS1 -0.3 0.86 0.82 6 clarityVVS2 -0.5 0.93 0.98 7 clarityVVS1 -0.7 0.92 1.03 8 clarityIF -0.62 1 1.11 9 carat NA 2.08 NA 10 log(carat) NA NA 1.81
モデル1は透明度のみのモデルであり、カラット数の影響を取り除いていないので、透明度と価格には負の関係があるように見えています。 その一方で、モデル2はカラット数を変数として加えることで影響を取り除いた上での透明度と価格の影響を見ています。モデル2の係数を見ると透明度が高いほど価格が高くなる傾向が見て取れ、こちらのほうがより尤もらしい結果であると言えそうです。
モデル3はカラット数に対数をとったモデルになります。 モデル2もモデル3も透明度と価格の関係には大きな変化はなさそうです。
glance()はモデルの性能をtidyなデータフレームとしてまとめてくれる関数です。
df_result %>% select(model_no, glanced) %>% unnest() %>% mutate_if(is.double, round, digits=2) # A tibble: 3 x 12 model_no r.squared adj.r.squared sigma statistic p.value df logLik <int> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> 1 1 0.05 0.05 0.99 415. 0 8 -75907. 2 2 0.87 0.87 0.37 43737. 0 9 -23024. 3 3 0.97 0.97 0.19 187918. 0 9 13378. # ... with 4 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>, # df.residual <int>
様々な指標がまとめて出力されますが、今回は自由度調整済み決定係数(adj.r.squared)を見ることにしましょう。
透明度のみモデル（モデル1）と比べると、カラット数を加えることでモデルの説明力は大幅に上昇しています（モデル2）。
カラット数に対数をとるとさらに説明力が改善されており（モデル3）、3つのモデルの中ではモデル3がベストであると言えそうです。
サブサンプルに分けて分析する場合
先程は複数のモデルを同じデータに当てはめましたが、同じモデルを複数のデータにまとめてフィットさせることもできます。 たとえば透明度によって価格とカラット数の関係が異なるのかを調べてみることにしましょう。 透明度の低いダイヤは大きくてもあまり価値は上がらないが、透明度が高い場合はサイズが大きくなると価値が大きく上昇する、といった関係があるかもしれません。 これは透明度でデータをサブサンプルに分割して、各データに対して回帰モデルを当てはめることで確認できます。
df_nested = df %>% group_by(clarity) %>% #nest(-clarity)でも同じ nest() %>% arrange(clarity) df_nested # A tibble: 8 x 2 clarity data <ord> <list> 1 I1 <tibble [741 × 9]> 2 SI2 <tibble [9,194 × 9]> 3 SI1 <tibble [13,065 × 9]> 4 VS2 <tibble [12,258 × 9]> 5 VS1 <tibble [8,171 × 9]> 6 VVS2 <tibble [5,066 × 9]> 7 VVS1 <tibble [3,655 × 9]> 8 IF <tibble [1,790 × 9]>
nest()はグループごとにデータフレームを分割してdata列に格納する関数です。
data列にmap()を適用することで各サブサンプルに対して同じモデルをまとめてフィットさせることができます。
df_result2 = df_nested %>% mutate(model = map(data, ~lm(log(price) ~ log(carat), data = .)), tidied = map(model, tidy))
モデルが作成できたので、回帰係数を見てみましょう。
df_result2 %>% select(clarity, tidied) %>% unnest() %>% filter(term != "(Intercept)") %>% #切片は無視 mutate_if(is.double, round, digits=2) # A tibble: 8 x 6 clarity term estimate std.error statistic p.value <ord> <chr> <dbl> <dbl> <dbl> <dbl> 1 I1 log(carat) 1.53 0.02 84.4 0 2 SI2 log(carat) 1.79 0 557. 0 3 SI1 log(carat) 1.82 0 676. 0 4 VS2 log(carat) 1.78 0 563. 0 5 VS1 log(carat) 1.83 0 473. 0 6 VVS2 log(carat) 1.85 0.01 338. 0 7 VVS1 log(carat) 1.83 0.01 245. 0 8 IF log(carat) 1.87 0.01 169. 0
I1だけは係数が小さくなっており、その他はほとんど同じ係数です。
どうやらと透明度が最低ランクの場合は、他のランクと比べてカラット数が価格に与える影響は小さいようです。
散布図も確認すると、確かにI1のみ傾きが少し緩やかになっているような気がします。
df %>% sample_frac(0.3) %>% #データが多いので減らす ggplot(aes(log(carat), log(price))) + geom_point(color = "gray", alpha = 0.5) + facet_wrap(~clarity) + geom_smooth(method = "lm", se = F, color = "#0072B2") + theme_scatter
まとめ
このように、探索的にデータを見ていく際には、purrrとbroomを使うと複数の回帰モデルを効率的に比較することができ、とても便利です。
